diff --git a/rust/Cargo.toml b/rust/Cargo.toml index 39123b1..bd34f74 100644 --- a/rust/Cargo.toml +++ b/rust/Cargo.toml @@ -3,13 +3,6 @@ name = "opencubes" version = "0.1.0" edition = "2021" -#feature flags to enable and disable features at compile time - -[features] -diagnostics = [] -size16 = [] -smallset = [] - [profile.release] debug = true diff --git a/rust/src/cli/cli.rs b/rust/src/cli/cli.rs index b298e5d..3b3bd2f 100644 --- a/rust/src/cli/cli.rs +++ b/rust/src/cli/cli.rs @@ -6,7 +6,7 @@ use std::{ use clap::{Args, Parser, Subcommand, ValueEnum}; use indicatif::{MultiProgress, ProgressBar, ProgressStyle}; -use opencubes::{naive_polycube::NaivePolyCube, pcube::PCubeFile}; +use opencubes::polycubes::{naive_polycube::NaivePolyCube, pcube::PCubeFile}; use rayon::prelude::{IntoParallelIterator, ParallelIterator}; mod enumerate; @@ -188,11 +188,11 @@ pub enum Compression { Gzip, } -impl From for opencubes::pcube::Compression { +impl From for opencubes::polycubes::pcube::Compression { fn from(value: Compression) -> Self { match value { - Compression::None => opencubes::pcube::Compression::None, - Compression::Gzip => opencubes::pcube::Compression::Gzip, + Compression::None => opencubes::polycubes::pcube::Compression::None, + Compression::Gzip => opencubes::polycubes::pcube::Compression::Gzip, } } } diff --git a/rust/src/cli/enumerate.rs b/rust/src/cli/enumerate.rs index 0901896..d20c2ec 100644 --- a/rust/src/cli/enumerate.rs +++ b/rust/src/cli/enumerate.rs @@ -3,10 +3,12 @@ use std::{io::ErrorKind, sync::Arc, time::Instant}; use opencubes::{ hashless::MapStore, iterator::{indicatif::PolycubeProgressBarIter, *}, - naive_polycube::NaivePolyCube, - pcube::{PCubeFile, RawPCube}, pointlist, - polycube_reps::CubeMapPos, + polycubes::{ + naive_polycube::NaivePolyCube, + pcube::{PCubeFile, RawPCube}, + point_list::CubeMapPos, + }, rotation_reduced, }; @@ -83,7 +85,7 @@ fn load_cache_file(n: usize) -> Option { } enum CacheOrbase { - Cache(opencubes::pcube::AllUnique), + Cache(opencubes::polycubes::pcube::AllUnique), Base(bool), } diff --git a/rust/src/hashless.rs b/rust/src/hashless.rs index 28d98f9..bae3816 100644 --- a/rust/src/hashless.rs +++ b/rust/src/hashless.rs @@ -1,110 +1,14 @@ -use std::cmp::max; - use hashbrown::HashSet; -use crate::{ - pointlist::{array_insert, array_shift}, - polycube_reps::{CubeMapPos, Dim}, - rotations::{rot_matrix_points, to_min_rot_points, MatrixCol}, +use crate::polycubes::{ + point_list::{CubeMapPos, PointListMeta}, + Dim, PolyCube, }; pub struct MapStore { inner: HashSet>, } -macro_rules! cube_map_pos_expand { - ($name:ident, $dim:ident, $shift:literal) => { - #[inline(always)] - pub fn $name<'a>( - &'a self, - shape: &'a Dim, - count: usize, - ) -> impl Iterator + 'a { - struct Iter<'a, const C: usize> { - inner: &'a CubeMapPos, - shape: &'a Dim, - count: usize, - i: usize, - stored: Option<(Dim, usize, CubeMapPos)>, - } - - impl<'a, const C: usize> Iterator for Iter<'a, C> { - type Item = (Dim, usize, CubeMapPos); - - fn next(&mut self) -> Option { - loop { - if let Some(stored) = self.stored.take() { - return Some(stored); - } - - let i = self.i; - - if i == self.count { - return None; - } - - self.i += 1; - let coord = *self.inner.cubes.get(i)?; - - let plus = coord + (1 << $shift); - let minus = coord - (1 << $shift); - - if !self.inner.cubes[(i + 1)..self.count].contains(&plus) { - let mut new_shape = *self.shape; - let mut new_map = *self.inner; - - array_insert(plus, &mut new_map.cubes[i..=self.count]); - new_shape.$dim = - max(new_shape.$dim, (((coord >> $shift) + 1) & 0x1f) as usize); - - self.stored = Some((new_shape, self.count + 1, new_map)); - } - - let mut new_map = *self.inner; - let mut new_shape = *self.shape; - - // If the coord is out of bounds for $dim, shift everything - // over and create the cube at the out-of-bounds position. - // If it is in bounds, check if the $dim - 1 value already - // exists. - let insert_coord = if (coord >> $shift) & 0x1f != 0 { - if !self.inner.cubes[0..i].contains(&minus) { - minus - } else { - continue; - } - } else { - new_shape.$dim += 1; - for i in 0..self.count { - new_map.cubes[i] += 1 << $shift; - } - coord - }; - - array_shift(&mut new_map.cubes[i..=self.count]); - array_insert(insert_coord, &mut new_map.cubes[0..=i]); - return Some((new_shape, self.count + 1, new_map)); - } - } - } - - Iter { - inner: self, - shape, - count, - i: 0, - stored: None, - } - } - }; -} - -impl CubeMapPos { - cube_map_pos_expand!(expand_x, x, 0); - cube_map_pos_expand!(expand_y, y, 5); - cube_map_pos_expand!(expand_z, z, 10); -} - impl MapStore { pub fn new() -> Self { Self { @@ -114,59 +18,13 @@ impl MapStore { /// helper function to not duplicate code for canonicalising polycubes /// and storing them in the hashset - fn insert_map(&mut self, dim: &Dim, map: &CubeMapPos, count: usize) { - if !self.inner.contains(map) { - let map = to_min_rot_points(map, dim, count); + fn insert_map(&mut self, dim: Dim, map: CubeMapPos, count: usize) { + if !self.inner.contains(&map) { + let map = map.to_min_rot_points(dim, count); self.inner.insert(map); } } - /// reduce number of expansions needing to be performed based on - /// X >= Y >= Z constraint on Dim - #[inline] - fn do_cube_expansion(&mut self, seed: &CubeMapPos, shape: &Dim, count: usize) { - let expand_ys = if shape.y < shape.x { - Some(seed.expand_y(shape, count)) - } else { - None - }; - - let expand_zs = if shape.z < shape.y { - Some(seed.expand_z(shape, count)) - } else { - None - }; - - seed.expand_x(shape, count) - .chain(expand_ys.into_iter().flatten()) - .chain(expand_zs.into_iter().flatten()) - .for_each(|(dim, new_count, map)| self.insert_map(&dim, &map, new_count)); - } - - /// perform the cube expansion for a given polycube - /// if perform extra expansions for cases where the dimensions are equal as - /// square sides may miss poly cubes otherwise - #[inline] - fn expand_cube_map(&mut self, seed: &CubeMapPos, shape: &Dim, count: usize) { - use MatrixCol::*; - - if shape.x == shape.y && shape.x > 0 { - let rotz = rot_matrix_points(seed, shape, count, YN, XN, ZN, 1025); - self.do_cube_expansion(&rotz, shape, count); - } - - if shape.y == shape.z && shape.y > 0 { - let rotx = rot_matrix_points(seed, shape, count, XN, ZP, YP, 1025); - self.do_cube_expansion(&rotx, shape, count); - } - if shape.x == shape.z && shape.x > 0 { - let roty = rot_matrix_points(seed, shape, count, ZP, YP, XN, 1025); - self.do_cube_expansion(&roty, shape, count); - } - - self.do_cube_expansion(seed, shape, count); - } - /// Calculate the amount of canonical children of size `target` /// that polycube `seed` of size `count` has. /// @@ -179,21 +37,33 @@ impl MapStore { count: usize, target: usize, ) -> usize { - let mut map = Self::new(); + let mut store = Self::new(); let shape = seed.extrapolate_dim(); - let seed = to_min_rot_points(seed, &shape, count); + let seed = seed.to_min_rot_points(shape, count); let shape = seed.extrapolate_dim(); - - map.expand_cube_map(&seed, &shape, count); - - map.inner + let meta = PointListMeta { + point_list: seed, + dim: shape, + count: count, + }; + meta.unique_expansions().for_each( + |PointListMeta { + point_list: map, + dim, + count, + }| store.insert_map(dim, map, count), + ); + + store + .inner .retain(|child| child.is_canonical_root(count, &seed)); if count + 1 == target { - map.inner.len() + store.inner.len() } else { - map.inner + store + .inner .iter() .map(|child| Self::enumerate_canonical_children_min_mem(child, count + 1, target)) .sum() diff --git a/rust/src/iterator.rs b/rust/src/iterator.rs index 4e8bde8..9865a6e 100644 --- a/rust/src/iterator.rs +++ b/rust/src/iterator.rs @@ -1,6 +1,6 @@ use std::collections::{HashMap, HashSet}; -use crate::pcube::RawPCube; +use crate::polycubes::pcube::RawPCube; /// An iterator over polycubes pub trait PolycubeIterator: Iterator diff --git a/rust/src/lib.rs b/rust/src/lib.rs index 2782875..9c269b9 100644 --- a/rust/src/lib.rs +++ b/rust/src/lib.rs @@ -2,12 +2,8 @@ mod test; pub mod iterator; -pub mod pcube; - -pub mod naive_polycube; +pub mod polycubes; pub mod hashless; pub mod pointlist; -pub mod polycube_reps; pub mod rotation_reduced; -pub mod rotations; diff --git a/rust/src/pointlist.rs b/rust/src/pointlist.rs index 5de9030..b260f99 100644 --- a/rust/src/pointlist.rs +++ b/rust/src/pointlist.rs @@ -1,309 +1,157 @@ -use std::{cmp::max, time::Instant}; +use std::time::Instant; -use crate::{ - polycube_reps::{CubeMapPos, CubeMapPosPart, Dim}, - rotations::{rot_matrix_points, to_min_rot_points, MatrixCol}, +use crate::polycubes::{ + pcube::RawPCube, + point_list::{CubeMapPos, PointListMeta}, + Dim, }; -use crate::pcube::RawPCube; use hashbrown::{HashMap, HashSet}; use indicatif::ProgressBar; use parking_lot::RwLock; use rayon::prelude::*; -///structure to store the polycubes -/// stores the key and fist block index as a key -/// to the set of 15 block tails that correspond to that shape and start -/// used for reducing mutex preasure on insertion -/// used as buckets for parallelising -/// however both of these give suboptomal performance due to the uneven distribution -type MapStore = HashMap<(Dim, u16), RwLock>>; - -/// helper function to not duplicate code for canonicalising polycubes -/// and storing them in the hashset -fn insert_map(store: &MapStore, dim: &Dim, map: &CubeMapPos<16>, count: usize) { - let map = to_min_rot_points(map, dim, count); - let mut body = CubeMapPosPart { cubes: [0; 15] }; - for i in 1..count { - body.cubes[i - 1] = map.cubes[i]; - } - match store.get(&(*dim, map.cubes[0])) { - Some(map) => { - map.write().insert(body); - } - None => { - panic!( - "shape {:?} data {} {:?} count {}", - dim, map.cubes[0], body, count - ); +/// Structure to store sets of polycubes +pub struct MapStore { + /// Stores the shape and fist block index as a key + /// to the set of 15 block tails that correspond to that shape and start. + /// used for reducing rwlock pressure on insertion + /// used as buckets for parallelising + /// however both of these give suboptomal performance due to the uneven distribution + inner: HashMap<(Dim, u16), RwLock>>>, +} + +impl MapStore { + pub fn new() -> Self { + Self { + inner: HashMap::new(), } } -} -///linearly scan backwards to insertion point overwrites end of slice -#[inline] -pub fn array_insert(val: u16, arr: &mut [u16]) { - for i in 1..(arr.len()) { - if arr[arr.len() - 1 - i] > val { - arr[arr.len() - i] = arr[arr.len() - 1 - i]; - } else { - arr[arr.len() - i] = val; + fn insert(&self, plm: &PointListMeta<16>) { + let map = &plm.point_list; + let dim = plm.dim; + let count = plm.count; + + // Check if we don't already happen to be in the minimum rotation position. + let mut body_maybemin = CubeMapPos::new(); + body_maybemin.cubes[0..count].copy_from_slice(&map.cubes[1..count + 1]); + let dim_maybe = map.extrapolate_dim(); + + // Weirdly enough, doing the copy and doing the lookup check this + // way is faster than only copying if `inner` has en entry for + // dim_maybe. + if self + .inner + .get(&(dim_maybe, map.cubes[0])) + .map(|v| v.read().contains(&body_maybemin)) + == Some(true) + { return; } - } - arr[0] = val; -} -/// moves contents of slice to index x+1, x==0 remains -#[inline] -pub fn array_shift(arr: &mut [u16]) { - for i in 1..(arr.len()) { - arr[arr.len() - i] = arr[arr.len() - 1 - i]; - } -} + let map = map.to_min_rot_points(dim, count); -/// try expaning each cube into both x+1 and x-1, calculating new dimension -/// and ensuring x is never negative -#[inline] -fn expand_xs(dst: &MapStore, seed: &CubeMapPos<16>, shape: &Dim, count: usize) { - for (i, coord) in seed.cubes[0..count].iter().enumerate() { - if !seed.cubes[(i + 1)..count].contains(&(coord + 1)) { - let mut new_shape = *shape; - let mut exp_map = *seed; - - array_insert(coord + 1, &mut exp_map.cubes[i..=count]); - new_shape.x = max(new_shape.x, ((coord + 1) & 0x1f) as usize); - insert_map(dst, &new_shape, &exp_map, count + 1) - } - if coord & 0x1f != 0 { - if !seed.cubes[0..i].contains(&(coord - 1)) { - let mut exp_map = *seed; - //faster move of top half hopefully - array_shift(&mut exp_map.cubes[i..=count]); - array_insert(coord - 1, &mut exp_map.cubes[0..=i]); - insert_map(dst, shape, &exp_map, count + 1) - } - } else { - let mut new_shape = *shape; - new_shape.x += 1; - let mut exp_map = *seed; - for i in 0..count { - exp_map.cubes[i] += 1; - } - array_shift(&mut exp_map.cubes[i..=count]); - array_insert(*coord, &mut exp_map.cubes[0..=i]); - insert_map(dst, &new_shape, &exp_map, count + 1) - } - } -} + let mut body = CubeMapPos::new(); + body.cubes[0..count].copy_from_slice(&map.cubes[1..count + 1]); -/// try expaning each cube into both y+1 and y-1, calculating new dimension -/// and ensuring y is never negative -#[inline] -fn expand_ys(dst: &MapStore, seed: &CubeMapPos<16>, shape: &Dim, count: usize) { - for (i, coord) in seed.cubes[0..count].iter().enumerate() { - if !seed.cubes[(i + 1)..count].contains(&(coord + (1 << 5))) { - let mut new_shape = *shape; - let mut exp_map = *seed; - array_insert(coord + (1 << 5), &mut exp_map.cubes[i..=count]); - new_shape.y = max(new_shape.y, (((coord >> 5) + 1) & 0x1f) as usize); - insert_map(dst, &new_shape, &exp_map, count + 1) - } - if (coord >> 5) & 0x1f != 0 { - if !seed.cubes[0..i].contains(&(coord - (1 << 5))) { - let mut exp_map = *seed; - //faster move of top half hopefully - array_shift(&mut exp_map.cubes[i..=count]); - array_insert(coord - (1 << 5), &mut exp_map.cubes[0..=i]); - insert_map(dst, shape, &exp_map, count + 1) - } - } else { - let mut new_shape = *shape; - new_shape.y += 1; - let mut exp_map = *seed; - for i in 0..count { - exp_map.cubes[i] += 1 << 5; - } - array_shift(&mut exp_map.cubes[i..=count]); - array_insert(*coord, &mut exp_map.cubes[0..=i]); - insert_map(dst, &new_shape, &exp_map, count + 1) - } - } -} + let entry = self + .inner + .get(&(plm.dim, map.cubes[0])) + .expect("Cube size does not have entry in destination map"); -/// try expaning each cube into both z+1 and z-1, calculating new dimension -/// and ensuring z is never negative -#[inline] -fn expand_zs(dst: &MapStore, seed: &CubeMapPos<16>, shape: &Dim, count: usize) { - for (i, coord) in seed.cubes[0..count].iter().enumerate() { - if !seed.cubes[(i + 1)..count].contains(&(coord + (1 << 10))) { - let mut new_shape = *shape; - let mut exp_map = *seed; - array_insert(coord + (1 << 10), &mut exp_map.cubes[i..=count]); - new_shape.z = max(new_shape.z, (((coord >> 10) + 1) & 0x1f) as usize); - insert_map(dst, &new_shape, &exp_map, count + 1) - } - if (coord >> 10) & 0x1f != 0 { - if !seed.cubes[0..i].contains(&(coord - (1 << 10))) { - let mut exp_map = *seed; - //faster move of top half hopefully - array_shift(&mut exp_map.cubes[i..=count]); - array_insert(coord - (1 << 10), &mut exp_map.cubes[0..=i]); - insert_map(dst, shape, &exp_map, count + 1) - } - } else { - let mut new_shape = *shape; - new_shape.z += 1; - let mut exp_map = *seed; - for i in 0..count { - exp_map.cubes[i] += 1 << 10; - } - array_shift(&mut exp_map.cubes[i..=count]); - array_insert(*coord, &mut exp_map.cubes[0..=i]); - insert_map(dst, &new_shape, &exp_map, count + 1) - } + entry.write().insert(body); } -} -/// reduce number of expansions needing to be performed based on -/// X >= Y >= Z constraint on Dim -#[inline] -fn do_cube_expansion(dst: &MapStore, seed: &CubeMapPos<16>, shape: &Dim, count: usize) { - if shape.y < shape.x { - expand_ys(dst, seed, shape, count); - } - if shape.z < shape.y { - expand_zs(dst, seed, shape, count); - } - expand_xs(dst, seed, shape, count); -} + /// helper for inner_exp in expand_cube_set it didnt like going directly in the closure + fn expand_cube_sub_set( + &self, + shape: Dim, + first_cube: u16, + body: impl Iterator>, + count: usize, + ) { + let mut seed = CubeMapPos { + cubes: [first_cube, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + }; -/// perform the cube expansion for a given polycube -/// if perform extra expansions for cases where the dimensions are equal as -/// square sides may miss poly cubes otherwise -#[inline] -fn expand_cube_map(dst: &MapStore, seed: &CubeMapPos<16>, shape: &Dim, count: usize) { - if shape.x == shape.y && shape.x > 0 { - let rotz = rot_matrix_points( - seed, - shape, - count, - MatrixCol::YN, - MatrixCol::XN, - MatrixCol::ZN, - 1025, - ); - do_cube_expansion(dst, &rotz, shape, count); - } - if shape.y == shape.z && shape.y > 0 { - let rotx = rot_matrix_points( - seed, - shape, - count, - MatrixCol::XN, - MatrixCol::ZP, - MatrixCol::YP, - 1025, - ); - do_cube_expansion(dst, &rotx, shape, count); - } - if shape.x == shape.z && shape.x > 0 { - let roty = rot_matrix_points( - seed, - shape, - count, - MatrixCol::ZP, - MatrixCol::YP, - MatrixCol::XN, - 1025, - ); - do_cube_expansion(dst, &roty, shape, count); - } - do_cube_expansion(dst, seed, shape, count); -} + for seed_body in body { + for i in 1..count { + seed.cubes[i] = seed_body.cubes[i - 1]; + } -/// helper for inner_exp in expand_cube_set it didnt like going directly in the closure -fn expand_cube_sub_set( - (shape, start): &(Dim, u16), - body: &RwLock>, - count: usize, - dst: &MapStore, -) { - let mut seed = CubeMapPos { - cubes: [*start, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], - }; - for seed_body in body.read().iter() { - for i in 1..count { - seed.cubes[i] = seed_body.cubes[i - 1]; + // body.cubes.copy_within(0..body.cubes.len() - 1, 1); + let seed_meta = PointListMeta { + point_list: seed, + dim: shape, + count, + }; + seed_meta.expand().for_each(|plm| self.insert(&plm)); } - expand_cube_map(dst, &seed, &shape, count); } -} -fn expand_cube_set( - seeds: &MapStore, - count: usize, - dst: &mut MapStore, - bar: &ProgressBar, - parallel: bool, -) { - // set up the dst sets before starting parallel processing so accessing doesnt block a global mutex - for x in 0..=count + 1 { - for y in 0..=(count + 1) / 2 { - for z in 0..=(count + 1) / 3 { - for i in 0..(y + 1) * 32 { - dst.insert((Dim { x, y, z }, i as u16), RwLock::new(HashSet::new())); + fn expand_cube_set(self, count: usize, dst: &mut MapStore, bar: &ProgressBar, parallel: bool) { + // set up the dst sets before starting parallel processing so accessing doesnt block a global mutex + for x in 0..=count + 1 { + for y in 0..=(count + 1) / 2 { + for z in 0..=(count + 1) / 3 { + for i in 0..(y + 1) * 32 { + dst.inner + .insert((Dim { x, y, z }, i as u16), RwLock::new(HashSet::new())); + } } } } - } - bar.set_message(format!("seed subsets expanded for N = {}...", count + 1)); + bar.set_position(0); + bar.set_length(self.inner.len() as u64); + bar.set_message(format!("seed subsets expanded for N = {}...", count + 1)); - let inner_exp = |(ss, body)| { - expand_cube_sub_set(ss, body, count, dst); - bar.inc(1); - }; + let inner_exp = |((shape, first_cube), body): (_, RwLock>)| { + dst.expand_cube_sub_set(shape, first_cube, body.into_inner().into_iter(), count); + bar.inc(1); + }; - //use parallel iterator or not to run expand_cube_set - if parallel { - seeds.par_iter().for_each(inner_exp); - } else { - seeds.iter().for_each(inner_exp); - } - //retain only subsets that have polycubes - dst.retain(|_, v| v.read().len() > 0); -} + // Use parallel iterator or not to run expand_cube_set + if parallel { + self.inner.into_par_iter().for_each(inner_exp); + } else { + self.inner.into_iter().for_each(inner_exp); + } -/// count the number of polycubes across all subsets -fn count_polycubes(maps: &MapStore) -> usize { - let mut total = 0; - #[cfg(feature = "diagnostics")] - for ((d, s), body) in maps.iter().rev() { - println!( - "({}, {}, {}) {} _> {}", - d.x + 1, - d.y + 1, - d.z + 1, - s, - body.len() - ); - } - for (_, body) in maps.iter() { - total += body.read().len() + //retain only subsets that have polycubes + dst.inner.retain(|_, v| v.read().len() > 0); + } + + /// Count the number of polycubes across all subsets + fn count_polycubes(&self) -> usize { + let mut total = 0; + #[cfg(feature = "diagnostics")] + for ((d, s), body) in maps.iter().rev() { + println!( + "({}, {}, {}) {} _> {}", + d.x + 1, + d.y + 1, + d.z + 1, + s, + body.len() + ); + } + for (_, body) in self.inner.iter() { + total += body.read().len() + } + + total } - total -} -/// distructively move the data from hashset to vector -fn move_polycubes_to_vec(maps: &mut MapStore) -> Vec> { - let mut v = Vec::new(); - while let Some(((dim, head), body)) = maps.iter().next() { - //extra scope to free lock and make borrow checker allow mutation of maps - { + /// Destructively move the data from hashset to vector + pub fn into_vec(self) -> Vec> { + let mut v = Vec::with_capacity(self.count_polycubes()); + + for ((_, head), body) in self.inner.into_iter() { let bod = body.read(); let mut cmp = CubeMapPos::new(); - cmp.cubes[0] = *head; + cmp.cubes[0] = head; for b in bod.iter() { for i in 0..15 { cmp.cubes[i + 1] = b.cubes[i]; @@ -311,20 +159,15 @@ fn move_polycubes_to_vec(maps: &mut MapStore) -> Vec> { v.push(cmp); } } - let dim = *dim; - let head = *head; - maps.remove(&(dim, head)); + + v } - v -} -/// distructively move the data from hashset to vector -fn _clone_polycubes_to_vec(maps: &mut MapStore) -> Vec> { - let mut v = Vec::new(); + /// Copy the data from hashset to vector + pub fn to_vec(&self) -> Vec> { + let mut v = Vec::with_capacity(self.count_polycubes()); - for ((_, head), body) in maps.iter() { - //extra scope to free lock and make borrow checker allow mutation of maps - { + for ((_, head), body) in self.inner.iter() { let bod = body.read(); let mut cmp = CubeMapPos::new(); cmp.cubes[0] = *head; @@ -335,8 +178,9 @@ fn _clone_polycubes_to_vec(maps: &mut MapStore) -> Vec> { v.push(cmp); } } + + v } - v } /// run pointlist based generation algorithm @@ -355,66 +199,39 @@ pub fn gen_polycubes( for seed in current.iter() { let seed: CubeMapPos<16> = seed.into(); let dim = seed.extrapolate_dim(); - if !seeds.contains_key(&(dim, seed.cubes[0])) { + if !seeds.inner.contains_key(&(dim, seed.cubes[0])) { for i in 0..(dim.y * 32 + dim.x + 1) { - seeds.insert((dim, i as u16), RwLock::new(HashSet::new())); + seeds + .inner + .insert((dim, i as u16), RwLock::new(HashSet::new())); } } - insert_map(&seeds, &dim, &seed, calculate_from - 1); + let seed_meta = PointListMeta { + point_list: seed, + dim, + count: calculate_from - 1, + }; + seeds.insert(&seed_meta); } drop(current); for i in calculate_from..=n as usize { bar.set_message(format!("seed subsets expanded for N = {}...", i)); let mut dst = MapStore::new(); - expand_cube_set(&mut seeds, i - 1, &mut dst, bar, parallel); + seeds.expand_cube_set(i - 1, &mut dst, bar, parallel); seeds = dst; - // if use_cache && i < n { - // let next = clone_polycubes_to_vec(&mut seeds); - // let name = &format!("cubes_{i}.pcube"); - // // if !std::fs::File::open(name).is_ok() { - // // println!("Saving {} cubes to cache file", next.len()); - // // PCubeFile::write_file( - // // false, - // // compression.into(), - // // next.iter().map(|v| v.into()), - // // name, - // // ) - // // .unwrap(); - // // } else { - // // println!("Cache file already exists for N = {i}. Not overwriting."); - // // } - // } - let t1_stop = Instant::now(); let time = t1_stop.duration_since(t1_start).as_micros(); bar.set_message(format!( - "Found {} unique expansions (N = {i}) in {}.{:06}s", - count_polycubes(&seeds), + "Found {} unique expansions (N = {i}) in {}.{:06}s", + seeds.count_polycubes(), time / 1000000, time % 1000000 )); bar.finish(); } - // exported eperately for memory concerns. already quite a lot more probably but not much I can do - let next = move_polycubes_to_vec(&mut seeds); - // if use_cache { - // let name = &format!("cubes_{n}.pcube"); - // if !std::fs::File::open(name).is_ok() { - // println!("Saving {} cubes to cache file", next.len()); - // PCubeFile::write_file( - // false, - // compression.into(), - // next.iter().map(|v| v.into()), - // name, - // ) - // .unwrap(); - // } else { - // println!("Cache file already exists for N = {n}. Not overwriting."); - // } - // } - next - //count_polycubes(&seeds); + + seeds.into_vec() } diff --git a/rust/src/polycube_reps.rs b/rust/src/polycube_reps.rs deleted file mode 100644 index ac82979..0000000 --- a/rust/src/polycube_reps.rs +++ /dev/null @@ -1,520 +0,0 @@ -use std::cmp::{max, min}; - -use crate::pcube::RawPCube; - -use crate::rotations::{map_coord, to_min_rot_points, MatrixCol::*}; - -/// Polycube representation -/// stores up to 16 blocks (number of cubes normally implicit or seperate in program state) -/// well formed polycubes are a sorted list of coordinates low to high -/// cordinates are group of packed 5 bit unsigned integers '-ZZZZZYYYYYXXXXX' -#[derive(PartialEq, Eq, Hash, Clone, Copy, Debug, PartialOrd, Ord)] -pub struct CubeMapPos { - pub cubes: [u16; N], -} - -/// Conversion function to assist with loading Polycubes from cache file to point-list implementation -/// Returned cubes may not be fully canonicalized (X >= Y >= Z guarenteed but not exact rotation) -impl From<&RawPCube> for CubeMapPos { - fn from(src: &RawPCube) -> Self { - let mut dst = CubeMapPos { cubes: [0; N] }; - let (x, y, z) = src.dims(); - let x = x as usize; - let y = y as usize; - let z = z as usize; - let dim = Dim { - x: x as usize - 1, - y: y as usize - 1, - z: z as usize - 1, - }; - - //correction matrix to convert to canonical dimension. I dont like it but it works - let (x_col, y_col, z_col, rdim) = if x >= y && y >= z { - (XP, YP, ZP, dim) - } else if x >= z && z >= y { - ( - XP, - ZP, - YN, - Dim { - x: x - 1, - y: z - 1, - z: y - 1, - }, - ) - } else if y >= x && x >= z { - ( - YP, - XP, - ZN, - Dim { - x: y - 1, - y: x - 1, - z: z - 1, - }, - ) - } else if y >= z && z >= x { - ( - YP, - ZP, - XP, - Dim { - x: y - 1, - y: z - 1, - z: x - 1, - }, - ) - } else if z >= x && x >= y { - ( - ZN, - XP, - YN, - Dim { - x: z - 1, - y: x - 1, - z: y - 1, - }, - ) - } else if z >= y && y >= x { - ( - ZN, - YN, - XP, - Dim { - x: z - 1, - y: y - 1, - z: x - 1, - }, - ) - } else { - panic!("imposible dimension of shape {:?}", dim) - }; - - let mut dst_index = 0; - for dz in 0..z as u16 { - for dy in 0..y as u16 { - for dx in 0..x as u16 { - if src.get(dx as u8, dy as u8, dz as u8) { - let cx = map_coord(dx, dy, dz, &dim, x_col); - let cy = map_coord(dx, dy, dz, &dim, y_col); - let cz = map_coord(dx, dy, dz, &dim, z_col); - if cx > rdim.x as u16 || cy > rdim.y as u16 || cz > rdim.z as u16 { - panic!("illegal block place {}, {}, {} {:?}", cx, cy, cz, dim) - } - let pack = ((cz << 10) | (cy << 5) | cx) as u16; - dst.cubes[dst_index] = pack; - dst_index += 1; - } - } - } - } - dst - } -} - -impl From<&'_ CubeMapPos> for RawPCube { - fn from(src: &'_ CubeMapPos) -> Self { - //cube is sorted numerically and then has trailing zeros - let count = src.extrapolate_count(); - let dim = src.extrapolate_dim(); - - let mut dst = RawPCube::new_empty(dim.x as u8 + 1, dim.y as u8 + 1, dim.z as u8 + 1); - for p in src.cubes[0..count].iter() { - let ix = *p & 0x1f; - let iy = (*p >> 5) & 0x1f; - let iz = (*p >> 10) & 0x1f; - dst.set(ix as u8, iy as u8, iz as u8, true); - } - dst - } -} - -impl From for CubeMapPos { - fn from(value: RawPCube) -> Self { - (&value).into() - } -} - -impl From> for RawPCube { - fn from(value: CubeMapPos) -> Self { - (&value).into() - } -} - -impl CubeMapPos { - pub fn new() -> Self { - CubeMapPos { cubes: [0; N] } - } - pub fn extrapolate_count(&self) -> usize { - let mut count = 1; - while self.cubes[count] > self.cubes[count - 1] { - count += 1; - } - count - } - pub fn extrapolate_dim(&self) -> Dim { - let count = self.extrapolate_count(); - let mut dim = Dim { x: 0, y: 0, z: 0 }; - for p in self.cubes[0..count].iter() { - let ix = *p & 0x1f; - let iy = (*p >> 5) & 0x1f; - let iz = (*p >> 10) & 0x1f; - dim.x = max(dim.x, ix as usize); - dim.y = max(dim.y, iy as usize); - dim.z = max(dim.z, iz as usize); - } - dim - } - - fn is_continuous(&self, len: usize) -> bool { - let start = self.cubes[0]; - let mut polycube2 = [start; 32]; - for i in 1..len { - polycube2[i] = self.cubes[i]; - } - let polycube = polycube2; - //sets were actually slower even when no allocating - let mut to_explore = [start; 32]; - let mut exp_head = 1; - let mut exp_tail = 0; - //to_explore[0] = start; - while exp_head > exp_tail { - let p = to_explore[exp_tail]; - exp_tail += 1; - if p & 0x1f != 0 && !to_explore.contains(&(p - 1)) && polycube.contains(&(p - 1)) { - to_explore[exp_head] = p - 1; - // unsafe {*to_explore.get_unchecked_mut(exp_head) = p - 1;} - exp_head += 1; - } - if p & 0x1f != 0x1f && !to_explore.contains(&(p + 1)) && polycube.contains(&(p + 1)) { - to_explore[exp_head] = p + 1; - // unsafe {*to_explore.get_unchecked_mut(exp_head) = p + 1;} - exp_head += 1; - } - if (p >> 5) & 0x1f != 0 - && !to_explore.contains(&(p - (1 << 5))) - && polycube.contains(&(p - (1 << 5))) - { - to_explore[exp_head] = p - (1 << 5); - // unsafe {*to_explore.get_unchecked_mut(exp_head) = p - (1 << 5);} - exp_head += 1; - } - if (p >> 5) & 0x1f != 0x1f - && !to_explore.contains(&(p + (1 << 5))) - && polycube.contains(&(p + (1 << 5))) - { - to_explore[exp_head] = p + (1 << 5); - // unsafe {*to_explore.get_unchecked_mut(exp_head) = p + (1 << 5);} - exp_head += 1; - } - if (p >> 10) & 0x1f != 0 - && !to_explore.contains(&(p - (1 << 10))) - && polycube.contains(&(p - (1 << 10))) - { - to_explore[exp_head] = p - (1 << 10); - // unsafe {*to_explore.get_unchecked_mut(exp_head) = p - (1 << 10);} - exp_head += 1; - } - if (p >> 10) & 0x1f != 0x1f - && !to_explore.contains(&(p + (1 << 10))) - && polycube.contains(&(p + (1 << 10))) - { - to_explore[exp_head] = p + (1 << 10); - // unsafe {*to_explore.get_unchecked_mut(exp_head) = p + (1 << 10);} - exp_head += 1; - } - } - exp_head == len - } - - fn renormalize(&self, dim: &Dim, count: usize) -> (Self, Dim) { - let mut dst = CubeMapPos::new(); - let x = dim.x; - let y = dim.y; - let z = dim.z; - let (x_col, y_col, z_col, rdim) = if x >= y && y >= z { - (XP, YP, ZP, Dim { x: x, y: y, z: z }) - } else if x >= z && z >= y { - (XP, ZP, YN, Dim { x: x, y: z, z: y }) - } else if y >= x && x >= z { - (YP, XP, ZN, Dim { x: y, y: x, z: z }) - } else if y >= z && z >= x { - (YP, ZP, XP, Dim { x: y, y: z, z: x }) - } else if z >= x && x >= y { - (ZN, XP, YN, Dim { x: z, y: x, z: y }) - } else if z >= y && y >= x { - (ZN, YN, XP, Dim { x: z, y: y, z: x }) - } else { - panic!("imposible dimension of shape {:?}", dim) - }; - for (i, d) in self.cubes[0..count].iter().enumerate() { - let dx = d & 0x1f; - let dy = (d >> 5) & 0x1f; - let dz = (d >> 10) & 0x1f; - let cx = map_coord(dx, dy, dz, &dim, x_col); - let cy = map_coord(dx, dy, dz, &dim, y_col); - let cz = map_coord(dx, dy, dz, &dim, z_col); - let pack = ((cz << 10) | (cy << 5) | cx) as u16; - dst.cubes[i] = pack; - // unsafe {*dst.cubes.get_unchecked_mut(i) = pack;} - } - //dst.cubes.sort(); - (dst, rdim) - } - - fn remove_cube(&self, point: usize, count: usize) -> (Self, Dim) { - let mut min_corner = Dim { - x: 0x1f, - y: 0x1f, - z: 0x1f, - }; - let mut max_corner = Dim { x: 0, y: 0, z: 0 }; - let mut root_candidate = CubeMapPos::new(); - let mut candidate_ptr = 0; - for i in 0..=count { - if i != point { - let pos = self.cubes[i]; - // let pos = unsafe {*exp.cubes.get_unchecked(i)}; - let x = pos as usize & 0x1f; - let y = (pos as usize >> 5) & 0x1f; - let z = (pos as usize >> 10) & 0x1f; - min_corner.x = min(min_corner.x, x); - min_corner.y = min(min_corner.y, y); - min_corner.z = min(min_corner.z, z); - max_corner.x = max(max_corner.x, x); - max_corner.y = max(max_corner.y, y); - max_corner.z = max(max_corner.z, z); - root_candidate.cubes[candidate_ptr] = pos; - // unsafe {*root_candidate.cubes.get_unchecked_mut(candidate_ptr) = pos;} - candidate_ptr += 1; - } - } - let offset = (min_corner.z << 10) | (min_corner.y << 5) | min_corner.x; - for i in 0..count { - root_candidate.cubes[i] -= offset as u16; - } - max_corner.x = max_corner.x - min_corner.x; - max_corner.y = max_corner.y - min_corner.y; - max_corner.z = max_corner.z - min_corner.z; - (root_candidate, max_corner) - } - - pub fn is_canonical_root(&self, count: usize, seed: &Self) -> bool { - for sub_cube_id in 0..=count { - let (mut root_candidate, mut dim) = self.remove_cube(sub_cube_id, count); - if !root_candidate.is_continuous(count) { - continue; - } - if dim.x < dim.y || dim.y < dim.z || dim.x < dim.z { - let (rroot_candidate, rdim) = root_candidate.renormalize(&dim, count); - root_candidate = rroot_candidate; - dim = rdim; - root_candidate.cubes[0..count].sort_unstable(); - } - let mrp = to_min_rot_points(&root_candidate, &dim, count); - if &mrp < seed { - return false; - } - } - true - } -} - -/// Partial Polycube representation for storage -/// the first block is stored seperately and used as a key to access the set -#[derive(PartialEq, Eq, Hash, Clone, Copy, Debug, PartialOrd, Ord)] -pub struct CubeMapPosPart { - pub cubes: [u16; 15], -} - -/// the "Dimension" or "Shape" of a poly cube -/// defines the maximum bounds of a polycube -/// X >= Y >= Z for efficiency reasons reducing the number of rotations needing to be performed -/// stores len() - 1 for each dimension so the unit cube has a size of (0, 0, 0) -/// and the 2x1x1 starting seed has a dimension of (1, 0, 0) -#[derive(PartialEq, Eq, Hash, Clone, Copy, Debug, PartialOrd, Ord)] -pub struct Dim { - pub x: usize, - pub y: usize, - pub z: usize, -} - -#[cfg(not(feature = "size16"))] -pub type CubeRow = u32; -#[cfg(feature = "size16")] -pub type CubeRow = u16; - -//CubeRow is an integer type either u16 or u32 based on flags -#[derive(PartialEq, Eq, Hash, Clone, Copy, Debug)] -pub struct CubeMap { - pub x: u32, //max x index (len(xs) - 1) - pub y: u32, //max y index (len(ys) - 1) - pub z: u32, //max z index (len(zs) - 1) - pub cube_map: [CubeRow; 6 * 6], -} - -impl PartialOrd for CubeMap { - fn partial_cmp(&self, other: &Self) -> Option { - match self.cube_map.partial_cmp(&other.cube_map) { - Some(core::cmp::Ordering::Equal) => {} - ord => return ord, - } - match self.x.partial_cmp(&other.x) { - Some(core::cmp::Ordering::Equal) => {} - ord => return ord, - } - match self.y.partial_cmp(&other.y) { - Some(core::cmp::Ordering::Equal) => {} - ord => return ord, - } - self.z.partial_cmp(&other.z) - } -} - -impl CubeMap { - /// returns 1 if it block at xyz exists - /// returns 0 if it doesnt - pub fn get_block(&self, x: usize, y: usize, z: usize) -> CubeRow { - (self.cube_map[z * (self.y as usize + 1) + y] >> x) & 1 - } - #[cfg(feature = "smallset")] - /// sets the block at xyz to exist - pub fn set_block(&mut self, x: usize, y: usize, z: usize) { - self.cube_map[z * (self.y as usize + 1) + y] |= 1 << x; - } - /// set a block to the bit v - /// IMORTANT: Sets, does not unset, performs an | on the vale will never clear even on set v = 0 - pub fn set_block_to(&mut self, x: usize, y: usize, z: usize, v: CubeRow) { - self.cube_map[z * (self.y as usize + 1) + y] |= v << x; - } - pub fn clear(&mut self) { - for i in 0..((self.z as usize + 1) * (self.y as usize + 1)) { - self.cube_map[i] = 0; - } - } - #[cfg(feature = "diagnostics")] - /// ensure expected number of cubes are set, only used as an integrity check - pub fn count_cubes(&self) -> usize { - let mut c = 0; - for i in 0..36 { - let mut x = self.cube_map[i]; - while x > 0 { - c += x as usize & 1; - x >>= 1; - } - } - c - } - #[cfg(feature = "diagnostics")] - /// ensure no blocks are set outside expected area, only used as an integrity check - pub fn validate_bounds(&self) -> bool { - for x in (self.x + 1)..MAX_X as u32 { - for y in 0..=self.y { - for z in 0..=self.z { - if self.get_block(x as usize, y as usize, z as usize) == 1 { - return false; - } - } - } - } - for i in ((self.z as usize + 1) * (self.y as usize + 1))..36 { - if self.cube_map[i] != 0 { - return false; - } - } - true - } - - #[cfg(feature = "diagnostics")] - /// find an existing block to seed continuity check - fn find_a_block(&self) -> Dim { - for y in 0..=self.y { - for z in 0..=self.z { - let mut x = 0; - let mut row = self.cube_map[(z * (self.y + 1) + y) as usize]; - if row != 0 { - while row > 0 { - if row & 1 == 1 { - return Dim { - x: x as usize, - y: y as usize, - z: z as usize, - }; - } - x += 1; - row >>= 1; - } - } - } - } - panic!("{:?} empty", self); - } - - #[cfg(feature = "diagnostics")] - /// ensure all blocks are connected, only used as an integrity check - pub fn validate_continuity(&self) -> bool { - let mut to_visit = HashSet::new(); - let mut map = *self; - let start = self.find_a_block(); - to_visit.insert(start); - while let Some(p) = to_visit.iter().next() { - let p = *p; - to_visit.remove(&p); - map.cube_map[p.z * (map.y as usize + 1) + p.y] &= !(1 << p.x); - if p.x > 0 && (map.cube_map[p.z * (map.y as usize + 1) + p.y] >> (p.x - 1)) & 1 == 1 { - to_visit.insert(Dim { - x: p.x - 1, - y: p.y, - z: p.z, - }); - } - if p.x < map.x as usize - && (map.cube_map[p.z * (map.y as usize + 1) + p.y] >> (p.x + 1)) & 1 == 1 - { - to_visit.insert(Dim { - x: p.x + 1, - y: p.y, - z: p.z, - }); - } - if p.y > 0 && (map.cube_map[p.z * (map.y as usize + 1) + (p.y - 1)] >> p.x) & 1 == 1 { - to_visit.insert(Dim { - x: p.x, - y: p.y - 1, - z: p.z, - }); - } - if p.y < map.y as usize - && (map.cube_map[p.z * (map.y as usize + 1) + (p.y + 1)] >> p.x) & 1 == 1 - { - to_visit.insert(Dim { - x: p.x, - y: p.y + 1, - z: p.z, - }); - } - if p.z > 0 && (map.cube_map[(p.z - 1) * (map.y as usize + 1) + p.y] >> p.x) & 1 == 1 { - to_visit.insert(Dim { - x: p.x, - y: p.y, - z: p.z - 1, - }); - } - if p.z < map.z as usize - && (map.cube_map[(p.z + 1) * (map.y as usize + 1) + p.y] >> p.x) & 1 == 1 - { - to_visit.insert(Dim { - x: p.x, - y: p.y, - z: p.z + 1, - }); - } - } - for i in 0..((map.y + 1) * (map.z + 1)) { - if map.cube_map[i as usize] != 0 { - return false; - } - } - true - } -} diff --git a/rust/src/polycubes/mod.rs b/rust/src/polycubes/mod.rs new file mode 100644 index 0000000..1039b80 --- /dev/null +++ b/rust/src/polycubes/mod.rs @@ -0,0 +1,34 @@ +// hopefully only expose pcube reps in the future not full modules +pub mod point_list; + +pub mod naive_polycube; + +pub mod pcube; +pub use self::pcube::RawPCube; + +pub mod rotation_reduced; + +/// the "Dimension" or "Shape" of a poly cube +/// defines the maximum bounds of a polycube +/// X >= Y >= Z for efficiency reasons reducing the number of rotations needing to be performed +/// stores len() for each dimension so the unit cube has a size of (1, 1, 1) +/// and the 2x1x1 starting seed has a dimension of (2, 1, 1) +#[derive(PartialEq, Eq, Hash, Clone, Copy, Debug, PartialOrd, Ord)] +pub struct Dim { + pub x: usize, + pub y: usize, + pub z: usize, +} + +pub trait PolyCube: From + Into + Sized { + type UniqueExpansionsIterator: Iterator; + + /// Produce an iterator that yields all unique n + 1 expansions of + /// `input`. + fn unique_expansions(&self) -> Self::UniqueExpansionsIterator; + + /// Return a copy of self in some "canonical" form. + fn canonical_form(&self) -> Self; + fn size(&self) -> usize; + fn dims(&self) -> Dim; +} diff --git a/rust/src/naive_polycube/expander.rs b/rust/src/polycubes/naive_polycube/expander.rs similarity index 100% rename from rust/src/naive_polycube/expander.rs rename to rust/src/polycubes/naive_polycube/expander.rs diff --git a/rust/src/naive_polycube/mod.rs b/rust/src/polycubes/naive_polycube/mod.rs similarity index 99% rename from rust/src/naive_polycube/mod.rs rename to rust/src/polycubes/naive_polycube/mod.rs index 1027907..bfcda47 100644 --- a/rust/src/naive_polycube/mod.rs +++ b/rust/src/polycubes/naive_polycube/mod.rs @@ -8,7 +8,7 @@ use crate::{ iterator::{ AllPolycubeIterator, AllUniquePolycubeIterator, PolycubeIterator, UniquePolycubeIterator, }, - pcube::RawPCube, + polycubes::pcube::RawPCube, }; mod expander; diff --git a/rust/src/naive_polycube/rotations.rs b/rust/src/polycubes/naive_polycube/rotations.rs similarity index 100% rename from rust/src/naive_polycube/rotations.rs rename to rust/src/polycubes/naive_polycube/rotations.rs diff --git a/rust/src/pcube/compression.rs b/rust/src/polycubes/pcube/compression.rs similarity index 100% rename from rust/src/pcube/compression.rs rename to rust/src/polycubes/pcube/compression.rs diff --git a/rust/src/pcube/mod.rs b/rust/src/polycubes/pcube/mod.rs similarity index 100% rename from rust/src/pcube/mod.rs rename to rust/src/polycubes/pcube/mod.rs diff --git a/rust/src/pcube/raw_pcube.rs b/rust/src/polycubes/pcube/raw_pcube.rs similarity index 100% rename from rust/src/pcube/raw_pcube.rs rename to rust/src/polycubes/pcube/raw_pcube.rs diff --git a/rust/src/polycubes/point_list/expand.rs b/rust/src/polycubes/point_list/expand.rs new file mode 100644 index 0000000..8b13789 --- /dev/null +++ b/rust/src/polycubes/point_list/expand.rs @@ -0,0 +1 @@ + diff --git a/rust/src/polycubes/point_list/mod.rs b/rust/src/polycubes/point_list/mod.rs new file mode 100644 index 0000000..a7a8686 --- /dev/null +++ b/rust/src/polycubes/point_list/mod.rs @@ -0,0 +1,594 @@ +pub mod expand; +pub mod rotate; + +use std::{ + cmp::{max, min}, + iter::{Chain, Flatten}, + option::IntoIter, +}; + +use super::{pcube::RawPCube, rotation_reduced::rotate::MatrixCol, Dim, PolyCube}; +use crate::polycubes::rotation_reduced::rotate::MatrixCol::*; + +/// Polycube representation +/// stores up to 16 blocks (number of cubes normally implicit or seperate in program state) +/// well formed polycubes are a sorted list of coordinates low to high +/// cordinates are group of packed 5 bit unsigned integers '-ZZZZZYYYYYXXXXX' +#[derive(PartialEq, Eq, Hash, Clone, Copy, Debug, PartialOrd, Ord)] +pub struct CubeMapPos { + pub cubes: [u16; N], +} + +/// Conversion function to assist with loading Polycubes from cache file to point-list implementation +/// Returned cubes may not be fully canonicalized (X >= Y >= Z guarenteed but not exact rotation) +impl From<&RawPCube> for CubeMapPos { + fn from(src: &RawPCube) -> Self { + let mut dst = CubeMapPos { cubes: [0; N] }; + let (x, y, z) = src.dims(); + let x = x as usize; + let y = y as usize; + let z = z as usize; + let dim = Dim { + x: x as usize - 1, + y: y as usize - 1, + z: z as usize - 1, + }; + + //correction matrix to convert to canonical dimension. I dont like it but it works + let (x_col, y_col, z_col, rdim) = if x >= y && y >= z { + (XP, YP, ZP, dim) + } else if x >= z && z >= y { + ( + XP, + ZP, + YN, + Dim { + x: x - 1, + y: z - 1, + z: y - 1, + }, + ) + } else if y >= x && x >= z { + ( + YP, + XP, + ZN, + Dim { + x: y - 1, + y: x - 1, + z: z - 1, + }, + ) + } else if y >= z && z >= x { + ( + YP, + ZP, + XP, + Dim { + x: y - 1, + y: z - 1, + z: x - 1, + }, + ) + } else if z >= x && x >= y { + ( + ZN, + XP, + YN, + Dim { + x: z - 1, + y: x - 1, + z: y - 1, + }, + ) + } else if z >= y && y >= x { + ( + ZN, + YN, + XP, + Dim { + x: z - 1, + y: y - 1, + z: x - 1, + }, + ) + } else { + panic!("imposible dimension of shape {:?}", dim) + }; + + let mut dst_index = 0; + for dz in 0..z as u16 { + for dy in 0..y as u16 { + for dx in 0..x as u16 { + if src.get(dx as u8, dy as u8, dz as u8) { + let cx = Self::map_coord(dx, dy, dz, &dim, x_col); + let cy = Self::map_coord(dx, dy, dz, &dim, y_col); + let cz = Self::map_coord(dx, dy, dz, &dim, z_col); + if cx > rdim.x as u16 || cy > rdim.y as u16 || cz > rdim.z as u16 { + panic!("illegal block place {}, {}, {} {:?}", cx, cy, cz, dim) + } + let pack = ((cz << 10) | (cy << 5) | cx) as u16; + dst.cubes[dst_index] = pack; + dst_index += 1; + } + } + } + } + dst + } +} + +impl From<&'_ CubeMapPos> for RawPCube { + fn from(src: &'_ CubeMapPos) -> Self { + //cube is sorted numerically and then has trailing zeros + let count = src.extrapolate_count(); + let dim = src.extrapolate_dim(); + + let mut dst = RawPCube::new_empty(dim.x as u8 + 1, dim.y as u8 + 1, dim.z as u8 + 1); + for p in src.cubes[0..count].iter() { + let ix = *p & 0x1f; + let iy = (*p >> 5) & 0x1f; + let iz = (*p >> 10) & 0x1f; + dst.set(ix as u8, iy as u8, iz as u8, true); + } + dst + } +} + +impl From for CubeMapPos { + fn from(value: RawPCube) -> Self { + (&value).into() + } +} + +impl From> for RawPCube { + fn from(value: CubeMapPos) -> Self { + (&value).into() + } +} + +/// Linearly scan backwards to insertion point overwrites end of slice +#[inline] +fn array_insert(val: u16, arr: &mut [u16]) { + for i in 1..(arr.len()) { + if arr[arr.len() - 1 - i] > val { + arr[arr.len() - i] = arr[arr.len() - 1 - i]; + } else { + arr[arr.len() - i] = val; + return; + } + } + arr[0] = val; +} + +/// Moves contents of slice to index x+1, x==0 remains +#[inline] +fn array_shift(arr: &mut [u16]) { + for i in 1..(arr.len()) { + arr[arr.len() - i] = arr[arr.len() - 1 - i]; + } +} + +impl CubeMapPos { + pub fn new() -> Self { + CubeMapPos { cubes: [0; N] } + } + + #[inline] + pub fn map_coord(x: u16, y: u16, z: u16, shape: &Dim, col: MatrixCol) -> u16 { + match col { + MatrixCol::XP => x, + MatrixCol::XN => shape.x as u16 - x, + MatrixCol::YP => y, + MatrixCol::YN => shape.y as u16 - y, + MatrixCol::ZP => z, + MatrixCol::ZN => shape.z as u16 - z, + } + } + + pub fn extrapolate_count(&self) -> usize { + let mut count = 1; + while self.cubes[count] > self.cubes[count - 1] { + count += 1; + } + count + } + + pub fn extrapolate_dim(&self) -> Dim { + let count = self.extrapolate_count(); + let mut dim = Dim { x: 0, y: 0, z: 0 }; + for p in self.cubes[0..count].iter() { + let ix = *p & 0x1f; + let iy = (*p >> 5) & 0x1f; + let iz = (*p >> 10) & 0x1f; + dim.x = max(dim.x, ix as usize); + dim.y = max(dim.y, iy as usize); + dim.z = max(dim.z, iz as usize); + } + dim + } + + fn is_continuous(&self, len: usize) -> bool { + let start = self.cubes[0]; + let mut polycube = [start; N]; + polycube[1..len].copy_from_slice(&self.cubes[1..len]); + + // sets were actually slower even when no allocating + let mut to_explore = [start; N]; + let mut exp_head = 1; + let mut exp_tail = 0; + + while exp_head > exp_tail { + let p = to_explore[exp_tail]; + exp_tail += 1; + + macro_rules! check { + ($([$check:ident, $shift:literal, $op:tt, $ignore_val:literal],)*) => { + $( + let value = (p $op (1 << $shift)); + if (p >> $shift) & 0x1F != $ignore_val + && !to_explore.contains(&value) + && polycube.contains(&value) + { + to_explore[exp_head] = value; + exp_head += 1; + } + )* + }; + } + + check!( + [x_min, 0, -, 0], + [x_plus, 0, +, 0x1f], + [y_min, 5, -, 0], + [y_plus, 5, +, 0x1f], + [z_min, 10, -, 0], + [z_plus, 10, +, 0x1f], + ); + } + exp_head == len + } + + fn renormalize(&self, dim: &Dim, count: usize) -> (Self, Dim) { + let mut dst = CubeMapPos::new(); + let x = dim.x; + let y = dim.y; + let z = dim.z; + let (x_col, y_col, z_col, rdim) = if x >= y && y >= z { + (XP, YP, ZP, Dim { x: x, y: y, z: z }) + } else if x >= z && z >= y { + (XP, ZP, YN, Dim { x: x, y: z, z: y }) + } else if y >= x && x >= z { + (YP, XP, ZN, Dim { x: y, y: x, z: z }) + } else if y >= z && z >= x { + (YP, ZP, XP, Dim { x: y, y: z, z: x }) + } else if z >= x && x >= y { + (ZN, XP, YN, Dim { x: z, y: x, z: y }) + } else if z >= y && y >= x { + (ZN, YN, XP, Dim { x: z, y: y, z: x }) + } else { + panic!("imposible dimension of shape {:?}", dim) + }; + for (i, d) in self.cubes[0..count].iter().enumerate() { + let dx = d & 0x1f; + let dy = (d >> 5) & 0x1f; + let dz = (d >> 10) & 0x1f; + let cx = Self::map_coord(dx, dy, dz, &dim, x_col); + let cy = Self::map_coord(dx, dy, dz, &dim, y_col); + let cz = Self::map_coord(dx, dy, dz, &dim, z_col); + let pack = ((cz << 10) | (cy << 5) | cx) as u16; + dst.cubes[i] = pack; + // unsafe {*dst.cubes.get_unchecked_mut(i) = pack;} + } + //dst.cubes.sort(); + (dst, rdim) + } + + fn remove_cube(&self, point: usize, count: usize) -> (Self, Dim) { + let mut min_corner = Dim { + x: 0x1f, + y: 0x1f, + z: 0x1f, + }; + let mut max_corner = Dim { x: 0, y: 0, z: 0 }; + let mut root_candidate = CubeMapPos::new(); + let mut candidate_ptr = 0; + for i in 0..=count { + if i != point { + let pos = self.cubes[i]; + // let pos = unsafe {*exp.cubes.get_unchecked(i)}; + let x = pos as usize & 0x1f; + let y = (pos as usize >> 5) & 0x1f; + let z = (pos as usize >> 10) & 0x1f; + min_corner.x = min(min_corner.x, x); + min_corner.y = min(min_corner.y, y); + min_corner.z = min(min_corner.z, z); + max_corner.x = max(max_corner.x, x); + max_corner.y = max(max_corner.y, y); + max_corner.z = max(max_corner.z, z); + root_candidate.cubes[candidate_ptr] = pos; + // unsafe {*root_candidate.cubes.get_unchecked_mut(candidate_ptr) = pos;} + candidate_ptr += 1; + } + } + let offset = (min_corner.z << 10) | (min_corner.y << 5) | min_corner.x; + for i in 0..count { + root_candidate.cubes[i] -= offset as u16; + } + max_corner.x = max_corner.x - min_corner.x; + max_corner.y = max_corner.y - min_corner.y; + max_corner.z = max_corner.z - min_corner.z; + (root_candidate, max_corner) + } + + pub fn is_canonical_root(&self, count: usize, seed: &Self) -> bool { + for sub_cube_id in 0..=count { + let (mut root_candidate, mut dim) = self.remove_cube(sub_cube_id, count); + if !root_candidate.is_continuous(count) { + continue; + } + if dim.x < dim.y || dim.y < dim.z || dim.x < dim.z { + let (rroot_candidate, rdim) = root_candidate.renormalize(&dim, count); + root_candidate = rroot_candidate; + dim = rdim; + root_candidate.cubes[0..count].sort_unstable(); + } + let mrp = root_candidate.to_min_rot_points(dim, count); + if &mrp < seed { + return false; + } + } + true + } +} + +#[derive(Copy, Clone)] +pub struct PointListMeta { + pub point_list: CubeMapPos, + pub dim: Dim, + pub count: usize, +} + +macro_rules! plm_expand { + ($iter_name:ident, $name:ident, $dim:ident, $shift:literal) => { + pub struct $iter_name { + inner: PointListMeta, + i: usize, + stored: Option>, + } + + impl<'a, const C: usize> Iterator for $iter_name { + type Item = PointListMeta; + + fn next(&mut self) -> Option { + loop { + if let Some(stored) = self.stored.take() { + return Some(stored); + } + + let i = self.i; + + if i == self.inner.count { + return None; + } + + self.i += 1; + let coord = *self.inner.point_list.cubes.get(i)?; + + let plus = coord + (1 << $shift); + let minus = coord - (1 << $shift); + + if !self.inner.point_list.cubes[(i + 1)..self.inner.count].contains(&plus) { + let mut new_shape = self.inner.dim; + let mut new_map = self.inner.point_list; + + array_insert(plus, &mut new_map.cubes[i..=self.inner.count]); + new_shape.$dim = + max(new_shape.$dim, (((coord >> $shift) + 1) & 0x1f) as usize); + + self.stored = Some(PointListMeta { + point_list: new_map, + dim: new_shape, + count: self.inner.count + 1, + }); + } + + let mut new_map = self.inner.point_list; + let mut new_shape = self.inner.dim; + + // If the coord is out of bounds for $dim, shift everything + // over and create the cube at the out-of-bounds position. + // If it is in bounds, check if the $dim - 1 value already + // exists. + let insert_coord = if (coord >> $shift) & 0x1f != 0 { + if !self.inner.point_list.cubes[0..i].contains(&minus) { + minus + } else { + continue; + } + } else { + new_shape.$dim += 1; + for i in 0..self.inner.count { + new_map.cubes[i] += 1 << $shift; + } + coord + }; + + array_shift(&mut new_map.cubes[i..=self.inner.count]); + array_insert(insert_coord, &mut new_map.cubes[0..=i]); + return Some(PointListMeta { + point_list: new_map, + dim: new_shape, + count: self.inner.count + 1, + }); + } + } + } + + impl PointListMeta { + #[inline(always)] + pub fn $name(self) -> $iter_name { + $iter_name { + inner: self, + i: 0, + stored: None, + } + } + } + }; +} + +plm_expand!(ExpandX, expand_x, x, 0); +plm_expand!(ExpandY, expand_y, y, 5); +plm_expand!(ExpandZ, expand_z, z, 10); + +pub type SubExpandIter = + Chain, Flatten>>>, Flatten>>>; + +pub type ExpandIter = Chain< + Chain< + Chain>>, Flatten>>>, + Flatten>>, + >, + SubExpandIter, +>; + +impl PointListMeta { + /// reduce number of expansions needing to be performed based on + /// X >= Y >= Z constraint on Dim + #[inline] + #[must_use] + fn do_expand(self) -> SubExpandIter { + let expand_ys = if self.dim.y < self.dim.x { + Some(self.expand_y()) + } else { + None + }; + + let expand_zs = if self.dim.z < self.dim.y { + Some(self.expand_z()) + } else { + None + }; + + let inner = self + .expand_x() + .chain(expand_ys.into_iter().flatten()) + .chain(expand_zs.into_iter().flatten()); + + inner + } + + /// perform the cube expansion for a given polycube + /// if perform extra expansions for cases where the dimensions are equal as + /// square sides may miss poly cubes otherwise + #[inline] + pub fn expand(&self) -> ExpandIter { + use MatrixCol::*; + + let z = if self.dim.x == self.dim.y && self.dim.x > 0 { + let rotz = self + .point_list + .rot_matrix_points(self.dim, self.count, YN, XN, ZN, 1025); + + let rotz = PointListMeta { + point_list: rotz, + dim: self.dim, + count: self.count, + }; + + Some(rotz.do_expand()) + } else { + None + }; + + let y = if self.dim.y == self.dim.z && self.dim.y > 0 { + let rotx = self + .point_list + .rot_matrix_points(self.dim, self.count, XN, ZP, YP, 1025); + + let rotx = PointListMeta { + point_list: rotx, + dim: self.dim, + count: self.count, + }; + + Some(rotx.do_expand()) + } else { + None + }; + + let x = if self.dim.x == self.dim.z && self.dim.x > 0 { + let roty = self + .point_list + .rot_matrix_points(self.dim, self.count, ZP, YP, XN, 1025); + + let roty = PointListMeta { + point_list: roty, + dim: self.dim, + count: self.count, + }; + + Some(roty.do_expand()) + } else { + None + }; + + let seed = self.do_expand(); + + let w = z + .into_iter() + .flatten() + .chain(y.into_iter().flatten()) + .chain(x.into_iter().flatten()) + .chain(seed); + w + } +} + +impl From for PointListMeta { + fn from(value: RawPCube) -> Self { + let (x, y, z) = value.dims(); + let dim = Dim { + x: x as usize, + y: y as usize, + z: z as usize, + }; + let point_list = value.into(); + PointListMeta { + point_list, + dim, + count: point_list.extrapolate_count(), + } + } +} + +impl From> for RawPCube { + fn from(value: PointListMeta) -> Self { + value.point_list.into() + } +} + +impl PolyCube for PointListMeta { + fn canonical_form(&self) -> Self { + PointListMeta { + point_list: self.point_list.to_min_rot_points(self.dims(), self.size()), + count: self.count, + dim: self.dim, + } + } + + fn size(&self) -> usize { + self.count + } + + fn dims(&self) -> Dim { + self.dim + } + + type UniqueExpansionsIterator = ExpandIter; + + fn unique_expansions(&self) -> Self::UniqueExpansionsIterator { + self.expand() + } +} diff --git a/rust/src/polycubes/point_list/rotate.rs b/rust/src/polycubes/point_list/rotate.rs new file mode 100644 index 0000000..2c56b35 --- /dev/null +++ b/rust/src/polycubes/point_list/rotate.rs @@ -0,0 +1,121 @@ +use std::cmp::min; + +use crate::polycubes::rotation_reduced::rotate::MatrixCol; + +use crate::polycubes::point_list::{CubeMapPos, Dim}; + +use MatrixCol::*; + +// NOTE: this could technically be an `fn`, but iterating over the tuples +// slows down the program by 2x. +macro_rules ! rot_matrix_points { + ($self:expr, $shape:expr, $count:expr, $res:expr, $(($x:expr, $y:expr, $z:expr),)*) => { + $( + $res = min($res, $self.rot_matrix_points($shape, $count, $x, $y, $z, $res.cubes[0])); + )* + } +} + +macro_rules! def_rot_matrix_points { + ($name:ident, $(($x:expr, $y:expr, $z:expr)),*) => { + #[inline(always)] + fn $name(&self, shape: Dim, count: usize, res: &mut CubeMapPos) { + rot_matrix_points!(self, shape, count, *res, $(($x, $y, $z),)*); + } + }; +} + +impl CubeMapPos { + #[inline] + pub fn rot_matrix_points( + &self, + shape: Dim, + count: usize, + x_col: MatrixCol, + y_col: MatrixCol, + z_col: MatrixCol, + pmin: u16, + ) -> CubeMapPos { + let mut res = CubeMapPos::new(); + let mut mmin = 1024; + for (i, coord) in self.cubes[0..count].iter().enumerate() { + let ix = coord & 0x1f; + let iy = (coord >> 5) & 0x1f; + let iz = (coord >> 10) & 0x1f; + let dx = Self::map_coord(ix, iy, iz, &shape, x_col); + let dy = Self::map_coord(ix, iy, iz, &shape, y_col); + let dz = Self::map_coord(ix, iy, iz, &shape, z_col); + let v = (dz << 10) | (dy << 5) | dx; + mmin = min(mmin, v); + res.cubes[i] = v; + } + //shorcut sorting because sort used to be >55% of runtime + if pmin < mmin { + res.cubes[0] = 1 << 10; + return res; + } + res.cubes[0..count].sort_unstable(); + res + } + + def_rot_matrix_points!( + xy_rots_points, + (YN, XN, ZN), + (YP, XP, ZN), + (YP, XN, ZP), + (YN, XP, ZP) + ); + + def_rot_matrix_points!( + yz_rots_points, + (XN, ZP, YP), + (XN, ZN, YN), + (XP, ZP, YN), + (XP, ZN, YP) + ); + + def_rot_matrix_points!( + xyz_rots_points, + // xz + (ZP, YP, XN), + (ZN, YN, XN), + (ZN, YP, XP), + (ZP, YN, XP), + // xyz + (ZP, XN, YN), + (YP, ZP, XP), + (YN, ZN, XP), + (ZN, XP, YN), + (YP, ZN, XN), + (YN, ZP, XN), + (ZN, XN, YP), + (ZP, XP, YP) + ); + + pub fn to_min_rot_points(&self, shape: Dim, count: usize) -> CubeMapPos { + let mut res = *self; + if shape.x == shape.y && shape.x != 0 { + self.xy_rots_points(shape, count, &mut res); + } + + if shape.y == shape.z && shape.y != 0 { + self.yz_rots_points(shape, count, &mut res); + } + + if shape.x == shape.y && shape.y == shape.z && shape.x != 0 { + self.xyz_rots_points(shape, count, &mut res); + } + + rot_matrix_points!( + self, + shape, + count, + res, + (XP, YN, ZN), + (XN, YP, ZN), + (XN, YN, ZP), + ); + + res + } +} diff --git a/rust/src/polycubes/rotation_reduced/expand.rs b/rust/src/polycubes/rotation_reduced/expand.rs new file mode 100644 index 0000000..8b13789 --- /dev/null +++ b/rust/src/polycubes/rotation_reduced/expand.rs @@ -0,0 +1 @@ + diff --git a/rust/src/polycubes/rotation_reduced/mod.rs b/rust/src/polycubes/rotation_reduced/mod.rs new file mode 100644 index 0000000..359c86c --- /dev/null +++ b/rust/src/polycubes/rotation_reduced/mod.rs @@ -0,0 +1,178 @@ +pub mod expand; +pub mod rotate; + +use hashbrown::HashSet; + +use super::Dim; + +//CubeRow is an integer type either u16 or u32 based on flags +#[derive(PartialEq, Eq, Hash, Clone, Copy, Debug)] +pub struct CubeMap { + pub x: u32, //max x index (len(xs) - 1) + pub y: u32, //max y index (len(ys) - 1) + pub z: u32, //max z index (len(zs) - 1) + pub cube_map: [u16; 6 * 6], +} + +impl PartialOrd for CubeMap { + fn partial_cmp(&self, other: &Self) -> Option { + match self.cube_map.partial_cmp(&other.cube_map) { + Some(core::cmp::Ordering::Equal) => {} + ord => return ord, + } + match self.x.partial_cmp(&other.x) { + Some(core::cmp::Ordering::Equal) => {} + ord => return ord, + } + match self.y.partial_cmp(&other.y) { + Some(core::cmp::Ordering::Equal) => {} + ord => return ord, + } + self.z.partial_cmp(&other.z) + } +} + +impl CubeMap { + /// returns 1 if it block at xyz exists + /// returns 0 if it doesnt + pub fn get_block(&self, x: usize, y: usize, z: usize) -> u16 { + (self.cube_map[z * (self.y as usize + 1) + y] >> x) & 1 + } + + /// sets the block at xyz to exist + pub fn set_block(&mut self, x: usize, y: usize, z: usize) { + self.cube_map[z * (self.y as usize + 1) + y] |= 1 << x; + } + /// set a block to the bit v + /// IMORTANT: Sets, does not unset, performs an | on the vale will never clear even on set v = 0 + pub fn set_block_to(&mut self, x: usize, y: usize, z: usize, v: u16) { + self.cube_map[z * (self.y as usize + 1) + y] |= v << x; + } + pub fn clear(&mut self) { + for i in 0..((self.z as usize + 1) * (self.y as usize + 1)) { + self.cube_map[i] = 0; + } + } + + /// ensure expected number of cubes are set, only used as an integrity check + pub fn count_cubes(&self) -> usize { + let mut c = 0; + for i in 0..36 { + let mut x = self.cube_map[i]; + while x > 0 { + c += x as usize & 1; + x >>= 1; + } + } + c + } + /// ensure no blocks are set outside expected area, only used as an integrity check + pub fn validate_bounds(&self) -> bool { + for x in (self.x + 1)..16 as u32 { + for y in 0..=self.y { + for z in 0..=self.z { + if self.get_block(x as usize, y as usize, z as usize) == 1 { + return false; + } + } + } + } + for i in ((self.z as usize + 1) * (self.y as usize + 1))..36 { + if self.cube_map[i] != 0 { + return false; + } + } + true + } + + /// find an existing block to seed continuity check + fn find_a_block(&self) -> Dim { + for y in 0..=self.y { + for z in 0..=self.z { + let mut x = 0; + let mut row = self.cube_map[(z * (self.y + 1) + y) as usize]; + if row != 0 { + while row > 0 { + if row & 1 == 1 { + return Dim { + x: x as usize, + y: y as usize, + z: z as usize, + }; + } + x += 1; + row >>= 1; + } + } + } + } + panic!("{:?} empty", self); + } + + /// ensure all blocks are connected, only used as an integrity check + pub fn validate_continuity(&self) -> bool { + let mut to_visit = HashSet::new(); + let mut map = *self; + let start = self.find_a_block(); + to_visit.insert(start); + while let Some(p) = to_visit.iter().next() { + let p = *p; + to_visit.remove(&p); + map.cube_map[p.z * (map.y as usize + 1) + p.y] &= !(1 << p.x); + if p.x > 0 && (map.cube_map[p.z * (map.y as usize + 1) + p.y] >> (p.x - 1)) & 1 == 1 { + to_visit.insert(Dim { + x: p.x - 1, + y: p.y, + z: p.z, + }); + } + if p.x < map.x as usize + && (map.cube_map[p.z * (map.y as usize + 1) + p.y] >> (p.x + 1)) & 1 == 1 + { + to_visit.insert(Dim { + x: p.x + 1, + y: p.y, + z: p.z, + }); + } + if p.y > 0 && (map.cube_map[p.z * (map.y as usize + 1) + (p.y - 1)] >> p.x) & 1 == 1 { + to_visit.insert(Dim { + x: p.x, + y: p.y - 1, + z: p.z, + }); + } + if p.y < map.y as usize + && (map.cube_map[p.z * (map.y as usize + 1) + (p.y + 1)] >> p.x) & 1 == 1 + { + to_visit.insert(Dim { + x: p.x, + y: p.y + 1, + z: p.z, + }); + } + if p.z > 0 && (map.cube_map[(p.z - 1) * (map.y as usize + 1) + p.y] >> p.x) & 1 == 1 { + to_visit.insert(Dim { + x: p.x, + y: p.y, + z: p.z - 1, + }); + } + if p.z < map.z as usize + && (map.cube_map[(p.z + 1) * (map.y as usize + 1) + p.y] >> p.x) & 1 == 1 + { + to_visit.insert(Dim { + x: p.x, + y: p.y, + z: p.z + 1, + }); + } + } + for i in 0..((map.y + 1) * (map.z + 1)) { + if map.cube_map[i as usize] != 0 { + return false; + } + } + true + } +} diff --git a/rust/src/rotations.rs b/rust/src/polycubes/rotation_reduced/rotate.rs similarity index 52% rename from rust/src/rotations.rs rename to rust/src/polycubes/rotation_reduced/rotate.rs index c9834a9..2f431ee 100644 --- a/rust/src/rotations.rs +++ b/rust/src/polycubes/rotation_reduced/rotate.rs @@ -1,6 +1,4 @@ -use std::cmp::min; - -use crate::polycube_reps::{CubeMap, CubeMapPos, CubeRow, Dim}; +use super::CubeMap; #[derive(PartialEq, Eq, Clone, Copy, Debug)] pub enum MatrixCol { @@ -13,7 +11,7 @@ pub enum MatrixCol { } #[inline] -pub fn reverse_bits(mut n: CubeRow, count: u32) -> CubeRow { +pub fn reverse_bits(mut n: u16, count: u32) -> u16 { let mut res = 0; for _ in 0..count { res = (res << 1) + (n & 1); @@ -325,393 +323,3 @@ pub fn to_min_rot(map: &CubeMap) -> CubeMap { } res } - -#[inline] -pub fn map_coord(x: u16, y: u16, z: u16, shape: &Dim, col: MatrixCol) -> u16 { - match col { - MatrixCol::XP => x, - MatrixCol::XN => shape.x as u16 - x, - MatrixCol::YP => y, - MatrixCol::YN => shape.y as u16 - y, - MatrixCol::ZP => z, - MatrixCol::ZN => shape.z as u16 - z, - } -} - -#[inline] -pub fn rot_matrix_points( - map: &CubeMapPos, - shape: &Dim, - count: usize, - x_col: MatrixCol, - y_col: MatrixCol, - z_col: MatrixCol, - pmin: u16, -) -> CubeMapPos { - let mut res = CubeMapPos::new(); - let mut mmin = 1024; - for (i, coord) in map.cubes[0..count].iter().enumerate() { - let ix = coord & 0x1f; - let iy = (coord >> 5) & 0x1f; - let iz = (coord >> 10) & 0x1f; - let dx = map_coord(ix, iy, iz, shape, x_col); - let dy = map_coord(ix, iy, iz, shape, y_col); - let dz = map_coord(ix, iy, iz, shape, z_col); - let v = (dz << 10) | (dy << 5) | dx; - mmin = min(mmin, v); - res.cubes[i] = v; - } - //shorcut sorting because sort used to be >55% of runtime - if pmin < mmin { - res.cubes[0] = 1 << 10; - return res; - } - res.cubes[0..count].sort_unstable(); - res -} - -#[inline] -fn xy_rots_points( - map: &CubeMapPos, - shape: &Dim, - count: usize, - res: &mut CubeMapPos, -) { - *res = min( - *res, - rot_matrix_points( - map, - shape, - count, - MatrixCol::YN, - MatrixCol::XN, - MatrixCol::ZN, - res.cubes[0], - ), - ); - - *res = min( - *res, - rot_matrix_points( - map, - shape, - count, - MatrixCol::YP, - MatrixCol::XP, - MatrixCol::ZN, - res.cubes[0], - ), - ); - - *res = min( - *res, - rot_matrix_points( - map, - shape, - count, - MatrixCol::YP, - MatrixCol::XN, - MatrixCol::ZP, - res.cubes[0], - ), - ); - - *res = min( - *res, - rot_matrix_points( - map, - shape, - count, - MatrixCol::YN, - MatrixCol::XP, - MatrixCol::ZP, - res.cubes[0], - ), - ); -} - -#[inline] -fn yz_rots_points( - map: &CubeMapPos, - shape: &Dim, - count: usize, - res: &mut CubeMapPos, -) { - *res = min( - *res, - rot_matrix_points( - map, - shape, - count, - MatrixCol::XN, - MatrixCol::ZP, - MatrixCol::YP, - res.cubes[0], - ), - ); - - *res = min( - *res, - rot_matrix_points( - map, - shape, - count, - MatrixCol::XN, - MatrixCol::ZN, - MatrixCol::YN, - res.cubes[0], - ), - ); - - *res = min( - *res, - rot_matrix_points( - map, - shape, - count, - MatrixCol::XP, - MatrixCol::ZP, - MatrixCol::YN, - res.cubes[0], - ), - ); - - *res = min( - *res, - rot_matrix_points( - map, - shape, - count, - MatrixCol::XP, - MatrixCol::ZN, - MatrixCol::YP, - res.cubes[0], - ), - ); -} - -#[inline] -fn xyz_rots_points( - map: &CubeMapPos, - shape: &Dim, - count: usize, - res: &mut CubeMapPos, -) { - //xz - *res = min( - *res, - rot_matrix_points( - map, - shape, - count, - MatrixCol::ZP, - MatrixCol::YP, - MatrixCol::XN, - res.cubes[0], - ), - ); - - *res = min( - *res, - rot_matrix_points( - map, - shape, - count, - MatrixCol::ZN, - MatrixCol::YN, - MatrixCol::XN, - res.cubes[0], - ), - ); - - *res = min( - *res, - rot_matrix_points( - map, - shape, - count, - MatrixCol::ZN, - MatrixCol::YP, - MatrixCol::XP, - res.cubes[0], - ), - ); - - *res = min( - *res, - rot_matrix_points( - map, - shape, - count, - MatrixCol::ZP, - MatrixCol::YN, - MatrixCol::XP, - res.cubes[0], - ), - ); - - //xyz - *res = min( - *res, - rot_matrix_points( - map, - shape, - count, - MatrixCol::ZP, - MatrixCol::XN, - MatrixCol::YN, - res.cubes[0], - ), - ); - - *res = min( - *res, - rot_matrix_points( - map, - shape, - count, - MatrixCol::YP, - MatrixCol::ZP, - MatrixCol::XP, - res.cubes[0], - ), - ); - - *res = min( - *res, - rot_matrix_points( - map, - shape, - count, - MatrixCol::YN, - MatrixCol::ZN, - MatrixCol::XP, - res.cubes[0], - ), - ); - - *res = min( - *res, - rot_matrix_points( - map, - shape, - count, - MatrixCol::ZN, - MatrixCol::XP, - MatrixCol::YN, - res.cubes[0], - ), - ); - - *res = min( - *res, - rot_matrix_points( - map, - shape, - count, - MatrixCol::YP, - MatrixCol::ZN, - MatrixCol::XN, - res.cubes[0], - ), - ); - - *res = min( - *res, - rot_matrix_points( - map, - shape, - count, - MatrixCol::YN, - MatrixCol::ZP, - MatrixCol::XN, - res.cubes[0], - ), - ); - - *res = min( - *res, - rot_matrix_points( - map, - shape, - count, - MatrixCol::ZN, - MatrixCol::XN, - MatrixCol::YP, - res.cubes[0], - ), - ); - - *res = min( - *res, - rot_matrix_points( - map, - shape, - count, - MatrixCol::ZP, - MatrixCol::XP, - MatrixCol::YP, - res.cubes[0], - ), - ); -} - -pub fn to_min_rot_points( - map: &CubeMapPos, - shape: &Dim, - count: usize, -) -> CubeMapPos { - let mut res = *map; - if shape.x == shape.y && shape.x != 0 { - xy_rots_points(map, shape, count, &mut res); - } - - if shape.y == shape.z && shape.y != 0 { - yz_rots_points(map, shape, count, &mut res); - } - - if shape.x == shape.y && shape.y == shape.z && shape.x != 0 { - xyz_rots_points(map, shape, count, &mut res); - } - - res = min( - res, - rot_matrix_points( - map, - shape, - count, - MatrixCol::XP, - MatrixCol::YN, - MatrixCol::ZN, - res.cubes[0], - ), - ); - - res = min( - res, - rot_matrix_points( - map, - shape, - count, - MatrixCol::XN, - MatrixCol::YP, - MatrixCol::ZN, - res.cubes[0], - ), - ); - - res = min( - res, - rot_matrix_points( - map, - shape, - count, - MatrixCol::XN, - MatrixCol::YN, - MatrixCol::ZP, - res.cubes[0], - ), - ); - - res -} diff --git a/rust/src/rotation_reduced.rs b/rust/src/rotation_reduced.rs index 22854cb..c43713e 100644 --- a/rust/src/rotation_reduced.rs +++ b/rust/src/rotation_reduced.rs @@ -1,26 +1,17 @@ -#[cfg(feature = "diagnostics")] -#[cfg(not(feature = "smallset"))] -use std::collections::HashMap; use std::{cmp::max, collections::HashSet, time::Instant}; use indicatif::ProgressBar; -use crate::{ - polycube_reps::CubeMap, - rotations::{self, rot_matrix, MatrixCol}, +use crate::polycubes::{ + point_list::CubeMapPos, + rotation_reduced::{ + rotate::{self, rot_matrix, MatrixCol}, + CubeMap, + }, }; -#[cfg(feature = "diagnostics")] -#[cfg(feature = "size16")] -pub static MAX_X: usize = 16; - -#[cfg(feature = "diagnostics")] -#[cfg(not(feature = "size16"))] -pub static MAX_X: usize = 32; - -#[cfg(feature = "smallset")] ///converts a cube map to map pos for hashset storage slow (+10% runtime combined with decode last measured) -fn cube_map_to_cube_map_pos(map: &CubeMap) -> CubeMapPos { +fn cube_map_to_cube_map_pos(map: &CubeMap) -> CubeMapPos<16> { let mut pos = CubeMapPos { cubes: [0; 16] }; let mut i = 0; for z in 0..=map.z as usize { @@ -34,19 +25,11 @@ fn cube_map_to_cube_map_pos(map: &CubeMap) -> CubeMapPos { } } pos.cubes[i - 1] |= 0x8000; - #[cfg(feature = "diagnostics")] - { - let a = cube_map_from_cube_map_pos(&pos); - if a != *map { - panic!("{:?} {:?} unequal", a, map); - } - } pos } -#[cfg(feature = "smallset")] ///converts a mappos from hashset storage to a cube map -fn cube_map_from_cube_map_pos(map: &CubeMapPos) -> CubeMap { +fn cube_map_from_cube_map_pos(map: &CubeMapPos<16>) -> CubeMap { let mut dst = CubeMap { x: 0, y: 0, @@ -81,31 +64,9 @@ fn cube_map_from_cube_map_pos(map: &CubeMapPos) -> CubeMap { dst } -#[cfg(feature = "smallset")] -type CubeEncoding = CubeMapPos; -#[cfg(not(feature = "smallset"))] -type CubeEncoding = CubeMap; - #[inline] -fn insert_map( - map: &CubeMap, - seen: &mut HashSet, - #[cfg(feature = "diagnostics")] depth: usize, -) { - let work_map = rotations::to_min_rot(map); - #[cfg(feature = "diagnostics")] - { - if map.count_cubes() != depth { - panic!("{:?} doesnt have {} cubes", map, depth) - } - if !map.validate_bounds() { - panic!("{:?} has blocks out of bounds", map) - } - if !map.validate_continuity() { - panic!("{:?} non continuous", map) - } - } - #[cfg(feature = "smallset")] +fn insert_map(map: &CubeMap, seen: &mut HashSet>) { + let work_map = rotate::to_min_rot(map); let work_map = cube_map_to_cube_map_pos(&work_map); seen.insert(work_map); } @@ -207,33 +168,19 @@ fn expand_cube_map_out(map: &CubeMap, y: usize, z: usize, offset: u32) -> CubeMa /// expand each cube +/-1 X where possible #[inline] -fn expand_xs( - map: &CubeMap, - seen: &mut HashSet, - #[cfg(feature = "diagnostics")] depth: usize, -) { +fn expand_xs(map: &CubeMap, seen: &mut HashSet>) { for yz in 0..(((map.y + 1) * (map.z + 1)) as usize) { let left_bits = ((map.cube_map[yz] << 1) | map.cube_map[yz]) ^ map.cube_map[yz]; let right_bits = ((map.cube_map[yz] << 1) | map.cube_map[yz]) ^ (map.cube_map[yz] << 1); for xoff in 1..(map.x + 2) { //start at 1 because shifting left cant be in the zero bit if left_bits & (1 << xoff) != 0 { - insert_map( - &expand_cube_map_left(map, yz, xoff), - seen, - #[cfg(feature = "diagnostics")] - depth, - ); + insert_map(&expand_cube_map_left(map, yz, xoff), seen); } } for xoff in 0..(map.x + 1) { if right_bits & (1 << xoff) != 0 { - insert_map( - &expand_cube_map_right(map, yz, xoff), - seen, - #[cfg(feature = "diagnostics")] - depth, - ); + insert_map(&expand_cube_map_right(map, yz, xoff), seen); } } } @@ -241,11 +188,7 @@ fn expand_xs( /// expand each cube +/-1 Y where possible #[inline] -fn expand_ys( - map: &CubeMap, - seen: &mut HashSet, - #[cfg(feature = "diagnostics")] depth: usize, -) { +fn expand_ys(map: &CubeMap, seen: &mut HashSet>) { for z in 0..=map.z as usize { for y in 0..=map.y as usize { let up_bits = if y == map.y as usize { @@ -263,20 +206,10 @@ fn expand_ys( for xoff in 0..=map.x { //start at 1 because shifting left cant be in the zero bit if up_bits & (1 << xoff) != 0 { - insert_map( - &expand_cube_map_up(map, y + 1, z, xoff), - seen, - #[cfg(feature = "diagnostics")] - depth, - ); + insert_map(&expand_cube_map_up(map, y + 1, z, xoff), seen); } if down_bits & (1 << xoff) != 0 { - insert_map( - &expand_cube_map_down(map, y, z, xoff), - seen, - #[cfg(feature = "diagnostics")] - depth, - ); + insert_map(&expand_cube_map_down(map, y, z, xoff), seen); } } } @@ -285,11 +218,7 @@ fn expand_ys( /// expand each cube +/-1 Z where possible #[inline] -fn expand_zs( - map: &CubeMap, - seen: &mut HashSet, - #[cfg(feature = "diagnostics")] depth: usize, -) { +fn expand_zs(map: &CubeMap, seen: &mut HashSet>) { for z in 0..=map.z as usize { for y in 0..=map.y as usize { let in_bits = map.cube_map[z * (map.y as usize + 1) + y] @@ -302,20 +231,10 @@ fn expand_zs( }; for xoff in 0..(map.x + 1) { if in_bits & (1 << xoff) != 0 { - insert_map( - &expand_cube_map_in(map, y, z + 1, xoff), - seen, - #[cfg(feature = "diagnostics")] - depth, - ); + insert_map(&expand_cube_map_in(map, y, z + 1, xoff), seen); } if out_bits & (1 << xoff) != 0 { - insert_map( - &expand_cube_map_out(map, y, z, xoff), - seen, - #[cfg(feature = "diagnostics")] - depth, - ); + insert_map(&expand_cube_map_out(map, y, z, xoff), seen); } } } @@ -324,52 +243,20 @@ fn expand_zs( /// expand in X, Y and Z abiding by the X >= Y >= Z constraint #[inline] -fn do_cube_expansion( - map: &CubeMap, - seen: &mut HashSet, - #[cfg(feature = "diagnostics")] depth: usize, -) { - expand_xs( - map, - seen, - #[cfg(feature = "diagnostics")] - depth, - ); +fn do_cube_expansion(map: &CubeMap, seen: &mut HashSet>) { + expand_xs(map, seen); if map.y < map.x { - expand_ys( - map, - seen, - #[cfg(feature = "diagnostics")] - depth, - ); + expand_ys(map, seen); } if map.z < map.y { - expand_zs( - map, - seen, - #[cfg(feature = "diagnostics")] - depth, - ); + expand_zs(map, seen); } } /// expand cube, rotate around square faces to catch adgecases that were getting missed due to the X >= Y >= Z constraint #[inline] -fn expand_cube_map( - map: &CubeMap, - seen: &mut HashSet, - #[cfg(feature = "diagnostics")] depth: usize, -) { - do_cube_expansion( - map, - seen, - #[cfg(feature = "diagnostics")] - depth, - ); - #[cfg(feature = "diagnostics")] - if map.count_cubes() != depth - 1 { - panic!("{:?} doesnt have {} cubes", map, depth - 1) - } +fn expand_cube_map(map: &CubeMap, seen: &mut HashSet>) { + do_cube_expansion(map, seen); if map.x == map.y && map.x > 0 { let mut rot = CubeMap { x: map.x, @@ -378,16 +265,7 @@ fn expand_cube_map( cube_map: [0; 36], }; rot_matrix(map, &mut rot, MatrixCol::YN, MatrixCol::XN, MatrixCol::ZN); - #[cfg(feature = "diagnostics")] - if rot.count_cubes() != depth - 1 { - panic!("{:?} doesnt have {} cubes", rot, depth - 1) - } - do_cube_expansion( - &rot, - seen, - #[cfg(feature = "diagnostics")] - depth, - ); + do_cube_expansion(&rot, seen); } if map.y == map.z && map.y > 0 { let mut rot = CubeMap { @@ -397,16 +275,7 @@ fn expand_cube_map( cube_map: [0; 36], }; rot_matrix(map, &mut rot, MatrixCol::XN, MatrixCol::ZP, MatrixCol::YP); - #[cfg(feature = "diagnostics")] - if rot.count_cubes() != depth - 1 { - panic!("{:?} doesnt have {} cubes", rot, depth - 1) - } - do_cube_expansion( - &rot, - seen, - #[cfg(feature = "diagnostics")] - depth, - ); + do_cube_expansion(&rot, seen); } if map.x == map.z && map.x > 0 { let mut rot = CubeMap { @@ -416,46 +285,20 @@ fn expand_cube_map( cube_map: [0; 36], }; rot_matrix(map, &mut rot, MatrixCol::ZP, MatrixCol::YP, MatrixCol::XN); - #[cfg(feature = "diagnostics")] - if rot.count_cubes() != depth - 1 { - panic!("{:?} doesnt have {} cubes", rot, depth - 1) - } - do_cube_expansion( - &rot, - seen, - #[cfg(feature = "diagnostics")] - depth, - ); - } -} - -#[cfg(feature = "diagnostics")] -#[cfg(not(feature = "smallset"))] -fn to_dim(cm: &CubeMap) -> Dim { - Dim { - x: cm.x as usize + 1, - y: cm.y as usize + 1, - z: cm.z as usize + 1, + do_cube_expansion(&rot, seen); } } /// expand all polycubes in set n-1 to get set n fn expand_cube_set( - in_set: &HashSet, - #[cfg(feature = "diagnostics")] depth: usize, - out_set: &mut HashSet, + in_set: &HashSet>, + out_set: &mut HashSet>, bar: &ProgressBar, ) { let mut i = 0; for map in in_set.iter() { - #[cfg(feature = "smallset")] let map = &cube_map_from_cube_map_pos(map); - expand_cube_map( - map, - out_set, - #[cfg(feature = "diagnostics")] - depth, - ); + expand_cube_map(map, out_set); i += 1; if i == 100 { bar.inc(100); @@ -463,27 +306,6 @@ fn expand_cube_set( } } bar.inc(i); - #[cfg(feature = "diagnostics")] - #[cfg(not(feature = "smallset"))] - { - let mut shape_map = HashMap::new(); - for map in out_set.iter() { - if map.count_cubes() != depth { - panic!("{:?} doesnt have {} cubes", map, depth) - } - let dim = to_dim(map); - shape_map.insert( - dim, - match shape_map.get(&dim) { - Some(res) => res + 1, - None => 1, - }, - ); - } - for (s, map) in shape_map.iter() { - println!("{}, {}, {} -> {:?}", s.x, s.y, s.z, map); - } - } } pub fn gen_polycubes(n: usize, bar: &ProgressBar) -> usize { @@ -499,23 +321,11 @@ pub fn gen_polycubes(n: usize, bar: &ProgressBar) -> usize { let t1_start = Instant::now(); let mut seeds = HashSet::new(); let mut dst = HashSet::new(); - insert_map( - &unit_cube, - &mut seeds, - #[cfg(feature = "diagnostics")] - 2, - ); + insert_map(&unit_cube, &mut seeds); for i in 3..=n as usize { bar.set_message(format!("seed subsets expanded for N = {}...", i)); - expand_cube_set( - &seeds, - #[cfg(feature = "diagnostics")] - i, - &mut dst, - bar, - ); - //if diagnostics enabled panic if the returned values are wrong - #[cfg(feature = "diagnostics")] + expand_cube_set(&seeds, &mut dst, bar); + // panic if the returned values are wrong if i == 3 && dst.len() != 2 { panic!("{} supposed to have {} elems not {}", i, 2, dst.len()) } else if i == 4 && dst.len() != 8 { diff --git a/rust/src/test.rs b/rust/src/test.rs index ad573e2..0a03da0 100644 --- a/rust/src/test.rs +++ b/rust/src/test.rs @@ -1,6 +1,6 @@ use std::collections::HashSet; -use crate::naive_polycube::NaivePolyCube; +use crate::polycubes::naive_polycube::NaivePolyCube; fn test_cube() -> NaivePolyCube { let mut cube = NaivePolyCube::new(2, 3, 4);