use rand::prelude::{SeedableRng, SliceRandom};
use std::cmp::{max, min};
use std::collections::HashMap;
use std::collections::VecDeque;
#[derive(Clone, Copy, Debug, PartialEq, Eq, Hash)]
pub enum Vertex {
Source,
Pup(usize), Pdown(usize), PZ(usize, usize), N(usize), Sink,
}
#[derive(Clone, Copy, Debug)]
pub struct FlowEdge {
cap: u64, flow: i64, dest: usize, rev: usize, }
#[derive(Clone, Copy, Debug)]
pub struct WeightedEdge {
w: i64, dest: usize,
}
pub trait Edge: Clone + Copy {}
impl Edge for FlowEdge {}
impl Edge for WeightedEdge {}
pub struct Graph<E: Edge> {
vertex_to_id: HashMap<Vertex, usize>,
id_to_vertex: Vec<Vertex>,
graph: Vec<Vec<E>>,
}
pub type CostFunction = HashMap<(Vertex, Vertex), i64>;
impl<E: Edge> Graph<E> {
pub fn new(vertices: &[Vertex]) -> Self {
let mut map = HashMap::<Vertex, usize>::new();
for (i, vert) in vertices.iter().enumerate() {
map.insert(*vert, i);
}
Graph::<E> {
vertex_to_id: map,
id_to_vertex: vertices.to_vec(),
graph: vec![Vec::<E>::new(); vertices.len()],
}
}
fn get_vertex_id(&self, v: &Vertex) -> Result<usize, String> {
self.vertex_to_id
.get(v)
.cloned()
.ok_or_else(|| format!("The graph does not contain vertex {:?}", v))
}
}
impl Graph<FlowEdge> {
pub fn add_edge(&mut self, u: Vertex, v: Vertex, c: u64) -> Result<(), String> {
let idu = self.get_vertex_id(&u)?;
let idv = self.get_vertex_id(&v)?;
if idu == idv {
return Err("Cannot add edge from vertex to itself in flow graph".into());
}
let rev_u = self.graph[idu].len();
let rev_v = self.graph[idv].len();
self.graph[idu].push(FlowEdge {
cap: c,
dest: idv,
flow: 0,
rev: rev_v,
});
self.graph[idv].push(FlowEdge {
cap: 0,
dest: idu,
flow: 0,
rev: rev_u,
});
Ok(())
}
pub fn get_positive_flow_from(&self, v: Vertex) -> Result<Vec<Vertex>, String> {
let idv = self.get_vertex_id(&v)?;
let mut result = Vec::<Vertex>::new();
for edge in self.graph[idv].iter() {
if edge.flow > 0 {
result.push(self.id_to_vertex[edge.dest]);
}
}
Ok(result)
}
pub fn get_inflow(&self, v: Vertex) -> Result<i64, String> {
let idv = self.get_vertex_id(&v)?;
let mut result = 0;
for edge in self.graph[idv].iter() {
result += max(0, self.graph[edge.dest][edge.rev].flow);
}
Ok(result)
}
pub fn get_outflow(&self, v: Vertex) -> Result<i64, String> {
let idv = self.get_vertex_id(&v)?;
let mut result = 0;
for edge in self.graph[idv].iter() {
result += max(0, edge.flow);
}
Ok(result)
}
pub fn get_flow_value(&mut self) -> Result<i64, String> {
self.get_outflow(Vertex::Source)
}
fn shuffle_edges(&mut self) {
let mut rng = rand::rngs::StdRng::from_seed([0x12u8; 32]);
for i in 0..self.graph.len() {
self.graph[i].shuffle(&mut rng);
for j in 0..self.graph[i].len() {
let target_v = self.graph[i][j].dest;
let target_rev = self.graph[i][j].rev;
self.graph[target_v][target_rev].rev = j;
}
}
}
pub fn flow_upper_bound(&self) -> Result<u64, String> {
let idsource = self.get_vertex_id(&Vertex::Source)?;
let mut flow_upper_bound = 0;
for edge in self.graph[idsource].iter() {
flow_upper_bound += edge.cap;
}
Ok(flow_upper_bound)
}
pub fn compute_maximal_flow(&mut self) -> Result<(), String> {
let idsource = self.get_vertex_id(&Vertex::Source)?;
let idsink = self.get_vertex_id(&Vertex::Sink)?;
let nb_vertices = self.graph.len();
let flow_upper_bound = self.flow_upper_bound()?;
self.shuffle_edges();
loop {
let mut level = vec![None; nb_vertices];
let mut fifo = VecDeque::new();
fifo.push_back((idsource, 0));
while let Some((id, lvl)) = fifo.pop_front() {
if level[id].is_none() {
level[id] = Some(lvl);
for edge in self.graph[id].iter() {
if edge.cap as i64 - edge.flow > 0 {
fifo.push_back((edge.dest, lvl + 1));
}
}
}
}
if level[idsink].is_none() {
break;
}
let mut next_nbd = vec![0; nb_vertices];
let mut lifo = Vec::new();
lifo.push((idsource, flow_upper_bound));
while let Some((id, f)) = lifo.last().cloned() {
if id == idsink {
lifo.pop();
while let Some((id, _)) = lifo.pop() {
let nbd = next_nbd[id];
self.graph[id][nbd].flow += f as i64;
let id_rev = self.graph[id][nbd].dest;
let nbd_rev = self.graph[id][nbd].rev;
self.graph[id_rev][nbd_rev].flow -= f as i64;
}
lifo.push((idsource, flow_upper_bound));
continue;
}
let nbd = next_nbd[id];
if nbd >= self.graph[id].len() {
lifo.pop();
if let Some((parent, _)) = lifo.last() {
next_nbd[*parent] += 1;
}
continue;
}
let new_flow = min(
f as i64,
self.graph[id][nbd].cap as i64 - self.graph[id][nbd].flow,
) as u64;
if new_flow == 0 {
next_nbd[id] += 1;
continue;
}
if let (Some(lvldest), Some(lvlid)) = (level[self.graph[id][nbd].dest], level[id]) {
if lvldest <= lvlid {
next_nbd[id] += 1;
continue;
}
}
lifo.push((self.graph[id][nbd].dest, new_flow));
}
}
Ok(())
}
pub fn optimize_flow_with_cost(
&mut self,
cost: &CostFunction,
path_length: usize,
) -> Result<(), String> {
let mut gf = self.build_cost_graph(cost)?;
let mut cycles = gf.list_negative_cycles(path_length);
while !cycles.is_empty() {
for c in cycles.iter() {
for i in 0..c.len() {
let idu = self.vertex_to_id[&c[i]];
let idv = self.vertex_to_id[&c[(i + 1) % c.len()]];
for j in 0..self.graph[idu].len() {
let edge = self.graph[idu][j];
if edge.dest == idv {
self.graph[idu][j].flow += 1;
self.graph[idv][edge.rev].flow -= 1;
break;
}
}
}
}
gf = self.build_cost_graph(cost)?;
cycles = gf.list_negative_cycles(path_length);
}
Ok(())
}
fn build_cost_graph(&self, cost: &CostFunction) -> Result<Graph<WeightedEdge>, String> {
let mut g = Graph::<WeightedEdge>::new(&self.id_to_vertex);
let nb_vertices = self.id_to_vertex.len();
for i in 0..nb_vertices {
for edge in self.graph[i].iter() {
if edge.cap as i64 - edge.flow > 0 {
let u = self.id_to_vertex[i];
let v = self.id_to_vertex[edge.dest];
if cost.contains_key(&(u, v)) {
g.add_edge(u, v, cost[&(u, v)])?;
} else if cost.contains_key(&(v, u)) {
g.add_edge(u, v, -cost[&(v, u)])?;
} else {
g.add_edge(u, v, 0)?;
}
}
}
}
Ok(g)
}
}
impl Graph<WeightedEdge> {
pub fn add_edge(&mut self, u: Vertex, v: Vertex, w: i64) -> Result<(), String> {
let idu = self.get_vertex_id(&u)?;
let idv = self.get_vertex_id(&v)?;
self.graph[idu].push(WeightedEdge { w, dest: idv });
Ok(())
}
fn list_negative_cycles(&self, path_length: usize) -> Vec<Vec<Vertex>> {
let nb_vertices = self.graph.len();
let mut distance = vec![0; nb_vertices];
let mut prev = vec![None; nb_vertices];
for _ in 0..path_length + 1 {
for id in 0..nb_vertices {
for e in self.graph[id].iter() {
if distance[id] + e.w < distance[e.dest] {
distance[e.dest] = distance[id] + e.w;
prev[e.dest] = Some(id);
}
}
}
}
let cycles_prev = cycles_of_1_forest(&prev);
return cycles_prev
.iter()
.map(|cycle| {
cycle
.iter()
.rev()
.map(|id| self.id_to_vertex[*id])
.collect()
})
.collect();
}
}
fn cycles_of_1_forest(forest: &[Option<usize>]) -> Vec<Vec<usize>> {
let mut cycles = Vec::<Vec<usize>>::new();
let mut time_of_discovery = vec![None; forest.len()];
for t in 0..forest.len() {
let mut id = t;
while time_of_discovery[id].is_none() {
time_of_discovery[id] = Some(t);
if let Some(i) = forest[id] {
id = i;
} else {
break;
}
}
if forest[id].is_some() && time_of_discovery[id] == Some(t) {
let mut cy = vec![id; 1];
let mut id2 = id;
while let Some(id_next) = forest[id2] {
id2 = id_next;
if id2 != id {
cy.push(id2);
} else {
break;
}
}
cycles.push(cy);
}
}
cycles
}