Google OR-Tools: ortools/graph/min_cost

Google OR-Tools: ortools/graph/min_cost_flow.cc Source File

16#include <algorithm>

17#include <cmath>

18#include <cstdint>

19#include <cstdlib>

20#include <limits>

21#include <memory>

22#include <string>

23#include <vector>

25#include "absl/base/attributes.h"

26#include "absl/flags/flag.h"

27#include "absl/log/check.h"

28#include "absl/log/log.h"

29#include "absl/strings/str_format.h"

30#include "absl/strings/string_view.h"

40 "Divide factor for epsilon at each refine step.");

41ABSL_FLAG(bool, min_cost_flow_check_feasibility, true,

42 "Check that the graph has enough capacity to send all supplies "

43 "and serve all demands. Also check that the sum of supplies "

44 "is equal to the sum of demands.");

45ABSL_FLAG(bool, min_cost_flow_check_result, true,

46 "Check that the result is valid.");

50template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>

54 alpha_(absl::GetFlag(FLAGS_min_cost_flow_alpha)),

55 stats_("MinCostFlow"),

56 check_feasibility_(absl::GetFlag(FLAGS_min_cost_flow_check_feasibility)) {

60 const NodeIndex max_num_nodes = graph_->node_capacity();

61 if (max_num_nodes > 0) {

62 first_admissible_arc_ = std::make_unique<ArcIndex[]>(max_num_nodes);

63 std::fill(first_admissible_arc_.get(),

64 first_admissible_arc_.get() + max_num_nodes, Graph::kNilArc);

65 node_potential_ = std::make_unique<CostValue[]>(max_num_nodes);

66 node_excess_ = std::make_unique<FlowQuantity[]>(max_num_nodes);

67 initial_node_excess_ = std::make_unique<FlowQuantity[]>(max_num_nodes);

68 }

69 const ArcIndex max_num_arcs = graph_->arc_capacity();

70 if (max_num_arcs > 0) {

71 residual_arc_capacity_.Reserve(-max_num_arcs, max_num_arcs - 1);

72 residual_arc_capacity_.SetAll(0);

73 scaled_arc_unit_cost_.Reserve(-max_num_arcs, max_num_arcs - 1);

74 scaled_arc_unit_cost_.SetAll(0);

75 }

76}

78template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>

81 DCHECK(graph_->IsNodeValid(node));

82 node_excess_[node] = supply;

83 initial_node_excess_[node] = supply;

85 feasibility_checked_ = false;

86}

88template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>

90 ArcIndex arc, ArcScaledCostType unit_cost) {

91 DCHECK(IsArcDirect(arc));

92 scaled_arc_unit_cost_.Set(arc, unit_cost);

93 scaled_arc_unit_cost_.Set(Opposite(arc), -scaled_arc_unit_cost_[arc]);

95 feasibility_checked_ = false;

96}

98template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>

100 ArcIndex arc, ArcFlowType new_capacity) {

101 DCHECK_LE(0, new_capacity);

102 DCHECK(IsArcDirect(arc));

103 const FlowQuantity free_capacity = residual_arc_capacity_[arc];

105 if (capacity_delta == 0) {

106 return;

107 }

109 feasibility_checked_ = false;

110 const FlowQuantity new_availability = free_capacity + capacity_delta;

111 if (new_availability >= 0) {

112

113

114

115

116

117 DCHECK((capacity_delta > 0) ||

118 (capacity_delta < 0 && new_availability >= 0));

119 residual_arc_capacity_.Set(arc, new_availability);

120 DCHECK_LE(0, residual_arc_capacity_[arc]);

121 } else {

122

123

124 const FlowQuantity flow = residual_arc_capacity_[Opposite(arc)];

125 const FlowQuantity flow_excess = flow - new_capacity;

126 residual_arc_capacity_.Set(arc, 0);

127 residual_arc_capacity_.Set(Opposite(arc), new_capacity);

128 node_excess_[Tail(arc)] += flow_excess;

129 node_excess_[Head(arc)] -= flow_excess;

130 DCHECK_LE(0, residual_arc_capacity_[arc]);

131 DCHECK_LE(0, residual_arc_capacity_[Opposite(arc)]);

132 }

133}

134

135template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>

137 ArcScaledCostType>::CheckInputConsistency() {

138 const FlowQuantity kMaxFlow = std::numeric_limits<FlowQuantity>::max();

139

140

141

142 FlowQuantity sum_of_supplies = 0;

143 FlowQuantity sum_of_demands = 0;

144 for (NodeIndex node = 0; node < graph_->num_nodes(); ++node) {

145 const FlowQuantity excess = node_excess_[node];

146 if (excess > 0) {

147 sum_of_supplies = CapAdd(sum_of_supplies, excess);

148 if (sum_of_supplies >= kMaxFlow) {

149 status_ = BAD_CAPACITY_RANGE;

150 LOG(ERROR) << "Input consistency error: sum of supplies overflow";

151 return false;

152 }

153 } else if (excess < 0) {

154 sum_of_demands = CapAdd(sum_of_demands, -excess);

155 if (sum_of_demands >= kMaxFlow) {

156 status_ = BAD_CAPACITY_RANGE;

157 LOG(ERROR) << "Input consistency error: sum of demands overflow";

158 return false;

159 }

160 }

161 }

162 if (sum_of_supplies != sum_of_demands) {

163 status_ = UNBALANCED;

164 LOG(ERROR) << "Input consistency error: unbalanced problem";

165 return false;

166 }

167

168 std::vector<FlowQuantity> max_node_excess(

169 node_excess_.get(), node_excess_.get() + graph_->num_nodes());

170 std::vector<FlowQuantity> min_node_excess = max_node_excess;

171 for (ArcIndex arc = 0; arc < graph_->num_arcs(); ++arc) {

172 const FlowQuantity capacity = residual_arc_capacity_[arc];

173 CHECK_GE(capacity, 0);

174 if (capacity == 0) continue;

175 const int tail = graph_->Tail(arc);

176 const int head = graph_->Head(arc);

177 min_node_excess[tail] = CapSub(min_node_excess[tail], capacity);

178 max_node_excess[head] = CapAdd(max_node_excess[head], capacity);

179 }

180

181 const int num_nodes = graph_->num_nodes();

182 for (NodeIndex node = 0; node < num_nodes; ++node) {

183 if (max_node_excess[node] >= kMaxFlow ||

184 min_node_excess[node] <= -kMaxFlow) {

185

186

187

188 if (max_node_excess[node] < std::numeric_limits<ArcFlowType>::max()) {

189

190 const ArcFlowType upper_bound =

191 std::max<ArcFlowType>(0, max_node_excess[node]);

192

193

194 min_node_excess[node] = node_excess_[node];

195 for (const ArcIndex arc : graph_->OutgoingArcs(node)) {

196 residual_arc_capacity_[arc] =

197 std::min(residual_arc_capacity_[arc], upper_bound);

198 min_node_excess[node] =

199 CapSub(min_node_excess[node], residual_arc_capacity_[arc]);

200 }

201 if (min_node_excess[node] > -kMaxFlow) continue;

202 }

203 if (min_node_excess[node] > -std::numeric_limits<ArcFlowType>::max()) {

204

205 const ArcFlowType upper_bound =

206 std::max<ArcFlowType>(0, -min_node_excess[node]);

207

208

209 max_node_excess[node] = node_excess_[node];

210 for (const ArcIndex arc : graph_->IncomingArcs(node)) {

211 residual_arc_capacity_[arc] =

212 std::min(residual_arc_capacity_[arc], upper_bound);

213 max_node_excess[node] =

214 CapAdd(max_node_excess[node], residual_arc_capacity_[arc]);

215 }

216 if (max_node_excess[node] < kMaxFlow) continue;

217 }

218

219 status_ = BAD_CAPACITY_RANGE;

220 LOG(ERROR) << "Maximum in or out flow of node + excess " << node

221 << " overflow the FlowQuantity type (int64_t).";

222 return false;

223 }

224 }

225 return true;

226}

227

228template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>

229bool GenericMinCostFlow<Graph, ArcFlowType, ArcScaledCostType>::CheckResult()

230 const {

231 for (NodeIndex node = 0; node < graph_->num_nodes(); ++node) {

232 if (node_excess_[node] != 0) {

233 LOG(DFATAL) << "node_excess_[" << node << "] != 0";

234 return false;

235 }

236 for (OutgoingOrOppositeIncomingArcIterator it(*graph_, node); it.Ok();

237 it.Next()) {

238 const ArcIndex arc = it.Index();

239 bool ok = true;

240 if (residual_arc_capacity_[arc] < 0) {

241 LOG(DFATAL) << "residual_arc_capacity_[" << arc << "] < 0";

242 ok = false;

243 }

244 if (residual_arc_capacity_[arc] > 0 && ReducedCost(arc) < -epsilon_) {

245 LOG(DFATAL) << "residual_arc_capacity_[" << arc

246 << "] > 0 && ReducedCost(" << arc << ") < " << -epsilon_

247 << ". (epsilon_ = " << epsilon_ << ").";

248 ok = false;

249 }

250 if (!ok) {

251 LOG(DFATAL) << DebugString("CheckResult ", arc);

252 return false;

253 }

254 }

255 }

256 return true;

257}

258

259template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>

262

263

264

265

266

267

268 DCHECK_GE(node_excess_[node], 0);

269 for (OutgoingOrOppositeIncomingArcIterator it(*graph_, node); it.Ok();

270 it.Next()) {

271 const ArcIndex arc = it.Index();

272 DCHECK(!IsAdmissible(arc)) << DebugString("CheckRelabelPrecondition:", arc);

273 }

274 return true;

275}

276

277template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>

278std::string

279GenericMinCostFlow<Graph, ArcFlowType, ArcScaledCostType>::DebugString(

280 absl::string_view context, ArcIndex arc) const {

281 const NodeIndex tail = Tail(arc);

282 const NodeIndex head = Head(arc);

283

284

285

286 const CostValue reduced_cost = scaled_arc_unit_cost_[arc] +

287 node_potential_[tail] - node_potential_[head];

288 return absl::StrFormat(

289 "%s Arc %d, from %d to %d, "

290 "Capacity = %d, Residual capacity = %d, "

291 "Flow = residual capacity for reverse arc = %d, "

292 "Height(tail) = %d, Height(head) = %d, "

293 "Excess(tail) = %d, Excess(head) = %d, "

294 "Cost = %d, Reduced cost = %d, ",

295 context, arc, tail, head, Capacity(arc),

296 static_cast<FlowQuantity>(residual_arc_capacity_[arc]), Flow(arc),

297 node_potential_[tail], node_potential_[head], node_excess_[tail],

298 node_excess_[head], static_cast<CostValue>(scaled_arc_unit_cost_[arc]),

299 reduced_cost);

300}

301

302template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>

304 ArcScaledCostType>::CheckFeasibility() {

306

307

308

309

310

311

312

313

314

315

316 feasibility_checked_ = false;

317 ArcIndex num_extra_arcs = 0;

318 for (NodeIndex node = 0; node < graph_->num_nodes(); ++node) {

319 if (initial_node_excess_[node] != 0) {

320 ++num_extra_arcs;

321 }

322 }

323 const NodeIndex num_nodes_in_max_flow = graph_->num_nodes() + 2;

324 const ArcIndex num_arcs_in_max_flow = graph_->num_arcs() + num_extra_arcs;

325 const NodeIndex source = num_nodes_in_max_flow - 2;

326 const NodeIndex sink = num_nodes_in_max_flow - 1;

327 using CheckerGraph = ::util::ReverseArcListGraph<>;

328 CheckerGraph checker_graph(num_nodes_in_max_flow, num_arcs_in_max_flow);

330

331 for (ArcIndex arc = 0; arc < graph_->num_arcs(); ++arc) {

332 const ArcIndex new_arc =

333 checker_graph.AddArc(graph_->Tail(arc), graph_->Head(arc));

334 DCHECK_EQ(arc, new_arc);

335 checker.SetArcCapacity(new_arc, Capacity(arc));

336 }

337 FlowQuantity total_demand = 0;

338 FlowQuantity total_supply = 0;

339

340 for (NodeIndex node = 0; node < graph_->num_nodes(); ++node) {

341 const FlowQuantity supply = initial_node_excess_[node];

342 if (supply > 0) {

343 const ArcIndex new_arc = checker_graph.AddArc(source, node);

344 checker.SetArcCapacity(new_arc, supply);

345 total_supply += supply;

346 } else if (supply < 0) {

347 const ArcIndex new_arc = checker_graph.AddArc(node, sink);

348 checker.SetArcCapacity(new_arc, -supply);

349 total_demand -= supply;

350 }

351 }

352 if (total_supply != total_demand) {

353 LOG(DFATAL) << "total_supply(" << total_supply << ") != total_demand("

354 << total_demand << ").";

355 return false;

356 }

357 if (!checker.Solve()) {

358 LOG(DFATAL) << "Max flow could not be computed.";

359 return false;

360 }

361 const FlowQuantity optimal_max_flow = checker.GetOptimalFlow();

362 feasibility_checked_ = true;

363 return optimal_max_flow == total_supply;

364}

365

366template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>

369 if (IsArcDirect(arc)) {

370 return residual_arc_capacity_[Opposite(arc)];

371 } else {

372 return -residual_arc_capacity_[arc];

373 }

374}

375

376

377template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>

380 if (IsArcDirect(arc)) {

381 return residual_arc_capacity_[arc] + residual_arc_capacity_[Opposite(arc)];

382 } else {

383 return 0;

384 }

385}

386

387template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>

390 DCHECK(IsArcValid(arc));

391 DCHECK_EQ(uint64_t{1}, cost_scaling_factor_);

392 return scaled_arc_unit_cost_[arc];

393}

394

395template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>

398 DCHECK(graph_->IsNodeValid(node));

399 return node_excess_[node];

400}

401

402template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>

403bool GenericMinCostFlow<Graph, ArcFlowType, ArcScaledCostType>::IsAdmissible(

404 ArcIndex arc) const {

405 return FastIsAdmissible(arc, node_potential_[Tail(arc)]);

406}

407

408template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>

409bool GenericMinCostFlow<Graph, ArcFlowType, ArcScaledCostType>::

410 FastIsAdmissible(ArcIndex arc, CostValue tail_potential) const {

411 DCHECK_EQ(node_potential_[Tail(arc)], tail_potential);

412 return residual_arc_capacity_[arc] > 0 &&

413 FastReducedCost(arc, tail_potential) < 0;

414}

415

416template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>

417bool GenericMinCostFlow<Graph, ArcFlowType, ArcScaledCostType>::IsActive(

419 return node_excess_[node] > 0;

420}

421

422template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>

423auto GenericMinCostFlow<Graph, ArcFlowType, ArcScaledCostType>::ReducedCost(

424 ArcIndex arc) const -> CostValue {

425 return FastReducedCost(arc, node_potential_[Tail(arc)]);

426}

427

428template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>

429auto GenericMinCostFlow<Graph, ArcFlowType, ArcScaledCostType>::FastReducedCost(

431 DCHECK_EQ(node_potential_[Tail(arc)], tail_potential);

432 DCHECK(graph_->IsNodeValid(Tail(arc)));

433 DCHECK(graph_->IsNodeValid(Head(arc)));

434 DCHECK_LE(node_potential_[Tail(arc)], 0) << DebugString("ReducedCost:", arc);

435 DCHECK_LE(node_potential_[Head(arc)], 0) << DebugString("ReducedCost:", arc);

436

437

438

439 return scaled_arc_unit_cost_[arc] + tail_potential -

440 node_potential_[Head(arc)];

441}

442

443template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>

446 OutgoingOrOppositeIncomingArcIterator arc_it(*graph_, node);

447 return arc_it.Index();

448}

449

450template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>

452 if (!CheckInputConsistency()) {

453 return false;

454 }

455 if (check_feasibility_ && !CheckFeasibility()) {

457 return false;

458 }

459

461 std::fill(node_potential_.get(),

462 node_potential_.get() + graph_->node_capacity(), 0);

463

464 ResetFirstAdmissibleArcs();

465 if (!ScaleCosts()) return false;

466 if (!Optimize()) return false;

469

470 if (absl::GetFlag(FLAGS_min_cost_flow_check_result) && !CheckResult()) {

472 UnscaleCosts();

473 return false;

474 }

475 UnscaleCosts();

476

478 return true;

479}

480

481template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>

484 if (status_ != OPTIMAL) {

485 return 0;

486 }

487

488

489

491 const CostValue kMaxCost = std::numeric_limits<CostValue>::max();

492 const CostValue kMinCost = std::numeric_limits<CostValue>::min();

493 for (ArcIndex arc = 0; arc < graph_->num_arcs(); ++arc) {

494 const CostValue flow_on_arc = residual_arc_capacity_[Opposite(arc)];

496 CapProd(scaled_arc_unit_cost_[arc], flow_on_arc);

497 if (flow_cost == kMaxCost || flow_cost == kMinCost) return kMaxCost;

498 total_flow_cost = CapAdd(flow_cost, total_flow_cost);

499 if (total_flow_cost == kMaxCost || total_flow_cost == kMinCost) {

500 return kMaxCost;

501 }

502 }

503 return total_flow_cost;

504}

505

506template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>

508 ArcScaledCostType>::ResetFirstAdmissibleArcs() {

509 for (NodeIndex node = 0; node < graph_->num_nodes(); ++node) {

510 first_admissible_arc_[node] = GetFirstOutgoingOrOppositeIncomingArc(node);

511 }

512}

513

514template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>

515bool GenericMinCostFlow<Graph, ArcFlowType, ArcScaledCostType>::ScaleCosts() {

517 cost_scaling_factor_ = scale_prices_ ? graph_->num_nodes() + 1 : 1;

518 epsilon_ = 1LL;

519 VLOG(3) << "Number of nodes in the graph = " << graph_->num_nodes();

520 VLOG(3) << "Number of arcs in the graph = " << graph_->num_arcs();

522 std::numeric_limits<CostValue>::max() / (2 * cost_scaling_factor_);

523 for (ArcIndex arc = 0; arc < graph_->num_arcs(); ++arc) {

524 if (scaled_arc_unit_cost_[arc] > threshold ||

525 scaled_arc_unit_cost_[arc] < -threshold) {

526 status_ = BAD_COST_RANGE;

527 return false;

528 }

529 const CostValue cost = scaled_arc_unit_cost_[arc] * cost_scaling_factor_;

530 scaled_arc_unit_cost_.Set(arc, cost);

531 scaled_arc_unit_cost_.Set(Opposite(arc), -cost);

532 epsilon_ = std::max(epsilon_, MathUtil::Abs(cost));

533 }

534

535

536

537

538 overflow_threshold_ = std::numeric_limits<CostValue>::min() + epsilon_;

539

540 VLOG(3) << "Initial epsilon = " << epsilon_;

541 VLOG(3) << "Cost scaling factor = " << cost_scaling_factor_;

542 return true;

543}

544

545template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>

546void GenericMinCostFlow<Graph, ArcFlowType, ArcScaledCostType>::UnscaleCosts() {

548 for (ArcIndex arc = 0; arc < graph_->num_arcs(); ++arc) {

549 const CostValue cost = scaled_arc_unit_cost_[arc] / cost_scaling_factor_;

550 scaled_arc_unit_cost_.Set(arc, cost);

551 scaled_arc_unit_cost_.Set(Opposite(arc), -cost);

552 }

553 cost_scaling_factor_ = 1;

554}

555

556template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>

557bool GenericMinCostFlow<Graph, ArcFlowType, ArcScaledCostType>::Optimize() {

558 const CostValue kEpsilonMin = 1LL;

559 num_relabels_since_last_price_update_ = 0;

560 do {

561

562 epsilon_ = std::max(epsilon_ / alpha_, kEpsilonMin);

563 VLOG(3) << "Epsilon changed to: " << epsilon_;

564 if (!Refine()) return false;

565 } while (epsilon_ != 1LL);

566 return true;

567}

568

569template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>

571 ArcScaledCostType>::SaturateAdmissibleArcs() {

573 for (NodeIndex node = 0; node < graph_->num_nodes(); ++node) {

574 const CostValue tail_potential = node_potential_[node];

575 for (OutgoingOrOppositeIncomingArcIterator it(*graph_, node,

576 first_admissible_arc_[node]);

577 it.Ok(); it.Next()) {

578 const ArcIndex arc = it.Index();

579 if (FastIsAdmissible(arc, tail_potential)) {

580 FastPushFlow(residual_arc_capacity_[arc], arc, node);

581 }

582 }

583

584

585

586

587

588

589

591 }

592}

593

594template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>

595void GenericMinCostFlow<Graph, ArcFlowType, ArcScaledCostType>::PushFlow(

596 FlowQuantity flow, ArcIndex arc) {

598 FastPushFlow(flow, arc, Tail(arc));

599}

600

601template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>

602void GenericMinCostFlow<Graph, ArcFlowType, ArcScaledCostType>::FastPushFlow(

603 FlowQuantity flow, ArcIndex arc, NodeIndex tail) {

605 DCHECK_EQ(Tail(arc), tail);

606 DCHECK_GT(residual_arc_capacity_[arc], 0);

607 DCHECK_LE(flow, residual_arc_capacity_[arc]);

608

609 residual_arc_capacity_.Set(arc, residual_arc_capacity_[arc] - flow);

610

611 const ArcIndex opposite = Opposite(arc);

612 residual_arc_capacity_.Set(opposite, residual_arc_capacity_[opposite] + flow);

613

614 node_excess_[tail] -= flow;

615 node_excess_[Head(arc)] += flow;

616}

617

618template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>

620 ArcScaledCostType>::InitializeActiveNodeStack() {

622 DCHECK(active_nodes_.empty());

623 for (NodeIndex node = 0; node < graph_->num_nodes(); ++node) {

624 if (IsActive(node)) {

625 active_nodes_.push(node);

626 }

627 }

628}

629

630template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>

631bool GenericMinCostFlow<Graph, ArcFlowType, ArcScaledCostType>::UpdatePrices() {

633

634

635

636

637

638

639

640

641

642

643

644

645

646

647

648

649

650 const NodeIndex num_nodes = graph_->num_nodes();

651 std::vector<NodeIndex> bfs_queue;

652 std::vector<bool> node_in_queue(num_nodes, false);

653

654

655 const CostValue kMinCostValue = std::numeric_limits<CostValue>::min();

656 std::vector<CostValue> min_non_admissible_potential(num_nodes, kMinCostValue);

657 std::vector<NodeIndex> nodes_to_process;

658

659

660 FlowQuantity remaining_excess = 0;

661

662

663 for (NodeIndex node = 0; node < num_nodes; ++node) {

664 if (node_excess_[node] < 0) {

665 bfs_queue.push_back(node);

666 node_in_queue[node] = true;

667

668

669 remaining_excess -= node_excess_[node];

670 }

671 }

672

673

674

675

676

678

679 int queue_index = 0;

680 while (remaining_excess > 0) {

681

682

683

684

685 for (; queue_index < bfs_queue.size(); ++queue_index) {

686 DCHECK_GE(num_nodes, bfs_queue.size());

687 const NodeIndex node = bfs_queue[queue_index];

688 for (OutgoingOrOppositeIncomingArcIterator it(*graph_, node); it.Ok();

689 it.Next()) {

690 const NodeIndex head = Head(it.Index());

691 if (node_in_queue[head]) continue;

692 const ArcIndex opposite_arc = Opposite(it.Index());

693 if (residual_arc_capacity_[opposite_arc] > 0) {

694 node_potential_[head] += potential_delta;

695 if (node_potential_[head] < overflow_threshold_) {

696 status_ = BAD_COST_RANGE;

697 return false;

698 }

699 if (ReducedCost(opposite_arc) < 0) {

700 DCHECK(IsAdmissible(opposite_arc));

701

702

703

704

705 remaining_excess -= node_excess_[head];

706 if (remaining_excess == 0) {

707 node_potential_[head] -= potential_delta;

708 break;

709 }

710 bfs_queue.push_back(head);

711 node_in_queue[head] = true;

712 if (potential_delta < 0) {

713 first_admissible_arc_[head] =

714 GetFirstOutgoingOrOppositeIncomingArc(head);

715 }

716 } else {

717

718

719 node_potential_[head] -= potential_delta;

720 if (min_non_admissible_potential[head] == kMinCostValue) {

721 nodes_to_process.push_back(head);

722 }

723 min_non_admissible_potential[head] = std::max(

724 min_non_admissible_potential[head],

725 node_potential_[node] - scaled_arc_unit_cost_[opposite_arc]);

726 }

727 }

728 }

729 if (remaining_excess == 0) break;

730 }

731 if (remaining_excess == 0) break;

732

733

734

735 CostValue max_potential_diff = kMinCostValue;

736 for (int i = 0; i < nodes_to_process.size(); ++i) {

737 const NodeIndex node = nodes_to_process[i];

738 if (node_in_queue[node]) continue;

739 max_potential_diff =

740 std::max(max_potential_diff,

741 min_non_admissible_potential[node] - node_potential_[node]);

742 if (max_potential_diff == potential_delta) break;

743 }

744 DCHECK_LE(max_potential_diff, potential_delta);

745 potential_delta = max_potential_diff - epsilon_;

746

747

748

749

750

751

752

753 int index = 0;

754 for (int i = 0; i < nodes_to_process.size(); ++i) {

755 const NodeIndex node = nodes_to_process[i];

756 if (node_in_queue[node]) continue;

757 if (node_potential_[node] + potential_delta <

758 min_non_admissible_potential[node]) {

759 node_potential_[node] += potential_delta;

760 if (node_potential_[node] < overflow_threshold_) {

761 status_ = BAD_COST_RANGE;

762 return false;

763 }

764 first_admissible_arc_[node] =

765 GetFirstOutgoingOrOppositeIncomingArc(node);

766 bfs_queue.push_back(node);

767 node_in_queue[node] = true;

768 remaining_excess -= node_excess_[node];

769 continue;

770 }

771

772

773 nodes_to_process[index] = node;

774 ++index;

775 }

776 nodes_to_process.resize(index);

777 }

778

779

780 if (potential_delta == 0) return true;

781 for (NodeIndex node = 0; node < num_nodes; ++node) {

782 if (!node_in_queue[node]) {

783 node_potential_[node] += potential_delta;

784 if (node_potential_[node] < overflow_threshold_) {

785 status_ = BAD_COST_RANGE;

786 return false;

787 }

788 first_admissible_arc_[node] = GetFirstOutgoingOrOppositeIncomingArc(node);

789 }

790 }

791 return true;

792}

793

794template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>

795bool GenericMinCostFlow<Graph, ArcFlowType, ArcScaledCostType>::Refine() {

797 SaturateAdmissibleArcs();

798 InitializeActiveNodeStack();

799

800 const NodeIndex num_nodes = graph_->num_nodes();

801 while (!active_nodes_.empty()) {

802

803 if (num_relabels_since_last_price_update_ >= num_nodes) {

804 num_relabels_since_last_price_update_ = 0;

805 if (use_price_update_) {

806 if (!UpdatePrices()) return false;

807 }

808 }

809 const NodeIndex node = active_nodes_.top();

810 active_nodes_.pop();

811 DCHECK(IsActive(node));

812 if (!Discharge(node)) return false;

813 }

814 return true;

815}

816

817template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>

818bool GenericMinCostFlow<Graph, ArcFlowType, ArcScaledCostType>::Discharge(

821 while (true) {

822

823

824 DCHECK(IsActive(node));

825 const CostValue tail_potential = node_potential_[node];

826 for (OutgoingOrOppositeIncomingArcIterator it(*graph_, node,

827 first_admissible_arc_[node]);

828 it.Ok(); it.Next()) {

829 const ArcIndex arc = it.Index();

830 if (!FastIsAdmissible(arc, tail_potential)) continue;

831

832

833

834

835 const NodeIndex head = Head(arc);

836 if (node_excess_[head] >= 0 && !NodeHasAdmissibleArc(head)) {

837

838 if (!Relabel(head)) return false;

839 if (!FastIsAdmissible(arc, tail_potential)) continue;

840 }

841

842 const FlowQuantity delta =

843 std::min(node_excess_[node],

844 static_cast<FlowQuantity>(residual_arc_capacity_[arc]));

845 const bool head_active_before_push = IsActive(head);

846 FastPushFlow(delta, arc, node);

847 if (IsActive(head) && !head_active_before_push) {

848 active_nodes_.push(head);

849 }

850

851 if (node_excess_[node] == 0) {

852

853 first_admissible_arc_[node] = arc;

854 return true;

855 }

856 }

857 if (!Relabel(node)) return false;

858 }

859}

860

861template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>

865 const CostValue tail_potential = node_potential_[node];

866 for (OutgoingOrOppositeIncomingArcIterator it(*graph_, node,

867 first_admissible_arc_[node]);

868 it.Ok(); it.Next()) {

869 const ArcIndex arc = it.Index();

870 if (FastIsAdmissible(arc, tail_potential)) {

871 first_admissible_arc_[node] = arc;

872 return true;

873 }

874 }

875 return false;

876}

877

878

879template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>

880bool GenericMinCostFlow<Graph, ArcFlowType, ArcScaledCostType>::Relabel(

883 DCHECK(CheckRelabelPrecondition(node));

884 ++num_relabels_since_last_price_update_;

885

886

887

888

889

890

891 const CostValue guaranteed_new_potential = node_potential_[node] - epsilon_;

892 if (guaranteed_new_potential < overflow_threshold_) {

893 status_ = BAD_COST_RANGE;

894 return false;

895 }

896

897

898

899

900

901

902 const CostValue kMinCostValue = std::numeric_limits<CostValue>::min();

903 CostValue min_non_admissible_potential = kMinCostValue;

904

905

906

907

908 CostValue previous_min_non_admissible_potential = kMinCostValue;

910

911 for (OutgoingOrOppositeIncomingArcIterator it(*graph_, node); it.Ok();

912 it.Next()) {

913 const ArcIndex arc = it.Index();

914 if (residual_arc_capacity_[arc] > 0) {

915 const CostValue min_non_admissible_potential_for_arc =

916 node_potential_[Head(arc)] - scaled_arc_unit_cost_[arc];

917

918

919

920 if (min_non_admissible_potential_for_arc > guaranteed_new_potential) {

921

922

923

924 node_potential_[node] = guaranteed_new_potential;

925 first_admissible_arc_[node] = arc;

926 return true;

927 }

928 if (min_non_admissible_potential_for_arc > min_non_admissible_potential) {

929 previous_min_non_admissible_potential = min_non_admissible_potential;

930 min_non_admissible_potential = min_non_admissible_potential_for_arc;

931 first_arc = arc;

932 }

933 }

934 }

935

936

937

938

939

940 if (min_non_admissible_potential == kMinCostValue) {

941 if (node_excess_[node] != 0) {

942

943

945 return false;

946 } else {

947

948

949

950

951 node_potential_[node] = guaranteed_new_potential;

953 }

954 return true;

955 }

956

957

958

959

960 if (min_non_admissible_potential < overflow_threshold_) {

961 status_ = BAD_COST_RANGE;

962 return false;

963 }

964 const CostValue new_potential = min_non_admissible_potential - epsilon_;

965 if (new_potential < overflow_threshold_) {

966 status_ = BAD_COST_RANGE;

967 return false;

968 }

969 node_potential_[node] = new_potential;

970 if (previous_min_non_admissible_potential <= new_potential) {

971 first_admissible_arc_[node] = first_arc;

972 } else {

973

974 first_admissible_arc_[node] = GetFirstOutgoingOrOppositeIncomingArc(node);

975 }

976 return true;

977}

978

979template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>

981GenericMinCostFlow<Graph, ArcFlowType, ArcScaledCostType>::Opposite(

982 ArcIndex arc) const {

983 return graph_->OppositeArc(arc);

984}

985

986template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>

987bool GenericMinCostFlow<Graph, ArcFlowType, ArcScaledCostType>::IsArcValid(

988 ArcIndex arc) const {

989 return graph_->IsArcValid(arc);

990}

991

992template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>

993bool GenericMinCostFlow<Graph, ArcFlowType, ArcScaledCostType>::IsArcDirect(

994 ArcIndex arc) const {

995 DCHECK(IsArcValid(arc));

996 return arc >= 0;

997}

998

999

1000

1001

1002

1006 ::util::ReverseArcStaticGraph<uint16_t, int32_t>>;

1008 int64_t, int64_t>;

1009

1010

1012 ::util::ReverseArcStaticGraph<uint16_t, int32_t>,

1013 int16_t,

1014 int32_t>;

1015

1017 ArcIndex reserve_num_arcs) {

1018 if (reserve_num_nodes > 0) {

1019 node_supply_.reserve(reserve_num_nodes);

1020 }

1021 if (reserve_num_arcs > 0) {

1022 arc_tail_.reserve(reserve_num_arcs);

1023 arc_head_.reserve(reserve_num_arcs);

1024 arc_capacity_.reserve(reserve_num_arcs);

1025 arc_cost_.reserve(reserve_num_arcs);

1026 arc_permutation_.reserve(reserve_num_arcs);

1027 arc_flow_.reserve(reserve_num_arcs);

1028 }

1029}

1030

1032 ResizeNodeVectors(node);

1033 node_supply_[node] = supply;

1034}

1035

1039 ResizeNodeVectors(std::max(tail, head));

1040 const ArcIndex arc = arc_tail_.size();

1041 arc_tail_.push_back(tail);

1042 arc_head_.push_back(head);

1043 arc_capacity_.push_back(capacity);

1044 arc_cost_.push_back(unit_cost);

1045 return arc;

1046}

1047

1049 arc_capacity_[arc] = capacity;

1050}

1051

1053 return arc < arc_permutation_.size() ? arc_permutation_[arc] : arc;

1054}

1055

1057 SupplyAdjustment adjustment) {

1058 optimal_cost_ = 0;

1059 maximum_flow_ = 0;

1060 arc_flow_.clear();

1061 const NodeIndex num_nodes = node_supply_.size();

1062 const ArcIndex num_arcs = arc_capacity_.size();

1063 if (num_nodes == 0) return OPTIMAL;

1064

1065 int supply_node_count = 0, demand_node_count = 0;

1066 FlowQuantity total_supply = 0, total_demand = 0;

1067 for (NodeIndex node = 0; node < num_nodes; ++node) {

1068 if (node_supply_[node] > 0) {

1069 ++supply_node_count;

1070 total_supply += node_supply_[node];

1071 } else if (node_supply_[node] < 0) {

1072 ++demand_node_count;

1073 total_demand -= node_supply_[node];

1074 }

1075 }

1076 if (adjustment == DONT_ADJUST && total_supply != total_demand) {

1078 }

1079

1080

1081

1082

1083

1084

1085

1086

1087

1088

1089

1090

1091 const ArcIndex augmented_num_arcs =

1092 num_arcs + supply_node_count + demand_node_count;

1093 const NodeIndex source = num_nodes;

1094 const NodeIndex sink = num_nodes + 1;

1095 const NodeIndex augmented_num_nodes = num_nodes + 2;

1096

1097 Graph graph(augmented_num_nodes, augmented_num_arcs);

1098 for (ArcIndex arc = 0; arc < num_arcs; ++arc) {

1099 graph.AddArc(arc_tail_[arc], arc_head_[arc]);

1100 }

1101

1102 for (NodeIndex node = 0; node < num_nodes; ++node) {

1103 if (node_supply_[node] > 0) {

1104 graph.AddArc(source, node);

1105 } else if (node_supply_[node] < 0) {

1106 graph.AddArc(node, sink);

1107 }

1108 }

1109

1110 graph.Build(&arc_permutation_);

1111

1112 {

1113 GenericMaxFlow<Graph> max_flow(&graph, source, sink);

1115 for (arc = 0; arc < num_arcs; ++arc) {

1116 max_flow.SetArcCapacity(PermutedArc(arc), arc_capacity_[arc]);

1117 }

1118 for (NodeIndex node = 0; node < num_nodes; ++node) {

1119 if (node_supply_[node] != 0) {

1120 max_flow.SetArcCapacity(PermutedArc(arc), std::abs(node_supply_[node]));

1121 ++arc;

1122 }

1123 }

1124 CHECK_EQ(arc, augmented_num_arcs);

1125 if (!max_flow.Solve()) {

1126 LOG(ERROR) << "Max flow could not be computed.";

1127 switch (max_flow.status()) {

1131 LOG(ERROR)

1132 << "Max flow failed but claimed to have an optimal solution";

1133 ABSL_FALLTHROUGH_INTENDED;

1134 default:

1136 }

1137 }

1138 maximum_flow_ = max_flow.GetOptimalFlow();

1139 }

1140

1141 if (adjustment == DONT_ADJUST && maximum_flow_ != total_supply) {

1143 }

1144

1145 GenericMinCostFlow<Graph> min_cost_flow(&graph);

1147 for (arc = 0; arc < num_arcs; ++arc) {

1148 ArcIndex permuted_arc = PermutedArc(arc);

1149 min_cost_flow.SetArcUnitCost(permuted_arc, arc_cost_[arc]);

1150 min_cost_flow.SetArcCapacity(permuted_arc, arc_capacity_[arc]);

1151 }

1152 for (NodeIndex node = 0; node < num_nodes; ++node) {

1153 if (node_supply_[node] != 0) {

1154 ArcIndex permuted_arc = PermutedArc(arc);

1155 min_cost_flow.SetArcCapacity(permuted_arc, std::abs(node_supply_[node]));

1156 min_cost_flow.SetArcUnitCost(permuted_arc, 0);

1157 ++arc;

1158 }

1159 }

1160 min_cost_flow.SetNodeSupply(source, maximum_flow_);

1161 min_cost_flow.SetNodeSupply(sink, -maximum_flow_);

1162 min_cost_flow.SetCheckFeasibility(false);

1163 min_cost_flow.SetPriceScaling(scale_prices_);

1164

1165 arc_flow_.resize(num_arcs);

1166 if (min_cost_flow.Solve()) {

1167 optimal_cost_ = min_cost_flow.GetOptimalCost();

1168 for (arc = 0; arc < num_arcs; ++arc) {

1169 arc_flow_[arc] = min_cost_flow.Flow(PermutedArc(arc));

1170 }

1171 }

1172

1173 return min_cost_flow.status();

1174}

1175

1177 return optimal_cost_;

1178}

1179

1181 return maximum_flow_;

1182}

1183

1185 return arc_flow_[arc];

1186}

1187

1189 const {

1190 return node_supply_.size();

1191}

1192

1194 return arc_tail_.size();

1195}

1196

1198 return arc_tail_[arc];

1199}

1200

1202 return arc_head_[arc];

1203}

1204

1207 return arc_capacity_[arc];

1208}

1209

1211 return arc_cost_[arc];

1212}

1213

1216 return node_supply_[node];

1217}

1218

1219void SimpleMinCostFlow::ResizeNodeVectors(NodeIndex node) {

1220 if (node < node_supply_.size()) return;

1221 node_supply_.resize(node + 1);

1222}

1223

1224}

FlowQuantity Flow(ArcIndex arc) const

Definition min_cost_flow.cc:367

void SetArcUnitCost(ArcIndex arc, ArcScaledCostType unit_cost)

Definition min_cost_flow.cc:89

GenericMinCostFlow(const Graph *graph)

Definition min_cost_flow.cc:51

bool Solve()

Definition min_cost_flow.cc:451

void SetArcCapacity(ArcIndex arc, ArcFlowType new_capacity)

Definition min_cost_flow.cc:99

void SetNodeSupply(NodeIndex node, FlowQuantity supply)

Definition min_cost_flow.cc:79

const Graph * graph() const

CostValue UnitCost(ArcIndex arc) const

Definition min_cost_flow.cc:388

FlowQuantity Supply(NodeIndex node) const

Definition min_cost_flow.cc:396

Graph::NodeIndex NodeIndex

FlowQuantity Capacity(ArcIndex arc) const

Definition min_cost_flow.cc:378

CostValue GetOptimalCost()

Definition min_cost_flow.cc:482

NodeIndex Tail(ArcIndex arc) const

Definition min_cost_flow.cc:1197

CostValue OptimalCost() const

Definition min_cost_flow.cc:1176

ArcIndex NumArcs() const

Definition min_cost_flow.cc:1193

void SetNodeSupply(NodeIndex node, FlowQuantity supply)

Definition min_cost_flow.cc:1031

void SetArcCapacity(ArcIndex arc, FlowQuantity capacity)

Definition min_cost_flow.cc:1048

FlowQuantity Flow(ArcIndex arc) const

Definition min_cost_flow.cc:1184

FlowQuantity Capacity(ArcIndex arc) const

Definition min_cost_flow.cc:1205

CostValue UnitCost(ArcIndex arc) const

Definition min_cost_flow.cc:1210

NodeIndex NumNodes() const

Definition min_cost_flow.cc:1188

ArcIndex AddArcWithCapacityAndUnitCost(NodeIndex tail, NodeIndex head, FlowQuantity capacity, CostValue unit_cost)

Definition min_cost_flow.cc:1036

FlowQuantity MaximumFlow() const

Definition min_cost_flow.cc:1180

NodeIndex Head(ArcIndex arc) const

Definition min_cost_flow.cc:1201

SimpleMinCostFlow(NodeIndex reserve_num_nodes=0, ArcIndex reserve_num_arcs=0)

Definition min_cost_flow.cc:1016

FlowQuantity Supply(NodeIndex node) const

Definition min_cost_flow.cc:1214

static constexpr int32_t kNilArc

static constexpr bool kHasNegativeReverseArcs

ABSL_FLAG(int64_t, min_cost_flow_alpha, 5, "Divide factor for epsilon at each refine step.")

int64_t CapAdd(int64_t x, int64_t y)

int64_t CapSub(int64_t x, int64_t y)

BlossomGraph::CostValue CostValue

BlossomGraph::NodeIndex NodeIndex

int64_t CapProd(int64_t x, int64_t y)

util::ReverseArcStaticGraph Graph

#define IF_STATS_ENABLED(instructions)

#define SCOPED_TIME_STAT(stats)