diff --git a/.gitignore b/.gitignore index 9d05126..f988492 100644 --- a/.gitignore +++ b/.gitignore @@ -13,3 +13,4 @@ target_wasm .lycheecache .vscode/* +*.patch diff --git a/Cargo.lock b/Cargo.lock index 2695099..0683ccd 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -2386,12 +2386,12 @@ dependencies = [ "futures-intrusive", "gif", "image", - "indexmap", "log", "opener", "palette", "pathfinding", "pollster", + "rayon", "rfd", "serde", "serde-wasm-bindgen", diff --git a/Cargo.toml b/Cargo.toml index 07216f9..f3dcf4f 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -28,8 +28,6 @@ bytemuck = { version = "1.16", features = ["derive"] } image = { version = "0.25.6", features = ["serde"] } palette = "0.7.6" -pathfinding = "4.13.0" -indexmap = "2.11.0" frand = "0.10.1" uuid = { version = "1.18.1", features = ["v4", "js", "serde"] } @@ -47,6 +45,7 @@ serde = { version = "1.0.219", features = ["derive"] } # native: [target.'cfg(not(target_arch = "wasm32"))'.dependencies] env_logger = "0.11.8" +rayon = "1.10" # multithreaded solvers (the web build stays single-threaded) # web: [target.'cfg(target_arch = "wasm32")'.dependencies] @@ -64,6 +63,9 @@ console_error_panic_hook = "0.1.7" # Override wgpu for WASM to use only WebGL backend wgpu = { version = "25.0", features = ["webgl"] } +[dev-dependencies] +pathfinding = "4.13.0" # reference Hungarian implementation for tests + [profile.release] opt-level = 3 # fast and small wasm diff --git a/src/app/calculate/genetic.rs b/src/app/calculate/genetic.rs new file mode 100644 index 0000000..6064e19 --- /dev/null +++ b/src/app/calculate/genetic.rs @@ -0,0 +1,484 @@ +//! The "genetic" algorithm: repeated random pair swaps that are kept if they lower +//! the total cost, with the allowed swap distance shrinking every generation. +//! +//! Native builds run each generation on several threads. All threads share one +//! grid of cells, and a swap only touches its two cells, so each cell carries a +//! tiny spin lock: a thread `try_lock`s both cells of a candidate swap and simply +//! skips the attempt if either is busy. Skipping is harmless for a randomised +//! search and, with thousands of cells and a handful of threads, rare. +//! +//! Because the threads interleave freely, parallel runs are not bit-for-bit +//! reproducible (the single-threaded path still is). + +#[cfg(not(target_arch = "wasm32"))] +use std::sync::{ + Arc, + atomic::{AtomicBool, AtomicI64, AtomicU32, Ordering}, +}; + +use super::util::{self, GenerationSettings, ProgressSink}; +use super::{ProgressMsg, SWAPS_PER_GENERATION_PER_PIXEL, heuristic, make_new_img}; +use crate::app::preset::{Preset, UnprocessedPreset}; + +const SEED: u64 = 12345; + +/// Picks the partner cell for `apos`, exactly like the original: a uniformly random +/// offset within `max_dist` in both axes, clamped to the image. +#[inline(always)] +fn pick_partner(rng: &mut frand::Rand, apos: usize, side: usize, max_dist: u32) -> (u16, u16) { + // 32-bit division: noticeably cheaper than usize on many CPUs, and this is the hot path. + let (ax, ay) = ((apos as u32 % side as u32) as i16, (apos as u32 / side as u32) as i16); + let d = max_dist as i16; + let bx = (ax + rng.gen_range(-d..(d + 1))).clamp(0, side as i16 - 1) as u16; + let by = (ay + rng.gen_range(-d..(d + 1))).clamp(0, side as i16 - 1) as u16; + (bx, by) +} + +/// What to do after a generation: keep going, or hand back the result. +fn finished(max_dist: u32, swaps_made: usize) -> bool { + max_dist < 4 && swaps_made < 10 +} + +fn next_max_dist(max_dist: u32) -> u32 { + (max_dist as f32 * 0.99).max(2.0) as u32 +} + +/// Shared by both paths: report one finished generation, or the final result. +/// Returns true when the run is over. +fn report( + tx: &mut S, + unprocessed: &UnprocessedPreset, + settings: &GenerationSettings, + source_pixels: &[(u8, u8, u8)], + assignments: Vec, + max_dist: u32, + swaps_made: usize, +) -> bool { + if finished(max_dist, swaps_made) { + tx.send(ProgressMsg::Done(Preset { + inner: UnprocessedPreset { + name: unprocessed.name.clone(), + width: settings.sidelen, + height: settings.sidelen, + source_img: source_pixels + .iter() + .flat_map(|(r, g, b)| [*r, *g, *b]) + .collect(), + }, + assignments, + })); + return true; + } + tx.send(ProgressMsg::UpdatePreview { + width: settings.sidelen, + height: settings.sidelen, + data: make_new_img(source_pixels, &assignments, settings.sidelen), + }); + tx.send(ProgressMsg::Progress( + 1.0 - max_dist as f32 / settings.sidelen as f32, + )); + false +} + +// --------------------------------------------------------------------------- +// Sequential path (wasm, or a single core). Same logic as before. +// --------------------------------------------------------------------------- + +#[derive(Clone, Copy)] +struct Pixel { + src_x: u16, + src_y: u16, + rgb: (u8, u8, u8), + h: i64, // current heuristic value +} + +fn process_sequential( + unprocessed: UnprocessedPreset, + settings: GenerationSettings, + tx: &mut S, + #[cfg(not(target_arch = "wasm32"))] cancel: &AtomicBool, +) -> Result<(), Box> { + let source_img = image::ImageBuffer::from_vec( + unprocessed.width, + unprocessed.height, + unprocessed.source_img.clone(), + ) + .unwrap(); + let (source_pixels, target_pixels, weights) = util::get_images(source_img, &settings)?; + let side = settings.sidelen as usize; + let prox = settings.proximity_importance; + + let mut pixels = source_pixels + .iter() + .enumerate() + .map(|(i, &rgb)| { + let x = (i % side) as u16; + let y = (i / side) as u16; + let h = heuristic((x, y), (x, y), rgb, target_pixels[i], weights[i], prox); + Pixel { + src_x: x, + src_y: y, + rgb, + h, + } + }) + .collect::>(); + + let mut rng = frand::Rand::with_seed(SEED); + let swaps_per_generation = SWAPS_PER_GENERATION_PER_PIXEL * pixels.len(); + let mut max_dist = settings.sidelen; + + loop { + let mut swaps_made = 0; + for _ in 0..swaps_per_generation { + let apos = rng.gen_range(0..pixels.len() as u32) as usize; + let (ax, ay) = ((apos as u32 % side as u32) as u16, (apos as u32 / side as u32) as u16); + let (bx, by) = pick_partner(&mut rng, apos, side, max_dist); + let bpos = by as usize * side + bx as usize; + + let (pa, pb) = (pixels[apos], pixels[bpos]); + let a_on_b = heuristic((pa.src_x, pa.src_y), (bx, by), pa.rgb, target_pixels[bpos], weights[bpos], prox); + let b_on_a = heuristic((pb.src_x, pb.src_y), (ax, ay), pb.rgb, target_pixels[apos], weights[apos], prox); + + if (pa.h - b_on_a) + (pb.h - a_on_b) > 0 { + pixels.swap(apos, bpos); + pixels[apos].h = b_on_a; + pixels[bpos].h = a_on_b; + swaps_made += 1; + } + } + + #[cfg(not(target_arch = "wasm32"))] + if cancel.load(Ordering::Relaxed) { + tx.send(ProgressMsg::Cancelled); + return Ok(()); + } + + let assignments = pixels + .iter() + .map(|p| p.src_y as usize * side + p.src_x as usize) + .collect::>(); + if report(tx, &unprocessed, &settings, &source_pixels, assignments, max_dist, swaps_made) { + return Ok(()); + } + max_dist = next_max_dist(max_dist); + } +} + +// --------------------------------------------------------------------------- +// Parallel path (native, more than one thread). +// --------------------------------------------------------------------------- + +/// One grid cell: which source pixel currently sits here, and its cost here. +/// The data fields are relaxed atomics only so the type is `Sync` without +/// `unsafe`; the `lock` (acquire/release) is what actually orders accesses. +#[cfg(not(target_arch = "wasm32"))] +struct Cell { + lock: AtomicBool, + /// Source position of the pixel in this cell: `x << 16 | y`. + src: AtomicU32, + h: AtomicI64, +} + +#[cfg(not(target_arch = "wasm32"))] +impl Cell { + #[inline(always)] + fn try_lock(&self) -> bool { + // Test first so a contended cell doesn't take the cache line exclusively. + !self.lock.load(Ordering::Relaxed) + && self + .lock + .compare_exchange(false, true, Ordering::Acquire, Ordering::Relaxed) + .is_ok() + } + + #[inline(always)] + fn unlock(&self) { + self.lock.store(false, Ordering::Release); + } +} + +#[cfg(not(target_arch = "wasm32"))] +struct Shared<'a> { + side: usize, + prox: i64, + cells: Vec, + source: &'a [(u8, u8, u8)], + target: &'a [(u8, u8, u8)], + weights: &'a [i64], +} + +#[cfg(not(target_arch = "wasm32"))] +impl Shared<'_> { + /// If swapping cells `a` and `b` lowers the total cost, the two new per-cell + /// costs `(b's pixel at a, a's pixel at b)`. Reads the cells as they are now; + /// only exact while the caller holds both locks. + #[inline(always)] + fn improvement( + &self, + ca: &Cell, + cb: &Cell, + (ax, ay): (u16, u16), + (bx, by): (u16, u16), + apos: usize, + bpos: usize, + ) -> Option<(i64, i64)> { + let side = self.side; + let (sa, sb) = (ca.src.load(Ordering::Relaxed), cb.src.load(Ordering::Relaxed)); + let (ha, hb) = (ca.h.load(Ordering::Relaxed), cb.h.load(Ordering::Relaxed)); + let (sax, say) = ((sa >> 16) as u16, (sa & 0xFFFF) as u16); + let (sbx, sby) = ((sb >> 16) as u16, (sb & 0xFFFF) as u16); + let rgb_a = self.source[say as usize * side + sax as usize]; + let rgb_b = self.source[sby as usize * side + sbx as usize]; + + let a_on_b = heuristic((sax, say), (bx, by), rgb_a, self.target[bpos], self.weights[bpos], self.prox); + let b_on_a = heuristic((sbx, sby), (ax, ay), rgb_b, self.target[apos], self.weights[apos], self.prox); + ((ha - b_on_a) + (hb - a_on_b) > 0).then_some((b_on_a, a_on_b)) + } + + #[inline(always)] + fn improves(&self, ca: &Cell, cb: &Cell, a: (u16, u16), b: (u16, u16), apos: usize, bpos: usize) -> bool { + self.improvement(ca, cb, a, b, apos, bpos).is_some() + } + + /// Try `attempts` random swaps; returns how many were made. + fn run(&self, rng: &mut frand::Rand, attempts: usize, max_dist: u32) -> usize { + let side = self.side; + let mut swaps_made = 0; + for _ in 0..attempts { + let apos = rng.gen_range(0..self.cells.len() as u32) as usize; + let (bx, by) = pick_partner(rng, apos, side, max_dist); + let bpos = by as usize * side + bx as usize; + if apos == bpos { + continue; // swapping a cell with itself never improves anything + } + + let (ax, ay) = ((apos as u32 % side as u32) as u16, (apos as u32 / side as u32) as u16); + let (ca, cb) = (&self.cells[apos], &self.cells[bpos]); + + // 1. Optimistic look, no locks. Most attempts don't improve anything, and + // for those this is all we do. The loads may be stale or even mix two + // moments in time (another thread may be mid-swap); that can only make + // us miss a swap or ask for a lock we then don't need, because step 2 + // re-checks everything under the lock. + if !self.improves(ca, cb, (ax, ay), (bx, by), apos, bpos) { + continue; + } + + // 2. Looks good: take both cells (or give up if someone else has one) + // and check again on data nobody else can change. + let (lo, hi) = if apos < bpos { (apos, bpos) } else { (bpos, apos) }; + if !self.cells[lo].try_lock() { + continue; + } + if !self.cells[hi].try_lock() { + self.cells[lo].unlock(); + continue; + } + if let Some((b_on_a, a_on_b)) = self.improvement(ca, cb, (ax, ay), (bx, by), apos, bpos) { + let (sa, sb) = (ca.src.load(Ordering::Relaxed), cb.src.load(Ordering::Relaxed)); + ca.src.store(sb, Ordering::Relaxed); + ca.h.store(b_on_a, Ordering::Relaxed); + cb.src.store(sa, Ordering::Relaxed); + cb.h.store(a_on_b, Ordering::Relaxed); + swaps_made += 1; + } + self.cells[hi].unlock(); + self.cells[lo].unlock(); + } + swaps_made + } +} + +#[cfg(not(target_arch = "wasm32"))] +fn process_parallel( + unprocessed: UnprocessedPreset, + settings: GenerationSettings, + tx: &mut S, + cancel: &AtomicBool, + threads: usize, +) -> Result<(), Box> { + use rayon::prelude::*; + + let source_img = image::ImageBuffer::from_vec( + unprocessed.width, + unprocessed.height, + unprocessed.source_img.clone(), + ) + .unwrap(); + let (source_pixels, target_pixels, weights) = util::get_images(source_img, &settings)?; + let side = settings.sidelen as usize; + let prox = settings.proximity_importance; + + let cells = source_pixels + .iter() + .enumerate() + .map(|(i, &rgb)| { + let (x, y) = ((i % side) as u16, (i / side) as u16); + Cell { + lock: AtomicBool::new(false), + src: AtomicU32::new((x as u32) << 16 | y as u32), + h: AtomicI64::new(heuristic((x, y), (x, y), rgb, target_pixels[i], weights[i], prox)), + } + }) + .collect::>(); + let shared = Shared { + side, + prox, + cells, + source: &source_pixels, + target: &target_pixels, + weights: &weights, + }; + + let n = source_pixels.len(); + let swaps_per_generation = SWAPS_PER_GENERATION_PER_PIXEL * n; + let mut max_dist = settings.sidelen; + let mut generation: u64 = 0; + + loop { + // Split the generation's attempts evenly; each task has its own RNG stream. + let swaps_made: usize = (0..threads) + .into_par_iter() + .map(|t| { + let attempts = swaps_per_generation / threads + + if t < swaps_per_generation % threads { 1 } else { 0 }; + let seed = SEED + .wrapping_add(generation.wrapping_mul(0x9E37_79B9_7F4A_7C15)) + .wrapping_add(t as u64); + shared.run(&mut frand::Rand::with_seed(seed), attempts, max_dist) + }) + .sum(); + generation += 1; + + if cancel.load(Ordering::Relaxed) { + tx.send(ProgressMsg::Cancelled); + return Ok(()); + } + + let assignments = shared + .cells + .iter() + .map(|c| { + let s = c.src.load(Ordering::Relaxed); + (s & 0xFFFF) as usize * side + (s >> 16) as usize + }) + .collect::>(); + if report(tx, &unprocessed, &settings, &source_pixels, assignments, max_dist, swaps_made) { + return Ok(()); + } + max_dist = next_max_dist(max_dist); + } +} + +// --------------------------------------------------------------------------- + +#[cfg(not(target_arch = "wasm32"))] +pub fn process_genetic( + unprocessed: UnprocessedPreset, + settings: GenerationSettings, + tx: &mut S, + cancel: Arc, +) -> Result<(), Box> { + let threads = rayon::current_num_threads(); + // Tiny images aren't worth the coordination. + if threads > 1 && settings.sidelen >= 32 { + process_parallel(unprocessed, settings, tx, &cancel, threads) + } else { + process_sequential(unprocessed, settings, tx, &cancel) + } +} + +#[cfg(target_arch = "wasm32")] +pub fn process_genetic( + unprocessed: UnprocessedPreset, + settings: GenerationSettings, + tx: &mut S, +) -> Result<(), Box> { + process_sequential(unprocessed, settings, tx) +} + +#[cfg(all(test, not(target_arch = "wasm32")))] +mod tests { + use super::*; + use std::sync::Mutex; + use uuid::Uuid; + + fn run(threads: usize, side: u32) -> Vec { + let mut rng = frand::Rand::with_seed(3); + let (w, h) = (48u32, 48u32); + let img: Vec = (0..w * h * 3).map(|_| rng.gen_range(0..256u32) as u8).collect(); + let unprocessed = UnprocessedPreset { + name: "t".into(), + width: w, + height: h, + source_img: img, + }; + let mut settings = GenerationSettings::default(Uuid::nil(), "t".into()); + settings.sidelen = side; + + let out = Arc::new(Mutex::new(None)); + let out2 = out.clone(); + let mut sink = move |m| { + if let ProgressMsg::Done(p) = m { + *out2.lock().unwrap() = Some(p.assignments); + } + }; + let cancel = AtomicBool::new(false); + let pool = rayon::ThreadPoolBuilder::new().num_threads(threads).build().unwrap(); + pool.install(|| { + if threads > 1 { + process_parallel(unprocessed, settings, &mut sink, &cancel, threads).unwrap(); + } else { + process_sequential(unprocessed, settings, &mut sink, &cancel).unwrap(); + } + }); + let result = out.lock().unwrap().take().unwrap(); + result + } + + /// However the threads interleave, every source pixel must end up in exactly one + /// cell: no pixel duplicated, none lost. + #[test] + fn parallel_run_yields_a_permutation() { + let side = 24; + for _ in 0..3 { + let mut a = run(4, side); + a.sort_unstable(); + assert!(a.iter().copied().eq(0..(side * side) as usize)); + } + } + + /// Parallel and sequential runs should land at about the same quality. + #[test] + fn parallel_quality_matches_sequential() { + let side = 24u32; + let cost = |assign: &[usize]| -> i64 { + // Rebuild the problem to score it with the shared heuristic. + let mut rng = frand::Rand::with_seed(3); + let img: Vec = (0..48 * 48 * 3).map(|_| rng.gen_range(0..256u32) as u8).collect(); + let mut settings = GenerationSettings::default(Uuid::nil(), "t".into()); + settings.sidelen = side; + let src = image::ImageBuffer::from_vec(48, 48, img).unwrap(); + let (sp, tp, wt) = util::get_images(src, &settings).unwrap(); + let s = side as usize; + assign + .iter() + .enumerate() + .map(|(i, &c)| { + heuristic( + ((c % s) as u16, (c / s) as u16), + ((i % s) as u16, (i / s) as u16), + sp[c], + tp[i], + wt[i], + settings.proximity_importance, + ) + }) + .sum() + }; + let (seq, par) = (cost(&run(1, side)), cost(&run(4, side))); + let ratio = par as f64 / seq as f64; + assert!((0.97..1.03).contains(&ratio), "seq={seq} par={par} ratio={ratio}"); + } +} diff --git a/src/app/calculate/mod.rs b/src/app/calculate/mod.rs index 717d98d..726effc 100644 --- a/src/app/calculate/mod.rs +++ b/src/app/calculate/mod.rs @@ -2,6 +2,8 @@ use std::sync::{Arc, atomic::AtomicBool}; #[cfg(not(target_arch = "wasm32"))] pub mod drawing_process; +mod genetic; +mod optimal; pub mod util; #[cfg(target_arch = "wasm32")] @@ -19,8 +21,6 @@ use crate::app::{ calculate::util::{GenerationSettings, ProgressSink}, preset::{Preset, UnprocessedPreset}, }; -use egui::ahash::AHasher; -use pathfinding::prelude::Weights; use serde::{Deserialize, Serialize}; #[inline(always)] @@ -39,48 +39,6 @@ fn heuristic( color * color_weight + (spatial * spatial_weight).pow(2) } -struct ImgDiffWeights<'a> { - source: Vec<(u8, u8, u8)>, - target: Vec<(u8, u8, u8)>, - weights: Vec, - sidelen: usize, - settings: &'a GenerationSettings, -} - -// const TARGET_IMAGE_PATH: &str = "./target.png"; -// const TARGET_WEIGHTS_PATH: &str = "./weights.png"; - -impl Weights for ImgDiffWeights<'_> { - fn rows(&self) -> usize { - self.target.len() - } - - fn columns(&self) -> usize { - self.source.len() - } - - #[inline(always)] - fn at(&self, row: usize, col: usize) -> i64 { - let (x1, y1) = (row % self.sidelen, row / self.sidelen); - let (x2, y2) = (col % self.sidelen, col / self.sidelen); - let (r1, g1, b1) = self.target[row]; - let (r2, g2, b2) = self.source[col]; - let weight = self.weights[row]; - -heuristic( - (x1 as u16, y1 as u16), - (x2 as u16, y2 as u16), - (r1, g1, b1), - (r2, g2, b2), - weight, - self.settings.proximity_importance, - ) - } - - fn neg(&self) -> Self { - todo!() - } -} - #[derive(Serialize, Deserialize)] pub enum ProgressMsg { Progress(f32), @@ -108,192 +66,6 @@ impl ProgressMsg { } } -type FxIndexSet = indexmap::IndexSet>; - -pub fn process_optimal( - unprocessed: UnprocessedPreset, - settings: GenerationSettings, - tx: &mut S, - #[cfg(not(target_arch = "wasm32"))] cancel: Arc, -) -> Result<(), Box> { - let source_img = image::ImageBuffer::from_vec( - unprocessed.width, - unprocessed.height, - unprocessed.source_img.clone(), - ) - .unwrap(); - // let start_time = std::time::Instant::now(); - let (source_pixels, target_pixels, weights) = util::get_images(source_img, &settings)?; - - let weights = ImgDiffWeights { - source: source_pixels.clone(), - target: target_pixels, - weights, - sidelen: settings.sidelen as usize, - settings: &settings, - }; - - // pathfinding::kuhn_munkres, inlined to allow for progress bar and cancelling - let (_total_diff, assignments) = { - // We call x the rows and y the columns. (nx, ny) is the size of the matrix. - let nx = weights.rows(); - let ny = weights.columns(); - assert!( - nx <= ny, - "number of rows must not be larger than number of columns" - ); - // xy represents matching for x, yz matching for y - let mut xy: Vec> = vec![None; nx]; - let mut yx: Vec> = vec![None; ny]; - // lx is the labelling for x nodes, ly the labelling for y nodes. We start - // with an acceptable labelling with the maximum possible values for lx - // and 0 for ly. - let mut lx: Vec = (0..nx) - .map(|row| (0..ny).map(|col| weights.at(row, col)).max().unwrap()) - .collect::>(); - let mut ly: Vec = vec![0; ny]; - // s, augmenting, and slack will be reset every time they are reused. augmenting - // contains Some(prev) when the corresponding node belongs to the augmenting path. - let mut s = FxIndexSet::::default(); - let mut alternating = Vec::with_capacity(ny); - let mut slack = vec![0; ny]; - let mut slackx = Vec::with_capacity(ny); - for root in 0..nx { - alternating.clear(); - alternating.resize(ny, None); - // Find y such that the path is augmented. This will be set when breaking for the - // loop below. Above the loop is some code to initialize the search. - let mut y = { - s.clear(); - s.insert(root); - // Slack for a vertex y is, initially, the margin between the - // sum of the labels of root and y, and the weight between root and y. - // As we add x nodes to the alternating path, we update the slack to - // represent the smallest margin between one of the x nodes and y. - for y in 0..ny { - slack[y] = lx[root] + ly[y] - weights.at(root, y); - } - slackx.clear(); - slackx.resize(ny, root); - Some(loop { - let mut delta = pathfinding::num_traits::Bounded::max_value(); - let mut x = 0; - let mut y = 0; - // Select one of the smallest slack delta and its edge (x, y) - // for y not in the alternating path already. - for yy in 0..ny { - if alternating[yy].is_none() && slack[yy] < delta { - delta = slack[yy]; - x = slackx[yy]; - y = yy; - } - } - // If some slack has been found, remove it from x nodes in the - // alternating path, and add it to y nodes in the alternating path. - // The slack of y nodes outside the alternating path will be reduced - // by this minimal slack as well. - if delta > 0 { - for &x in &s { - lx[x] -= delta; - } - for y in 0..ny { - if alternating[y].is_some() { - ly[y] += delta; - } else { - slack[y] -= delta; - } - } - } - // Add (x, y) to the alternating path. - alternating[y] = Some(x); - if yx[y].is_none() { - // We have found an augmenting path. - break y; - } - // This y node had a predecessor, add it to the set of x nodes - // in the augmenting path. - let x = yx[y].unwrap(); - s.insert(x); - // Update slack because of the added vertex in s might contain a - // greater slack than with previously inserted x nodes in the augmenting - // path. - for y in 0..ny { - if alternating[y].is_none() { - let alternate_slack = lx[x] + ly[y] - weights.at(x, y); - if slack[y] > alternate_slack { - slack[y] = alternate_slack; - slackx[y] = x; - } - } - } - }) - }; - // Inverse edges along the augmenting path. - while y.is_some() { - let x = alternating[y.unwrap()].unwrap(); - let prec = xy[x]; - yx[y.unwrap()] = Some(x); - xy[x] = y; - y = prec; - } - if root % 100 == 0 { - // send progress - #[cfg(not(target_arch = "wasm32"))] - { - if cancel.load(std::sync::atomic::Ordering::Relaxed) { - tx.send(ProgressMsg::Cancelled); - return Ok(()); - } - } - - tx.send(ProgressMsg::Progress(root as f32 / nx as f32)); - - let data = make_new_img( - &source_pixels, - &xy.clone() - .into_iter() - .map(|a| a.unwrap_or(0)) - .collect::>(), - settings.sidelen, - ); - - tx.send(ProgressMsg::UpdatePreview { - width: settings.sidelen, - height: settings.sidelen, - data, - }); - } - } - ( - lx.into_iter().sum::() + ly.into_iter().sum::(), - xy.into_iter().map(Option::unwrap).collect::>(), - ) - }; - - //let img = make_new_img(&source_pixels, &assignments, target.width()); - - //let dir_name = util::save_result(target, "todo".to_string(), source, assignments, img)?; - - tx.send(ProgressMsg::Done(Preset { - inner: UnprocessedPreset { - name: unprocessed.name, - width: settings.sidelen, - height: settings.sidelen, - source_img: source_pixels - .into_iter() - .flat_map(|(r, g, b)| [r, g, b]) - .collect(), - }, - assignments: assignments.clone(), - })); - - // println!( - // "finished in {:.2?} seconds", - // std::time::Instant::now().duration_since(start_time) - // ); - Ok(()) -} - fn make_new_img(source_pixels: &[(u8, u8, u8)], assignments: &[usize], sidelen: u32) -> Vec { let mut img = vec![0; (sidelen * sidelen * 3) as usize]; for (target_idx, source_idx) in assignments.iter().enumerate() { @@ -306,191 +78,43 @@ fn make_new_img(source_pixels: &[(u8, u8, u8)], assignments: &[usize], sidelen: img } -#[derive(Clone, Copy)] -struct Pixel { - src_x: u16, - src_y: u16, - rgb: (u8, u8, u8), - h: i64, // current heuristic value -} - -impl Pixel { - fn new(src_x: u16, src_y: u16, rgb: (u8, u8, u8), h: i64) -> Self { - Self { - src_x, - src_y, - rgb, - h, - } - } - - fn update_heuristic(&mut self, new_h: i64) { - self.h = new_h; - } - - #[inline(always)] - fn calc_heuristic( - &self, - target_pos: (u16, u16), - target_col: (u8, u8, u8), - weight: i64, - proximity_importance: i64, - ) -> i64 { - heuristic( - (self.src_x, self.src_y), - target_pos, - self.rgb, - target_col, - weight, - proximity_importance, - ) - } -} - const SWAPS_PER_GENERATION_PER_PIXEL: usize = 128; -pub fn process_genetic( - unprocessed: UnprocessedPreset, - settings: GenerationSettings, - tx: &mut S, - #[cfg(not(target_arch = "wasm32"))] cancel: Arc, -) -> Result<(), Box> { - let source_img = image::ImageBuffer::from_vec( - unprocessed.width, - unprocessed.height, - unprocessed.source_img.clone(), - ) - .unwrap(); - // let start_time = std::time::Instant::now(); - let (source_pixels, target_pixels, weights) = util::get_images(source_img, &settings)?; - - let mut pixels = source_pixels - .iter() - .enumerate() - .map(|(i, &(r, g, b))| { - let x = (i as u32 % settings.sidelen) as u16; - let y = (i as u32 / settings.sidelen) as u16; - let mut p = Pixel::new(x, y, (r, g, b), 0); - let h = p.calc_heuristic( - (x, y), - target_pixels[i], - weights[i], - settings.proximity_importance, - ); - p.update_heuristic(h); - p - }) - .collect::>(); - - let mut rng = frand::Rand::with_seed(12345); - let swaps_per_generation = SWAPS_PER_GENERATION_PER_PIXEL * pixels.len(); - - let mut max_dist = settings.sidelen; - loop { - let mut swaps_made = 0; - for _ in 0..swaps_per_generation { - let apos = rng.gen_range(0..pixels.len() as u32) as usize; - let ax = apos as u16 % settings.sidelen as u16; - let ay = apos as u16 / settings.sidelen as u16; - let bx = (ax as i16 + rng.gen_range(-(max_dist as i16)..(max_dist as i16 + 1))) - .clamp(0, settings.sidelen as i16 - 1) as u16; - let by = (ay as i16 + rng.gen_range(-(max_dist as i16)..(max_dist as i16 + 1))) - .clamp(0, settings.sidelen as i16 - 1) as u16; - let bpos = by as usize * settings.sidelen as usize + bx as usize; - - let t_a = target_pixels[apos]; - let t_b = target_pixels[bpos]; - - let a_on_b_h = pixels[apos].calc_heuristic( - (bx, by), - t_b, - weights[bpos], - settings.proximity_importance, - ); - - let b_on_a_h = pixels[bpos].calc_heuristic( - (ax, ay), - t_a, - weights[apos], - settings.proximity_importance, - ); - - let improvement_a = pixels[apos].h - b_on_a_h; - let improvement_b = pixels[bpos].h - a_on_b_h; - if improvement_a + improvement_b > 0 { - // swap - pixels.swap(apos, bpos); - pixels[apos].update_heuristic(b_on_a_h); - pixels[bpos].update_heuristic(a_on_b_h); - swaps_made += 1; - } - } - - #[cfg(not(target_arch = "wasm32"))] - { - if cancel.load(std::sync::atomic::Ordering::Relaxed) { - println!("cancelled"); - tx.send(ProgressMsg::Cancelled); - return Ok(()); - } - } - - let assignments = pixels - .iter() - .map(|p| p.src_y as usize * settings.sidelen as usize + p.src_x as usize) - .collect::>(); - //debug_print(format!("max_dist = {max_dist}, swaps made = {swaps_made}")); - if max_dist < 4 && swaps_made < 10 { - //let dir_name = util::save_result(target, base_name, source, assignments, img)?; - tx.send(ProgressMsg::Done(Preset { - inner: UnprocessedPreset { - name: unprocessed.name, - width: settings.sidelen, - height: settings.sidelen, - source_img: source_pixels - .iter() - .flat_map(|(r, g, b)| [*r, *g, *b]) - .collect(), - }, - assignments: assignments.clone(), - })); - return Ok(()); - } - let data = make_new_img(&source_pixels, &assignments, settings.sidelen); - tx.send(ProgressMsg::UpdatePreview { - width: settings.sidelen, - height: settings.sidelen, - data, - }); - tx.send(ProgressMsg::Progress( - 1.0 - max_dist as f32 / settings.sidelen as f32, - )); - - max_dist = (max_dist as f32 * 0.99).max(2.0) as u32; - } +/// How many threads the native solvers may use: all cores by default, or whatever +/// `OBAMIFY_THREADS` says (`1` forces the single-threaded code paths). +#[cfg(not(target_arch = "wasm32"))] +fn worker_threads() -> usize { + std::env::var("OBAMIFY_THREADS") + .ok() + .and_then(|v| v.parse::().ok()) + .filter(|&n| n > 0) + .unwrap_or_else(|| std::thread::available_parallelism().map_or(1, |n| n.get())) } -// fn serialize_assignments(assignments: Vec) -> String { -// format!( -// "[{}]", -// assignments -// .iter() -// .map(|a| a.to_string()) -// .collect::>() -// .join(",") -// ) -// } #[cfg(not(target_arch = "wasm32"))] -pub fn process( +pub fn process( unprocessed: UnprocessedPreset, settings: GenerationSettings, tx: &mut S, cancel: Arc, ) -> Result<(), Box> { - match settings.algorithm { - Algorithm::Optimal => process_optimal(unprocessed, settings, tx, cancel), - Algorithm::Genetic => process_genetic(unprocessed, settings, tx, cancel), - } + // A private pool, so the solvers' worker count doesn't depend on (or fight with) + // anything else in the app. `install` runs the closure *inside* the pool, which + // makes the solvers' fork/join calls much cheaper than entering from outside. + let pool = rayon::ThreadPoolBuilder::new() + .num_threads(worker_threads()) + .thread_name(|i| format!("obamify-worker-{i}")) + .build()?; + + // Box isn't Send, so only a String crosses the pool boundary. + pool.install(move || { + match settings.algorithm { + Algorithm::Optimal => optimal::process_optimal(unprocessed, settings, tx, cancel), + Algorithm::Genetic => genetic::process_genetic(unprocessed, settings, tx, cancel), + } + .map_err(|e| e.to_string()) + }) + .map_err(Into::into) } #[cfg(target_arch = "wasm32")] @@ -499,8 +123,9 @@ pub fn process( settings: GenerationSettings, tx: &mut S, ) -> Result<(), Box> { + // The web build runs in a single worker without shared-memory threads. match settings.algorithm { - Algorithm::Optimal => process_optimal(unprocessed, settings, tx), - Algorithm::Genetic => process_genetic(unprocessed, settings, tx), + Algorithm::Optimal => optimal::process_optimal(unprocessed, settings, tx), + Algorithm::Genetic => genetic::process_genetic(unprocessed, settings, tx), } } diff --git a/src/app/calculate/optimal.rs b/src/app/calculate/optimal.rs new file mode 100644 index 0000000..57f865f --- /dev/null +++ b/src/app/calculate/optimal.rs @@ -0,0 +1,525 @@ +//! The "optimal" algorithm: an exact minimum-cost assignment between source and +//! target cells (the Hungarian method, in its shortest-augmenting-path form). +//! +//! Compared to the previous inlined `kuhn_munkres`, this version +//! +//! * does a **single pass** over the columns per search step (cost evaluation, +//! slack update and arg-min fused) instead of three, by keeping the dual +//! variables *lazy* (they are only brought up to date once per augmentation), +//! * stores the columns in struct-of-arrays form and hoists everything that only +//! depends on the row out of the inner loop, so the hot loop has no division, +//! modulo or bounds-check-heavy indexing, +//! * needs no O(n²) initialisation pass (all costs are >= 0, so zero duals are +//! already feasible), +//! * checks for cancellation after every row instead of every 100th, and +//! * on native targets, splits the per-step scan across a rayon thread pool. +//! +//! The result is still an exact optimum. Ties are broken towards the lowest +//! column index both sequentially and in parallel, so the output does not depend +//! on how many threads are used. + +#[cfg(not(target_arch = "wasm32"))] +use std::sync::{Arc, atomic::AtomicBool, atomic::Ordering}; + +use super::util::{self, GenerationSettings, ProgressSink}; +use super::{ProgressMsg, make_new_img}; +use crate::app::preset::{Preset, UnprocessedPreset}; + +const INF: i64 = i64::MAX / 4; +/// Marks "no row" / "no column" / "the root of the search tree". +const NONE: u32 = u32::MAX; + +/// Below this many columns, splitting a step across threads costs more than it saves. +#[cfg(not(target_arch = "wasm32"))] +const MIN_COLS_FOR_PARALLEL: usize = 4096; +/// Smallest slice of columns handed to one thread. +#[cfg(not(target_arch = "wasm32"))] +const MIN_CHUNK: usize = 1024; + +/// Everything that is fixed for the whole run. Rows are target cells, columns are +/// source cells. The cost of putting source cell `j` at target cell `i` is +/// `heuristic(..)`, which we minimise. +struct Problem<'a> { + side: usize, + spatial_weight: i64, + target: &'a [(u8, u8, u8)], + weights: &'a [i64], + // Source cells, struct-of-arrays so the hot loop streams through memory. + col_x: Vec, + col_y: Vec, + col_r: Vec, + col_g: Vec, + col_b: Vec, +} + +/// The parts of a row that the inner loop needs. +#[derive(Clone, Copy)] +struct RowCtx { + x: i32, + y: i32, + r: i32, + g: i32, + b: i32, + color_weight: i64, +} + +impl<'a> Problem<'a> { + fn new( + source: &[(u8, u8, u8)], + target: &'a [(u8, u8, u8)], + weights: &'a [i64], + side: usize, + spatial_weight: i64, + ) -> Self { + let n = source.len(); + let mut p = Problem { + side, + spatial_weight, + target, + weights, + col_x: Vec::with_capacity(n), + col_y: Vec::with_capacity(n), + col_r: Vec::with_capacity(n), + col_g: Vec::with_capacity(n), + col_b: Vec::with_capacity(n), + }; + for (j, &(r, g, b)) in source.iter().enumerate() { + p.col_x.push((j % side) as i32); + p.col_y.push((j / side) as i32); + p.col_r.push(r as i32); + p.col_g.push(g as i32); + p.col_b.push(b as i32); + } + p + } + + #[inline] + fn row_ctx(&self, row: usize) -> RowCtx { + let (r, g, b) = self.target[row]; + RowCtx { + x: (row % self.side) as i32, + y: (row / self.side) as i32, + r: r as i32, + g: g as i32, + b: b as i32, + color_weight: self.weights[row], + } + } + +} + +/// Same value as `heuristic(..)` for (row, column), with the row-dependent parts +/// already in `c`. Kept in i64 so results are bit-identical to the old code. +#[inline(always)] +fn cost_of(c: &RowCtx, spatial_weight: i64, x: i32, y: i32, r: i32, g: i32, b: i32) -> i64 { + let dx = (c.x - x) as i64; + let dy = (c.y - y) as i64; + let dr = (c.r - r) as i64; + let dg = (c.g - g) as i64; + let db = (c.b - b) as i64; + let spatial = dx * dx + dy * dy; + let color = dr * dr + dg * dg + db * db; + color * c.color_weight + (spatial * spatial_weight).pow(2) +} + +/// Relax the columns `start .. start + minv.len()` against one new row. +/// +/// For every column not yet in the search tree this lowers `minv[j]` (and records +/// where it came from in `way[j]`) if this row offers a cheaper reduced cost, and +/// returns the smallest `(minv, column)` seen, preferring the lowest column on ties. +/// +/// `offset` is `acc - u[row]`; see the lazy-dual notes in `solve`. +#[inline] +fn relax_range( + prob: &Problem<'_>, + ctx: &RowCtx, + offset: i64, + prev_col: u32, + v: &[i64], + used: &[bool], + start: usize, + minv: &mut [i64], + way: &mut [u32], +) -> (i64, usize) { + // Re-slice everything to exactly this chunk so the loop has no bounds checks. + let len = minv.len(); + let end = start + len; + let (cx, cy) = (&prob.col_x[start..end], &prob.col_y[start..end]); + let (cr, cg, cb) = ( + &prob.col_r[start..end], + &prob.col_g[start..end], + &prob.col_b[start..end], + ); + let (v, used, way) = (&v[start..end], &used[start..end], &mut way[..len]); + let sw = prob.spatial_weight; + + let mut best = (INF, usize::MAX); + for k in 0..len { + if used[k] { + continue; + } + let cur = cost_of(ctx, sw, cx[k], cy[k], cr[k], cg[k], cb[k]) - v[k] + offset; + if cur < minv[k] { + minv[k] = cur; + way[k] = prev_col; + } + // Strict `<` keeps the first (lowest) column among equal minima. + if minv[k] < best.0 { + best = (minv[k], start + k); + } + } + best +} + +/// Scratch state for one run of the solver, plus the (optional) thread count. +struct Solver<'a> { + prob: &'a Problem<'a>, + /// Column potentials (only brought up to date at the end of each augmentation). + v: Vec, + /// Row potentials. + u: Vec, + col_row: Vec, + row_col: Vec, + // Per-search scratch. + minv: Vec, + way: Vec, + used: Vec, + #[cfg(not(target_arch = "wasm32"))] + chunk: usize, // 0 = run sequentially +} + +impl<'a> Solver<'a> { + fn new(prob: &'a Problem<'a>, n_rows: usize) -> Self { + let n_cols = prob.col_x.len(); + #[cfg(not(target_arch = "wasm32"))] + let chunk = { + let threads = rayon::current_num_threads(); + if threads > 1 && n_cols >= MIN_COLS_FOR_PARALLEL { + // Round up to a multiple of 64 so chunks don't share cache lines. + (n_cols.div_ceil(threads).max(MIN_CHUNK)).next_multiple_of(64) + } else { + 0 + } + }; + Solver { + prob, + v: vec![0; n_cols], + u: vec![0; n_rows], + col_row: vec![NONE; n_cols], + row_col: vec![NONE; n_rows], + minv: vec![INF; n_cols], + way: vec![NONE; n_cols], + used: vec![false; n_cols], + #[cfg(not(target_arch = "wasm32"))] + chunk, + } + } + + /// One relaxation step over every column, in parallel when worthwhile. + fn relax(&mut self, ctx: &RowCtx, offset: i64, prev_col: u32) -> (i64, usize) { + #[cfg(not(target_arch = "wasm32"))] + if self.chunk > 0 { + use rayon::prelude::*; + let chunk = self.chunk; + let prob = self.prob; + let v = &self.v[..]; + let used = &self.used[..]; + return self + .minv + .par_chunks_mut(chunk) + .zip(self.way.par_chunks_mut(chunk)) + .enumerate() + .map(|(ci, (minv, way))| { + relax_range(prob, ctx, offset, prev_col, v, used, ci * chunk, minv, way) + }) + // Lexicographic min on (value, column) == what the sequential scan picks. + .reduce(|| (INF, usize::MAX), |a, b| if b < a { b } else { a }); + } + + relax_range( + self.prob, + ctx, + offset, + prev_col, + &self.v, + &self.used, + 0, + &mut self.minv, + &mut self.way, + ) + } + + /// Find a minimum-cost augmenting path from the free row `root` and apply it. + /// + /// This is Dijkstra over reduced costs `c(i,j) - u[i] - v[j]`. The textbook + /// version shifts the potentials of every column on each step; instead we keep + /// `acc`, the total shift so far, and store `minv[j]` as an absolute value + /// (`true_minv[j] + acc`). A column's potential then only needs fixing once, + /// at the end, by the amount of shift that happened after it joined the tree. + fn augment(&mut self, root: usize) { + self.minv.fill(INF); + self.used.fill(false); + + // Columns that joined the tree, with the value of `acc` when they did. + let mut visited: Vec<(usize, i64)> = Vec::new(); + + let mut acc: i64 = 0; + let mut row = root; + let mut prev_col = NONE; + let end_col = loop { + let ctx = self.prob.row_ctx(row); + let offset = acc - self.u[row]; + let (m, j) = self.relax(&ctx, offset, prev_col); + debug_assert!(j != usize::MAX, "no unused column left to relax"); + + acc = m; // everything in the tree is shifted by (m - acc_old) + if self.col_row[j] == NONE { + break j; // free column: augmenting path found + } + self.used[j] = true; + visited.push((j, acc)); + row = self.col_row[j] as usize; + prev_col = j as u32; + }; + + // Bring the potentials of the tree up to date. + self.u[root] += acc; + for &(j, a) in &visited { + let shift = acc - a; + self.v[j] -= shift; + self.u[self.col_row[j] as usize] += shift; + } + + // Flip the matching along the path end_col <- ... <- root. + let mut j = end_col; + loop { + let prev = self.way[j]; + let r = if prev == NONE { + root + } else { + self.col_row[prev as usize] as usize + }; + self.col_row[j] = r as u32; + self.row_col[r] = j as u32; + if prev == NONE { + break; + } + j = prev as usize; + } + } +} + +pub fn process_optimal( + unprocessed: UnprocessedPreset, + settings: GenerationSettings, + tx: &mut S, + #[cfg(not(target_arch = "wasm32"))] cancel: Arc, +) -> Result<(), Box> { + let source_img = image::ImageBuffer::from_vec( + unprocessed.width, + unprocessed.height, + unprocessed.source_img.clone(), + ) + .unwrap(); + let (source_pixels, target_pixels, weights) = util::get_images(source_img, &settings)?; + + let side = settings.sidelen as usize; + let nx = target_pixels.len(); + let prob = Problem::new( + &source_pixels, + &target_pixels, + &weights, + side, + settings.proximity_importance, + ); + assert!(nx <= prob.col_x.len(), "more target cells than source cells"); + + let mut solver = Solver::new(&prob, nx); + + for root in 0..nx { + #[cfg(not(target_arch = "wasm32"))] + if cancel.load(Ordering::Relaxed) { + tx.send(ProgressMsg::Cancelled); + return Ok(()); + } + + solver.augment(root); + + if root % 100 == 0 { + tx.send(ProgressMsg::Progress(root as f32 / nx as f32)); + + // Rows after `root` are not matched yet; show them as column 0 like before. + let assignments = solver + .row_col + .iter() + .map(|&c| if c == NONE { 0 } else { c as usize }) + .collect::>(); + tx.send(ProgressMsg::UpdatePreview { + width: settings.sidelen, + height: settings.sidelen, + data: make_new_img(&source_pixels, &assignments, settings.sidelen), + }); + } + } + + let assignments = solver + .row_col + .iter() + .map(|&c| c as usize) + .collect::>(); + + tx.send(ProgressMsg::Done(Preset { + inner: UnprocessedPreset { + name: unprocessed.name, + width: settings.sidelen, + height: settings.sidelen, + source_img: source_pixels + .into_iter() + .flat_map(|(r, g, b)| [r, g, b]) + .collect(), + }, + assignments, + })); + Ok(()) +} + +#[cfg(test)] +mod tests { + use super::*; + use crate::app::calculate::heuristic; + use pathfinding::prelude::{Matrix, kuhn_munkres}; + + type Rgb = (u8, u8, u8); + + fn random_instance(side: usize, seed: u64) -> (Vec, Vec, Vec) { + let mut rng = frand::Rand::with_seed(seed); + let mut px = |n: usize| -> Vec { + (0..n) + .map(|_| { + ( + rng.gen_range(0..256u32) as u8, + rng.gen_range(0..256u32) as u8, + rng.gen_range(0..256u32) as u8, + ) + }) + .collect() + }; + let n = side * side; + let (source, target) = (px(n), px(n)); + let weights = (0..n).map(|_| rng.gen_range(1..256u32) as i64).collect(); + (source, target, weights) + } + + fn total_cost( + assignment: &[usize], + source: &[Rgb], + target: &[Rgb], + weights: &[i64], + side: usize, + prox: i64, + ) -> i64 { + assignment + .iter() + .enumerate() + .map(|(row, &col)| { + heuristic( + ((row % side) as u16, (row / side) as u16), + ((col % side) as u16, (col / side) as u16), + target[row], + source[col], + weights[row], + prox, + ) + }) + .sum() + } + + fn solve(prob: &Problem, n: usize, chunk: Option) -> Vec { + let mut solver = Solver::new(prob, n); + #[cfg(not(target_arch = "wasm32"))] + if let Some(c) = chunk { + solver.chunk = c; + } + let _ = chunk; + for root in 0..n { + solver.augment(root); + } + solver.row_col.iter().map(|&c| c as usize).collect() + } + + fn is_permutation(a: &[usize]) -> bool { + let mut seen = vec![false; a.len()]; + a.iter().all(|&c| c < a.len() && !std::mem::replace(&mut seen[c], true)) + } + + /// Same optimum as the old pathfinding-based implementation. + #[test] + fn matches_reference_kuhn_munkres() { + for (side, prox) in [(2, 13), (3, 13), (5, 1), (6, 13), (8, 40), (10, 0)] { + for seed in 0..4 { + let (source, target, weights) = random_instance(side, seed * 31 + side as u64); + let n = side * side; + let prob = Problem::new(&source, &target, &weights, side, prox); + let ours = solve(&prob, n, None); + assert!(is_permutation(&ours)); + + let matrix = Matrix::from_vec( + n, + n, + (0..n * n) + .map(|k| { + let (row, col) = (k / n, k % n); + -heuristic( + ((row % side) as u16, (row / side) as u16), + ((col % side) as u16, (col / side) as u16), + target[row], + source[col], + weights[row], + prox, + ) + }) + .collect(), + ) + .unwrap(); + let (best, _) = kuhn_munkres(&matrix); + + assert_eq!( + total_cost(&ours, &source, &target, &weights, side, prox), + -best, + "side={side} prox={prox} seed={seed}" + ); + } + } + } + + /// Splitting the scan over threads must not change the answer at all. + #[cfg(not(target_arch = "wasm32"))] + #[test] + fn parallel_scan_is_identical_to_sequential() { + let pool = rayon::ThreadPoolBuilder::new().num_threads(4).build().unwrap(); + for side in [7, 12, 20] { + let (source, target, weights) = random_instance(side, 99 + side as u64); + let n = side * side; + let prob = Problem::new(&source, &target, &weights, side, 13); + let seq = solve(&prob, n, None); + for chunk in [1, 5, 16, 64] { + let par = pool.install(|| solve(&prob, n, Some(chunk))); + assert_eq!(seq, par, "side={side} chunk={chunk}"); + } + } + } + + /// Rectangular case from the old code's assert: fewer rows than columns is fine. + #[test] + fn fewer_rows_than_columns() { + let side = 4; + let (source, target, weights) = random_instance(side, 5); + let prob = Problem::new(&source, &target[..10], &weights, side, 13); + let a = solve(&prob, 10, None); + assert!(a.iter().all(|&c| c < 16)); + let mut s = a.clone(); + s.sort(); + s.dedup(); + assert_eq!(s.len(), 10); + } +}