path: root/src/util/bipartite.rs

/*
 * This module deals with graph algorithm in complete bipartite
 * graphs. It is used in layout.rs to build the partition to node
 * assignation.
 * */

use rand::prelude::SliceRandom;
use std::cmp::{max, min};
use std::collections::VecDeque;

//Graph data structure for the flow algorithm.
#[derive(Clone, Copy, Debug)]
struct EdgeFlow {
	c: i32,
	flow: i32,
	v: usize,
	rev: usize,
}

//Graph data structure for the detection of positive cycles.
#[derive(Clone, Copy, Debug)]
struct WeightedEdge {
	w: i32,
	u: usize,
	v: usize,
}

/* This function takes two matchings (old_match and new_match) in a
 * complete bipartite graph. It returns a matching that has the
 * same degree as new_match at every vertex, and that is as close
 * as possible to old_match.
 * */
pub fn optimize_matching(
	old_match: &Vec<Vec<usize>>,
	new_match: &Vec<Vec<usize>>,
	nb_right: usize,
) -> Vec<Vec<usize>> {
	let nb_left = old_match.len();
	let ed = WeightedEdge { w: -1, u: 0, v: 0 };
	let mut edge_vec = vec![ed; nb_left * nb_right];

	//We build the complete bipartite graph structure, represented
	//by the list of all edges.
	for i in 0..nb_left {
		for j in 0..nb_right {
			edge_vec[i * nb_right + j].u = i;
			edge_vec[i * nb_right + j].v = nb_left + j;
		}
	}

	for i in 0..edge_vec.len() {
		//We add the old matchings
		if old_match[edge_vec[i].u].contains(&(edge_vec[i].v - nb_left)) {
			edge_vec[i].w *= -1;
		}
		//We add the new matchings
		if new_match[edge_vec[i].u].contains(&(edge_vec[i].v - nb_left)) {
			(edge_vec[i].u, edge_vec[i].v) = (edge_vec[i].v, edge_vec[i].u);
			edge_vec[i].w *= -1;
		}
	}
	//Now edge_vec is a graph where edges are oriented LR if we
	//can add them to new_match, and RL otherwise. If
	//adding/removing them makes the matching closer to old_match
	//they have weight 1; and -1 otherwise.

	//We shuffle the edge list so that there is no bias depending in
	//partitions/zone label in the triplet dispersion
	let mut rng = rand::thread_rng();
	edge_vec.shuffle(&mut rng);

	//Discovering and flipping a cycle with positive weight in this
	//graph will make the matching closer to old_match.
	//We use Bellman Ford algorithm to discover positive cycles
	loop {
		if let Some(cycle) = positive_cycle(&edge_vec, nb_left, nb_right) {
			for i in cycle {
				//We flip the edges of the cycle.
				(edge_vec[i].u, edge_vec[i].v) = (edge_vec[i].v, edge_vec[i].u);
				edge_vec[i].w *= -1;
			}
		} else {
			//If there is no cycle, we return the optimal matching.
			break;
		}
	}

	//The optimal matching is build from the graph structure.
	let mut matching = vec![Vec::<usize>::new(); nb_left];
	for e in edge_vec {
		if e.u > e.v {
			matching[e.v].push(e.u - nb_left);
		}
	}
	matching
}

//This function finds a positive cycle in a bipartite wieghted graph.
fn positive_cycle(
	edge_vec: &Vec<WeightedEdge>,
	nb_left: usize,
	nb_right: usize,
) -> Option<Vec<usize>> {
	let nb_side_min = min(nb_left, nb_right);
	let nb_vertices = nb_left + nb_right;
	let weight_lowerbound = -((nb_left + nb_right) as i32) - 1;
	let mut accessed = vec![false; nb_left];

	//We try to find a positive cycle accessible from the left
	//vertex i.
	for i in 0..nb_left {
		if accessed[i] {
			continue;
		}
		let mut weight = vec![weight_lowerbound; nb_vertices];
		let mut prev = vec![edge_vec.len(); nb_vertices];
		weight[i] = 0;
		//We compute largest weighted paths from i.
		//Since the graph is bipartite, any simple cycle has length
		//at most 2*nb_side_min. In the general Bellman-Ford
		//algorithm, the bound here is the number of vertices. Since
		//the number of partitions can be much larger than the
		//number of nodes, we optimize that.
		for _ in 0..(2 * nb_side_min) {
			for j in 0..edge_vec.len() {
				let e = edge_vec[j];
				if weight[e.v] < weight[e.u] + e.w {
					weight[e.v] = weight[e.u] + e.w;
					prev[e.v] = j;
				}
			}
		}
		//We update the accessed table
		for i in 0..nb_left {
			if weight[i] > weight_lowerbound {
				accessed[i] = true;
			}
		}
		//We detect positive cycle
		for e in edge_vec {
			if weight[e.v] < weight[e.u] + e.w {
				//it means e is on a path branching from a positive cycle
				let mut was_seen = vec![false; nb_vertices];
				let mut curr = e.u;
				//We track back with prev until we reach the cycle.
				while !was_seen[curr] {
					was_seen[curr] = true;
					curr = edge_vec[prev[curr]].u;
				}
				//Now curr is on the cycle. We collect the edges ids.
				let mut cycle = Vec::<usize>::new();
				cycle.push(prev[curr]);
				let mut cycle_vert = edge_vec[prev[curr]].u;
				while cycle_vert != curr {
					cycle.push(prev[cycle_vert]);
					cycle_vert = edge_vec[prev[cycle_vert]].u;
				}

				return Some(cycle);
			}
		}
	}

	None
}

// This function takes two arrays of capacity and computes the
// maximal matching in the complete bipartite graph such that the
// left vertex i is matched to left_cap_vec[i] right vertices, and
// the right vertex j is matched to right_cap_vec[j] left vertices.
// To do so, we use Dinic's maximum flow algorithm.
pub fn dinic_compute_matching(left_cap_vec: Vec<u32>, right_cap_vec: Vec<u32>) -> Vec<Vec<usize>> {
	let mut graph = Vec::<Vec<EdgeFlow>>::new();
	let ed = EdgeFlow {
		c: 0,
		flow: 0,
		v: 0,
		rev: 0,
	};

	// 0 will be the source
	graph.push(vec![ed; left_cap_vec.len()]);
	for i in 0..left_cap_vec.len() {
		graph[0][i].c = left_cap_vec[i] as i32;
		graph[0][i].v = i + 2;
		graph[0][i].rev = 0;
	}

	//1 will be the sink
	graph.push(vec![ed; right_cap_vec.len()]);
	for i in 0..right_cap_vec.len() {
		graph[1][i].c = right_cap_vec[i] as i32;
		graph[1][i].v = i + 2 + left_cap_vec.len();
		graph[1][i].rev = 0;
	}

	//we add left vertices
	for i in 0..left_cap_vec.len() {
		graph.push(vec![ed; 1 + right_cap_vec.len()]);
		graph[i + 2][0].c = 0; //directed
		graph[i + 2][0].v = 0;
		graph[i + 2][0].rev = i;

		for j in 0..right_cap_vec.len() {
			graph[i + 2][j + 1].c = 1;
			graph[i + 2][j + 1].v = 2 + left_cap_vec.len() + j;
			graph[i + 2][j + 1].rev = i + 1;
		}
	}

	//we add right vertices
	for i in 0..right_cap_vec.len() {
		let lft_ln = left_cap_vec.len();
		graph.push(vec![ed; 1 + lft_ln]);
		graph[i + lft_ln + 2][0].c = graph[1][i].c;
		graph[i + lft_ln + 2][0].v = 1;
		graph[i + lft_ln + 2][0].rev = i;

		for j in 0..left_cap_vec.len() {
			graph[i + 2 + lft_ln][j + 1].c = 0; //directed
			graph[i + 2 + lft_ln][j + 1].v = j + 2;
			graph[i + 2 + lft_ln][j + 1].rev = i + 1;
		}
	}

	//To ensure the dispersion of the triplets generated by the
	//assignation, we shuffle the neighbours of the nodes. Hence,
	//left vertices do not consider the right ones in the same order.
	let mut rng = rand::thread_rng();
	for i in 0..graph.len() {
		graph[i].shuffle(&mut rng);
		//We need to update the ids of the reverse edges.
		for j in 0..graph[i].len() {
			let target_v = graph[i][j].v;
			let target_rev = graph[i][j].rev;
			graph[target_v][target_rev].rev = j;
		}
	}

	let nb_vertices = graph.len();

	//We run Dinic's max flow algorithm
	loop {
		//We build the level array from Dinic's algorithm.
		let mut level = vec![-1; nb_vertices];

		let mut fifo = VecDeque::new();
		fifo.push_back((0, 0));
		while !fifo.is_empty() {
			if let Some((id, lvl)) = fifo.pop_front() {
				if level[id] == -1 {
					level[id] = lvl;
					for e in graph[id].iter() {
						if e.c - e.flow > 0 {
							fifo.push_back((e.v, lvl + 1));
						}
					}
				}
			}
		}
		if level[1] == -1 {
			//There is no residual flow
			break;
		}

		//Now we run DFS respecting the level array
		let mut next_nbd = vec![0; nb_vertices];
		let mut lifo = VecDeque::new();

		let flow_upper_bound;
		if let Some(x) = left_cap_vec.iter().max() {
			flow_upper_bound = *x as i32;
		} else {
			flow_upper_bound = 0;
			assert!(false);
		}

		lifo.push_back((0, flow_upper_bound));

		loop {
			if let Some((id_tmp, f_tmp)) = lifo.back() {
				let id = *id_tmp;
				let f = *f_tmp;
				if id == 1 {
					//The DFS reached the sink, we can add a
					//residual flow.
					lifo.pop_back();
					while !lifo.is_empty() {
						if let Some((id, _)) = lifo.pop_back() {
							let nbd = next_nbd[id];
							graph[id][nbd].flow += f;
							let id_v = graph[id][nbd].v;
							let nbd_v = graph[id][nbd].rev;
							graph[id_v][nbd_v].flow -= f;
						}
					}
					lifo.push_back((0, flow_upper_bound));
					continue;
				}
				//else we did not reach the sink
				let nbd = next_nbd[id];
				if nbd >= graph[id].len() {
					//There is nothing to explore from id anymore
					lifo.pop_back();
					if let Some((parent, _)) = lifo.back() {
						next_nbd[*parent] += 1;
					}
					continue;
				}
				//else we can try to send flow from id to its nbd
				let new_flow = min(f, graph[id][nbd].c - graph[id][nbd].flow);
				if level[graph[id][nbd].v] <= level[id] || new_flow == 0 {
					//We cannot send flow to nbd.
					next_nbd[id] += 1;
					continue;
				}
				//otherwise, we send flow to nbd.
				lifo.push_back((graph[id][nbd].v, new_flow));
			} else {
				break;
			}
		}
	}

	//We return the association
	let assoc_table = (0..left_cap_vec.len())
		.map(|id| {
			graph[id + 2]
				.iter()
				.filter(|e| e.flow > 0)
				.map(|e| e.v - 2 - left_cap_vec.len())
				.collect()
		})
		.collect();

	//consistency check

	//it is a flow
	for i in 3..graph.len() {
		assert!(graph[i].iter().map(|e| e.flow).sum::<i32>() == 0);
		for e in graph[i].iter() {
			assert!(e.flow + graph[e.v][e.rev].flow == 0);
		}
	}

	//it solves the matching problem
	for i in 0..left_cap_vec.len() {
		assert!(left_cap_vec[i] as i32 == graph[i + 2].iter().map(|e| max(0, e.flow)).sum::<i32>());
	}
	for i in 0..right_cap_vec.len() {
		assert!(
			right_cap_vec[i] as i32
				== graph[i + 2 + left_cap_vec.len()]
					.iter()
					.map(|e| max(0, e.flow))
					.sum::<i32>()
		);
	}

	assoc_table
}

#[cfg(test)]
mod tests {
	use super::*;

	#[test]
	fn test_flow() {
		let left_vec = vec![3; 8];
		let right_vec = vec![0, 4, 8, 4, 8];
		//There are asserts in the function that computes the flow
		let _ = dinic_compute_matching(left_vec, right_vec);
	}

	//maybe add tests relative to the matching optilization ?
}