diff --git a/ColumnGeneration/base.cpp b/ColumnGeneration/base.cpp index 5d779a2d4f157f28017b554b49041d0ee5e2e7ff..0ea60116e31c71a039b66f8dcf053431fbee7f5e 100644 --- a/ColumnGeneration/base.cpp +++ b/ColumnGeneration/base.cpp @@ -1,3 +1,5 @@ +#include "base.h" +#include "base.h" /******************** Include eigenes Modul ********************************/ #include "base.h" @@ -42,70 +44,85 @@ c_Controller c_Controller::c_Controller( string settings_file ) : c_SettingsReader( settings_file ), - next_node_id( 0 ), - ub( 10.0 * INFTY ), - aspiration_lb( INFTY ), - aspiration_ub( -INFTY ), - env( NULL ), - p_rmp( NULL ), - p_pool( NULL ), - p_root_node( NULL ), - p_conv_info( NULL ), - cols_in_rmp(0), - last_node_terminated( false ), - last_node_integer( false ), + env(nullptr), + p_pool(nullptr), + p_red_costs_fixing(nullptr), + p_currentNode(nullptr), + p_conv_info(nullptr), + next_node_id(0), + ub(100.0 * INFTY), + aspiration_lb(INFTY), + aspiration_ub(-INFTY), + last_node_terminated(false), + last_node_integer(false), + // protected members + p_root_node(nullptr), + p_rmp(nullptr), + p_best_solution(nullptr), + info_number_of_rmp_iterations( 0 ), info_number_of_generated_columns( 0 ), info_number_of_solved_pricing_problems( 0 ), info_number_of_generated_cuts( 0 ), info_number_of_solved_nodes( 0 ), + + first_integer_solution(INFTY), info_time_first_integer_solution( INFTY ), info_time_best_integer_solution( INFTY ), - info_time_pricing( 0.0 ), - info_time_separation( 0.0 ), - info_time_reoptimization( 0.0 ), - info_time_root_branching( 0.0 ), - info_time_Init(0.0), - info_time_AddColumn(0.0), - info_time_NodeTerm(0.0), - info_time_NodeUpdateUB(0.0), - info_timer_nodeSolve(0.0), - info_timer_ContUpdate(0.0), - info_timer_CreateDual(0.0), - info_timer_OrgPool(0.0), - info_timer_GetDual(0.0), - info_timer_RMPUpdate(0.0), - info_timer_PricingUpdate(0.0), - info_time_MIPForUB(0.0), + + info_time_overall(0.0), + info_time_pricing(0.0), + info_time_separation(0.0), + info_time_reoptimization(0.0), + info_time_root_lb1(0.0), + info_time_root_lb2(0.0), + info_time_root_branching(0.0), + info_time_AddColumn(0.0), + info_time_Init(0.0), + info_time_NodeTerm(0.0), + info_time_NodeRun(0.0), + info_time_NodeUpdateUB(0.0), + info_time_MIPForUB(0.0), + info_timer_nodeSolve(0.0), + info_timer_ContUpdate(0.0), + info_timer_CreateDual(0.0), + info_timer_OrgPool(0.0), + info_timer_GetDual(0.0), + info_timer_RMPUpdate(0.0), + info_timer_PricingUpdate(0.0), + info_timer_StrongBranching(0.0), + info_timer_last_PP_root(0.0), + info_1( &cout ), info_2( &cout ), info_3( &cout ), info_err( &cout ), info_statistics( &cout ), - p_best_solution( NULL ), - first_integer_solution(INFTY), // Parameter DEBUG( "Debug", false, false ), UseLPOptimizer( "UseLPOptimizer", "dual", false ), ReloadSettings( "ReloadSettings", false, false ), Stabilization( "Stabilization", false, false ), + InfoLevel( "InfoLevel", 1, false ), InfoStatisticsLevel( "InfoStatisticsLevel", 1, false), InfoDottyOutput( "InfoDottyOutput", false, false ), InfoGnuplotOutput( "InfoGnuplotOutput", false, false ), InfoGraphMLOutput( "InfoGraphMLOutput", false, false ), + InfoTikzOutput("InfoTikzOutput", false, false), InfoDottySuppressIntermediateNodes( "InfoDottySuppressIntermediateNodes", false, false ), + InfoInstance("InfoInstance", "Any Instance", false), InfoPrintOverallBestSolution( "InfoPrintOverallBestSolution", true, false ), InfoPrintBestSolution( "InfoPrintBestSolution", true, false ), InfoOutputFileName( "InfoOutputFileName", "Solutions.csv", false ), InfoConvergency( "InfoConvergency", false, false ), InfoStatisticsColumns( "InfoStatisticsColumns", false, false ), + MaxSolutionTimeSec( "MaxSolutionTimeSec", INFTY, false ), MaxColumnsToAddInPricing( "MaxColumnsToAddInPricing", 500, false ), PricingWithColumnSelection( "PricingWithColumnSelection", true, false ), PricingHierarchyMaxNumFailuresBeforeSwitchOff( "PricingHierarchyMaxNumFailuresBeforeSwitchOff", 999999, false ), PricingHierarchyMinNumColsToGenerate( "PricingHierarchyMinNumColsToGenerate", 1, false ), - PricingHierarchyStandardOrder("PricingHierarchyStandardOrder", true, false), MaxRatioColsToRows( "MaxRatioColsToRows", 5.0, false ), MaxRatioColsToRowsMIP( "MaxRatioColsToRowsMIP", 5.0, false ), Branching( "Branching", true, false ), @@ -116,20 +133,31 @@ c_Controller::c_Controller( string settings_file ) StrongBranchingNumCandidates( "StrongBranchingNumCandidates", 10 , false ), StrongBranchingMinNumCandidates( "StrongBranchingMinNumCandidates", 5, false ), StrongBranchingNumCandidatesDecreasePerLevel( "StrongBranchingNumCandidatesDecreasePerLevel", 0, false ), - StrongBranchingSolveExact("StrongBranchingSolveExact", true, false ), + CuttingPlanes( "CuttingPlanes", true, false ), - InfoInstance( "InfoInstance", "Any Instance", false ), + UseMIPSolverForUBs( "UseMIPSolverForUBs", false, false ), UseMIPSolverUptoLevel( "UseMIPSolverUptoLevel", 0, false ), + UseMIPSolverUpFromLevel("UseMIPSolverUpFromLevel", 0, false), MIPMaxSolutionTimeSec( "MIPMaxSolutionTimeSec", 60, false ), + MIPScreenOutput("MIPScreenOutput", false, false), + MIPAccelerate("MIPAccelerate", false, false), + SolveAsMipAfterTimeOut("SolveAsMipAfterTimeOut", false, false), + UseLagrangeanLB( "UseLagrangeanLB", false, false ), ReducedCostVariableFixing( "ReducedCostVariableFixing", true, false ), StabFrequence("StabFrequence",5,false), FeasyCutFrequence ("FeasyCutFrequence",0,false), FeasyCutGap("FeasyCutGap",1.00,false), NumThreadsRMP("NumThreadsRMP",1,false), - info_timer_last_PP_root(0.0), - i_startHierarchyId(0) + NameDottyOutput("NameDottyOutput", "bbtree.dty", false), + CheckForInftyAndLBError("CheckForInftyAndLBError",false,false), + + StrongBranchingRule("StrongBranchingRule", "BestSmall", false), + StrongBranchingMu("StrongBranchingMu", 0.0, false), + PricersUsedForStrongBranching("PricersUsedForStrongBranching", MaxNumPricers, false), + LocalCG_EPS("LocalCG_EPS", CG_EPS, false), + LocalInfty("LocalInfty", INFTY, false) { env = new CPLEXEnvironment(); p_pool = new c_ColumnPool( this ); @@ -144,6 +172,7 @@ c_Controller::c_Controller( string settings_file ) AddParameter( &InfoDottyOutput ); AddParameter( &InfoGnuplotOutput ); AddParameter( &InfoGraphMLOutput ); + AddParameter(&InfoTikzOutput); AddParameter( &InfoDottySuppressIntermediateNodes ); AddParameter( &InfoPrintOverallBestSolution ); AddParameter( &InfoPrintBestSolution ); @@ -155,7 +184,6 @@ c_Controller::c_Controller( string settings_file ) AddParameter( &PricingWithColumnSelection ); AddParameter( &PricingHierarchyMaxNumFailuresBeforeSwitchOff ); AddParameter( &PricingHierarchyMinNumColsToGenerate ); - AddParameter( &PricingHierarchyStandardOrder ); AddParameter( &MaxRatioColsToRows ); AddParameter( &MaxRatioColsToRowsMIP ); AddParameter( &Branching ); @@ -166,7 +194,6 @@ c_Controller::c_Controller( string settings_file ) AddParameter( &StrongBranchingNumCandidates ); AddParameter( &StrongBranchingMinNumCandidates ); AddParameter( &StrongBranchingNumCandidatesDecreasePerLevel ); - AddParameter( &StrongBranchingSolveExact ); AddParameter( &CuttingPlanes ); AddParameter( &InfoInstance ); AddParameter( &UseMIPSolverForUBs ); @@ -178,6 +205,17 @@ c_Controller::c_Controller( string settings_file ) AddParameter( &FeasyCutFrequence); AddParameter( &FeasyCutGap); AddParameter( &NumThreadsRMP); + AddParameter( &NameDottyOutput); + AddParameter( &PricersUsedForStrongBranching); + AddParameter( &StrongBranchingRule); + AddParameter( &StrongBranchingMu); + AddParameter( &LocalInfty); + AddParameter( &LocalCG_EPS); + AddParameter( &SolveAsMipAfterTimeOut); + AddParameter(&CheckForInftyAndLBError); + AddParameter(&UseMIPSolverUpFromLevel); + AddParameter(&MIPScreenOutput); + AddParameter(&MIPAccelerate); // AddParameter( ); } @@ -185,49 +223,48 @@ c_Controller::c_Controller( string settings_file ) c_Controller::~c_Controller() { // clean memory at the end - for ( auto hierarchy = pricing_problem_hierarchy.begin(); hierarchy != pricing_problem_hierarchy.end(); ++hierarchy ) - delete *hierarchy; + for ( auto hierarchy : pricing_problem_hierarchy ) + delete hierarchy; pricing_problem_hierarchy.clear(); - + // colums for ( int i=0; i<(int)cols_in_rmp.size(); i++ ) delete cols_in_rmp[i]; - // separation problems ... - for ( auto sep_prob=separation_problems.begin(); sep_prob!=separation_problems.end(); ++sep_prob ) - delete *sep_prob; + for ( auto sep_prob : separation_problems ) + delete sep_prob; separation_problems.clear(); - // nodes of the B&B-Tree DeleteTree(); - // constraints - for ( auto con_iter=l_constraints.begin(); con_iter!=l_constraints.end(); ++con_iter ) - delete *con_iter; - + for ( auto con_iter : l_constraints ) + delete con_iter; // basic objects // NOT: "delete p_root_node" delete p_rmp; - delete env; delete p_pool; delete p_best_solution; + delete p_currentNode; + env->terminate(); + delete env; } void c_Controller::DeleteTree() { // remove existing tree - for ( auto node=active_nodes.begin(); node !=active_nodes.end(); ++node ) - delete *node; + for ( auto node : active_nodes ) + delete node; active_nodes.clear(); - for ( auto node=inactive_nodes.begin(); node != inactive_nodes.end(); ++node ) - delete *node; + for ( auto node : inactive_nodes ) + delete node; inactive_nodes.clear(); p_root_node = nullptr; next_node_id = 0; // reset bounds - ub = 10.0 * INFTY; + ub = 100.0 * INFTY; // remove solution information delete p_best_solution; p_best_solution = nullptr; + p_currentNode = nullptr; first_integer_solution = INFTY; // reset infos about the solution process last_node_terminated = false; @@ -243,6 +280,12 @@ void c_Controller::DeleteTree() info_time_separation = 0.0 ; info_time_reoptimization = 0.0 ; info_time_root_branching = 0.0 ; + for ( int i=0; i<(int)v_info_time_pricing_level.size(); i++ ) + v_info_time_pricing_level[i] = 0; + for (int i = 0; i<(int)v_info_calls_pricing_level.size(); i++) + v_info_calls_pricing_level[i] = 0; + for (int i = 0; i<(int)v_info_success_pricing_level.size(); i++) + v_info_success_pricing_level[i] = 0; info_time_Init = 0.0; info_time_AddColumn = 0.0; info_time_NodeTerm = 0.0; @@ -290,11 +333,9 @@ c_Controller::e_branching c_Controller::BranchingRule() const return best_first; } - void c_Controller::TerminateNode( c_BranchAndBoundNode* node ) { - list<c_BranchAndBoundNode*>::iterator tester; - for ( tester=active_nodes.begin(); tester!=active_nodes.end(); ++tester ) + for ( auto tester=active_nodes.begin(); tester!=active_nodes.end(); ++tester ) if ( (*tester) == node ) { active_nodes.erase( tester ); @@ -303,15 +344,11 @@ void c_Controller::TerminateNode( c_BranchAndBoundNode* node ) } } - void c_Controller::AddPricingProblemHierarchy( c_PricingProblemHierarchy* hierarchy ) { pricing_problem_hierarchy.push_back( hierarchy ); - pricing_levels_in_hierarchies.push_back(0); } - - void c_Controller::AddPricingProblem( c_PricingProblem* prob ) { int max_num_failures = PricingHierarchyMaxNumFailuresBeforeSwitchOff(); @@ -324,11 +361,9 @@ void c_Controller::AddPricingProblem( c_PricingProblem* prob ) void c_Controller::FixPartOfSolution() {} - void c_Controller::UnfixSolution() {} - class c_GnuPlotNodeInfo { public: double from_bound; @@ -356,18 +391,13 @@ void GnuPlotTree( double from_bound, double from_pos, c_BranchAndBoundNode* node int num_sons = node->NumSons(); double offset = ( interval_rhs - interval_lhs ) / double( num_sons ); double startpos = interval_lhs; - const list<c_BranchAndBoundNode*>& sons = node->Sons(); - list<c_BranchAndBoundNode*>::const_iterator i_son; - for ( i_son=sons.begin(); i_son !=sons.end(); ++i_son ) + for ( auto son : node->Sons()) { - c_BranchAndBoundNode* son = *i_son; GnuPlotTree( gnu_node.bound, gnu_node.pos, son, startpos, startpos+offset, result ); startpos += offset; } } - - void c_Controller::UpdateAllNodeInformations( c_BranchAndBoundNode* selected_node ) { /////////////////// @@ -399,47 +429,47 @@ void c_Controller::UpdateAllNodeInformations( c_BranchAndBoundNode* selected_nod else gnu << "set ytics add (\"opt\" " << lb << ") " << endl; int as = 1; - for ( set<int>::const_iterator ii=line_styles.begin(); ii!=line_styles.end(); ++ii ) + for ( auto ii : line_styles ) { - gnu << "set style arrow " << as++ << " nohead nofilled front lt " << *ii << " lw 1 " << endl; + gnu << "set style arrow " << as++ << " nohead nofilled front lt " << ii << " lw 1 " << endl; } gnu << "set style rect back fc rgb \"white\" fs solid 1.0 noborder " << endl; gnu << "set object 1 rect from -0.02," << lb << " to 1.02," << ub << " fc rgb \"gray\" " << endl; gnu << "plot [-.02:1.02] [" << root_lb-0.1*(ub-root_lb) << ":" << ub+0.3*(ub-root_lb) << "] "; as = 1; - for ( set<int>::const_iterator ii=line_styles.begin(); ii!=line_styles.end(); ++ii ) + for ( int as=0; as<(int)line_styles.size(); as++ ) { - gnu << " \"-\" using 2:1:($4-$2):($3-$1) with vectors as " << as++ << " notitle, \\" << endl; + gnu << " \"-\" using 2:1:($4-$2):($3-$1) with vectors as " << (as+1) << " notitle, \\" << endl; } gnu << " \"-\" using 4:3 with points pt 13 ps .5 notitle, \\" << endl; gnu << " \"-\" using 4:3 with points pt 10 ps .5 notitle " << endl; // now the data - for ( set<int>::const_iterator ii=line_styles.begin(); ii!=line_styles.end(); ++ii ) + for ( const auto ii : line_styles ) { - for ( int i=0; i<(int) result.size(); i++ ) - if ( result[i].line_style == *ii ) - gnu << result[i].from_bound - << " " << result[i].from_pos - << " " << result[i].bound - << " " << result[i].pos + for ( const auto& res : result ) + if ( res.line_style == ii ) + gnu << res.from_bound + << " " << res.from_pos + << " " << res.bound + << " " << res.pos << endl; gnu << "e" << endl; } - for ( int i=0; i<(int) result.size(); i++ ) - if ( result[i].bound < ub ) - gnu << result[i].from_bound - << " " << result[i].from_pos - << " " << result[i].bound - << " " << result[i].pos + for ( const auto& res : result ) + if ( res.bound < ub ) + gnu << res.from_bound + << " " << res.from_pos + << " " << res.bound + << " " << res.pos << endl; gnu << "e" << endl; - for ( int i=0; i<(int) result.size(); i++ ) - if ( result[i].bound >= ub ) - gnu << result[i].from_bound - << " " << result[i].from_pos - << " " << result[i].bound - << " " << result[i].pos + for ( auto res : result ) + if ( res.bound >= ub ) + gnu << res.from_bound + << " " << res.from_pos + << " " << res.bound + << " " << res.pos << endl; gnu << "e" << endl; } @@ -450,173 +480,171 @@ void c_Controller::UpdateAllNodeInformations( c_BranchAndBoundNode* selected_nod { // Get all nodes... list<c_BranchAndBoundNode*> all_nodes; - list<c_BranchAndBoundNode*>::const_iterator node; - for ( node=active_nodes.begin(); node!=active_nodes.end(); ++node ) - all_nodes.push_back( *node ); - for ( node=inactive_nodes.begin(); node!=inactive_nodes.end(); ++node ) - all_nodes.push_back( *node ); + for ( auto node : active_nodes ) + all_nodes.push_back( node ); + for ( auto node : inactive_nodes ) + all_nodes.push_back( node ); // Hier auch die Ausgabe des B&B Suchbaums als .dty file - ofstream bbout( "bbtree.dty" ); + ofstream bbout( NameDottyOutput()); bbout << "digraph bbtree { " << endl; if ( InfoInstance().size() != 0 ) bbout << "title = \"" << InfoInstance() << "\"\n"; // alle Knoten beschriften - for ( node=all_nodes.begin(); node!=all_nodes.end(); ++node ) + for ( auto node : all_nodes ) { ostringstream buffer; // alle Knoten beschriften; Fallunterscheidung bool inactive= true; - list<c_BranchAndBoundNode*>::const_iterator act; - for ( act=active_nodes.begin(); act!=active_nodes.end(); ++act ) - if ( *node == *act ) + for ( auto act : active_nodes ) + if ( node == act ) { inactive = false; break; } - // do not display intermediate nodes - if ( InfoDottySuppressIntermediateNodes() && (*node)->NumSons() == 1 ) - continue; - if ( (*node)->NodeNumber() < 0 ) - continue; - // determine (proper) father - c_BranchAndBoundNode* son = *node; - c_BranchAndBoundNode* father = (*node)->Father(); - // do not display intermediate nodes - while( InfoDottySuppressIntermediateNodes() && father && father->Level() == (*node)->Level() ) - { - son = father; - father = father->Father(); - } - - if ( (*node) == selected_node ) - { - // current node - buffer << "node" << (*node)->NodeNumber() - << " [color=green, shape=record, label= \"{num=" << (*node)->NodeNumber() - << "|}|{LB=" << (*node)->LB() - << "|*current*}\"];"; - } - else + // do not display intermediate nodes + if ( InfoDottySuppressIntermediateNodes() && node->NumSons() == 1 ) + continue; + if ( node->NodeNumber() < 0 ) + continue; + // determine (proper) father + c_BranchAndBoundNode* son = node; + c_BranchAndBoundNode* father = node->Father(); + // do not display intermediate nodes + while( InfoDottySuppressIntermediateNodes() && father && father->Level() == node->Level() ) + { + son = father; + father = father->Father(); + } + if ( node == selected_node ) + { + // current node + buffer << "node" << node->NodeNumber() + << " [color=green, shape=record, label= \"{num=" << node->NodeNumber() + << "|}|{LB=" << node->LB() + << "|*current*}\"];"; + } + else + { + // already solved or terminated? + if ( inactive ) { - // already solved or terminated? - if ( inactive ) + string text = ""; + string color = "black"; + ostringstream sol_text2; + sol_text2 << "sol=" << node->NodeNumberSolved(); + string sol_text = sol_text2.str(); + if ( node->Terminated() ) { - string text = ""; - string color = "black"; - ostringstream sol_text2; - sol_text2 << "sol=" << (*node)->NodeNumberSolved(); - string sol_text = sol_text2.str(); - if ( (*node)->Terminated() ) - { - text = "terminated"; - sol_text = ""; - } - else + text = "terminated"; + sol_text = ""; + } + else + { + if ( node->IsFeasible() ) { - if ( (*node)->IsFeasible() ) + if (node->IsIntegral()) { - if ( (*node)->IsIntegral() ) - if ( (*node)->LB() == UB() ) - { - text = "best integer"; - color = "blue"; - } - else - { - text = "integer"; - color = "lightblue"; - } - } - else - { - text = "infeasible"; - color = "yellow"; + if (node->LB() == UB()) + { + text = "best integer"; + color = "blue"; + } + else + { + text = "integer"; + color = "lightblue"; + } } } - - if ( (*node)->LB() > son->LB() + CG_EPS ) - buffer << "node" << (*node)->NodeNumber() - << " [color=" << color - << ", shape=record, label= \"{num="<< (*node)->NodeNumber() - << "|" << sol_text - << "}|{LB=" << son->LB() - << "\\n" << (*node)->LB() - << "|" << text - << "}\"];"; else - buffer << "node" << (*node)->NodeNumber() - << " [color=" << color - << ", shape=record, label= \"{num=" << (*node)->NodeNumber() - << "|" << sol_text - << "}|{LB=" << (*node)->LB() - << "|" << text - << "}\"];"; + { + text = "infeasible"; + color = "yellow"; + } } - // active node - if ( !inactive ) - buffer << "node" << (*node)->NodeNumber() << " [color=red, shape=circle, label=\"active\"];"; + + if (node->LB() > son->LB() + LocalCG_EPS()) + buffer << "node" << node->NodeNumber() + << " [color=" << color + << ", shape=record, label= \"{num="<< node->NodeNumber() + << "|" << sol_text + << "}|{LB=" << son->LB() + << "\\n" << node->LB() + << "|" << text + << "}\"];"; + else + buffer << "node" << node->NodeNumber() + << " [color=" << color + << ", shape=record, label= \"{num=" << node->NodeNumber() + << "|" << sol_text + << "}|{LB=" << node->LB() + << "|" << text + << "}\"];"; } - bbout << buffer.str() << endl; + // active node + if ( !inactive ) + buffer << "node" << node->NodeNumber() << " [color=red, shape=circle, label=\"active\"];"; + } + bbout << buffer.str() << endl; - // arcs - // connection to father node - if ( father ) + // arcs + // connection to father node + if ( father ) + { + ostringstream buffer; + buffer << "node" << father->NodeNumber() + << " -> node" << node->NodeNumber() << "[label=\""; + bbout << buffer.str(); + bool first = true; + for ( auto con : son->AddedConstraints()) { - ostringstream buffer; - buffer << "node" << father->NodeNumber() - << " -> node" << (*node)->NodeNumber() << "[label=\""; - bbout << buffer.str(); - const list<c_BranchAndBoundConstraint*>& constraints = son->AddedConstraints(); - list<c_BranchAndBoundConstraint*>::const_iterator con; - bool first = true; - for ( con=constraints.begin(); con !=constraints.end(); ++con ) - { - if ( !first ) - bbout << "\\n"; - bbout << *(*con); - first = false; - } - bbout << "\"]; " << endl; + if ( !first ) + bbout << "\\n"; + bbout << *con; + first = false; } + bbout << "\"]; " << endl; + } } - bbout << "}" << endl; } } - void c_Controller::SetBestSolution( c_Solution* sol ) { // Store Solution delete p_best_solution; p_best_solution = sol; - if ( sol->Objective() < UB() ) - SetUB( sol->Objective() ); + if (sol->Objective() < UB()) + { + SetUB(sol->Objective()); + //new set first int + if (info_time_first_integer_solution > timer.Seconds()) + { + info_time_first_integer_solution = timer.Seconds(); + first_integer_solution = sol->Objective(); + } + //end new first int + } // Best Integer Solution // and First Integer Solution info_time_best_integer_solution = timer.Seconds(); } - void c_Controller::Go( void ) { // initialization - timer.Start(); - int status = 0; - c_TimeInfo gc_timer; - gc_timer.Start(); - + //timer.Start(); + DeleteTree(); c_TimeInfo timer_init; timer_init.Restart(); - // creation of pricers and separators i_phase_in_bap = initialization; if ( pricing_problem_hierarchy.empty() ) AddPricingProblems(); if ( separation_problems.empty() ) AddSeparationProblems(); - // create RMP if ( !p_rmp ) { @@ -627,9 +655,9 @@ void c_Controller::Go( void ) p_rmp->Reset(); //CPXENVptr env = (CPXENVptr) this->CPLEXEnv()->getEnv(); - ////useful debugging parameter to check input to cplex more carefully and show problems in the structure of the model immediately, leads to earlier CplexExceptions + //useful debugging parameter to check input to cplex more carefully and show problems in the structure of the model immediately, leads to earlier CplexExceptions //CPXsetintparam(env, CPX_PARAM_DATACHECK, CPX_ON); - ////prints cplex logging in the console + //prints cplex logging in the console //CPXsetintparam(env, CPX_PARAM_SCRIND, CPX_ON); ////to write cplex messages in file "CplexOutput.txt" ////get pointers of type CPXCHANNELptr to hold the address of the channel corresponding cpxresults, cpxwarning, cpxerror and cpxlog. May be NULL. @@ -656,7 +684,6 @@ void c_Controller::Go( void ) p_red_costs_fixing = CreateReducedCostsFixing(); // create root node of b&b-tree - DeleteTree(); i_phase_in_bap = root_node; p_root_node = CreateRootNode(); AddNode( p_root_node ); @@ -667,36 +694,28 @@ void c_Controller::Go( void ) while ( ( akt_node = SelectNextNode() ) && !TimeOut() ) { + p_currentNode = akt_node; + if ( akt_node->Father() ) i_phase_in_bap = tree; if ( ReloadSettings() ) ReadSettings(); - // garbage collection - gc_timer.Stop(); - if(gc_timer.Minutes() > 10) - gc_timer.Reset(); - gc_timer.Restart(); - c_TimeInfo timer_Node; timer_Node.Restart(); // check if node can be terminated - if ( p_rmp->RoundLBTo( akt_node->LB() ) >= ub - CG_EPS - || ( p_root_node->SubtreeLB() >= AspirationLB() - CG_EPS ) - || ( UB() < AspirationUB() + CG_EPS ) ) + if (p_rmp->RoundLBTo(akt_node->LB()) >= ub - LocalCG_EPS() + || (p_root_node->SubtreeLB() >= AspirationLB() - LocalCG_EPS()) + || (UB() < AspirationUB() + LocalCG_EPS())) { // yes, terminate... - if ( InfoLevel() >= 2 ) + if ( InfoLevel() >= 2 ) { - if ( akt_node->LB() >= ub ){ - Info2() << "\nNew Node: Termination (LB[" - << akt_node->LB() - << "] >= UB[" << ub << "]).\n" << flush; - } else{ - Info2() << "\nNew Node: Termination because of aspiration level\n" - << flush; - } + if ( akt_node->LB() >= ub ) + Info2() << "\nNew Node: Termination (LB[" << akt_node->LB() << "] >= UB[" << ub << "]).\n" << flush; + else + Info2() << "\nNew Node: Termination because of aspiration level\n" << flush; } akt_node->Terminate(); timer_Node.Stop(); @@ -712,12 +731,10 @@ void c_Controller::Go( void ) { Info2() << "\nLB = " << p_root_node->SubtreeLB() << ", UB = " << ub - << ", open nodes: " << (int)active_nodes.size() << "\n" << flush; + << ", open nodes: " << NumOpenBranchAndBoundNodes() << "\n" << flush; } - // here comes the code for solving c_DualVariables* p_dual = akt_node->Go(); - node_time.Stop(); info_timer_nodeSolve += node_time.Seconds(); @@ -735,6 +752,7 @@ void c_Controller::Go( void ) c_TimeInfo timer_updateUB; timer_updateUB.Restart(); bool is_integral = akt_node->IsIntegral(); + bool branching_or_cutting = false; if ( is_integral ) { // solution is integral @@ -772,7 +790,7 @@ void c_Controller::Go( void ) else { // solution is not integral - if ( ( akt_node->LB() < ub - CG_EPS ) + if ((p_rmp->RoundLBTo(akt_node->LB()) < ub - LocalCG_EPS()) && ( akt_node->IsFeasible() ) ) { if ( !TimeOut() ) @@ -782,9 +800,11 @@ void c_Controller::Go( void ) ReducedCostsFixing()->ComputeLowerBounds( p_rmp, p_dual, akt_node->RMPValue(), true ); // branching - if ( Branching() || CuttingPlanes() ) - PerformBranching( akt_node ); - + if (Branching() || CuttingPlanes()) + { + PerformBranching(akt_node); + branching_or_cutting = true; + } last_node_terminated = false; } } @@ -797,10 +817,11 @@ void c_Controller::Go( void ) } last_node_integer = false; - // Heuristik, try to solve as MIP - if ( UseMIPSolverForUBs() - && akt_node->DetailedLevel() <= UseMIPSolverUptoLevel() - && !TimeOut() ) + // Heuristic, try to solve as MIP + if ( ( UseMIPSolverForUBs() + && akt_node->Level() <= UseMIPSolverUptoLevel() && akt_node->Level() >= UseMIPSolverUpFromLevel() + && !TimeOut() && akt_node->IsFeasible()) + ) { akt_node->SolveAsMIP( MIPMaxSolutionTimeSec() ); if ( ReducedCostsFixing() ) @@ -829,9 +850,9 @@ void c_Controller::Go( void ) InfoStatistics() << "Node " << setw(5) << akt_node->NodeNumber() << "(" << setw(2) << akt_node->Level() << ")\t"; InfoStatistics() << "(" - << setw(5) << (int)active_nodes.size() << "/" + << setw(5) << NumOpenBranchAndBoundNodes() << "/" << setw(5) << (int)inactive_nodes.size() << "/" - << setw(5) << (int)(active_nodes.size()+this->inactive_nodes.size() ) << ")\t"; + << setw(5) << NumBranchAndBoundNodes() << ")\t"; InfoStatistics() << setw(8) << p_root_node->SubtreeLB() << setw(8) << ub @@ -839,7 +860,7 @@ void c_Controller::Go( void ) if ( akt_node->IsIntegral() ) { - InfoStatistics() << " int"; + InfoStatistics() << " int"; if ( ub == akt_node->LB() ) { InfoStatistics() << ": currently best solution "; @@ -847,17 +868,47 @@ void c_Controller::Go( void ) InfoStatistics() << endl << *BestSolution() << flush; } } - if ( !akt_node->IsFeasible() ) - InfoStatistics() << " infeasible"; - if ( !akt_node->Sons().empty() ) + else + { + if (!akt_node->IsFeasible()) + InfoStatistics() << " infeasible"; + else + { + if (branching_or_cutting) + { + if (akt_node->NumSons() == 1 + && !(akt_node->Sons().front()->AddedConstraints().size() ==1 && akt_node->Sons().front()->AddedConstraints().front()->IsStrongBranchingCandidates())) + { + InfoStatistics() << " cutting: " << akt_node->Sons().front()->AddedConstraints().front()->TextInBranching(); + } + else + { + InfoStatistics() << " branching: "; + if (akt_node->NumSons() > 0) + InfoStatistics() << akt_node->Sons().front()->AddedConstraints().front()->TextInBranching(); + else + InfoStatistics() << "strong br"; + } + } + else + InfoStatistics() << " bounding"; + } + } + /*if ( !akt_node->Sons().empty() ) InfoStatistics() << " #" << akt_node->Sons().size() - << ":" << akt_node->Sons().front()->AddedConstraints().front()->TextInBranching(); + << ":" << akt_node->Sons().front()->AddedConstraints().front()->TextInBranching();*/ + InfoStatistics() << endl << flush; } delete p_dual; } } // post processing + if ((TimeOut() && SolveAsMipAfterTimeOut())) + { + p_rmp->Update(p_root_node); + p_root_node->SolveAsMIP(MIPMaxSolutionTimeSec()); + } timer.Stop(); UpdateAllNodeInformations( NULL ); if ( InfoLevel() >= 1 ) @@ -887,17 +938,16 @@ void c_Controller::Go( void ) StatisticsBranching( counter_for_branching ); map<string,int>::const_iterator ii; int num_branchings = 0; - for ( ii=counter_for_branching.begin(); ii!=counter_for_branching.end(); ++ii ) - num_branchings += ii->second; - for ( ii=counter_for_branching.begin(); ii!=counter_for_branching.end(); ++ii ) + for ( auto ii : counter_for_branching ) + num_branchings += ii.second; + for ( auto ii : counter_for_branching ) Info1() << " | " - << setw(5) << ii->second << " branchings on " << setw(10) << ii->first << " (" - << int( double(ii->second)/double(num_branchings)*100.0 +.5 ) + << setw(5) << ii.second << " branchings on " << setw(10) << ii.first << " (" + << int( double(ii.second)/double(num_branchings)*100.0 +.5 ) << "%)" << endl; // cutting statistics - list<c_SeparationProblem*>::const_iterator prob; - for ( prob=SeparationProblems().begin(); prob!=SeparationProblems().end(); ++prob ) - (*prob)->OutputInOStream(cout); + for ( auto prob : SeparationProblems() ) + prob->OutputInOStream(cout); // overall best solution if ( p_best_solution && InfoPrintOverallBestSolution() ) Info1() << "Best integer solution: \n" @@ -916,22 +966,22 @@ void c_Controller::PrintColumnStatistics() { // compute aggregated statistics (like in pivot tables) map<list<int>,int> agg_statistics_columns; - for ( auto ii=o_info_statistics_columns.begin(); ii!=o_info_statistics_columns.end(); ++ii ) + for ( auto ii : o_info_statistics_columns ) { - list<int> col_info = ii->first; + list<int> col_info = ii.first; while ( col_info.size() > 2 ) { col_info.pop_back(); - agg_statistics_columns[ col_info ] += ii->second; + agg_statistics_columns[ col_info ] += ii.second; } } - for ( auto ii=agg_statistics_columns.begin(); ii!=agg_statistics_columns.end(); ++ii ) - o_info_statistics_columns[ ii->first ] += ii->second; + for ( auto ii : agg_statistics_columns ) + o_info_statistics_columns[ ii.first ] += ii.second; // output if ( InfoLevel() >= 1 ) { Info1() << "Columns Statistics: " << endl; - for ( int i=0; i<3; i++ ) // init, root, tree + for (int i : {0, 1, 2}) // init, root, tree { if ( i==0 ) Info1() << " initialization: " << endl; @@ -939,15 +989,15 @@ void c_Controller::PrintColumnStatistics() Info1() << " linear relax.: " << endl; else if ( i==2 ) Info1() << " search tree: " << endl; - for ( auto ii=o_info_statistics_columns.begin(); ii!=o_info_statistics_columns.end(); ++ii ) - if ( ii->first.front() == i ) + for ( auto ii : o_info_statistics_columns ) + if ( ii.first.front() == i ) { Info1() << " "; - list<int> col_info = ii->first; + list<int> col_info = ii.first; Info1() << "#"; - for ( auto jj=col_info.begin(); jj!=col_info.end(); ++jj ) - Info1() << *jj << "."; + for ( auto jj : col_info ) + Info1() << jj << "."; col_info.pop_front(); int col_type = col_info.empty() ? -1 : col_info.front(); @@ -956,31 +1006,35 @@ void c_Controller::PrintColumnStatistics() Info1() << s_info_column_type_names[col_type]; else Info1() << col_type; - for ( auto jj=col_info.begin(); jj!=col_info.end(); ++jj ) - Info1() << "~" << *jj; - Info1() << ": " << ii->second << endl; + for ( auto jj : col_info ) + Info1() << "~" << jj; + Info1() << ": " << ii.second << endl; } } } } } +bool pairsorterNodes(const pair<double, list<c_BranchAndBoundNode*> >& i, const pair<double, list<c_BranchAndBoundNode*> >& j) +{ + return i.first < j.first; +} void c_Controller::PerformBranching( c_BranchAndBoundNode* akt_node ) { - // Was Branching performed earlier? + GetBranchingInfo.Reset(); + GetBranchingInfo.Start(); + // was branching performed earlier? if ( akt_node->NumSons() ) return; - list< list<c_BranchAndBoundConstraint*> > cons_list; - - // Sind Schnittebenen erlaubt und werden auch ggf. welche separiert? + // cutting planes if ( CuttingPlanes() ) { list<c_BranchAndBoundConstraint*> cons; - akt_node->GetSeparationConstraints( cons ); - if ( !cons.empty() ) - cons_list.push_back( cons ); + akt_node->GetSeparationConstraints(cons); + if (!cons.empty()) + cons_list.push_back(cons); } // time for root node if ( !akt_node->Father() ) @@ -993,139 +1047,30 @@ void c_Controller::PerformBranching( c_BranchAndBoundNode* akt_node ) } if ( Branching() && cons_list.empty() ) { - // ********** strong branching *********** // - int max_k = StrongBranchingNumCandidates(); - // check whether UB is a good estimate - if ( UB() / p_root_node->LB() < 1.2 ) // 20% gap + if ( StrongBranching() ) { - double exp_alpha = 1.5; - double alpha = pow( ( UB() - akt_node->LB() )/( UB() - p_root_node->LB() ), exp_alpha ); - // alpha == 0 <=> close to opt - // alpha == 1 <=> at root - max_k = (int) floor( 0.5 + alpha * StrongBranchingNumCandidates() + (1-alpha) * StrongBranchingMinNumCandidates() ); - } - else - max_k = max( StrongBranchingMinNumCandidates(), StrongBranchingNumCandidates() - int( StrongBranchingNumCandidatesDecreasePerLevel() * akt_node->Level()) ); - - if ( StrongBranching() && akt_node->Level() < StrongBranchingUptoLevel() && max_k > 1 ) - { - if ( InfoLevel() >=2 ) - Info2() << "Strong Branching..." << endl; - - GetBranchingInfo.Restart(); - vector< pair<double, list<c_BranchAndBoundNode*> > > nodes; - for ( int k=0; k<max_k; k++ ) + // strong branching + int max_k = StrongBranchingNumCandidates(); + max_k = max(StrongBranchingMinNumCandidates(), StrongBranchingNumCandidates() - int(StrongBranchingNumCandidatesDecreasePerLevel() * akt_node->Level())); + if (max_k > 1) { - list< list<c_BranchAndBoundConstraint*> > local_cons_list; - - akt_node->Get_k_th_best_BranchAndBoundContraints( k, local_cons_list ); - if ( k==0 ) - cons_list = local_cons_list; - if ( local_cons_list.size() == 0 ) - { - max_k = k; - break; - } - nodes.push_back( pair<double, list<c_BranchAndBoundNode*> >() ); - - list< list<c_BranchAndBoundConstraint*> >::iterator cons; - for ( cons=local_cons_list.begin(); cons!=local_cons_list.end(); ++cons ) - { - // create node - c_BranchAndBoundNode* new_node = akt_node->CreateNewNode( *cons ); - nodes[k].second.push_back( new_node ); - } - if ( local_cons_list.size() == 1 ) - { - max_k = k+1; - break; - } - + c_StrongBranchingCandidates* candidates = akt_node->ComputeStrongBranchingCandidates(max_k); + cons_list.push_back({ candidates }); } - double best_lb_of_test_sons = akt_node->LB(); - list<c_BranchAndBoundNode*>::iterator node; - for ( int k=0; k<max_k; k++ ) - { - double min_lb_of_test_sons = INFTY; - for ( node=nodes[k].second.begin(); node!=nodes[k].second.end(); ++node ) - { - list<c_BranchAndBoundConstraint*>::const_iterator ii; - if ( InfoLevel() >=2 ) - { - Info2() << "Probing { "; - for ( ii=(*node)->AddedConstraints().begin(); ii!=(*node)->AddedConstraints().end(); ++ii ) - Info2() << *(*ii) << "; "; - Info2() << "} "; - } - if( StrongBranchingSolveExact() ){ - // exact - c_DualVariables* dual = (*node)->Go( numeric_limits<double>::infinity(), best_lb_of_test_sons ); - delete dual; - (*node)->SetUnsolved(); - if ( (*node)->LB() < min_lb_of_test_sons ) - min_lb_of_test_sons = (*node)->LB(); - } else{ - // heuristic - Update( *node ); - RMP()->Optimize( c_RMP::ePRIMOPT ); - //cout << " :: " << RMP()->ActualObjectiveValue(); - if ( RMP()->ActualObjectiveValue() < min_lb_of_test_sons ) - min_lb_of_test_sons = RMP()->ActualObjectiveValue(); - } - - if ( min_lb_of_test_sons < best_lb_of_test_sons ) - break; - } - nodes[k].first = -min_lb_of_test_sons; - - if ( InfoLevel() >=2 ) - Info2() << "Pair gives " << min_lb_of_test_sons; - if ( best_lb_of_test_sons < min_lb_of_test_sons ) - { - best_lb_of_test_sons = min_lb_of_test_sons; - if ( InfoLevel() >=2 ) - Info2() << " *(curr best)"; - } - if ( InfoLevel() >=2 ) - Info2() << endl; - } - - // According to improvement - sort( nodes.begin(), nodes.end() ); - - // undo changes; very important for a following MIP run - Update( akt_node ); - - for ( node=nodes[0].second.begin(); node!=nodes[0].second.end(); ++node ) - { - akt_node->AddSon( *node ); - AddNode( *node ); - } - if ( max_k > 1 ) - { - const int kk=0; - for ( int k=kk+1; k<max_k; k++ ) - for ( node=nodes[k].second.begin(); node!=nodes[k].second.end(); ++node ) - delete *node; - } - GetBranchingInfo.Stop(); - info_time_root_branching += GetBranchingInfo.Seconds(); - return; + else // standard + akt_node->GetBranchAndBoundConstraints(cons_list); } else { - // ********** nicht strong branching *********** // - GetBranchingInfo.Restart(); + // standard branching akt_node->GetBranchAndBoundConstraints( cons_list ); - GetBranchingInfo.Stop(); - info_time_root_branching += GetBranchingInfo.Seconds(); } + GetBranchingInfo.Stop(); + info_time_root_branching += GetBranchingInfo.Seconds(); } - - list< list<c_BranchAndBoundConstraint*> >::iterator cons; - for ( cons=cons_list.begin(); cons!=cons_list.end(); ++cons ) + for ( auto cons : cons_list ) { - c_BranchAndBoundNode* new_node = akt_node->CreateNewNode( *cons ); + c_BranchAndBoundNode* new_node = akt_node->CreateNewNode( cons ); akt_node->AddSon( new_node ); AddNode( new_node ); } @@ -1133,76 +1078,97 @@ void c_Controller::PerformBranching( c_BranchAndBoundNode* akt_node ) c_BranchAndBoundNode* c_Controller::SelectNextNode() { - // any node available - if ( active_nodes.size() == 0 ) - return NULL; - c_BranchAndBoundNode* erg = NULL; e_branching branching = BranchingRule(); - switch ( branching ) + switch (branching) { - case depth_first: + case depth_first: { // depth-first erg = *active_nodes.rbegin(); list<c_BranchAndBoundNode*>::reverse_iterator brother = active_nodes.rbegin(); brother++; - while ( brother!=active_nodes.rend() - && (*brother)->Father() == erg->Father() ) + while (brother != active_nodes.rend() && (*brother)->Father() == erg->Father()) { erg = *brother; brother++; } break; } - case best_first: + case best_first: { // best first, return node with smallest LB // break ties by selecting nodes from the tree from left to right... - double best_lb = INFTY * INFTY; + double best_lb = LocalInfty() * LocalInfty(); + int best_lvl = LocalInfty() ; list<c_BranchAndBoundNode*>::iterator i; - for ( i=active_nodes.begin(); i!=active_nodes.end(); ++i ) - if ( (*i)->LB() < best_lb + CG_EPS*CG_EPS ) + for (i = active_nodes.begin(); i != active_nodes.end();) + { + if (RMP()->RoundLBTo((*i)->LB()) > UB() - (LocalCG_EPS() * LocalCG_EPS())) //pruning + { + inactive_nodes.push_back(*i); + auto new_i = i; + new_i++; + active_nodes.erase(i); + i = new_i; + } + else { - erg = *i; - best_lb = (*i)->LB(); + if ( (*i)->LB() < best_lb || ( (int) (100* (*i)->LB() ) == (int)(100* best_lb) && (*i)->GetLevel_new() < best_lvl ) ) + { + erg = *i; + best_lb = (*i)->LB(); + best_lvl = (*i)->GetLevel_new(); + } + + + /*if ((*i)->LB() < best_lb + (LocalCG_EPS() * LocalCG_EPS())) + { + erg = *i; + best_lb = (*i)->LB(); + }*/ + ++i; } - break; + } + break; } - case fifo: + case fifo: { erg = *active_nodes.begin(); break; } } // this should not happen... - if ( erg == NULL ) + if (erg == NULL) return NULL; else { - active_nodes.remove( erg ); - inactive_nodes.push_back( erg ); + // is this node a strong banching node with candidates not evaluated yet + if (erg->AddedConstraints().size() == 1 && erg->AddedConstraints().front()->IsStrongBranchingCandidates()) + { + erg->DecideOnStrongBranchingCandidate(); + } + // now we have a node to process + active_nodes.remove(erg); + inactive_nodes.push_back(erg); return erg; } } - void c_Controller::AddColumns( const list<c_Column*>& col_list ) { p_rmp->AddColumns( col_list ); // store in controller int position = (int)cols_in_rmp.size(); cols_in_rmp.resize( cols_in_rmp.size() + col_list.size() ); - list<c_Column*>::const_iterator col; - for ( col = col_list.begin(); col!=col_list.end(); ++col ) + for ( auto col : col_list ) { - cols_in_rmp[position] = *col; + cols_in_rmp[position] = col; cols_in_rmp[position]->SetNumInRMP( position ); position++; } } - void c_Controller::AddNode( c_BranchAndBoundNode* node ) { node->SetRemoved( false ); @@ -1211,32 +1177,26 @@ void c_Controller::AddNode( c_BranchAndBoundNode* node ) { inactive_nodes.push_back( node ); // recursively for the sons... - list<c_BranchAndBoundNode*>::const_iterator son_iter; - for ( son_iter = node->Sons().begin(); son_iter != node->Sons().end(); ++son_iter ) - AddNode( *son_iter ); - + for ( auto son_iter : node->Sons() ) + AddNode( son_iter ); return; } - if ( node->NodeNumber() < 0 ) node->SetNodeNumber( next_node_id++ ); active_nodes.push_back( node ); // globally valid constraints can be added directly to the RMP - const list<c_BranchAndBoundConstraint*>& cons = node->AddedConstraints(); - for ( auto bb_con=cons.begin(); bb_con!=cons.end(); ++bb_con ) + for ( auto bb_con : node->AddedConstraints()) { // loop over all rmp constraints for current B&B constraint - const list<c_Constraint*>& rmp_constraints = (*bb_con)->RMPConstraints(); list<c_Constraint*> constraintsToAdd; - for ( auto rmp_con=rmp_constraints.begin(); rmp_con!=rmp_constraints.end(); ++rmp_con ) - if ( (*rmp_con)->IsActiveAtNode( p_root_node ) && !p_rmp->IsConstraintInRMP( *rmp_con ) ) - constraintsToAdd.push_back(*rmp_con); + for ( auto rmp_con : bb_con->RMPConstraints()) + if ( rmp_con->IsActiveAtNode( p_root_node ) && !p_rmp->IsConstraintInRMP( rmp_con ) ) + constraintsToAdd.push_back(rmp_con); p_rmp->AddConstraints( constraintsToAdd ); } } - void c_Controller::RemoveNode( c_BranchAndBoundNode* node ) { node->SetRemoved( true ); @@ -1255,19 +1215,16 @@ void c_Controller::RemoveNode( c_BranchAndBoundNode* node ) inactive_nodes.erase( i ); // recursively for the sons... - list<c_BranchAndBoundNode*>::const_iterator son_iter; - for ( son_iter = node->Sons().begin(); son_iter != node->Sons().end(); ++son_iter ) - RemoveNode( *son_iter ); - + for ( auto son_iter : node->Sons() ) + RemoveNode( son_iter ); if ( InfoLevel() >=3 ) Info3() << "Removed node " << node->NodeNumber() << endl; } - void c_Controller::Update( c_BranchAndBoundNode* node ) { c_TimeInfo timer_RMPUpdate; - timer_RMPUpdate.Restart(); + timer_RMPUpdate.Start(); // Update of RMP RMP()->Update( node ); timer_RMPUpdate.Stop(); @@ -1277,23 +1234,13 @@ void c_Controller::Update( c_BranchAndBoundNode* node ) timer_PricingUpdate.Restart(); // loop over all Pricing Problems and update vector<c_PricingProblem*>::iterator prob; - for( int h = 0; h < pricing_problem_hierarchy.size(); h++){ - pricing_problem_hierarchy[h]->Update( node ); - pricing_levels_in_hierarchies[h] = 0; - } - ResetNumFailures(); - - /*vector<c_PricingProblemHierarchy*>::iterator hierarchy; - for ( hierarchy = pricing_problem_hierarchy.begin(); hierarchy != pricing_problem_hierarchy.end(); ++hierarchy ){ - (*hierarchy)->Update( node ); - }*/ + for ( auto hierarchy : pricing_problem_hierarchy ) + hierarchy->Update( node ); timer_PricingUpdate.Stop(); info_timer_PricingUpdate += timer_PricingUpdate.Seconds(); } - - -void c_Controller::Pricing( int iteration, c_DualVariables* dual, list<c_Column*>& col_list ) +void c_Controller::Pricing( int iteration, c_DualVariables* dual, list<c_Column*>& col_list, int NumPricersUsed ) { // Initialisierung list<c_ColumnPricePair> column_price_list; @@ -1303,125 +1250,15 @@ void c_Controller::Pricing( int iteration, c_DualVariables* dual, list<c_Column* Info3() << "Pricing " << flush; c_TimeInfo pricing_time_info; + pricing_time_info.Start(); int generated_cols = 0; - vector<c_PricingProblemHierarchy*>& hierarchies = PricingProblemHierarchies(); - - if( PricingHierarchyStandardOrder() ){ - vector<c_PricingProblemHierarchy*>::iterator hh; - for ( hh=hierarchies.begin(); hh!=hierarchies.end(); ++hh ) - { - if( TimeOut() ) - return; - c_PricingProblemHierarchy* hierarchy = *hh; - hierarchy->Update( dual ); - generated_cols += hierarchy->Go( dual, column_price_list ); - //info_number_of_solved_pricing_problems++; - } - } else{ - //NowTODO - int startHierarchy = GetStartHierarchy(); - int numHierarchies = (int)hierarchies.size(); - int maxPricingLevel = hierarchies[0]->GetNumPricingLevels()-1; - for ( int h = 0; h < numHierarchies; h++ ) { - hierarchies[h]->SetBestReducedCost(-10*INFTY); //Initialize reduced cost to -infinity, so that the lagrangean lower bound is also valid, if not all pricing hierarchies are started - } - - bool continuePricing = true; - vector<bool> v_failed = vector<bool>(numHierarchies, false); - //vector<bool> v_failed = vector<bool>(numHierarchies, false); - do{ - if( TimeOut() ) - return; - - c_PricingProblemHierarchy* hierarchy = hierarchies[startHierarchy]; - hierarchy->Update( dual ); - - if( all_of(v_failed.begin(), v_failed.end(), [](bool b){return b;}) ){ - IncreaseNumFailures( min(hierarchy->GetPricingStartLevel(), maxPricingLevel) ); - if( hierarchy->GetPricingStartLevel() >= maxPricingLevel ) //no more improving columns for this node - break; //Other behaviour when MinNumCols > 0? - else{ - //increase pricing level of all hierarchies (at least temporarily; if maxNumFailures is not reached, it is decreased again) - int minPricingStartLevel = maxPricingLevel; - for( int p = 0; p < numHierarchies; p++){ - if ( hierarchies[p]->GetPricingStartLevel() < minPricingStartLevel ) - minPricingStartLevel = hierarchies[p]->GetPricingStartLevel(); - } - for( int p = 0; p < v_failed.size(); p++){ - hierarchies[p]->SetPricingStartLevel( minPricingStartLevel + 1 ); //hierarchies[p]->GetPricingStartLevel() +1 ); - v_failed[p] = false; - } - } - } - - int added_cols = hierarchy->Go( dual, column_price_list ); - generated_cols += added_cols; - - if( added_cols >= hierarchy->GetMinNumCols() ){ //stay at the same hierarchy until no more columns for it can be found - continuePricing = false; - SetStartHierarchy(startHierarchy); - } else{ - v_failed[startHierarchy] = true; - startHierarchy = (startHierarchy + 1)%numHierarchies; - } - - } while( continuePricing ); - - ////New test: always start with pricing hierarchy which has the minimal pricingStartLevel - //int maxPricingLevel = hierarchies[0]->GetNumPricingLevels(); - //int minPricingStartLevel = maxPricingLevel+1; - //int bestH = -1; - //int numHierarchies = (int)hierarchies.size(); - //for ( int h = 0; h < numHierarchies; h++ ) { - // hierarchies[h]->SetBestReducedCost(-10*INFTY); //Initialize reduced cost to -infinity, so that the lagrangean lower bound is also valid, if not all pricing hierarchies are started - // if( pricing_levels_in_hierarchies[h] < minPricingStartLevel ){ - // minPricingStartLevel = hierarchies[h]->GetPricingStartLevel(); - // bestH = h; - // } - //} - - ////stay at one pricing hierarchy until no more columns can be found. Then switch first to other hierarchies with same pricing level. - //int hNum = bestH; - //bool continuePricing = true; - //vector<bool> hierarchyRunWithMaxPricingLevel = vector<bool>(numHierarchies); - //for( int p = 0; p < numHierarchies; p++) - // hierarchyRunWithMaxPricingLevel[p] = false; - //bool pricingLevelIncreased = false; - //do{ - // if( TimeOut() ) - // return; - // c_PricingProblemHierarchy* hierarchy = hierarchies[hNum]; - // hierarchy->Update( dual ); - // hierarchyRunWithMaxPricingLevel[hNum] = (pricing_levels_in_hierarchies[hNum] >= maxPricingLevel-1)? true : false; - // generated_cols += hierarchy->Go( dual, column_price_list ); - // int curPricingLevel = hierarchy->GetPricingStartLevel(); - // if( curPricingLevel > pricing_levels_in_hierarchies[hNum] ) - // pricingLevelIncreased = true; - // if( generated_cols > hierarchies[hNum]->GetMinNumCols() ){ //stay at the same hierarchy until no more columns for it can be found - // //if the maximal number of failures is not reached and the pricing level below, go back - // if( curPricingLevel > 0 && hierarchy->GetNumFailures(curPricingLevel-1) < hierarchy->GetMaxNumFailures() ){ - // hierarchy->SetPricingStartLevel( curPricingLevel-1 ); - // } - // pricing_levels_in_hierarchies[hNum] = hierarchy->GetPricingStartLevel(); - // break; - // } - // pricing_levels_in_hierarchies[hNum] = curPricingLevel; - // hNum= hNum + 1; - // hNum = hNum % numHierarchies; - - // //if maxNumFailures is greater 1, ensure that the pricing level is although increased if no hierarchy found new columns and the same pricing problem would be resolved - // if( hNum == bestH && generated_cols == 0 && !pricingLevelIncreased){ //GetMinPricingStartLevel() != maxPricingLevel - // //all hierarchies were tested and found no column, but maximal pricing level not reached, such that pricing will continue: - // //increase all minimal pricing levels (by ignoring a maxNumFailures > 1), since nothing has changed - // for ( int h = 0; h < numHierarchies; h++ ) { - // if( hierarchies[h]->GetPricingStartLevel() == GetMinPricingStartLevel() ){ - // hierarchies[h]->SetPricingStartLevel(pricing_levels_in_hierarchies[h]+1); - // pricing_levels_in_hierarchies[h] += 1; - // } - // } - // } - // continuePricing = ( find(begin(hierarchyRunWithMaxPricingLevel), end(hierarchyRunWithMaxPricingLevel), false) != end(hierarchyRunWithMaxPricingLevel) ); // not all true - //} while( continuePricing ); //hNum != bestH || GetMinPricingStartLevel() != maxPricingLevel + for ( auto hh : PricingProblemHierarchies() ) + { + if( TimeOut() ) + return; + hh->Update( dual ); + generated_cols += hh->Go( dual, column_price_list, NumPricersUsed ); + info_number_of_solved_pricing_problems++; } // Select an appropriate subset of columns to include into the @@ -1431,9 +1268,8 @@ void c_Controller::Pricing( int iteration, c_DualVariables* dual, list<c_Column* if ( PricingWithColumnSelection() ) { ColumnSelectionInPricing( column_price_list, to_remove_list ); - list<c_ColumnPricePair>::iterator col_price; - for ( col_price=to_remove_list.begin(); col_price!=to_remove_list.end(); ++col_price) - delete (*col_price).first; + for ( auto col_price : to_remove_list) + delete col_price.first; } pricing_time_info.Stop(); @@ -1451,22 +1287,22 @@ void c_Controller::Pricing( int iteration, c_DualVariables* dual, list<c_Column* Info3() << " (" << pricing_time_info.MilliSeconds() << "ms)" << flush; // check, if too many columns generated - double threshold = INFTY * INFTY; + int max_cols_to_add = MaxColumnsToAddInPricing(); + double threshold = LocalInfty() * LocalInfty(); if ( (int)column_price_list.size() >= MaxColumnsToAddInPricing() ) { // take only the best MaxColumnsToAddInPricing() of c_FindKthElement k_of_n_alg( (int)column_price_list.size() ); - list<c_ColumnPricePair>::iterator col_price; int idx=0; - for ( col_price=column_price_list.begin(); col_price!=column_price_list.end(); ++col_price ) - k_of_n_alg.SetElement( idx++, (*col_price).second ); - threshold = k_of_n_alg.FindKthElementOfN( MaxColumnsToAddInPricing() - 1, (int)column_price_list.size() ); + for ( auto col_price : column_price_list ) + k_of_n_alg.SetElement( idx++, col_price.second ); + threshold = k_of_n_alg.FindKthElementOfN(max_cols_to_add - 1, (int)column_price_list.size() ); if ( InfoLevel() >= 3 ) { int to_add=0; - for ( col_price=column_price_list.begin(); col_price!=column_price_list.end(); ++col_price ) - if ( (*col_price).second <= threshold ) + for ( auto col_price : column_price_list ) + if ( col_price.second <= threshold && to_add < max_cols_to_add ) to_add++; Info3() << ", only " << to_add << " added.\n"; } @@ -1477,21 +1313,17 @@ void c_Controller::Pricing( int iteration, c_DualVariables* dual, list<c_Column* // Transfer Columns from pairs to single columns col_list.clear(); - - list<c_ColumnPricePair>::iterator col_price; - for ( col_price= column_price_list.begin(); - col_price!=column_price_list.end(); - ++col_price ) - if ( (*col_price).second <= threshold ) - col_list.push_back( (*col_price).first ); + for ( auto col_price : column_price_list ) + if ( col_price.second <= threshold && (int)col_list.size() < max_cols_to_add ) + col_list.push_back( col_price.first ); else - delete (*col_price).first; + delete col_price.first; column_price_list.clear(); } void c_Controller::Separation( c_RMP* rmp, list<c_Constraint*>& cuts, int node_level ) { - const list<c_SeparationProblem*>& separators = SeparationProblems(); + auto& separators = SeparationProblems(); if ( separators.size() == 0 ) return; if ( InfoLevel() >= 3 && separators.size() ) @@ -1499,9 +1331,9 @@ void c_Controller::Separation( c_RMP* rmp, list<c_Constraint*>& cuts, int node_l // loop over all separation problems c_TimeInfo separation_time_info; - list<c_SeparationProblem*>::const_iterator prob; - for ( prob=separators.begin(); prob!=separators.end(); ++prob ) - (*prob)->Go( rmp, cuts, node_level ); + separation_time_info.Start(); + for ( auto prob : separators ) + prob->Go( rmp, cuts, node_level ); separation_time_info.Stop(); int num_of_cuts = (int)cuts.size(); @@ -1512,7 +1344,6 @@ void c_Controller::Separation( c_RMP* rmp, list<c_Constraint*>& cuts, int node_l Info3() << " (" << separation_time_info.MilliSeconds() << "ms)." << flush; } - void c_Controller::AddSeparationProblem( c_SeparationProblem* prob ) { separation_problems.push_back( prob ); @@ -1523,15 +1354,36 @@ void c_Controller::AddFeasibilityCuttingProblem( c_SeparationProblem* prob ) feasibility_cutting_problems.push_back( prob ); } +int c_Controller::NumBranchAndBoundNodes() const +{ + return NumOpenBranchAndBoundNodes() + (int)inactive_nodes.size(); +} + +int c_Controller::NumOpenBranchAndBoundNodes() const +{ + return (int)active_nodes.size() + NumOpenStrongBranchingNodes(); +} + +int c_Controller::NumOpenStrongBranchingNodes() const +{ + int num_strong_branching_nodes = 0; + for (auto node : active_nodes) + if (node->AddedConstraints().size() == 1 && node->AddedConstraints().front()->IsStrongBranchingCandidates()) + { + c_StrongBranchingCandidates* cand = (c_StrongBranchingCandidates*)node->AddedConstraints().front(); + num_strong_branching_nodes += (cand->EstimatedNumBranchAndBoundNodes() - 1); + } + return num_strong_branching_nodes; +} + bool c_Controller::FeasibilityCutting() { bool result = false; - list<c_SeparationProblem*>::iterator iter; - for ( iter=feasibility_cutting_problems.begin(); iter!=feasibility_cutting_problems.end(); ++iter ) + for ( auto iter : feasibility_cutting_problems ) { // call each cutting problem list<c_Constraint*> cuts; - (*iter)->Go( RMP(), cuts, -1 ); + iter->Go( RMP(), cuts, -1 ); if ( !cuts.empty() ) { RMP()->AddConstraints( cuts ); @@ -1545,47 +1397,40 @@ double c_Controller::LB() const { if ( p_root_node ) return p_root_node->SubtreeLB(); - return -INFTY; + return -LocalInfty(); } double c_Controller::LB1() const // lower bound at root node { // Get all nodes... list<c_BranchAndBoundNode*> all_nodes; - list<c_BranchAndBoundNode*>::const_iterator node; - for ( node=active_nodes.begin(); node!=active_nodes.end(); ++node ) - all_nodes.insert( all_nodes.end(), *node ); - for ( node=inactive_nodes.begin(); node!=inactive_nodes.end(); ++node ) - all_nodes.insert( all_nodes.end(), *node ); - - for ( node=all_nodes.begin(); node!=all_nodes.end(); ++node ) - if ( !(*node)->Father() ) - return (*node)->LB(); - return -INFTY; + for ( auto node : active_nodes ) + all_nodes.insert( all_nodes.end(), node ); + for ( auto node : inactive_nodes ) + all_nodes.insert( all_nodes.end(), node ); + for ( auto node : all_nodes ) + if ( !node->Father() ) + return node->LB(); + return -LocalInfty(); } double c_Controller::LB2() const // lower bound at root node + cutting planes { // Get all nodes... list<c_BranchAndBoundNode*> all_nodes; - list<c_BranchAndBoundNode*>::const_iterator node; - for ( node=active_nodes.begin(); node!=active_nodes.end(); ++node ) - all_nodes.insert( all_nodes.end(), *node ); - for ( node=inactive_nodes.begin(); node!=inactive_nodes.end(); ++node ) - all_nodes.insert( all_nodes.end(), *node ); - - double result = -INFTY; - for ( node=all_nodes.begin(); node!=all_nodes.end(); ++node ) - if ( (*node)->IsRootNode() && result < (*node)->LB() ) - result = (*node)->LB(); + for ( auto node : active_nodes ) + all_nodes.insert( all_nodes.end(), node ); + for ( auto node : inactive_nodes ) + all_nodes.insert( all_nodes.end(), node ); + double result = -LocalInfty(); + for ( auto node : all_nodes ) + if ( node->IsRootNode() && result < node->LB() ) + result = node->LB(); return result; } - void c_Controller::OrganizePool() { - // - c_TimeInfo pool_time_info; double max_ratio_cols_rows = MaxRatioColsToRows(); long max_cols_allowed = (long) ceil( 1.333 * max_ratio_cols_rows * double( p_rmp->NumRows() ) ); @@ -1601,6 +1446,8 @@ void c_Controller::OrganizePool() Info3() << "\nColumnPool: " << flush; Info3() << "sort... " << flush; } + c_TimeInfo pool_time_info; + pool_time_info.Start(); double* c_red = new double[n]; p_rmp->GetReducedCosts( c_red ); @@ -1613,56 +1460,57 @@ void c_Controller::OrganizePool() if ( ColumnInRMP(i)->IsTransferableToPool() ) k_of_n_alg.SetElement( i, c_red[i] ); else - k_of_n_alg.SetElement( i, -INFTY ); + k_of_n_alg.SetElement(i, -LocalInfty()); double threshold = k_of_n_alg.FindKthElementOfN( max_cols_allowed2 - 1, n ); - for ( i=0; i<n; i++ ) - if ( ColumnInRMP(i)->IsTransferableToPool() - && ( c_red[i] > threshold ) - && ( c_red[i] >= CG_EPS ) ) + for (i = 0; i < n; i++) + { + if (ColumnInRMP(i)->IsTransferableToPool() + && (c_red[i] > threshold) + && (c_red[i] >= LocalCG_EPS())) { // test value double val; - p_rmp->getX( &val, i, i ); - if ( val == 0.0 ) - add_to_pool.insert( add_to_pool.end(), ColumnInRMP(i) ); + p_rmp->getX(&val, i, i); + if (val == 0.0) + add_to_pool.insert(add_to_pool.end(), ColumnInRMP(i)); } - delete[] c_red; + + } + delete[] c_red; - if ( InfoLevel() >= 3 ) - Info3() << "transfer... " << flush; - int pool_plus = (int)add_to_pool.size(); - TransferColumnsToPool( add_to_pool ); - add_to_pool.clear(); // no longer valid + if (InfoLevel() >= 3) + Info3() << "transfer... " << flush; + int pool_plus = (int)add_to_pool.size(); + TransferColumnsToPool(add_to_pool); + add_to_pool.clear(); // no longer valid - if ( InfoLevel() >= 3 ) - Info3() << "reopt... " << flush; - p_rmp->Optimize( c_RMP::eDUALOPT ); + if (InfoLevel() >= 3) + Info3() << "reopt... " << flush; + p_rmp->Optimize(c_RMP::eDUALOPT); - // Test - if ( fabs( test_rmp_val - p_rmp->objval() ) > CG_EPS ) - { - InfoErr() << "#ERROR: reopt. error after pool transfer" << endl; - InfoErr() << "#ERROR: " << test_rmp_val << " <> " << p_rmp->objval() << endl; - } + // Test + if (fabs(test_rmp_val - p_rmp->objval()) > 100.0 * LocalCG_EPS()) + { + InfoErr() << "#ERROR: reopt. error after pool transfer" << endl; + InfoErr() << "#ERROR: " << test_rmp_val << " <> " << p_rmp->objval() << endl; + } - pool_time_info.Stop(); - if ( InfoLevel() >= 3 ) - { - Info3() << pool_plus - << " columns added to pool (" - << pool_time_info.MilliSeconds() << " ms).\n" << flush; - } + pool_time_info.Stop(); + if (InfoLevel() >= 3) + { + Info3() << pool_plus + << " columns added to pool (" + << pool_time_info.MilliSeconds() << " ms).\n" << flush; + } } } - void c_Controller::TransferColumnsToPool( const list<c_Column*>& col_list ) { - list<c_Column*>::const_iterator col; - for ( col=col_list.begin(); col!= col_list.end(); ++col ) - if ( (*col)->IsInRMP() && (*col)->IsTransferableToPool() ) + for ( auto col : col_list ) + if ( col->IsInRMP() && col->IsTransferableToPool() ) { - int pos = (*col)->NumInRMP(); + int pos = col->NumInRMP(); cols_in_rmp[ pos ] = NULL; } else @@ -1674,19 +1522,19 @@ void c_Controller::TransferColumnsToPool( const list<c_Column*>& col_list ) p_rmp->delcols( i, i ); int pos = 0; - for ( i=0; i<(int)cols_in_rmp.size(); i++ ) - if ( cols_in_rmp[i] ) + for (i = 0; i < (int)cols_in_rmp.size(); i++) + { + if (cols_in_rmp[i]) { cols_in_rmp[pos] = cols_in_rmp[i]; - cols_in_rmp[pos]->SetNumInRMP( pos ); + cols_in_rmp[pos]->SetNumInRMP(pos); pos++; } - cols_in_rmp.resize(pos); - - p_pool->AddColumns( col_list ); + } + cols_in_rmp.resize(pos); + p_pool->AddColumns(col_list); } - void c_Controller::OutputInOStream( ostream& s ) const { s << "Controller: \n" @@ -1698,18 +1546,15 @@ void c_Controller::StatisticsBranching( map<string,int>& counter_for_branching ) { counter_for_branching.clear(); list<c_BranchAndBoundNode*> all_nodes; - list<c_BranchAndBoundNode*>::const_iterator node; - for ( node=active_nodes.begin(); node!=active_nodes.end(); ++node ) - all_nodes.insert( all_nodes.end(), *node ); - for ( node=inactive_nodes.begin(); node!=inactive_nodes.end(); ++node ) - all_nodes.insert( all_nodes.end(), *node ); - int num_nodes_with_sons = 0; - for ( node=all_nodes.begin(); node!=all_nodes.end(); ++node ) - if ( !(*node)->Sons().empty() ) + for ( auto node : active_nodes ) + all_nodes.insert( all_nodes.end(), node ); + for ( auto node : inactive_nodes ) + all_nodes.insert( all_nodes.end(), node ); + for ( auto node : all_nodes ) + if ( !node->Sons().empty() ) { - string text = (*node)->Sons().front()->AddedConstraints().front()->TextInBranching(); + string text = node->Sons().front()->AddedConstraints().front()->TextInBranching(); counter_for_branching[text]++; - num_nodes_with_sons++; } } @@ -1753,6 +1598,18 @@ int c_Controller::StatisticsColumns( const list<int>& col_info ) const return 0; } +bool c_Controller::local_IsFractional(double x) +{ + return (((ceil(x) - x > x - floor(x)) ? x - floor(x) : ceil(x) - x) > LocalCG_EPS()); +} + +c_Solution* c_Controller::CreateSolutionFromMIP(c_RMP* theRMP, vector<pair<int, double> >& columns_in_Mip_sol) +{ + cerr << "Serious Error! If you use SolveAsMip, you have to implement the function virtual c_Solution* CreateSolutionFromMIP(c_RMP* theRMP) in your problem-specific controller" << endl; + return nullptr; +} + + /**************************************************************************** c_ColumnPool @@ -1763,15 +1620,13 @@ c_ColumnPool::~c_ColumnPool() { } - void c_ColumnPool::AddColumns( const list<c_Column*>& col_list ) { bool do_nothing = true; if ( do_nothing ) { - list<c_Column*>::const_iterator col; - for ( col=col_list.begin(); col!=col_list.end(); ++col ) - delete (*col); + for ( auto col : col_list ) + delete col; return; } else @@ -1793,20 +1648,17 @@ void c_ColumnPool::AddColumns( const list<c_Column*>& col_list ) // in pool? // pairwise comparison - list<c_Column*>::const_iterator col1; - for ( col1=col_list.begin(); col1!=col_list.end(); ++col1 ) + for ( auto col1 : col_list ) { // nzcnt1 = (*col1)->CPXnzcnt(); - nzcnt1 = (*col1)->GetCPXColumn( cost1, ind1, val1, lb1, ub1, NULL ); + nzcnt1 = col1->GetCPXColumn( cost1, ind1, val1, lb1, ub1, NULL ); bool vorhanden = false; - - list<c_Column*>::iterator col2; - for ( col2=cols_in_pool.begin(); col2!=cols_in_pool.end(); ++col2 ) + for ( auto col2 : cols_in_pool ) { // Teste *col1 und *col2 auf Gleichheit bool equal = true; // nzcnt2 = (*col2)->CPXnzcnt(); - nzcnt2 = (*col2)->GetCPXColumn( cost2, ind2, val2, lb2, ub2, NULL ); + nzcnt2 = col2->GetCPXColumn( cost2, ind2, val2, lb2, ub2, NULL ); if ( nzcnt1!=nzcnt2 || cost1 != cost2 ) equal = false; else @@ -1824,9 +1676,8 @@ void c_ColumnPool::AddColumns( const list<c_Column*>& col_list ) } } if ( !vorhanden ) - cols_in_pool.insert( cols_in_pool.end(), *col1 ); + cols_in_pool.insert( cols_in_pool.end(), col1 ); } - // clean up everything delete[] ind1; delete[] ind2; @@ -1837,8 +1688,7 @@ void c_ColumnPool::AddColumns( const list<c_Column*>& col_list ) void c_ColumnPool::RemoveColumns( c_DualVariables* dual, double threshold ) { - list<c_Column*>::iterator col; - for ( col=cols_in_pool.begin(); col!=cols_in_pool.end(); ) + for ( auto col=cols_in_pool.begin(); col!=cols_in_pool.end(); ) if ( (*col)->ReducedCosts( *dual ) >= threshold ) { list<c_Column*>::iterator merk = col; @@ -1849,7 +1699,6 @@ void c_ColumnPool::RemoveColumns( c_DualVariables* dual, double threshold ) } else ++col; - } void c_ColumnPool::ReducePool( c_DualVariables* dual, int max_number_of_columns ) @@ -1857,34 +1706,25 @@ void c_ColumnPool::ReducePool( c_DualVariables* dual, int max_number_of_columns int size = (int)cols_in_pool.size(); if ( size <= max_number_of_columns ) return; - c_FindKthElement find_k_of_n( (int)cols_in_pool.size() ); - list<c_Column*>::iterator col; int pos = 0; - for ( col=cols_in_pool.begin(); col!=cols_in_pool.end(); ++col ) - find_k_of_n.SetElement( pos++, (*col)->ReducedCosts( *dual ) ); - - double threshold = find_k_of_n.FindKthElementOfN - ( max_number_of_columns, size, c_FindKthElement::increasing_order ); + for ( auto col : cols_in_pool ) + find_k_of_n.SetElement( pos++, col->ReducedCosts( *dual ) ); + double threshold = find_k_of_n.FindKthElementOfN( max_number_of_columns, size, c_FindKthElement::increasing_order ); RemoveColumns( dual, threshold ); } - - - void c_ColumnPool::MoveAllColumnsToRMP() { Controller()->AddColumns( cols_in_pool ); cols_in_pool.clear(); } - void c_ColumnPool::Pricing( c_DualVariables* dual, list<c_Column*>& col_list, list<double>& prices ) { } - /**************************************************************************** c_BranchAndBoundNode @@ -1892,73 +1732,55 @@ c_BranchAndBoundNode ****************************************************************************/ c_BranchAndBoundNode::c_BranchAndBoundNode( c_Controller* controller ) - : p_controller( controller ), - father( NULL ), - is_removed( false ), - is_integral( false ), - is_feasible( true ), - terminated( false ), - lb( -numeric_limits<double>::max() ), - rmp_value( -INFTY ), - iteration( NO_ITERATION ), - node_number( -1 ) +: p_controller( controller ), + p_father( NULL ), + i_node_number(-1), + i_iteration(NO_ITERATION), + i_node_number_solved(0), + i_num_iter_wo_feasibility_cut(0), + b_is_removed( false ), + b_is_integral( false ), + b_is_feasible( true ), + b_terminated( false ), + d_LB( -numeric_limits<double>::max() ), + d_RMP_value( -INFTY ), + i_level(0) {} c_BranchAndBoundNode::c_BranchAndBoundNode( c_BranchAndBoundNode* p_father, list<c_BranchAndBoundConstraint*>& add_constraints ) - : p_controller( p_father->Controller() ), - father( p_father ), - // alternative_bab_nodes( NULL ), - lb( p_father->LB() ), - is_removed( false ), - is_integral( false ), - is_feasible( true ), - terminated( false ), - rmp_value( -INFTY ), - iteration( NO_ITERATION ), - node_number( -1 ), - NumIterWithoutFeasyCut(0) +: p_controller( p_father->Controller() ), + p_father( p_father ), + i_node_number(-1), + i_iteration( NO_ITERATION ), + i_num_iter_wo_feasibility_cut(0), + b_is_removed(false), + b_is_integral(false), + b_is_feasible(true), + b_terminated(false), + d_LB(p_father->LB()), + d_RMP_value(-INFTY), + i_level(1+p_father->GetLevel_new()) { // B&B constraints from father - list<c_BranchAndBoundConstraint*>::const_iterator con; - // add new B&B constraints - for ( con=add_constraints.begin(); con!=add_constraints.end(); con++ ) + for ( auto con : add_constraints ) { - added_constraints.push_back( *con ); - + l_added_constraints.push_back( con ); // set the RMP constraint active - list<c_Constraint*>::const_iterator rmp_con; - const list<c_Constraint*>& rmp_constraints = (*con)->RMPConstraints(); - for ( rmp_con=rmp_constraints.begin(); rmp_con!=rmp_constraints.end(); ++rmp_con ) + for ( auto rmp_con : con->RMPConstraints()) { - if ( (*rmp_con)->IsGloballyValid() ) - (*rmp_con)->SetActiveAtNode( NULL ); + if ( rmp_con->IsGloballyValid() ) + rmp_con->SetActiveAtNode( NULL ); else - (*rmp_con)->SetActiveAtNode( this ); + rmp_con->SetActiveAtNode( this ); } } } c_BranchAndBoundNode::~c_BranchAndBoundNode() { - list<c_BranchAndBoundConstraint*>::const_iterator con; - for ( con=added_constraints.begin(); con!=added_constraints.end(); con++ ) - delete *con; -} - - - -c_BranchAndBoundNode* c_BranchAndBoundNode::Brother() const -{ - if ( !Father() ) - return NULL; - if ( Father()->NumSons() != 2 ) - return NULL; - list<c_BranchAndBoundNode*>::const_iterator ii; - for ( ii=Father()->sons.begin(); ii!=Father()->sons.end(); ++ii ) - if ( *ii != this ) - return *ii; - return NULL; + for (auto con : l_added_constraints ) + delete con; } string c_BranchAndBoundNode::UniquePattern() const @@ -1967,124 +1789,251 @@ string c_BranchAndBoundNode::UniquePattern() const return string("0"); string result = Father()->UniquePattern(); list<c_BranchAndBoundNode*>::const_iterator ii; - int pos; - for ( ii=Father()->sons.begin(), pos=0; ii!=Father()->sons.end(); ++ii, pos++ ) - if ( *ii != this ) + int pos = 0; + for (ii = Father()->l_sons.begin(); ii != Father()->l_sons.end(); ++ii, pos++) + { + if (*ii != this) { ostringstream buffer; buffer << pos; - result.append( buffer.str() ); + result.append(buffer.str()); break; } - return result; + } + return result; } void c_BranchAndBoundNode::AddSon( c_BranchAndBoundNode* son ) { - sons.insert( sons.end(), son ); + /* + if (_eldestChild == nullptr) + _eldestChild = son; + if (_youngestChild != nullptr) + son->setElderSibling(_youngestChild); + _youngestChild = son; + */ + l_sons.insert( l_sons.end(), son ); } - int c_BranchAndBoundNode::Level() const { - if ( !father ) + if ( !p_father ) return 0; - if ( father->sons.size() == 1 ) - return father->Level(); - return 1 + father->Level(); + if ( p_father->l_sons.size() == 1 ) + return p_father->Level(); + return 1 + p_father->Level(); } - int c_BranchAndBoundNode::DetailedLevel() const { - if ( !father ) + if ( !p_father ) return 0; - return 1 + father->DetailedLevel(); + return 1 + p_father->DetailedLevel(); } +int CGBase::c_BranchAndBoundNode::GetLevel_new() const +{ + return i_level; +} bool c_BranchAndBoundNode::IsIntegral() const { - return is_integral; + return b_is_integral; } bool c_BranchAndBoundNode::IsFeasible() const { - return is_feasible; + return b_is_feasible; } void c_BranchAndBoundNode::SetLB( double val ) { - lb = val; - //TODO: remove - if( p_controller->InfoLevel() >= 2 ) - cout << "Node-LB set to " << val << endl; + d_LB = val; p_controller->LBUpdateAllert(); } +//double c_BranchAndBoundNode::SubtreeLB() const +//{ +// // if node is infeasible +// if (!IsFeasible() && !p_controller->TimeOut()) +// return INFTY; +// +// double sLB = LB(); +// for (const auto s : sons) // consider all LBs from all sons and determine minimum +// { +// if (s->IsFeasible()) +// { +// double sub_lb = s->SubtreeLB(); +// if (sub_lb < sLB) sLB = sub_lb; +// } +// } +// return sLB; +//} + double c_BranchAndBoundNode::SubtreeLB() const { // if node is infeasible - if ( !IsFeasible() && !p_controller->TimeOut() ) + if (!IsFeasible() && !p_controller->TimeOut()) return INFTY; - // if there are no sons - if ( sons.size() == 0 ) + if (l_sons.size() == 0) return LB(); // consider all LBs from all sons and determine minimum - double erg = INFTY; - list<c_BranchAndBoundNode*>::const_iterator son; - for ( son=sons.begin(); son!=sons.end(); ++son ) + double result = INFTY; + for (auto son : l_sons) { - if ( (*son)->IsFeasible() ) + if (son->IsFeasible()) { - double sub_lb = (*son)->SubtreeLB(); - if ( sub_lb < erg ) - erg = sub_lb; + double sub_lb = son->SubtreeLB(); + if (sub_lb < result) + result = sub_lb; } } - return erg; + return result; } - double c_BranchAndBoundNode::ComputeLagrangeanLB() { - double result = -INFTY; + double result = -Controller()->LocalInfty(); if ( Controller()->UseLagrangeanLB() ) { // determine actual LP lower bound result = Controller()->RMP()->ActualObjectiveValue(); - // loop over all pricing problems double sum = 0.0; - const vector<c_PricingProblemHierarchy*>& hierarchies = Controller()->PricingProblemHierarchies(); - vector<c_PricingProblemHierarchy*>::const_iterator hier; - for ( hier=hierarchies.begin(); hier != hierarchies.end(); ++hier ) - { - c_PricingProblemHierarchy* hierachy = *hier; - //cout << "best rdc = " << hierachy->BestReducedCosts() << ", rhs = " << hierachy->RHSofConvexityConstraint() << endl; - sum += hierachy->BestReducedCosts() * hierachy->RHSofConvexityConstraint(); - } + for ( auto hier : Controller()->PricingProblemHierarchies()) + sum += hier->BestReducedCosts() * hier->RHSofConvexityConstraint(); result += sum; - // Safe new lower bound, if bound is good if ( LB() < result) SetLB( result ); } - return result; } +/* void c_BranchAndBoundNode::Get_k_th_best_BranchAndBoundContraints - ( int k, list< list<c_BranchAndBoundConstraint*> >& con_list ) + ( int k, list< list<c_BranchAndBoundConstraint*> >& con_list ) { cout << "throw - c_BranchAndBoundNode::Get_k_th_best_BranchAndBoundContraints -> Strong Branching is not yet implemented" << endl; throw; } +*/ + +void c_BranchAndBoundNode::DecideOnStrongBranchingCandidate() +{ + if (AddedConstraints().size() != 1) + { + cerr << "Something wrong with strong branching. (1)\n"; + throw; + } + c_BranchAndBoundConstraint* con = AddedConstraints().front(); + if (!con->IsStrongBranchingCandidates()) + { + cerr << "Something wrong with strong branching. (2)\n"; + throw; + } + if (Father()->NumSons() != 1) + { + cerr << "Something wrong with strong branching. (3)\n"; + throw; + } + if (Controller()->InfoLevel() >= 2) + Controller()->Info2() << "Strong Branching..." << endl; + // statistics + c_TimeInfo Timer; + Timer.Start(); + // OK, lets evaluate the candidates + auto& info2(Controller()->Info2()); + c_BranchAndBoundNode* father = Father(); + double best_lb_of_test_sons = LB(); + int idx = 0; + c_StrongBranchingCandidates* candidates = (c_StrongBranchingCandidates*)con; + vector< pair<double, list<c_BranchAndBoundNode*> > > nodes(candidates->Candidates().size()); + for (auto cand : candidates->Candidates()) + { + nodes[idx].first = 0.0; + double min_lb_of_test_sons = Controller()->LocalInfty(); + double max_lb_of_test_sons = -Controller()->LocalInfty(); + double product = 1.0; + for (auto cons_list : cand) + { + c_BranchAndBoundNode* node = Father()->CreateNewNode(cons_list); + nodes[idx].second.push_back(node); + if (Controller()->InfoLevel() >= 2) + { + info2 << "Probing { "; + for (auto ii : node->AddedConstraints()) + info2 << *ii << "; "; + info2 << "} "; + } + // exact + int old_info_level = Controller()->InfoLevel(); + Controller()->InfoLevel.SetValue(0); + c_DualVariables* dual = node->Go(numeric_limits<double>::infinity(), best_lb_of_test_sons, Controller()->PricersUsedForStrongBranching()); + delete dual; + node->SetUnsolved(); + Controller()->InfoLevel.SetValue(old_info_level); + double value = (Controller()->PricersUsedForStrongBranching() == MaxNumPricers) ? node->LB() : node->RMPValue(); + if (value < min_lb_of_test_sons) + min_lb_of_test_sons = value; + if (value > max_lb_of_test_sons) + max_lb_of_test_sons = value; + product *= (std::max(value, Controller()->LocalCG_EPS())); // ??? delta? + if (min_lb_of_test_sons < best_lb_of_test_sons) + break; // ??? + } + // Depending on Strategy: product rule or weighted sum + if (Controller()->StrongBranchingRule() == "Product") + nodes[idx].first = -product; + else // weighted sum + nodes[idx].first = -(((1 - Controller()->StrongBranchingMu()) * min_lb_of_test_sons) + (Controller()->StrongBranchingMu() * max_lb_of_test_sons)); + if (Controller()->InfoLevel() >= 2) + info2 << "Pairs MinLB " << min_lb_of_test_sons << " and MaxLB " << max_lb_of_test_sons << endl; + // increment counter + idx++; + } + // Select the best candidate + sort(nodes.begin(), nodes.end() ); + bool first = true; + for (auto node : nodes[0].second) + { + if (first) + { + l_added_constraints = node->AddedConstraints(); + } + else + { + father->AddSon(node); + Controller()->AddNode(node); + } + first = false; + } + // Throw away the other B&B nodes + for (int idx = 1; idx < (int)nodes.size(); idx++) + for (auto node : nodes[idx].second) + delete node; // remark: the B&B constraints are deleted in the destructor of c_BranchAndBoundNode + // clear memory + delete candidates; + // undo changes; very important for a following MIP run + Controller()->Update(this); + // statistics + Timer.Stop(); + Controller()->SetTimeStrongBranching(Timer.Seconds() + Controller()->InfoTimeStrongBranching()); +} + +c_StrongBranchingCandidates* c_BranchAndBoundNode::ComputeStrongBranchingCandidates(int max_k) +{ + cerr << "Serious error: Please provide member function\n" + << " c_StrongBranchingCandidates* c_BranchAndBoundNode::ComputeStrongBranchingCandidates(int k)\n" + << "in the derived class of\n" + << " c_BranchAndBoundNode\n" + << "if you intend to use strong branching.\n"; + throw; +} - -c_DualVariables* c_BranchAndBoundNode::Go( double termination_lb, double termination_ub ) +c_DualVariables* c_BranchAndBoundNode::Go( double termination_lb, double termination_ub, int NumPricersUsed ) // terminate if LB(node) goes above termination_lb // or UB(node) falls below termination_ub { @@ -2104,7 +2053,7 @@ c_DualVariables* c_BranchAndBoundNode::Go( double termination_lb, double termina c_TimeInfo timer_dual; timer_dual.Restart(); - iteration = 0; + i_iteration = 0; c_DualVariables* dual = p_controller->CreateDualVariables(); c_RMP* rmp = p_controller->RMP(); if ( p_controller->Stabilization() ) @@ -2112,8 +2061,6 @@ c_DualVariables* c_BranchAndBoundNode::Go( double termination_lb, double termina rmp->UpdateStabilization( -1, 0, dual ); //rmp->write( "mylp.lp", "LP" ); } - - //rmp->write( "mylp.lp", "LP" ); //TODO: remove double old_lb = LB(); // text if ( info_level >= 2 ) @@ -2129,17 +2076,17 @@ c_DualVariables* c_BranchAndBoundNode::Go( double termination_lb, double termina p_controller->SetInfoTimeCreateDual(p_controller->InfoTimeCreateDual() + timer_dual.Seconds()); bool end = false; bool early_termination = false; - double LLB = -INFTY; + double LLB = -Controller()->LocalInfty(); while ( !end ) { - LLB = -INFTY; + LLB = -Controller()->LocalInfty(); // count iterations - iteration++; - NumIterWithoutFeasyCut++; + i_iteration++; + i_num_iter_wo_feasibility_cut++; if ( info_level >= 3 ) { - info_3 << "Iteration " << iteration << ":\n" << flush; + info_3 << "Iteration " << i_iteration << ":\n" << flush; //rmp->write( "lp", "LP" ); info_3 << "Opt[" << p_controller->UseLPOptimizer() << "](cols=" << rmp->NumCols() << "; rows=" << rmp->NumRows() @@ -2150,10 +2097,10 @@ c_DualVariables* c_BranchAndBoundNode::Go( double termination_lb, double termina // dualen Simplex-Algorithmus zum Reoptimieren zu Beginn eines B&B-Knotens // primalen Simplex-Algorithmus danach, ist meist schneller c_TimeInfo reopt_time_info; + reopt_time_info.Start(); c_TimeInfo timer_rmp; - timer_rmp.Restart(); - + if ( p_controller->UseLPOptimizer().compare("primal" ) == 0 ) rmp->Optimize( c_RMP::ePRIMOPT ); else if ( p_controller->UseLPOptimizer().compare("dual" ) == 0 ) @@ -2168,12 +2115,10 @@ c_DualVariables* c_BranchAndBoundNode::Go( double termination_lb, double termina reopt_time_info.Stop(); timer_rmp.Stop(); - p_controller->SetInfoTimeReoptimization - ( p_controller->InfoTimeReoptimization() + reopt_time_info.Seconds() ); - + p_controller->SetInfoTimeReoptimization( p_controller->InfoTimeReoptimization() + reopt_time_info.Seconds() ); if ( info_level >= 3 ) info_3 << "(" << reopt_time_info.MilliSeconds() << "ms) "; - + double act_rmp = rmp->ActualObjectiveValue(); int status = rmp->stat(); #if CPX_VERSION >= 900 @@ -2190,27 +2135,28 @@ c_DualVariables* c_BranchAndBoundNode::Go( double termination_lb, double termina } } - if ( status == CPX_STAT_OPTIMAL_INFEAS ) + if ( status == CPX_STAT_OPTIMAL_INFEAS ) { + int helpcounter = 0; int i = 0; - while (status == CPX_STAT_OPTIMAL_INFEAS){ + while (status == CPX_STAT_OPTIMAL_INFEAS && helpcounter < 12) + { + helpcounter++; timer_rmp.Restart(); if (i == 0) - rmp->Optimize(c_RMP::ePRIMOPT); + rmp->Optimize(c_RMP::ePRIMOPT); if (i == 1) - rmp->Optimize(c_RMP::eDUALOPT); + rmp->Optimize(c_RMP::eDUALOPT); if (i == 2) - rmp->Optimize(c_RMP::eBAROPT); + rmp->Optimize(c_RMP::eBAROPT); timer_rmp.Stop(); status = rmp->stat(); - ++i; i = i % 3; - // cout << "i: " << i << endl; } - + cout << "Optimal Infeasible: Runs " << helpcounter+1 << endl; } - //rmp->write( "RMP.lp", "LP" ); + // rmp->write( "RMP.lp", "LP" ); #else if ( status != CPX_OPTIMAL && status != CPX_OPTIMAL_INFEAS ) { @@ -2225,6 +2171,58 @@ c_DualVariables* c_BranchAndBoundNode::Go( double termination_lb, double termina } } #endif + if (Controller()->CheckForInftyAndLBError()) + { + int counter = 0; //10 Iterations + if (status == CPX_STAT_INFEASIBLE || act_rmp < old_lb - 0.1) + { + int i = 0; + while ((status == CPX_STAT_INFEASIBLE || act_rmp < old_lb - 0.1) && counter < 10) + { + bool infeasy = false; + if (status == CPX_STAT_INFEASIBLE) + infeasy = true; + + if (infeasy) + cout << "Reoptimize (infeasible) iteration: " << counter << endl; + else + cout << "Reoptimize (LB-Error) iteration: " << counter << endl; + + timer_rmp.Restart(); + if (i == 0) + { + Controller()->CPLEXEnv()->setparam(CPX_PARAM_ADVIND, 0); + rmp->primopt(); + //rmp->Optimize(c_RMP::ePRIMOPT); + } + if (i == 1) + { + Controller()->CPLEXEnv()->setparam(CPX_PARAM_ADVIND, 0); + Controller()->CPLEXEnv()->setparam(CPX_PARAM_DPRIIND, CPX_DPRIIND_AUTO); + rmp->dualopt(); + //rmp->Optimize(c_RMP::eDUALOPT); + } + if (i == 2) + { + rmp->baropt(); + // rmp->Optimize(c_RMP::eBAROPT); + } + timer_rmp.Stop(); + status = rmp->stat(); + + if (status != CPX_STAT_INFEASIBLE) + { + rmp->SetObjval(rmp->objval()); + act_rmp = rmp->ActualObjectiveValue(); + cout << "Val-rmp: " << act_rmp << endl; + //rmp->OutputInOStream(cout); + } + ++i; + i = i % 3; + counter++; + } + } + } if (Controller()->InfoNumberOfRMPIterations()>1) { c_TimeInfo timer_OrgPool; @@ -2233,22 +2231,19 @@ c_DualVariables* c_BranchAndBoundNode::Go( double termination_lb, double termina timer_OrgPool.Stop(); p_controller->SetInfoTimeOrgPool(p_controller->InfoTimeOrgPool() + timer_OrgPool.Seconds()); } - - //WriteFile(); if ( act_rmp < old_lb - 0.1 && !p_controller->Stabilization() ) { // here is something wrong info_err << "#ERROR: new LB= " << act_rmp << " is smaller than LB before.\n"; + info_err << "Current Node number is " << i_node_number << " father number is " << p_father->NodeNumber() << "\n"; info_err << *rmp << "\n" << flush; - ofstream solSmallerLB( "SolutionWithSmallerLBThanBefore.txt"); - if(solSmallerLB.is_open()) - { - solSmallerLB << *rmp << endl; - } - solSmallerLB.close(); + rmp->GetDualVariables(*dual); + info_err << *dual << "\n" << flush; + rmp->write("lbsmaller_err.lp"); WriteFile(); - throw; + cin.get(); + //throw; } // text if ( info_level >= 3 ) @@ -2278,13 +2273,13 @@ c_DualVariables* c_BranchAndBoundNode::Go( double termination_lb, double termina timer_pricing.Restart(); bool feasibility_cut_added = false; do { - p_controller->Pricing( iteration, dual, new_cols ); + p_controller->Pricing( i_iteration, dual, new_cols, NumPricersUsed ); // if there are no more new columns, but the RMP is infeasible because of missing constraints feasibility_cut_added = false; bool smallGAP=false; if (LB()/rmp->ActualObjectiveValue()>Controller()->FeasyCutGap()) smallGAP=true; - if ( (new_cols.empty() && NodeNumber() >= 0 ) || NumIterWithoutFeasyCut == p_controller->FeasyCutFrequence() || iteration ==1 || smallGAP ) + if ( (new_cols.empty() && NodeNumber() >= 0 ) || i_num_iter_wo_feasibility_cut == p_controller->FeasyCutFrequence() || i_iteration ==1 || smallGAP ) { feasibility_cut_added = p_controller->FeasibilityCutting(); if ( feasibility_cut_added ) @@ -2299,19 +2294,19 @@ c_DualVariables* c_BranchAndBoundNode::Go( double termination_lb, double termina dual = p_controller->CreateDualVariables(); rmp->GetDualVariables( *dual ); } - NumIterWithoutFeasyCut=0; + i_num_iter_wo_feasibility_cut=0; } } while ( feasibility_cut_added ); timer_pricing.Stop(); - if (!father && new_cols.empty()) + if (!p_father && new_cols.empty()) p_controller->SetInfoTimeLastPPRoot(timer_pricing.Seconds()); p_controller->SetInfoTimePricing(p_controller->InfoTimePricing() + timer_pricing.Seconds()); // update stabilization - if ( p_controller->Stabilization() && ( new_cols.size() ==0 || iteration % p_controller->StabFrequence() ==0 ) ) + if ( p_controller->Stabilization() && ( new_cols.size() ==0 || i_iteration % p_controller->StabFrequence() ==0 ) ) { - last_update = rmp->UpdateStabilization( iteration, (int)new_cols.size(), dual ); + last_update = rmp->UpdateStabilization( i_iteration, (int)new_cols.size(), dual ); } // LLB if ( last_update ) @@ -2324,7 +2319,7 @@ c_DualVariables* c_BranchAndBoundNode::Go( double termination_lb, double termina #else if ( status == CPX_OPTIMAL || status == CPX_OPTIMAL_INFEAS || status == CPX_STAT_NUM_BEST ) { - p_controller->Pricing( iteration, dual, new_cols ); + p_controller->Pricing( iteration, dual, new_cols, NumPricersUsed ); // update stabilization if ( p_controller->Stabilization() && ( new_cols.size() ==0 || iteration % stab_freq ==0 ) ) @@ -2343,10 +2338,10 @@ c_DualVariables* c_BranchAndBoundNode::Go( double termination_lb, double termina timer_AddColumn.Restart(); // complex conditions for termination... end = ( last_update && ( new_cols.size() == 0 ) ) - || ( last_update && ( early_termination = ( Controller()->EarlyBranching() && this->EarlyBranching() && NumIterWithoutFeasyCut == 0) ) ) + || ( last_update && ( early_termination = ( Controller()->EarlyBranching() && this->EarlyBranching() && i_num_iter_wo_feasibility_cut == 0) ) ) || ( p_controller->TimeOut() ) - || ( rmp->RoundLBTo(LLB) >= Controller()->UB() - CG_EPS ) - || ( LB() >= 0.1 * INFTY ) + || (rmp->RoundLBTo(LLB) >= Controller()->UB() - Controller()->LocalCG_EPS()) + || (LB() >= 1 * Controller()->LocalInfty()) || ( LLB >= termination_lb ) || ( act_rmp <= termination_ub ); // add cols @@ -2354,9 +2349,8 @@ c_DualVariables* c_BranchAndBoundNode::Go( double termination_lb, double termina p_controller->AddColumns( new_cols ); else { - list<c_Column*>::const_iterator col; - for ( col=new_cols.begin(); col!=new_cols.end(); ++col ) - delete *col; + for ( auto col : new_cols ) + delete col; new_cols.clear(); } timer_AddColumn.Stop(); @@ -2392,11 +2386,11 @@ c_DualVariables* c_BranchAndBoundNode::Go( double termination_lb, double termina throw; } } - is_integral = rmp->IsIntegral() && !p_controller->TimeOut(); - is_feasible = rmp->IsFeasible() || p_controller->TimeOut(); + b_is_integral = rmp->IsIntegral() && !p_controller->TimeOut(); + b_is_feasible = rmp->IsFeasible() || p_controller->TimeOut(); // Debug: - //rmp->write( "rmp.lp", "LP" ); + // rmp->write( "lp", "LP" ); // node done; print out results if ( info_level >= 2 ) @@ -2409,16 +2403,15 @@ c_DualVariables* c_BranchAndBoundNode::Go( double termination_lb, double termina if ( LB() != rmp->ActualObjectiveValue() ) info_2 << " RMP-value: " << rmp->ActualObjectiveValue() << " (gap = " << RelativeDifference( LB(), rmp->ActualObjectiveValue() ) << ")\n"; - info_2 << " iterations: " << iteration << "\n"; - if ( is_feasible ) - info_2 << " integral: " << ( is_integral ? "yes" :"no" ) << "\n" << flush; + info_2 << " iterations: " << i_iteration << "\n"; + if ( b_is_feasible ) + info_2 << " integral: " << ( b_is_integral ? "yes" :"no" ) << "\n" << flush; else info_2 << " feasible: no\n" << flush; - if ( info_level >= 2 ){ + if ( info_level >= 2 ) + { info_2 << *rmp << flush; - } - if ( info_level >= 3 ){ info_2 << *dual << flush; rmp->write( "mylp.lp", "LP" ); } @@ -2432,26 +2425,23 @@ void c_BranchAndBoundNode::WriteFile() bool c_BranchAndBoundNode::AnySonSolved() const { - list<c_BranchAndBoundNode*>::const_iterator son_iter; - for ( son_iter = sons.begin(); son_iter != sons.end(); ++son_iter ) - if ( (*son_iter)->IsSolved() ) + for ( auto son_iter :l_sons ) + if ( son_iter->IsSolved() ) return true; return false; } - - bool c_BranchAndBoundNode::GetSeparationConstraints( list<c_BranchAndBoundConstraint*>& con_list ) { int num_cuts = 0; c_RMP* rmp = Controller()->RMP(); c_TimeInfo sep_time_info; + sep_time_info.Start(); list<c_Constraint*> cuts; p_controller->Separation( rmp, cuts, Level() ); sep_time_info.Stop(); - p_controller->SetInfoTimeSeparation - (p_controller->InfoTimeSeparation() + sep_time_info.Seconds()); + p_controller->SetInfoTimeSeparation(p_controller->InfoTimeSeparation() + sep_time_info.Seconds()); num_cuts = (int)cuts.size(); if ( num_cuts ) @@ -2502,7 +2492,7 @@ void *cbhandle; glob_controller->SetUB( best_integer ); } // Check Timeout - if ( glob_controller->TimeOut() || mip_timer.Seconds() >= mip_time_limit ) + if ((glob_controller->TimeOut() && !glob_controller->SolveAsMipAfterTimeOut()) || mip_timer.Seconds() >= mip_time_limit) { glob_controller->Info1() << "... MIP timeout." << endl; return (status=1) ; @@ -2512,92 +2502,136 @@ void *cbhandle; // solving the node with a MIP solver is a heuristic (to improve UB) -double c_BranchAndBoundNode::SolveAsMIP( long time_limit ) -{ - double result = INFTY*INFTY; +//double c_BranchAndBoundNode::SolveAsMIP( long time_limit ) +//{ +// double result = Controller()->LocalInfty()*Controller()->LocalInfty(); +// c_RMP* rmp = Controller()->RMP(); +// c_Controller* ctlr = Controller(); +// //rmp->write("RootReady.lp", "LP"); +// cout << "Solve as MIP..." << endl; +// mip_time_limit = time_limit; +// mip_timer.Start(); +// //ctlr->RMP()->Optimize( c_RMP::eDUALOPT ); +// ctlr->FixPartOfSolution(); +// // set my callback function +// CPXENVptr env = (CPXENVptr) ctlr->CPLEXEnv()->getEnv(); +// //CPXsetdblparam(env, CPX_PARAM_EPGAP, 0.001); +// CPXsetintparam(env, CPX_PARAM_THREADS, ctlr->NumThreadsRMP()); +// CPXsetintparam(env, CPX_PARAM_MIPSEARCH, 2); +// glob_controller = ctlr; +// //CPXsetmipcallbackfunc( env, mymipcallback, NULL); +// // remove redundant cols (if any) +// if ( ctlr->MaxRatioColsToRowsMIP() < ctlr->MaxRatioColsToRows() ) +// { +// if ( ctlr->InfoLevel() >= 1 ) +// { +// ctlr->Info1() << "Invoking MIP Solver of cplex... " << endl; +// ctlr->Info1() << "... Reducing columns for MIP-Solver... " << endl; +// } +// double merk = Controller()->MaxRatioColsToRows(); +// ctlr->SetMaxRatioColsToRows( ctlr->MaxRatioColsToRowsMIP() ); +// Controller()->RMP()->Optimize( c_RMP::eDUALOPT ); +// ctlr->OrganizePool(); +// ctlr->SetMaxRatioColsToRows( merk ); +// } +// // Do optimization +// double old_UB = ctlr->UB(); +// rmp->OptimizeMIP( c_RMP::ePRIMOPT ); +// result = rmp->mipobjval(); +// mip_timer.Stop(); +// ctlr->SetInfoTimeMIPForUB(ctlr->InfoTimeMIPForUB() + mip_timer.Seconds()); +// // quality? +// //if (rmp->IsFeasible() && (result - ctlr->UB()) < -ctlr->LocalCG_EPS()) // new best solution found +// if (rmp->IsFeasible() && (result - old_UB) < -ctlr->LocalCG_EPS()) // new best solution found +// { +// // Store MIP-Solution +// ctlr->SetBestSolution( +// ctlr->CreateSolution( Controller()->RMP(), +// NULL /* no dual info available */ ) ); +// if ( ctlr->InfoLevel() >= 1 ) +// { +// ctlr->Info1() << "... MIP time: " << mip_timer.Seconds() << "s" << endl; +// if ( ctlr->BestSolution() && ctlr->InfoPrintBestSolution() ) +// ctlr->Info1() << *ctlr->BestSolution() << flush; +// } +// } +// ctlr->UnfixSolution(); +// // change back to LP +// rmp->chgprobtype( CPLEXProblem::LP ); +// // turn off callback +// CPXsetmipcallbackfunc( env, NULL, NULL); +// glob_controller = NULL; +// return result; +//} + +double c_BranchAndBoundNode::SolveAsMIP(long time_limit) +{ + double result = Controller()->LocalInfty() * Controller()->LocalInfty(); c_RMP* rmp = Controller()->RMP(); c_Controller* ctlr = Controller(); //rmp->write("RootReady.lp", "LP"); cout << "Solve as MIP..." << endl; mip_time_limit = time_limit; mip_timer.Start(); - ctlr->RMP()->Optimize( c_RMP::eDUALOPT ); - ctlr->FixPartOfSolution(); - // set my callback function - CPXENVptr env = (CPXENVptr) ctlr->CPLEXEnv()->getEnv(); - //CPXsetdblparam(env, CPX_PARAM_EPGAP, 0.001); - CPXsetintparam(env, CPX_PARAM_THREADS, 1); - glob_controller = ctlr; - CPXsetmipcallbackfunc( env, mymipcallback, NULL); + + // remove redundant cols (if any) - if ( ctlr->MaxRatioColsToRowsMIP() < ctlr->MaxRatioColsToRows() ) + if (ctlr->MaxRatioColsToRowsMIP() < ctlr->MaxRatioColsToRows()) { - if ( ctlr->InfoLevel() >= 1 ) + if (ctlr->InfoLevel() >= 1) { ctlr->Info1() << "Invoking MIP Solver of cplex... " << endl; ctlr->Info1() << "... Reducing columns for MIP-Solver... " << endl; } double merk = Controller()->MaxRatioColsToRows(); - ctlr->SetMaxRatioColsToRows( ctlr->MaxRatioColsToRowsMIP() ); - Controller()->RMP()->Optimize( c_RMP::eDUALOPT ); + ctlr->SetMaxRatioColsToRows(ctlr->MaxRatioColsToRowsMIP()); + Controller()->RMP()->Optimize(c_RMP::eDUALOPT); ctlr->OrganizePool(); - ctlr->SetMaxRatioColsToRows( merk ); + ctlr->SetMaxRatioColsToRows(merk); } + // Do optimization - rmp->OptimizeMIP( c_RMP::ePRIMOPT ); + rmp->OptimizeMIP(c_RMP::ePRIMOPT); result = rmp->mipobjval(); mip_timer.Stop(); ctlr->SetInfoTimeMIPForUB(ctlr->InfoTimeMIPForUB() + mip_timer.Seconds()); - // quality? - if ( rmp->IsFeasible() && ( result - ctlr->UB() ) < CG_EPS ) // new best solution found - { - // Store MIP-Solution - ctlr->SetBestSolution( - ctlr->CreateSolution( Controller()->RMP(), - NULL /* no dual info available */ ) ); - if ( ctlr->InfoLevel() >= 1 ) - { - ctlr->Info1() << "... MIP time: " << mip_timer.Seconds() << "s" << endl; - if ( ctlr->BestSolution() && ctlr->InfoPrintBestSolution() ) - ctlr->Info1() << *ctlr->BestSolution() << flush; - } - } - ctlr->UnfixSolution(); + + // change back to LP - rmp->chgprobtype( CPLEXProblem::LP ); + rmp->chgprobtype(CPLEXProblem::LP); // turn off callback - CPXsetmipcallbackfunc( env, NULL, NULL); - glob_controller = NULL; return result; } - void c_BranchAndBoundNode::GetConstraints( list<c_BranchAndBoundConstraint*>& con_list ) const { - con_list = added_constraints; + con_list = l_added_constraints; c_BranchAndBoundNode* p_father = Father(); while ( p_father ) { - list<c_BranchAndBoundConstraint*>::const_iterator iter; const list<c_BranchAndBoundConstraint*>& added_list = p_father->AddedConstraints(); - for( iter=added_list.begin(); iter!=added_list.end(); ++iter ) - con_list.push_back( *iter ); - p_father = p_father->father; + for(auto iter : added_list ) + con_list.push_back( iter ); + p_father = p_father->p_father; } } void c_BranchAndBoundNode::OutputInOStream( ostream& s ) const { - s << "node LB = " << lb << "\n" + s << "node LB = " << d_LB << "\n" << "added node constraints: "; - if ( !added_constraints.empty() ) + if ( !l_added_constraints.empty() ) { - list<c_BranchAndBoundConstraint*>::const_iterator con; - for ( con=added_constraints.begin(); con!=added_constraints.end(); ++con ) - if ( con !=added_constraints.begin() ) - s << ", " << **con; - else - s << **con; + bool first = true; + for (auto con : l_added_constraints) + { + if (first) + { + s << ", "; + first = false; + } + s << *con; + } } else s << "(none)"; @@ -2609,10 +2643,9 @@ c_Constraint ****************************************************************************/ -c_Constraint::c_Constraint( c_Controller* controller, - char sense, double rhs, double estimator_dual, - int type, int n1, int n2, int n3 ) - : _type( type ), +c_Constraint::c_Constraint( c_Controller* controller, char sense, double rhs, double estimator_dual, int type, int n1, int n2, int n3 ) +: p_controller(controller), + _type( type ), _n1( n1 ), _n2( n2 ), _n3( n3 ), @@ -2621,10 +2654,10 @@ c_Constraint::c_Constraint( c_Controller* controller, _estimator_dual( estimator_dual ), col_yps_plus( NULL ), col_yps_minus( NULL ), - p_controller( controller ), active_at_node( NULL ) { controller->AddConstraint( this ); + } @@ -2697,7 +2730,8 @@ c_RMP::c_RMP( c_Controller* controller ) actual_obj_value( INFTY ), info_number_of_simplex_iterations( 0 ), num_removed_constraints(0), - num_added_constraints(0) + num_added_constraints(0), + help_mipname("") { Controller()->CPLEXEnv()->setparam( CPX_PARAM_THREADS, Controller()->NumThreadsRMP()); } @@ -2706,17 +2740,16 @@ c_RMP::~c_RMP() { // do not delete c_Constraint* in constraints // they are deleted by the controller + remove(help_mipname.c_str()); } - void c_RMP::Reset() { actual_iteration = -1; - actual_obj_value = INFTY; + actual_obj_value = Controller()->LocalInfty(); info_number_of_simplex_iterations = 0; } - void c_RMP::GetDualVariables( c_DualVariables& dual ) const { getPi( dual.CPXDual() ); @@ -2742,6 +2775,7 @@ void c_RMP::Optimize( e_solver solver ) break; case eBAROPT: // barrier solver + Controller()->CPLEXEnv()->setparam(CPX_PARAM_SOLUTIONTYPE, 2); // hybbaropt(PRIMOPT); if (p_controller->InfoNumberOfRMPIterations()<200) Controller()->CPLEXEnv()->setparam(CPX_PARAM_BAREPCOMP, 0.000001); @@ -2769,18 +2803,12 @@ void c_RMP::Optimize( e_solver solver ) info_err << "#what = " << e.what() << endl; info_err << "#writing LP to error.lp" << endl; write("error.lp"); - if( stat() != CPX_STAT_OPTIMAL ){ - Controller()->CPLEXEnv()->setparam(CPX_PARAM_NUMERICALEMPHASIS, true); - primopt(); - if( stat() != CPX_STAT_OPTIMAL){ - cout << "Reoptimization with emphasis on numerical issues did not help. RMP cannot be solved to optimality." << endl; - } - } } // count simplex iterations info_number_of_simplex_iterations += getitcnt(); #if CPX_VERSION >= 900 + if ( stat() == CPX_STAT_OPTIMAL || stat() == CPX_STAT_OPTIMAL_INFEAS ) actual_obj_value = objval(); else @@ -2788,18 +2816,19 @@ void c_RMP::Optimize( e_solver solver ) cout << "#ERROR: Optimizer errorcode " << stat() << endl; switch ( stat() ) { - case CPX_STAT_UNBOUNDED: actual_obj_value = -INFTY; break; - case CPX_STAT_INFEASIBLE: actual_obj_value = INFTY * INFTY; break; - case CPX_STAT_NUM_BEST: actual_obj_value = -INFTY; break; - //case CPX_STAT_FEASIBLE_RELAXED: actual_obj_value = -INFTY; break; - //case CPX_STAT_OPTIMAL_RELAXED: actual_obj_value = -INFTY; break; + case CPX_STAT_UNBOUNDED: actual_obj_value = -Controller()->LocalInfty(); break; + case CPX_STAT_INFEASIBLE: actual_obj_value = Controller()->LocalInfty() * Controller()->LocalInfty(); break; + + case CPX_STAT_NUM_BEST: actual_obj_value = -Controller()->LocalInfty(); break; + //case CPX_STAT_FEASIBLE_RELAXED: actual_obj_value = -Controller()->LocalInfty(); break; + //case CPX_STAT_OPTIMAL_RELAXED: actual_obj_value = -Controller()->LocalInfty(); break; // Abbruch case CPX_STAT_ABORT_IT_LIM: case CPX_STAT_ABORT_TIME_LIM: case CPX_STAT_ABORT_OBJ_LIM: - case CPX_STAT_ABORT_USER: actual_obj_value = INFTY; break; + case CPX_STAT_ABORT_USER: actual_obj_value = Controller()->LocalInfty(); break; // XXX - case CPX_STAT_INForUNBD: actual_obj_value = INFTY; break; + case CPX_STAT_INForUNBD: actual_obj_value = Controller()->LocalInfty(); break; } write( "error.lp", "LP" ); } @@ -2811,15 +2840,15 @@ void c_RMP::Optimize( e_solver solver ) cout << "#ERROR: Optimizer errorcode " << stat() << endl; switch ( stat() ) { - case CPX_INFEASIBLE: actual_obj_value = INFTY * INFTY; break; - case CPX_UNBOUNDED: actual_obj_value = -INFTY; break; - case CPX_OBJ_LIM: actual_obj_value = INFTY; break; + case CPX_INFEASIBLE: actual_obj_value = Controller()->LocalInfty() * Controller()->LocalInfty(); break; + case CPX_UNBOUNDED: actual_obj_value = -Controller()->LocalInfty(); break; + case CPX_OBJ_LIM: actual_obj_value = Controller()->LocalInfty(); break; case CPX_IT_LIM_FEAS: actual_obj_value = objval(); break; - case CPX_IT_LIM_INFEAS: actual_obj_value = INFTY; break; + case CPX_IT_LIM_INFEAS: actual_obj_value = Controller()->LocalInfty(); break; case CPX_TIME_LIM_FEAS: actual_obj_value = objval(); break; - case CPX_TIME_LIM_INFEAS: actual_obj_value = INFTY; break; + case CPX_TIME_LIM_INFEAS: actual_obj_value = Controller()->LocalInfty(); break; case CPX_NUM_BEST_FEAS: actual_obj_value = objval(); break; - case CPX_NUM_BEST_INFEAS: actual_obj_value = INFTY; break; + case CPX_NUM_BEST_INFEAS: actual_obj_value = Controller()->LocalInfty(); break; case CPX_OPTIMAL_INFEAS: actual_obj_value = objval(); break; // Abbruch case CPX_ABORT_FEAS: @@ -2828,9 +2857,9 @@ void c_RMP::Optimize( e_solver solver ) case CPX_ABORT_PRIM_INFEAS: case CPX_ABORT_PRIM_DUAL_INFEAS: case CPX_ABORT_PRIM_DUAL_FEAS: - case CPX_ABORT_CROSSOVER: actual_obj_value = INFTY; break; + case CPX_ABORT_CROSSOVER: actual_obj_value = Controller()->LocalInfty(); break; // XXX - case CPX_INForUNBD: actual_obj_value = INFTY; break; + case CPX_INForUNBD: actual_obj_value = Controller()->LocalInfty(); break; } write( "error.lp", "LP" ); } @@ -2845,27 +2874,126 @@ void c_RMP::SetRequiredVarsIntegral() } +void c_RMP::SetConstraintSense(int contraintIndex, char sense) { + + c_Controller* ctlr = Controller(); + + char _sense[] = { sense }; + int _indices[] = { contraintIndex }; + + //CPXchgsense(ctlr->CPLEXEnv()->getEnv(), (CPXLPptr)this, 1, _indices, _sense); + chgsense(1, _indices, _sense); + +} + + // try to solve the LP as a MIP (->UB) void c_RMP::OptimizeMIP( e_solver solver ) { + c_Controller* ctlr = Controller(); + const bool accelerate_MIP_solver = ctlr->MIPAccelerate(); + chgprobtype( CPLEXProblem::MIP ); char integral[] = { 'I' }; for ( int i=0; (i<numcols()); i++ ) chgctype( 1, &i, integral ); try { - d_mip_objval = INFTY*INFTY; + d_mip_objval = Controller()->LocalInfty() * Controller()->LocalInfty(); + vector<pair<int, double> > columns_in_MIPsol; + if ( !accelerate_MIP_solver ) + { + CPXENVptr env = (CPXENVptr) ctlr->CPLEXEnv()->getEnv(); + //set parameters + CPXsetintparam(env, CPX_PARAM_THREADS, ctlr->NumThreadsRMP()); + if (ctlr->MIPScreenOutput()) + CPXsetintparam(env, CPX_PARAM_SCRIND, CPX_ON); + CPXsetdblparam(env, CPX_PARAM_TILIM, ctlr->MIPMaxSolutionTimeSec()); + CPXsetdblparam(env, CPX_PARAM_CUTUP, ctlr->UB() + 1); + + // We use the same CPLEX problem mipopt(); + //Reset parameters for lp + CPXsetintparam(env, CPX_PARAM_SCRIND, CPX_OFF); + CPXsetdblparam(env, CPX_PARAM_TILIM, ctlr->MaxSolutionTimeSec()); + CPXsetdblparam(env, CPX_PARAM_CUTUP, INFTY); d_mip_objval = CPLEXProblem::mipobjval(); // Call the function of the base class + if ((d_mip_objval - ctlr->UB()) < -ctlr->LocalCG_EPS()) // new best solution found + { + for (int i = 0; i < numcols(); i++) + { + double val_x; + getX(&val_x, i, i); + if (val_x > CG_EPS) + columns_in_MIPsol.push_back(make_pair(i, val_x)); + } + } + } + else + { + // We write out a new problem, read it in, and solve the new problem + // CPLEX uses all the sophisticated preprocessing etc. here + // + // Warning: this causes a segmentation fault on the MOGON cluster for unknown reasons + // + if (help_mipname == "") + { + stringstream mipname; + srand((int)time(nullptr)); + int help = rand(); + //mipname << "./tmp/MIP" << (int)(10000 * ctlr->InfoTimeReoptimization()) << (int)(10000 * ctlr->InfoTimePricing()) << help << ".mps"; + mipname << "./tmp/MIP" << ".mps"; + help_mipname = mipname.str().c_str(); + } + write(help_mipname.c_str()); + CPLEXEnvironment env_mip; + if (ctlr->MIPScreenOutput()) + env_mip.setparam(CPX_PARAM_SCRIND, CPX_ON); + env_mip.setparam(CPX_PARAM_THREADS, ctlr->NumThreadsRMP()); + double soltimemip = ctlr->MIPMaxSolutionTimeSec(); + env_mip.setparam(CPX_PARAM_TILIM, soltimemip); + env_mip.setparam(CPX_PARAM_CUTUP, ctlr->UB()); + CPLEXProblem the_mip(env_mip); + the_mip.read(help_mipname.c_str()); + //the_mip.write("MIP_write.lp"); + the_mip.mipopt(); + d_mip_objval = the_mip.mipobjval(); + if ((d_mip_objval - ctlr->UB()) < -ctlr->LocalCG_EPS()) // new best solution found + { + for (int i = 0; i < the_mip.numcols(); i++) + { + double val_x; + the_mip.getX(&val_x, i, i); + if (val_x > CG_EPS) + columns_in_MIPsol.push_back(make_pair(i, val_x)); + } + } + env_mip.terminate(); + } + if (ctlr->InfoLevel() >= 1) + ctlr->Info1() << "... MIP time: " << mip_timer.Seconds() << "s" << endl; + + + c_Solution* sol = ctlr->CreateSolutionFromMIP(Controller()->RMP(), columns_in_MIPsol); + if (sol != nullptr) //if no dummys in sol + { + if (d_mip_objval < ctlr->UB()) //if new UB + { + ctlr->SetBestSolution(sol); + ctlr->SetUB(d_mip_objval); + if (ctlr->InfoLevel() >= 1) + ctlr->Info1() << "... New best solution with value: " << d_mip_objval << endl; + } + } } catch(CPLEXException e) { - cout << e.what() << endl; + if(e.errcode() != 1217) //No onscreen Info if no solution exists + cout << e.what() << endl; } // LATER: // chgprobtype( CPLEXProblem::LP ); } - double c_RMP::mipobjval() const { return d_mip_objval; @@ -2881,20 +3009,18 @@ bool c_RMP::IsFeasible() const if ( stat() != CPX_OPTIMAL ) return false; #endif - int i; double x; double obj; - for ( i=0; (i<numcols()); i++ ) + for (int i = 0, maxi = numcols(); i < maxi; ++i) { getX( &x, i, i ); getobj( &obj, i, i ); - if ( ( x > CG_EPS ) && ( obj >= INFTY ) ) + if ((x > Controller()->LocalCG_EPS()) && (obj >= Controller()->LocalInfty())) return false; } return true; } - bool c_RMP::IsIntegral() const { #if CPX_VERSION >= 900 @@ -2904,23 +3030,19 @@ bool c_RMP::IsIntegral() const if ( stat() != CPX_OPTIMAL ) return false; #endif - int i; double x; - for ( i=0; (i<numcols()); i++ ) + for (int i = 0, maxi = numcols(); i < maxi; ++i) { getX( &x, i, i ); - if ( ( ( x - floor( x ) ) > CG_EPS ) - && ( ( ceil( x ) - x ) > CG_EPS ) ) + if (((x - floor(x)) > Controller()->LocalCG_EPS()) + && ((ceil(x) - x) > Controller()->LocalCG_EPS())) return false; } return true; } - - void c_RMP::AddColumns( const list<c_Column*>& col_list ) { - class Y{ public: Y():ccnt_max(0),obj(NULL),cmatbeg(NULL),lb(NULL),ub(NULL),nzcnt_max(0),cmatind(NULL),cmatval(NULL){}; @@ -2939,8 +3061,6 @@ void c_RMP::AddColumns( const list<c_Column*>& col_list ) static Y y; // determine needed space int ccnt = (int)col_list.size(); - int col_cnt = 0; - list<c_Column*>::const_iterator col; int nzcnt_estimation = ccnt * NumRows(); /* @@ -2971,25 +3091,23 @@ void c_RMP::AddColumns( const list<c_Column*>& col_list ) // do insertion int col_no=0; int nzcnt = 0; - for ( col=col_list.begin(); col!=col_list.end(); ++col ) + for ( auto col : col_list ) { y.cmatbeg[col_no] = nzcnt; - nzcnt += (*col)->GetCPXColumn( y.obj[col_no], &y.cmatind[nzcnt], &y.cmatval[nzcnt], y.lb[col_no], y.ub[col_no], NULL ); + nzcnt += col->GetCPXColumn( y.obj[col_no], &y.cmatind[nzcnt], &y.cmatval[nzcnt], y.lb[col_no], y.ub[col_no], NULL ); col_no++; // do statistics if ( Controller()->InfoStatisticsColumns() ) - Controller()->RecordStatisticsColumns( (*col)->InfoColumnType() ); + Controller()->RecordStatisticsColumns( col->InfoColumnType() ); } addcols( ccnt, nzcnt, y.obj, y.cmatbeg, y.cmatind, y.cmatval, y.lb, y.ub, NULL ); } - void c_RMP::GetReducedCosts( double c_red[] ) { getDj(c_red, 0, -1 ); } - void c_RMP::GetDualPrices( double pi[] ) { getPi(pi, 0, -1 ); @@ -3000,7 +3118,6 @@ void c_RMP::GetRHS( double rhs[] ) getRHS(rhs, 0, -1 ); } - void c_RMP::RemoveAllBranchAndBoundConstraints() { } @@ -3082,19 +3199,12 @@ bool c_RMP::UpdateStabilization( int iteration, int num_cols_pricing, c_DualVari obj_cnt++; } } - - - chgbds( col_cnt, col_ind, col_ul, col_bnd ); for ( int i=0;i<obj_cnt;i++ ) chgcoef( -1, obj_ind[i], obj_val[i] ); - - return last_update; } - - void c_RMP::Update( c_BranchAndBoundNode* node ) { // Info @@ -3105,10 +3215,9 @@ void c_RMP::Update( c_BranchAndBoundNode* node ) // loop over all constraints and eliminate those that are not valid here list<c_Constraint*> non_active; - vector<c_Constraint*>::iterator con; - for ( con=constraints.begin(); con!=constraints.end(); ++con ) - if ( !(*con)->IsActiveAtNode( node ) ) - non_active.push_back( *con ); + for (auto con:constraints ) + if ( !con->IsActiveAtNode( node ) ) + non_active.push_back( con ); if ( !non_active.empty() ) RemoveConstraints( non_active ); @@ -3121,7 +3230,6 @@ void c_RMP::Update( c_BranchAndBoundNode* node ) num_added_constraints += (int)constraint_list.size(); } - void c_RMP::AddBranchAndBoundConstraints( const list<c_BranchAndBoundConstraint*>& constraint_list ) { class X{ @@ -3171,17 +3279,26 @@ void c_RMP::AddBranchAndBoundConstraints( const list<c_BranchAndBoundConstraint* double new_lb = x.lb[pos]; double new_ub = x.ub[pos]; list<c_BranchAndBoundConstraint*>::const_iterator con; - for ( con=constraint_list.begin(); con!=constraint_list.end(); ++con ) + for ( auto con : constraint_list ) { - (*con)->ComputeBounds( col, new_lb, new_ub ); + con->ComputeBounds( col, new_lb, new_ub ); + if (new_lb >= new_ub) + { + x.lb[pos] = new_lb; + x.ub[pos] = new_lb; + break; + } + if ( x.lb[pos] < new_lb ) x.lb[pos] = new_lb; if ( x.ub[pos] > new_ub ) x.ub[pos] = new_ub; + /* if ( new_lb >= new_ub ) break; + */ } - pos++; + ++pos; } } @@ -3194,8 +3311,8 @@ void c_RMP::AddBranchAndBoundConstraints( const list<c_BranchAndBoundConstraint* { if ( x.ub[i] == 0.0 && x.lb[i]==0.0 ) { - x.cost[i] = INFTY; - x.ub[i] = INFTY; + x.cost[i] = Controller()->LocalInfty(); + x.ub[i] = Controller()->LocalInfty(); //x.ub[i] = 0; //TODO -> TEST } @@ -3208,22 +3325,11 @@ void c_RMP::AddBranchAndBoundConstraints( const list<c_BranchAndBoundConstraint* // RMP constraints in RMP list<c_Constraint*> rmp_constraints_to_add; - list<c_BranchAndBoundConstraint*>::const_iterator bb_con; - for ( bb_con=constraint_list.begin(); bb_con!=constraint_list.end(); ++bb_con ) - { - // loop over all rmp constraints - list<c_Constraint*>::const_iterator rmp_con; - const list<c_Constraint*>& rmp_constraints = (*bb_con)->RMPConstraints(); - for ( rmp_con=rmp_constraints.begin(); - rmp_con!=rmp_constraints.end(); - ++rmp_con ){ - if ( !IsConstraintInRMP( *rmp_con ) ){ - rmp_constraints_to_add.push_back( *rmp_con ); - } - } - } + for ( auto bb_con:constraint_list ) + for ( auto rmp_con : bb_con->RMPConstraints()) + if ( !IsConstraintInRMP( rmp_con ) ) + rmp_constraints_to_add.push_back( rmp_con ); AddConstraints( rmp_constraints_to_add ); - } bool c_RMP::IsConstraintInRMP( const c_Constraint* con ) const @@ -3262,15 +3368,13 @@ void c_RMP::RemoveConstraint( c_Constraint* con ) } } - void c_RMP::RemoveConstraints( list<c_Constraint*>& con_list ) { list<c_Constraint*>::const_iterator con; - for ( con=con_list.begin(); con!=con_list.end(); ++con ) - RemoveConstraint( *con ); + for ( auto con : con_list) + RemoveConstraint( con ); } - int c_RMP::AddConstraint( c_Constraint* con ) { list<c_Constraint*> con_list; @@ -3279,7 +3383,6 @@ int c_RMP::AddConstraint( c_Constraint* con ) return con->Index(); } - void c_RMP::AddConstraints( list<c_Constraint*>& con_list ) { class Z{ @@ -3310,10 +3413,8 @@ void c_RMP::AddConstraints( list<c_Constraint*>& con_list ) if ( con_list.size() == 0 ) return; - list<c_Constraint*>::iterator constr; - for ( constr=con_list.begin(); constr!=con_list.end(); ++constr ) + for ( auto con : con_list ) { - c_Constraint* con = *constr; #ifdef _DEBUG if ( IndexExists( con->_type, con->_n1, con->_n2, con->_n3 ) ) { @@ -3324,7 +3425,6 @@ void c_RMP::AddConstraints( list<c_Constraint*>& con_list ) constraints.push_back( con ); NewIndex( con->_type, con->_n1, con->_n2, con->_n3 ); } - // Check, if there are already integrated columns int n = NumCols(); if ( n < 1 ) @@ -3375,9 +3475,8 @@ void c_RMP::AddConstraints( list<c_Constraint*>& con_list ) // read new constraints rcnt = 0; - for ( constr=con_list.begin(); constr!=con_list.end(); ++constr ) + for ( auto con : con_list ) { - c_Constraint* con = *constr; z.rhs[rcnt] = con->Rhs(); z.sense[rcnt] = con->Sense(); z.anz_entries_row[rcnt] = 0; @@ -3392,9 +3491,8 @@ void c_RMP::AddConstraints( list<c_Constraint*>& con_list ) z.ind[i] = -1; int col_nzcnt = col->GetCPXColumn(cost,z.ind,z.val,lb,ub,NULL); rcnt = 0; - for ( constr=con_list.begin(); constr!=con_list.end(); ++constr ) + for ( auto con:con_list ) { - c_Constraint* con = *constr; int idx = con->Index(); // integrate additional values for ( int i=0; i<col_nzcnt; i++ ) @@ -3438,12 +3536,9 @@ void c_RMP::AddConstraints( list<c_Constraint*>& con_list ) } - - int c_RMP::ConstraintIndex( int type, int n1, int n2, int n3 ) const { - int index = Index( type, n1, n2, n3 ); - return index; + return Index( type, n1, n2, n3 ); } bool c_RMP::ConstraintIndexExists( int i1, int i2, int i3, int i4 ) const @@ -3451,7 +3546,7 @@ bool c_RMP::ConstraintIndexExists( int i1, int i2, int i3, int i4 ) const return IndexExists( i1, i2, i3, i4 ) ; } -int c_RMP::NumConstraints( void ) const +int c_RMP::NumConstraints(void) const { return NumIndices(); } @@ -3464,12 +3559,10 @@ c_Constraint* c_RMP::Constraint( int i1, int i2, int i3, int i4 ) const return constraints[Index( i1, i2, i3, i4)]; } - void c_RMP::ConstructEmptyRMP( CPLEXProblem::ObjSense objsen ) { int num_rows = NumConstraints(); int num_cols = 0; // add cols later - double* obj = new double[1]; double* rhs = new double[ num_rows ]; char* sense = new char[ num_rows ]; @@ -3479,52 +3572,44 @@ void c_RMP::ConstructEmptyRMP( CPLEXProblem::ObjSense objsen ) double* cmatval = new double[1]; double* lb = new double[1]; double* ub = new double[1]; - - vector< c_Constraint* >::iterator constraint; - for ( constraint=constraints.begin(); constraint != constraints.end(); ++constraint ) + for ( auto constraint : constraints ) { - rhs[ (*constraint)->Index() ] = (*constraint)->Rhs(); - sense[ (*constraint)->Index() ] = (*constraint)->Sense(); + rhs[ constraint->Index() ] = constraint->Rhs(); + sense[ constraint->Index() ] = constraint->Sense(); } - // format: ...beg cnt ind val lb ub copylpdata( num_cols, num_rows, objsen, obj, rhs, sense, cmatbeg, cmatcnt, cmatind, cmatval, lb, ub ); - // add surplus and slack column for stabilized column generation list< c_Column* > col_list; - for ( constraint=constraints.begin(); constraint != constraints.end(); ++constraint ) - if ( (*constraint)->Stabilization() ) + for (auto constraint : constraints) + { + if (constraint->Stabilization()) { - (*constraint)->CreateYpsColumns(); - if ( (*constraint)->YpsMinusColumn() ) - col_list.insert( col_list.end(), (*constraint)->YpsMinusColumn() ); - if ( (*constraint)->YpsPlusColumn() ) - col_list.insert( col_list.end(), (*constraint)->YpsPlusColumn() ); + constraint->CreateYpsColumns(); + if (constraint->YpsMinusColumn()) + col_list.insert(col_list.end(), constraint->YpsMinusColumn()); + if (constraint->YpsPlusColumn()) + col_list.insert(col_list.end(), constraint->YpsPlusColumn()); } - Controller()->AddColumns( col_list ); - // clean memory - delete[] obj; - delete[] rhs; - delete[] sense; - delete[] cmatbeg; - delete[] cmatcnt; - delete[] cmatind; - delete[] cmatval; - delete[] lb; - delete[] ub; + } + Controller()->AddColumns( col_list ); + // clean memory + delete[] obj; + delete[] rhs; + delete[] sense; + delete[] cmatbeg; + delete[] cmatcnt; + delete[] cmatind; + delete[] cmatval; + delete[] lb; + delete[] ub; } - - - void c_RMP::OutputInOStream( ostream& s ) const { s << "RMP: \n"; - - vector< c_Constraint*>::const_iterator constraint; - for ( constraint = constraints.begin(); constraint != constraints.end(); ++constraint ) + for ( auto con : constraints ) { - c_Constraint* con = *constraint; s << con->Index() << ": " << *con << endl; if ( !IsConstraintInRMP( con ) ) s << "#Error: constraint not in RMP." << endl; @@ -3592,14 +3677,11 @@ c_Column::c_Column( c_Controller* controller ) num_in_rmp( -1 ) {} - - void c_Column::OutputInOStream( ostream& s ) const { s << "Column: undefined Type"; } - double c_Column::ReducedCosts( const c_DualVariables& dual ) { class X{ @@ -3636,7 +3718,45 @@ double c_Column::ReducedCosts( const c_DualVariables& dual ) return cost; } +double c_Column::CoutReducedCosts(const c_DualVariables& dual) +{ + class X { + public: + X() :ind(NULL), val(NULL), nzcnt_max(0) {}; + ~X() { delete[] ind; delete[] val; }; + int* ind; + double* val; + int nzcnt_max; + }; + // static init + static X x; + //int nzcnt_max = -1; + //static int* ind = NULL; + //static double* val = NULL; + // enough memory available? + double cost = 0.0; + int nzcnt_estimation = Controller()->RMP()->NumRows(); + if (x.nzcnt_max < nzcnt_estimation) + { + x.nzcnt_max = 2 * nzcnt_estimation; + delete[] x.ind; + delete[] x.val; + x.ind = new int[x.nzcnt_max]; + x.val = new double[x.nzcnt_max]; + } + double lb = 0.0; + double ub = 0.0; + char buffer[256]; + int nzcnt = GetCPXColumn(cost, x.ind, x.val, lb, ub, buffer); + for (int i = 0; i < nzcnt; i++) + { + cost -= dual[x.ind[i]] * x.val[i]; + cout << "Dual: " << dual[x.ind[i]] << " * " << x.val[i] << endl; + } + + return cost; +} const list<int> c_Column::InfoColumnType() const { @@ -3652,19 +3772,16 @@ void c_Column::ChangeObjCoeff( double obj ) } - /*************************************************************************** c_BranchAndBoundConstraint ***************************************************************************/ - c_BranchAndBoundConstraint::c_BranchAndBoundConstraint( c_Controller* controller ) : p_controller( controller ) {} - c_BranchAndBoundConstraint::c_BranchAndBoundConstraint( c_Controller* controller, const list<c_Constraint*>& rmp_constraints ) : p_controller( controller ), _rmp_constraints( rmp_constraints ) @@ -3676,15 +3793,59 @@ c_BranchAndBoundConstraint::~c_BranchAndBoundConstraint() // they are deleted by the controller } - void c_BranchAndBoundConstraint::OutputInOStream( ostream& s ) const { s << "B&B-Constraint: "; list<c_Constraint*>::const_iterator rmp_con; - for ( rmp_con=_rmp_constraints.begin(); rmp_con!=_rmp_constraints.end(); ++rmp_con ) - s << *(*rmp_con) << "\\n"; + for ( auto rmp_con : _rmp_constraints ) + s << *rmp_con /*<< "\n"*/; +} + +/*************************************************************************** + +c_BranchAndBoundConstraintSeparation + +***************************************************************************/ + +inline string c_BranchAndBoundConstraintSeparation::TextInBranching() const +{ + stringstream ss; + if (_rmp_constraints.size() > 0) + { + map<string, int> cut_names; + for (auto con : _rmp_constraints) + cut_names[con->TextInBranching()]++; + bool first = true; + for (auto elem : cut_names) + { + ss << (first ? "" : "+") << elem.first; + first = false; + } + } + else + ss << "cut " << "Unknown"; + return ss.str(); +} + +/*************************************************************************** + +c_StrongBranchingCandidates + +***************************************************************************/ + +c_StrongBranchingCandidates::c_StrongBranchingCandidates(c_Controller* controller) + : c_BranchAndBoundConstraint(controller, list<c_Constraint*>() ) +{} + +c_StrongBranchingCandidates::~c_StrongBranchingCandidates() +{ + // TODO: clean memory } +void c_StrongBranchingCandidates::AddCandidate(const list<list< c_BranchAndBoundConstraint*> >& candidate) +{ + l_candidates.push_back(candidate); +} /*************************************************************************** @@ -3698,7 +3859,6 @@ c_YpsMinusColumn::c_YpsMinusColumn( c_Controller* controller, c_Constraint* my_c constraint( my_constraint ) {} - int c_YpsMinusColumn::GetCPXColumn( double& cost, int* ind, double* val, double& lb, double& ub, char* name ) { cost = -constraint->DeltaMinus(); @@ -3737,14 +3897,12 @@ c_YpsPlusColumn ***************************************************************************/ - // column for slack in stabilized constraint c_YpsPlusColumn::c_YpsPlusColumn( c_Controller* controller, c_Constraint* my_constraint ) : c_Column( controller ), constraint( my_constraint ) {} - int c_YpsPlusColumn::GetCPXColumn( double& cost, int* ind, double* val, double& lb, double& ub, char* name ) { cost = constraint->DeltaPlus(); @@ -3815,7 +3973,8 @@ c_Solution ***************************************************************************/ c_Solution::c_Solution( c_Controller* ctlr, c_RMP*, c_DualVariables* ) - : p_controller( ctlr ) +: p_controller( ctlr ), + obj(numeric_limits<double>::max()) {} c_Solution::~c_Solution() @@ -3834,15 +3993,14 @@ c_PricingProblemHierarchy ***************************************************************************/ c_PricingProblemHierarchy::c_PricingProblemHierarchy( c_Controller* controller, int max_num_failures, int min_num_cols ) - : p_controller( controller ), +: p_controller( controller ), + i_pricing_start_level(0), i_max_num_failures( max_num_failures ), i_min_num_columns( min_num_cols ), i_num_failures(0), - i_pricing_start_level(0), d_best_reduced_costs(0.0) {} - void c_PricingProblemHierarchy::Add( c_PricingProblem* pricer ) { if ( !pricer ) @@ -3854,18 +4012,17 @@ void c_PricingProblemHierarchy::Add( c_PricingProblem* pricer ) throw; } v_hierarchy.push_back( pricer ); + p_controller->IncreasePricingStatisticsVectors(); v_num_failures.push_back(0); pricer->p_my_hierarchy = this; } - c_PricingProblemHierarchy::~c_PricingProblemHierarchy() { for ( int i=0; i<(int)v_hierarchy.size(); i++ ) delete v_hierarchy[i]; } - void c_PricingProblemHierarchy::Update( c_BranchAndBoundNode* node ) { for ( int i=0; i<(int)v_hierarchy.size(); i++ ){ @@ -3884,63 +4041,53 @@ void c_PricingProblemHierarchy::Update( c_DualVariables* dual ) // by all pricers in the hierarchy. } -int c_PricingProblemHierarchy::Go( c_DualVariables* dual, list<c_ColumnPricePair>& cols_and_prices ) +int c_PricingProblemHierarchy::Go( c_DualVariables* dual, list<c_ColumnPricePair>& cols_and_prices, int NumPricersUsed ) { - //NowTODO int num_generated_cols = 0; - d_best_reduced_costs = -10*INFTY; - int start_level = min( i_pricing_start_level, (int)v_hierarchy.size()-1); - for ( int pricing_level= start_level; pricing_level<(int)v_hierarchy.size(); pricing_level++ ) + d_best_reduced_costs = -1000*Controller()->LocalInfty(); + + int counter = 0; + for ( int pricing_level=i_pricing_start_level; pricing_level<(int)v_hierarchy.size() && counter<NumPricersUsed ; pricing_level++,counter++ ) { if( p_controller->TimeOut() ) return 0; - // call pricer + c_TimeInfo timer_pricer_level; + timer_pricer_level.Start(); c_PricingProblem* pricer = v_hierarchy[pricing_level]; pricer->Update(dual); double value = 0.0; int cols_before = (int) cols_and_prices.size(); if ( p_controller->InfoLevel() >= 3 ) p_controller->Info3() << "[" << pricing_level << "] "; - int pricer_return_val=pricer->Go( dual, cols_and_prices, value ); + pricer->Go( dual, cols_and_prices, value ); num_generated_cols += ((int)cols_and_prices.size() - cols_before); - p_controller->Set_number_of_solved_pricing_problems(p_controller->InfoNumberOfSolvedPricingProblems()+1); + timer_pricer_level.Stop(); + p_controller->SetInfoTimePricingLevel(pricing_level, timer_pricer_level.Seconds()); + p_controller->IncreaseInfoCallPricingLevel(pricing_level); + + //cout << p_controller->InfoTimePricingLevel(pricing_level) << endl; + + if(((int)cols_and_prices.size() - cols_before)>0) + p_controller->IncreaseInfoSuccessPricingLevel(pricing_level); // update d_best_reduced_costs if ( pricing_level == (int)v_hierarchy.size()-1 ) { d_best_reduced_costs = value; - if ( num_generated_cols == 0 ) + if ( num_generated_cols == 0 && !p_controller->TimeOut() ) d_best_reduced_costs = 0.0; } - - // update failures and start level - v_num_failures[pricing_level] = (num_generated_cols < i_min_num_columns) ? v_num_failures[pricing_level]+1 : v_num_failures[pricing_level]; //else 0 if consecutive failures are counted, otherwise v_num_failures[pricing_level] - - if( !p_controller->PricingHierarchyStandardOrder() ){ //pricing_level > p_controller->GetMinPricingStartLevel() - - //if not all failed attempts of lower level are used, go back - if( num_generated_cols >= i_min_num_columns && pricing_level > 0 && p_controller->GetNumFailures(pricing_level-1) < i_max_num_failures ){ - i_pricing_start_level = pricing_level - 1; - } - - return num_generated_cols; - } - - if( v_num_failures[pricing_level] >= i_max_num_failures ){ //in not standard order, the pricing level is increased in function Pricing() - i_pricing_start_level = pricing_level+1; - } - // sufficiently many cols generated? + v_num_failures[pricing_level] = (num_generated_cols == 0) ? v_num_failures[pricing_level]+1 : 0; if ( num_generated_cols >= i_min_num_columns) break; - } - //// update failures and start level - //for( int i = 0; i < (int)v_hierarchy.size()-1; i++){ - // if( v_num_failures[i] >= i_max_num_failures ){ - // i_pricing_start_level = i+1; - // } - //} + // update failures and start level + for( int i = 0; i < (int)v_hierarchy.size()-1; i++){ + if( v_num_failures[i] >= i_max_num_failures ){ + i_pricing_start_level = i+1; + } + } return num_generated_cols; } @@ -3951,49 +4098,30 @@ double c_PricingProblemHierarchy::RHSofConvexityConstraint() const /*************************************************************************** -class c_Quadrupel - -***************************************************************************/ - - -bool operator<( const c_Quadrupel& first, const c_Quadrupel& second ) -{ - return (first.i1<second.i1) - || (first.i1==second.i1 && first.i2<second.i2) - || (first.i1==second.i1 && first.i2==second.i2 && first.i3<second.i3 ) - || (first.i1==second.i1 && first.i2==second.i2 && first.i3==second.i3 && first.i4<second.i4); -} - -/*************************************************************************** - class c_IndexInfo ***************************************************************************/ c_IndexInfo::c_IndexInfo( int first_index ) : act_index( first_index ), - indices(), - inverse() + indices() {} - c_IndexInfo::~c_IndexInfo( ) {} - - -int c_IndexInfo::NewIndex ( int i1, int i2, int i3, int i4 ) -{ - indices[ c_Quadrupel( i1, i2, i3, i4 ) ] = act_index; +int c_IndexInfo::NewIndex(int i1, int i2, int i3, int i4) +{ + array<int, 4> key = { i1,i2,i3,i4 }; + indices[key] = act_index; return act_index++; } - -int c_IndexInfo::Index( int i1, int i2, int i3, int i4 ) const +int c_IndexInfo::Index(int i1, int i2, int i3, int i4) const { - map<c_Quadrupel,int>::const_iterator i; - i = indices.find( c_Quadrupel( i1, i2, i3, i4 ) ); - if ( i!=indices.end() ) + array<int, 4> key = { i1,i2,i3,i4 }; + unordered_map<array<int, 4>, int>::const_iterator i = indices.find(key); + if (i != indices.end()) return (*i).second; else { @@ -4004,26 +4132,25 @@ int c_IndexInfo::Index( int i1, int i2, int i3, int i4 ) const bool c_IndexInfo::IndexExists( int i1, int i2, int i3, int i4 ) const { - map<c_Quadrupel,int>::const_iterator i; - i = indices.find( c_Quadrupel( i1, i2, i3, i4 ) ); - if ( i!=indices.end() ) + array<int, 4> key = { i1,i2,i3,i4 }; + unordered_map<array<int, 4>, int>::const_iterator i = indices.find(key); + if (i != indices.end()) return true; else return false; } - -void c_IndexInfo::DeleteIndex( int index_to_del ) +void c_IndexInfo::DeleteIndex(int index_to_del) { - map<c_Quadrupel,int>::iterator i; - for ( i=indices.begin(); i!=indices.end(); ) + unordered_map<array<int, 4>, int>::const_iterator i; + for (i = indices.begin(); i != indices.end(); ) { - if ( (*i).second == index_to_del ) - i = indices.erase( i ); + if ((*i).second == index_to_del) + i = indices.erase(i); else { - if ( (*i).second > index_to_del ) - (*i).second--; + if ((*i).second > index_to_del) + indices[i->first]--; ++i; } } @@ -4042,9 +4169,8 @@ c_ConvergencyInfo::c_ConvergencyInfo() c_ConvergencyInfo::~c_ConvergencyInfo() { - list<c_ConvergencyInfoPoint*>::iterator iter; - for ( iter=l_values.begin(); iter!=l_values.end(); ++iter ) - delete *iter; + for ( auto iter : l_values ) + delete iter; l_values.clear(); } @@ -4067,16 +4193,15 @@ void c_ConvergencyInfo::AddPoint( double time_ms, double lb, double ub ) l_values.push_back( point ); } - void c_ConvergencyInfo::WriteFile( const char* filename ) { std::ofstream out( filename ); list<c_ConvergencyInfoPoint*>::iterator iter; - for ( iter=l_values.begin(); iter!=l_values.end(); ++iter ) + for ( auto iter : l_values) { - out << (*iter)->_time_ms << " "; - out << (*iter)->_lb << " "; - out << (*iter)->_ub << endl; + out << iter->_time_ms << " "; + out << iter->_lb << " "; + out << iter->_ub << endl; } } @@ -4115,15 +4240,12 @@ c_FindKthElement::~c_FindKthElement() delete[] elem; } - - void c_FindKthElement::SetElement( int index, double value ) { assert( (index>=0) && (index<max_size) ); elem[index] = value; } - int compare_less( const void *arg1, const void *arg2 ) { /* Compare all of both strings: */ @@ -4148,7 +4270,6 @@ int compare_greater( const void *arg1, const void *arg2 ) return 0; } - double c_FindKthElement::FindKthElementOfN( int k, int n, order ord ) { assert( (k>=0) && (k<n) && (n<=max_size) ); @@ -4162,7 +4283,6 @@ double c_FindKthElement::FindKthElementOfN( int k, int n, order ord ) } - /*************************************************************************** globale Funktionen @@ -4194,7 +4314,6 @@ double FractionalValue( double x ) return ( ceil(x)-x > x-floor(x) ) ? x-floor(x) : ceil(x)-x; } - double RelativeDifference( double x_100_percent, double value ) { if ( value > x_100_percent ) @@ -4265,4 +4384,5 @@ ostream& operator<<( ostream& s, const c_Constraint& con ) return s; } + } // end of namespace CGBase diff --git a/ColumnGeneration/base.h b/ColumnGeneration/base.h index 29dbfb982f11b91c6dc9cfab501f8f28af02cb2d..bcc72add03e8deaf80a61f5d562f3d0a23856108 100644 --- a/ColumnGeneration/base.h +++ b/ColumnGeneration/base.h @@ -3,6 +3,45 @@ /**************************************************************************** Major Revisions: +- 2023-03-15: New and very clean mechanism for strong branching + * provide: virtual c_StrongBranchingCandidates* ComputeStrongBranchingCandidates(int max_k); + in class derived from c_BranchAndBoundNode + the new class c_StrongBranchingCandidates is a container for the candidate branching decisions + * deprecated: GetkthStrongBranchingCandidate + DecideonBranching + * todo: we may integrate more parameters into the class c_StrongBranchingCandidates to + control the way in which the RMPs are evaluated (number of iterations, which type of pricing algorithms etc.) + +- 2021-09-23: *SolveAsMip have been heavily modified such that preprocessing and dynamic search are used in CPLEX. + Therefore, the problem is written in file "MIP.lp" and read-in again with another instance of CPLEX and the Environment + *If you want to use the SolveAsMip Functionality, you have to implement the function 'virtual c_Solution* CreateSolutionFromMIP(c_RMP* theRMP, vector<int,double>& columns_in_MIPsol)' + in your problem-specific controller. + *To do so, the parameter columns_in_MIPsol contains the column_numbers and the value of all columns with positive value in the solution found by the MIP. + * If the solution is infeasible (e.g., due to containing dummy columns), the function CreateSolutionFromMIP must return nullptr! + * Parameter MIPScreenOutput controls the screen output of cplex for the MIP + * The MIP .lp-file has now a randomized name and this .lp-file is deleted in the destructor of the MIP with 'remove ' + +- 2021-03-30: Unordered_map instead of map in c_IndexInfo. + +- 2019-02-05: Strong Branching Changed - All projects using Strong Branching must be recoded: + + *c_Controller has a vector v_undecidedBranching_nodes that saves completed node for which the branching candidate is not determined yet + info_timer_StrongBranching counts the time that strong branching takes + + * c_BranchAndBoundNodeMVBP has two new functions that needs to be implemented + * and the derived class needs a member that saves the branching candidates + * (e.g., a vector<pair<init,int>> containing possible arc branching decisions) + 1) ComputeStrongBranchingCandidates -> in that function the vector with candidates need to be filled + 2) GetkthStrongBranchingCandidate -> returns the k-th candidate for branching + + *The new function DecideonBranching is called in SelectNextNode to decide on the branching candidates as late as possible + + + +- 2017-06-07: + * INFTY and CG_EPS are now member of the class controller to deal with Nested CG + * New Function "isFractional" also as a member of the controller + - 2015-02-27: Memberfunctions of c_Column: @@ -32,9 +71,6 @@ Major Revisions: ****************************************************************************/ -// turn of message that identifier (in template) get too long -#pragma warning(disable: 4786 ) - #include <math.h> #include <list> @@ -42,86 +78,96 @@ Major Revisions: #include <vector> #include <iostream> #include <list> +#include <limits> +#include <unordered_map> +#include <array> using namespace std; #include "cplexcpp.h" -#include "../SharedFiles/timeinfo.h" +//#include "../SharedFiles/timeinfo.h" +#include "../SharedFiles/stopwatch.h" #include "../SharedFiles/settings_reader.h" +#include <assert.h> + +typedef c_Stopwatch c_TimeInfo; + namespace CGBase { /******************* constants *********************************************/ const double INFTY = 10E8; -const double CG_EPS = 1e-3; +const double CG_EPS = 1e-4; const int NO_ITERATION = -1; - +const int MaxNumPricers = 100000; /******************* class declarations ************************************/ -class c_Controller; +class c_BranchAndBoundNode; +class c_BranchAndBoundConstraint; class c_Column; class c_ColumnPool; +class c_Constraint; +class c_Controller; +class c_ConvergencyInfo; class c_DualVariables; -class c_RMP; -class c_BranchAndBoundNode; -class c_BranchAndBoundConstraint; +class c_FindKthElement; +class c_IndexInfo; class c_PricingProblemHierarchy; class c_PricingProblem; +class c_ReducedCostsFixing; +class c_RMP; class c_SeparationProblem; -class c_Constraint; +class c_Solution; +class c_StrongBranchingCandidates; class c_YpsPlusColumn; class c_YpsMinusColumn; -class c_Solution; -class c_ReducedCostsFixing; -class c_ConvergencyInfo; -class c_IndexInfo; -class c_FindKthElement; -/****************** classes ***********************************************/ +/****************** typedefs ***********************************************/ typedef pair<c_Column*,double> c_ColumnPricePair; +/****************** classes ************************************************/ class c_Controller : public c_SettingsReader { public: enum e_phase_in_bap { initialization, root_node, tree }; private: CPLEXEnvironment* env; - - c_RMP* p_rmp; c_ColumnPool* p_pool; c_ReducedCostsFixing* p_red_costs_fixing; - vector<c_PricingProblemHierarchy*> pricing_problem_hierarchy; - list<c_SeparationProblem*> separation_problems; - vector<c_Column*> cols_in_rmp; + c_BranchAndBoundNode* p_currentNode; + c_ConvergencyInfo* p_conv_info; int next_node_id; - list<c_BranchAndBoundNode*> active_nodes; - list<c_BranchAndBoundNode*> inactive_nodes; - list<c_Constraint*> l_constraints; - c_ConvergencyInfo* p_conv_info; - map<list<int>,int> o_info_statistics_columns; - map<int,string> s_info_column_type_names; - double ub; double aspiration_lb; double aspiration_ub; bool last_node_terminated; bool last_node_integer; - // Info - double first_integer_solution; + map<list<int>, int> o_info_statistics_columns; + map<int, string> s_info_column_type_names; + list<c_BranchAndBoundNode*> active_nodes; + list<c_BranchAndBoundNode*> inactive_nodes; + vector<c_Column*> cols_in_rmp; +protected: + c_BranchAndBoundNode* p_root_node; + c_RMP* p_rmp; c_Solution* p_best_solution; + // Info int info_number_of_rmp_iterations; int info_number_of_generated_columns; int info_number_of_solved_pricing_problems; int info_number_of_generated_cuts; int info_number_of_solved_nodes; + + double first_integer_solution; + double info_time_first_integer_solution; + double info_time_best_integer_solution; + double info_time_overall; double info_time_pricing; double info_time_separation; double info_time_reoptimization; - double info_time_first_integer_solution; - double info_time_best_integer_solution; double info_time_root_lb1; double info_time_root_lb2; double info_time_root_branching; @@ -138,33 +184,32 @@ private: double info_timer_GetDual; double info_timer_RMPUpdate; double info_timer_PricingUpdate; + double info_timer_StrongBranching; + double info_timer_last_PP_root; ostream* info_1; ostream* info_2; ostream* info_3; ostream* info_err; ostream* info_statistics; - c_TimeInfo timer; - c_TimeInfo GetBranchingInfo; e_phase_in_bap i_phase_in_bap; - - double info_timer_last_PP_root; - int i_startHierarchyId; + vector<double> v_info_time_pricing_level; + vector<int> v_info_calls_pricing_level; + vector<int> v_info_success_pricing_level; + c_TimeInfo timer; + c_TimeInfo GetBranchingInfo; protected: - c_BranchAndBoundNode* p_root_node; + list<c_Constraint*> l_constraints; + list<c_PricingProblemHierarchy*> pricing_problem_hierarchy; + list<c_SeparationProblem*> separation_problems; bool LastNodeInteger() const { return last_node_integer; } bool LastNodeTerminated() const { return last_node_terminated; } list<c_SeparationProblem*> feasibility_cutting_problems; // setter to overwrite pricing method - + void Set_number_of_solved_pricing_problems(int new_number){info_number_of_solved_pricing_problems=new_number;} void Set_number_of_generated_columns(int new_number){info_number_of_generated_columns=new_number;} void Set_number_of_rmp_iterations(int new_number){info_number_of_rmp_iterations=new_number;} - void Set_number_of_rmp_update(double new_number){info_timer_RMPUpdate=new_number;} void Set_timer_pricing_update(double new_number){info_timer_PricingUpdate=new_number;} - // getter to overwrite virtuals - vector<c_PricingProblemHierarchy*> Get_pricing_problem_hierarchy(){return pricing_problem_hierarchy;} - vector<int> pricing_levels_in_hierarchies; - vector<int> v_numFailures; public: c_Controller( string settings_file = "settings.txt" ); virtual ~c_Controller(); @@ -187,12 +232,10 @@ public: void AddNode( c_BranchAndBoundNode* node ); void RemoveNode( c_BranchAndBoundNode* node ); - void DeleteTree(); void AddConstraint( c_Constraint* con ) { l_constraints.push_back( con ); } - void AddPricingProblem( c_PricingProblem* prob ); void AddPricingProblemHierarchy( c_PricingProblemHierarchy* hierarchy ); - vector<c_PricingProblemHierarchy*>& PricingProblemHierarchies() { return pricing_problem_hierarchy;} + list<c_PricingProblemHierarchy*>& PricingProblemHierarchies() { return pricing_problem_hierarchy;} void AddSeparationProblem( c_SeparationProblem* ); const list<c_SeparationProblem*>& SeparationProblems() const { return separation_problems; } @@ -200,7 +243,12 @@ public: void AddFeasibilityCuttingProblem( c_SeparationProblem* ); const list<c_SeparationProblem*>& FeasibilityCuttingProblems() const { return feasibility_cutting_problems; } - int NumBranchAndBoundNodes() const { return (int)active_nodes.size() + (int)inactive_nodes.size(); } + int NumBranchAndBoundNodes() const; + int NumOpenBranchAndBoundNodes() const; + int NumOpenStrongBranchingNodes() const; + + const list<c_BranchAndBoundNode*>& GetInactiveNodes() const { return inactive_nodes; } + void TerminateNode( c_BranchAndBoundNode* node ); enum e_branching { depth_first, best_first, fifo }; void StatisticsBranching( map<string,int>& counter_for_branching ); @@ -211,11 +259,13 @@ public: c_StringParameter UseLPOptimizer; c_BoolParameter ReloadSettings; c_BoolParameter Stabilization; + c_IntParameter InfoLevel; c_IntParameter InfoStatisticsLevel; c_BoolParameter InfoDottyOutput; c_BoolParameter InfoGnuplotOutput; c_BoolParameter InfoGraphMLOutput; + c_BoolParameter InfoTikzOutput; c_BoolParameter InfoDottySuppressIntermediateNodes; c_StringParameter InfoInstance; c_BoolParameter InfoPrintOverallBestSolution; @@ -223,12 +273,12 @@ public: c_StringParameter InfoOutputFileName; c_BoolParameter InfoConvergency; c_BoolParameter InfoStatisticsColumns; + c_DoubleParameter MaxSolutionTimeSec; c_IntParameter MaxColumnsToAddInPricing; c_BoolParameter PricingWithColumnSelection; c_IntParameter PricingHierarchyMaxNumFailuresBeforeSwitchOff; c_IntParameter PricingHierarchyMinNumColsToGenerate; - c_BoolParameter PricingHierarchyStandardOrder; c_DoubleParameter MaxRatioColsToRows; c_DoubleParameter MaxRatioColsToRowsMIP; c_BoolParameter Branching; @@ -239,17 +289,31 @@ public: c_IntParameter StrongBranchingNumCandidates; c_IntParameter StrongBranchingMinNumCandidates; c_DoubleParameter StrongBranchingNumCandidatesDecreasePerLevel; - c_BoolParameter StrongBranchingSolveExact; + c_BoolParameter CuttingPlanes; + c_BoolParameter UseMIPSolverForUBs; c_IntParameter UseMIPSolverUptoLevel; + c_IntParameter UseMIPSolverUpFromLevel; c_IntParameter MIPMaxSolutionTimeSec; + c_BoolParameter MIPScreenOutput; + c_BoolParameter MIPAccelerate; // Warning: This still causes a segmentation fault on Mogon II + c_BoolParameter SolveAsMipAfterTimeOut; + c_BoolParameter UseLagrangeanLB; c_BoolParameter ReducedCostVariableFixing; c_IntParameter StabFrequence; c_IntParameter FeasyCutFrequence; c_DoubleParameter FeasyCutGap; c_IntParameter NumThreadsRMP; + c_StringParameter NameDottyOutput; + c_BoolParameter CheckForInftyAndLBError; + + c_StringParameter StrongBranchingRule; + c_DoubleParameter StrongBranchingMu; + c_IntParameter PricersUsedForStrongBranching; + c_DoubleParameter LocalCG_EPS; + c_DoubleParameter LocalInfty; // setters for parameters void SetDEBUG( bool val ) { DEBUG.SetValue( val ); } void SetUseLPOptimizer( const char* txt) { UseLPOptimizer.SetValue( txt ); } @@ -260,7 +324,7 @@ public: void SetInfoDottyOutput(bool flag) { InfoDottyOutput.SetValue( flag ); } void SetInfoGnuplotOutput(bool flag) { InfoGnuplotOutput.SetValue( flag ); } void SetInfoGraphMLOutput(bool flag) { InfoGraphMLOutput.SetValue( flag ); } - void SetInfoDottySuppressIntermediateNodes(bool f) { InfoDottySuppressIntermediateNodes.SetValue(f); } + void SetInfoDottySuppressIntermediateNodes(bool f) { InfoDottySuppressIntermediateNodes.SetValue(f); } void SetInfoPrintOverallBestSolution(bool f) { InfoPrintOverallBestSolution.SetValue(f); } void SetInfoPrintBestSolution(bool f) { InfoPrintBestSolution.SetValue(f); } void SetInfoOutputFileName(const char* txt) { InfoOutputFileName.SetValue( txt ); } @@ -303,7 +367,6 @@ public: double InfoTimePricing() const { return info_time_pricing; } double InfoTimeSeparation() const { return info_time_separation; } double InfoTimeReoptimization() const { return info_time_reoptimization; } - int InfoNumberOfSimplexIterations() const; int InfoNumberOfRMPIterations() const { return info_number_of_rmp_iterations; } int InfoNumberOfGeneratedColumns() const { return info_number_of_generated_columns; } @@ -312,6 +375,9 @@ public: int InfoNumberOfSolvedNodes() const { return info_number_of_solved_nodes; } double InfoTimeFirstIntegerSolution() const { return info_time_first_integer_solution; } double InfoTimeBestIntegerSolution() const {return info_time_best_integer_solution; } + double InfoTimePricingLevel(int level) { double x=0.0; (level < (int)v_info_time_pricing_level.size()) ? x=v_info_time_pricing_level[level] : x=0.0; return x; } + double InfoCallPricingLevel(int level) { int x = 0; (level < (int)v_info_calls_pricing_level.size()) ? x = v_info_calls_pricing_level[level] : x = 0; return x; } + double InfoSuccessPricingLevel(int level) { int x = 0; (level < (int)v_info_success_pricing_level.size()) ? x = v_info_success_pricing_level[level] : x = 0; return x; } double InfoTimeRootLB1() const { return info_time_root_lb1; } double InfoTimeRootLB2() const { return info_time_root_lb2; } double FirstIntegerSolution() const { return first_integer_solution; } @@ -329,6 +395,9 @@ public: double InfoTimeGetDual() const {return info_timer_GetDual; } double InfoTimeRMPUpdate() const {return info_timer_RMPUpdate; } double InfoTimePricingUpdate() const {return info_timer_PricingUpdate; } + double InfoTimeStrongBranching() const { return info_timer_StrongBranching; } + void SetTimeStrongBranching(double new_number) { info_timer_StrongBranching = new_number; } + void RecordStatisticsColumns( const list<int>& col_info_type ); void SetStatisticsColumnTypeName( int i, string s ) { s_info_column_type_names[i] = s; } e_phase_in_bap PhaseInBaP() const; @@ -344,6 +413,7 @@ public: bool SolutionIsFound() const { return ( p_best_solution != NULL ); } double SolutionTimeSec() const { return timer.Seconds(); } bool TimeOut() const { return ( timer.Seconds() > MaxSolutionTimeSec() ); } + c_BranchAndBoundNode* GetCurrentNode() const { return p_currentNode; } void SetInfo1( ostream& out ) { info_1 = &out; } void SetInfo2( ostream& out ) { info_2 = &out; } @@ -360,19 +430,12 @@ public: void SetInfoTimeOrgPool(double i) { info_timer_OrgPool = i; } void SetInfoTimeGetDual(double i) { info_timer_GetDual = i; } void SetInfoTimeMIPForUB(double i) { info_time_MIPForUB = i; } - - void Set_number_of_solved_pricing_problems(int new_number){info_number_of_solved_pricing_problems=new_number;} - - int GetMinPricingStartLevel(){ - int minPricingLevel = pricing_levels_in_hierarchies[0]; - for (auto l= pricing_levels_in_hierarchies.begin();l!= pricing_levels_in_hierarchies.end();l++){ - if( *l < minPricingLevel ){ - minPricingLevel = *l; - } - } - return minPricingLevel; - //return *(std::min_element(pricing_levels_in_hierarchies.begin(), pricing_levels_in_hierarchies.end()) ); - } + void IncreasePricingStatisticsVectors() { + v_info_time_pricing_level.push_back(0.0); v_info_calls_pricing_level.push_back(0); v_info_success_pricing_level.push_back(0); } + void SetInfoTimePricingLevel(int level, double i) { v_info_time_pricing_level[level] += i; } + void IncreaseInfoCallPricingLevel(int level) { v_info_calls_pricing_level[level] ++; } + void IncreaseInfoSuccessPricingLevel(int level) { v_info_success_pricing_level[level] ++; } + bool local_IsFractional(double x); // we may overwrite these methods virtual void Go( void ); @@ -380,34 +443,26 @@ public: virtual c_BranchAndBoundNode* SelectNextNode(); virtual void OutputInOStream( ostream& s ) const; virtual void Update( c_BranchAndBoundNode* node ); - virtual void Pricing( int iteration, c_DualVariables* dual, list<c_Column*>& col_list ); + virtual void Pricing(int iteration, c_DualVariables* dual, list<c_Column*>& col_list, int NumPricersUsed = MaxNumPricers); virtual void ColumnSelectionInPricing( list<c_ColumnPricePair>& cols_and_prices, list<c_ColumnPricePair>& others ) { /*do nothing*/ }; virtual void OrganizePool(); - virtual c_ReducedCostsFixing* CreateReducedCostsFixing() { return NULL; } + virtual c_ReducedCostsFixing* CreateReducedCostsFixing() { return nullptr; } virtual void Separation( c_RMP* rmp, list<c_Constraint*>& cuts, int BaB_node_level ); virtual void PerformBranching( c_BranchAndBoundNode* node ); virtual void UpdateAllNodeInformations( c_BranchAndBoundNode* selected_node ); virtual bool FeasibilityCutting(); virtual void FixPartOfSolution(); virtual void UnfixSolution(); + virtual void DeleteTree(); // these methods must be provided by derived class virtual c_RMP* CreateNewRMP() = 0; virtual c_BranchAndBoundNode* CreateRootNode() = 0; virtual c_DualVariables* CreateDualVariables() = 0; virtual c_Solution* CreateSolution( c_RMP*, c_DualVariables* ) = 0; + virtual c_Solution* CreateSolutionFromMIP(c_RMP* theRMP, vector<pair<int, double> >& columns_in_Mip_sol); virtual void AddPricingProblems() = 0; virtual void AddSeparationProblems() = 0; - - int GetStartHierarchy(){return i_startHierarchyId;} - void SetStartHierarchy(int hierarchyId){i_startHierarchyId = hierarchyId;} - - int GetNumFailures(int pricingLevel){return v_numFailures[pricingLevel];} - void IncreaseNumFailures(int pricingLevel){ v_numFailures[pricingLevel]++;} - void ResetNumFailures(){ - for( int i = 0; i < v_numFailures.size(); i++) - v_numFailures[i] = 0; - } }; @@ -427,8 +482,8 @@ public: bool IsInRMP() const { return is_in_rmp; } long NumInRMP() const { return num_in_rmp; } double ReducedCosts( const c_DualVariables& dual ); - void ChangeObjCoeff( double obj ); - // double PrintReducedCosts( const c_DualVariables& dual ); + double CoutReducedCosts(const c_DualVariables& dual); + virtual void ChangeObjCoeff( double obj ); // we may overwrite these methods virtual void OutputInOStream( ostream& s ) const; @@ -440,12 +495,10 @@ public: virtual bool IsSystemColumn() const { return false; }; // CPLEX format of columns virtual double GetCPXObj() = 0; - //virtual int CPXnzcnt() = 0; virtual int GetCPXColumn( double& cost, int* ind, double* val, double& lb, double& ub, char* name ) = 0; }; - class c_DualVariables { c_Controller* p_controller; @@ -456,6 +509,7 @@ public: c_Controller* Controller() const { return p_controller; } double* CPXDual() { return dual_val; } + double* CPXDual() const { return dual_val; } double& operator[](int idx) { return dual_val[idx]; } double operator[](int idx) const { return dual_val[idx]; } @@ -524,28 +578,32 @@ public: virtual double DefaultLowerBound(); virtual double DefaultUpperBound(); - virtual void OutputInOStream( ostream& s ) const; + virtual void OutputInOStream( ostream& s ) const; }; -class c_Quadrupel { +class hasher_fct { public: - int i1, i2, i3, i4; - c_Quadrupel( int v1=0, int v2=0, int v3=0, int v4=0 ) - : i1(v1), i2(v2), i3(v3), i4(v4) {} -}; + int operator()(const array<int, 4>& a) const + { + int h = 0; + for (int i = 0; i < 4; ++i) + { + h = h * 31 + a[i]; + } + return h; + } +}; -bool operator<( const c_Quadrupel& first, const c_Quadrupel& second ); class c_IndexInfo { int act_index; - map<c_Quadrupel,int> indices; - map<int,c_Quadrupel> inverse; + unordered_map<array<int, 4>, int, hasher_fct> indices; public: c_IndexInfo( int first_index = 0 ); virtual ~c_IndexInfo( void ); - typedef enum { NO_INDEX = -1 }; + // typedef enum { NO_INDEX = -1 }; int NewIndex( int i1, int i2=0, int i3=0, int i4=0 ); int Index( int i1, int i2=0, int i3=0, int i4=0 ) const; bool IndexExists( int i1, int i2=0, int i3=0, int i4=0 ) const; @@ -554,7 +612,6 @@ public: }; - class c_Constraint { friend class c_RMP; c_Controller* p_controller; @@ -562,10 +619,9 @@ friend class c_RMP; int _sense; double _rhs; double _estimator_dual; - c_BranchAndBoundNode* active_at_node; - c_YpsPlusColumn* col_yps_plus; c_YpsMinusColumn* col_yps_minus; + c_BranchAndBoundNode* active_at_node; protected: // use only in the construction phase of derived classes void SetRhs( double rhs ) { _rhs = rhs; } @@ -592,46 +648,36 @@ public: // we may overwrite these methods virtual bool IsGloballyValid() const { return false; } // when added in B&B-tree - + virtual int Type() const { return _type; } virtual bool Stabilization() const { return false; } virtual c_YpsPlusColumn* CreateYpsPlusColumn() { return NULL; } virtual c_YpsMinusColumn* CreateYpsMinusColumn() { return NULL; } - virtual bool UpdateStabilization( int iteration, int num_cols_pricing, double dual ) { return true; } virtual double EpsMinus() const { return 1.0; } virtual double EpsPlus() const { return 1.0; } virtual double DeltaMinus() const { return 1.0; } virtual double DeltaPlus() const { return 1.0; } - + virtual string TextInBranching() const { return string("constraint"); } virtual void OutputInOStream( ostream& s ) const; }; -class c_GlobalValidConstraint : public c_Constraint{ -public: - c_GlobalValidConstraint(c_Controller* controller, char sense, double rhs, double estimator_dual, int type, int n1=0, int n2=0, int n3=0 ) - : c_Constraint( controller, sense, rhs, estimator_dual, type, n1, n2, n3) {}; - virtual ~c_GlobalValidConstraint() {}; - virtual bool IsGloballyValid() const { return true;} -}; - class c_RMP: public CPLEXProblem, private c_IndexInfo { c_Controller* p_controller; - +protected: + int actual_iteration; double actual_obj_value; double d_mip_objval; -protected: - int actual_iteration; vector< c_Constraint* > constraints; int info_number_of_simplex_iterations; - int num_removed_constraints; int num_added_constraints; - + string help_mipname; public: c_RMP( c_Controller* controller ); virtual ~c_RMP(); - void Reset(); + + void Reset(); c_Controller* Controller() const { return p_controller; } @@ -659,9 +705,11 @@ public: bool ConstraintIndexExists( int i1, int i2=0, int i3=0, int i4=0 ) const; int NumConstraints( void ) const; c_Constraint* Constraint( int i1, int i2=0, int i3=0, int i4=0 ) const; - void ConstructEmptyRMP( CPLEXProblem::ObjSense objsen = CPLEXProblem::MINIMIZE ); const vector<c_Constraint*> Constraints() const { return constraints; } + int GetNumRemCon() const { return num_removed_constraints; } + int GetNumAddCon() const { return num_added_constraints; } + void SetObjval(double the_val) { actual_obj_value = the_val; } // we may overwrite these methods virtual void Optimize( e_solver solver ); @@ -669,119 +717,117 @@ public: virtual bool IsIntegral() const; virtual bool IsFeasible() const; virtual bool UpdateStabilization( int iteration, int num_cols_pricing, c_DualVariables* dual ); - virtual double RoundLBTo( double lb ) { return ceil( lb - CG_EPS * CG_EPS ); } + virtual double RoundLBTo(double lb) { return ceil(lb - Controller()->LocalCG_EPS() * Controller()->LocalCG_EPS()); } virtual void RemoveAllBranchAndBoundConstraints(); virtual void Update( c_BranchAndBoundNode* node ); virtual void OutputInOStream( ostream& s ) const; virtual void AddBranchAndBoundConstraints( const list<c_BranchAndBoundConstraint*>& constraint_list ); virtual void SetRequiredVarsIntegral(); + virtual void SetConstraintSense(int contraintIndex, char sense); + // these methods must be provided by derived class virtual void InitializeFirstRMP( ) = 0; - - int GetNumRemCon() { return num_removed_constraints; } - int GetNumAddCon() { return num_added_constraints; } }; - class c_BranchAndBoundNode { c_Controller* p_controller; - int node_number; - int iteration; - int node_number_solved; - list<c_BranchAndBoundConstraint*> added_constraints; - bool is_removed; - bool is_integral; - bool is_feasible; - bool terminated; - double lb; - double rmp_value; - c_BranchAndBoundNode* father; - list<c_BranchAndBoundNode*> sons; + c_BranchAndBoundNode* p_father; + int i_node_number; + int i_iteration; + int i_node_number_solved; + int i_num_iter_wo_feasibility_cut; + bool b_is_removed; + bool b_is_integral; + bool b_is_feasible; + bool b_terminated; + double d_LB; + double d_RMP_value; + list<c_BranchAndBoundConstraint*> l_added_constraints; + list<c_BranchAndBoundNode*> l_sons; // private methods - void SetRMPValue( double val ) { rmp_value = val; } + void SetRMPValue( double val ) { d_RMP_value = val; } double ComputeLagrangeanLB(); - int NumIterWithoutFeasyCut; -protected: - void SetIsIntegral() { is_integral = true; } + + int i_level; public: c_BranchAndBoundNode( c_Controller* controller ); c_BranchAndBoundNode( c_BranchAndBoundNode* p_father, list<c_BranchAndBoundConstraint*>& add_constraints ); virtual ~c_BranchAndBoundNode(); - + // getter c_Controller* Controller() const { return p_controller; } - c_DualVariables* Go( double termination_lb = numeric_limits<double>::infinity(), double termination_ub = -numeric_limits<double>::infinity() ); // returns the dual optimal solution - // we may overwrite these methods - int NodeNumber() const { return node_number; } - void SetNodeNumber( int num ) { node_number = num; } - int NodeNumberSolved() const { return node_number_solved; } - void SetNodeNumberSolved( int num ) { node_number_solved = num; } - c_BranchAndBoundNode* Father() const { return father; } - c_BranchAndBoundNode* Brother() const; - const list<c_BranchAndBoundNode*>& Sons() const { return sons; } + int NodeNumber() const { return i_node_number; } + int NodeNumberSolved() const { return i_node_number_solved; } + c_BranchAndBoundNode* Father() const { return p_father; } + const list<c_BranchAndBoundNode*>& Sons() const { return l_sons; } string UniquePattern() const; - - int NumSons() const { return (int)sons.size(); } - bool IsRootNode() const { return ( father == NULL ) || ( father->NumSons() == 1 && father->IsRootNode() ); } + int NumSons() const { return (int)l_sons.size(); } + bool IsRootNode() const { return ( p_father == NULL ) || ( p_father->NumSons() == 1 && p_father->IsRootNode() ); } int Level() const; int DetailedLevel() const; - double LB() const { return lb; } - void SetLB( double val ); + int GetLevel_new() const; + double LB() const { return d_LB; } double SubtreeLB() const; - double RMPValue() const { return rmp_value; } - void Terminate() { p_controller->TerminateNode( this ); terminated = true; } - bool Terminated() const { return terminated; } - bool IsSolved() const { return (iteration > NO_ITERATION); } - void SetUnsolved() { iteration = NO_ITERATION; } - bool IsRemoved() const { return is_removed; } - void SetRemoved( bool removed ) { is_removed = removed; } + double RMPValue() const { return d_RMP_value; } + bool Terminated() const { return b_terminated; } + bool IsSolved() const { return (i_iteration > NO_ITERATION); } + bool IsRemoved() const { return b_is_removed; } bool AnySonSolved() const; - void AddSon( c_BranchAndBoundNode* son ); - virtual void GetConstraints( list<c_BranchAndBoundConstraint*>& ) const; //TODO: virtual added by Anka to ignore Aggregated SRCuts - const list<c_BranchAndBoundConstraint*>& AddedConstraints() const { return added_constraints; } + void GetConstraints( list<c_BranchAndBoundConstraint*>& ) const; + const list<c_BranchAndBoundConstraint*>& AddedConstraints() const { return l_added_constraints; } + // setters + void SetNodeNumber(int num) { i_node_number = num; } + void SetNodeNumberSolved(int num) { i_node_number_solved = num; } + void SetIsIntegral(bool val) { b_is_integral = val; } + void SetLB(double val); + void SetRemoved(bool removed) { b_is_removed = removed; } + void SetUnsolved() { i_iteration = NO_ITERATION; } + // non-const methods + c_DualVariables* Go(double termination_lb = numeric_limits<double>::infinity(), double termination_ub = -numeric_limits<double>::infinity(), int NumPricersUsed = MaxNumPricers); // returns the dual optimal solution + void Terminate() { p_controller->TerminateNode(this); b_terminated = true; } + void AddSon(c_BranchAndBoundNode* son); // we may overwrite these methods - //virtual void ComputeBounds( c_Column* col, double& lb, double& ub ); + virtual void WriteFile(); virtual double SolveAsMIP( long time_limit ); virtual bool IsIntegral() const; virtual bool IsFeasible() const; virtual void OutputInOStream( ostream& s ) const; virtual bool EarlyBranching() { return false; } virtual bool GetSeparationConstraints( list<c_BranchAndBoundConstraint*>& con_list ); - virtual void Get_k_th_best_BranchAndBoundContraints( int k, list< list<c_BranchAndBoundConstraint*> >& con_list ); - + virtual void DecideOnStrongBranchingCandidate(); // these methods must be provided by derived class - virtual c_BranchAndBoundNode* - CreateNewNode( list<c_BranchAndBoundConstraint*>& new_constraint_list )= 0; + virtual c_BranchAndBoundNode* CreateNewNode( list<c_BranchAndBoundConstraint*>& new_constraint_list )= 0; virtual void GetBranchAndBoundConstraints( list< list<c_BranchAndBoundConstraint*> >& con_list ) = 0; - - virtual void WriteFile(); + // this method must be provided by derived class if SB is used + virtual c_StrongBranchingCandidates* ComputeStrongBranchingCandidates(int max_k); }; - class c_BranchAndBoundConstraint { +protected: c_Controller* p_controller; list<c_Constraint*> _rmp_constraints; public: - c_BranchAndBoundConstraint( c_Controller* controller ); + c_BranchAndBoundConstraint( c_Controller* controller ); c_BranchAndBoundConstraint( c_Controller* controller, const list<c_Constraint*>& rmp_constraints ); virtual ~c_BranchAndBoundConstraint() ; c_Controller* Controller() const { return p_controller; } const list<c_Constraint*>& RMPConstraints() const { return _rmp_constraints; } - void AddConstraint( c_Constraint* con ) { _rmp_constraints.push_back(con); } - virtual bool IsSeparationConstraint() const { return false; } - + void AddConstraint( c_Constraint* con ) { _rmp_constraints.push_back(con); } + // we may overwrite these methods - virtual void OutputInOStream( ostream& s ) const; - virtual void ComputeBounds( c_Column* col, double& lb, double& ub ) {} // mache nichts + virtual bool IsSeparationConstraint() const { return false; } + virtual bool IsStrongBranchingCandidates() const { return false; } + virtual void ComputeBounds( c_Column* col, double& lb, double& ub ) {} // do nothing + // output + virtual void OutputInOStream(ostream& s) const; virtual int GnuPlotLineStyle() { return 12; } virtual string TextInBranching() const { return string("branch"); } - virtual void Serialize( ostream& ) = 0; - virtual void Serialize( istream& ) = 0; }; + class c_BranchAndBoundConstraintSeparation : public c_BranchAndBoundConstraint { public: c_BranchAndBoundConstraintSeparation( c_Controller* controller, list<c_Constraint*>& rmp_constraints ) @@ -789,9 +835,22 @@ public: virtual ~c_BranchAndBoundConstraintSeparation() {} virtual bool IsSeparationConstraint() const { return true; } - virtual void Serialize( ostream& ) {}; - virtual void Serialize( istream& ) {}; - virtual string TextInBranching() const { return string("cut"); } + virtual string TextInBranching() const; +}; + + +class c_StrongBranchingCandidates : public c_BranchAndBoundConstraint { + vector<list<list< c_BranchAndBoundConstraint*> > > l_candidates; +public: + c_StrongBranchingCandidates(c_Controller* controller); + virtual ~c_StrongBranchingCandidates(); + + void AddCandidate(const list<list< c_BranchAndBoundConstraint*> >& candidate); + bool IsStrongBranchingCandidates() const { return true; } + auto Candidates() const { return l_candidates; } + // we may overwrite these methods + virtual int EstimatedNumBranchAndBoundNodes() const { return 2; } + string TextInBranching() const { return string("strong br"); } }; /* @@ -812,51 +871,44 @@ the parameters have the following meaning: */ class c_PricingProblemHierarchy { c_Controller* p_controller; - double d_best_reduced_costs; +protected: + int i_pricing_start_level; int i_max_num_failures; int i_min_num_columns; - - int i_pricing_start_level; - int i_num_failures; - vector<int> v_num_failures; -protected: - vector<c_PricingProblem*> v_hierarchy; + int i_num_failures; + double d_best_reduced_costs; + vector<int> v_num_failures; + vector<c_PricingProblem*> v_hierarchy; public: c_PricingProblemHierarchy( c_Controller*, int max_num_failures, int min_num_columns ); virtual ~c_PricingProblemHierarchy(); c_Controller* Controller() const { return p_controller; } - void Add( c_PricingProblem* pricer ); double BestReducedCosts() const { return d_best_reduced_costs; } double RHSofConvexityConstraint() const; int GetNumFailures() const {return i_num_failures;} - int GetNumFailures(int level) const {return v_num_failures[level];} - int GetPricingStartLevel() const { return i_pricing_start_level;} - void SetPricingStartLevel(int val) { i_pricing_start_level = val;} - void SetBestReducedCost(double val) { d_best_reduced_costs = val;} - int GetNumPricingLevels() const {return (int)v_hierarchy.size();} - const int GetMinNumCols() const { return i_min_num_columns;} - const int GetMaxNumFailures() const { return i_max_num_failures;} + int Size() const { return (int) v_hierarchy.size(); } + c_PricingProblem* PricingProblem(int i) const { return v_hierarchy[i]; }; // we may overwrite these methods virtual void Update( c_BranchAndBoundNode* node ); virtual void Update( c_DualVariables* dual ); - virtual int Go( c_DualVariables* dual, list<c_ColumnPricePair>& cols_and_prices ); + virtual int Go(c_DualVariables* dual, list<c_ColumnPricePair>& cols_and_prices, int NumPricersUsed = MaxNumPricers); }; class c_PricingProblem { - friend class c_PricingProblemHierarchy; +friend class c_PricingProblemHierarchy; c_Controller* p_controller; - c_PricingProblemHierarchy* p_my_hierarchy; + c_PricingProblemHierarchy* p_my_hierarchy; public: c_PricingProblem( c_Controller* ); virtual ~c_PricingProblem() {}; c_Controller* Controller() const { return p_controller; } c_PricingProblemHierarchy* MyPricingProblemHierarchy() const { return p_my_hierarchy; } - // we may overwrite these methods + // we may overwrite these methods virtual void OutputInOStream( ostream& s ) const; // these methods must be provided by derived class virtual void Update( c_BranchAndBoundNode* node ) = 0; @@ -872,12 +924,9 @@ public: c_SeparationProblem( c_Controller* ); virtual ~c_SeparationProblem() {} - // we may overwrite these methods - c_Controller* Controller() const { return p_controller; } - + c_Controller* Controller() const { return p_controller; } // we may overwrite these methods virtual void OutputInOStream( ostream& s ) const; - // these methods must be provided by derived class virtual int Go( c_RMP* rmp, list<c_Constraint*>& cuts, int node_level ) = 0; }; @@ -947,7 +996,7 @@ public: c_FindKthElement( int _max_size ); ~c_FindKthElement(); - typedef enum order { increasing_order, decreasing_order }; + typedef enum order { increasing_order, decreasing_order } order; void SetElement( int index, double value ); double FindKthElementOfN( int k, int n, order ord=increasing_order ); }; @@ -976,4 +1025,4 @@ ostream& operator<<( ostream& s, const c_Solution& ); ostream& operator<<( ostream& s, const c_Constraint& ); } // end of namespace CGBase -#endif // of BASE_H +#endif // of BASE_H \ No newline at end of file diff --git a/bcp_cvrptw.cpp b/bcp_cvrptw.cpp index b840e9a56d0d6724b2072617579cf5eb635ad9d8..9eaca8ae38a1795449fb1f0dc03d7f04c28df756 100644 --- a/bcp_cvrptw.cpp +++ b/bcp_cvrptw.cpp @@ -386,7 +386,6 @@ void c_ControllerCVRPTW::AddPricingProblems() maxNumPricingLevels = hierarchies[d]->GetNumPricingLevels(); AddPricingProblemHierarchy( hierarchies[d] ); } - v_numFailures = vector<int>(maxNumPricingLevels, 0); //setNumEasierPricingProblems(easierPricingProblems); } @@ -395,7 +394,203 @@ void c_ControllerCVRPTW::AddSeparationProblems() AddSeparationProblem(new c_SeparationProblemCVRPTW(this, *this)); } +void c_ControllerCVRPTW::Pricing(int iteration, c_DualVariables* dual, list<c_Column*>& col_list, int NumPricersUsed) +{ + // Initialisierung + list<c_ColumnPricePair> column_price_list; + list<c_ColumnPricePair> to_remove_list; + + if (InfoLevel() >= 3) + Info3() << "Pricing " << flush; + + c_TimeInfo pricing_time_info; + int generated_cols = 0; + vector<c_PricingProblemHierarchy*>& hierarchies = PricingProblemHierarchies(); + + if (PricingHierarchyStandardOrder()) { + vector<c_PricingProblemHierarchy*>::iterator hh; + for (hh = hierarchies.begin(); hh != hierarchies.end(); ++hh) + { + if (TimeOut()) + return; + c_PricingProblemHierarchy* hierarchy = *hh; + hierarchy->Update(dual); + generated_cols += hierarchy->Go(dual, column_price_list); + //info_number_of_solved_pricing_problems++; + } + } + else { + //NowTODO + int startHierarchy = GetStartHierarchy(); + int numHierarchies = (int)hierarchies.size(); + int maxPricingLevel = hierarchies[0]->GetNumPricingLevels() - 1; + for (int h = 0; h < numHierarchies; h++) { + hierarchies[h]->SetBestReducedCost(-10 * INFTY); //Initialize reduced cost to -infinity, so that the lagrangean lower bound is also valid, if not all pricing hierarchies are started + } + + bool continuePricing = true; + vector<bool> v_failed = vector<bool>(numHierarchies, false); + //vector<bool> v_failed = vector<bool>(numHierarchies, false); + do { + if (TimeOut()) + return; + c_PricingProblemHierarchy* hierarchy = hierarchies[startHierarchy]; + hierarchy->Update(dual); + + if (all_of(v_failed.begin(), v_failed.end(), [](bool b) {return b; })) { + IncreaseNumFailures(min(hierarchy->GetPricingStartLevel(), maxPricingLevel)); + if (hierarchy->GetPricingStartLevel() >= maxPricingLevel) //no more improving columns for this node + break; //Other behaviour when MinNumCols > 0? + else { + //increase pricing level of all hierarchies (at least temporarily; if maxNumFailures is not reached, it is decreased again) + int minPricingStartLevel = maxPricingLevel; + for (int p = 0; p < numHierarchies; p++) { + if (hierarchies[p]->GetPricingStartLevel() < minPricingStartLevel) + minPricingStartLevel = hierarchies[p]->GetPricingStartLevel(); + } + for (int p = 0; p < v_failed.size(); p++) { + hierarchies[p]->SetPricingStartLevel(minPricingStartLevel + 1); //hierarchies[p]->GetPricingStartLevel() +1 ); + v_failed[p] = false; + } + } + } + + int added_cols = hierarchy->Go(dual, column_price_list); + generated_cols += added_cols; + + if (added_cols >= hierarchy->GetMinNumCols()) { //stay at the same hierarchy until no more columns for it can be found + continuePricing = false; + SetStartHierarchy(startHierarchy); + } + else { + v_failed[startHierarchy] = true; + startHierarchy = (startHierarchy + 1) % numHierarchies; + } + + } while (continuePricing); + + ////New test: always start with pricing hierarchy which has the minimal pricingStartLevel + //int maxPricingLevel = hierarchies[0]->GetNumPricingLevels(); + //int minPricingStartLevel = maxPricingLevel+1; + //int bestH = -1; + //int numHierarchies = (int)hierarchies.size(); + //for ( int h = 0; h < numHierarchies; h++ ) { + // hierarchies[h]->SetBestReducedCost(-10*INFTY); //Initialize reduced cost to -infinity, so that the lagrangean lower bound is also valid, if not all pricing hierarchies are started + // if( pricing_levels_in_hierarchies[h] < minPricingStartLevel ){ + // minPricingStartLevel = hierarchies[h]->GetPricingStartLevel(); + // bestH = h; + // } + //} + + ////stay at one pricing hierarchy until no more columns can be found. Then switch first to other hierarchies with same pricing level. + //int hNum = bestH; + //bool continuePricing = true; + //vector<bool> hierarchyRunWithMaxPricingLevel = vector<bool>(numHierarchies); + //for( int p = 0; p < numHierarchies; p++) + // hierarchyRunWithMaxPricingLevel[p] = false; + //bool pricingLevelIncreased = false; + //do{ + // if( TimeOut() ) + // return; + // c_PricingProblemHierarchy* hierarchy = hierarchies[hNum]; + // hierarchy->Update( dual ); + // hierarchyRunWithMaxPricingLevel[hNum] = (pricing_levels_in_hierarchies[hNum] >= maxPricingLevel-1)? true : false; + // generated_cols += hierarchy->Go( dual, column_price_list ); + // int curPricingLevel = hierarchy->GetPricingStartLevel(); + // if( curPricingLevel > pricing_levels_in_hierarchies[hNum] ) + // pricingLevelIncreased = true; + // if( generated_cols > hierarchies[hNum]->GetMinNumCols() ){ //stay at the same hierarchy until no more columns for it can be found + // //if the maximal number of failures is not reached and the pricing level below, go back + // if( curPricingLevel > 0 && hierarchy->GetNumFailures(curPricingLevel-1) < hierarchy->GetMaxNumFailures() ){ + // hierarchy->SetPricingStartLevel( curPricingLevel-1 ); + // } + // pricing_levels_in_hierarchies[hNum] = hierarchy->GetPricingStartLevel(); + // break; + // } + // pricing_levels_in_hierarchies[hNum] = curPricingLevel; + // hNum= hNum + 1; + // hNum = hNum % numHierarchies; + + // //if maxNumFailures is greater 1, ensure that the pricing level is although increased if no hierarchy found new columns and the same pricing problem would be resolved + // if( hNum == bestH && generated_cols == 0 && !pricingLevelIncreased){ //GetMinPricingStartLevel() != maxPricingLevel + // //all hierarchies were tested and found no column, but maximal pricing level not reached, such that pricing will continue: + // //increase all minimal pricing levels (by ignoring a maxNumFailures > 1), since nothing has changed + // for ( int h = 0; h < numHierarchies; h++ ) { + // if( hierarchies[h]->GetPricingStartLevel() == GetMinPricingStartLevel() ){ + // hierarchies[h]->SetPricingStartLevel(pricing_levels_in_hierarchies[h]+1); + // pricing_levels_in_hierarchies[h] += 1; + // } + // } + // } + // continuePricing = ( find(begin(hierarchyRunWithMaxPricingLevel), end(hierarchyRunWithMaxPricingLevel), false) != end(hierarchyRunWithMaxPricingLevel) ); // not all true + //} while( continuePricing ); //hNum != bestH || GetMinPricingStartLevel() != maxPricingLevel + } + + // Select an appropriate subset of columns to include into the + // RMP. Especially in situations, where the number of generated columns is + // large this method should take care of a good choice of columns. + int size_bofore_sel = (int)column_price_list.size(); + if (PricingWithColumnSelection()) + { + ColumnSelectionInPricing(column_price_list, to_remove_list); + list<c_ColumnPricePair>::iterator col_price; + for (col_price = to_remove_list.begin(); col_price != to_remove_list.end(); ++col_price) + delete (*col_price).first; + } + pricing_time_info.Stop(); + + //info_time_pricing += pricing_time_info.Seconds(); + info_number_of_generated_columns += generated_cols; + info_number_of_rmp_iterations++; + + if (generated_cols && (InfoLevel() >= 3)) + { + Info3() << " " << generated_cols << " col's generated" << flush; + if (size_bofore_sel != column_price_list.size()) + Info3() << ", " << (int)column_price_list.size() << " col's selected" << flush; + } + if (InfoLevel() >= 3) + Info3() << " (" << pricing_time_info.MilliSeconds() << "ms)" << flush; + + // check, if too many columns generated + double threshold = INFTY * INFTY; + if ((int)column_price_list.size() >= MaxColumnsToAddInPricing()) + { + // take only the best MaxColumnsToAddInPricing() of + c_FindKthElement k_of_n_alg((int)column_price_list.size()); + + list<c_ColumnPricePair>::iterator col_price; + int idx = 0; + for (col_price = column_price_list.begin(); col_price != column_price_list.end(); ++col_price) + k_of_n_alg.SetElement(idx++, (*col_price).second); + threshold = k_of_n_alg.FindKthElementOfN(MaxColumnsToAddInPricing() - 1, (int)column_price_list.size()); + if (InfoLevel() >= 3) + { + int to_add = 0; + for (col_price = column_price_list.begin(); col_price != column_price_list.end(); ++col_price) + if ((*col_price).second <= threshold) + to_add++; + Info3() << ", only " << to_add << " added.\n"; + } + } + else + if (InfoLevel() >= 3) + Info3() << ".\n"; + + // Transfer Columns from pairs to single columns + col_list.clear(); + + list<c_ColumnPricePair>::iterator col_price; + for (col_price = column_price_list.begin(); + col_price != column_price_list.end(); + ++col_price) + if ((*col_price).second <= threshold) + col_list.push_back((*col_price).first); + else + delete (*col_price).first; + column_price_list.clear(); +} void c_ControllerCVRPTW::createREFs() { diff --git a/bcp_cvrptw.h b/bcp_cvrptw.h index 7814a2af1020478221e47e2a38ab5bbdd206baf5..b056882b5297e29d21b7b965f7d1bb6a97558a31 100644 --- a/bcp_cvrptw.h +++ b/bcp_cvrptw.h @@ -47,6 +47,7 @@ public: // these methods can be modified virtual void Go( void ); virtual void WriteSolutionInFile( string FileName ); + virtual void Pricing(int iteration, c_DualVariables* dual, list<c_Column*>& col_list, int NumPricersUsed = MaxNumPricers); // these methods must be provided by derived class virtual c_RMP* CreateNewRMP(); virtual c_BranchAndBoundNode* CreateRootNode(); @@ -378,7 +379,7 @@ public: virtual ~c_Constraint_PVRPTW() {}; virtual void OutputInOStream(ostream& s) const; // = 0; //virtual int Type() = 0; - virtual const int Type() const { return i_type; } + }; @@ -648,7 +649,7 @@ public: const vector<int>& GetSubset(){return v_Subset;} int getCutId() const {return i_cutNumber;} virtual void OutputInOStream( ostream& s ) const; - const int Type() const { return c_RMP_CVRPTW::eSubsetRowCut; } + int Type() const { return c_RMP_CVRPTW::eSubsetRowCut; } }; @@ -664,7 +665,7 @@ public: const vector<int>& GetSubset(int day) const { return v_SubsetsPerDay[day];} int getCutId() const {return i_cutNumber;} virtual void OutputInOStream( ostream& s ) const; - const int Type() const { return c_RMP_CVRPTW::eAggregatedSubsetRowCut; } + int Type() const { return c_RMP_CVRPTW::eAggregatedSubsetRowCut; } }; class c_NgIncreasementCVRPTW : public c_GlobalValidConstraint_PVRPTW{ @@ -678,7 +679,7 @@ public: const vector<int>& GetLocations(){return v_changedLocations;} int getCutId() const {return i_cutNumber;} virtual void OutputInOStream( ostream& s ) const; - const int Type() const { return c_RMP_CVRPTW::eNgIncreasementCut; } + int Type() const { return c_RMP_CVRPTW::eNgIncreasementCut; } }; @@ -697,7 +698,7 @@ public: int GetRHS(){return i_rhs;} int getCutId() const {return i_cutNumber;} virtual void OutputInOStream( ostream& s ) const; - const int Type() const { return c_RMP_CVRPTW::eCapacityCut; } + int Type() const { return c_RMP_CVRPTW::eCapacityCut; } }; @@ -715,7 +716,7 @@ public: // moeglicherweise zu spezialisieren virtual bool IsGloballyValid() const { return true; } virtual void OutputInOStream(ostream& s) const = 0; - const int Type() const { return c_RMP_CVRPTW::eARCSET; } + int Type() const { return c_RMP_CVRPTW::eARCSET; } }; class c_TwoPathConstraint : public c_ArcSetConstraint {