From 5aebf8d9e7a590f3073a8cddac22a37e7438e2b6 Mon Sep 17 00:00:00 2001 From: Hubert Hirtz Date: Thu, 17 Mar 2022 09:06:27 +0100 Subject: [PATCH 1/6] Rework trait family This simplifies the family of traits from before into one "Partition" trait that is generic over the input of the algorithm. Thus, algorithms can partition points, weights or graphs without implementing several traits. The partition function can now mutate algorithm state, which is useful for PRNG state inside coupe::Random. Algorithms can now return errors instead of panicking, which is more FFI-friendly. Algorithms can now return arbitrary information (e.g. coupe::VnBest returns its number of iterations), which removes the need for coupe::RunInfo. == Limitations == Some implementations cannot be as generic as we want (e.g. Kmeans) At first, I wanted k-means to implement `Partition` for any `P: AsRef<[PointND]>` and `W: AsRef<[f64]>` (so that it can also be used with a &Vec<_>), but rust doesn't accept those kind of "loose" constraints on const-generic parameters. I suppose this will do for now. Surprisingly, RCB, which has looser constraints, is fine. == Misc changes == All algorithm structs have been moved in their respective module, to declutter lib.rs. Constructors/Builder patterns have been removed since they didn't add much (and a fn new with 5 parameters is not that readable). It's simpler to have a Default impl and let the user do Algo { specific_param, ..Algo::default() }. --- Cargo.lock | 1 + benches/benchmark.rs | 72 +- benches/generator.rs | 18 - src/algorithms.rs | 99 +- src/algorithms/ckk.rs | 55 +- src/algorithms/fiduccia_mattheyses.rs | 118 ++- src/algorithms/graph_growth.rs | 85 +- src/algorithms/greedy.rs | 39 +- src/algorithms/hilbert_curve.rs | 82 +- src/algorithms/k_means.rs | 328 ++---- src/algorithms/kernighan_lin.rs | 137 ++- src/algorithms/kk.rs | 39 +- src/algorithms/multi_jagged.rs | 82 +- src/algorithms/recursive_bisection.rs | 145 ++- src/algorithms/vn/best.rs | 80 +- src/algorithms/vn/mod.rs | 2 + src/algorithms/z_curve.rs | 85 +- src/geometry.rs | 67 +- src/imbalance.rs | 8 +- src/lib.rs | 1330 +------------------------ src/partition.rs | 442 -------- src/run_info.rs | 16 - tools/Cargo.toml | 1 + tools/num-part/src/main.rs | 3 +- tools/src/bin/mesh-part/main.rs | 330 +++--- 25 files changed, 1277 insertions(+), 2387 deletions(-) delete mode 100644 src/partition.rs delete mode 100644 src/run_info.rs diff --git a/Cargo.lock b/Cargo.lock index 8f2d51a2..fb9b67b1 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -161,6 +161,7 @@ dependencies = [ "num", "rand", "rand_pcg", + "rayon", ] [[package]] diff --git a/benches/benchmark.rs b/benches/benchmark.rs index dc5d6367..e57c19da 100644 --- a/benches/benchmark.rs +++ b/benches/benchmark.rs @@ -1,8 +1,7 @@ mod generator; -use coupe::algorithms::k_means::simplified_k_means; -use coupe::algorithms::recursive_bisection::{axis_sort, rcb}; use coupe::geometry::Point2D; +use coupe::Partition as _; use criterion::{criterion_group, criterion_main}; use criterion::{Criterion, Throughput}; use rayon::prelude::*; @@ -10,21 +9,6 @@ use rayon::prelude::*; const SAMPLE_SIZE: usize = 5000; const NUM_ITER: usize = 2; -fn bench_axis_sort_random(c: &mut Criterion) { - c.benchmark_group("axis_sort_random") - .throughput(Throughput::Elements(SAMPLE_SIZE as u64)) - .bench_function("axis_sort_random", move |b| { - let sample_points = generator::uniform_rectangle( - Point2D::from([0., 0.]), - Point2D::from([30., 10.]), - SAMPLE_SIZE, - ); - let mut permutation: Vec<_> = (0..SAMPLE_SIZE).collect(); - - b.iter(|| axis_sort(&sample_points, &mut permutation, 0)) - }); -} - fn bench_raw_pdqsort_random(c: &mut Criterion) { c.benchmark_group("raw_pdqsort_random") .throughput(Throughput::Elements(SAMPLE_SIZE as u64)) @@ -87,20 +71,6 @@ fn bench_parallel_raw_pdqsort_sorted(c: &mut Criterion) { }); } -fn bench_axis_sort_sorted(c: &mut Criterion) { - c.benchmark_group("axis_sort_sorted") - .throughput(Throughput::Elements(SAMPLE_SIZE as u64)) - .bench_function("axis_sort_sorted", move |b| { - let sample_points = generator::already_x_sorted_rectangle( - Point2D::from([0., 0.]), - Point2D::from([30., 10.]), - SAMPLE_SIZE, - ); - let mut permutation: Vec<_> = (0..SAMPLE_SIZE).collect(); - b.iter(|| axis_sort(&sample_points, &mut permutation, 0)) - }); -} - fn bench_rcb_random(c: &mut Criterion) { c.benchmark_group("rcb_random") .throughput(Throughput::Elements(SAMPLE_SIZE as u64)) @@ -110,52 +80,54 @@ fn bench_rcb_random(c: &mut Criterion) { Point2D::from([30., 10.]), SAMPLE_SIZE, ); - let weights = vec![1.0; sample_points.len()]; - let mut ids = vec![0; sample_points.len()]; + let weights = vec![1; SAMPLE_SIZE]; + let mut ids = vec![0; SAMPLE_SIZE]; b.iter(|| { use rayon::iter::IntoParallelRefIterator as _; use rayon::iter::ParallelIterator as _; let p = sample_points.par_iter().cloned(); let w = weights.par_iter().cloned(); - rcb(&mut ids, p, w, NUM_ITER) + coupe::Rcb { + iter_count: NUM_ITER, + } + .partition(&mut ids, (p, w)) + .unwrap() }) }); } -fn bench_simplified_k_means(c: &mut Criterion) { - c.benchmark_group("simplified_k_means") +fn bench_k_means(c: &mut Criterion) { + c.benchmark_group("k_means") .throughput(Throughput::Elements(SAMPLE_SIZE as u64)) - .bench_function("simplified_k_means", move |b| { + .bench_function("k_means", move |b| { let sample_points = generator::uniform_rectangle( Point2D::new(0., 0.), Point2D::new(30., 10.), SAMPLE_SIZE, ); - let ids: Vec<_> = (0..SAMPLE_SIZE).collect(); - let weights: Vec<_> = ids.iter().map(|_| 1.).collect(); + let weights = vec![1.0; SAMPLE_SIZE]; b.iter(|| { - simplified_k_means( - &sample_points, - &weights, - 2usize.pow(NUM_ITER as u32), - 5., - 1000, - true, - ) + coupe::KMeans { + part_count: 2_usize.pow(NUM_ITER as u32), + imbalance_tol: 5.0, + delta_threshold: 0.0, + hilbert: true, + ..Default::default() + } + .partition(&mut [0; SAMPLE_SIZE], (&sample_points, &weights)) + .unwrap() }) }); } criterion_group!( benches, - bench_axis_sort_random, - bench_axis_sort_sorted, bench_raw_pdqsort_random, bench_parallel_raw_pdqsort_random, bench_raw_pdqsort_sorted, bench_parallel_raw_pdqsort_sorted, bench_rcb_random, - bench_simplified_k_means + bench_k_means ); criterion_main!(benches); diff --git a/benches/generator.rs b/benches/generator.rs index 66bfa5d5..34ed9987 100644 --- a/benches/generator.rs +++ b/benches/generator.rs @@ -1,5 +1,4 @@ use coupe::geometry::Point2D; -use itertools::Itertools; use rand::{self, Rng}; pub fn uniform_f64(min: f64, max: f64, num_vals: usize) -> Vec { @@ -18,20 +17,3 @@ pub fn uniform_rectangle(p_min: Point2D, p_max: Point2D, num_points: usize) -> V }) .collect() } - -pub fn already_x_sorted_rectangle( - p_min: Point2D, - p_max: Point2D, - num_points: usize, -) -> Vec { - let mut rng = rand::thread_rng(); - (0..num_points) - .map(|_| { - Point2D::from([ - rng.gen_range(p_min.x..p_max.x), - rng.gen_range(p_min.y..p_max.y), - ]) - }) - .sorted_by(|p1, p2| p1.x.partial_cmp(&p2.x).unwrap_or(std::cmp::Ordering::Equal)) - .collect() -} diff --git a/src/algorithms.rs b/src/algorithms.rs index 45d1df40..3a7c470b 100644 --- a/src/algorithms.rs +++ b/src/algorithms.rs @@ -1,12 +1,87 @@ -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; +use std::fmt; + +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; + +/// Common errors thrown by algorithms. +#[derive(Debug)] +#[non_exhaustive] +pub enum Error { + /// No partition that matches the given criteria could been found. + NotFound, + + /// Input sets don't have matching lengths. + InputLenMismatch { expected: usize, actual: usize }, +} + +impl fmt::Display for Error { + fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { + match self { + Error::NotFound => write!(f, "no partition found"), + Error::InputLenMismatch { expected, actual } => write!( + f, + "input sets don't have the same length (expected {expected} items, got {actual})", + ), + } + } +} + +impl std::error::Error for Error {} + +/// Map elements to parts randomly. +/// +/// # Example +/// +/// ```rust +/// use coupe::Partition as _; +/// use rand; +/// +/// let mut partition = [0; 12]; +/// +/// coupe::Random { rng: rand::thread_rng(), part_count: 3 } +/// .partition(&mut partition, ()) +/// .unwrap(); +/// ``` +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..9959dea8 100644 --- a/src/algorithms/ckk.rs +++ b/src/algorithms/ckk.rs @@ -1,3 +1,4 @@ +use super::Error; use num::FromPrimitive; use num::ToPrimitive; use std::iter::Sum; @@ -92,7 +93,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) -> Result<(), Error> where I: IntoIterator, T: Sum + Add + Sub, @@ -100,9 +101,14 @@ where T: Ord + Default + Copy, { let mut weights: Vec<(T, usize)> = weights.into_iter().zip(0..).collect(); - debug_assert_eq!(weights.len(), partition.len()); + if weights.len() != partition.len() { + return Err(Error::InputLenMismatch { + expected: partition.len(), + actual: weights.len(), + }); + } if weights.is_empty() { - return true; + return Ok(()); } weights.sort_unstable(); @@ -111,7 +117,48 @@ where let mut steps = Vec::new(); - ckk_bipart_rec(partition, &mut weights, tolerance, &mut steps) + if ckk_bipart_rec(partition, &mut weights, tolerance, &mut steps) { + Ok(()) + } else { + Err(Error::NotFound) + } +} + +/// The exact variant of the [Karmarkar-Karp][crate::KarmarkarKarp] algorithm. +/// +/// # Example +/// +/// ```rust +/// use coupe::Partition as _; +/// +/// let weights: [i32; 4] = [3, 5, 3, 9]; +/// let mut partition = [0; 4]; +/// +/// coupe::CompleteKarmarkarKarp { tolerance: 0.1 } +/// .partition(&mut partition, weights) +/// .unwrap(); +/// ``` +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 = Error; + + fn partition( + &mut self, + part_ids: &mut [usize], + weights: W, + ) -> Result { + ckk_bipart(part_ids, weights, self.tolerance) + } } // TODO tests diff --git a/src/algorithms/fiduccia_mattheyses.rs b/src/algorithms/fiduccia_mattheyses.rs index 088313ba..8cc31a40 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( + 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,107 @@ 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::Partition as _; +/// use coupe::Point2D; +/// use sprs::CsMat; +/// +/// // swap +/// // 0 1 0 1 +/// // +--+--+--+ +/// // | | | | +/// // +--+--+--+ +/// // 0 0 1 1 +/// let points = [ +/// 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 = [1.0; 8]; +/// let mut partition = [0, 0, 1, 1, 0, 1, 0, 1]; +/// +/// 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.); +/// +/// coupe::FiducciaMattheyses { max_bad_move_in_a_row: 1, ..Default::default() } +/// .partition(&mut partition, (adjacency.view(), &weights)) +/// .unwrap(); +/// +/// assert_eq!(partition[5], 0); +/// assert_eq!(partition[6], 1); +/// ``` +#[derive(Debug, Clone, Copy, Default)] +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..e6e558fc 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,81 @@ 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::Partition as _; +/// use coupe::Point2D; +/// use sprs::CsMat; +/// +/// // +--+--+--+ +/// // | | | | +/// // +--+--+--+ +/// +/// let weights = [1.0; 8]; +/// let mut partition = [0; 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.); +/// +/// coupe::GraphGrowth { part_count: 2 } +/// .partition(&mut partition, (adjacency.view(), &weights)) +/// .unwrap(); +/// ``` +#[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..e58fdfa4 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::Partition as _; +/// use coupe::Real; +/// +/// let weights = [3.2, 6.8, 10.0, 7.5].map(Real::from); +/// let mut partition = [0; 4]; +/// +/// coupe::Greedy { part_count: 2 } +/// .partition(&mut partition, weights) +/// .unwrap(); +/// ``` +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..5f91d359 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,73 @@ 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::Partition as _; +/// use coupe::Point2D; +/// +/// let points = [ +/// 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 = [1.0; 8]; +/// let mut partition = [0; 8]; +/// +/// // generate a partition of 4 parts +/// coupe::HilbertCurve { part_count: 4, order: 5 } +/// .partition(&mut partition, (points, weights)) +/// .unwrap(); +/// +/// assert_eq!(partition[0], partition[1]); +/// assert_eq!(partition[2], partition[3]); +/// assert_eq!(partition[4], partition[5]); +/// assert_eq!(partition[6], partition[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..bc8c26d4 100644 --- a/src/algorithms/k_means.rs +++ b/src/algorithms/k_means.rs @@ -15,9 +15,6 @@ use rayon::prelude::*; use std::cmp::Ordering; use std::sync::atomic::{self, AtomicPtr}; -// use super::hilbert_curve; -use super::z_curve; - use itertools::iproduct; use itertools::Itertools; @@ -26,143 +23,6 @@ use itertools::Itertools; /// for the k-means algorithm and not a partition id type ClusterId = usize; -/// A simplified implementation of the algorithm described in the paper -/// by Moritz von Looz et al. that follows the same idea but without the small -/// optimizations that would improve the efficiency of the algorithm. In particular, -/// 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( - points: &[PointND], - weights: &[f64], - num_partitions: usize, - imbalance_tol: f64, - mut n_iter: isize, - hilbert: bool, -) -> Vec -where - Const: DimSub>, - DefaultAllocator: Allocator, Const, Buffer = ArrayStorage> - + Allocator, Const<1>>>, -{ - let sfc_order = (f64::from(std::u32::MAX)).log(f64::from(2u32.pow(D as u32))) as u32; - let permu = if hilbert { - unimplemented!("hilbert curve currently not implemented for n-dimension"); - // hilbert_curve::hilbert_curve_reorder(points, 15) - } else { - z_curve::z_curve_reorder(points, sfc_order) - }; - - let points_per_center = points.len() / num_partitions; - - let mut centers: Vec<_> = permu - .iter() - .cloned() - .skip(points_per_center / 2) - .step_by(points_per_center) - .map(|idx| points[idx]) - .collect(); - - let center_ids: Vec = centers.par_iter().map(|_| crate::uid()).collect(); - - let mut influences = centers.par_iter().map(|_| 1.).collect::>(); - - let dummy_id = crate::uid(); - let mut assignments: Vec = permu.par_iter().map(|_| dummy_id).collect(); - let atomic_handle = AtomicPtr::from(assignments.as_mut_ptr()); - permu - .par_chunks(points_per_center) - .zip(center_ids.par_iter()) - .for_each(|(chunk, id)| { - let ptr = atomic_handle.load(atomic::Ordering::Relaxed); - for idx in chunk { - unsafe { std::ptr::write(ptr.add(*idx), *id) } - } - }); - - let mut imbalance = std::f64::MAX; - - let target_weight = weights.par_iter().sum::() / num_partitions as f64; - - while imbalance > imbalance_tol && n_iter > 0 { - n_iter -= 1; - - // find new assignments - permu.par_iter().for_each(|&idx| { - // find closest center - let mut distances = ::std::iter::repeat(points[idx]) - .zip(centers.iter()) - .zip(center_ids.iter()) - .zip(influences.iter()) - .map(|(((p, center), id), influence)| (*id, (p - center).norm() * influence)) - .collect::>(); - - distances - .as_mut_slice() - .sort_unstable_by(|(_, d1), (_, d2)| d1.partial_cmp(d2).unwrap_or(Ordering::Equal)); - - // update assignment with closest center found - - let new_assignment = distances.into_iter().next().unwrap().0; - let ptr = atomic_handle.load(atomic::Ordering::Relaxed); - unsafe { std::ptr::write(ptr.add(idx), new_assignment) } - }); - - // update centers position from new assignments - centers - .par_iter_mut() - .zip(center_ids.par_iter()) - .for_each(|(center, id)| { - let new_center = geometry::center( - &permu - .iter() - // TODO: maybe cloning is avoidable here - .map(|&idx| (points[idx], assignments[idx])) - .filter(|(_, point_id)| *id == *point_id) - .map(|(p, _)| p) - .collect::>(), - ); - - *center = new_center; - }); - - // compute cluster weights - let cluster_weights = center_ids - .par_iter() - .map(|id| { - assignments - .iter() - .zip(weights.iter()) - .filter(|(point_id, _)| *id == **point_id) - .fold(0., |acc, (_, weight)| acc + weight) - }) - .collect::>(); - - // update influence - cluster_weights - .par_iter() - .zip(influences.par_iter_mut()) - .for_each(|(cluster_weight, influence)| { - let ratio = target_weight / *cluster_weight as f64; - - let new_influence = *influence / ratio.sqrt(); - let max_diff = 0.05 * *influence; - if (*influence - new_influence).abs() < max_diff { - *influence = new_influence; - } else if new_influence > *influence { - *influence += max_diff; - } else { - *influence -= max_diff; - } - }); - - // update imbalance - imbalance = self::imbalance(&cluster_weights); - } - - assignments -} - fn imbalance(weights: &[f64]) -> f64 { match ( weights @@ -216,89 +76,7 @@ impl Default for BalancedKmeansSettings { } } -pub fn balanced_k_means( - points: &[PointND], - weights: &[f64], - settings: impl Into>, -) -> Vec -where - Const: DimSub>, - DefaultAllocator: Allocator, Const, Buffer = ArrayStorage> - + Allocator, Const<1>>>, -{ - let settings = settings.into().unwrap_or_default(); - - // sort points with space filling curve - let sfc_order = (f64::from(std::u32::MAX)).log(f64::from(2u32.pow(D as u32))) as u32; - let mut permu = if settings.hilbert { - unimplemented!("hilbert curve currently not implemented for n-dimension"); - // hilbert_curve::hilbert_curve_reorder(points, 15) - } else { - z_curve::z_curve_reorder(points, sfc_order) - }; - - // Compute how many points will be initially assigned to each cluster - let points_per_center = points.len() / settings.num_partitions; - - // select num_partitions initial centers from the ordered points - let centers: Vec<_> = permu - .iter() - // for each partition yielded by the Z-curve reordering - // we select the median point to be the initial cluster center - // because it is in most cases in the middle of the partition - .skip(points_per_center / 2) - .step_by(points_per_center) - .par_bridge() - .map(|&idx| points[idx]) - .collect(); - - // generate unique ids for each initial partition that will live throughout - // the algorithm (no new id is generated afterwards) - let center_ids: Vec = centers.par_iter().map(|_| crate::uid()).collect(); - - let dummy_id = crate::uid(); - let mut assignments: Vec<_> = permu.par_iter().map(|_| dummy_id).collect(); - let atomic_handle = AtomicPtr::from(assignments.as_mut_ptr()); - permu - .par_chunks(points_per_center) - .zip(center_ids.par_iter()) - .for_each(|(chunk, id)| { - let ptr = atomic_handle.load(atomic::Ordering::Relaxed); - for idx in chunk { - unsafe { std::ptr::write(ptr.add(*idx), *id) } - } - }); - - debug_assert_eq!(assignments.len(), points.len()); - - // Generate initial influences (to 1) - let mut influences: Vec<_> = centers.par_iter().map(|_| 1.).collect(); - - // Generate initial lower and upper bounds. These two variables represent bounds on - // the effective distance between an point and the cluster it is assigned to. - let mut lbs: Vec<_> = points.par_iter().map(|_| 0.).collect(); - let mut ubs: Vec<_> = points.par_iter().map(|_| std::f64::MAX).collect(); // we use f64::MAX to represent infinity - - balanced_k_means_iter( - Inputs { points, weights }, - Clusters { - centers, - center_ids: ¢er_ids, - }, - &mut permu, - AlgorithmState { - assignments: &mut assignments, - influences: &mut influences, - lbs: &mut lbs, - ubs: &mut ubs, - }, - &settings, - settings.max_iter, - ); - 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 +515,107 @@ 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::Partition as _; +/// use coupe::Point2D; +/// +/// let points = [ +/// 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 = [1.; 9]; +/// +/// // create an unbalanced partition: +/// // - p1: total weight = 1 +/// // - p2: total weight = 7 +/// // - p3: total weight = 1 +/// let mut partition = [1, 2, 2, 2, 2, 2, 2, 2, 3]; +/// +/// coupe::KMeans { part_count: 3, delta_threshold: 0.0, ..Default::default() } +/// .partition(&mut partition, (&points, &weights)) +/// .unwrap(); +/// +/// assert_eq!(partition[0], partition[1]); +/// assert_eq!(partition[0], partition[2]); +/// +/// assert_eq!(partition[3], partition[4]); +/// assert_eq!(partition[3], partition[5]); +/// +/// assert_eq!(partition[6], partition[7]); +/// assert_eq!(partition[6], partition[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, weights, settings, part_ids); + Ok(()) + } +} diff --git a/src/algorithms/kernighan_lin.rs b/src/algorithms/kernighan_lin.rs index 31654f6d..fe041e0e 100644 --- a/src/algorithms/kernighan_lin.rs +++ b/src/algorithms/kernighan_lin.rs @@ -3,18 +3,16 @@ //! 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>, +fn kernighan_lin( + part_ids: &mut [usize], + weights: &[f64], adjacency: CsMatView, - max_passes: impl Into>, - max_flips_per_pass: impl Into>, - max_imbalance_per_flip: impl Into>, + max_passes: Option, + max_flips_per_pass: Option, + max_imbalance_per_flip: Option, max_bad_move_in_a_row: usize, ) { // To adapt Kernighan-Lin to a partition of more than 2 parts, @@ -22,15 +20,10 @@ pub(crate) fn kernighan_lin<'a, const D: usize>( // are adjacent if there exists an element in one part that is linked to // an element in the other part). - 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, + adjacency, max_passes, max_flips_per_pass, max_imbalance_per_flip, @@ -38,10 +31,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, @@ -57,7 +50,7 @@ fn kernighan_lin_2_impl( unimplemented!(); } - let mut cut_size = crate::topology::cut_size(adjacency.view(), initial_partition); + let mut cut_size = crate::topology::cut_size(adjacency, initial_partition); println!("Initial cut size: {}", cut_size); let mut new_cut_size = cut_size; @@ -180,3 +173,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::Partition as _; +/// use coupe::Point2D; +/// use sprs::CsMat; +/// +/// // swap +/// // 0 1 0 1 +/// // +--+--+--+ +/// // | | | | +/// // +--+--+--+ +/// // 0 0 1 1 +/// let points = [ +/// 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 = [1.; 8]; +/// let mut partition = [0, 0, 1, 1, 0, 1, 0, 1]; +/// +/// 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 +/// coupe::KernighanLin { +/// max_passes: Some(1), +/// max_flips_per_pass: Some(1), +/// max_imbalance_per_flip: None, +/// max_bad_move_in_a_row: 1, +/// } +/// .partition(&mut partition, (adjacency.view(), &weights)) +/// .unwrap(); +/// +/// assert_eq!(partition[5], 0); +/// assert_eq!(partition[6], 1); +/// ``` +#[derive(Debug, Clone, Copy, Default)] +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..944fe701 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,40 @@ where parts.truncate(partition.len()); partition.copy_from_slice(&parts); } + +/// The Karmarkar-Karp algorithm, also called the Largest Differencing Method. +/// +/// # Example +/// +/// ```rust +/// use coupe::Partition as _; +/// +/// let weights = [3, 5, 3, 9]; +/// let mut partition = [0; 4]; +/// +/// coupe::KarmarkarKarp { part_count: 3 } +/// .partition(&mut partition, weights) +/// .unwrap(); +/// ``` +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..cfd3acaf 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,73 @@ 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::Partition as _; +/// use coupe::Point2D; +/// +/// 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 = [4.2; 9]; +/// let mut partition = [0; 9]; +/// +/// // generate a partition of 4 parts +/// coupe::MultiJagged { part_count: 9, max_iter: 4 } +/// .partition(&mut partition, (&points, &weights)) +/// .unwrap(); +/// +/// for i in 0..9 { +/// for j in 0..9 { +/// if j == i { +/// continue +/// } +/// assert_ne!(partition[i], partition[j]) +/// } +/// } +/// ``` +pub struct MultiJagged { + pub part_count: 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, weights, self.part_count, 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 f4a8ccb9..326af565 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, @@ -547,6 +548,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::Partition as _; +/// use coupe::Point2D; +/// +/// let points = [ +/// Point2D::new(1., 1.), +/// Point2D::new(-1., 1.), +/// Point2D::new(1., -1.), +/// Point2D::new(-1., -1.), +/// ]; +/// let weights = [1; 4]; +/// let mut partition = [0; 4]; +/// +/// // generate a partition of 4 parts +/// coupe::Rcb { iter_count: 2 } +/// .partition(&mut partition, (points, weights)) +/// .unwrap(); +/// +/// for i in 0..4 { +/// for j in 0..4 { +/// if j == i { +/// continue +/// } +/// assert_ne!(partition[i], partition[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], @@ -580,12 +647,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>>>, @@ -608,6 +671,74 @@ 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::Partition as _; +/// use coupe::Point2D; +/// +/// // Here, the inertia axis is the y axis. +/// // We thus expect Rib to split horizontally first. +/// let points = [ +/// Point2D::new(1., 10.), +/// Point2D::new(-1., 10.), +/// Point2D::new(1., -10.), +/// Point2D::new(-1., -10.), +/// ]; +/// let weights = [1; 4]; +/// let mut partition = [0; 4]; +/// +/// // generate a partition of 2 parts (1 split) +/// coupe::Rib { iter_count: 1 } +/// .partition(&mut partition, (&points, weights)) +/// .unwrap(); +/// +/// // the two points at the top are in the same part +/// assert_eq!(partition[0], partition[1]); +/// +/// // the two points at the bottom are in the same part +/// assert_eq!(partition[2], partition[3]); +/// +/// // there are two different parts +/// assert_ne!(partition[1], partition[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 3c519413..e8054533 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,10 +7,10 @@ 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, T: Ord + Copy, { @@ -21,22 +20,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 +117,52 @@ 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::Partition as _; +/// use rand; +/// +/// let part_count = 2; +/// let mut partition = [0; 4]; +/// let weights = [4, 6, 2, 9]; +/// +/// coupe::Random { rng: rand::thread_rng(), part_count } +/// .partition(&mut partition, ()) +/// .unwrap(); +/// coupe::VnBest { part_count } +/// .partition(&mut partition, weights) +/// .unwrap(); +/// ``` +pub struct VnBest { + pub part_count: usize, +} + +impl crate::Partition for VnBest +where + W: IntoIterator, + W::Item: AddAssign + Sub + Div + Mul, + 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 +178,7 @@ mod tests { let w = [1, 2, 3, 4, 5, 6]; let mut part = vec![0; w.len()]; let imb_ini = imbalance::imbalance(2, &part, &w); - vn_best_mono::(&mut part, &w, 2); + vn_best_mono(&mut part, w, 2); let imb_end = imbalance::imbalance(2, &part, &w); assert!(imb_end < imb_ini); println!("imbalance : {} < {}", imb_end, imb_ini); @@ -144,9 +196,9 @@ mod tests { prop::collection::vec(0..1usize, num_weights)) }) ) { - let imb_ini = imbalance::max_imbalance(2, &partition, &weights); - vn_best_mono::(&mut partition, &weights, 2); - let imb_end = imbalance::max_imbalance(2, &partition, &weights); + let imb_ini = imbalance::max_imbalance(2, &partition, weights.iter().cloned()); + vn_best_mono::<_, u64>(&mut partition, weights.iter().cloned(), 2); + let imb_end = imbalance::max_imbalance(2, &partition, weights.iter().cloned()); // 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..98492a84 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,63 @@ 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::Partition as _; +/// use coupe::Point2D; +/// +/// let points = [ +/// 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 mut partition = [0; 8]; +/// +/// // generate a partition of 4 parts +/// coupe::ZCurve { part_count: 4, order: 5 } +/// .partition(&mut partition, &points) +/// .unwrap(); +/// +/// assert_eq!(partition[0], partition[1]); +/// assert_eq!(partition[2], partition[3]); +/// assert_eq!(partition[4], partition[5]); +/// assert_eq!(partition[6], partition[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 +250,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 +261,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/geometry.rs b/src/geometry.rs index 67ea48b0..6bd09e11 100644 --- a/src/geometry.rs +++ b/src/geometry.rs @@ -1,7 +1,6 @@ //! A few useful geometric types use approx::Ulps; -use itertools::Itertools; use nalgebra::allocator::Allocator; use nalgebra::ArrayStorage; use nalgebra::Const; @@ -125,24 +124,6 @@ impl Aabb { Self { p_min, p_max } } - pub fn aspect_ratio(&self) -> f64 { - use itertools::MinMaxResult::*; - - let (min, max) = match self - .p_min() - .iter() - .zip(self.p_max().iter()) - .map(|(min, max)| (max - min).abs()) - .minmax() - { - MinMax(min, max) => (min, max), - _ => unimplemented!(), - }; - - // TODO: What if min == 0.0 ? - max / min - } - pub fn contains(&self, point: &PointND) -> bool { let eps = 10. * std::f64::EPSILON; self.p_min @@ -261,14 +242,6 @@ impl Mbr { } } - /// Computes the aspect ratio of the Mbr. - /// - /// It is defined as follows: - /// `ratio = long_side / short_side` - pub fn aspect_ratio(&self) -> f64 { - self.aabb.aspect_ratio() - } - /// Computes the distance between a point and the current mbr. pub fn distance_to_point(&self, point: &PointND) -> f64 { self.aabb.distance_to_point(&(self.mbr_to_aabb * point)) @@ -377,6 +350,7 @@ mod tests { use super::*; use approx::assert_relative_eq; use approx::assert_ulps_eq; + use itertools::Itertools as _; use nalgebra::Matrix2; #[test] @@ -390,11 +364,9 @@ mod tests { ]; let aabb = Aabb::from_points(&points); - let aspect_ratio = aabb.aspect_ratio(); assert_ulps_eq!(aabb.p_min, Point2D::from([0., 0.])); assert_ulps_eq!(aabb.p_max, Point2D::from([5., 5.])); - assert_ulps_eq!(aspect_ratio, 1.); } #[test] @@ -411,21 +383,6 @@ mod tests { let _aabb = Aabb::from_points(&points); } - #[test] - fn test_mbr_2d() { - let points = vec![ - Point2D::from([5., 3.]), - Point2D::from([0., 0.]), - Point2D::from([1., -1.]), - Point2D::from([4., 4.]), - ]; - - let mbr = Mbr::from_points(&points); - let aspect_ratio = mbr.aspect_ratio(); - - assert_relative_eq!(aspect_ratio, 4., epsilon = 1e-14); - } - #[test] fn test_inertia_matrix() { let points = vec![ @@ -567,31 +524,9 @@ mod tests { ]; let aabb = Aabb::from_points(&points); - let aspect_ratio = aabb.aspect_ratio(); assert_ulps_eq!(aabb.p_min, Point3D::from([0., 0., -2.])); assert_ulps_eq!(aabb.p_max, Point3D::from([5., 5., 5.])); - assert_ulps_eq!(aspect_ratio, 7. / 5.); - } - - #[test] - fn test_mbr_3d() { - let points = vec![ - Point3D::from([5., 3., 0.]), - Point3D::from([0., 0., 0.]), - Point3D::from([1., -1., 0.]), - Point3D::from([4., 4., 0.]), - Point3D::from([5., 3., 2.]), - Point3D::from([0., 0., 2.]), - Point3D::from([1., -1., 2.]), - Point3D::from([4., 4., 2.]), - ]; - - let mbr = Mbr::from_points(&points); - println!("mbr = {:?}", mbr); - let aspect_ratio = mbr.aspect_ratio(); - - assert_relative_eq!(aspect_ratio, 4., epsilon = 1e-14); } #[test] diff --git a/src/imbalance.rs b/src/imbalance.rs index 1fad5f4f..ba4ef189 100644 --- a/src/imbalance.rs +++ b/src/imbalance.rs @@ -12,14 +12,14 @@ 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() .zip(weights) .fold(vec![T::zero(); num_parts], |mut acc, (&part, w)| { - acc[part] += w.clone(); + acc[part] += w; acc }) } @@ -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..8e3e2b66 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,1324 +25,39 @@ //! - [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; mod real; -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; } diff --git a/src/partition.rs b/src/partition.rs deleted file mode 100644 index d7b57e12..00000000 --- a/src/partition.rs +++ /dev/null @@ -1,442 +0,0 @@ -//! Utilities to manipulate partitions. - -use itertools::Itertools; -use nalgebra::allocator::Allocator; -use nalgebra::ArrayStorage; -use nalgebra::Const; -use nalgebra::DefaultAllocator; -use nalgebra::DimDiff; -use nalgebra::DimSub; -use num::{Num, Signed}; -use sprs::CsMatView; - -use std::cmp::PartialOrd; -use std::iter::Sum; - -use crate::geometry::Mbr; -use crate::PointND; - -/// Represents a partition. -/// -/// This struct is usually created by a partitioning algorithm. It internally uses unique IDs to represent a partition of a set of points and weights. -/// The `Partition` object exposes a convenient interface to manipulate a partition: -/// - Iterate over parts -/// - Analyze quality (e.g. weight imbalance) -/// -/// # Example -/// -/// ```rust -/// use coupe::Point2D; -/// use coupe::partition::Partition; -/// use approx::*; -/// -/// -/// let ids = vec![1, 2, 3, 1, 2, 3, 1, 2, 3]; -/// let weights = vec![1.; 9]; -/// let points = vec![ -/// Point2D::new(2., 2.), -/// Point2D::new(3., 3.), -/// Point2D::new(7., 7.), -/// Point2D::new(8., 8.), -/// Point2D::new(9., 9.), -/// Point2D::new(5., 5.), -/// Point2D::new(1., 1.), -/// Point2D::new(6., 6.), -/// Point2D::new(4., 4.), -/// ]; -/// -/// let partition = Partition::from_ids(&points, &weights, ids); -/// -/// // iterate over each part corresponding to a unique ID -/// for part in partition.parts() { -/// assert_ulps_eq!(part.total_weight(), 3.) -/// } -/// ``` -#[derive(Clone)] -pub struct Partition<'a, P, W> { - points: &'a [P], - weights: &'a [W], - ids: Vec, -} - -impl<'a, P, W> Partition<'a, P, W> { - /// Constructs a default partition from a set of points and weights and initialize it with default IDs (with a single part). - /// - /// # Panics - /// - /// Panics if `points` and `weights` have different lenghts. - pub fn new(points: &'a [P], weights: &'a [W]) -> Self { - if points.len() != weights.len() { - panic!( - "Cannot initialize a partition with points and weights of different sizes. Found {} points and {} weights", - points.len(), weights.len() - ); - } - - Self { - points, - weights, - ids: vec![0; points.len()], - } - } - - /// Constructs a partition from a set of points, weights and IDs. - /// - /// # Panics - /// - /// Panics if `points`, `weights` and `ids` do not have the same length. - pub fn from_ids(points: &'a [P], weights: &'a [W], ids: Vec) -> Self { - if points.len() != weights.len() { - panic!( - "Cannot initialize a partition with points and weights of different sizes. Found {} points and {} weights", - points.len(), weights.len() - ); - } - - if points.len() != ids.len() { - panic!( - "Cannot initialize a partition with points and ids of different sizes. Found {} points and {} ids", - points.len(), ids.len() - ); - } - - Self { - points, - weights, - ids, - } - } - - /// Consumes the partition, returning the internal array of ids representing the partition. - pub fn into_ids(self) -> Vec { - self.ids - } - - /// Consumes the partition, returning the internal components representing the partition - /// - /// # Example - /// - /// ```rust - /// use coupe::partition::Partition; - /// use coupe::Point2D; - /// - /// fn consume_partition(partition: Partition) { - /// let (_points, _weights, mut _ids) = partition.into_raw(); - /// // do something with the fields - /// } - /// ``` - pub fn into_raw(self) -> (&'a [P], &'a [W], Vec) { - (self.points, self.weights, self.ids) - } - - pub fn as_raw(&self) -> (&'a [P], &'a [W], &[usize]) { - (self.points, self.weights, &self.ids) - } - - pub fn as_raw_mut(&mut self) -> (&'a [P], &'a [W], &mut [usize]) { - (self.points, self.weights, &mut self.ids) - } - - /// Returns a slice of the IDs used to represent the partition. - pub fn ids(&self) -> &[usize] { - &self.ids - } - - /// Returns of mutable slice of the IDs used to represent the partition. - pub fn ids_mut(&mut self) -> &mut [usize] { - &mut self.ids - } - - /// Returns a slice of the points used for the partition. - pub fn points(&self) -> &[P] { - self.points - } - - /// Returns a slice of the weights used for the partition. - pub fn weights(&self) -> &[W] { - self.weights - } - - /// Returns an iterator over each part of the partition. - /// - /// Each part is a set of points and weights that share the same ID. - pub fn parts(&'a self) -> impl Iterator> { - let indices = (0..self.points.len()).collect::>(); - self.ids.iter().unique().map(move |id| { - let indices = indices - .iter() - .filter(|idx| self.ids[**idx] == *id) - .cloned() - .collect::>(); - Part::<'a, P, W> { - partition: self, - indices, - } - }) - } - - /// Computes the maximum imbalance of the partition. - /// It is defined as the maximum difference between the weights of two different parts. - /// - /// # Example - /// - /// ```rust - /// use coupe::partition::Partition; - /// - /// let dummy_point = coupe::Point2D::new(0., 0.); - /// let points = [dummy_point; 5]; - /// - /// let weights = [1, 2, 3, 4, 5]; - /// - /// let ids = vec![1, 2, 2, 2, 3]; - /// - /// let partition = Partition::from_ids(&points[..], &weights[..], ids); - /// - /// let max_imbalance = partition.max_imbalance(); - /// assert_eq!(max_imbalance, 8); - /// ``` - pub fn max_imbalance(&'a self) -> W - where - W: Num + PartialOrd + Signed + Sum + Clone, - { - let total_weights = self - .parts() - .map(|part| part.total_weight()) - .collect::>(); - total_weights - .iter() - .flat_map(|w1| { - total_weights - .iter() - .map(move |w2| (w1.clone() - w2.clone()).abs()) - }) - .max_by(|a, b| a.partial_cmp(b).unwrap()) - // if the partition is empty, then there is the imbalance is null - .unwrap_or_else(W::zero) - } - - /// Computes the relative imbalance of the partition. - /// It is defined as follows: `relative_imbalance = maximum_imbalance / total_weight_in_partition`. - /// It expresses the imbalance in terms of percentage of the total weight. It may be more meaningful than raw numbers in some cases. - /// - /// # Example - /// - /// ```rust - /// use coupe::partition::Partition; - /// - /// let dummy_point = coupe::Point2D::new(0., 0.); - /// let points = [dummy_point; 5]; - /// - /// let weights = [1000, 1000, 1000, 1000, 6000]; - /// - /// let ids = vec![1, 2, 2, 2, 3]; - /// - /// let partition = Partition::from_ids(&points[..], &weights[..], ids); - /// - /// let relative_imbalance = partition.relative_imbalance(); - /// assert_eq!(relative_imbalance, 0.5); // = 50% of total weight - /// ``` - pub fn relative_imbalance(&'a self) -> f64 - where - W: Num + PartialOrd + Signed + Sum + Clone, - f64: std::convert::From, - { - let total_weights = self - .parts() - .map(|part| part.total_weight()) - .collect::>(); - let max_imbalance = total_weights - .iter() - .flat_map(|w1| { - total_weights - .iter() - .map(move |w2| (w1.clone() - w2.clone()).abs()) - }) - .max_by(|a, b| a.partial_cmp(b).unwrap()) - .unwrap_or_else(W::zero); - - f64::from(max_imbalance) / f64::from(total_weights.into_iter().sum()) - } - - /// Merge two parts of the partition (identified by unique id). - /// - /// If either `id1` or `id2` is not contained in the partition, does nothing. - pub fn merge_parts(&mut self, id1: usize, id2: usize) { - for id in self.ids.iter_mut() { - if *id == id2 { - *id = id1; - } - } - } - - /// Construct pairs of adjacent parts. - /// Two parts are adjacent if there exist one element in the first - /// part and one element in the second part that are linked together - /// in the adjacency graph. - pub fn adjacent_parts( - &'a self, - adjacency: CsMatView, - ) -> Vec<(Part<'a, P, W>, Part<'a, P, W>)> { - // build parts adjacency graph - let parts = self.parts().collect::>(); - - let mut adjacent_parts = Vec::new(); - - for (i, p) in parts.iter().enumerate() { - for (j, q) in parts.iter().enumerate() { - if j < i - && adjacency - .iter() - .any(|(_, (i, j))| p.indices().contains(&i) && q.indices().contains(&j)) - { - adjacent_parts.push((p.clone(), q.clone())) - } - } - } - - adjacent_parts - } - - /// Swaps the parts of two elements - /// - /// Panics if `idx1` or `idx2` is out of range - pub fn swap(&mut self, idx1: usize, idx2: usize) { - self.ids.swap(idx1, idx2); - } -} - -/// Represent a part of a partition. -/// -/// This struct is not meaningful on its own, and is usually constructed by -/// [iterating over a partition](struct.Partition.html#method.parts). -pub struct Part<'a, P, W> { - partition: &'a Partition<'a, P, W>, - indices: Vec, -} - -// Custom Clone impl since we juste clone the indices and not the partition itself -impl<'a, P, W> Clone for Part<'a, P, W> { - fn clone(&self) -> Self { - Self { - partition: self.partition, - indices: self.indices.clone(), - } - } -} - -impl<'a, P, W> Part<'a, P, W> { - /// Computes the total weight (i.e. the sum of all the weights) of the part. - pub fn total_weight(&self) -> W - where - W: Sum + Clone, - { - self.indices - .iter() - .map(|idx| &self.partition.weights[*idx]) - .cloned() - .sum() - } - - pub fn indices(&self) -> &[usize] { - &self.indices - } - - pub fn into_indices(self) -> Vec { - self.indices - } - - pub fn points(&self) -> &'a [P] { - self.partition.points() - } - - pub fn weights(&self) -> &'a [W] { - self.partition.weights() - } - - /// Iterate over the points and weights of the part. - /// - /// # Example - /// - /// ```rust - /// use coupe::Point2D; - /// use coupe::partition::Partition; - /// - /// fn iterate(partition: &Partition) { - /// // iterate over each part - /// for part in partition.parts() { - /// //iterate over each point/weight of the current part - /// for (p, w) in part.iter() { - /// // do something with p, w - /// } - /// } - /// } - /// ``` - pub fn iter(&self) -> PartIter { - PartIter { - partition: self.partition, - indices: &self.indices, - } - } -} - -impl<'a, W, const D: usize> Part<'a, PointND, W> { - /// Computes the aspect ratio of a part. It is defined as the aspect ratio of a minimal bounding rectangle - /// of the set of points contained in the part. - /// - /// # Example - /// - /// ```rust - /// use coupe::Point2D; - /// use coupe::partition::Partition; - /// use approx::*; - /// - /// let points = [ - /// Point2D::new(0., 0.), - /// Point2D::new(6., 0.), - /// Point2D::new(0., 2.), - /// Point2D::new(6., 2.), - /// ]; - /// let weights = [1.; 4]; - /// - /// let partition = Partition::new(&points[..], &weights[..]); // initialized with a single part - /// - /// // only 1 iteration - /// for part in partition.parts() { - /// assert_ulps_eq!(part.aspect_ratio(), 3.); - /// } - /// ``` - pub fn aspect_ratio(&self) -> f64 - where - Const: DimSub>, - DefaultAllocator: Allocator, Const, Buffer = ArrayStorage> - + Allocator, Const<1>>>, - { - if self.indices.len() <= 2 { - panic!("Cannot compute the aspect ratio of a part of less than 2 points"); - } - let points = self - .indices - .iter() - .map(|&idx| self.partition.points()[idx]) - .collect::>(); - Mbr::from_points(&points).aspect_ratio() - } -} - -/// An iterator over points and weights of a part. -/// -/// Is struct is usually created when [iterating over a part](struct.Part.html#method.iter). -pub struct PartIter<'a, P, W> { - partition: &'a Partition<'a, P, W>, - indices: &'a [usize], -} - -impl<'a, P, W> Iterator for PartIter<'a, P, W> { - type Item = (&'a P, &'a W); - fn next(&mut self) -> Option { - self.indices.get(0).map(|idx| { - self.indices = &self.indices[1..]; - (&self.partition.points[*idx], &self.partition.weights[*idx]) - }) - } -} diff --git a/src/run_info.rs b/src/run_info.rs deleted file mode 100644 index 457157d3..00000000 --- a/src/run_info.rs +++ /dev/null @@ -1,16 +0,0 @@ -/// Information on an algorithm run. -/// -/// Filled in by algorithms when run on a given input. Gives information about how the run went. -#[derive(Default, Copy, Clone)] -pub struct RunInfo { - /// Number of iterations the algorithm underwent to provide a partition. - pub algo_iterations: Option, -} - -impl RunInfo { - pub fn skip() -> RunInfo { - RunInfo { - algo_iterations: Some(0), - } - } -} diff --git a/tools/Cargo.toml b/tools/Cargo.toml index 4c1133d0..77a41c52 100644 --- a/tools/Cargo.toml +++ b/tools/Cargo.toml @@ -30,3 +30,4 @@ anyhow = { version = "1", default-features = false, features = ["std"] } # Other utilities itertools = { version = "0.10", default-features = false } num = { version = "0.4", default-features = false } +rayon = "1.3" diff --git a/tools/num-part/src/main.rs b/tools/num-part/src/main.rs index db40ade4..9a7a88aa 100644 --- a/tools/num-part/src/main.rs +++ b/tools/num-part/src/main.rs @@ -1,6 +1,5 @@ use anyhow::Context as _; use anyhow::Result; -use coupe::RunInfo; use rand::SeedableRng as _; use std::env; use std::path::PathBuf; @@ -73,7 +72,7 @@ where /// - `Err(msg)`, where `msg` is the reason the algorithm cannot be run with the given arguments /// (e.g. it's a bi-partitioning algorithm but the caller passed a number of parts that is /// different than 2). -type Algorithm = Box Result>; +type Algorithm = Box Result<()>>; fn main() -> Result<()> { let mut options = getopts::Options::new(); diff --git a/tools/src/bin/mesh-part/main.rs b/tools/src/bin/mesh-part/main.rs index 6672aa12..c8fd4ce6 100644 --- a/tools/src/bin/mesh-part/main.rs +++ b/tools/src/bin/mesh-part/main.rs @@ -1,77 +1,127 @@ use anyhow::Context as _; use anyhow::Result; -use coupe::RunInfo; +use coupe::Partition as _; +use coupe::PointND; +use mesh_io::medit::Mesh; use mesh_io::weight; +use rayon::iter::IntoParallelRefIterator as _; +use rayon::iter::ParallelIterator as _; use std::env; use std::fs; use std::io; +use std::mem; -macro_rules! coupe_part { - ( $algo:expr, $dimension:expr, $points:expr, $weights:expr ) => {{ - match $dimension { - 2 => { - coupe::Partitioner::<2>::partition(&$algo, $points.as_slice(), $weights).into_ids() +struct Problem { + points: Vec>, + weights: weight::Array, +} + +trait Algorithm { + fn run(&mut self, partition: &mut [usize], problem: &Problem) -> Result<()>; +} + +impl Algorithm for coupe::Random +where + R: rand::Rng, +{ + fn run(&mut self, partition: &mut [usize], _: &Problem) -> Result<()> { + self.partition(partition, ())?; + Ok(()) + } +} + +impl Algorithm for coupe::Greedy { + fn run(&mut self, partition: &mut [usize], problem: &Problem) -> Result<()> { + use weight::Array::*; + match &problem.weights { + Integers(is) => { + let weights = is.iter().map(|weight| weight[0]); + self.partition(partition, weights)?; } - 3 => { - coupe::Partitioner::<3>::partition(&$algo, $points.as_slice(), $weights).into_ids() + Floats(fs) => { + let weights = fs.iter().map(|weight| coupe::Real::from(weight[0])); + self.partition(partition, weights)?; } - _ => anyhow::bail!("only 2D and 3D meshes are supported"), } - }}; - ( $algo:expr, $problem:expr ) => {{ - let weights: Vec = match &$problem.weights { - weight::Array::Floats(ws) => ws.iter().map(|weight| weight[0]).collect(), - weight::Array::Integers(_) => anyhow::bail!("integers are not supported by coupe yet"), - }; - coupe_part!($algo, $problem.dimension, $problem.points, &weights) - }}; + Ok(()) + } } -macro_rules! coupe_part_improve { - ( $algo:expr, $partition:expr, $dimension:expr, $points:expr, $weights:expr ) => {{ - let ids = match $dimension { - 2 => { - let points = coupe::PointsView::to_points_nd($points.as_slice()); - let partition = - coupe::partition::Partition::from_ids(points, $weights, $partition.to_vec()); - coupe::PartitionImprover::<2>::improve_partition(&$algo, partition).into_ids() +impl Algorithm for coupe::KarmarkarKarp { + fn run(&mut self, partition: &mut [usize], problem: &Problem) -> Result<()> { + use weight::Array::*; + match &problem.weights { + Integers(is) => { + let weights = is.iter().map(|weight| weight[0]); + self.partition(partition, weights)?; } - 3 => { - let points = coupe::PointsView::to_points_nd($points.as_slice()); - let partition = - coupe::partition::Partition::from_ids(points, $weights, $partition.to_vec()); - coupe::PartitionImprover::<3>::improve_partition(&$algo, partition).into_ids() + Floats(fs) => { + let weights = fs.iter().map(|weight| coupe::Real::from(weight[0])); + self.partition(partition, weights)?; } - _ => anyhow::bail!("only 2D and 3D meshes are supported"), - }; - $partition.copy_from_slice(&ids); - }}; - ( $algo:expr, $partition:expr, $problem:expr ) => {{ - let weights: Vec = match &$problem.weights { - weight::Array::Floats(ws) => ws.iter().map(|weight| weight[0]).collect(), - weight::Array::Integers(_) => anyhow::bail!("integers are not supported by coupe yet"), - }; - coupe_part_improve!( - $algo, - $partition, - $problem.dimension, - $problem.points, - &weights - ) - }}; + } + Ok(()) + } } -struct Problem { - dimension: usize, - points: Vec, - weights: weight::Array, +impl Algorithm for coupe::VnBest { + fn run(&mut self, partition: &mut [usize], problem: &Problem) -> Result<()> { + use weight::Array::*; + match &problem.weights { + Integers(is) => { + let weights = is.iter().map(|weight| weight[0]); + self.partition(partition, weights)?; + } + Floats(fs) => { + let weights = fs.iter().map(|weight| coupe::Real::from(weight[0])); + self.partition(partition, weights)?; + } + } + Ok(()) + } } -type Algorithm = Box Result>; +impl Algorithm for coupe::Rcb { + fn run(&mut self, partition: &mut [usize], problem: &Problem) -> Result<()> { + use weight::Array::*; + let points = problem.points.par_iter().cloned(); + match &problem.weights { + Integers(is) => { + let weights = is.par_iter().map(|weight| weight[0]); + self.partition(partition, (points, weights))?; + } + Floats(fs) => { + let weights = fs.par_iter().map(|weight| weight[0]); + self.partition(partition, (points, weights))?; + } + } + Ok(()) + } +} + +impl Algorithm for coupe::HilbertCurve { + fn run(&mut self, partition: &mut [usize], problem: &Problem) -> Result<()> { + use weight::Array::*; + if D != 2 { + anyhow::bail!("hilbert is only implemented for 2D meshes"); + } + // SAFETY: is a noop since D == 2 + let points = + unsafe { mem::transmute::<&Vec>, &Vec>>(&problem.points) }; + match &problem.weights { + Integers(_) => anyhow::bail!("hilbert is only implemented for floats"), + Floats(fs) => { + let weights: Vec = fs.iter().map(|weight| weight[0]).collect(); + self.partition(partition, (points, weights))?; + } + } + Ok(()) + } +} -fn parse_algorithm(spec: String) -> Result { +fn parse_algorithm(spec: &str) -> Result>> { let mut args = spec.split(','); - let name = args.next().context("empty algorithm spec")?; + let name = args.next().context("it's empty")?; fn optional(maybe_arg: Option>, default: T) -> Result { Ok(maybe_arg.transpose()?.unwrap_or(default)) @@ -100,77 +150,75 @@ fn parse_algorithm(spec: String) -> Result { bytes.resize(32_usize, 0_u8); bytes.try_into().unwrap() }; - let mut rng = rand_pcg::Pcg64::from_seed(seed); - Box::new(move |partition, problem| { - let algo = coupe::Random::new(&mut rng, part_count); - let weights = vec![0.0; problem.points.len()]; - let ids = coupe_part!(algo, problem.dimension, problem.points, &weights); - partition.copy_from_slice(&ids); - Ok(RunInfo::default()) - }) - } - "greedy" => { - let part_count = required(usize_arg(args.next()))?; - Box::new(move |partition, problem| { - let ids = coupe_part!(coupe::Greedy::new(part_count), problem); - partition.copy_from_slice(&ids); - Ok(RunInfo::default()) - }) - } - "kk" => { - let part_count = required(usize_arg(args.next()))?; - Box::new(move |partition, problem| { - let ids = coupe_part!(coupe::KarmarkarKarp::new(part_count), problem); - partition.copy_from_slice(&ids); - Ok(RunInfo::default()) - }) - } - "vn-best" => { - let part_count = required(usize_arg(args.next()))?; - Box::new(move |partition, problem| { - coupe_part_improve!(coupe::VnBest::new(part_count), partition, problem); - Ok(RunInfo::default()) - }) - } - "rcb" => { - let iter_count = required(usize_arg(args.next()))?; - Box::new(move |partition, problem| { - let ids = coupe_part!(coupe::Rcb::new(iter_count), problem); - partition.copy_from_slice(&ids); - Ok(RunInfo::default()) - }) - } - "hilbert" => { - let num_partitions = required(usize_arg(args.next()))?; - let order = optional(usize_arg(args.next()), 12)? as u32; - Box::new(move |partition, problem| { - let algo = coupe::HilbertCurve { - num_partitions, - order, - }; - let weights: Vec = match &problem.weights { - weight::Array::Floats(ws) => ws.iter().map(|weight| weight[0]).collect(), - weight::Array::Integers(_) => { - anyhow::bail!("hilbert is not implemented for integers") - } - }; - let res = match problem.dimension { - 2 => coupe::Partitioner::<2>::partition( - &algo, - problem.points.as_slice(), - &weights, - ) - .into_ids(), - _ => anyhow::bail!("hilbert is only implemented in 2D"), - }; - partition.copy_from_slice(&res); - Ok(RunInfo::default()) - }) + let rng = rand_pcg::Pcg64::from_seed(seed); + Box::new(coupe::Random { rng, part_count }) } + "greedy" => Box::new(coupe::Greedy { + part_count: required(usize_arg(args.next()))?, + }), + "kk" => Box::new(coupe::KarmarkarKarp { + part_count: required(usize_arg(args.next()))?, + }), + "vn-best" => Box::new(coupe::VnBest { + part_count: required(usize_arg(args.next()))?, + }), + "rcb" => Box::new(coupe::Rcb { + iter_count: required(usize_arg(args.next()))?, + }), + "hilbert" => Box::new(coupe::HilbertCurve { + part_count: required(usize_arg(args.next()))?, + order: optional(usize_arg(args.next()), 12)? as u32, + }), _ => anyhow::bail!("unknown algorithm {:?}", name), }) } +fn main_d( + matches: getopts::Matches, + mesh: Mesh, + weights: weight::Array, +) -> Result> { + let points: Vec<_> = mesh + .elements() + .filter_map(|(element_type, nodes, _element_ref)| { + if element_type.dimension() != mesh.dimension() { + return None; + } + let mut barycentre = [0.0; D]; + for node_idx in nodes { + let node_coordinates = mesh.node(*node_idx); + for (bc_coord, node_coord) in barycentre.iter_mut().zip(node_coordinates) { + *bc_coord += node_coord; + } + } + for bc_coord in &mut barycentre { + *bc_coord /= nodes.len() as f64; + } + Some(PointND::from(barycentre)) + }) + .collect(); + + let mut partition = vec![0; points.len()]; + let problem = Problem { points, weights }; + + let algorithms: Vec<_> = matches + .opt_strs("a") + .into_iter() + .map(|algorithm_spec| { + parse_algorithm(&algorithm_spec) + .with_context(|| format!("invalid algorithm {:?}", algorithm_spec)) + }) + .collect::>()?; + + for mut algorithm in algorithms { + algorithm + .run(&mut partition, &problem) + .context("failed to apply algorithm")?; + } + + Ok(partition) +} + fn main() -> Result<()> { let mut options = getopts::Options::new(); options.optflag("h", "help", "print this help menu"); @@ -191,37 +239,10 @@ fn main() -> Result<()> { return Ok(()); } - let algorithms: Vec<_> = matches - .opt_strs("a") - .into_iter() - .map(parse_algorithm) - .collect::>()?; - let mesh_file = matches .opt_str("m") .context("missing required option 'mesh'")?; - let mesh = mesh_io::medit::Mesh::from_file(mesh_file).context("failed to read mesh file")?; - - let points: Vec<_> = mesh - .elements() - .filter_map(|(element_type, nodes, _element_ref)| { - if element_type.dimension() != mesh.dimension() { - return None; - } - let mut barycentre = vec![0.0; mesh.dimension()]; - for node_idx in nodes { - let node_coordinates = mesh.node(*node_idx); - for (bc_coord, node_coord) in barycentre.iter_mut().zip(node_coordinates) { - *bc_coord += node_coord; - } - } - for bc_coord in &mut barycentre { - *bc_coord /= nodes.len() as f64; - } - Some(barycentre) - }) - .flatten() - .collect(); + let mesh = Mesh::from_file(mesh_file).context("failed to read mesh file")?; let weight_file = matches .opt_str("w") @@ -230,16 +251,11 @@ fn main() -> Result<()> { let weight_file = io::BufReader::new(weight_file); let weights = weight::read(weight_file).context("failed to read weight file")?; - let problem = Problem { - dimension: mesh.dimension(), - points, - weights, + let partition = match mesh.dimension() { + 2 => main_d::<2>(matches, mesh, weights)?, + 3 => main_d::<3>(matches, mesh, weights)?, + n => anyhow::bail!("expected 2D or 3D mesh, got a {n}D mesh"), }; - let mut partition = vec![0; problem.points.len() / problem.dimension]; - - for mut algorithm in algorithms { - algorithm(&mut partition, &problem).context("failed to apply algorithm")?; - } let stdout = io::stdout(); let stdout = stdout.lock(); From 586f1d4a0c71631a352905b234e0c0fc2dd832f4 Mon Sep 17 00:00:00 2001 From: Hubert Hirtz Date: Wed, 23 Mar 2022 13:58:35 +0100 Subject: [PATCH 2/6] Hide geometry module All of its items were exported anyway --- benches/benchmark.rs | 2 +- benches/generator.rs | 2 +- src/lib.rs | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/benches/benchmark.rs b/benches/benchmark.rs index e57c19da..09253483 100644 --- a/benches/benchmark.rs +++ b/benches/benchmark.rs @@ -1,7 +1,7 @@ mod generator; -use coupe::geometry::Point2D; use coupe::Partition as _; +use coupe::Point2D; use criterion::{criterion_group, criterion_main}; use criterion::{Criterion, Throughput}; use rayon::prelude::*; diff --git a/benches/generator.rs b/benches/generator.rs index 34ed9987..2191fafe 100644 --- a/benches/generator.rs +++ b/benches/generator.rs @@ -1,4 +1,4 @@ -use coupe::geometry::Point2D; +use coupe::Point2D; use rand::{self, Rng}; pub fn uniform_f64(min: f64, max: f64, num_vals: usize) -> Vec { diff --git a/src/lib.rs b/src/lib.rs index 8e3e2b66..c0470515 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -29,7 +29,7 @@ //! - [Kernighan-Lin][KernighanLin] mod algorithms; -pub mod geometry; +mod geometry; pub mod imbalance; mod real; pub mod topology; From 90033fc26833f2867d53d0816d8348579d9500dd Mon Sep 17 00:00:00 2001 From: Hubert Hirtz Date: Wed, 23 Mar 2022 14:32:01 +0100 Subject: [PATCH 3/6] hilbert: replace panics with errors --- src/algorithms.rs | 1 + src/algorithms/hilbert_curve.rs | 63 +++++++++++++++++++++------------ 2 files changed, 41 insertions(+), 23 deletions(-) diff --git a/src/algorithms.rs b/src/algorithms.rs index 3a7c470b..906328cf 100644 --- a/src/algorithms.rs +++ b/src/algorithms.rs @@ -17,6 +17,7 @@ pub use ckk::CompleteKarmarkarKarp; pub use fiduccia_mattheyses::FiducciaMattheyses; pub use graph_growth::GraphGrowth; pub use greedy::Greedy; +pub use hilbert_curve::Error as HilbertCurveError; pub use hilbert_curve::HilbertCurve; pub use k_means::KMeans; pub use kernighan_lin::KernighanLin; diff --git a/src/algorithms/hilbert_curve.rs b/src/algorithms/hilbert_curve.rs index 5f91d359..b7d730fc 100644 --- a/src/algorithms/hilbert_curve.rs +++ b/src/algorithms/hilbert_curve.rs @@ -12,6 +12,7 @@ use crate::geometry::{Mbr, Point2D}; use rayon::prelude::*; +use std::fmt; fn hilbert_curve_partition( partition: &mut [usize], @@ -20,10 +21,7 @@ fn hilbert_curve_partition( part_count: usize, order: usize, ) { - assert!( - order < 32, - "Cannot construct a Hilbert curve of order >= 32 because 2^32 would overflow u32 capacity." - ); + debug_assert!(order < 32); let compute_hilbert_index = hilbert_index_computer(points, order); let mut permutation = (0..points.len()).into_par_iter().collect::>(); @@ -72,15 +70,8 @@ pub(crate) fn hilbert_curve_reorder(points: &[Point2D], order: usize) -> Vec= 32 because 2^32 would overflow u32 capacity." - ); +fn hilbert_curve_reorder_permu(points: &[Point2D], permutation: &mut [usize], order: usize) { + debug_assert!(order < 32); let compute_hilbert_index = hilbert_index_computer(points, order); @@ -113,17 +104,14 @@ fn hilbert_index_computer(points: &[Point2D], order: usize) -> impl Fn((f64, f64 } fn encode(x: u32, y: u32, order: usize) -> u32 { - assert!( - order < 32, - "Cannot construct a Hilbert curve of order >= 32 because 2^32 would overflow u32 capacity." - ); - assert!( + debug_assert!(order < 32); + debug_assert!( x < 2u32.pow(order as u32), "Cannot encode the point {:?} on an hilbert curve of order {} because x >= 2^order.", (x, y), order, ); - assert!( + debug_assert!( y < 2u32.pow(order as u32), "Cannot encode the point {:?} on an hilbert curve of order {} because y >= 2^order.", (x, y), @@ -174,13 +162,13 @@ fn interleave_bits(odd: u32, even: u32) -> u32 { // Compute a mapping from [a_min; a_max] to [b_min; b_max] fn segment_to_segment(a_min: f64, a_max: f64, b_min: f64, b_max: f64) -> impl Fn(f64) -> f64 { - assert!( + debug_assert!( a_min <= a_max, "Cannot construct a segment to segment mapping because a_max < a_min. a_min = {}, a_max = {}.", a_min, a_max, ); - assert!( + debug_assert!( b_min <= b_max, "Cannot construct a segment to segment mapping because b_max < b_min. b_min = {}, b_max = {}.", b_min, @@ -193,7 +181,7 @@ fn segment_to_segment(a_min: f64, a_max: f64, b_min: f64, b_max: f64) -> impl Fn let beta = b_min - a_min * alpha; move |x| { - assert!( + debug_assert!( a_min <= x && x <= a_max, "Called a mapping from [{}, {}] to [{}, {}] with the invalid value {}.", a_min, @@ -206,6 +194,29 @@ fn segment_to_segment(a_min: f64, a_max: f64, b_min: f64, b_max: f64) -> impl Fn } } +#[derive(Debug)] +#[non_exhaustive] +pub enum Error { + /// Invalid space filling curve order. + InvalidOrder { max: u32, actual: u32 }, +} + +impl fmt::Display for Error { + fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { + match self { + Error::InvalidOrder { max, actual } => { + write!( + f, + "given hilbert curve order too high. Got {}, max={}.", + actual, max + ) + } + } + } +} + +impl std::error::Error for Error {} + /// # Hilbert space-filling curve algorithm /// /// An implementation of the Hilbert curve based on @@ -255,13 +266,19 @@ where W: AsRef<[f64]>, { type Metadata = (); - type Error = std::convert::Infallible; + type Error = Error; fn partition( &mut self, part_ids: &mut [usize], (points, weights): (P, W), ) -> Result { + if self.order >= 32 { + return Err(Error::InvalidOrder { + max: 31, + actual: self.order, + }); + } hilbert_curve_partition( part_ids, points.as_ref(), From cca044ab9ec5f48df1c6ba62e1b0f5820d360d1a Mon Sep 17 00:00:00 2001 From: Hubert Hirtz Date: Wed, 23 Mar 2022 14:32:20 +0100 Subject: [PATCH 4/6] greedy: replace panics with errors --- src/algorithms/greedy.rs | 30 ++++++++++++++++++++---------- 1 file changed, 20 insertions(+), 10 deletions(-) diff --git a/src/algorithms/greedy.rs b/src/algorithms/greedy.rs index e58fdfa4..db35643b 100644 --- a/src/algorithms/greedy.rs +++ b/src/algorithms/greedy.rs @@ -1,3 +1,4 @@ +use super::Error; use num::Zero; use std::ops::AddAssign; @@ -5,10 +6,11 @@ use std::ops::AddAssign; fn greedy( partition: &mut [usize], weights: impl IntoIterator, - num_parts: usize, -) { - if num_parts < 2 { - return; + part_count: usize, +) -> Result<(), Error> { + if part_count < 2 { + partition.fill(0); + return Ok(()); } // Initialization: make the partition and record the weight of each part in another vector. @@ -16,9 +18,16 @@ fn greedy( .into_iter() .zip(0..) // Keep track of the weights' indicies .collect(); - assert_eq!(partition.len(), weights.len()); + + if weights.len() != partition.len() { + return Err(Error::InputLenMismatch { + expected: partition.len(), + actual: weights.len(), + }); + } + weights.sort_unstable(); - let mut part_weights = vec![T::zero(); num_parts]; + let mut part_weights = vec![T::zero(); part_count]; // Put each weight in the lightweightest part. for (weight, weight_id) in weights.into_iter().rev() { @@ -26,10 +35,12 @@ fn greedy( .iter() .enumerate() .min_by_key(|(_idx, part_weight)| *part_weight) - .unwrap(); + .unwrap(); // Will not panic because !part_weights.is_empty() partition[weight_id] = min_part_weight_idx; part_weights[min_part_weight_idx] += weight; } + + Ok(()) } /// Greedily assign weights to each part. @@ -57,14 +68,13 @@ where W::Item: Ord + Zero + Clone + AddAssign, { type Metadata = (); - type Error = std::convert::Infallible; + type Error = Error; fn partition( &mut self, part_ids: &mut [usize], weights: W, ) -> Result { - greedy(part_ids, weights, self.part_count); - Ok(()) + greedy(part_ids, weights, self.part_count) } } From 0b8bf5c534694c91f152c461d3e5fe61f6ca33df Mon Sep 17 00:00:00 2001 From: Hubert Hirtz Date: Wed, 23 Mar 2022 14:32:32 +0100 Subject: [PATCH 5/6] rcb: replace panics with errors --- src/algorithms/recursive_bisection.rs | 59 ++++++++++++++++++--------- 1 file changed, 39 insertions(+), 20 deletions(-) diff --git a/src/algorithms/recursive_bisection.rs b/src/algorithms/recursive_bisection.rs index 326af565..eaf1acaa 100644 --- a/src/algorithms/recursive_bisection.rs +++ b/src/algorithms/recursive_bisection.rs @@ -1,3 +1,4 @@ +use super::Error; use crate::geometry::Mbr; use crate::geometry::PointND; use async_lock::Mutex; @@ -162,7 +163,7 @@ where .map(|item| item.point[coord]) .minmax() .into_option() - .unwrap(); + .unwrap(); // Won't panic because items has at least two elements. mem::drop(enter); @@ -501,7 +502,12 @@ fn rcb_thread( futures_lite::future::block_on(task); } -fn rcb(partition: &mut [usize], points: P, weights: W, iter_count: usize) +fn rcb( + partition: &mut [usize], + points: P, + weights: W, + iter_count: usize, +) -> Result<(), Error> where P: rayon::iter::IntoParallelIterator>, P::Iter: rayon::iter::IndexedParallelIterator, @@ -514,8 +520,18 @@ where let points = points.into_par_iter(); let weights = weights.into_par_iter(); - assert_eq!(points.len(), weights.len()); - assert_eq!(points.len(), partition.len()); + if weights.len() != partition.len() { + return Err(Error::InputLenMismatch { + expected: partition.len(), + actual: weights.len(), + }); + } + if points.len() != partition.len() { + return Err(Error::InputLenMismatch { + expected: partition.len(), + actual: points.len(), + }); + } let init_span = tracing::info_span!("convert input and make initial data structures"); let enter = init_span.enter(); @@ -546,6 +562,8 @@ where s.spawn(move |_| rcb_thread(iteration_ctxs, chunk, iter_count, 0.05)); } }); + + Ok(()) } /// # Recursive Coordinate Bisection algorithm @@ -602,15 +620,18 @@ where W::Iter: rayon::iter::IndexedParallelIterator, { type Metadata = (); - type Error = std::convert::Infallible; + type Error = Error; fn partition( &mut self, part_ids: &mut [usize], (points, weights): (P, W), ) -> Result { - rcb(part_ids, points, weights, self.iter_count); - Ok(()) + if part_ids.len() < 2 { + // Would make Itertools::minmax().into_option() return None. + return Ok(()); + } + rcb(part_ids, points, weights, self.iter_count) } } @@ -647,7 +668,12 @@ 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. -fn rib(partition: &mut [usize], points: &[PointND], weights: W, n_iter: usize) +fn rib( + partition: &mut [usize], + points: &[PointND], + weights: W, + n_iter: usize, +) -> Result<(), Error> where Const: DimSub>, DefaultAllocator: Allocator, Const, Buffer = ArrayStorage> @@ -658,15 +684,8 @@ where W::Item: num::ToPrimitive, W::Iter: rayon::iter::IndexedParallelIterator, { - let weights = weights.into_par_iter(); - - assert_eq!(points.len(), weights.len()); - assert_eq!(points.len(), partition.len()); - let mbr = Mbr::from_points(points); - let points = points.par_iter().map(|p| mbr.mbr_to_aabb(p)); - // When the rotation is done, we just apply RCB rcb(partition, points, weights, n_iter) } @@ -727,15 +746,14 @@ where W::Iter: rayon::iter::IndexedParallelIterator, { type Metadata = (); - type Error = std::convert::Infallible; + type Error = Error; fn partition( &mut self, part_ids: &mut [usize], (points, weights): (&'a [PointND], W), ) -> Result { - rib(part_ids, points, weights, self.iter_count); - Ok(()) + rib(part_ids, points, weights, self.iter_count) } } @@ -795,7 +813,8 @@ mod tests { .num_threads(1) // make the test deterministic .build() .unwrap() - .install(|| rcb(&mut partition, points, weights, 2)); + .install(|| rcb(&mut partition, points, weights, 2)) + .unwrap(); assert_eq!(partition[0], partition[6]); assert_eq!(partition[1], partition[7]); @@ -825,7 +844,7 @@ mod tests { let weights: Vec = (0..points.len()).map(|_| rand::random()).collect(); let mut partition = vec![0; points.len()]; - rcb(&mut partition, points, weights.par_iter().cloned(), 3); + rcb(&mut partition, points, weights.par_iter().cloned(), 3).unwrap(); let mut loads: HashMap = HashMap::new(); let mut sizes: HashMap = HashMap::new(); From e6b14ec10ccb6a861706e275e36a5fc0ff68ae34 Mon Sep 17 00:00:00 2001 From: Hubert Hirtz Date: Wed, 23 Mar 2022 14:32:48 +0100 Subject: [PATCH 6/6] kk: replace panics with errors --- src/algorithms/kk.rs | 52 +++++++++++++++++++++----------------------- 1 file changed, 25 insertions(+), 27 deletions(-) diff --git a/src/algorithms/kk.rs b/src/algorithms/kk.rs index 944fe701..60b970fc 100644 --- a/src/algorithms/kk.rs +++ b/src/algorithms/kk.rs @@ -1,3 +1,4 @@ +use super::Error; use std::collections::BinaryHeap; use std::ops::Sub; use std::ops::SubAssign; @@ -9,7 +10,7 @@ use num::Zero; /// # Differences with the k-partitioning implementation /// /// This function has better performance than [kk] called with `num_parts == 2`. -fn kk_bipart(partition: &mut [usize], weights: impl IntoIterator) +fn kk_bipart(partition: &mut [usize], weights: impl Iterator) where T: Ord + Sub, { @@ -17,10 +18,6 @@ where .into_iter() .zip(0..) // Keep track of the weights' indicies .collect(); - assert_eq!(partition.len(), weights.len()); - if weights.is_empty() { - return; - } // Core algorithm: find the imbalance of the partition. // "opposites" is built in this loop to backtrack the solution. It tracks weights that must end @@ -54,30 +51,15 @@ where fn kk(partition: &mut [usize], weights: I, num_parts: usize) where T: Zero + Ord + Sub + SubAssign + Copy, - I: IntoIterator, - ::IntoIter: ExactSizeIterator, + I: Iterator + ExactSizeIterator, { - let weights = weights.into_iter(); - let num_weights = weights.len(); - assert_eq!(partition.len(), num_weights); - if num_weights == 0 { - return; - } - if num_parts < 2 { - return; - } - if num_parts == 2 { - // The bi-partitioning is a special case that can be handled faster than - // the general case. - return kk_bipart(partition, weights); - } - // Initialize "m", a "k*num_weights" matrix whose first column is "weights". + let weight_count = weights.len(); let mut m: BinaryHeap> = weights .zip(0..) .map(|(w, id)| { let mut v: Vec<(T, usize)> = (0..num_parts) - .map(|p| (T::zero(), num_weights * p + id)) + .map(|p| (T::zero(), weight_count * p + id)) .collect(); v[0].0 = w; v @@ -88,7 +70,7 @@ where // largest weights in two different parts, the largest weight of each row is put into the same // part as the smallest one, and so on. - let mut opposites = Vec::with_capacity(num_weights); + let mut opposites = Vec::with_capacity(weight_count); while 2 <= m.len() { let a = m.pop().unwrap(); let b = m.pop().unwrap(); @@ -119,7 +101,7 @@ where // Backtracking. Same as the bi-partitioning case. // parts = [ [m0i] for m0i in m[0] ] - let mut parts: Vec = vec![0; num_parts * num_weights]; + let mut parts: Vec = vec![0; num_parts * weight_count]; let imbalance = m.pop().unwrap(); // first and last element of "m". for (i, w) in imbalance.into_iter().enumerate() { // Put each remaining element in a different part. @@ -160,14 +142,30 @@ where W::Item: Zero + Ord + Sub + SubAssign + Copy, { type Metadata = (); - type Error = std::convert::Infallible; + type Error = Error; fn partition( &mut self, part_ids: &mut [usize], weights: W, ) -> Result { - kk(part_ids, weights, self.part_count); + if self.part_count < 2 || part_ids.len() < 2 { + return Ok(()); + } + let weights = weights.into_iter(); + if weights.len() != part_ids.len() { + return Err(Error::InputLenMismatch { + expected: part_ids.len(), + actual: weights.len(), + }); + } + if self.part_count == 2 { + // The bi-partitioning is a special case that can be handled faster + // than the general case. + kk_bipart(part_ids, weights); + } else { + kk(part_ids, weights, self.part_count); + } Ok(()) } }