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..2e66ff3 100644 --- a/rust/src/cli/enumerate.rs +++ b/rust/src/cli/enumerate.rs @@ -3,11 +3,8 @@ 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, - rotation_reduced, + rotation_reduced, polycubes::{pcube::{RawPCube, PCubeFile}, naive_polycube::NaivePolyCube, point_list::CubeMapPos}, }; use crate::{finish_bar, make_bar, unknown_bar, Compression, EnumerateOpts, EnumerationMode}; @@ -83,7 +80,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..50518a1 100644 --- a/rust/src/hashless.rs +++ b/rust/src/hashless.rs @@ -1,110 +1,11 @@ -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, Dim}; 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 +15,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 +34,24 @@ 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); + seed.expand(shape, count) + .for_each(|(dim, count, map)| store.insert_map(dim, map, count)); - map.inner + 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..064a29a 100644 --- a/rust/src/lib.rs +++ b/rust/src/lib.rs @@ -2,12 +2,9 @@ mod test; pub mod iterator; -pub mod pcube; +pub mod polycubes; -pub mod naive_polycube; 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..bc08732 100644 --- a/rust/src/pointlist.rs +++ b/rust/src/pointlist.rs @@ -1,309 +1,130 @@ -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, 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>>>, } -///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; - return; +impl MapStore { + pub fn new() -> Self { + Self { + inner: HashMap::new(), } } - 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]; - } -} + fn insert(&self, dim: Dim, map: CubeMapPos<16>, count: usize) { + let map = map.to_min_rot_points(dim, count); + let mut body = CubeMapPos::new(); -/// 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; + body.cubes[0..count].copy_from_slice(&map.cubes[1..count + 1]); - 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 entry = self + .inner + .get(&(dim, map.cubes[0])) + .expect("Cube size does not have entry in destination map"); -/// 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) - } + entry.write().insert(body); } -} -/// 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) - } - } -} + /// 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], + }; -/// 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); -} + for seed_body in body { + for i in 1..count { + seed.cubes[i] = seed_body.cubes[i - 1]; + } -/// 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); -} + // body.cubes.copy_within(0..body.cubes.len() - 1, 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]; + seed.expand(shape, count) + .for_each(|(dim, count, map)| self.insert(dim, map, count)); } - 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_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 +132,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 +151,9 @@ fn _clone_polycubes_to_vec(maps: &mut MapStore) -> Vec> { v.push(cmp); } } + + v } - v } /// run pointlist based generation algorithm @@ -355,66 +172,34 @@ 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); + seeds.insert(dim, seed, calculate_from - 1); } 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/polycubes/mod.rs b/rust/src/polycubes/mod.rs new file mode 100644 index 0000000..2a39175 --- /dev/null +++ b/rust/src/polycubes/mod.rs @@ -0,0 +1,5 @@ + +pub mod point_list; +pub mod naive_polycube; +pub mod pcube; +pub mod rotation_reduced; \ No newline at end of file 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..e69de29 diff --git a/rust/src/polycube_reps.rs b/rust/src/polycubes/point_list/mod.rs similarity index 62% rename from rust/src/polycube_reps.rs rename to rust/src/polycubes/point_list/mod.rs index ac82979..f91ab0c 100644 --- a/rust/src/polycube_reps.rs +++ b/rust/src/polycubes/point_list/mod.rs @@ -1,8 +1,22 @@ +pub mod expand; +pub mod rotate; + use std::cmp::{max, min}; -use crate::pcube::RawPCube; +use super::{pcube::RawPCube, rotation_reduced::rotate::MatrixCol}; +use crate::polycubes::rotation_reduced::rotate::MatrixCol::*; -use crate::rotations::{map_coord, to_min_rot_points, MatrixCol::*}; +/// 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, +} /// Polycube representation /// stores up to 16 blocks (number of cubes normally implicit or seperate in program state) @@ -95,9 +109,9 @@ impl From<&RawPCube> for CubeMapPos { 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); + 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) } @@ -141,10 +155,47 @@ impl From> for RawPCube { } } +impl CubeMapPos<0> { + ///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] { @@ -152,6 +203,7 @@ impl CubeMapPos { } count } + pub fn extrapolate_dim(&self) -> Dim { let count = self.extrapolate_count(); let mut dim = Dim { x: 0, y: 0, z: 0 }; @@ -251,9 +303,9 @@ impl CubeMapPos { 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 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;} @@ -311,7 +363,7 @@ impl CubeMapPos { dim = rdim; root_candidate.cubes[0..count].sort_unstable(); } - let mrp = to_min_rot_points(&root_candidate, &dim, count); + let mrp = root_candidate.to_min_rot_points(dim, count); if &mrp < seed { return false; } @@ -320,201 +372,154 @@ impl CubeMapPos { } } -/// 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], -} +macro_rules! cube_map_pos_expand { + ($name:ident, $dim:ident, $shift:literal) => { + #[inline(always)] + pub fn $name(self, shape: Dim, count: usize) -> impl Iterator { + struct Iter { + inner: CubeMapPos, + shape: Dim, + count: usize, + i: usize, + stored: Option<(Dim, usize, CubeMapPos)>, + } -/// 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, -} + impl<'a, const C: usize> Iterator for Iter { + type Item = (Dim, usize, CubeMapPos); -#[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], -} + fn next(&mut self) -> Option { + loop { + if let Some(stored) = self.stored.take() { + return Some(stored); + } -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) - } -} + let i = self.i; -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 - } + 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; + + CubeMapPos::array_insert(plus, &mut new_map.cubes[i..=self.count]); + new_shape.$dim = + max(new_shape.$dim, (((coord >> $shift) + 1) & 0x1f) as usize); - #[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, - }; + self.stored = Some((new_shape, self.count + 1, new_map)); } - x += 1; - row >>= 1; + + 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 + }; + + CubeMapPos::array_shift(&mut new_map.cubes[i..=self.count]); + CubeMapPos::array_insert(insert_coord, &mut new_map.cubes[0..=i]); + return Some((new_shape, self.count + 1, new_map)); } } } - } - 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, - }); + Iter { + inner: self, + shape, + count, + i: 0, + stored: None, } } - for i in 0..((map.y + 1) * (map.z + 1)) { - if map.cube_map[i as usize] != 0 { - return false; - } - } - true + }; +} + +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); + + /// reduce number of expansions needing to be performed based on + /// X >= Y >= Z constraint on Dim + #[inline] + #[must_use] + fn do_expand(self, shape: Dim, count: usize) -> impl Iterator { + let expand_ys = if shape.y < shape.x { + Some(self.expand_y(shape, count)) + } else { + None + }; + + let expand_zs = if shape.z < shape.y { + Some(self.expand_z(shape, count)) + } else { + None + }; + + self.expand_x(shape, count) + .chain(expand_ys.into_iter().flatten()) + .chain(expand_zs.into_iter().flatten()) + } + + /// 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, + shape: Dim, + count: usize, + ) -> impl Iterator + '_ { + use MatrixCol::*; + + let z = if shape.x == shape.y && shape.x > 0 { + let rotz = self.rot_matrix_points(shape, count, YN, XN, ZN, 1025); + Some(rotz.do_expand(shape, count)) + } else { + None + }; + + let y = if shape.y == shape.z && shape.y > 0 { + let rotx = self.rot_matrix_points(shape, count, XN, ZP, YP, 1025); + Some(rotx.do_expand(shape, count)) + } else { + None + }; + + let x = if shape.x == shape.z && shape.x > 0 { + let roty = self.rot_matrix_points(shape, count, ZP, YP, XN, 1025); + Some(roty.do_expand(shape, count)) + } else { + None + }; + + let seed = self.do_expand(shape, count); + + z.into_iter() + .flatten() + .chain(y.into_iter().flatten()) + .chain(x.into_iter().flatten()) + .chain(seed) } } 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..e69de29 diff --git a/rust/src/polycubes/rotation_reduced/mod.rs b/rust/src/polycubes/rotation_reduced/mod.rs new file mode 100644 index 0000000..613eb2a --- /dev/null +++ b/rust/src/polycubes/rotation_reduced/mod.rs @@ -0,0 +1,179 @@ +pub mod rotate; +pub mod expand; + +use hashbrown::HashSet; + +use super::point_list::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..f8eae72 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); @@ -326,392 +324,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..b8c2b93 100644 --- a/rust/src/rotation_reduced.rs +++ b/rust/src/rotation_reduced.rs @@ -1,26 +1,12 @@ -#[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}, -}; - -#[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; +use crate::polycubes::{rotation_reduced::{CubeMap, rotate::{self, MatrixCol, rot_matrix}}, point_list::CubeMapPos}; -#[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 +20,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 +59,12 @@ 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, + seen: &mut HashSet> ) { - 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")] + let work_map = rotate::to_min_rot(map); let work_map = cube_map_to_cube_map_pos(&work_map); seen.insert(work_map); } @@ -209,8 +168,7 @@ fn expand_cube_map_out(map: &CubeMap, y: usize, z: usize, offset: u32) -> CubeMa #[inline] fn expand_xs( map: &CubeMap, - seen: &mut HashSet, - #[cfg(feature = "diagnostics")] depth: usize, + 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]; @@ -221,8 +179,6 @@ fn expand_xs( insert_map( &expand_cube_map_left(map, yz, xoff), seen, - #[cfg(feature = "diagnostics")] - depth, ); } } @@ -231,8 +187,6 @@ fn expand_xs( insert_map( &expand_cube_map_right(map, yz, xoff), seen, - #[cfg(feature = "diagnostics")] - depth, ); } } @@ -243,8 +197,7 @@ fn expand_xs( #[inline] fn expand_ys( map: &CubeMap, - seen: &mut HashSet, - #[cfg(feature = "diagnostics")] depth: usize, + seen: &mut HashSet>, ) { for z in 0..=map.z as usize { for y in 0..=map.y as usize { @@ -266,16 +219,12 @@ fn expand_ys( insert_map( &expand_cube_map_up(map, y + 1, z, xoff), seen, - #[cfg(feature = "diagnostics")] - depth, ); } if down_bits & (1 << xoff) != 0 { insert_map( &expand_cube_map_down(map, y, z, xoff), seen, - #[cfg(feature = "diagnostics")] - depth, ); } } @@ -287,8 +236,7 @@ fn expand_ys( #[inline] fn expand_zs( map: &CubeMap, - seen: &mut HashSet, - #[cfg(feature = "diagnostics")] depth: usize, + seen: &mut HashSet>, ) { for z in 0..=map.z as usize { for y in 0..=map.y as usize { @@ -305,16 +253,12 @@ fn expand_zs( insert_map( &expand_cube_map_in(map, y, z + 1, xoff), seen, - #[cfg(feature = "diagnostics")] - depth, ); } if out_bits & (1 << xoff) != 0 { insert_map( &expand_cube_map_out(map, y, z, xoff), seen, - #[cfg(feature = "diagnostics")] - depth, ); } } @@ -326,29 +270,22 @@ fn expand_zs( #[inline] fn do_cube_expansion( map: &CubeMap, - seen: &mut HashSet, - #[cfg(feature = "diagnostics")] depth: usize, + seen: &mut HashSet>, ) { expand_xs( map, seen, - #[cfg(feature = "diagnostics")] - depth, ); if map.y < map.x { expand_ys( map, seen, - #[cfg(feature = "diagnostics")] - depth, ); } if map.z < map.y { expand_zs( map, seen, - #[cfg(feature = "diagnostics")] - depth, ); } } @@ -357,19 +294,12 @@ fn do_cube_expansion( #[inline] fn expand_cube_map( map: &CubeMap, - seen: &mut HashSet, - #[cfg(feature = "diagnostics")] depth: usize, + seen: &mut HashSet>, ) { 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) - } if map.x == map.y && map.x > 0 { let mut rot = CubeMap { x: map.x, @@ -378,15 +308,9 @@ 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, ); } if map.y == map.z && map.y > 0 { @@ -397,15 +321,9 @@ 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, ); } if map.x == map.z && map.x > 0 { @@ -416,45 +334,25 @@ 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, - } -} - /// 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, ); i += 1; if i == 100 { @@ -463,27 +361,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 { @@ -502,20 +379,15 @@ pub fn gen_polycubes(n: usize, bar: &ProgressBar) -> usize { insert_map( &unit_cube, &mut seeds, - #[cfg(feature = "diagnostics")] - 2, ); 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")] + // 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);