29 #include "scip/cons_linear.h"
44 const int p_num_nodes,
46 const vector< int >& p_demand,
47 const vector< vector<int> >& p_distance,
48 const vector< vector<SCIP_VAR*> >& p_arc_var,
49 const vector< vector<SCIP_CONS*> >& p_arc_con,
50 const vector<SCIP_CONS* >& p_part_con
52 ObjPricer(scip, p_name,
"Finds tour with negative reduced cost.", 0, TRUE),
53 _num_nodes(p_num_nodes),
54 _capacity(p_capacity),
56 _distance(p_distance),
76 for (
int i = 0; i < num_nodes(); ++i)
78 for (
int j = 0; j < i; ++j)
80 SCIP_CALL( SCIPgetTransformedVar(scip, _arc_var[i][j], &_arc_var[i][j]) );
81 SCIP_CALL( SCIPgetTransformedCons(scip, _arc_con[i][j], &_arc_con[i][j]) );
84 for (
int i = 1; i < num_nodes(); ++i)
86 SCIP_CALL( SCIPgetTransformedCons(scip, _part_con[i], &_part_con[i]) );
103 vector< vector<SCIP_Real> > red_length(
num_nodes());
105 red_length[i].resize(i, 0.0);
112 assert( i == 0 ||
part_con(i) != 0 );
113 for (
int j = 0; j < i; ++j)
118 r -= SCIPgetDualfarkasLinear(scip,
arc_con(i,j));
120 r -= 0.5 * SCIPgetDualfarkasLinear(scip,
part_con(j));
122 r -= 0.5 * SCIPgetDualfarkasLinear(scip,
part_con(i));
123 red_length[i][j] = r;
131 assert( i == 0 ||
part_con(i) != 0 );
132 for (
int j = 0; j < i; ++j)
137 r -= SCIPgetDualsolLinear(scip,
arc_con(i,j));
139 r -= 0.5 * SCIPgetDualsolLinear(scip,
part_con(j));
141 r -= 0.5 * SCIPgetDualsolLinear(scip,
part_con(i));
142 red_length[i][j] = r;
150 SCIPinfoMessage(scip, NULL,
"dual ray solution:\n");
153 for (
int j = 0; j < i; ++j)
154 SCIPinfoMessage(scip, NULL,
"arc_%d_%d: %g\n", i, j, SCIPgetDualfarkasLinear(scip,
arc_con(i,j)));
158 SCIPinfoMessage(scip, NULL,
"part_%d: %g\n", i, SCIPgetDualfarkasLinear(scip,
part_con(i)));
162 for (
int j = 0; j < i; ++j)
163 SCIPinfoMessage(scip, NULL,
"length_%d_%d: %g\n", i, j, red_length[i][j]);
168 SCIPinfoMessage(scip, NULL,
"dual solution:\n");
171 for (
int j = 0; j < i; ++j)
172 SCIPinfoMessage(scip, NULL,
"arc_%d_%d: %g\n", i, j, SCIPgetDualsolLinear(scip,
arc_con(i,j)));
176 SCIPinfoMessage(scip, NULL,
"part_%d: %g\n", i, SCIPgetDualsolLinear(scip,
part_con(i)));
180 for (
int j = 0; j < i; ++j)
181 SCIPinfoMessage(scip, NULL,
"length_%d_%d: %g\n", i, j, red_length[i][j]);
191 if ( SCIPisNegative(scip, reduced_cost) )
197 SCIP_CALL( SCIPwriteTransProblem(scip,
"vrp.lp",
"lp", FALSE) );
218 SCIPdebugMessage(
"call scip_redcost ...\n");
221 *result = SCIP_SUCCESS;
224 SCIP_CALL( pricing(scip,
false) );
239 SCIPdebugMessage(
"call scip_farkas ...\n");
242 SCIP_CALL( pricing(scip,
true) );
251 const list<int>& tour
257 SCIPsnprintf(var_name, 255,
"T");
258 for (list<int>::const_iterator it = tour.begin(); it != tour.end(); ++it)
260 strncpy(tmp_name, var_name, 255);
261 SCIPsnprintf(var_name, 255,
"%s_%d", tmp_name, *it);
263 SCIPdebugMessage(
"new variable <%s>\n", var_name);
270 SCIP_CALL( SCIPcreateVar(scip, &var, var_name,
274 SCIP_VARTYPE_CONTINUOUS,
275 false,
false, 0, 0, 0, 0, 0) );
278 SCIP_CALL( SCIPaddPricedVar(scip, var, 1.0) );
281 for (list<int>::const_iterator it = tour.begin(); it != tour.end(); ++it)
284 SCIP_CALL( SCIPaddCoefLinear(scip,
part_con(*it), var, 1.0) );
289 for (list<int>::const_iterator it = tour.begin(); it != tour.end(); ++it)
292 SCIP_CALL( SCIPaddCoefLinear(scip,
arc_con(last, *it), var, 1.0) );
295 SCIP_CALL( SCIPaddCoefLinear(scip,
arc_con(last, 0), var, 1.0 ) );
298 SCIP_CALL( SCIPreleaseVar(scip, &var) );
313 static const SCIP_Real eps = 1e-9;
320 PQUEUE_KEY() : demand(0), length(0.0) {}
323 bool operator< (
const PQUEUE_KEY& l1,
const PQUEUE_KEY& l2)
325 if ( l1.demand < l2.demand )
327 if ( l1.demand > l2.demand )
329 if ( l1.length < l2.length-eps )
338 typedef int PQUEUE_DATA;
340 typedef PQUEUE::pqueue_item PQUEUE_ITEM;
344 struct NODE_TABLE_DATA
348 PQUEUE::pqueue_item queue_item;
350 NODE_TABLE_DATA( ) : length(0.0), predecessor(-1), queue_item( NULL ) {}
353 typedef int NODE_TABLE_KEY;
354 typedef std::map< NODE_TABLE_KEY, NODE_TABLE_DATA > NODE_TABLE;
365 const vector< vector<SCIP_Real> >& length,
371 SCIPdebugMessage(
"Enter RSP - capacity: %d\n",
capacity());
378 PQUEUE_KEY queue_key;
379 PQUEUE_DATA queue_data = 0;
380 PQUEUE_ITEM queue_item = PQ.insert(queue_key, queue_data);
382 NODE_TABLE_KEY table_key = 0;
383 NODE_TABLE_DATA table_entry;
386 while ( ! PQ.empty() )
389 queue_item = PQ.top();
390 queue_key = PQ.get_key (queue_item);
391 queue_data = PQ.get_data(queue_item);
395 const int curr_node = queue_data;
396 const SCIP_Real curr_length = queue_key.length;
397 const int curr_demand = queue_key.demand;
400 if ( curr_node == 0 && curr_length < -eps )
404 if ( curr_node == 0 && curr_demand != 0 )
408 for (
int next_node = 0; next_node <
num_nodes(); ++next_node)
410 if ( next_node == curr_node )
412 if (
have_edge( next_node, curr_node ) == false )
415 const int next_demand = curr_demand +
demand(next_node);
420 const SCIP_Real next_length = curr_length + ( curr_node > next_node ?
421 length[curr_node][next_node] :
422 length[next_node][curr_node] );
424 NODE_TABLE& next_table = table[next_node];
428 list<NODE_TABLE::iterator> dominated;
430 for (NODE_TABLE::iterator it = next_table.begin(); it != next_table.end() && ! skip; ++it)
432 if ( next_demand >= it->first && next_length >= it->second.length - eps )
435 if ( next_demand <= it->first && next_length <= it->second.length + eps )
436 dominated.push_front( it );
442 for (list<NODE_TABLE::iterator>::iterator it = dominated.begin(); it != dominated.end(); ++it)
444 PQ.remove( (*it)->second.queue_item );
445 next_table.erase( *it );
449 queue_key.demand = next_demand;
450 queue_key.length = next_length;
451 queue_data = next_node;
453 queue_item = PQ.insert(queue_key, queue_data);
455 table_key = next_demand;
456 table_entry.length = next_length;
457 table_entry.predecessor = curr_node;
458 table_entry.queue_item = queue_item;
460 next_table[table_key] = table_entry;
463 printf(
"new entry node = %d demand = %d length = %g pref = %d\n", next_node, next_demand, next_length, curr_node);
468 SCIPdebugMessage(
"Done RSP DP.\n");
470 table_entry.predecessor = -1;
471 table_entry.length = 0;
475 for (NODE_TABLE::iterator it = table[0].begin(); it != table[0].end(); ++it)
477 if ( it->second.length < table_entry.length )
479 table_key = it->first;
480 table_entry = it->second;
483 SCIP_Real tour_length = table_entry.length;
485 while ( table_entry.predecessor > 0 )
487 table_key -=
demand(curr_node);
488 curr_node = table_entry.predecessor;
489 tour.push_front(curr_node);
490 table_entry = table[curr_node][table_key];
493 SCIPdebugMessage(
"Leave RSP tour length = %g\n", tour_length);