From c14aa3b61c133ce1bb2abfb0e7463f36ab3bd82a Mon Sep 17 00:00:00 2001 From: janw20 Date: Thu, 18 Apr 2024 11:09:39 +0200 Subject: [PATCH 01/14] Add new type `PackedArray` --- pineappl/src/lib.rs | 1 + pineappl/src/packed_array.rs | 825 +++++++++++++++++++++++++++++++++++ 2 files changed, 826 insertions(+) create mode 100644 pineappl/src/packed_array.rs diff --git a/pineappl/src/lib.rs b/pineappl/src/lib.rs index 42879684f..6c1479b13 100644 --- a/pineappl/src/lib.rs +++ b/pineappl/src/lib.rs @@ -8,6 +8,7 @@ pub mod evolution; pub mod fk_table; pub mod grid; pub mod import_only_subgrid; +pub mod packed_array; pub mod lagrange_subgrid; pub mod lumi; pub mod ntuple_subgrid; diff --git a/pineappl/src/packed_array.rs b/pineappl/src/packed_array.rs new file mode 100644 index 000000000..061666563 --- /dev/null +++ b/pineappl/src/packed_array.rs @@ -0,0 +1,825 @@ +use std::ops::{Index, IndexMut, MulAssign}; +use std::{mem, vec}; + +use ndarray::ArrayView3; +use serde::{Deserialize, Serialize}; + +#[derive(Clone, Deserialize, Serialize)] +pub struct PackedArray { + entries: Vec, + start_indices: Vec, + lengths: Vec, + shape: Vec, +} + +impl PackedArray { + #[must_use] + pub fn new(shape: [usize; D]) -> Self { + Self { + // entries: BTreeMap::new(), + entries: vec![], + start_indices: vec![], + lengths: vec![], + shape: shape.to_vec(), + } + } + + #[must_use] + pub fn is_empty(&self) -> bool { + self.entries.is_empty() + } + + pub fn shape(&self) -> &[usize] { + &self.shape + } + + pub fn clear(&mut self) { + self.entries.clear(); + self.start_indices.clear(); + self.lengths.clear(); + } + + #[must_use] + pub fn overhead(&self) -> usize { + ((self.start_indices.len() + self.lengths.len()) * mem::size_of::()) + / mem::size_of::() + } + + #[must_use] + pub fn explicit_zeros(&self) -> usize { + self.entries.iter().filter(|x| **x == T::default()).count() + } + + pub fn non_zeros(&self) -> usize { + self.entries.iter().filter(|x| **x != T::default()).count() + } + + pub fn indexed_iter(&self) -> impl Iterator + '_ { + self.start_indices + .iter() + .zip(&self.lengths) + .flat_map(|(&start_index, &length)| { + (start_index..(start_index + length)).map(|i| unravel_index(i, &self.shape)) + }) + .zip(&self.entries) + .filter(|&(_, entry)| *entry != Default::default()) + .map(|(indices, entry)| (indices, *entry)) + } +} + +impl, const D: usize> MulAssign for PackedArray { + fn mul_assign(&mut self, rhs: T) { + self.entries.iter_mut().for_each(|x| *x *= rhs); + } +} + +impl PackedArray { + pub fn from_ndarray(array: ArrayView3, xstart: usize, xsize: usize) -> Self { + let shape = array.shape(); + + let mut result = Self::new([xsize, shape[1], shape[2]]); + + for ((i, j, k), &entry) in array + .indexed_iter() + .filter(|(_, &entry)| entry != Default::default()) + { + result[[i + xstart, j, k]] = entry; + } + + result + } +} + +fn ravel_multi_index(indices: &[usize], dimensions: &[usize]) -> usize { + assert_eq!(indices.len(), dimensions.len()); + + indices + .iter() + .skip(1) + .zip(dimensions.iter().skip(1)) + .fold(indices[0], |acc, (i, d)| acc * d + i) +} + +fn unravel_index(index: usize, dimensions: &[usize]) -> [usize; D] { + assert!(index < dimensions.iter().product()); + let mut indices = [0; D]; + indices + .iter_mut() + .rev() + .zip(dimensions.iter().rev()) + .fold(index, |acc, (i, d)| { + *i = acc % d; + acc / d + }); + indices +} + +impl Index<[usize; D]> for PackedArray { + type Output = T; + + fn index(&self, index: [usize; D]) -> &Self::Output { + debug_assert_eq!(index.len(), self.shape.len()); + assert!( + index.iter().zip(self.shape.iter()).all(|(&i, &d)| i < d), + "index {:?} is out of bounds for array of shape {:?}", + index, + self.shape + ); + + let raveled_index = ravel_multi_index::(&index, &self.shape); + let point = self.start_indices.partition_point(|&i| i <= raveled_index); + + assert!( + point > 0, + "entry at index {index:?} is implicitly set to the default value" + ); + + let start_index = self.start_indices[point - 1]; + let length = self.lengths[point - 1]; + + let point_entries = + self.lengths.iter().take(point - 1).sum::() + raveled_index - start_index; + + assert!( + raveled_index < (start_index + length), + "entry at index {index:?} is implicitly set to the default value" + ); + + &self.entries[point_entries] + } +} + +// TODO: implement axis swapping optimization +impl IndexMut<[usize; D]> + for PackedArray +{ + fn index_mut(&mut self, index: [usize; D]) -> &mut Self::Output { + debug_assert_eq!(index.len(), self.shape.len()); + assert!( + index.iter().zip(self.shape.iter()).all(|(&i, &d)| i < d), + "index {:?} is out of bounds for array of shape {:?}", + index, + self.shape + ); + + let raveled_index = ravel_multi_index::(&index, &self.shape); + + // if self.entries.is_empty() { + // self.start_indices.push(raveled_index); + // self.lengths.push(1); + // self.entries.push(Default::default()); + + // return &mut self.entries[0]; + // } + + // the closest start_index that is greater than raveled_index + let point = self.start_indices.partition_point(|&i| i <= raveled_index); + + // the index that is stored in `point`, translated to the indices of the entries + let point_entries = self.lengths.iter().take(point).sum::(); + + // println!("test 1"); + + // maximum distance for merging regions + let threshold_distance = 2; + + if point > 0 { + let start_index = self.start_indices[point - 1]; + let length = self.lengths[point - 1]; + + if raveled_index < start_index + length { + // println!("test 2"); + return &mut self.entries[point_entries - length + raveled_index - start_index]; + } else if raveled_index < start_index + length + threshold_distance { + // println!("test 3"); + + let distance = raveled_index - (start_index + length) + 1; + // println!("distance: {}", distance); + self.lengths[point - 1] += distance; + for _ in 0..(distance) { + self.entries.insert(point_entries, Default::default()); + } + + if let Some(start_index_next) = self.start_indices.get(point) { + if raveled_index + threshold_distance >= *start_index_next { + // println!("test 4"); + let distance_next = start_index_next - raveled_index; + + self.lengths[point - 1] += distance_next - 1 + self.lengths[point]; + self.lengths.remove(point); + self.start_indices.remove(point); + for _ in 0..(distance_next - 1) { + self.entries.insert(point_entries, Default::default()); + } + } + } + + return &mut self.entries[point_entries - 1 + distance]; + } else if let Some(start_index_next) = self.start_indices.get(point) { + if raveled_index + threshold_distance >= *start_index_next { + // println!("test 5"); + let distance = start_index_next - raveled_index; + + self.start_indices[point] = raveled_index; + self.lengths[point] += distance; + for _ in 0..distance { + self.entries.insert(point_entries, Default::default()); + } + return &mut self.entries[point_entries]; + } + } + } + + // println!("test 6"); + + // in the most general case we have to insert a new start_index + self.start_indices.insert(point, raveled_index); + self.lengths.insert(point, 1); + self.entries.insert(point_entries, Default::default()); + + &mut self.entries[point_entries] + } +} + +// impl IntoIterator for &JaggedArray { +// type Item = (Vec, T); +// type IntoIter = Iterator, T)>; + +// fn into_iter(self) -> Self::IntoIter { + +// } +// } + +#[cfg(test)] +mod tests { + use ndarray::Array3; + use std::mem; + + use super::*; + + #[test] + fn unravel_index() { + assert_eq!(super::unravel_index(0, &[3, 2]), [0, 0]); + assert_eq!(super::unravel_index(1, &[3, 2]), [0, 1]); + assert_eq!(super::unravel_index(2, &[3, 2]), [1, 0]); + assert_eq!(super::unravel_index(3, &[3, 2]), [1, 1]); + assert_eq!(super::unravel_index(4, &[3, 2]), [2, 0]); + assert_eq!(super::unravel_index(5, &[3, 2]), [2, 1]); + } + + #[test] + fn ravel_multi_index() { + assert_eq!(super::ravel_multi_index::<2>(&[0, 0], &[3, 2]), 0); + assert_eq!(super::ravel_multi_index::<2>(&[0, 1], &[3, 2]), 1); + assert_eq!(super::ravel_multi_index::<2>(&[1, 0], &[3, 2]), 2); + assert_eq!(super::ravel_multi_index::<2>(&[1, 1], &[3, 2]), 3); + assert_eq!(super::ravel_multi_index::<2>(&[2, 0], &[3, 2]), 4); + assert_eq!(super::ravel_multi_index::<2>(&[2, 1], &[3, 2]), 5); + } + + #[test] + fn index() { + let mut a = PackedArray::::new([4, 2]); + + a[[0, 0]] = 1.0; + println!( + "entries: {:?}, start_indices: {:?}, lengths: {:?}", + a.entries, a.start_indices, a.lengths + ); + assert_eq!(a[[0, 0]], 1.0); + assert_eq!(a.entries, vec![1.0]); + assert_eq!(a.start_indices, vec![0]); + assert_eq!(a.lengths, vec![1]); + + a[[3, 0]] = 2.0; + println!( + "entries: {:?}, start_indices: {:?}, lengths: {:?}", + a.entries, a.start_indices, a.lengths + ); + assert_eq!(a[[0, 0]], 1.0); + assert_eq!(a[[3, 0]], 2.0); + assert_eq!(a.entries, vec![1.0, 2.0]); + assert_eq!(a.start_indices, vec![0, 6]); + assert_eq!(a.lengths, vec![1, 1]); + + a[[3, 1]] = 3.0; + println!( + "entries: {:?}, start_indices: {:?}, lengths: {:?}", + a.entries, a.start_indices, a.lengths + ); + + assert_eq!(a[[0, 0]], 1.0); + assert_eq!(a[[3, 0]], 2.0); + assert_eq!(a[[3, 1]], 3.0); + assert_eq!(a.entries, vec![1.0, 2.0, 3.0]); + assert_eq!(a.start_indices, vec![0, 6]); + assert_eq!(a.lengths, vec![1, 2]); + + a[[2, 0]] = 3.5; + println!( + "entries: {:?}, start_indices: {:?}, lengths: {:?}", + a.entries, a.start_indices, a.lengths + ); + assert_eq!(a[[0, 0]], 1.0); + assert_eq!(a[[3, 0]], 2.0); + assert_eq!(a[[3, 1]], 3.0); + assert_eq!(a[[2, 0]], 3.5); + // assert_eq!(a[[2, 1]], 0.0); + assert_eq!(a.entries, vec![1.0, 3.5, 0.0, 2.0, 3.0]); + assert_eq!(a.start_indices, vec![0, 4]); + assert_eq!(a.lengths, vec![1, 4]); + + a[[2, 0]] = 4.0; + println!( + "entries: {:?}, start_indices: {:?}, lengths: {:?}", + a.entries, a.start_indices, a.lengths + ); + + assert_eq!(a[[0, 0]], 1.0); + assert_eq!(a[[3, 0]], 2.0); + assert_eq!(a[[3, 1]], 3.0); + assert_eq!(a[[2, 0]], 4.0); + // assert_eq!(a[[2, 1]], 0.0); + assert_eq!(a.entries, vec![1.0, 4.0, 0.0, 2.0, 3.0]); + assert_eq!(a.start_indices, vec![0, 4]); + assert_eq!(a.lengths, vec![1, 4]); + + a[[1, 0]] = 5.0; + assert_eq!(a[[0, 0]], 1.0); + assert_eq!(a[[3, 0]], 2.0); + assert_eq!(a[[3, 1]], 3.0); + assert_eq!(a[[2, 0]], 4.0); + // assert_eq!(a[[2, 1]], 0.0); + assert_eq!(a[[1, 0]], 5.0); + assert_eq!(a.entries, vec![1.0, 0.0, 5.0, 0.0, 4.0, 0.0, 2.0, 3.0]); + assert_eq!(a.start_indices, vec![0]); + assert_eq!(a.lengths, vec![8]); + + println!( + "entries: {:?}, start_indices: {:?}, lengths: {:?}", + a.entries, a.start_indices, a.lengths + ); + } + + #[test] + fn iter() { + let mut a = PackedArray::::new([6, 5]); + a[[2, 2]] = 1; + a[[2, 4]] = 2; + a[[4, 1]] = 3; + a[[4, 4]] = 4; + println!( + "entries: {:?}, start_indices: {:?}, lengths: {:?}", + a.entries, a.start_indices, a.lengths + ); + a[[5, 0]] = 5; + println!( + "entries: {:?}, start_indices: {:?}, lengths: {:?}", + a.entries, a.start_indices, a.lengths + ); + assert_eq!( + a.indexed_iter().collect::>(), + &[ + ([2, 2], 1), + ([2, 4], 2), + ([4, 1], 3), + ([4, 4], 4), + ([5, 0], 5), + ] + ); + } + + fn index_access() { + let mut array = PackedArray::new([40, 50, 50]); + + // after creation the array must be empty + assert_eq!(array.overhead(), 2); + assert!(array.is_empty()); + + // insert the first element + array[[5, 10, 10]] = 1.0; + assert_eq!(array[[5, 10, 10]], 1.0); + assert_eq!(array.non_zeros(), 1); + assert_eq!(array.explicit_zeros(), 0); + assert_eq!(array.overhead(), 102); + assert!(!array.is_empty()); + + // insert an element after the first one + array[[8, 10, 10]] = 2.0; + assert_eq!(array[[8, 10, 10]], 2.0); + assert_eq!(array[[5, 10, 10]], 1.0); + assert_eq!(array.non_zeros(), 2); + assert_eq!(array.explicit_zeros(), 0); + assert_eq!(array.overhead(), 402); + assert!(!array.is_empty()); + + // insert an element before the first one + array[[1, 10, 10]] = 3.0; + assert_eq!(array[[1, 10, 10]], 3.0); + assert_eq!(array[[8, 10, 10]], 2.0); + assert_eq!(array[[5, 10, 10]], 1.0); + assert_eq!(array.non_zeros(), 3); + assert_eq!(array.explicit_zeros(), 0); + assert_eq!(array.overhead(), 802); + assert!(!array.is_empty()); + + array[[1, 10, 11]] = 4.0; + assert_eq!(array[[1, 10, 11]], 4.0); + assert_eq!(array[[1, 10, 10]], 3.0); + assert_eq!(array[[8, 10, 10]], 2.0); + assert_eq!(array[[5, 10, 10]], 1.0); + assert_eq!(array.non_zeros(), 4); + assert_eq!(array.explicit_zeros(), 0); + assert_eq!(array.overhead(), 802); + assert!(!array.is_empty()); + + array[[1, 10, 9]] = 5.0; + assert_eq!(array[[1, 10, 9]], 5.0); + assert_eq!(array[[1, 10, 11]], 4.0); + assert_eq!(array[[1, 10, 10]], 3.0); + assert_eq!(array[[8, 10, 10]], 2.0); + assert_eq!(array[[5, 10, 10]], 1.0); + assert_eq!(array.non_zeros(), 5); + assert_eq!(array.explicit_zeros(), 0); + assert_eq!(array.overhead(), 802); + assert!(!array.is_empty()); + + array[[1, 10, 0]] = 6.0; + assert_eq!(array[[1, 10, 0]], 6.0); + assert_eq!(array[[1, 10, 9]], 5.0); + assert_eq!(array[[1, 10, 11]], 4.0); + assert_eq!(array[[1, 10, 10]], 3.0); + assert_eq!(array[[8, 10, 10]], 2.0); + assert_eq!(array[[5, 10, 10]], 1.0); + assert_eq!(array.non_zeros(), 6); + assert_eq!(array.overhead(), 802); + assert!(!array.is_empty()); + + // check zeros + assert_eq!(array[[1, 10, 1]], 0.0); + assert_eq!(array[[1, 10, 2]], 0.0); + assert_eq!(array[[1, 10, 3]], 0.0); + assert_eq!(array[[1, 10, 4]], 0.0); + assert_eq!(array[[1, 10, 5]], 0.0); + assert_eq!(array[[1, 10, 6]], 0.0); + assert_eq!(array[[1, 10, 7]], 0.0); + assert_eq!(array[[1, 10, 8]], 0.0); + assert_eq!(array.explicit_zeros(), 8); + + // insert where previously a zero was + array[[1, 10, 2]] = 7.0; + assert_eq!(array[[1, 10, 2]], 7.0); + assert_eq!(array[[1, 10, 0]], 6.0); + assert_eq!(array[[1, 10, 9]], 5.0); + assert_eq!(array[[1, 10, 11]], 4.0); + assert_eq!(array[[1, 10, 10]], 3.0); + assert_eq!(array[[8, 10, 10]], 2.0); + assert_eq!(array[[5, 10, 10]], 1.0); + assert_eq!(array.non_zeros(), 7); + assert_eq!(array.overhead(), 802); + assert!(!array.is_empty()); + + // check zeros + assert_eq!(array[[1, 10, 1]], 0.0); + assert_eq!(array[[1, 10, 3]], 0.0); + assert_eq!(array[[1, 10, 4]], 0.0); + assert_eq!(array[[1, 10, 5]], 0.0); + assert_eq!(array[[1, 10, 6]], 0.0); + assert_eq!(array[[1, 10, 7]], 0.0); + assert_eq!(array[[1, 10, 8]], 0.0); + assert_eq!(array.explicit_zeros(), 7); + + array[[1, 15, 2]] = 8.0; + assert_eq!(array[[1, 15, 2]], 8.0); + assert_eq!(array[[1, 10, 2]], 7.0); + assert_eq!(array[[1, 10, 0]], 6.0); + assert_eq!(array[[1, 10, 9]], 5.0); + assert_eq!(array[[1, 10, 11]], 4.0); + assert_eq!(array[[1, 10, 10]], 3.0); + assert_eq!(array[[8, 10, 10]], 2.0); + assert_eq!(array[[5, 10, 10]], 1.0); + assert_eq!(array.non_zeros(), 8); + assert_eq!(array.overhead(), 802); + assert!(!array.is_empty()); + + // check zeros + assert_eq!(array[[1, 10, 1]], 0.0); + assert_eq!(array[[1, 10, 3]], 0.0); + assert_eq!(array[[1, 10, 4]], 0.0); + assert_eq!(array[[1, 10, 5]], 0.0); + assert_eq!(array[[1, 10, 6]], 0.0); + assert_eq!(array[[1, 10, 7]], 0.0); + assert_eq!(array[[1, 10, 8]], 0.0); + assert_eq!(array.explicit_zeros(), 7); + + array[[1, 15, 4]] = 9.0; + assert_eq!(array[[1, 15, 4]], 9.0); + assert_eq!(array[[1, 15, 2]], 8.0); + assert_eq!(array[[1, 10, 2]], 7.0); + assert_eq!(array[[1, 10, 0]], 6.0); + assert_eq!(array[[1, 10, 9]], 5.0); + assert_eq!(array[[1, 10, 11]], 4.0); + assert_eq!(array[[1, 10, 10]], 3.0); + assert_eq!(array[[8, 10, 10]], 2.0); + assert_eq!(array[[5, 10, 10]], 1.0); + assert_eq!(array.non_zeros(), 9); + assert_eq!(array.overhead(), 802); + assert!(!array.is_empty()); + + // check zeros + assert_eq!(array[[1, 15, 3]], 0.0); + assert_eq!(array[[1, 10, 1]], 0.0); + assert_eq!(array[[1, 10, 3]], 0.0); + assert_eq!(array[[1, 10, 4]], 0.0); + assert_eq!(array[[1, 10, 5]], 0.0); + assert_eq!(array[[1, 10, 6]], 0.0); + assert_eq!(array[[1, 10, 7]], 0.0); + assert_eq!(array[[1, 10, 8]], 0.0); + assert_eq!(array.explicit_zeros(), 8); + + array[[1, 15, 0]] = 10.0; + assert_eq!(array[[1, 15, 0]], 10.0); + assert_eq!(array[[1, 15, 4]], 9.0); + assert_eq!(array[[1, 15, 2]], 8.0); + assert_eq!(array[[1, 10, 2]], 7.0); + assert_eq!(array[[1, 10, 0]], 6.0); + assert_eq!(array[[1, 10, 9]], 5.0); + assert_eq!(array[[1, 10, 11]], 4.0); + assert_eq!(array[[1, 10, 10]], 3.0); + assert_eq!(array[[8, 10, 10]], 2.0); + assert_eq!(array[[5, 10, 10]], 1.0); + assert_eq!(array.non_zeros(), 10); + assert_eq!(array.overhead(), 802); + assert!(!array.is_empty()); + + // check zeros + assert_eq!(array[[1, 15, 1]], 0.0); + assert_eq!(array[[1, 15, 3]], 0.0); + assert_eq!(array[[1, 10, 1]], 0.0); + assert_eq!(array[[1, 10, 3]], 0.0); + assert_eq!(array[[1, 10, 4]], 0.0); + assert_eq!(array[[1, 10, 5]], 0.0); + assert_eq!(array[[1, 10, 6]], 0.0); + assert_eq!(array[[1, 10, 7]], 0.0); + assert_eq!(array[[1, 10, 8]], 0.0); + assert_eq!(array.explicit_zeros(), 9); + } + + #[test] + #[should_panic(expected = "index [40, 0, 50] is out of bounds for array of shape [40, 50, 50]")] + fn index_mut_panic_dim0() { + let mut array = PackedArray::new([40, 50, 50]); + + array[[40, 0, 50]] = 1.0; + } + + #[test] + #[should_panic(expected = "index [0, 50, 0] is out of bounds for array of shape [40, 50, 50]")] + fn index_mut_panic_dim1() { + let mut array = PackedArray::new([40, 50, 50]); + + array[[0, 50, 0]] = 1.0; + } + + #[test] + #[should_panic(expected = "index [0, 0, 50] is out of bounds for array of shape [40, 50, 50]")] + fn index_mut_panic_dim2() { + let mut array = PackedArray::new([40, 50, 50]); + + array[[0, 0, 50]] = 1.0; + } + + #[test] + #[should_panic(expected = "entry at index [0, 0, 0] is implicitly set to the default value")] + fn index_panic_dim0_0() { + let mut array = PackedArray::new([40, 50, 50]); + + array[[1, 0, 0]] = 1.0; + + assert_eq!(array[[0, 0, 0]], 0.0); + } + + #[test] + #[should_panic(expected = "entry at index [2, 0, 0] is implicitly set to the default value")] + fn index_panic_dim0_1() { + let mut array = PackedArray::new([40, 50, 50]); + + array[[1, 0, 0]] = 1.0; + + assert_eq!(array[[2, 0, 0]], 0.0); + } + + #[test] + #[should_panic(expected = "index [1, 50, 0] is out of bounds for array of shape [40, 50, 50]")] + fn index_panic_dim1() { + let mut array = PackedArray::new([40, 50, 50]); + + array[[1, 0, 0]] = 1.0; + + assert_eq!(array[[1, 50, 0]], 0.0); + } + + #[test] + #[should_panic(expected = "entry at index [0, 0, 0] is implicitly set to the default value")] + fn index_panic_dim2_0() { + let mut array = PackedArray::new([40, 50, 50]); + + array[[0, 0, 1]] = 1.0; + + assert_eq!(array[[0, 0, 0]], 0.0); + } + + #[test] + #[should_panic(expected = "entry at index [0, 0, 2] is implicitly set to the default value")] + fn index_panic_dim2_1() { + let mut array = PackedArray::new([40, 50, 50]); + + array[[0, 0, 1]] = 1.0; + + assert_eq!(array[[0, 0, 2]], 0.0); + } + + #[test] + fn indexed_iter() { + let mut array = PackedArray::new([40, 50, 50]); + + // check empty iterator + assert_eq!(array.indexed_iter().next(), None); + + // insert an element + array[[2, 3, 4]] = 1.0; + + let mut iter = array.indexed_iter(); + + // check iterator with one element + assert_eq!(iter.next(), Some(([2, 3, 4], 1.0))); + assert_eq!(iter.next(), None); + + mem::drop(iter); + + // insert another element + array[[2, 3, 6]] = 2.0; + + let mut iter = array.indexed_iter(); + + assert_eq!(iter.next(), Some(([2, 3, 4], 1.0))); + assert_eq!(iter.next(), Some(([2, 3, 6], 2.0))); + assert_eq!(iter.next(), None); + + mem::drop(iter); + + // insert yet another element + array[[4, 5, 7]] = 3.0; + + let mut iter = array.indexed_iter(); + + assert_eq!(iter.next(), Some(([2, 3, 4], 1.0))); + assert_eq!(iter.next(), Some(([2, 3, 6], 2.0))); + assert_eq!(iter.next(), Some(([4, 5, 7], 3.0))); + assert_eq!(iter.next(), None); + + mem::drop(iter); + + // insert at the very first position + array[[2, 0, 0]] = 4.0; + + let mut iter = array.indexed_iter(); + + assert_eq!(iter.next(), Some(([2, 0, 0], 4.0))); + assert_eq!(iter.next(), Some(([2, 3, 4], 1.0))); + assert_eq!(iter.next(), Some(([2, 3, 6], 2.0))); + assert_eq!(iter.next(), Some(([4, 5, 7], 3.0))); + assert_eq!(iter.next(), None); + } + + #[test] + fn clear() { + let mut array = PackedArray::new([40, 50, 50]); + + array[[3, 5, 1]] = 1.0; + array[[7, 8, 9]] = 2.0; + array[[9, 1, 4]] = 3.0; + + assert!(!array.is_empty()); + assert_eq!(array.non_zeros(), 3); + assert_eq!(array.explicit_zeros(), 0); + + array.clear(); + + assert!(array.is_empty()); + assert_eq!(array.non_zeros(), 0); + assert_eq!(array.explicit_zeros(), 0); + } + + #[test] + fn from_ndarray() { + let mut ndarray = Array3::zeros((2, 50, 50)); + + ndarray[[0, 4, 3]] = 1.0; + ndarray[[0, 4, 4]] = 2.0; + ndarray[[0, 4, 6]] = 3.0; + ndarray[[0, 5, 1]] = 4.0; + ndarray[[0, 5, 7]] = 5.0; + ndarray[[1, 3, 9]] = 6.0; + + let array = PackedArray::from_ndarray(ndarray.view(), 3, 40); + + assert_eq!(array[[3, 4, 3]], 1.0); + assert_eq!(array[[3, 4, 4]], 2.0); + // assert_eq!(array[[3, 4, 5]], 0.0); + assert_eq!(array[[3, 4, 6]], 3.0); + assert_eq!(array[[3, 5, 1]], 4.0); + // assert_eq!(array[[3, 5, 2]], 0.0); + // assert_eq!(array[[3, 5, 3]], 0.0); + // assert_eq!(array[[3, 5, 4]], 0.0); + // assert_eq!(array[[3, 5, 5]], 0.0); + // assert_eq!(array[[3, 5, 6]], 0.0); + assert_eq!(array[[3, 5, 7]], 5.0); + assert_eq!(array[[4, 3, 9]], 6.0); + + // assert_eq!(array.len(), 6); + // assert_eq!(array.zeros(), 6); + } + + // #[test] + // fn test_index_swap() { + // let mut array = JaggedArray::new([5, 50, 2]); + + // array[[0, 0, 0]] = 1.0; + // array[[0, 0, 1]] = 2.0; + // array[[1, 2, 1]] = 3.0; + // array[[1, 5, 1]] = 4.0; + // array[[1, 6, 0]] = 5.0; + // array[[1, 8, 0]] = 6.0; + // array[[1, 9, 0]] = 7.0; + // array[[2, 0, 0]] = 8.0; + // array[[3, 2, 1]] = 9.0; + // array[[3, 4, 0]] = 10.0; + // array[[3, 4, 1]] = 11.0; + // array[[4, 0, 0]] = 12.0; + // array[[4, 0, 1]] = 13.0; + + // assert_eq!(array[[0, 0, 0]], 1.0); + // assert_eq!(array[[0, 0, 1]], 2.0); + // assert_eq!(array[[1, 2, 1]], 3.0); + // assert_eq!(array[[1, 5, 1]], 4.0); + // assert_eq!(array[[1, 6, 0]], 5.0); + // assert_eq!(array[[1, 8, 0]], 6.0); + // assert_eq!(array[[1, 9, 0]], 7.0); + // assert_eq!(array[[2, 0, 0]], 8.0); + // assert_eq!(array[[3, 2, 1]], 9.0); + // assert_eq!(array[[3, 4, 0]], 10.0); + // assert_eq!(array[[3, 4, 1]], 11.0); + // assert_eq!(array[[4, 0, 0]], 12.0); + // assert_eq!(array[[4, 0, 1]], 13.0); + + // let mut iter = array.indexed_iter(); + + // assert_eq!(iter.next(), Some(([0, 0, 0], 1.0))); + // assert_eq!(iter.next(), Some(([0, 0, 1], 2.0))); + // assert_eq!(iter.next(), Some(([1, 6, 0], 5.0))); + // assert_eq!(iter.next(), Some(([1, 8, 0], 6.0))); + // assert_eq!(iter.next(), Some(([1, 9, 0], 7.0))); + // assert_eq!(iter.next(), Some(([1, 2, 1], 3.0))); + // assert_eq!(iter.next(), Some(([1, 5, 1], 4.0))); + // assert_eq!(iter.next(), Some(([2, 0, 0], 8.0))); + // assert_eq!(iter.next(), Some(([3, 4, 0], 10.0))); + // assert_eq!(iter.next(), Some(([3, 2, 1], 9.0))); + // assert_eq!(iter.next(), Some(([3, 4, 1], 11.0))); + // assert_eq!(iter.next(), Some(([4, 0, 0], 12.0))); + // assert_eq!(iter.next(), Some(([4, 0, 1], 13.0))); + // assert_eq!(iter.next(), None); + + // let mut ndarray = Array3::zeros((5, 50, 2)); + + // ndarray[[0, 0, 0]] = 1.0; + // ndarray[[0, 0, 1]] = 2.0; + // ndarray[[1, 2, 1]] = 3.0; + // ndarray[[1, 5, 1]] = 4.0; + // ndarray[[1, 6, 0]] = 5.0; + // ndarray[[1, 8, 0]] = 6.0; + // ndarray[[1, 9, 0]] = 7.0; + // ndarray[[2, 0, 0]] = 8.0; + // ndarray[[3, 2, 1]] = 9.0; + // ndarray[[3, 4, 0]] = 10.0; + // ndarray[[3, 4, 1]] = 11.0; + // ndarray[[4, 0, 0]] = 12.0; + // ndarray[[4, 0, 1]] = 13.0; + + // let mut other = JaggedArray::from_ndarray(ndarray.view(), 0, 5); + + // assert_eq!(other[[0, 0, 0]], 1.0); + // assert_eq!(other[[0, 0, 1]], 2.0); + // assert_eq!(other[[1, 2, 1]], 3.0); + // assert_eq!(other[[1, 5, 1]], 4.0); + // assert_eq!(other[[1, 6, 0]], 5.0); + // assert_eq!(other[[1, 8, 0]], 6.0); + // assert_eq!(other[[1, 9, 0]], 7.0); + // assert_eq!(other[[2, 0, 0]], 8.0); + // assert_eq!(other[[3, 2, 1]], 9.0); + // assert_eq!(other[[3, 4, 0]], 10.0); + // assert_eq!(other[[3, 4, 1]], 11.0); + // assert_eq!(other[[4, 0, 0]], 12.0); + // assert_eq!(other[[4, 0, 1]], 13.0); + // } +} From e94437ec717be9786a32d0bccab919f67723c6fc Mon Sep 17 00:00:00 2001 From: janw20 Date: Thu, 18 Apr 2024 14:30:23 +0200 Subject: [PATCH 02/14] Fix `index_access` test and one edge case --- pineappl/src/packed_array.rs | 175 +++++++---------------------------- 1 file changed, 32 insertions(+), 143 deletions(-) diff --git a/pineappl/src/packed_array.rs b/pineappl/src/packed_array.rs index 061666563..bb2a32626 100644 --- a/pineappl/src/packed_array.rs +++ b/pineappl/src/packed_array.rs @@ -215,18 +215,20 @@ impl IndexMut<[usize; D]> } return &mut self.entries[point_entries - 1 + distance]; - } else if let Some(start_index_next) = self.start_indices.get(point) { - if raveled_index + threshold_distance >= *start_index_next { - // println!("test 5"); - let distance = start_index_next - raveled_index; - - self.start_indices[point] = raveled_index; - self.lengths[point] += distance; - for _ in 0..distance { - self.entries.insert(point_entries, Default::default()); - } - return &mut self.entries[point_entries]; + } + } + + if let Some(start_index_next) = self.start_indices.get(point) { + if raveled_index + threshold_distance >= *start_index_next { + // println!("test 5"); + let distance = start_index_next - raveled_index; + + self.start_indices[point] = raveled_index; + self.lengths[point] += distance; + for _ in 0..distance { + self.entries.insert(point_entries, Default::default()); } + return &mut self.entries[point_entries]; } } @@ -389,11 +391,12 @@ mod tests { ); } + #[test] fn index_access() { let mut array = PackedArray::new([40, 50, 50]); // after creation the array must be empty - assert_eq!(array.overhead(), 2); + assert_eq!(array.overhead(), 0); assert!(array.is_empty()); // insert the first element @@ -401,7 +404,7 @@ mod tests { assert_eq!(array[[5, 10, 10]], 1.0); assert_eq!(array.non_zeros(), 1); assert_eq!(array.explicit_zeros(), 0); - assert_eq!(array.overhead(), 102); + assert_eq!(array.overhead(), 2); assert!(!array.is_empty()); // insert an element after the first one @@ -410,7 +413,7 @@ mod tests { assert_eq!(array[[5, 10, 10]], 1.0); assert_eq!(array.non_zeros(), 2); assert_eq!(array.explicit_zeros(), 0); - assert_eq!(array.overhead(), 402); + assert_eq!(array.overhead(), 4); assert!(!array.is_empty()); // insert an element before the first one @@ -420,7 +423,7 @@ mod tests { assert_eq!(array[[5, 10, 10]], 1.0); assert_eq!(array.non_zeros(), 3); assert_eq!(array.explicit_zeros(), 0); - assert_eq!(array.overhead(), 802); + assert_eq!(array.overhead(), 6); assert!(!array.is_empty()); array[[1, 10, 11]] = 4.0; @@ -430,7 +433,7 @@ mod tests { assert_eq!(array[[5, 10, 10]], 1.0); assert_eq!(array.non_zeros(), 4); assert_eq!(array.explicit_zeros(), 0); - assert_eq!(array.overhead(), 802); + assert_eq!(array.overhead(), 6); assert!(!array.is_empty()); array[[1, 10, 9]] = 5.0; @@ -441,7 +444,9 @@ mod tests { assert_eq!(array[[5, 10, 10]], 1.0); assert_eq!(array.non_zeros(), 5); assert_eq!(array.explicit_zeros(), 0); - assert_eq!(array.overhead(), 802); + // dbg!(&array.start_indices); + // dbg!(&array.lengths); + assert_eq!(array.overhead(), 6); assert!(!array.is_empty()); array[[1, 10, 0]] = 6.0; @@ -452,21 +457,10 @@ mod tests { assert_eq!(array[[8, 10, 10]], 2.0); assert_eq!(array[[5, 10, 10]], 1.0); assert_eq!(array.non_zeros(), 6); - assert_eq!(array.overhead(), 802); + assert_eq!(array.explicit_zeros(), 0); + assert_eq!(array.overhead(), 8); assert!(!array.is_empty()); - // check zeros - assert_eq!(array[[1, 10, 1]], 0.0); - assert_eq!(array[[1, 10, 2]], 0.0); - assert_eq!(array[[1, 10, 3]], 0.0); - assert_eq!(array[[1, 10, 4]], 0.0); - assert_eq!(array[[1, 10, 5]], 0.0); - assert_eq!(array[[1, 10, 6]], 0.0); - assert_eq!(array[[1, 10, 7]], 0.0); - assert_eq!(array[[1, 10, 8]], 0.0); - assert_eq!(array.explicit_zeros(), 8); - - // insert where previously a zero was array[[1, 10, 2]] = 7.0; assert_eq!(array[[1, 10, 2]], 7.0); assert_eq!(array[[1, 10, 0]], 6.0); @@ -476,18 +470,12 @@ mod tests { assert_eq!(array[[8, 10, 10]], 2.0); assert_eq!(array[[5, 10, 10]], 1.0); assert_eq!(array.non_zeros(), 7); - assert_eq!(array.overhead(), 802); + assert_eq!(array.overhead(), 8); assert!(!array.is_empty()); // check zeros assert_eq!(array[[1, 10, 1]], 0.0); - assert_eq!(array[[1, 10, 3]], 0.0); - assert_eq!(array[[1, 10, 4]], 0.0); - assert_eq!(array[[1, 10, 5]], 0.0); - assert_eq!(array[[1, 10, 6]], 0.0); - assert_eq!(array[[1, 10, 7]], 0.0); - assert_eq!(array[[1, 10, 8]], 0.0); - assert_eq!(array.explicit_zeros(), 7); + assert_eq!(array.explicit_zeros(), 1); array[[1, 15, 2]] = 8.0; assert_eq!(array[[1, 15, 2]], 8.0); @@ -499,18 +487,12 @@ mod tests { assert_eq!(array[[8, 10, 10]], 2.0); assert_eq!(array[[5, 10, 10]], 1.0); assert_eq!(array.non_zeros(), 8); - assert_eq!(array.overhead(), 802); + assert_eq!(array.overhead(), 10); assert!(!array.is_empty()); // check zeros assert_eq!(array[[1, 10, 1]], 0.0); - assert_eq!(array[[1, 10, 3]], 0.0); - assert_eq!(array[[1, 10, 4]], 0.0); - assert_eq!(array[[1, 10, 5]], 0.0); - assert_eq!(array[[1, 10, 6]], 0.0); - assert_eq!(array[[1, 10, 7]], 0.0); - assert_eq!(array[[1, 10, 8]], 0.0); - assert_eq!(array.explicit_zeros(), 7); + assert_eq!(array.explicit_zeros(), 1); array[[1, 15, 4]] = 9.0; assert_eq!(array[[1, 15, 4]], 9.0); @@ -523,19 +505,13 @@ mod tests { assert_eq!(array[[8, 10, 10]], 2.0); assert_eq!(array[[5, 10, 10]], 1.0); assert_eq!(array.non_zeros(), 9); - assert_eq!(array.overhead(), 802); + assert_eq!(array.overhead(), 10); assert!(!array.is_empty()); // check zeros assert_eq!(array[[1, 15, 3]], 0.0); assert_eq!(array[[1, 10, 1]], 0.0); - assert_eq!(array[[1, 10, 3]], 0.0); - assert_eq!(array[[1, 10, 4]], 0.0); - assert_eq!(array[[1, 10, 5]], 0.0); - assert_eq!(array[[1, 10, 6]], 0.0); - assert_eq!(array[[1, 10, 7]], 0.0); - assert_eq!(array[[1, 10, 8]], 0.0); - assert_eq!(array.explicit_zeros(), 8); + assert_eq!(array.explicit_zeros(), 2); array[[1, 15, 0]] = 10.0; assert_eq!(array[[1, 15, 0]], 10.0); @@ -549,20 +525,14 @@ mod tests { assert_eq!(array[[8, 10, 10]], 2.0); assert_eq!(array[[5, 10, 10]], 1.0); assert_eq!(array.non_zeros(), 10); - assert_eq!(array.overhead(), 802); + assert_eq!(array.overhead(), 10); assert!(!array.is_empty()); // check zeros assert_eq!(array[[1, 15, 1]], 0.0); assert_eq!(array[[1, 15, 3]], 0.0); assert_eq!(array[[1, 10, 1]], 0.0); - assert_eq!(array[[1, 10, 3]], 0.0); - assert_eq!(array[[1, 10, 4]], 0.0); - assert_eq!(array[[1, 10, 5]], 0.0); - assert_eq!(array[[1, 10, 6]], 0.0); - assert_eq!(array[[1, 10, 7]], 0.0); - assert_eq!(array[[1, 10, 8]], 0.0); - assert_eq!(array.explicit_zeros(), 9); + assert_eq!(array.explicit_zeros(), 3); } #[test] @@ -741,85 +711,4 @@ mod tests { // assert_eq!(array.zeros(), 6); } - // #[test] - // fn test_index_swap() { - // let mut array = JaggedArray::new([5, 50, 2]); - - // array[[0, 0, 0]] = 1.0; - // array[[0, 0, 1]] = 2.0; - // array[[1, 2, 1]] = 3.0; - // array[[1, 5, 1]] = 4.0; - // array[[1, 6, 0]] = 5.0; - // array[[1, 8, 0]] = 6.0; - // array[[1, 9, 0]] = 7.0; - // array[[2, 0, 0]] = 8.0; - // array[[3, 2, 1]] = 9.0; - // array[[3, 4, 0]] = 10.0; - // array[[3, 4, 1]] = 11.0; - // array[[4, 0, 0]] = 12.0; - // array[[4, 0, 1]] = 13.0; - - // assert_eq!(array[[0, 0, 0]], 1.0); - // assert_eq!(array[[0, 0, 1]], 2.0); - // assert_eq!(array[[1, 2, 1]], 3.0); - // assert_eq!(array[[1, 5, 1]], 4.0); - // assert_eq!(array[[1, 6, 0]], 5.0); - // assert_eq!(array[[1, 8, 0]], 6.0); - // assert_eq!(array[[1, 9, 0]], 7.0); - // assert_eq!(array[[2, 0, 0]], 8.0); - // assert_eq!(array[[3, 2, 1]], 9.0); - // assert_eq!(array[[3, 4, 0]], 10.0); - // assert_eq!(array[[3, 4, 1]], 11.0); - // assert_eq!(array[[4, 0, 0]], 12.0); - // assert_eq!(array[[4, 0, 1]], 13.0); - - // let mut iter = array.indexed_iter(); - - // assert_eq!(iter.next(), Some(([0, 0, 0], 1.0))); - // assert_eq!(iter.next(), Some(([0, 0, 1], 2.0))); - // assert_eq!(iter.next(), Some(([1, 6, 0], 5.0))); - // assert_eq!(iter.next(), Some(([1, 8, 0], 6.0))); - // assert_eq!(iter.next(), Some(([1, 9, 0], 7.0))); - // assert_eq!(iter.next(), Some(([1, 2, 1], 3.0))); - // assert_eq!(iter.next(), Some(([1, 5, 1], 4.0))); - // assert_eq!(iter.next(), Some(([2, 0, 0], 8.0))); - // assert_eq!(iter.next(), Some(([3, 4, 0], 10.0))); - // assert_eq!(iter.next(), Some(([3, 2, 1], 9.0))); - // assert_eq!(iter.next(), Some(([3, 4, 1], 11.0))); - // assert_eq!(iter.next(), Some(([4, 0, 0], 12.0))); - // assert_eq!(iter.next(), Some(([4, 0, 1], 13.0))); - // assert_eq!(iter.next(), None); - - // let mut ndarray = Array3::zeros((5, 50, 2)); - - // ndarray[[0, 0, 0]] = 1.0; - // ndarray[[0, 0, 1]] = 2.0; - // ndarray[[1, 2, 1]] = 3.0; - // ndarray[[1, 5, 1]] = 4.0; - // ndarray[[1, 6, 0]] = 5.0; - // ndarray[[1, 8, 0]] = 6.0; - // ndarray[[1, 9, 0]] = 7.0; - // ndarray[[2, 0, 0]] = 8.0; - // ndarray[[3, 2, 1]] = 9.0; - // ndarray[[3, 4, 0]] = 10.0; - // ndarray[[3, 4, 1]] = 11.0; - // ndarray[[4, 0, 0]] = 12.0; - // ndarray[[4, 0, 1]] = 13.0; - - // let mut other = JaggedArray::from_ndarray(ndarray.view(), 0, 5); - - // assert_eq!(other[[0, 0, 0]], 1.0); - // assert_eq!(other[[0, 0, 1]], 2.0); - // assert_eq!(other[[1, 2, 1]], 3.0); - // assert_eq!(other[[1, 5, 1]], 4.0); - // assert_eq!(other[[1, 6, 0]], 5.0); - // assert_eq!(other[[1, 8, 0]], 6.0); - // assert_eq!(other[[1, 9, 0]], 7.0); - // assert_eq!(other[[2, 0, 0]], 8.0); - // assert_eq!(other[[3, 2, 1]], 9.0); - // assert_eq!(other[[3, 4, 0]], 10.0); - // assert_eq!(other[[3, 4, 1]], 11.0); - // assert_eq!(other[[4, 0, 0]], 12.0); - // assert_eq!(other[[4, 0, 1]], 13.0); - // } } From bbd50966c74c5a6be5f609a1f2cf779e420b58c3 Mon Sep 17 00:00:00 2001 From: Christopher Schwan Date: Thu, 18 Apr 2024 11:48:56 +0200 Subject: [PATCH 03/14] Add missing `must_use` attributes --- pineappl/src/packed_array.rs | 3 +++ 1 file changed, 3 insertions(+) diff --git a/pineappl/src/packed_array.rs b/pineappl/src/packed_array.rs index bb2a32626..f0b119f59 100644 --- a/pineappl/src/packed_array.rs +++ b/pineappl/src/packed_array.rs @@ -29,6 +29,7 @@ impl PackedArray { self.entries.is_empty() } + #[must_use] pub fn shape(&self) -> &[usize] { &self.shape } @@ -50,6 +51,7 @@ impl PackedArray { self.entries.iter().filter(|x| **x == T::default()).count() } + #[must_use] pub fn non_zeros(&self) -> usize { self.entries.iter().filter(|x| **x != T::default()).count() } @@ -74,6 +76,7 @@ impl, const D: usize> MulAssign for PackedArray } impl PackedArray { + #[must_use] pub fn from_ndarray(array: ArrayView3, xstart: usize, xsize: usize) -> Self { let shape = array.shape(); From 9023fa71600cb2f9935ba4e1f6e494bc090897c0 Mon Sep 17 00:00:00 2001 From: Christopher Schwan Date: Thu, 18 Apr 2024 11:52:44 +0200 Subject: [PATCH 04/14] Avoid having to specify dimension `D` explicitly --- pineappl/src/packed_array.rs | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/pineappl/src/packed_array.rs b/pineappl/src/packed_array.rs index f0b119f59..bc874617d 100644 --- a/pineappl/src/packed_array.rs +++ b/pineappl/src/packed_array.rs @@ -93,7 +93,7 @@ impl PackedArray { } } -fn ravel_multi_index(indices: &[usize], dimensions: &[usize]) -> usize { +fn ravel_multi_index(indices: &[usize; D], dimensions: &[usize]) -> usize { assert_eq!(indices.len(), dimensions.len()); indices @@ -129,7 +129,7 @@ impl Index<[usize; D]> for Packed self.shape ); - let raveled_index = ravel_multi_index::(&index, &self.shape); + let raveled_index = ravel_multi_index(&index, &self.shape); let point = self.start_indices.partition_point(|&i| i <= raveled_index); assert!( @@ -165,7 +165,7 @@ impl IndexMut<[usize; D]> self.shape ); - let raveled_index = ravel_multi_index::(&index, &self.shape); + let raveled_index = ravel_multi_index(&index, &self.shape); // if self.entries.is_empty() { // self.start_indices.push(raveled_index); @@ -274,12 +274,12 @@ mod tests { #[test] fn ravel_multi_index() { - assert_eq!(super::ravel_multi_index::<2>(&[0, 0], &[3, 2]), 0); - assert_eq!(super::ravel_multi_index::<2>(&[0, 1], &[3, 2]), 1); - assert_eq!(super::ravel_multi_index::<2>(&[1, 0], &[3, 2]), 2); - assert_eq!(super::ravel_multi_index::<2>(&[1, 1], &[3, 2]), 3); - assert_eq!(super::ravel_multi_index::<2>(&[2, 0], &[3, 2]), 4); - assert_eq!(super::ravel_multi_index::<2>(&[2, 1], &[3, 2]), 5); + assert_eq!(super::ravel_multi_index(&[0, 0], &[3, 2]), 0); + assert_eq!(super::ravel_multi_index(&[0, 1], &[3, 2]), 1); + assert_eq!(super::ravel_multi_index(&[1, 0], &[3, 2]), 2); + assert_eq!(super::ravel_multi_index(&[1, 1], &[3, 2]), 3); + assert_eq!(super::ravel_multi_index(&[2, 0], &[3, 2]), 4); + assert_eq!(super::ravel_multi_index(&[2, 1], &[3, 2]), 5); } #[test] From 5d40cc75a150f924d1a05e3b721e6576ccfb3dd3 Mon Sep 17 00:00:00 2001 From: janw20 Date: Thu, 18 Apr 2024 15:18:23 +0200 Subject: [PATCH 05/14] Add documentation and do cleanup --- pineappl/src/packed_array.rs | 110 ++++++++++------------------------- 1 file changed, 30 insertions(+), 80 deletions(-) diff --git a/pineappl/src/packed_array.rs b/pineappl/src/packed_array.rs index bc874617d..6cb1ca2ac 100644 --- a/pineappl/src/packed_array.rs +++ b/pineappl/src/packed_array.rs @@ -4,15 +4,21 @@ use std::{mem, vec}; use ndarray::ArrayView3; use serde::{Deserialize, Serialize}; +/// `D`-dimensional array similar to ndarray::ArrayBase, except that `T::default()` is not stored to save space. Instead, adjacent elements are grouped together and only the index of their first element (`start_index`) and the length of the group (`lengths`) is stored. #[derive(Clone, Deserialize, Serialize)] pub struct PackedArray { + /// The actual values stored in the array. The length of `entries` is always the sum of the elements in `lengths`. entries: Vec, + /// The indices of the first elements in each group. start_indices: Vec, + /// The length of each group. lengths: Vec, + /// The shape (dimensions) of the array. shape: Vec, } impl PackedArray { + /// Constructs a new and empty `PackedArray` of shape `shape`. #[must_use] pub fn new(shape: [usize; D]) -> Self { Self { @@ -24,38 +30,45 @@ impl PackedArray { } } + /// Returns `true` if the array contains no element. #[must_use] pub fn is_empty(&self) -> bool { self.entries.is_empty() } + /// Returns the shape of the array. #[must_use] pub fn shape(&self) -> &[usize] { &self.shape } + /// Clears the contents of the array. pub fn clear(&mut self) { self.entries.clear(); self.start_indices.clear(); self.lengths.clear(); } + /// Returns the overhead of storing the `start_indices` and the `lengths` of the groups, in units of `f64`. #[must_use] pub fn overhead(&self) -> usize { ((self.start_indices.len() + self.lengths.len()) * mem::size_of::()) / mem::size_of::() } + /// Returns the number of default (zero) elements that are explicitly stored in `entries`. If there is one default elements between adjacent groups, it is more economical to store the one default element explicitly and merge the two groups, than to store the `start_indices` and `lengths` of those groups. #[must_use] pub fn explicit_zeros(&self) -> usize { self.entries.iter().filter(|x| **x == T::default()).count() } + /// Returns the number of non-default (non-zero) elements stored in the array. #[must_use] pub fn non_zeros(&self) -> usize { self.entries.iter().filter(|x| **x != T::default()).count() } + /// Returns an `Iterator` over the non-default (non-zero) elements of this array. The type of an iterator element is `([usize; D], T)`. pub fn indexed_iter(&self) -> impl Iterator + '_ { self.start_indices .iter() @@ -70,12 +83,16 @@ impl PackedArray { } impl, const D: usize> MulAssign for PackedArray { + + /// Perform `self *= rhs` as elementwise multiplication (in place). fn mul_assign(&mut self, rhs: T) { self.entries.iter_mut().for_each(|x| *x *= rhs); } } impl PackedArray { + + /// Converts `array` into a `PackedArray` of dimension 3. #[must_use] pub fn from_ndarray(array: ArrayView3, xstart: usize, xsize: usize) -> Self { let shape = array.shape(); @@ -93,16 +110,18 @@ impl PackedArray { } } -fn ravel_multi_index(indices: &[usize; D], dimensions: &[usize]) -> usize { - assert_eq!(indices.len(), dimensions.len()); - - indices - .iter() - .skip(1) - .zip(dimensions.iter().skip(1)) - .fold(indices[0], |acc, (i, d)| acc * d + i) +/// Converts a `multi_index` into an flat index. +fn ravel_multi_index(multi_index: &[usize; D], dimensions: &[usize]) -> usize { + assert_eq!(multi_index.len(), dimensions.len()); + + multi_index + .iter() + .skip(1) + .zip(dimensions.iter().skip(1)) + .fold(multi_index[0], |acc, (i, d)| acc * d + i) } +/// Converts a flat `index` into a multi_index. fn unravel_index(index: usize, dimensions: &[usize]) -> [usize; D] { assert!(index < dimensions.iter().product()); let mut indices = [0; D]; @@ -152,7 +171,6 @@ impl Index<[usize; D]> for Packed } } -// TODO: implement axis swapping optimization impl IndexMut<[usize; D]> for PackedArray { @@ -167,22 +185,12 @@ impl IndexMut<[usize; D]> let raveled_index = ravel_multi_index(&index, &self.shape); - // if self.entries.is_empty() { - // self.start_indices.push(raveled_index); - // self.lengths.push(1); - // self.entries.push(Default::default()); - - // return &mut self.entries[0]; - // } - // the closest start_index that is greater than raveled_index let point = self.start_indices.partition_point(|&i| i <= raveled_index); // the index that is stored in `point`, translated to the indices of the entries let point_entries = self.lengths.iter().take(point).sum::(); - // println!("test 1"); - // maximum distance for merging regions let threshold_distance = 2; @@ -191,14 +199,11 @@ impl IndexMut<[usize; D]> let length = self.lengths[point - 1]; if raveled_index < start_index + length { - // println!("test 2"); return &mut self.entries[point_entries - length + raveled_index - start_index]; } else if raveled_index < start_index + length + threshold_distance { - // println!("test 3"); - let distance = raveled_index - (start_index + length) + 1; - // println!("distance: {}", distance); self.lengths[point - 1] += distance; + for _ in 0..(distance) { self.entries.insert(point_entries, Default::default()); } @@ -223,7 +228,6 @@ impl IndexMut<[usize; D]> if let Some(start_index_next) = self.start_indices.get(point) { if raveled_index + threshold_distance >= *start_index_next { - // println!("test 5"); let distance = start_index_next - raveled_index; self.start_indices[point] = raveled_index; @@ -235,8 +239,6 @@ impl IndexMut<[usize; D]> } } - // println!("test 6"); - // in the most general case we have to insert a new start_index self.start_indices.insert(point, raveled_index); self.lengths.insert(point, 1); @@ -246,14 +248,6 @@ impl IndexMut<[usize; D]> } } -// impl IntoIterator for &JaggedArray { -// type Item = (Vec, T); -// type IntoIter = Iterator, T)>; - -// fn into_iter(self) -> Self::IntoIter { - -// } -// } #[cfg(test)] mod tests { @@ -287,20 +281,12 @@ mod tests { let mut a = PackedArray::::new([4, 2]); a[[0, 0]] = 1.0; - println!( - "entries: {:?}, start_indices: {:?}, lengths: {:?}", - a.entries, a.start_indices, a.lengths - ); assert_eq!(a[[0, 0]], 1.0); assert_eq!(a.entries, vec![1.0]); assert_eq!(a.start_indices, vec![0]); assert_eq!(a.lengths, vec![1]); a[[3, 0]] = 2.0; - println!( - "entries: {:?}, start_indices: {:?}, lengths: {:?}", - a.entries, a.start_indices, a.lengths - ); assert_eq!(a[[0, 0]], 1.0); assert_eq!(a[[3, 0]], 2.0); assert_eq!(a.entries, vec![1.0, 2.0]); @@ -308,11 +294,6 @@ mod tests { assert_eq!(a.lengths, vec![1, 1]); a[[3, 1]] = 3.0; - println!( - "entries: {:?}, start_indices: {:?}, lengths: {:?}", - a.entries, a.start_indices, a.lengths - ); - assert_eq!(a[[0, 0]], 1.0); assert_eq!(a[[3, 0]], 2.0); assert_eq!(a[[3, 1]], 3.0); @@ -321,30 +302,19 @@ mod tests { assert_eq!(a.lengths, vec![1, 2]); a[[2, 0]] = 3.5; - println!( - "entries: {:?}, start_indices: {:?}, lengths: {:?}", - a.entries, a.start_indices, a.lengths - ); assert_eq!(a[[0, 0]], 1.0); assert_eq!(a[[3, 0]], 2.0); assert_eq!(a[[3, 1]], 3.0); assert_eq!(a[[2, 0]], 3.5); - // assert_eq!(a[[2, 1]], 0.0); assert_eq!(a.entries, vec![1.0, 3.5, 0.0, 2.0, 3.0]); assert_eq!(a.start_indices, vec![0, 4]); assert_eq!(a.lengths, vec![1, 4]); a[[2, 0]] = 4.0; - println!( - "entries: {:?}, start_indices: {:?}, lengths: {:?}", - a.entries, a.start_indices, a.lengths - ); - assert_eq!(a[[0, 0]], 1.0); assert_eq!(a[[3, 0]], 2.0); assert_eq!(a[[3, 1]], 3.0); assert_eq!(a[[2, 0]], 4.0); - // assert_eq!(a[[2, 1]], 0.0); assert_eq!(a.entries, vec![1.0, 4.0, 0.0, 2.0, 3.0]); assert_eq!(a.start_indices, vec![0, 4]); assert_eq!(a.lengths, vec![1, 4]); @@ -354,16 +324,10 @@ mod tests { assert_eq!(a[[3, 0]], 2.0); assert_eq!(a[[3, 1]], 3.0); assert_eq!(a[[2, 0]], 4.0); - // assert_eq!(a[[2, 1]], 0.0); assert_eq!(a[[1, 0]], 5.0); assert_eq!(a.entries, vec![1.0, 0.0, 5.0, 0.0, 4.0, 0.0, 2.0, 3.0]); assert_eq!(a.start_indices, vec![0]); assert_eq!(a.lengths, vec![8]); - - println!( - "entries: {:?}, start_indices: {:?}, lengths: {:?}", - a.entries, a.start_indices, a.lengths - ); } #[test] @@ -373,15 +337,7 @@ mod tests { a[[2, 4]] = 2; a[[4, 1]] = 3; a[[4, 4]] = 4; - println!( - "entries: {:?}, start_indices: {:?}, lengths: {:?}", - a.entries, a.start_indices, a.lengths - ); a[[5, 0]] = 5; - println!( - "entries: {:?}, start_indices: {:?}, lengths: {:?}", - a.entries, a.start_indices, a.lengths - ); assert_eq!( a.indexed_iter().collect::>(), &[ @@ -699,19 +655,13 @@ mod tests { assert_eq!(array[[3, 4, 3]], 1.0); assert_eq!(array[[3, 4, 4]], 2.0); - // assert_eq!(array[[3, 4, 5]], 0.0); + assert_eq!(array[[3, 4, 5]], 0.0); assert_eq!(array[[3, 4, 6]], 3.0); assert_eq!(array[[3, 5, 1]], 4.0); - // assert_eq!(array[[3, 5, 2]], 0.0); - // assert_eq!(array[[3, 5, 3]], 0.0); - // assert_eq!(array[[3, 5, 4]], 0.0); - // assert_eq!(array[[3, 5, 5]], 0.0); - // assert_eq!(array[[3, 5, 6]], 0.0); assert_eq!(array[[3, 5, 7]], 5.0); assert_eq!(array[[4, 3, 9]], 6.0); - // assert_eq!(array.len(), 6); - // assert_eq!(array.zeros(), 6); + assert_eq!(array.explicit_zeros(), 1); } } From ced69fef39602cf222ac619598bc3b2a03db3317 Mon Sep 17 00:00:00 2001 From: janw20 Date: Thu, 18 Apr 2024 17:42:51 +0200 Subject: [PATCH 06/14] Add small changes to the documentation --- pineappl/src/packed_array.rs | 22 +++++++++------------- 1 file changed, 9 insertions(+), 13 deletions(-) diff --git a/pineappl/src/packed_array.rs b/pineappl/src/packed_array.rs index 6cb1ca2ac..cf4dbedda 100644 --- a/pineappl/src/packed_array.rs +++ b/pineappl/src/packed_array.rs @@ -4,7 +4,7 @@ use std::{mem, vec}; use ndarray::ArrayView3; use serde::{Deserialize, Serialize}; -/// `D`-dimensional array similar to ndarray::ArrayBase, except that `T::default()` is not stored to save space. Instead, adjacent elements are grouped together and only the index of their first element (`start_index`) and the length of the group (`lengths`) is stored. +/// `D`-dimensional array similar to ndarray::ArrayBase, except that `T::default()` is not stored to save space. Instead, adjacent non-default elements are grouped together and the index of their first element (`start_index`) and the length of the group (`lengths`) is stored. #[derive(Clone, Deserialize, Serialize)] pub struct PackedArray { /// The actual values stored in the array. The length of `entries` is always the sum of the elements in `lengths`. @@ -56,7 +56,7 @@ impl PackedArray { / mem::size_of::() } - /// Returns the number of default (zero) elements that are explicitly stored in `entries`. If there is one default elements between adjacent groups, it is more economical to store the one default element explicitly and merge the two groups, than to store the `start_indices` and `lengths` of those groups. + /// Returns the number of default (zero) elements that are explicitly stored in `entries`. If there is one default element between adjacent groups, it is more economical to store the one default element explicitly and merge the two groups, than to store the `start_indices` and `lengths` of both groups. #[must_use] pub fn explicit_zeros(&self) -> usize { self.entries.iter().filter(|x| **x == T::default()).count() @@ -68,7 +68,7 @@ impl PackedArray { self.entries.iter().filter(|x| **x != T::default()).count() } - /// Returns an `Iterator` over the non-default (non-zero) elements of this array. The type of an iterator element is `([usize; D], T)`. + /// Returns an `Iterator` over the non-default (non-zero) elements of this array. The type of an iterator element is `([usize; D], T)` where the first element of the tuple is the index and the second element is the value. pub fn indexed_iter(&self) -> impl Iterator + '_ { self.start_indices .iter() @@ -83,7 +83,6 @@ impl PackedArray { } impl, const D: usize> MulAssign for PackedArray { - /// Perform `self *= rhs` as elementwise multiplication (in place). fn mul_assign(&mut self, rhs: T) { self.entries.iter_mut().for_each(|x| *x *= rhs); @@ -91,7 +90,6 @@ impl, const D: usize> MulAssign for PackedArray } impl PackedArray { - /// Converts `array` into a `PackedArray` of dimension 3. #[must_use] pub fn from_ndarray(array: ArrayView3, xstart: usize, xsize: usize) -> Self { @@ -110,15 +108,15 @@ impl PackedArray { } } -/// Converts a `multi_index` into an flat index. +/// Converts a `multi_index` into a flat index. fn ravel_multi_index(multi_index: &[usize; D], dimensions: &[usize]) -> usize { assert_eq!(multi_index.len(), dimensions.len()); - + multi_index - .iter() - .skip(1) - .zip(dimensions.iter().skip(1)) - .fold(multi_index[0], |acc, (i, d)| acc * d + i) + .iter() + .skip(1) + .zip(dimensions.iter().skip(1)) + .fold(multi_index[0], |acc, (i, d)| acc * d + i) } /// Converts a flat `index` into a multi_index. @@ -248,7 +246,6 @@ impl IndexMut<[usize; D]> } } - #[cfg(test)] mod tests { use ndarray::Array3; @@ -663,5 +660,4 @@ mod tests { assert_eq!(array.explicit_zeros(), 1); } - } From 3ca681fcce2def5ef378e7a41dc7c5aadc4ec6d4 Mon Sep 17 00:00:00 2001 From: janw20 Date: Thu, 18 Apr 2024 17:59:02 +0200 Subject: [PATCH 07/14] Fix ordering --- pineappl/src/lib.rs | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pineappl/src/lib.rs b/pineappl/src/lib.rs index 6c1479b13..7a115d085 100644 --- a/pineappl/src/lib.rs +++ b/pineappl/src/lib.rs @@ -8,10 +8,10 @@ pub mod evolution; pub mod fk_table; pub mod grid; pub mod import_only_subgrid; -pub mod packed_array; pub mod lagrange_subgrid; pub mod lumi; pub mod ntuple_subgrid; +pub mod packed_array; pub mod pids; pub mod sparse_array3; pub mod subgrid; From 20580028290f366bb90351ece8a4fc44e6d2c839 Mon Sep 17 00:00:00 2001 From: Christopher Schwan Date: Tue, 23 Apr 2024 13:40:04 +0200 Subject: [PATCH 08/14] Adjust documentation a bit --- pineappl/src/packed_array.rs | 27 ++++++++++++++++++--------- 1 file changed, 18 insertions(+), 9 deletions(-) diff --git a/pineappl/src/packed_array.rs b/pineappl/src/packed_array.rs index cf4dbedda..e15a67615 100644 --- a/pineappl/src/packed_array.rs +++ b/pineappl/src/packed_array.rs @@ -1,13 +1,18 @@ +//! Provides the [`PackedArray`] struct. + use std::ops::{Index, IndexMut, MulAssign}; use std::{mem, vec}; use ndarray::ArrayView3; use serde::{Deserialize, Serialize}; -/// `D`-dimensional array similar to ndarray::ArrayBase, except that `T::default()` is not stored to save space. Instead, adjacent non-default elements are grouped together and the index of their first element (`start_index`) and the length of the group (`lengths`) is stored. +/// `D`-dimensional array similar to [`ndarray::ArrayBase`], except that `T::default()` is not +/// stored to save space. Instead, adjacent non-default elements are grouped together and the index +/// of their first element (`start_index`) and the length of the group (`lengths`) is stored. #[derive(Clone, Deserialize, Serialize)] pub struct PackedArray { - /// The actual values stored in the array. The length of `entries` is always the sum of the elements in `lengths`. + /// The actual values stored in the array. The length of `entries` is always the sum of the + /// elements in `lengths`. entries: Vec, /// The indices of the first elements in each group. start_indices: Vec, @@ -22,7 +27,6 @@ impl PackedArray { #[must_use] pub fn new(shape: [usize; D]) -> Self { Self { - // entries: BTreeMap::new(), entries: vec![], start_indices: vec![], lengths: vec![], @@ -49,14 +53,18 @@ impl PackedArray { self.lengths.clear(); } - /// Returns the overhead of storing the `start_indices` and the `lengths` of the groups, in units of `f64`. + /// Returns the overhead of storing the `start_indices` and the `lengths` of the groups, in + /// units of `f64`. #[must_use] pub fn overhead(&self) -> usize { ((self.start_indices.len() + self.lengths.len()) * mem::size_of::()) / mem::size_of::() } - /// Returns the number of default (zero) elements that are explicitly stored in `entries`. If there is one default element between adjacent groups, it is more economical to store the one default element explicitly and merge the two groups, than to store the `start_indices` and `lengths` of both groups. + /// Returns the number of default (zero) elements that are explicitly stored in `entries`. If + /// there is one default element between adjacent groups, it is more economical to store the + /// one default element explicitly and merge the two groups, than to store the `start_indices` + /// and `lengths` of both groups. #[must_use] pub fn explicit_zeros(&self) -> usize { self.entries.iter().filter(|x| **x == T::default()).count() @@ -68,7 +76,9 @@ impl PackedArray { self.entries.iter().filter(|x| **x != T::default()).count() } - /// Returns an `Iterator` over the non-default (non-zero) elements of this array. The type of an iterator element is `([usize; D], T)` where the first element of the tuple is the index and the second element is the value. + /// Returns an `Iterator` over the non-default (non-zero) elements of this array. The type of + /// an iterator element is `([usize; D], T)` where the first element of the tuple is the index + /// and the second element is the value. pub fn indexed_iter(&self) -> impl Iterator + '_ { self.start_indices .iter() @@ -83,14 +93,13 @@ impl PackedArray { } impl, const D: usize> MulAssign for PackedArray { - /// Perform `self *= rhs` as elementwise multiplication (in place). fn mul_assign(&mut self, rhs: T) { self.entries.iter_mut().for_each(|x| *x *= rhs); } } impl PackedArray { - /// Converts `array` into a `PackedArray` of dimension 3. + /// Converts `array` into a `PackedArray`. #[must_use] pub fn from_ndarray(array: ArrayView3, xstart: usize, xsize: usize) -> Self { let shape = array.shape(); @@ -119,7 +128,7 @@ fn ravel_multi_index(multi_index: &[usize; D], dimensions: &[usi .fold(multi_index[0], |acc, (i, d)| acc * d + i) } -/// Converts a flat `index` into a multi_index. +/// Converts a flat `index` into a `multi_index`. fn unravel_index(index: usize, dimensions: &[usize]) -> [usize; D] { assert!(index < dimensions.iter().product()); let mut indices = [0; D]; From 1f74315b938478728f226110b8bc71dbb9e98006 Mon Sep 17 00:00:00 2001 From: Christopher Schwan Date: Tue, 23 Apr 2024 13:41:24 +0200 Subject: [PATCH 09/14] Group `use` statements together --- pineappl/src/packed_array.rs | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/pineappl/src/packed_array.rs b/pineappl/src/packed_array.rs index e15a67615..d1d4ed680 100644 --- a/pineappl/src/packed_array.rs +++ b/pineappl/src/packed_array.rs @@ -1,10 +1,9 @@ //! Provides the [`PackedArray`] struct. -use std::ops::{Index, IndexMut, MulAssign}; -use std::{mem, vec}; - use ndarray::ArrayView3; use serde::{Deserialize, Serialize}; +use std::mem; +use std::ops::{Index, IndexMut, MulAssign}; /// `D`-dimensional array similar to [`ndarray::ArrayBase`], except that `T::default()` is not /// stored to save space. Instead, adjacent non-default elements are grouped together and the index From 68f66ff0572bb882a676faa55bf2fefd9628ea47 Mon Sep 17 00:00:00 2001 From: Christopher Schwan Date: Tue, 23 Apr 2024 14:52:34 +0200 Subject: [PATCH 10/14] Use `zip` early In this way all necessary modifications are made uniformly (with late `zip` it's possible to forget the first/second `skip`/`rev` --- pineappl/src/packed_array.rs | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pineappl/src/packed_array.rs b/pineappl/src/packed_array.rs index d1d4ed680..437b68c20 100644 --- a/pineappl/src/packed_array.rs +++ b/pineappl/src/packed_array.rs @@ -122,8 +122,8 @@ fn ravel_multi_index(multi_index: &[usize; D], dimensions: &[usi multi_index .iter() + .zip(dimensions) .skip(1) - .zip(dimensions.iter().skip(1)) .fold(multi_index[0], |acc, (i, d)| acc * d + i) } @@ -133,8 +133,8 @@ fn unravel_index(index: usize, dimensions: &[usize]) -> [usize; let mut indices = [0; D]; indices .iter_mut() + .zip(dimensions) .rev() - .zip(dimensions.iter().rev()) .fold(index, |acc, (i, d)| { *i = acc % d; acc / d From 4227d48ac07ea85254d4f8f6b59f42ddb6a1aed4 Mon Sep 17 00:00:00 2001 From: Christopher Schwan Date: Wed, 24 Apr 2024 10:24:37 +0200 Subject: [PATCH 11/14] Add some cosmetic changes --- pineappl/src/packed_array.rs | 18 ++++++++---------- 1 file changed, 8 insertions(+), 10 deletions(-) diff --git a/pineappl/src/packed_array.rs b/pineappl/src/packed_array.rs index 437b68c20..a83eb3e1c 100644 --- a/pineappl/src/packed_array.rs +++ b/pineappl/src/packed_array.rs @@ -117,23 +117,23 @@ impl PackedArray { } /// Converts a `multi_index` into a flat index. -fn ravel_multi_index(multi_index: &[usize; D], dimensions: &[usize]) -> usize { - assert_eq!(multi_index.len(), dimensions.len()); +fn ravel_multi_index(multi_index: &[usize; D], shape: &[usize]) -> usize { + assert_eq!(multi_index.len(), shape.len()); multi_index .iter() - .zip(dimensions) + .zip(shape) .skip(1) .fold(multi_index[0], |acc, (i, d)| acc * d + i) } /// Converts a flat `index` into a `multi_index`. -fn unravel_index(index: usize, dimensions: &[usize]) -> [usize; D] { - assert!(index < dimensions.iter().product()); +fn unravel_index(index: usize, shape: &[usize]) -> [usize; D] { + assert!(index < shape.iter().product()); let mut indices = [0; D]; indices .iter_mut() - .zip(dimensions) + .zip(shape) .rev() .fold(index, |acc, (i, d)| { *i = acc % d; @@ -210,13 +210,12 @@ impl IndexMut<[usize; D]> let distance = raveled_index - (start_index + length) + 1; self.lengths[point - 1] += distance; - for _ in 0..(distance) { + for _ in 0..distance { self.entries.insert(point_entries, Default::default()); } if let Some(start_index_next) = self.start_indices.get(point) { if raveled_index + threshold_distance >= *start_index_next { - // println!("test 4"); let distance_next = start_index_next - raveled_index; self.lengths[point - 1] += distance_next - 1 + self.lengths[point]; @@ -256,11 +255,10 @@ impl IndexMut<[usize; D]> #[cfg(test)] mod tests { + use super::*; use ndarray::Array3; use std::mem; - use super::*; - #[test] fn unravel_index() { assert_eq!(super::unravel_index(0, &[3, 2]), [0, 0]); From cac89150e10ca4075ff0a9e451049203b65ea277 Mon Sep 17 00:00:00 2001 From: Christopher Schwan Date: Thu, 25 Apr 2024 08:45:32 +0200 Subject: [PATCH 12/14] Simplify `ravel_multi_index` and `unravel_index` --- pineappl/src/packed_array.rs | 17 ++++++----------- 1 file changed, 6 insertions(+), 11 deletions(-) diff --git a/pineappl/src/packed_array.rs b/pineappl/src/packed_array.rs index a83eb3e1c..7458475f6 100644 --- a/pineappl/src/packed_array.rs +++ b/pineappl/src/packed_array.rs @@ -123,22 +123,17 @@ fn ravel_multi_index(multi_index: &[usize; D], shape: &[usize]) multi_index .iter() .zip(shape) - .skip(1) - .fold(multi_index[0], |acc, (i, d)| acc * d + i) + .fold(0, |acc, (i, d)| acc * d + i) } /// Converts a flat `index` into a `multi_index`. -fn unravel_index(index: usize, shape: &[usize]) -> [usize; D] { +fn unravel_index(mut index: usize, shape: &[usize]) -> [usize; D] { assert!(index < shape.iter().product()); let mut indices = [0; D]; - indices - .iter_mut() - .zip(shape) - .rev() - .fold(index, |acc, (i, d)| { - *i = acc % d; - acc / d - }); + for (i, d) in indices.iter_mut().zip(shape).rev() { + *i = index % d; + index /= d; + } indices } From ca8f09c29d04baa18c9951f1a354db7013a70c22 Mon Sep 17 00:00:00 2001 From: Christopher Schwan Date: Thu, 25 Apr 2024 08:46:17 +0200 Subject: [PATCH 13/14] Use `assert_eq` instead of `debug_assert_eq` --- pineappl/src/packed_array.rs | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pineappl/src/packed_array.rs b/pineappl/src/packed_array.rs index 7458475f6..bbb3ec3a9 100644 --- a/pineappl/src/packed_array.rs +++ b/pineappl/src/packed_array.rs @@ -141,7 +141,7 @@ impl Index<[usize; D]> for Packed type Output = T; fn index(&self, index: [usize; D]) -> &Self::Output { - debug_assert_eq!(index.len(), self.shape.len()); + assert_eq!(index.len(), self.shape.len()); assert!( index.iter().zip(self.shape.iter()).all(|(&i, &d)| i < d), "index {:?} is out of bounds for array of shape {:?}", @@ -176,7 +176,7 @@ impl IndexMut<[usize; D]> for PackedArray { fn index_mut(&mut self, index: [usize; D]) -> &mut Self::Output { - debug_assert_eq!(index.len(), self.shape.len()); + assert_eq!(index.len(), self.shape.len()); assert!( index.iter().zip(self.shape.iter()).all(|(&i, &d)| i < d), "index {:?} is out of bounds for array of shape {:?}", From 270aeddd5c4f432e6e030882c2ad47408e9ebeda Mon Sep 17 00:00:00 2001 From: Christopher Schwan Date: Thu, 25 Apr 2024 08:53:01 +0200 Subject: [PATCH 14/14] Replace loop over `insert` with `splice` This should be more efficient since Rust can calculate the final size and thus allocate only once --- pineappl/src/packed_array.rs | 23 +++++++++++++---------- 1 file changed, 13 insertions(+), 10 deletions(-) diff --git a/pineappl/src/packed_array.rs b/pineappl/src/packed_array.rs index bbb3ec3a9..80109c4ab 100644 --- a/pineappl/src/packed_array.rs +++ b/pineappl/src/packed_array.rs @@ -2,6 +2,7 @@ use ndarray::ArrayView3; use serde::{Deserialize, Serialize}; +use std::iter; use std::mem; use std::ops::{Index, IndexMut, MulAssign}; @@ -204,10 +205,10 @@ impl IndexMut<[usize; D]> } else if raveled_index < start_index + length + threshold_distance { let distance = raveled_index - (start_index + length) + 1; self.lengths[point - 1] += distance; - - for _ in 0..distance { - self.entries.insert(point_entries, Default::default()); - } + self.entries.splice( + point_entries..point_entries, + iter::repeat(Default::default()).take(distance), + ); if let Some(start_index_next) = self.start_indices.get(point) { if raveled_index + threshold_distance >= *start_index_next { @@ -216,9 +217,10 @@ impl IndexMut<[usize; D]> self.lengths[point - 1] += distance_next - 1 + self.lengths[point]; self.lengths.remove(point); self.start_indices.remove(point); - for _ in 0..(distance_next - 1) { - self.entries.insert(point_entries, Default::default()); - } + self.entries.splice( + point_entries..point_entries, + iter::repeat(Default::default()).take(distance_next - 1), + ); } } @@ -232,9 +234,10 @@ impl IndexMut<[usize; D]> self.start_indices[point] = raveled_index; self.lengths[point] += distance; - for _ in 0..distance { - self.entries.insert(point_entries, Default::default()); - } + self.entries.splice( + point_entries..point_entries, + iter::repeat(Default::default()).take(distance), + ); return &mut self.entries[point_entries]; } }