From 5521e9c9945281b45fc390137ce77559e3e8c8d9 Mon Sep 17 00:00:00 2001 From: Hubert Hirtz Date: Thu, 17 Mar 2022 09:06:27 +0100 Subject: [PATCH] WIP: rework Partitioner trait family This simplifies the family of traits from before into one "Partition" trait. --- src/algorithms.rs | 74 +- src/algorithms/ckk.rs | 41 +- src/algorithms/fiduccia_mattheyses.rs | 124 ++- src/algorithms/graph_growth.rs | 96 +- src/algorithms/greedy.rs | 39 +- src/algorithms/hilbert_curve.rs | 86 +- src/algorithms/k_means.rs | 123 ++- src/algorithms/kernighan_lin.rs | 123 ++- src/algorithms/kk.rs | 40 +- src/algorithms/multi_jagged.rs | 93 +- src/algorithms/recursive_bisection.rs | 146 ++- src/algorithms/vn/best.rs | 77 +- src/algorithms/vn/mod.rs | 2 + src/algorithms/z_curve.rs | 90 +- src/imbalance.rs | 8 +- src/lib.rs | 1327 +------------------------ 16 files changed, 1091 insertions(+), 1398 deletions(-) diff --git a/src/algorithms.rs b/src/algorithms.rs index 45d1df40..37b9795b 100644 --- a/src/algorithms.rs +++ b/src/algorithms.rs @@ -1,12 +1,62 @@ -pub mod ckk; -pub mod fiduccia_mattheyses; -pub mod graph_growth; -pub mod greedy; -pub mod hilbert_curve; -pub mod k_means; -pub mod kernighan_lin; -pub mod kk; -pub mod multi_jagged; -pub mod recursive_bisection; -pub mod vn; -pub mod z_curve; +mod ckk; +mod fiduccia_mattheyses; +mod graph_growth; +mod greedy; +mod hilbert_curve; +mod k_means; +mod kernighan_lin; +mod kk; +mod multi_jagged; +mod recursive_bisection; +mod vn; +mod z_curve; + +pub use ckk::CompleteKarmarkarKarp; +pub use fiduccia_mattheyses::FiducciaMattheyses; +pub use graph_growth::GraphGrowth; +pub use greedy::Greedy; +pub use hilbert_curve::HilbertCurve; +pub use k_means::KMeans; +pub use kernighan_lin::KernighanLin; +pub use kk::KarmarkarKarp; +pub use multi_jagged::MultiJagged; +pub use recursive_bisection::Rcb; +pub use recursive_bisection::Rib; +pub use vn::VnBest; +pub use z_curve::ZCurve; + +/// Map mesh elements to parts randomly. +/// +/// # Example +/// +/// ```rust +/// use coupe::Partitioner; +/// use coupe::Point2D; +/// use coupe::Random; +/// use rand; +/// +/// let num_parts = 3; +/// let points: &[Point2D] = &[Point2D::new(1., 0.), Point2D::new(0., 1.)]; +/// +/// let r = Random::new(rand::thread_rng(), num_parts); +/// r.partition(points, &[1., 1.]); +/// ``` +pub struct Random { + pub rng: R, + pub part_count: usize, +} + +impl crate::Partition<()> for Random +where + R: rand::Rng, +{ + type Metadata = (); + type Error = std::convert::Infallible; + + fn partition(&mut self, part_ids: &mut [usize], _: ()) -> Result { + for part_id in part_ids { + *part_id = self.rng.gen_range(0..self.part_count); + } + Ok(()) + } +} diff --git a/src/algorithms/ckk.rs b/src/algorithms/ckk.rs index a3cec67c..cf6b6bdb 100644 --- a/src/algorithms/ckk.rs +++ b/src/algorithms/ckk.rs @@ -92,7 +92,7 @@ where false } -pub fn ckk_bipart(partition: &mut [usize], weights: I, tolerance: f64) -> bool +fn ckk_bipart(partition: &mut [usize], weights: I, tolerance: f64) -> bool where I: IntoIterator, T: Sum + Add + Sub, @@ -114,4 +114,43 @@ where ckk_bipart_rec(partition, &mut weights, tolerance, &mut steps) } +/// The exact variant of the [Karmarkar-Karp][crate::KarmarkarKarp] algorithm. +/// +/// # Example +/// +/// ```rust +/// use coupe::Partitioner; +/// use coupe::Point2D; +/// use coupe::CompleteKarmarkarKarp; +/// +/// let tolerance = 0.1; +/// let points: &[Point2D] = &[Point2D::new(1., 0.), Point2D::new(0., 1.)]; +/// +/// let ckk = CompleteKarmarkarKarp::new(tolerance); +/// ckk.partition(points, &[1., 1.]); +/// ``` +pub struct CompleteKarmarkarKarp { + pub tolerance: f64, +} + +impl crate::Partition for CompleteKarmarkarKarp +where + W: IntoIterator, + W::Item: Sum + Add + Sub, + W::Item: FromPrimitive + ToPrimitive, + W::Item: Ord + Default + Copy, +{ + type Metadata = (); + type Error = std::convert::Infallible; + + fn partition( + &mut self, + part_ids: &mut [usize], + weights: W, + ) -> Result { + ckk_bipart(part_ids, weights, self.tolerance); + Ok(()) + } +} + // TODO tests diff --git a/src/algorithms/fiduccia_mattheyses.rs b/src/algorithms/fiduccia_mattheyses.rs index 088313ba..9db5c484 100644 --- a/src/algorithms/fiduccia_mattheyses.rs +++ b/src/algorithms/fiduccia_mattheyses.rs @@ -1,12 +1,11 @@ use itertools::Itertools; use sprs::CsMatView; -use crate::partition::Partition; -use crate::PointND; use std::collections::HashMap; -pub fn fiduccia_mattheyses<'a, const D: usize>( - initial_partition: &mut Partition<'a, PointND, f64>, +fn fiduccia_mattheyses<'a>( + part_ids: &mut [usize], + weights: &[f64], adjacency: CsMatView, max_passes: impl Into>, max_flips_per_pass: impl Into>, @@ -16,12 +15,11 @@ pub fn fiduccia_mattheyses<'a, const D: usize>( let max_passes = max_passes.into(); let max_flips_per_pass = max_flips_per_pass.into(); let max_imbalance_per_flip = max_imbalance_per_flip.into(); - let (_points, weights, ids) = initial_partition.as_raw_mut(); fiduccia_mattheyses_impl( + part_ids, weights, - adjacency.view(), - ids, + adjacency, max_passes, max_flips_per_pass, max_imbalance_per_flip, @@ -30,9 +28,9 @@ pub fn fiduccia_mattheyses<'a, const D: usize>( } fn fiduccia_mattheyses_impl( + initial_partition: &mut [usize], weights: &[f64], adjacency: CsMatView, - initial_partition: &mut [usize], max_passes: Option, max_flips_per_pass: Option, max_imbalance_per_flip: Option, @@ -250,3 +248,113 @@ fn fiduccia_mattheyses_impl( println!("final cut size: {}", new_cut_size); } + +/// FiducciaMattheyses +/// +/// An implementation of the Fiduccia Mattheyses topologic algorithm +/// for graph partitioning. This implementation is an extension of the +/// original algorithm to handle partitioning into more than two parts. +/// +/// This algorithm repeats an iterative pass during which a set of graph nodes are assigned to +/// a new part, reducing the overall cutsize of the partition. As opposed to the +/// Kernighan-Lin algorithm, during each pass iteration, only one node is flipped at a time. +/// The algorithm thus does not preserve partition weights balance and may produce an unbalanced +/// partition. +/// +/// Original algorithm from "A Linear-Time Heuristic for Improving Network Partitions" +/// by C.M. Fiduccia and R.M. Mattheyses. +/// +/// # Example +/// +/// ```rust +/// use coupe::Point2D; +/// use coupe::TopologicPartitionImprover; +/// use coupe::partition::Partition; +/// use sprs::CsMat; +/// +/// // swap +/// // 0 1 0 1 +/// // +--+--+--+ +/// // | | | | +/// // +--+--+--+ +/// // 0 0 1 1 +/// let points = vec![ +/// Point2D::new(0., 0.), +/// Point2D::new(1., 0.), +/// Point2D::new(2., 0.), +/// Point2D::new(3., 0.), +/// Point2D::new(0., 1.), +/// Point2D::new(1., 1.), +/// Point2D::new(2., 1.), +/// Point2D::new(3., 1.), +/// ]; +/// +/// let ids = vec![0, 0, 1, 1, 0, 1, 0, 1]; +/// let weights = vec![1.; 8]; +/// +/// let mut partition = Partition::from_ids(&points, &weights, ids); +/// +/// let mut adjacency = CsMat::empty(sprs::CSR, 8); +/// adjacency.reserve_outer_dim(8); +/// eprintln!("shape: {:?}", adjacency.shape()); +/// adjacency.insert(0, 1, 1.); +/// adjacency.insert(1, 2, 1.); +/// adjacency.insert(2, 3, 1.); +/// adjacency.insert(4, 5, 1.); +/// adjacency.insert(5, 6, 1.); +/// adjacency.insert(6, 7, 1.); +/// adjacency.insert(0, 4, 1.); +/// adjacency.insert(1, 5, 1.); +/// adjacency.insert(2, 6, 1.); +/// adjacency.insert(3, 7, 1.); +/// +/// // symmetry +/// adjacency.insert(1, 0, 1.); +/// adjacency.insert(2, 1, 1.); +/// adjacency.insert(3, 2, 1.); +/// adjacency.insert(5, 4, 1.); +/// adjacency.insert(6, 5, 1.); +/// adjacency.insert(7, 6, 1.); +/// adjacency.insert(4, 0, 1.); +/// adjacency.insert(5, 1, 1.); +/// adjacency.insert(6, 2, 1.); +/// adjacency.insert(7, 3, 1.); +/// +/// // 1 iteration +/// let algo = coupe::FiducciaMattheyses::new(None, None, None, 1); +/// +/// let partition = algo.improve_partition(partition, adjacency.view()); +/// +/// let new_ids = partition.into_ids(); +/// assert_eq!(new_ids[5], 0); +/// assert_eq!(new_ids[6], 1); +/// ``` +#[derive(Debug, Clone, Copy)] +pub struct FiducciaMattheyses { + pub max_passes: Option, + pub max_flips_per_pass: Option, + pub max_imbalance_per_flip: Option, + pub max_bad_move_in_a_row: usize, +} + +impl<'a> crate::Partition<(CsMatView<'a, f64>, &'a [f64])> for FiducciaMattheyses { + type Metadata = (); + type Error = std::convert::Infallible; + + fn partition( + &mut self, + part_ids: &mut [usize], + (adjacency, weights): (CsMatView, &'a [f64]), + ) -> Result { + fiduccia_mattheyses( + part_ids, + weights, + adjacency, + self.max_passes, + self.max_flips_per_pass, + self.max_imbalance_per_flip, + self.max_bad_move_in_a_row, + ); + Ok(()) + } +} diff --git a/src/algorithms/graph_growth.rs b/src/algorithms/graph_growth.rs index b9872195..6538adb0 100644 --- a/src/algorithms/graph_growth.rs +++ b/src/algorithms/graph_growth.rs @@ -1,12 +1,16 @@ use rand::seq::SliceRandom; use sprs::CsMatView; -pub fn graph_growth(weights: &[f64], adjacency: CsMatView, num_parts: usize) -> Vec { +fn graph_growth( + initial_ids: &mut [usize], + weights: &[f64], + adjacency: CsMatView, + num_parts: usize, +) { let (shape_x, shape_y) = adjacency.shape(); assert_eq!(shape_x, shape_y); assert_eq!(weights.len(), shape_x); - let mut initial_ids = vec![0; weights.len()]; // let total_weight = weights.iter().sum::(); // let weight_per_part = total_weight / num_parts as f64; let max_expansion_per_pass = 20; @@ -54,6 +58,92 @@ pub fn graph_growth(weights: &[f64], adjacency: CsMatView, num_parts: usize } } } +} + +/// Graph Growth algorithm +/// +/// A topologic algorithm that generates a partition from a topologic mesh. +/// Given a number k of parts, the algorithm selects k nodes randomly and assigns them to a different part. +/// Then, at each iteration, each part is expanded to neighbor nodes that are not yet assigned to a part +/// +/// # Example +/// +/// ```rust +/// use coupe::Point2D; +/// use coupe::TopologicPartitioner; +/// use coupe::partition::Partition; +/// use sprs::CsMat; +/// +/// // +--+--+--+ +/// // | | | | +/// // +--+--+--+ +/// +/// let points = vec![ +/// Point2D::new(0., 0.), +/// Point2D::new(1., 0.), +/// Point2D::new(2., 0.), +/// Point2D::new(3., 0.), +/// Point2D::new(0., 1.), +/// Point2D::new(1., 1.), +/// Point2D::new(2., 1.), +/// Point2D::new(3., 1.), +/// ]; +/// +/// let weights = vec![1.; 8]; +/// +/// let mut adjacency = CsMat::empty(sprs::CSR, 8); +/// adjacency.reserve_outer_dim(8); +/// eprintln!("shape: {:?}", adjacency.shape()); +/// adjacency.insert(0, 1, 1.); +/// adjacency.insert(1, 2, 1.); +/// adjacency.insert(2, 3, 1.); +/// adjacency.insert(4, 5, 1.); +/// adjacency.insert(5, 6, 1.); +/// adjacency.insert(6, 7, 1.); +/// adjacency.insert(0, 4, 1.); +/// adjacency.insert(1, 5, 1.); +/// adjacency.insert(2, 6, 1.); +/// adjacency.insert(3, 7, 1.); +/// +/// // symmetry +/// adjacency.insert(1, 0, 1.); +/// adjacency.insert(2, 1, 1.); +/// adjacency.insert(3, 2, 1.); +/// adjacency.insert(5, 4, 1.); +/// adjacency.insert(6, 5, 1.); +/// adjacency.insert(7, 6, 1.); +/// adjacency.insert(4, 0, 1.); +/// adjacency.insert(5, 1, 1.); +/// adjacency.insert(6, 2, 1.); +/// adjacency.insert(7, 3, 1.); +/// +/// let gg = coupe::GraphGrowth::new(2); +/// +/// let _partition = gg.partition(points.as_slice(), &weights, adjacency.view()); +/// ``` +#[derive(Debug, Clone, Copy)] +pub struct GraphGrowth { + pub part_count: usize, +} - initial_ids +impl<'a, W> crate::Partition<(CsMatView<'a, f64>, W)> for GraphGrowth +where + W: AsRef<[f64]>, +{ + type Metadata = (); + type Error = std::convert::Infallible; + + fn partition( + &mut self, + part_ids: &mut [usize], + (adjacency, weights): (CsMatView, W), + ) -> Result { + graph_growth( + part_ids, + weights.as_ref(), + adjacency.view(), + self.part_count, + ); + Ok(()) + } } diff --git a/src/algorithms/greedy.rs b/src/algorithms/greedy.rs index 4f778d55..3604f8ac 100644 --- a/src/algorithms/greedy.rs +++ b/src/algorithms/greedy.rs @@ -2,7 +2,7 @@ use num::Zero; use std::ops::AddAssign; /// Implementation of the greedy algorithm. -pub fn greedy( +fn greedy( partition: &mut [usize], weights: impl IntoIterator, num_parts: usize, @@ -31,3 +31,40 @@ pub fn greedy( part_weights[min_part_weight_idx] += weight; } } + +/// Greedily assign weights to each part. +/// +/// # Example +/// +/// ```rust +/// use coupe::Partitioner; +/// use coupe::Point2D; +/// use coupe::Greedy; +/// +/// let num_parts = 3; +/// let points: &[Point2D] = &[Point2D::new(1., 0.), Point2D::new(0., 1.)]; +/// +/// let g = Greedy::new(num_parts); +/// g.partition(points, &[1., 1.]); +/// ``` +pub struct Greedy { + pub part_count: usize, +} + +impl crate::Partition for Greedy +where + W: IntoIterator, + W::Item: Ord + Zero + Clone + AddAssign, +{ + type Metadata = (); + type Error = std::convert::Infallible; + + fn partition( + &mut self, + part_ids: &mut [usize], + weights: W, + ) -> Result { + greedy(part_ids, weights, self.part_count); + Ok(()) + } +} diff --git a/src/algorithms/hilbert_curve.rs b/src/algorithms/hilbert_curve.rs index ee41a7ef..da142354 100644 --- a/src/algorithms/hilbert_curve.rs +++ b/src/algorithms/hilbert_curve.rs @@ -13,12 +13,13 @@ use crate::geometry::{Mbr, Point2D}; use rayon::prelude::*; -pub fn hilbert_curve_partition( +fn hilbert_curve_partition( + partition: &mut [usize], points: &[Point2D], weights: &[f64], - num_partitions: usize, + part_count: usize, order: usize, -) -> Vec { +) { assert!( order < 32, "Cannot construct a Hilbert curve of order >= 32 because 2^32 would overflow u32 capacity." @@ -36,14 +37,12 @@ pub fn hilbert_curve_partition( .par_sort_by_key(|idx| hilbert_indices[*idx]); // dummy modifiers to use directly the routine from multi_jagged - let modifiers = vec![1. / num_partitions as f64; num_partitions]; - - let mut partition = vec![0; points.len()]; + let modifiers = vec![1. / part_count as f64; part_count]; let split_positions = super::multi_jagged::compute_split_positions( weights, &permutation, - num_partitions - 1, + part_count - 1, &modifiers, ); @@ -60,8 +59,6 @@ pub fn hilbert_curve_partition( unsafe { std::ptr::write(ptr.add(*i), part_id) } } }); - - partition } #[allow(unused)] @@ -209,6 +206,77 @@ fn segment_to_segment(a_min: f64, a_max: f64, b_min: f64, b_max: f64) -> impl Fn } } +/// # Hilbert space-filling curve algorithm +/// +/// An implementation of the Hilbert curve based on +/// "Encoding and Decoding the Hilbert Order" by XIAN LIU and GÜNTHER SCHRACK. +/// +/// This algorithm uses space hashing to reorder points alongside the Hilbert curve ov a giver order. +/// See [wikipedia](https://en.wikipedia.org/wiki/Hilbert_curve) for more details. +/// +/// # Example +/// +/// ```rust +/// use coupe::Point2D; +/// use coupe::Partitioner; +/// +/// let points = vec![ +/// Point2D::new(0., 0.), +/// Point2D::new(1., 1.), +/// Point2D::new(0., 10.), +/// Point2D::new(1., 9.), +/// Point2D::new(9., 1.), +/// Point2D::new(10., 0.), +/// Point2D::new(10., 10.), +/// Point2D::new(9., 9.), +/// ]; +/// +/// let weights = vec![1.; 8]; +/// +/// let part_count = 4; +/// let order = 5; +/// +/// // generate a partition of 4 parts +/// let hilbert = coupe::HilbertCurve { part_count, order }; +/// +/// let partition = hilbert.partition(points.as_slice(), &weights); +/// let ids = partition.ids(); +/// +/// assert_eq!(ids[0], ids[1]); +/// assert_eq!(ids[2], ids[3]); +/// assert_eq!(ids[4], ids[5]); +/// assert_eq!(ids[6], ids[7]); +/// ``` +pub struct HilbertCurve { + pub part_count: usize, + pub order: u32, +} + +// hilbert curve is only implemented in 2d for now +impl crate::Partition<(P, W)> for HilbertCurve +where + P: AsRef<[Point2D]>, + W: AsRef<[f64]>, +{ + type Metadata = (); + type Error = std::convert::Infallible; + + fn partition( + &mut self, + part_ids: &mut [usize], + (points, weights): (P, W), + ) -> Result { + hilbert_curve_partition( + part_ids, + points.as_ref(), + weights.as_ref(), + self.part_count, + self.order as usize, + ); + Ok(()) + } +} + #[cfg(test)] mod tests { use super::*; diff --git a/src/algorithms/k_means.rs b/src/algorithms/k_means.rs index 6e78814e..6a2a5e01 100644 --- a/src/algorithms/k_means.rs +++ b/src/algorithms/k_means.rs @@ -32,7 +32,7 @@ type ClusterId = usize; /// this version shows some noticeable oscillations when imposing a restrictive balance constraint. /// It also skips the bounding boxes optimization which would slightly reduce the complexity of the /// algorithm. -pub fn simplified_k_means( +fn simplified_k_means( points: &[PointND], weights: &[f64], num_partitions: usize, @@ -216,7 +216,7 @@ impl Default for BalancedKmeansSettings { } } -pub fn balanced_k_means( +fn balanced_k_means( points: &[PointND], weights: &[f64], settings: impl Into>, @@ -298,7 +298,7 @@ where assignments } -pub fn balanced_k_means_with_initial_partition( +fn balanced_k_means_with_initial_partition( points: &[PointND], weights: &[f64], settings: impl Into>, @@ -737,3 +737,120 @@ fn max_distance(points: &[PointND]) -> f64 { .max_by(|a, b| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal)) .unwrap() } + +/// K-means algorithm +/// +/// An implementation of the balanced k-means algorithm inspired from +/// "Balanced k-means for Parallel Geometric Partitioning" by Moritz von Looz, +/// Charilaos Tzovas and Henning Meyerhenke (2018, University of Cologne). +/// +/// From an initial partition, the K-means algorithm will generate points clusters that will, +/// at each iteration, exchage points with other clusters that are "closer", and move by recomputing the clusters position (defined as +/// the centroid of the points assigned to the cluster). Eventually the clusters will stop moving, yielding a new partition. +/// +/// # Example +/// +/// ```rust +/// use coupe::Point2D; +/// use coupe::PartitionImprover; +/// use coupe::partition::Partition; +/// +/// // create ids for initial partition +/// +/// let points = vec![ +/// Point2D::new(0., 0.), +/// Point2D::new(1., 0.), +/// Point2D::new(2., 0.), +/// Point2D::new(0., 5.), +/// Point2D::new(1., 5.), +/// Point2D::new(2., 5.), +/// Point2D::new(0., 10.), +/// Point2D::new(1., 10.), +/// Point2D::new(2., 10.), +/// ]; +/// +/// let weights = vec![1.; 9]; +/// +/// // create an unbalanced partition: +/// // - p1: total weight = 1 +/// // - p2: total weight = 7 +/// // - p3: total weight = 1 +/// let ids = vec![1, 2, 2, 2, 2, 2, 2, 2, 3]; +/// let partition = Partition::from_ids(&points, &weights, ids); +/// +/// let mut k_means = coupe::KMeans::default(); +/// k_means.num_partitions = 3; +/// k_means.delta_threshold = 0.; +/// +/// let partition = k_means.improve_partition(partition); +/// let ids = partition.ids(); +/// +/// assert_eq!(ids[0], ids[1]); +/// assert_eq!(ids[0], ids[2]); +/// +/// assert_eq!(ids[3], ids[4]); +/// assert_eq!(ids[3], ids[5]); +/// +/// assert_eq!(ids[6], ids[7]); +/// assert_eq!(ids[6], ids[8]); +/// ``` +#[derive(Debug, Clone, Copy)] +pub struct KMeans { + pub part_count: usize, + pub imbalance_tol: f64, + pub delta_threshold: f64, + pub max_iter: usize, + pub max_balance_iter: usize, + pub erode: bool, + pub hilbert: bool, + pub mbr_early_break: bool, +} + +impl Default for KMeans { + fn default() -> Self { + Self { + part_count: 7, + imbalance_tol: 5., + delta_threshold: 0.01, + max_iter: 500, + max_balance_iter: 20, // for now, `max_balance_iter > 1` yields poor convergence time + erode: false, // for now, `erode` yields` enabled yields wrong results + hilbert: true, + mbr_early_break: false, // for now, `mbr_early_break` enabled yields wrong results + } + } +} + +impl<'a, const D: usize> crate::Partition<(&'a [PointND], &'a [f64])> for KMeans +where + Const: DimSub>, + DefaultAllocator: Allocator, Const, Buffer = ArrayStorage> + + Allocator, Const<1>>>, +{ + type Metadata = (); + type Error = std::convert::Infallible; + + fn partition( + &mut self, + part_ids: &mut [usize], + (points, weights): (&'a [PointND], &'a [f64]), + ) -> Result { + let settings = BalancedKmeansSettings { + num_partitions: self.part_count, + imbalance_tol: self.imbalance_tol, + delta_threshold: self.delta_threshold, + max_iter: self.max_iter, + max_balance_iter: self.max_balance_iter, + erode: self.erode, + hilbert: self.hilbert, + mbr_early_break: self.mbr_early_break, + }; + balanced_k_means_with_initial_partition( + points.as_ref(), + weights.as_ref(), + settings, + part_ids, + ); + Ok(()) + } +} diff --git a/src/algorithms/kernighan_lin.rs b/src/algorithms/kernighan_lin.rs index 31654f6d..16a54646 100644 --- a/src/algorithms/kernighan_lin.rs +++ b/src/algorithms/kernighan_lin.rs @@ -3,14 +3,12 @@ //! At each iteration, two nodes of different partition will be swapped, decreasing the overall cutsize //! of the partition. The swap is performed in such a way that the added partition imbalanced is controlled. -use crate::geometry::PointND; -use crate::partition::Partition; - use itertools::Itertools; use sprs::CsMatView; -pub(crate) fn kernighan_lin<'a, const D: usize>( - initial_partition: &mut Partition<'a, PointND, f64>, +pub(crate) fn kernighan_lin<'a>( + part_ids: &mut [usize], + weights: &[f64], adjacency: CsMatView, max_passes: impl Into>, max_flips_per_pass: impl Into>, @@ -25,12 +23,11 @@ pub(crate) fn kernighan_lin<'a, const D: usize>( let max_passes = max_passes.into(); let max_flips_per_pass = max_flips_per_pass.into(); let max_imbalance_per_flip = max_imbalance_per_flip.into(); - let (_points, weights, ids) = initial_partition.as_raw_mut(); - kernighan_lin_2_impl::( + kernighan_lin_2_impl( + part_ids, weights, adjacency.view(), - ids, max_passes, max_flips_per_pass, max_imbalance_per_flip, @@ -38,10 +35,10 @@ pub(crate) fn kernighan_lin<'a, const D: usize>( ); } -fn kernighan_lin_2_impl( +fn kernighan_lin_2_impl( + initial_partition: &mut [usize], weights: &[f64], adjacency: CsMatView, - initial_partition: &mut [usize], max_passes: Option, max_flips_per_pass: Option, _max_imbalance_per_flip: Option, @@ -180,3 +177,109 @@ fn kernighan_lin_2_impl( } println!("final cut size: {}", new_cut_size); } + +/// KernighanLin algorithm +/// +/// An implementation of the Kernighan Lin topologic algorithm +/// for graph partitioning. The current implementation currently only handles +/// partitioning a graph into two parts, as described in the original algorithm in +/// "An efficient heuristic procedure for partitioning graphs" by W. Kernighan and S. Lin. +/// +/// The algorithms repeats an iterative pass during which several pairs of nodes have +/// their part assignment swapped in order to reduce the cutsize of the partition. +/// If all the nodes are equally weighted, the algorithm preserves the partition balance. +/// +/// # Example +/// +/// ```rust +/// use coupe::Point2D; +/// use coupe::TopologicPartitionImprover; +/// use coupe::partition::Partition; +/// use sprs::CsMat; +/// +/// // swap +/// // 0 1 0 1 +/// // +--+--+--+ +/// // | | | | +/// // +--+--+--+ +/// // 0 0 1 1 +/// let points = vec![ +/// Point2D::new(0., 0.), +/// Point2D::new(1., 0.), +/// Point2D::new(2., 0.), +/// Point2D::new(3., 0.), +/// Point2D::new(0., 1.), +/// Point2D::new(1., 1.), +/// Point2D::new(2., 1.), +/// Point2D::new(3., 1.), +/// ]; +/// +/// let ids = vec![0, 0, 1, 1, 0, 1, 0, 1]; +/// let weights = vec![1.; 8]; +/// +/// let mut partition = Partition::from_ids(&points, &weights, ids); +/// +/// let mut adjacency = CsMat::empty(sprs::CSR, 8); +/// adjacency.reserve_outer_dim(8); +/// eprintln!("shape: {:?}", adjacency.shape()); +/// adjacency.insert(0, 1, 1.); +/// adjacency.insert(1, 2, 1.); +/// adjacency.insert(2, 3, 1.); +/// adjacency.insert(4, 5, 1.); +/// adjacency.insert(5, 6, 1.); +/// adjacency.insert(6, 7, 1.); +/// adjacency.insert(0, 4, 1.); +/// adjacency.insert(1, 5, 1.); +/// adjacency.insert(2, 6, 1.); +/// adjacency.insert(3, 7, 1.); +/// +/// // symmetry +/// adjacency.insert(1, 0, 1.); +/// adjacency.insert(2, 1, 1.); +/// adjacency.insert(3, 2, 1.); +/// adjacency.insert(5, 4, 1.); +/// adjacency.insert(6, 5, 1.); +/// adjacency.insert(7, 6, 1.); +/// adjacency.insert(4, 0, 1.); +/// adjacency.insert(5, 1, 1.); +/// adjacency.insert(6, 2, 1.); +/// adjacency.insert(7, 3, 1.); +/// +/// // 1 iteration +/// let algo = coupe::KernighanLin::new(1, 1, None, 1); +/// +/// let partition = algo.improve_partition(partition, adjacency.view()); +/// +/// let new_ids = partition.into_ids(); +/// assert_eq!(new_ids[5], 0); +/// assert_eq!(new_ids[6], 1); +/// ``` +#[derive(Debug, Clone, Copy)] +pub struct KernighanLin { + pub max_passes: Option, + pub max_flips_per_pass: Option, + pub max_imbalance_per_flip: Option, + pub max_bad_move_in_a_row: usize, +} + +impl<'a> crate::Partition<(CsMatView<'a, f64>, &'a [f64])> for KernighanLin { + type Metadata = (); + type Error = std::convert::Infallible; + + fn partition( + &mut self, + part_ids: &mut [usize], + (adjacency, weights): (CsMatView, &'a [f64]), + ) -> Result { + kernighan_lin( + part_ids, + weights, + adjacency, + self.max_passes, + self.max_flips_per_pass, + self.max_imbalance_per_flip, + self.max_bad_move_in_a_row, + ); + Ok(()) + } +} diff --git a/src/algorithms/kk.rs b/src/algorithms/kk.rs index ce091aa3..7b74244a 100644 --- a/src/algorithms/kk.rs +++ b/src/algorithms/kk.rs @@ -51,7 +51,7 @@ where } /// Implementation of the Karmarkar-Karp algorithm -pub fn kk(partition: &mut [usize], weights: I, num_parts: usize) +fn kk(partition: &mut [usize], weights: I, num_parts: usize) where T: Zero + Ord + Sub + SubAssign + Copy, I: IntoIterator, @@ -134,3 +134,41 @@ where parts.truncate(partition.len()); partition.copy_from_slice(&parts); } + +/// The Karmarkar-Karp algorithm, also called the Largest Differencing Method. +/// +/// # Example +/// +/// ```rust +/// use coupe::Partitioner; +/// use coupe::Point2D; +/// use coupe::KarmarkarKarp; +/// +/// let num_parts = 3; +/// let points: &[Point2D] = &[Point2D::new(1., 0.), Point2D::new(0., 1.)]; +/// +/// let kk = KarmarkarKarp::new(num_parts); +/// kk.partition(points, &[1., 1.]); +/// ``` +pub struct KarmarkarKarp { + pub part_count: usize, +} + +impl crate::Partition for KarmarkarKarp +where + W: IntoIterator, + W::IntoIter: ExactSizeIterator, + W::Item: Zero + Ord + Sub + SubAssign + Copy, +{ + type Metadata = (); + type Error = std::convert::Infallible; + + fn partition( + &mut self, + part_ids: &mut [usize], + weights: W, + ) -> Result { + kk(part_ids, weights, self.part_count); + Ok(()) + } +} diff --git a/src/algorithms/multi_jagged.rs b/src/algorithms/multi_jagged.rs index 26018af3..746541d4 100644 --- a/src/algorithms/multi_jagged.rs +++ b/src/algorithms/multi_jagged.rs @@ -145,35 +145,34 @@ fn compute_modifiers( .collect() } -pub fn multi_jagged( +fn multi_jagged( + partition: &mut [usize], points: &[PointND], weights: &[f64], num_parts: usize, max_iter: usize, -) -> Vec { +) { let partition_scheme = partition_scheme(num_parts, max_iter); - multi_jagged_with_scheme(points, weights, partition_scheme) + multi_jagged_with_scheme(partition, points, weights, partition_scheme); } fn multi_jagged_with_scheme( + partition: &mut [usize], points: &[PointND], weights: &[f64], partition_scheme: PartitionScheme, -) -> Vec { +) { let len = points.len(); let mut permutation = (0..len).into_par_iter().collect::>(); - let mut initial_partition = vec![0; points.len()]; multi_jagged_recurse( points, weights, &mut permutation, - &AtomicPtr::new(initial_partition.as_mut_ptr()), + &AtomicPtr::new(partition.as_mut_ptr()), 0, partition_scheme, ); - - initial_partition } fn multi_jagged_recurse( @@ -322,6 +321,84 @@ pub(crate) fn split_at_mut_many<'a, T>( head } +/// # Multi-Jagged algorithm +/// +/// This algorithm is inspired by Multi-Jagged: A Scalable Parallel Spatial Partitioning Algorithm" +/// by Mehmet Deveci, Sivasankaran Rajamanickam, Karen D. Devine, Umit V. Catalyurek. +/// +/// It improves over [RCB](struct.Rcb.html) by following the same idea but by creating more than two subparts +/// in each iteration which leads to decreasing recursion depth. It also allows to generate a partition +/// of any number of parts. +/// +/// More precisely, given a number of parts, the algorithm will generate a "partition scheme", which describes how +/// to perform splits at each iteration, such that the total number of iteration is less than `max_iter`. +/// +/// More iteration does not necessarily result in a better partition. +/// +/// # Example +/// +/// ```rust +/// use coupe::Point2D; +/// use coupe::Partitioner; +/// +/// let points = vec![ +/// Point2D::new(0., 0.), +/// Point2D::new(1., 0.), +/// Point2D::new(2., 0.), +/// Point2D::new(0., 1.), +/// Point2D::new(1., 1.), +/// Point2D::new(2., 1.), +/// Point2D::new(0., 2.), +/// Point2D::new(1., 2.), +/// Point2D::new(2., 2.), +/// ]; +/// +/// let weights = vec![1.; 9]; +/// +/// let num_partitions = 9; +/// let max_iter = 4; +/// +/// // generate a partition of 4 parts +/// let multi_jagged = coupe::MultiJagged::new(num_partitions, max_iter); +/// +/// let partition = multi_jagged.partition(points.as_slice(), &weights); +/// +/// let ids = partition.ids(); +/// +/// for i in 0..9 { +/// for j in 0..9 { +/// if j == i { +/// continue +/// } +/// assert_ne!(ids[i], ids[j]) +/// } +/// } +/// ``` +pub struct MultiJagged { + pub num_partitions: usize, + pub max_iter: usize, +} + +impl<'a, const D: usize> crate::Partition<(&'a [PointND], &'a [f64])> for MultiJagged { + type Metadata = (); + type Error = std::convert::Infallible; + + fn partition( + &mut self, + part_ids: &mut [usize], + (points, weights): (&'a [PointND], &'a [f64]), + ) -> Result { + multi_jagged( + part_ids, + points.as_ref(), + weights.as_ref(), + self.num_partitions, + self.max_iter, + ); + Ok(()) + } +} + #[cfg(test)] mod tests { use super::*; diff --git a/src/algorithms/recursive_bisection.rs b/src/algorithms/recursive_bisection.rs index f0671614..4ad1750c 100644 --- a/src/algorithms/recursive_bisection.rs +++ b/src/algorithms/recursive_bisection.rs @@ -9,6 +9,7 @@ use nalgebra::Const; use nalgebra::DefaultAllocator; use nalgebra::DimDiff; use nalgebra::DimSub; +use nalgebra::ToTypenum; use rayon::prelude::*; use std::cmp; use std::future::Future; @@ -500,7 +501,7 @@ fn rcb_thread( futures_lite::future::block_on(task); } -pub fn rcb(partition: &mut [usize], points: P, weights: W, iter_count: usize) +fn rcb(partition: &mut [usize], points: P, weights: W, iter_count: usize) where P: rayon::iter::IntoParallelIterator>, P::Iter: rayon::iter::IndexedParallelIterator, @@ -546,6 +547,72 @@ where }); } +/// # Recursive Coordinate Bisection algorithm +/// +/// Partitions a mesh based on the nodes coordinates and coresponding weights. +/// +/// This is the most simple and straightforward geometric algorithm. It operates as follows for a N-dimensional set of points: +/// +/// At each iteration, select a vector `n` of the canonical basis `(e_0, ..., e_{n-1})`. Then, split the set of points with an hyperplane orthogonal +/// to `n`, such that the two parts of the splits are evenly weighted. Finally, recurse by reapplying the algorithm to the two parts with an other +/// normal vector selection. +/// +/// # Example +/// +/// ```rust +/// use coupe::Point2D; +/// use coupe::Partitioner; +/// +/// let points = vec![ +/// Point2D::new(1., 1.), +/// Point2D::new(-1., 1.), +/// Point2D::new(1., -1.), +/// Point2D::new(-1., -1.), +/// ]; +/// +/// let weights = vec![1., 1., 1., 1.]; +/// +/// // generate a partition of 4 parts +/// let rcb = coupe::Rcb::new(2); +/// let partition = rcb.partition(points.as_slice(), &weights); +/// +/// let ids = partition.ids(); +/// for i in 0..4 { +/// for j in 0..4 { +/// if j == i { +/// continue +/// } +/// assert_ne!(ids[i], ids[j]) +/// } +/// } +/// ``` +pub struct Rcb { + pub iter_count: usize, +} + +impl crate::Partition<(P, W)> for Rcb +where + P: rayon::iter::IntoParallelIterator>, + P::Iter: rayon::iter::IndexedParallelIterator, + W: rayon::iter::IntoParallelIterator, + W::Item: Copy + std::fmt::Debug + Default, + W::Item: Add + AddAssign + Sub + Sum + PartialOrd, + W::Item: num::ToPrimitive, + W::Iter: rayon::iter::IndexedParallelIterator, +{ + type Metadata = (); + type Error = std::convert::Infallible; + + fn partition( + &mut self, + part_ids: &mut [usize], + (points, weights): (P, W), + ) -> Result { + rcb(part_ids, points, weights, self.iter_count); + Ok(()) + } +} + // pub because it is also useful for multijagged and required for benchmarks pub fn axis_sort( points: &[PointND], @@ -579,12 +646,8 @@ pub fn axis_sort( /// The global shape of the data is first considered and the separator is computed to /// be parallel to the inertia axis of the global shape, which aims to lead to better shaped /// partitions. -pub fn rib( - partition: &mut [usize], - points: &[PointND], - weights: W, - n_iter: usize, -) where +fn rib(partition: &mut [usize], points: &[PointND], weights: W, n_iter: usize) +where Const: DimSub>, DefaultAllocator: Allocator, Const, Buffer = ArrayStorage> + Allocator, Const<1>>>, @@ -607,6 +670,75 @@ pub fn rib( rcb(partition, points, weights, n_iter) } +/// # Recursive Inertial Bisection algorithm +/// +/// Partitions a mesh based on the nodes coordinates and coresponding weights +/// +/// This is a variant of the [`Rcb`](struct.Rcb.html) algorithm, where a basis change is performed beforehand so that +/// the first coordinate of the new basis is colinear to the inertia axis of the set of points. This has the goal +/// of producing better shaped partition than [Rcb](struct.Rcb.html). +/// +/// # Example +/// +/// ```rust +/// use coupe::Point2D; +/// use coupe::Partitioner; +/// +/// // Here, the inertia axis is the y axis. +/// // We thus expect Rib to split horizontally first. +/// let points = vec![ +/// Point2D::new(1., 10.), +/// Point2D::new(-1., 10.), +/// Point2D::new(1., -10.), +/// Point2D::new(-1., -10.), +/// ]; +/// +/// let weights = vec![1., 1., 1., 1.]; +/// +/// // generate a partition of 2 parts (1 split) +/// let rib = coupe::Rib::new(1); +/// let partition = rib.partition(points.as_slice(), &weights); +/// +/// let ids = partition.ids(); +/// +/// // the two points at the top are in the same partition +/// assert_eq!(ids[0], ids[1]); +/// +/// // the two points at the bottom are in the same partition +/// assert_eq!(ids[2], ids[3]); +/// +/// // there are two different partition +/// assert_ne!(ids[1], ids[2]); +/// ``` +pub struct Rib { + /// The number of iterations of the algorithm. This will yield a partition of `2^num_iter` parts. + pub iter_count: usize, +} + +impl<'a, const D: usize, W> crate::Partition<(&'a [PointND], W)> for Rib +where + Const: DimSub> + ToTypenum, + DefaultAllocator: Allocator, Const, Buffer = ArrayStorage> + + Allocator, Const<1>>>, + W: rayon::iter::IntoParallelIterator, + W::Item: Copy + std::fmt::Debug + Default, + W::Item: Add + AddAssign + Sub + Sum + PartialOrd, + W::Item: num::ToPrimitive, + W::Iter: rayon::iter::IndexedParallelIterator, +{ + type Metadata = (); + type Error = std::convert::Infallible; + + fn partition( + &mut self, + part_ids: &mut [usize], + (points, weights): (&'a [PointND], W), + ) -> Result { + rib(part_ids, points, weights, self.iter_count); + Ok(()) + } +} + #[cfg(test)] mod tests { use super::*; diff --git a/src/algorithms/vn/best.rs b/src/algorithms/vn/best.rs index 490960e5..be45d80d 100644 --- a/src/algorithms/vn/best.rs +++ b/src/algorithms/vn/best.rs @@ -1,5 +1,4 @@ use crate::imbalance::compute_parts_load; -use crate::RunInfo; use itertools::Itertools; use num::One; use num::Zero; @@ -8,8 +7,9 @@ use std::ops::Div; use std::ops::Mul; use std::ops::Sub; -pub fn vn_best_mono(partition: &mut [usize], criterion: &[T], nb_parts: usize) -> RunInfo +fn vn_best_mono(partition: &mut [usize], criterion: W, nb_parts: usize) -> usize where + W: IntoIterator, T: AddAssign + Sub + Div + Mul, for<'a> &'a T: Sub, T: Zero + One, @@ -21,22 +21,31 @@ where two }; + let mut criterion: Vec<(T, usize)> = criterion.into_iter().zip(0..).collect(); + assert_eq!(partition.len(), criterion.len()); // We expect weights to be non-negative values - assert!(criterion.iter().all(|weight| *weight >= T::zero())); + assert!(criterion + .iter() + .all(|(weight, _weight_id)| *weight >= T::zero())); // We check if all weights are 0 because this screw the gain table // initialization and a quick fix could interfere with the algorithm. if partition.is_empty() || criterion.is_empty() - || criterion.iter().all(|item| item.is_zero()) + || criterion + .iter() + .all(|(weight, _weight_id)| weight.is_zero()) || nb_parts < 2 { - return RunInfo::skip(); + return 0; } - let mut part_loads = compute_parts_load(partition, nb_parts, criterion); - let mut criterion: Vec<(T, usize)> = criterion.iter().copied().zip(0..).collect(); + let mut part_loads = compute_parts_load( + partition, + nb_parts, + criterion.iter().map(|(weight, _weight_id)| *weight), + ); criterion.sort_unstable(); let mut algo_iterations = 0; @@ -109,8 +118,54 @@ where algo_iterations += 1; } - RunInfo { - algo_iterations: Some(algo_iterations), + algo_iterations +} + +/// Greedy Vector-of-Numbers algorithms. +/// +/// This algorithm greedily moves weights from parts to parts in such a way that +/// the imbalance gain is maximized on each step. +/// +/// # Example +/// +/// ```rust +/// use coupe::Partitioner; +/// use coupe::PartitionImprover; +/// use coupe::Point2D; +/// use coupe::Random; +/// use coupe::VnBest; +/// use rand; +/// +/// let num_parts = 3; +/// let points: &[Point2D] = &[Point2D::new(1., 0.), Point2D::new(0., 1.)]; +/// +/// let r = Random::new(rand::thread_rng(), num_parts); +/// let partition = r.partition(points, &[1., 1.]); +/// let vn = VnBest::new(num_parts); +/// vn.improve_partition(partition); +/// ``` +pub struct VnBest { + pub part_count: usize, +} + +impl crate::Partition for VnBest +where + W: IntoIterator, + W::Item: AddAssign + Sub + Div + Mul, + for<'a> &'a W::Item: Sub, + W::Item: Zero + One, + W::Item: Ord + Copy, +{ + type Metadata = usize; + type Error = std::convert::Infallible; + + fn partition( + &mut self, + part_ids: &mut [usize], + weights: W, + ) -> Result { + let algo_iterations = vn_best_mono::(part_ids, weights, self.part_count); + Ok(algo_iterations) } } @@ -126,7 +181,7 @@ mod tests { let w = [1, 2, 3, 4, 5, 6]; let mut part = vec![0 as usize; w.len()]; let imb_ini = imbalance::imbalance(2, &part, &w); - vn_best_mono::(&mut part, &w, 2); + vn_best_mono::<_, u32>(&mut part, w, 2); let imb_end = imbalance::imbalance(2, &part, &w); assert!(imb_end < imb_ini); println!("imbalance : {} < {}", imb_end, imb_ini); @@ -145,7 +200,7 @@ mod tests { }) ) { let imb_ini = imbalance::max_imbalance(2, &partition, &weights); - vn_best_mono::(&mut partition, &weights, 2); + vn_best_mono::<_, u64>(&mut partition, weights.iter().cloned(), 2); let imb_end = imbalance::max_imbalance(2, &partition, &weights); // Not sure if it is true for max_imbalance (i.e. weighter - lighter) proptest::prop_assert!(imb_end <= imb_ini, "{} < {}", imb_ini, imb_end); diff --git a/src/algorithms/vn/mod.rs b/src/algorithms/vn/mod.rs index 3e42fd2c..22f4211c 100644 --- a/src/algorithms/vn/mod.rs +++ b/src/algorithms/vn/mod.rs @@ -1 +1,3 @@ pub mod best; + +pub use best::VnBest; diff --git a/src/algorithms/z_curve.rs b/src/algorithms/z_curve.rs index c8b3ebbb..2305b48e 100644 --- a/src/algorithms/z_curve.rs +++ b/src/algorithms/z_curve.rs @@ -25,6 +25,7 @@ use nalgebra::Const; use nalgebra::DefaultAllocator; use nalgebra::DimDiff; use nalgebra::DimSub; +use nalgebra::ToTypenum; use rayon::prelude::*; use std::cmp::Ordering; @@ -36,13 +37,13 @@ use std::sync::atomic::{self, AtomicPtr}; type HashType = u128; const HASH_TYPE_MAX: HashType = std::u128::MAX; -pub fn z_curve_partition( +fn z_curve_partition( + partition: &mut [usize], points: &[PointND], - num_partitions: usize, + part_count: usize, order: u32, -) -> Vec -where - Const: DimSub>, +) where + Const: DimSub> + ToTypenum, DefaultAllocator: Allocator, Const, Buffer = ArrayStorage> + Allocator, Const<1>>>, { @@ -58,15 +59,13 @@ where let mut permutation: Vec<_> = (0..points.len()).into_par_iter().collect(); - let mut ininial_partition = vec![0; points.len()]; - // reorder points z_curve_partition_recurse(points, order, &mbr, &mut permutation); - let points_per_partition = points.len() / num_partitions; - let remainder = points.len() % num_partitions; + let points_per_partition = points.len() / part_count; + let remainder = points.len() % part_count; - let atomic_handle = AtomicPtr::from(ininial_partition.as_mut_ptr()); + let atomic_handle = AtomicPtr::from(partition.as_mut_ptr()); // give an id to each partition // @@ -85,8 +84,6 @@ where unsafe { std::ptr::write(ptr.add(*idx), id) } } }); - - ininial_partition } // reorders `permu` to sort points by increasing z-curve hash @@ -189,6 +186,68 @@ fn compute_hash(point: &PointND, order: u32, mbr: &Mbr) -> } } +/// # Z space-filling curve algorithm +/// +/// The Z-curve uses space hashing to partition points. The points in the same part of a partition +/// have the same Z-hash. This hash is computed by recursively constructing a N-dimensional region tree. +/// +/// # Example +/// +/// ```rust +/// use coupe::Point2D; +/// use coupe::Partitioner; +/// +/// let points = vec![ +/// Point2D::new(0., 0.), +/// Point2D::new(1., 1.), +/// Point2D::new(0., 10.), +/// Point2D::new(1., 9.), +/// Point2D::new(9., 1.), +/// Point2D::new(10., 0.), +/// Point2D::new(10., 10.), +/// Point2D::new(9., 9.), +/// ]; +/// +/// let weights = vec![1.; 8]; +/// +/// let part_count = 4; +/// let order = 5; +/// +/// // generate a partition of 4 parts +/// let z_curve = coupe::ZCurve::new(part_count, order); +/// +/// let partition = z_curve.partition(points.as_slice(), &weights); +/// let ids = partition.ids(); +/// +/// assert_eq!(ids[0], ids[1]); +/// assert_eq!(ids[2], ids[3]); +/// assert_eq!(ids[4], ids[5]); +/// assert_eq!(ids[6], ids[7]); +/// ``` +pub struct ZCurve { + pub part_count: usize, + pub order: u32, +} + +impl<'a, const D: usize> crate::Partition<&'a [PointND]> for ZCurve +where + Const: DimSub> + ToTypenum, + DefaultAllocator: Allocator, Const, Buffer = ArrayStorage> + + Allocator, Const<1>>>, +{ + type Metadata = (); + type Error = std::convert::Infallible; + + fn partition( + &mut self, + part_ids: &mut [usize], + points: &'a [PointND], + ) -> Result { + z_curve_partition(part_ids, points, self.part_count, self.order); + Ok(()) + } +} + #[cfg(test)] mod tests { use super::*; @@ -196,7 +255,7 @@ mod tests { #[test] fn test_partition() { - let points = vec![ + let points = [ Point2D::from([0., 0.]), Point2D::from([20., 10.]), Point2D::from([0., 10.]), @@ -207,8 +266,9 @@ mod tests { Point2D::from([4., 2.]), ]; - let ids = z_curve_partition(&points, 4, 1); - for id in ids.iter() { + let mut ids = [0; 8]; + z_curve_partition(&mut ids, &points, 4, 1); + for id in ids { println!("{}", id); } assert_eq!(ids[0], ids[7]); diff --git a/src/imbalance.rs b/src/imbalance.rs index 1fad5f4f..48571ec4 100644 --- a/src/imbalance.rs +++ b/src/imbalance.rs @@ -12,11 +12,11 @@ use num::Zero; pub fn compute_parts_load( partition: &[usize], num_parts: usize, - weights: &[T], + weights: impl IntoIterator, ) -> Vec { debug_assert!(*partition.iter().max().unwrap_or(&0) < num_parts); partition - .iter() + .into_iter() .zip(weights) .fold(vec![T::zero(); num_parts], |mut acc, (&part, w)| { acc[part] += w.clone(); @@ -64,7 +64,7 @@ where pub fn imbalance_target + PartialOrd + Copy>( targets: &[T], partition: &[usize], - weights: &[T], + weights: impl IntoIterator, ) -> T { let num_parts = targets.len(); debug_assert!(*partition.iter().max().unwrap_or(&0) < num_parts); @@ -81,7 +81,7 @@ pub fn imbalance_target + Pa pub fn max_imbalance>( num_parts: usize, partition: &[usize], - weights: &[T], + weights: impl IntoIterator, ) -> T { compute_parts_load(partition, num_parts, weights) .iter() diff --git a/src/lib.rs b/src/lib.rs index dc39216b..e20390be 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -2,10 +2,9 @@ //! //! # Crate Layout //! -//! Coupe exposes each algorithm with a struct that implements a trait. There are currently two traits available: -//! -//! - A [`Partitioner`] represents an algorithm that will generate a partition given a set of geometric points and weights. -//! - A [`PartitionImprover`] represents an algorithm that will improve an existing partition (previously generated with a [`Partitioner`]). +//! Coupe exposes a [`Partition`] trait, which is in turn implemented by +//! algorithms. See its documentation for more details. The trait is generic around its input, which means algorithms +//! can partition different type of collections (e.g. 2D and 3D meshes). //! //! # Available algorithms //! @@ -26,8 +25,10 @@ //! - [K-means][KMeans] //! - Number partitioning: //! + [VN-Best][VnBest] +//! - [Fiduccia-Mattheyses][FiducciaMattheyses] +//! - [Kernighan-Lin][KernighanLin] -pub mod algorithms; +mod algorithms; pub mod geometry; pub mod imbalance; pub mod partition; @@ -36,1314 +37,30 @@ mod run_info; pub mod topology; mod uid; -// API - -// SUBMODULES REEXPORT +pub use crate::algorithms::*; pub use crate::geometry::{Point2D, Point3D, PointND}; pub use crate::real::Real; pub use crate::run_info::RunInfo; pub use crate::uid::uid; -pub mod dimension { - pub use nalgebra::base::dimension::*; -} - -use nalgebra::allocator::Allocator; -use nalgebra::ArrayStorage; -use nalgebra::Const; -use nalgebra::DefaultAllocator; -use nalgebra::DimDiff; -use nalgebra::DimSub; -use ndarray::ArrayView1; -use ndarray::ArrayView2; -use rayon::iter::IntoParallelIterator as _; -use rayon::iter::ParallelIterator as _; -use sprs::CsMatView; - -use crate::partition::*; - -// Trait that allows conversions from/to different kinds of -// points views representation as partitioner inputs -// e.g. &[f64], &[PointND], slice from ndarray, ... -pub trait PointsView<'a, const D: usize> { - fn to_points_nd(self) -> &'a [PointND]; -} - -impl<'a, const D: usize> PointsView<'a, D> for &'a [f64] { - fn to_points_nd(self) -> &'a [PointND] { - if self.len() % D != 0 { - panic!("error: tried to convert a &[f64] to a &[PointND] with D = {}, but input slice has len {}", D, self.len()); - } - unsafe { std::slice::from_raw_parts(self.as_ptr() as *const PointND, self.len() / D) } - } -} - -impl<'a, const D: usize> PointsView<'a, D> for ArrayView1<'a, f64> { - fn to_points_nd(self) -> &'a [PointND] { - let slice = self.to_slice().expect( - "Cannot convert an ArrayView1 with dicontiguous storage repr to a slice. Try cloning the data into a contiguous array first" - ); - slice.to_points_nd() - } -} - -impl<'a, const D: usize> PointsView<'a, D> for ArrayView2<'a, f64> { - fn to_points_nd(self) -> &'a [PointND] { - let slice = self.to_slice().expect( - "Cannot convert an ArrayView2 with dicontiguous storage repr to a slice. Try cloning the data into a contiguous array first" - ); - slice.to_points_nd() - } -} - -impl<'a, const D: usize> PointsView<'a, D> for &'a [PointND] { - fn to_points_nd(self) -> &'a [PointND] { - // nothing to do - self - } -} - -/// # A geometric partitioning algorithm. -/// -/// Algorithms that implement [`Partitioner`](trait.Partitioner.html) operate on a set of geometric points and associated weights and generate -/// a partition. A partition is described by an array of ids, each unique id represents a part of the partition. -/// -/// These algorithms generate a partition from scratch, as opposed to those which implement [`PartitionImprover`](trait.PartitionImprover.html), which work on -/// an existing partition to improve it. -/// -/// See the [implementors](trait.Partitioner.html#implementors) for more information about the currently available algorithms. -pub trait Partitioner { - fn partition<'a>( - &self, - points: impl PointsView<'a, D>, - weights: &'a [f64], - ) -> Partition<'a, PointND, f64>; -} - -/// A geometric algorithm to improve a partition. -/// -/// Algorithms that implement [`PartitionImprover`](trait.PartitionImprover.html) operate on a set of geometric -/// points and associated weights to modify and improve an existing partition (typically generated by a [`Partitioner`](trait.Partitioner.html)). -/// -/// See the [implementors](trait.PartitionImprover.html#implementors) for more information about the currently available algorithms. -pub trait PartitionImprover { - fn improve_partition<'a>( - &self, - partition: Partition<'a, PointND, f64>, - ) -> Partition<'a, PointND, f64>; -} - -pub trait TopologicPartitioner { - fn partition<'a>( - &self, - points: impl PointsView<'a, D>, - weights: &'a [f64], - adjacency: CsMatView, - ) -> Partition<'a, PointND, f64>; -} - -pub trait TopologicPartitionImprover { - fn improve_partition<'a>( - &self, - partition: Partition<'a, PointND, f64>, - adjacency: CsMatView, - ) -> Partition<'a, PointND, f64>; -} - -/// # Recursive Coordinate Bisection algorithm -/// -/// Partitions a mesh based on the nodes coordinates and coresponding weights. -/// -/// This is the most simple and straightforward geometric algorithm. It operates as follows for a N-dimensional set of points: -/// -/// At each iteration, select a vector `n` of the canonical basis `(e_0, ..., e_{n-1})`. Then, split the set of points with an hyperplane orthogonal -/// to `n`, such that the two parts of the splits are evenly weighted. Finally, recurse by reapplying the algorithm to the two parts with an other -/// normal vector selection. -/// -/// # Example -/// -/// ```rust -/// use coupe::Point2D; -/// use coupe::Partitioner; -/// -/// let points = vec![ -/// Point2D::new(1., 1.), -/// Point2D::new(-1., 1.), -/// Point2D::new(1., -1.), -/// Point2D::new(-1., -1.), -/// ]; -/// -/// let weights = vec![1., 1., 1., 1.]; -/// -/// // generate a partition of 4 parts -/// let rcb = coupe::Rcb::new(2); -/// let partition = rcb.partition(points.as_slice(), &weights); -/// -/// let ids = partition.ids(); -/// for i in 0..4 { -/// for j in 0..4 { -/// if j == i { -/// continue -/// } -/// assert_ne!(ids[i], ids[j]) -/// } -/// } -/// ``` -pub struct Rcb { - pub num_iter: usize, -} - -impl Rcb { - pub fn new(num_iter: usize) -> Self { - Self { num_iter } - } -} - -impl Partitioner for Rcb { - fn partition<'a>( - &self, - points: impl PointsView<'a, D>, - weights: &'a [f64], - ) -> Partition<'a, PointND, f64> { - let points = points.to_points_nd(); - let mut ids = vec![0; points.len()]; - let p = points.into_par_iter().cloned(); - let w = weights.into_par_iter().cloned(); - crate::algorithms::recursive_bisection::rcb(&mut ids, p, w, self.num_iter); - Partition::from_ids(points, weights, ids) - } -} - -/// # Recursive Inertial Bisection algorithm -/// -/// Partitions a mesh based on the nodes coordinates and coresponding weights -/// -/// This is a variant of the [`Rcb`](struct.Rcb.html) algorithm, where a basis change is performed beforehand so that -/// the first coordinate of the new basis is colinear to the inertia axis of the set of points. This has the goal -/// of producing better shaped partition than [Rcb](struct.Rcb.html). -/// -/// # Example -/// -/// ```rust -/// use coupe::Point2D; -/// use coupe::Partitioner; -/// -/// // Here, the inertia axis is the y axis. -/// // We thus expect Rib to split horizontally first. -/// let points = vec![ -/// Point2D::new(1., 10.), -/// Point2D::new(-1., 10.), -/// Point2D::new(1., -10.), -/// Point2D::new(-1., -10.), -/// ]; -/// -/// let weights = vec![1., 1., 1., 1.]; -/// -/// // generate a partition of 2 parts (1 split) -/// let rib = coupe::Rib::new(1); -/// let partition = rib.partition(points.as_slice(), &weights); -/// -/// let ids = partition.ids(); -/// -/// // the two points at the top are in the same partition -/// assert_eq!(ids[0], ids[1]); -/// -/// // the two points at the bottom are in the same partition -/// assert_eq!(ids[2], ids[3]); -/// -/// // there are two different partition -/// assert_ne!(ids[1], ids[2]); -/// ``` -pub struct Rib { - /// The number of iterations of the algorithm. This will yield a partition of `2^num_iter` parts. - pub num_iter: usize, -} - -impl Rib { - pub fn new(num_iter: usize) -> Self { - Self { num_iter } - } -} - -impl Partitioner for Rib -where - Const: DimSub>, - DefaultAllocator: Allocator, Const, Buffer = ArrayStorage> - + Allocator, Const<1>>>, -{ - fn partition<'a>( - &self, - points: impl PointsView<'a, D>, - weights: &'a [f64], - ) -> Partition<'a, PointND, f64> { - let points = points.to_points_nd(); - let mut ids = vec![0; points.len()]; - let w = weights.into_par_iter().cloned(); - crate::algorithms::recursive_bisection::rib(&mut ids, points, w, self.num_iter); - Partition::from_ids(points, weights, ids) - } -} - -/// # Multi-Jagged algorithm -/// -/// This algorithm is inspired by Multi-Jagged: A Scalable Parallel Spatial Partitioning Algorithm" -/// by Mehmet Deveci, Sivasankaran Rajamanickam, Karen D. Devine, Umit V. Catalyurek. -/// -/// It improves over [RCB](struct.Rcb.html) by following the same idea but by creating more than two subparts -/// in each iteration which leads to decreasing recursion depth. It also allows to generate a partition -/// of any number of parts. -/// -/// More precisely, given a number of parts, the algorithm will generate a "partition scheme", which describes how -/// to perform splits at each iteration, such that the total number of iteration is less than `max_iter`. -/// -/// More iteration does not necessarily result in a better partition. -/// -/// # Example -/// -/// ```rust -/// use coupe::Point2D; -/// use coupe::Partitioner; -/// -/// let points = vec![ -/// Point2D::new(0., 0.), -/// Point2D::new(1., 0.), -/// Point2D::new(2., 0.), -/// Point2D::new(0., 1.), -/// Point2D::new(1., 1.), -/// Point2D::new(2., 1.), -/// Point2D::new(0., 2.), -/// Point2D::new(1., 2.), -/// Point2D::new(2., 2.), -/// ]; -/// -/// let weights = vec![1.; 9]; -/// -/// let num_partitions = 9; -/// let max_iter = 4; -/// -/// // generate a partition of 4 parts -/// let multi_jagged = coupe::MultiJagged::new(num_partitions, max_iter); -/// -/// let partition = multi_jagged.partition(points.as_slice(), &weights); -/// -/// let ids = partition.ids(); -/// -/// for i in 0..9 { -/// for j in 0..9 { -/// if j == i { -/// continue -/// } -/// assert_ne!(ids[i], ids[j]) -/// } -/// } -/// ``` -pub struct MultiJagged { - pub num_partitions: usize, - pub max_iter: usize, -} - -impl MultiJagged { - pub fn new(num_partitions: usize, max_iter: usize) -> Self { - Self { - num_partitions, - max_iter, - } - } -} - -impl Partitioner for MultiJagged { - fn partition<'a>( - &self, - points: impl PointsView<'a, D>, - weights: &'a [f64], - ) -> Partition<'a, PointND, f64> { - let points = points.to_points_nd(); - let ids = crate::algorithms::multi_jagged::multi_jagged( - points, - weights, - self.num_partitions, - self.max_iter, - ); - - Partition::from_ids(points, weights, ids) - } -} - -/// # Z space-filling curve algorithm -/// -/// The Z-curve uses space hashing to partition points. The points in the same part of a partition -/// have the same Z-hash. This hash is computed by recursively constructing a N-dimensional region tree. -/// -/// # Example -/// -/// ```rust -/// use coupe::Point2D; -/// use coupe::Partitioner; -/// -/// let points = vec![ -/// Point2D::new(0., 0.), -/// Point2D::new(1., 1.), -/// Point2D::new(0., 10.), -/// Point2D::new(1., 9.), -/// Point2D::new(9., 1.), -/// Point2D::new(10., 0.), -/// Point2D::new(10., 10.), -/// Point2D::new(9., 9.), -/// ]; -/// -/// let weights = vec![1.; 8]; -/// -/// let num_partitions = 4; -/// let order = 5; -/// -/// // generate a partition of 4 parts -/// let z_curve = coupe::ZCurve::new(num_partitions, order); -/// -/// let partition = z_curve.partition(points.as_slice(), &weights); -/// let ids = partition.ids(); -/// -/// assert_eq!(ids[0], ids[1]); -/// assert_eq!(ids[2], ids[3]); -/// assert_eq!(ids[4], ids[5]); -/// assert_eq!(ids[6], ids[7]); -/// ``` -pub struct ZCurve { - pub num_partitions: usize, - pub order: u32, -} - -impl ZCurve { - pub fn new(num_partitions: usize, order: u32) -> Self { - Self { - num_partitions, - order, - } - } -} - -impl Partitioner for ZCurve -where - Const: DimSub>, - DefaultAllocator: Allocator, Const, Buffer = ArrayStorage> - + Allocator, Const<1>>>, -{ - fn partition<'a>( - &self, - points: impl PointsView<'a, D>, - weights: &'a [f64], - ) -> Partition<'a, PointND, f64> { - let points = points.to_points_nd(); - let ids = - crate::algorithms::z_curve::z_curve_partition(points, self.num_partitions, self.order); - Partition::from_ids(points, weights, ids) - } -} - -/// # Hilbert space-filling curve algorithm -/// -/// An implementation of the Hilbert curve based on -/// "Encoding and Decoding the Hilbert Order" by XIAN LIU and GÜNTHER SCHRACK. -/// -/// This algorithm uses space hashing to reorder points alongside the Hilbert curve ov a giver order. -/// See [wikipedia](https://en.wikipedia.org/wiki/Hilbert_curve) for more details. -/// -/// # Example -/// -/// ```rust -/// use coupe::Point2D; -/// use coupe::Partitioner; -/// -/// let points = vec![ -/// Point2D::new(0., 0.), -/// Point2D::new(1., 1.), -/// Point2D::new(0., 10.), -/// Point2D::new(1., 9.), -/// Point2D::new(9., 1.), -/// Point2D::new(10., 0.), -/// Point2D::new(10., 10.), -/// Point2D::new(9., 9.), -/// ]; -/// -/// let weights = vec![1.; 8]; -/// -/// let num_partitions = 4; -/// let order = 5; -/// -/// // generate a partition of 4 parts -/// let hilbert = coupe::HilbertCurve::new(num_partitions, order); -/// -/// let partition = hilbert.partition(points.as_slice(), &weights); -/// let ids = partition.ids(); -/// -/// assert_eq!(ids[0], ids[1]); -/// assert_eq!(ids[2], ids[3]); -/// assert_eq!(ids[4], ids[5]); -/// assert_eq!(ids[6], ids[7]); -/// ``` -pub struct HilbertCurve { - pub num_partitions: usize, - pub order: u32, -} - -impl HilbertCurve { - pub fn new(num_partitions: usize, order: u32) -> Self { - Self { - num_partitions, - order, - } - } -} - -// hilbert curve is only implemented in 2d for now -impl Partitioner<2> for HilbertCurve { - fn partition<'a>( - &self, - points: impl PointsView<'a, 2>, - _weights: &'a [f64], - ) -> Partition<'a, PointND<2>, f64> { - let points = points.to_points_nd(); - let ids = crate::algorithms::hilbert_curve::hilbert_curve_partition( - points, - _weights, - self.num_partitions, - self.order as usize, - ); - - Partition::from_ids(points, _weights, ids) - } -} - -use std::cell::RefCell; - -/// Map mesh elements to partitions randomly. -/// -/// # Example -/// -/// ```rust -/// use coupe::Partitioner; -/// use coupe::Point2D; -/// use coupe::Random; -/// use rand; -/// -/// let num_parts = 3; -/// let points: &[Point2D] = &[Point2D::new(1., 0.), Point2D::new(0., 1.)]; -/// -/// let r = Random::new(rand::thread_rng(), num_parts); -/// r.partition(points, &[1., 1.]); -/// ``` -pub struct Random { - rng: RefCell, - num_parts: usize, -} - -impl Random { - pub fn new(rng: R, num_parts: usize) -> Random { - Random { - rng: RefCell::new(rng), - num_parts, - } - } -} - -impl Partitioner for Random -where - R: rand::Rng, -{ - fn partition<'a>( - &self, - points: impl PointsView<'a, D>, - weights: &'a [f64], - ) -> Partition<'a, PointND, f64> { - let mut rng = self.rng.borrow_mut(); - let ids = (0..weights.len()) - .map(|_| rng.gen_range(0..self.num_parts)) - .collect(); - Partition::from_ids(points.to_points_nd(), weights, ids) - } -} - -/// Greedily assign weights to each part. -/// -/// # Example -/// -/// ```rust -/// use coupe::Partitioner; -/// use coupe::Point2D; -/// use coupe::Greedy; -/// -/// let num_parts = 3; -/// let points: &[Point2D] = &[Point2D::new(1., 0.), Point2D::new(0., 1.)]; -/// -/// let g = Greedy::new(num_parts); -/// g.partition(points, &[1., 1.]); -/// ``` -pub struct Greedy { - num_parts: usize, -} - -impl Greedy { - pub fn new(num_parts: usize) -> Greedy { - Greedy { num_parts } - } -} - -impl Partitioner for Greedy { - fn partition<'a>( - &self, - points: impl PointsView<'a, D>, - weights: &'a [f64], - ) -> Partition<'a, PointND, f64> { - let mut ids = vec![0; weights.len()]; - { - let weights = weights.iter().map(Real::from); - algorithms::greedy::greedy(&mut ids, weights, self.num_parts); - } - Partition::from_ids(points.to_points_nd(), weights, ids) - } -} - -/// The Karmarkar-Karp algorithm, also called the Largest Differencing Method. -/// -/// # Example -/// -/// ```rust -/// use coupe::Partitioner; -/// use coupe::Point2D; -/// use coupe::KarmarkarKarp; -/// -/// let num_parts = 3; -/// let points: &[Point2D] = &[Point2D::new(1., 0.), Point2D::new(0., 1.)]; -/// -/// let kk = KarmarkarKarp::new(num_parts); -/// kk.partition(points, &[1., 1.]); -/// ``` -pub struct KarmarkarKarp { - num_parts: usize, -} - -impl KarmarkarKarp { - pub fn new(num_parts: usize) -> KarmarkarKarp { - KarmarkarKarp { num_parts } - } -} - -impl Partitioner for KarmarkarKarp { - fn partition<'a>( - &self, - points: impl PointsView<'a, D>, - weights: &'a [f64], - ) -> Partition<'a, PointND, f64> { - let mut ids = vec![0; weights.len()]; - { - let weights = weights.iter().map(Real::from); - algorithms::kk::kk(&mut ids, weights, self.num_parts); - } - Partition::from_ids(points.to_points_nd(), weights, ids) - } -} - -/// The exact variant of the [Karmarkar-Karp][KarmarkarKarp] algorithm. -/// -/// # Example -/// -/// ```rust -/// use coupe::Partitioner; -/// use coupe::Point2D; -/// use coupe::CompleteKarmarkarKarp; -/// -/// let tolerance = 0.1; -/// let points: &[Point2D] = &[Point2D::new(1., 0.), Point2D::new(0., 1.)]; -/// -/// let ckk = CompleteKarmarkarKarp::new(tolerance); -/// ckk.partition(points, &[1., 1.]); -/// ``` -pub struct CompleteKarmarkarKarp { - tolerance: f64, -} - -impl CompleteKarmarkarKarp { - pub fn new(tolerance: f64) -> CompleteKarmarkarKarp { - CompleteKarmarkarKarp { tolerance } - } -} - -impl Partitioner for CompleteKarmarkarKarp { - fn partition<'a>( - &self, - points: impl PointsView<'a, D>, - weights: &'a [f64], - ) -> Partition<'a, PointND, f64> { - let mut ids = vec![0; weights.len()]; - { - let weights = weights.iter().map(Real::from); - algorithms::ckk::ckk_bipart(&mut ids, weights, self.tolerance); - } - Partition::from_ids(points.to_points_nd(), weights, ids) - } -} - -/// Greedy Vector-of-Numbers algorithms. -/// -/// This algorithm greedily moves weights from parts to parts in such a way that -/// the imbalance gain is maximized on each step. -/// -/// # Example -/// -/// ```rust -/// use coupe::Partitioner; -/// use coupe::PartitionImprover; -/// use coupe::Point2D; -/// use coupe::Random; -/// use coupe::VnBest; -/// use rand; -/// -/// let num_parts = 3; -/// let points: &[Point2D] = &[Point2D::new(1., 0.), Point2D::new(0., 1.)]; -/// -/// let r = Random::new(rand::thread_rng(), num_parts); -/// let partition = r.partition(points, &[1., 1.]); -/// let vn = VnBest::new(num_parts); -/// vn.improve_partition(partition); -/// ``` -pub struct VnBest { - num_parts: usize, -} - -impl VnBest { - pub fn new(num_parts: usize) -> VnBest { - VnBest { num_parts } - } -} - -impl PartitionImprover for VnBest { - fn improve_partition<'a>( - &self, - partition: Partition<'a, PointND, f64>, - ) -> Partition<'a, PointND, f64> { - let (points, weights, mut ids) = partition.into_raw(); - let weights_real: Vec = weights.iter().map(Real::from).collect(); - algorithms::vn::best::vn_best_mono::(&mut ids, &weights_real, self.num_parts); - Partition::from_ids(points, weights, ids) - } -} - -/// K-means algorithm -/// -/// An implementation of the balanced k-means algorithm inspired from -/// "Balanced k-means for Parallel Geometric Partitioning" by Moritz von Looz, -/// Charilaos Tzovas and Henning Meyerhenke (2018, University of Cologne). -/// -/// From an initial partition, the K-means algorithm will generate points clusters that will, -/// at each iteration, exchage points with other clusters that are "closer", and move by recomputing the clusters position (defined as -/// the centroid of the points assigned to the cluster). Eventually the clusters will stop moving, yielding a new partition. -/// -/// # Example -/// -/// ```rust -/// use coupe::Point2D; -/// use coupe::PartitionImprover; -/// use coupe::partition::Partition; -/// -/// // create ids for initial partition -/// -/// let points = vec![ -/// Point2D::new(0., 0.), -/// Point2D::new(1., 0.), -/// Point2D::new(2., 0.), -/// Point2D::new(0., 5.), -/// Point2D::new(1., 5.), -/// Point2D::new(2., 5.), -/// Point2D::new(0., 10.), -/// Point2D::new(1., 10.), -/// Point2D::new(2., 10.), -/// ]; -/// -/// let weights = vec![1.; 9]; -/// -/// // create an unbalanced partition: -/// // - p1: total weight = 1 -/// // - p2: total weight = 7 -/// // - p3: total weight = 1 -/// let ids = vec![1, 2, 2, 2, 2, 2, 2, 2, 3]; -/// let partition = Partition::from_ids(&points, &weights, ids); +/// The `Partition` trait allows for partitioning data. /// -/// let mut k_means = coupe::KMeans::default(); -/// k_means.num_partitions = 3; -/// k_means.delta_threshold = 0.; +/// Partitioning algorithms implement this trait. /// -/// let partition = k_means.improve_partition(partition); -/// let ids = partition.ids(); +/// The generic argument `M` defines the input of the algorithms (e.g. an +/// adjacency matrix or a 2D set of points). /// -/// assert_eq!(ids[0], ids[1]); -/// assert_eq!(ids[0], ids[2]); -/// -/// assert_eq!(ids[3], ids[4]); -/// assert_eq!(ids[3], ids[5]); -/// -/// assert_eq!(ids[6], ids[7]); -/// assert_eq!(ids[6], ids[8]); -/// ``` -#[derive(Debug, Clone, Copy)] -pub struct KMeans { - pub num_partitions: usize, - pub imbalance_tol: f64, - pub delta_threshold: f64, - pub max_iter: usize, - pub max_balance_iter: usize, - pub erode: bool, - pub hilbert: bool, - pub mbr_early_break: bool, -} - -// KMeans builder pattern -// to reduce construction boilerplate -// e.g. -// ```rust -// let k_means = KMeansBuilder::default() -// .imbalance_tol(5.) -// .max_balance_iter(12) -// .build(); -// ``` -#[derive(Debug, Default, Clone, Copy)] -pub struct KMeansBuilder { - inner: KMeans, -} - -impl KMeansBuilder { - pub fn build(self) -> KMeans { - self.inner - } - - pub fn num_partitions(mut self, num_partitions: usize) -> Self { - self.inner.num_partitions = num_partitions; - self - } - - pub fn imbalance_tol(mut self, imbalance_tol: f64) -> Self { - self.inner.imbalance_tol = imbalance_tol; - self - } +/// The input partition must be of the correct size and its contents may or may +/// not be used by the algorithms. +pub trait Partition { + /// Diagnostic data returned for a specific run of the algorithm. + type Metadata; - pub fn delta_threshold(mut self, delta_threshold: f64) -> Self { - self.inner.delta_threshold = delta_threshold; - self - } - - pub fn max_iter(mut self, max_iter: usize) -> Self { - self.inner.max_iter = max_iter; - self - } - - pub fn max_balance_iter(mut self, max_balance_iter: usize) -> Self { - self.inner.max_balance_iter = max_balance_iter; - self - } - - pub fn erode(mut self, erode: bool) -> Self { - self.inner.erode = erode; - self - } - - pub fn hilbert(mut self, hilbert: bool) -> Self { - self.inner.hilbert = hilbert; - self - } - - pub fn mbr_early_break(mut self, mbr_early_break: bool) -> Self { - self.inner.mbr_early_break = mbr_early_break; - self - } -} - -impl KMeans { - #[allow(clippy::too_many_arguments)] - pub fn new( - num_partitions: usize, - imbalance_tol: f64, - delta_threshold: f64, - max_iter: usize, - max_balance_iter: usize, - erode: bool, - hilbert: bool, - mbr_early_break: bool, - ) -> Self { - Self { - num_partitions, - imbalance_tol, - delta_threshold, - max_iter, - max_balance_iter, - erode, - hilbert, - mbr_early_break, - } - } -} - -impl Default for KMeans { - fn default() -> Self { - Self { - num_partitions: 7, - imbalance_tol: 5., - delta_threshold: 0.01, - max_iter: 500, - max_balance_iter: 20, // for now, `max_balance_iter > 1` yields poor convergence time - erode: false, // for now, `erode` yields` enabled yields wrong results - hilbert: true, - mbr_early_break: false, // for now, `mbr_early_break` enabled yields wrong results - } - } -} - -impl PartitionImprover for KMeans -where - Const: DimSub>, - DefaultAllocator: Allocator, Const, Buffer = ArrayStorage> - + Allocator, Const<1>>>, -{ - fn improve_partition<'a>( - &self, - partition: Partition<'a, PointND, f64>, - ) -> Partition<'a, PointND, f64> { - let settings = crate::algorithms::k_means::BalancedKmeansSettings { - num_partitions: self.num_partitions, - imbalance_tol: self.imbalance_tol, - delta_threshold: self.delta_threshold, - max_iter: self.max_iter, - max_balance_iter: self.max_balance_iter, - erode: self.erode, - hilbert: self.hilbert, - mbr_early_break: self.mbr_early_break, - }; - let (points, weights, mut ids) = partition.into_raw(); - crate::algorithms::k_means::balanced_k_means_with_initial_partition( - points, weights, settings, &mut ids, - ); - Partition::from_ids(points, weights, ids) - } -} - -/// KernighanLin algorithm -/// -/// An implementation of the Kernighan Lin topologic algorithm -/// for graph partitioning. The current implementation currently only handles -/// partitioning a graph into two parts, as described in the original algorithm in -/// "An efficient heuristic procedure for partitioning graphs" by W. Kernighan and S. Lin. -/// -/// The algorithms repeats an iterative pass during which several pairs of nodes have -/// their part assignment swapped in order to reduce the cutsize of the partition. -/// If all the nodes are equally weighted, the algorithm preserves the partition balance. -/// -/// # Example -/// -/// ```rust -/// use coupe::Point2D; -/// use coupe::TopologicPartitionImprover; -/// use coupe::partition::Partition; -/// use sprs::CsMat; -/// -/// // swap -/// // 0 1 0 1 -/// // +--+--+--+ -/// // | | | | -/// // +--+--+--+ -/// // 0 0 1 1 -/// let points = vec![ -/// Point2D::new(0., 0.), -/// Point2D::new(1., 0.), -/// Point2D::new(2., 0.), -/// Point2D::new(3., 0.), -/// Point2D::new(0., 1.), -/// Point2D::new(1., 1.), -/// Point2D::new(2., 1.), -/// Point2D::new(3., 1.), -/// ]; -/// -/// let ids = vec![0, 0, 1, 1, 0, 1, 0, 1]; -/// let weights = vec![1.; 8]; -/// -/// let mut partition = Partition::from_ids(&points, &weights, ids); -/// -/// let mut adjacency = CsMat::empty(sprs::CSR, 8); -/// adjacency.reserve_outer_dim(8); -/// eprintln!("shape: {:?}", adjacency.shape()); -/// adjacency.insert(0, 1, 1.); -/// adjacency.insert(1, 2, 1.); -/// adjacency.insert(2, 3, 1.); -/// adjacency.insert(4, 5, 1.); -/// adjacency.insert(5, 6, 1.); -/// adjacency.insert(6, 7, 1.); -/// adjacency.insert(0, 4, 1.); -/// adjacency.insert(1, 5, 1.); -/// adjacency.insert(2, 6, 1.); -/// adjacency.insert(3, 7, 1.); -/// -/// // symmetry -/// adjacency.insert(1, 0, 1.); -/// adjacency.insert(2, 1, 1.); -/// adjacency.insert(3, 2, 1.); -/// adjacency.insert(5, 4, 1.); -/// adjacency.insert(6, 5, 1.); -/// adjacency.insert(7, 6, 1.); -/// adjacency.insert(4, 0, 1.); -/// adjacency.insert(5, 1, 1.); -/// adjacency.insert(6, 2, 1.); -/// adjacency.insert(7, 3, 1.); -/// -/// // 1 iteration -/// let algo = coupe::KernighanLin::new(1, 1, None, 1); -/// -/// let partition = algo.improve_partition(partition, adjacency.view()); -/// -/// let new_ids = partition.into_ids(); -/// assert_eq!(new_ids[5], 0); -/// assert_eq!(new_ids[6], 1); -/// ``` -#[derive(Debug, Clone, Copy)] -pub struct KernighanLin { - max_passes: Option, - max_flips_per_pass: Option, - max_imbalance_per_flip: Option, - max_bad_move_in_a_row: usize, -} - -impl KernighanLin { - pub fn new( - max_passes: impl Into>, - max_flips_per_pass: impl Into>, - max_imbalance_per_flip: impl Into>, - max_bad_move_in_a_row: usize, - ) -> Self { - Self { - max_passes: max_passes.into(), - max_flips_per_pass: max_flips_per_pass.into(), - max_imbalance_per_flip: max_imbalance_per_flip.into(), - max_bad_move_in_a_row, - } - } -} - -impl TopologicPartitionImprover for KernighanLin { - fn improve_partition<'a>( - &self, - mut partition: Partition<'a, PointND, f64>, - adjacency: CsMatView, - ) -> Partition<'a, PointND, f64> { - crate::algorithms::kernighan_lin::kernighan_lin( - &mut partition, - adjacency, - self.max_passes, - self.max_flips_per_pass, - self.max_imbalance_per_flip, - self.max_bad_move_in_a_row, - ); - partition - } -} - -/// FiducciaMattheyses -/// -/// An implementation of the Fiduccia Mattheyses topologic algorithm -/// for graph partitioning. This implementation is an extension of the -/// original algorithm to handle partitioning into more than two parts. -/// -/// This algorithm repeats an iterative pass during which a set of graph nodes are assigned to -/// a new part, reducing the overall cutsize of the partition. As opposed to the -/// Kernighan-Lin algorithm, during each pass iteration, only one node is flipped at a time. -/// The algorithm thus does not preserve partition weights balance and may produce an unbalanced -/// partition. -/// -/// Original algorithm from "A Linear-Time Heuristic for Improving Network Partitions" -/// by C.M. Fiduccia and R.M. Mattheyses. -/// -/// # Example -/// -/// ```rust -/// use coupe::Point2D; -/// use coupe::TopologicPartitionImprover; -/// use coupe::partition::Partition; -/// use sprs::CsMat; -/// -/// // swap -/// // 0 1 0 1 -/// // +--+--+--+ -/// // | | | | -/// // +--+--+--+ -/// // 0 0 1 1 -/// let points = vec![ -/// Point2D::new(0., 0.), -/// Point2D::new(1., 0.), -/// Point2D::new(2., 0.), -/// Point2D::new(3., 0.), -/// Point2D::new(0., 1.), -/// Point2D::new(1., 1.), -/// Point2D::new(2., 1.), -/// Point2D::new(3., 1.), -/// ]; -/// -/// let ids = vec![0, 0, 1, 1, 0, 1, 0, 1]; -/// let weights = vec![1.; 8]; -/// -/// let mut partition = Partition::from_ids(&points, &weights, ids); -/// -/// let mut adjacency = CsMat::empty(sprs::CSR, 8); -/// adjacency.reserve_outer_dim(8); -/// eprintln!("shape: {:?}", adjacency.shape()); -/// adjacency.insert(0, 1, 1.); -/// adjacency.insert(1, 2, 1.); -/// adjacency.insert(2, 3, 1.); -/// adjacency.insert(4, 5, 1.); -/// adjacency.insert(5, 6, 1.); -/// adjacency.insert(6, 7, 1.); -/// adjacency.insert(0, 4, 1.); -/// adjacency.insert(1, 5, 1.); -/// adjacency.insert(2, 6, 1.); -/// adjacency.insert(3, 7, 1.); -/// -/// // symmetry -/// adjacency.insert(1, 0, 1.); -/// adjacency.insert(2, 1, 1.); -/// adjacency.insert(3, 2, 1.); -/// adjacency.insert(5, 4, 1.); -/// adjacency.insert(6, 5, 1.); -/// adjacency.insert(7, 6, 1.); -/// adjacency.insert(4, 0, 1.); -/// adjacency.insert(5, 1, 1.); -/// adjacency.insert(6, 2, 1.); -/// adjacency.insert(7, 3, 1.); -/// -/// // 1 iteration -/// let algo = coupe::FiducciaMattheyses::new(None, None, None, 1); -/// -/// let partition = algo.improve_partition(partition, adjacency.view()); -/// -/// let new_ids = partition.into_ids(); -/// assert_eq!(new_ids[5], 0); -/// assert_eq!(new_ids[6], 1); -/// ``` -#[derive(Debug, Clone, Copy)] -pub struct FiducciaMattheyses { - max_passes: Option, - max_flips_per_pass: Option, - max_imbalance_per_flip: Option, - max_bad_move_in_a_row: usize, -} - -impl FiducciaMattheyses { - pub fn new( - max_passes: impl Into>, - max_flips_per_pass: impl Into>, - max_imbalance_per_flip: impl Into>, - max_bad_move_in_a_row: usize, - ) -> Self { - Self { - max_passes: max_passes.into(), - max_flips_per_pass: max_flips_per_pass.into(), - max_imbalance_per_flip: max_imbalance_per_flip.into(), - max_bad_move_in_a_row, - } - } -} - -impl TopologicPartitionImprover for FiducciaMattheyses { - fn improve_partition<'a>( - &self, - mut partition: Partition<'a, PointND, f64>, - adjacency: CsMatView, - ) -> Partition<'a, PointND, f64> { - crate::algorithms::fiduccia_mattheyses::fiduccia_mattheyses( - &mut partition, - adjacency, - self.max_passes, - self.max_flips_per_pass, - self.max_imbalance_per_flip, - self.max_bad_move_in_a_row, - ); - partition - } -} - -/// Graph Growth algorithm -/// -/// A topologic algorithm that generates a partition from a topologic mesh. -/// Given a number k of parts, the algorithm selects k nodes randomly and assigns them to a different part. -/// Then, at each iteration, each part is expanded to neighbor nodes that are not yet assigned to a part -/// -/// # Example -/// -/// ```rust -/// use coupe::Point2D; -/// use coupe::TopologicPartitioner; -/// use coupe::partition::Partition; -/// use sprs::CsMat; -/// -/// // +--+--+--+ -/// // | | | | -/// // +--+--+--+ -/// -/// let points = vec![ -/// Point2D::new(0., 0.), -/// Point2D::new(1., 0.), -/// Point2D::new(2., 0.), -/// Point2D::new(3., 0.), -/// Point2D::new(0., 1.), -/// Point2D::new(1., 1.), -/// Point2D::new(2., 1.), -/// Point2D::new(3., 1.), -/// ]; -/// -/// let weights = vec![1.; 8]; -/// -/// let mut adjacency = CsMat::empty(sprs::CSR, 8); -/// adjacency.reserve_outer_dim(8); -/// eprintln!("shape: {:?}", adjacency.shape()); -/// adjacency.insert(0, 1, 1.); -/// adjacency.insert(1, 2, 1.); -/// adjacency.insert(2, 3, 1.); -/// adjacency.insert(4, 5, 1.); -/// adjacency.insert(5, 6, 1.); -/// adjacency.insert(6, 7, 1.); -/// adjacency.insert(0, 4, 1.); -/// adjacency.insert(1, 5, 1.); -/// adjacency.insert(2, 6, 1.); -/// adjacency.insert(3, 7, 1.); -/// -/// // symmetry -/// adjacency.insert(1, 0, 1.); -/// adjacency.insert(2, 1, 1.); -/// adjacency.insert(3, 2, 1.); -/// adjacency.insert(5, 4, 1.); -/// adjacency.insert(6, 5, 1.); -/// adjacency.insert(7, 6, 1.); -/// adjacency.insert(4, 0, 1.); -/// adjacency.insert(5, 1, 1.); -/// adjacency.insert(6, 2, 1.); -/// adjacency.insert(7, 3, 1.); -/// -/// let gg = coupe::GraphGrowth::new(2); -/// -/// let _partition = gg.partition(points.as_slice(), &weights, adjacency.view()); -/// ``` -#[derive(Debug, Clone, Copy)] -pub struct GraphGrowth { - num_partitions: usize, -} - -impl GraphGrowth { - pub fn new(num_partitions: usize) -> Self { - Self { num_partitions } - } -} - -impl TopologicPartitioner for GraphGrowth { - fn partition<'a>( - &self, - points: impl PointsView<'a, D>, - weights: &'a [f64], - adjacency: CsMatView, - ) -> Partition<'a, PointND, f64> { - let ids = crate::algorithms::graph_growth::graph_growth( - weights, - adjacency.view(), - self.num_partitions, - ); - Partition::from_ids(points.to_points_nd(), weights, ids) - } -} - -/// # Represents the composition algorithm. -/// -/// This structure is created by calling the [`compose`](trait.Compose.html#tymethod.compose) -/// of the [`Compose`](trait.Compose.html) trait. -pub struct Composition { - first: T, - second: U, -} - -impl Composition { - pub fn new(first: T, second: U) -> Self { - Self { first, second } - } -} - -impl Partitioner for Composition -where - T: Partitioner, - U: PartitionImprover, -{ - fn partition<'a>( - &self, - points: impl PointsView<'a, D>, - weights: &'a [f64], - ) -> Partition<'a, PointND, f64> { - let points = points.to_points_nd(); - let partition = self.first.partition(points, weights); - self.second.improve_partition(partition) - } -} - -impl PartitionImprover for Composition -where - T: PartitionImprover, - U: PartitionImprover, -{ - fn improve_partition<'a>( - &self, - partition: Partition<'a, PointND, f64>, - ) -> Partition<'a, PointND, f64> { - self.second - .improve_partition(self.first.improve_partition(partition)) - } -} - -impl TopologicPartitioner for Composition -where - T: Partitioner, - U: TopologicPartitionImprover, -{ - fn partition<'a>( - &self, - points: impl PointsView<'a, D>, - weights: &'a [f64], - adjacency: CsMatView, - ) -> Partition<'a, PointND, f64> { - let partition = self.first.partition(points, weights); - self.second.improve_partition(partition, adjacency) - } -} - -/// # Compose two algorithms. -/// -/// This trait enables algorithms to be composed together. Doing so allows for instance to improve the partition -/// generated by a [`Partitioner`](trait.Partitioner.html) with a [`PartitionImprover`](trait.PartitionImprover.html). -/// -/// The resulting composition implements either [`Partitioner`](trait.Partitioner.html) or [`PartitionImprover`](trait.PartitionImprover.html) -/// based on the input algorithms. -/// -/// # Example -/// ```rust -/// use coupe::{Compose, Partitioner, PartitionImprover}; -/// use coupe::dimension::U3; -/// -/// let num_partitions = 7; -/// let max_iter = 3; -/// -/// let mut k_means = coupe::KMeans::default(); -/// k_means.num_partitions = 7; -/// let multi_jagged_then_k_means = coupe::MultiJagged::new( -/// num_partitions, -/// max_iter -/// ).compose(k_means); -/// -/// ``` -pub trait Compose { - type Composed; - fn compose(self, other: T) -> Self::Composed; -} + /// Error details, should the algorithm fail to run. + type Error; -impl Compose for U { - type Composed = Composition; - fn compose(self, other: T) -> Self::Composed { - Composition::new(self, other) - } + /// Partition the given data and output the part ID of each element in + /// `part_ids`. + fn partition(&mut self, part_ids: &mut [usize], data: M) + -> Result; }