q1tsim/
permutation.rs

1// Copyright 2019 Q1t BV
2//
3// Licensed under the Apache License, Version 2.0 (the "License");
4// you may not use this file except in compliance with the License.
5// You may obtain a copy of the License at
6//
7//    http://www.apache.org/licenses/LICENSE-2.0
8//
9// Unless required by applicable law or agreed to in writing, software
10// distributed under the License is distributed on an "AS IS" BASIS,
11// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
12// See the License for the specific language governing permissions and
13// limitations under the License.
14
15
16extern crate ndarray;
17extern crate num_traits;
18
19/// Struct for permutations
20///
21/// Struct Permutation is used to represents permutations. It can be used to
22/// shuffle the elements in a vector, or rows and columns in a matrix.
23pub struct Permutation
24{
25    /// The permuted indexes.
26    idxs: Vec<usize>
27}
28
29impl Permutation
30{
31    /// Create a new permutation.
32    ///
33    /// Create a new permutation from `idxs`. If `n` is the length of `idxs`,
34    /// the array should contain all elements in the range `[0, n-1]` (inclusive).
35    /// When `idxs[i] = j`, application of the permutation on a vector means that
36    /// the the `j`th element of the original vector is taken to position `i` in
37    /// the permuted vector.
38    pub fn new(idxs: Vec<usize>) -> crate::error::Result<Self>
39    {
40        if idxs.is_empty()
41        {
42            return Err(crate::error::Error::EmptyPermutation);
43        }
44
45        let n = idxs.len();
46        let max = *idxs.iter().max().unwrap();
47        if max >= n
48        {
49            return Err(crate::error::Error::InvalidPermutationElement(max, n));
50        }
51
52        let mut seen = vec![false; n];
53        for &pi in idxs.iter()
54        {
55            if seen[pi]
56            {
57                return Err(crate::error::Error::DoublePermutationElement(pi));
58            }
59            seen[pi] = true;
60        }
61
62        Ok(Permutation { idxs: idxs })
63    }
64
65    /// Return the permutation index vector
66    pub fn indices(&self) -> &[usize]
67    {
68        &self.idxs
69    }
70
71    /// Return the size of the permutation
72    pub fn size(&self) -> usize
73    {
74        self.idxs.len()
75    }
76
77    /// The inverse permutation
78    ///
79    /// Compute the inverse of this permutation. If this permutation maps
80    /// element `i` to `j`, then then in the inverse `j` is mapped to `i`.
81    pub fn inverse(&self) -> Self
82    {
83        let mut invidxs = vec![0; self.size()];
84        for (i, &pi) in self.idxs.iter().enumerate()
85        {
86            invidxs[pi] = i;
87        }
88        Permutation::new(invidxs).unwrap()
89    }
90
91    /// Matrix representation.
92    ///
93    /// Return a matrix representation of the permutation. For a matrix
94    /// representation `P` of permutation `perm`, the result of `P.dot(v)` for a
95    /// vector `v` is the same as `perm.apply_vec(&mut v)`.
96    pub fn matrix<T>(&self) -> ndarray::Array2<T>
97    where T: Clone + num_traits::Zero + num_traits::One
98    {
99        let n = self.size();
100        let mut res = ndarray::Array::zeros((n, n));
101        for (i, &pi) in self.idxs.iter().enumerate()
102        {
103            res[[i, pi]] = T::one();
104        }
105        res
106    }
107
108    fn apply<F>(&self, mut swap: F)
109    where F: FnMut(usize, usize)
110    {
111        let mut in_place = vec![false; self.size()];
112        for i in 0..self.size()
113        {
114            if !in_place[i]
115            {
116                let mut j = i;
117                let mut pj = self.idxs[j];
118                while pj != i
119                {
120                    swap(j, pj);
121                    in_place[j] = true;
122                    j = pj;
123                    pj = self.idxs[pj];
124                }
125                in_place[j] = true;
126            }
127        }
128    }
129
130    /// Permute a vector
131    ///
132    /// Apply this permutation to vector `v`.
133    pub fn apply_vec_in_place<T>(&self, v: &mut ndarray::Array1<T>)
134    {
135        self.apply(|i, j| v.swap(i, j));
136    }
137
138    /// Permute a vector
139    ///
140    /// Apply this permutation to vector `v`, and store the result in `perm_v`.
141    pub fn apply_vec_into<T>(&self, v: ndarray::ArrayView1<T>,
142        mut perm_v: ndarray::ArrayViewMut1<T>)
143    where T: Copy
144    {
145        for (new_idx, &old_idx) in self.idxs.iter().enumerate()
146        {
147            perm_v[new_idx] = v[old_idx];
148        }
149    }
150
151    /// Permute a vector
152    ///
153    /// Apply the inverse of this permutation to vector `v`, and store the
154    /// result in `inv_perm_v`.
155    pub fn apply_inverse_vec_into<T>(&self, v: ndarray::ArrayView1<T>,
156        mut inv_perm_v: ndarray::ArrayViewMut1<T>)
157    where T: Copy
158    {
159        for (new_idx, &old_idx) in self.idxs.iter().enumerate()
160        {
161            inv_perm_v[old_idx] = v[new_idx];
162        }
163    }
164
165    /// Transform a matrix
166    ///
167    /// Transform matrix `a`, permuting both its rows and columns by this
168    /// permutation.
169    pub fn transform<T>(&self, a: &ndarray::Array2<T>) -> ndarray::Array2<T>
170    where T: Copy
171    {
172        a.select(ndarray::Axis(0), &self.idxs).select(ndarray::Axis(1), &self.idxs)
173    }
174}
175
176#[cfg(test)]
177mod tests
178{
179    use super::Permutation;
180
181    #[test]
182    fn test_new()
183    {
184        let perm = Permutation::new(vec![0, 1, 2, 3]).unwrap();
185        assert_eq!(perm.indices(), &[0, 1, 2, 3]);
186        let perm = Permutation::new(vec![3, 1, 0, 2]).unwrap();
187        assert_eq!(perm.indices(), &[3, 1, 0, 2]);
188
189        assert!(matches!(Permutation::new(vec![]), Err(crate::error::Error::EmptyPermutation)));
190        assert!(matches!(Permutation::new(vec![0, 2, 3, 4]), Err(crate::error::Error::InvalidPermutationElement(4, 4))));
191        assert!(matches!(Permutation::new(vec![0, 2, 3, 2]), Err(crate::error::Error::DoublePermutationElement(2))));
192    }
193
194    #[test]
195    fn test_inverse()
196    {
197        let iperm = Permutation::new(vec![0, 1, 2, 3]).unwrap().inverse();
198        assert_eq!(iperm.indices(), &[0, 1, 2, 3]);
199        let iperm = Permutation::new(vec![3, 0, 1, 2]).unwrap().inverse();
200        assert_eq!(iperm.indices(), &[1, 2, 3, 0]);
201        let iperm = Permutation::new(vec![6, 2, 7, 1, 4, 3, 0, 5]).unwrap().inverse();
202        assert_eq!(iperm.indices(), &[6, 3, 1, 5, 4, 7, 0, 2]);
203    }
204
205    #[test]
206    fn test_matrix()
207    {
208        let perm = Permutation::new(vec![1, 3, 0, 2]).unwrap();
209        assert_eq!(perm.matrix::<usize>(), array![
210            [0, 1, 0, 0],
211            [0, 0, 0, 1],
212            [1, 0, 0, 0],
213            [0, 0, 1, 0]
214        ]);
215        let perm = Permutation::new(vec![1, 3, 5, 2, 7, 0, 4, 6]).unwrap();
216        assert_eq!(perm.matrix::<f64>(), array![
217            [0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
218            [0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0],
219            [0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0],
220            [0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0],
221            [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0],
222            [1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
223            [0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0],
224            [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0]
225        ]);
226
227        let perm = Permutation::new(vec![1, 3, 0, 2]).unwrap();
228        let v = array![1, 3, 5, 7];
229        assert_eq!(perm.matrix().dot(&v), array![3, 7, 1, 5]);
230
231        let perm = Permutation::new(vec![1, 3, 5, 2, 7, 0, 4, 6]).unwrap();
232        let v = array![1.2, 2.3, 3.4, 4.5, 5.6, 6.7, 7.8, 8.9];
233        assert_eq!(perm.matrix().dot(&v), array![2.3, 4.5, 6.7, 3.4, 8.9, 1.2, 5.6, 7.8]);
234    }
235
236    #[test]
237    fn test_apply_vec_in_place()
238    {
239        let perm = Permutation::new(vec![1, 3, 0, 2]).unwrap();
240        let mut v = array![1, 3, 5, 7];
241        perm.apply_vec_in_place(&mut v);
242        assert_eq!(v, array![3, 7, 1, 5]);
243
244        let perm = Permutation::new(vec![1, 3, 5, 2, 7, 0, 4, 6]).unwrap();
245        let mut v = array![1.2, 2.3, 3.4, 4.5, 5.6, 6.7, 7.8, 8.9];
246        perm.apply_vec_in_place(&mut v);
247        assert_eq!(v, array![2.3, 4.5, 6.7, 3.4, 8.9, 1.2, 5.6, 7.8]);
248    }
249
250    #[test]
251    fn test_apply_vec()
252    {
253        let perm = Permutation::new(vec![1, 3, 0, 2]).unwrap();
254        let v = array![1, 3, 5, 7];
255        let mut work = array![0, 0, 0, 0];
256        perm.apply_vec_into(v.view(), work.view_mut());
257        assert_eq!(work, array![3, 7, 1, 5]);
258
259        let perm = Permutation::new(vec![1, 3, 5, 2, 7, 0, 4, 6]).unwrap();
260        let v = array![1.2, 2.3, 3.4, 4.5, 5.6, 6.7, 7.8, 8.9];
261        let mut work = array![0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0];
262        perm.apply_vec_into(v.view(), work.view_mut());
263        assert_eq!(work, array![2.3, 4.5, 6.7, 3.4, 8.9, 1.2, 5.6, 7.8]);
264    }
265
266    #[test]
267    fn test_transform()
268    {
269        let a = array![[1, 2, 3], [4, 5, 6], [7, 8, 9]];
270        let perm = Permutation::new(vec![2, 0, 1]).unwrap();
271        let b = perm.transform(&a);
272        assert_eq!(b, array![[9, 7, 8], [3, 1, 2], [6, 4, 5]]);
273    }
274}
q1tsim/permutation.rs

q1tsim/
permutation.rs
