diff --git a/src/cartesian/mod.rs b/src/cartesian/mod.rs index 9b50f33f..d512465a 100644 --- a/src/cartesian/mod.rs +++ b/src/cartesian/mod.rs @@ -6,6 +6,10 @@ use rayon::iter::IndexedParallelIterator; use rayon::iter::IntoParallelRefIterator; use rayon::iter::IntoParallelRefMutIterator; use rayon::iter::ParallelIterator; +use std::collections::BTreeMap; +use std::collections::HashMap; +use std::collections::VecDeque; +use std::fmt; use std::iter::Sum; use std::marker::PhantomData; use std::num::NonZeroUsize; @@ -131,6 +135,7 @@ impl Grid<2> { where W: Send + Sync + PartialOrd + Num + Sum + AsPrimitive, f64: AsPrimitive, + W: num_traits::NumAssign, { let total_weight: W = weights.par_iter().cloned().sum(); let iters = rcb::recurse_2d( @@ -145,6 +150,226 @@ impl Grid<2> { let pos = self.position_of(i); *p = iters.part_of(pos, 1); }); + + println!("{}", iters.fmt_svg(self, 1)); + let mut segs = iters.segments_2d(self, 1); + + loop { + #[derive(Copy, Clone, Debug)] + enum Direction { + Lower, + Higher, + } + + #[derive(Copy, Clone, Debug)] + enum Axis { + X = 0, + Y = 1, + } + + #[derive(Debug)] + struct Move { + orientation: Axis, + gain: f64, + seg: Segment, + direction: Direction, + } + + eprintln!("\nNew pass"); + let mut best_lambda_move = Move { + orientation: Axis::X, + gain: 0.0, + seg: segs.c[0][0], + direction: Direction::Lower, + }; + + let lambda = |partition: &[usize], pos: [usize; 2]| -> W { + let i = self.index_of(pos); + let w = weights[i]; + let p = partition[i]; + self.neighbors(i) + .filter(|(n, _): &(usize, i32)| partition[*n] != p) + .map(|_| w) + .sum() + }; + let seg_lambda = |partition: &[usize], seg: Segment, orientation: Axis| -> W { + (seg.at..seg.at + 2) + .flat_map(|at| { + (seg.start..seg.end).map(move |se| match orientation { + Axis::X => lambda(partition, [se, at]), + Axis::Y => lambda(partition, [at, se]), + }) + }) + .sum() + }; + let check_move = |partition: &mut [usize], + seg: Segment, + orientation: Axis, + direction: Direction| + -> f64 { + let src_at; // position of the parts that will expand + let dst_at; // position of the parts that will shrink + let lambda_check_at; // start of the 3-wide region to check for lambda + match direction { + Direction::Lower => { + src_at = seg.at; + dst_at = seg.at - 1; + lambda_check_at = seg.at - 2; + } + Direction::Higher => { + src_at = seg.at - 1; + dst_at = seg.at; + lambda_check_at = seg.at - 1; + } + } + let src_cells = (seg.start..seg.end).map(|a| match orientation { + Axis::X => [a, src_at], + Axis::Y => [src_at, a], + }); + let dst_cells = (seg.start..seg.end).map(|a| match orientation { + Axis::X => [a, dst_at], + Axis::Y => [dst_at, a], + }); + let prev = seg_lambda( + partition, + Segment { + at: lambda_check_at, + ..seg + }, + orientation, + ); + let old_cells: Vec = dst_cells + .clone() + .map(|p| partition[self.index_of(p)]) + .collect(); + for (src, dst) in src_cells.zip(dst_cells.clone()) { + let src = partition[self.index_of(src)]; + let dst = &mut partition[self.index_of(dst)]; + debug_assert_ne!(*dst, src); + *dst = src; + } + let next = seg_lambda( + partition, + Segment { + at: lambda_check_at, + ..seg + }, + orientation, + ); + for (dst, old) in dst_cells.zip(old_cells) { + partition[self.index_of(dst)] = old; + } + prev.as_() - next.as_() + }; + + // Testing horizontal segments for lambda cut. + for seg in segs.moves(0) { + if seg.at >= 2 { + // Test if we can move segment down. + let gain = check_move(partition, seg, Axis::X, Direction::Lower); + if gain > best_lambda_move.gain { + best_lambda_move = Move { + orientation: Axis::X, + gain, + seg, + direction: Direction::Lower, + }; + } + } + if seg.at + 2 >= usize::from(self.size[1]) { + // Test if we can move segment up. + let gain = check_move(partition, seg, Axis::X, Direction::Higher); + if gain > best_lambda_move.gain { + best_lambda_move = Move { + orientation: Axis::X, + gain, + seg, + direction: Direction::Higher, + }; + } + } + } + + // Testing vertical segments for lambda cut. + for seg in segs.moves(1) { + if seg.at >= 2 { + // Test if we can move segment down. + let gain = check_move(partition, seg, Axis::Y, Direction::Lower); + if gain > best_lambda_move.gain { + best_lambda_move = Move { + orientation: Axis::Y, + gain, + seg, + direction: Direction::Lower, + }; + } + } + if seg.at + 2 >= usize::from(self.size[0]) { + // Test if we can move segment up. + let gain = check_move(partition, seg, Axis::Y, Direction::Higher); + if gain > best_lambda_move.gain { + best_lambda_move = Move { + orientation: Axis::Y, + gain, + seg, + direction: Direction::Higher, + }; + } + } + } + + let Move { + orientation, + gain, + seg, + direction, + } = best_lambda_move; + + if gain == 0.0 { + eprintln!("no gain, no pain"); + break; + } + + eprintln!(" Best lambda move: {best_lambda_move:?}"); + + let src_at; // position of the parts that will expand + let dst_at; // position of the parts that will shrink + let new_seg_at; // new position of the segment + match direction { + Direction::Lower => { + src_at = seg.at; + dst_at = seg.at - 1; + new_seg_at = seg.at - 1; + } + Direction::Higher => { + src_at = seg.at - 1; + dst_at = seg.at; + new_seg_at = seg.at + 1; + } + } + let src_cells = (seg.start..seg.end).map(|a| match orientation { + Axis::X => [a, src_at], + Axis::Y => [src_at, a], + }); + let dst_cells = (seg.start..seg.end).map(|a| match orientation { + Axis::X => [a, dst_at], + Axis::Y => [dst_at, a], + }); + for (src, dst) in src_cells.zip(dst_cells) { + partition[self.index_of(dst)] = partition[self.index_of(src)]; + } + for comp in segs.components(orientation as usize, &seg) { + comp.at = new_seg_at; + } + for ortho in &mut segs.c[1 - (orientation as usize)] { + // Update orthogonal segments. + if ortho.start == seg.at { + ortho.start = new_seg_at; + } else if ortho.end == seg.at { + ortho.end = new_seg_at; + } + } + } } } @@ -281,6 +506,330 @@ where } } +fn transpose(p: [T; 2]) -> [T; 2] { + let [a, b] = p; + [b, a] +} + +#[derive(Debug)] +pub enum SplitTree { + Whole, + Split { + position: usize, + left: Box, + right: Box, + }, +} + +impl SplitTree { + fn part_of(&self, pos: [usize; D], mut start_coord: usize) -> usize { + let mut it = self; + let mut part_id = 0; + while let Self::Split { + position, + left, + right, + } = it + { + if pos[start_coord] < *position { + part_id *= 2; + it = left; + } else { + part_id = 2 * part_id + 1; + it = right; + } + start_coord = (start_coord + 1) % D; + } + part_id + } + + pub fn fmt_svg(&self, grid: Grid<2>, start_coord: usize) -> impl fmt::Display + '_ { + struct ShowSvg<'a> { + tree: &'a SplitTree, + grid: Grid<2>, + start_coord: usize, + } + + fn print_splits( + f: &mut fmt::Formatter<'_>, + sg: SubGrid<2>, + tree: &SplitTree, + coord: usize, + iter: usize, + ) -> fmt::Result { + let SplitTree::Split { position, left, right } = tree + else { return Ok(()) }; + + // Recurse before so that lines from first iterations are shown + // above lines from the next ones. + let (sg_left, sg_right) = sg.split_at(coord, *position); + print_splits(f, sg_left, left, (coord + 1) % 2, iter + 1)?; + print_splits(f, sg_right, right, (coord + 1) % 2, iter + 1)?; + + let Range { start, end } = sg.axis(1 - coord); + let mut p1 = [*position, start]; + let mut p2 = [*position, end]; + if coord == 1 { + p1 = transpose(p1); + p2 = transpose(p2); + } + let color = match iter % 10 { + 0 => "maroon", + 1 => "green", + 2 => "red", + 3 => "lime", + 4 => "purple", + 5 => "olive", + 6 => "fuchsia", + 7 => "yellow", + 8 => "navy", + _ => "blue", + }; + writeln!( + f, + r#""#, + p1[0], p1[1], p2[0], p2[1], color, + ) + } + + impl fmt::Display for ShowSvg<'_> { + fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { + let [gwidth, gheight] = self.grid.size; + let sg = self.grid.into_subgrid(); + + writeln!( + f, + r#""# + )?; + + print_splits(f, sg, self.tree, self.start_coord, 0)?; + + for [x, y] in self.tree.joints_2d(self.grid, self.start_coord) { + writeln!(f, r#""#)?; + } + + writeln!(f, "") + } + } + + ShowSvg { + tree: self, + grid, + start_coord, + } + } + + pub fn joints_2d(&self, grid: Grid<2>, start_coord: usize) -> impl Iterator { + fn aux( + joints: &mut HashMap<[usize; 2], usize>, + tree: &SplitTree, + sg: SubGrid<2>, + coord: usize, + ) { + let SplitTree::Split { position, left, right } = tree + else { return }; + + let Range { start, end } = sg.axis(1 - coord); + let mut p1 = [*position, start]; + let mut p2 = [*position, end]; + if coord == 1 { + p1 = transpose(p1); + p2 = transpose(p2); + } + + *joints.entry(p1).or_default() += 1; + *joints.entry(p2).or_default() += 1; + + let (sg_left, sg_right) = sg.split_at(coord, *position); + aux(joints, left, sg_left, (coord + 1) % 2); + aux(joints, right, sg_right, (coord + 1) % 2); + } + + let mut joints = HashMap::new(); + + aux(&mut joints, self, grid.into_subgrid(), start_coord); + joints.retain(|_, occ| *occ > 1); + + joints.into_keys() + } + + pub fn segments_2d(&self, grid: Grid<2>, start_coord: usize) -> Segments<2> { + fn aux(c: &mut [Vec; 2], tree: &SplitTree, sg: SubGrid<2>, coord: usize) { + let SplitTree::Split { position, left, right } = tree + else { return }; + + let Range { start, end } = sg.axis(1 - coord); + c[1 - coord].push(Segment { + start, + end, + at: *position, + }); + + let (sg_left, sg_right) = sg.split_at(coord, *position); + aux(c, left, sg_left, (coord + 1) % 2); + aux(c, right, sg_right, (coord + 1) % 2); + } + + let mut c = [Vec::new(), Vec::new()]; + aux(&mut c, self, grid.into_subgrid(), start_coord); + + for [x, y] in self.joints_2d(grid, start_coord) { + for i in 0..c[0].len() { + let seg = &c[0][i]; + if y != seg.at { + continue; + } + let Some((left, right)) = seg.split_at(x) + else { continue }; + c[0][i] = left; + c[0].push(right); + } + for i in 0..c[1].len() { + let seg = &c[1][i]; + if x != seg.at { + continue; + } + let Some((left, right)) = seg.split_at(y) + else { continue }; + c[1][i] = left; + c[1].push(right); + } + } + + Segments { c, grid } + } +} + +#[derive(Clone, Copy, Debug, PartialEq)] +pub struct Segment { + start: usize, + end: usize, + at: usize, +} + +impl Segment { + pub fn split_at(&self, pos: usize) -> Option<(Segment, Segment)> { + if pos <= self.start || self.end <= pos { + return None; + } + Some(( + Segment { + start: self.start, + end: pos, + at: self.at, + }, + Segment { + start: pos, + end: self.end, + at: self.at, + }, + )) + } + + pub fn p_start_2d(&self, coord: usize) -> [usize; 2] { + if coord == 0 { + [self.at, self.start] + } else { + [self.start, self.at] + } + } + + pub fn p_end_2d(&self, coord: usize) -> [usize; 2] { + if coord == 0 { + [self.at, self.end] + } else { + [self.end, self.at] + } + } + + pub fn contains(&self, other: &Self) -> bool { + self.at == other.at && self.start <= other.start && other.end <= self.end + } +} + +#[derive(Debug)] +pub struct Segments { + c: [Vec; D], + grid: Grid, +} + +impl Segments<2> { + pub fn moves(&self, coord: usize) -> impl IntoIterator { + let mut occs: BTreeMap<[usize; 2], Vec<&Segment>> = BTreeMap::new(); + + for seg in &self.c[coord] { + let p1 = seg.p_start_2d(coord); + let p2 = seg.p_end_2d(coord); + occs.entry(p1).or_default().push(seg); + occs.entry(p2).or_default().push(seg); + } + + occs.retain(|_, segs| segs.len() > 1); + + let mut moves = self.c[coord].clone(); + + while let Some((joint, segs)) = occs.pop_first() { + debug_assert_eq!(segs.len(), 2); + let seg1 = segs[0]; + let seg2 = segs[1]; + + let mut multi_seg = VecDeque::new(); + + if seg1.p_end_2d(coord) == joint { + // seg1 is before seg2 + multi_seg.push_back(seg1); + multi_seg.push_back(seg2); + } else { + // seg1 is after seg2 + multi_seg.push_back(seg2); + multi_seg.push_back(seg1); + } + + loop { + // Add segments to the left of the multi-segment. + let joint = multi_seg.front().unwrap().p_start_2d(coord); + let Some(segs) = occs.remove(&joint) else { break }; + let seg = *segs + .iter() + .find(|seg| seg.p_end_2d(coord) == joint) + .unwrap(); + multi_seg.push_front(seg); + } + loop { + // Add segments to the right of the multi-segment. + let joint = multi_seg.back().unwrap().p_end_2d(coord); + let Some(segs) = occs.remove(&joint) else { break }; + let seg = *segs + .iter() + .find(|seg| seg.p_start_2d(coord) == joint) + .unwrap(); + multi_seg.push_back(seg); + } + + for i in 0..multi_seg.len() { + // i+1 because "moves" already contains individual segments. + for j in i + 1..multi_seg.len() { + moves.push(Segment { + start: multi_seg[i].start, + end: multi_seg[j].end, + at: multi_seg[i].at, + }); + } + } + } + + moves + } + + pub fn components<'a>( + &'a mut self, + coord: usize, + seg: &'a Segment, + ) -> impl Iterator + 'a { + self.c[coord].iter_mut().filter(|s| seg.contains(s)) + } +} + #[cfg(test)] mod tests { use super::*; diff --git a/src/cartesian/rcb.rs b/src/cartesian/rcb.rs index 5734ce1f..9382021e 100644 --- a/src/cartesian/rcb.rs +++ b/src/cartesian/rcb.rs @@ -1,4 +1,5 @@ use super::Grid; +use super::SplitTree; use super::SubGrid; use num_traits::AsPrimitive; use num_traits::Num; @@ -8,39 +9,6 @@ use rayon::iter::IntoParallelRefIterator; use rayon::iter::ParallelIterator; use std::iter::Sum; -#[derive(Debug)] -pub enum IterationResult { - Whole, - Split { - position: usize, - left: Box, - right: Box, - }, -} - -impl IterationResult { - pub fn part_of(&self, pos: [usize; D], mut start_coord: usize) -> usize { - let mut it = self; - let mut part_id = 0; - while let Self::Split { - position, - left, - right, - } = it - { - if pos[start_coord] < *position { - part_id *= 2; - it = left; - } else { - part_id = 2 * part_id + 1; - it = right; - } - start_coord = (start_coord + 1) % D; - } - part_id - } -} - const TOLERANCE: f64 = 0.01; #[derive(Debug)] @@ -105,13 +73,13 @@ pub(super) fn recurse_2d( total_weight: W, iter_count: usize, coord: usize, -) -> IterationResult +) -> SplitTree where W: Send + Sync + PartialOrd + Num + Sum + AsPrimitive, f64: AsPrimitive, { if subgrid.size[coord] == 0 || iter_count == 0 { - return IterationResult::Whole; + return SplitTree::Whole; } let axis_weights: Vec = if coord == 0 { @@ -170,7 +138,7 @@ where }, ); - IterationResult::Split { + SplitTree::Split { position: split_position, left: Box::new(left), right: Box::new(right), @@ -184,13 +152,13 @@ pub(super) fn recurse_3d( total_weight: W, iter_count: usize, coord: usize, -) -> IterationResult +) -> SplitTree where W: Send + Sync + PartialOrd + Num + Sum + AsPrimitive, f64: AsPrimitive, { if subgrid.size[coord] == 0 || iter_count == 0 { - return IterationResult::Whole; + return SplitTree::Whole; } let axis_weights: Vec = if coord == 0 { @@ -273,7 +241,7 @@ where }, ); - IterationResult::Split { + SplitTree::Split { position: split_position, left: Box::new(left), right: Box::new(right), diff --git a/src/main.rs b/src/main.rs index ca4ac288..a43c0dd0 100644 --- a/src/main.rs +++ b/src/main.rs @@ -9,7 +9,9 @@ fn main() { eprintln!("grid size: ({x},{y}); rcb iters: {iter}"); let grid = coupe::Grid::new_2d(x, y); let n = usize::from(x) * usize::from(y); - let weights: Vec = (0..n).map(|i| i as f64).collect(); + let weights: Vec = (0..n) + .map(|i| if i % x < 50 && i / y < 50 { 2 } else { 3 } as f64) + .collect(); let mut partition = vec![0; n]; let domain = ittapi::Domain::new("MyIncredibleDomain");