hilbert 0.1.2

Hilbert curve transform and inverse for points having two to thousands of dimensions, using Skilling's algorithm.
Documentation
//! The `permutation` module contains the `Permutation` struct, which is used to 
//! reorder the coordinates in one or more `Points` in a consistent way.
//! This permits the generation of multiple **Hilbert Curves**, each rendering
//! the points in a different order. With multiple curves, you can improve the 
//! accuracy of K-Nearest Nighbor searches and clustering algorithms.
use rand::thread_rng;
use rand::seq::SliceRandom;

/// Defines a one-to-one mapping from one ordering of values to another.
/// This can be used to rearrange the coordinates in a Point.
/// 
/// No coordinate value is lost or duplicated, so long as every value from 0 to dimensions - 1 is present in the `moves_to` collection.
#[derive(Clone, PartialEq, Debug)]
pub struct Permutation {
    /// Indicates to which dimension a give dimension should be moved. 
    /// 
    /// Example: 
    /// If permutation [1,2,0,3] is applied to Point [6,7,8,9] the result is Point [8,6,7,9].
    moves_to : Vec<usize>,

    /// Inverted index formed from `moves_to`. 
    moves_from : Vec<usize>
}

impl Permutation {
    /// Create a Permutation that causes no changes to the order of dimensions if applied to a Point.
    pub fn identity(dimensions : usize) -> Self {
        let mut p = Permutation { 
            moves_to : Vec::with_capacity(dimensions),
            moves_from : Vec::with_capacity(dimensions),
        };
        for dimension in 0 .. dimensions {
            p.moves_to[dimension] = dimension;
            p.moves_from[dimension] = dimension;
        }
        p
    }

    /// Create a Permutation that performs the given set of moves. 
    pub fn new(moves : &[usize]) -> Self {
        let mut p = Permutation { 
            moves_to : moves.to_vec(), 
            moves_from : vec![0; moves.len()]
        };
        if !p.is_valid() { panic!("Permutation has missing, duplicate and/or out of range coordinate indices"); }
        p.invert();
        p
    }

    /// Create a random permutation which will cause almost every dimension to be shifted to a new position. 
    pub fn random(dimensions : usize) -> Self {
        let mut p = Self::identity(dimensions);
        let mut rng = thread_rng();
        p.moves_to.shuffle(&mut rng);
        p.invert();
        p
    }

    /// Create a new Permutation where the two given dimensions are swapped.
    /// 
    /// This is useful for controlled permutation, unlike the chaotically scrambled 
    /// `Permutations` generated by the `random` method. 
    pub fn swap(&self, first_dimension : usize, second_dimension : usize) -> Self {
        let mut p = self.clone();
        let i_inverse = self.moves_to[first_dimension];
        let j_inverse = self.moves_to[second_dimension];
        p.moves_to.swap(first_dimension, second_dimension);
        p.moves_from.swap(i_inverse, j_inverse);
        p
    }

    /// Number of dimensions of points that can be handled by the permutation.
    pub fn dimensions(&self) -> usize {
        self.moves_to.len()
    }

    /// Verify that all the numbers from zero to `dimensions() - 1` are present in moves. 
    /// If there are missing, duplicate, or out of range values, the Permutation is invalid. 
    /// 
    /// These checks are performed, where `N = dimensions() - 1`: 
    /// 
    ///   1. Sum of values compared to N(N+1)/2
    ///   2. Sum of squares compared to N(N+1)(N+2)/6
    ///   3. Cumulative XOR
    ///   4. Range test 0 <= x <= N
    fn is_valid(&self) -> bool {
        // I do not know if this algorithm is guaranteed to spot all errors, but it is highly unlikely to miss any.
        let n : u64 = (self.dimensions() - 1) as u64;
        let mut sum : u64 = 0;
        let mut squares : u64 = 0;
        let mut expected_xor : usize = 0;
        let mut xor : usize = 0;
        for (i, dim) in self.moves_to.iter().enumerate() {
            if dim.clone() as u64 > n { return false; }
            xor ^= dim;
            expected_xor ^= i;
            let dim64 : u64 = dim.clone() as u64;
            sum += dim64;
            squares += dim64 * dim64;
        }
        (xor == expected_xor)
          && (sum == n * (n + 1) / 2)
          && (squares == n * (n + 1) * (2*n + 1) / 6)
    }

    /// Initialize the `moves_from` collection with an inverse of the `moves_to` collection.
    ///  
    /// This allows the apply method to fill the output vector in ascending order using push(),
    /// so that we do not need to fill the vector with default values first. 
    /// (In C#, would have initialized `moves_from` with nulls, but do not want to use Options here; too tedious.)
    fn invert(&mut self) {
        for i in 0 .. self.dimensions() {
            self.moves_from[self.moves_to[i]] = i;
        }
    }

    /// Apply the `Permutation` to data, cloning and reordering the input values appropriately.
    pub fn apply<T>(&self, unpermuted_data : &[T]) -> Vec<T>
    where T : Clone
    {
        let mut permuted_values = Vec::with_capacity(self.dimensions());
        for dim in 0 .. self.dimensions() {
            permuted_values.push(unpermuted_data[self.moves_from[dim]].clone());
        }
        permuted_values
    }
}

#[cfg(test)]
/// Tests of the Permutation methods.
mod tests {
    #[allow(unused_imports)]
    use spectral::prelude::*;
    use crate::permutation::Permutation;

    #[test]
    fn apply() {
        let data : [i32; 5] = [2,4,6,8,10];
        let moves : [usize; 5] = [2,0,1,4,3];
        let permutation = Permutation::new(&moves);
        let actual_permuted_data = permutation.apply(&data);
        let expected_permuted_data = vec![4,6,2,10,8];
        asserting("Correct permutation").that(&actual_permuted_data).is_equal_to(expected_permuted_data);
    }

    #[test]
    #[should_panic]
    fn is_valid_fails() {
        let data : [i32; 5] = [2,4,6,8,10];
        let moves : [usize; 5] = [2,0,1,4,999];
        // This line should panic with: 'Permutation has missing, duplicate and/or out of range coordinate indices'
        let permutation = Permutation::new(&moves);
        let _actual_permuted_data = permutation.apply(&data);
    }

    #[test]
    fn swap() {
        let data : [i32; 5] = [2,4,6,8,10];
        let permutation = Permutation::new(&[2,0,1,4,3]);
        let swapped = permutation.swap(2, 3); // -> [2,0,4,1,3]
        let actual_permuted_data = swapped.apply(&data);
        let expected_permuted_data = vec![4,8,2,10,6];
        asserting("swapped permutation").that(&actual_permuted_data).is_equal_to(expected_permuted_data);
    }

}