steiner-tree 0.0.1

Fast construction of rectilinear steiner minimal trees (RSMT) in two dimensions.
Documentation
// SPDX-FileCopyrightText: 2022 Thomas Kramer <code@tkramer.ch>
//
// SPDX-License-Identifier: GPL-3.0-or-later

//! Helper functions related to permuations and rank of permuations (group index).

use itertools::Itertools;

use super::position_sequence::PositionSequence;
use std::cell::RefCell;

// fn permutation_rank(permutation: &[usize]) -> usize {
//     let n = permutation.len();
//
//     // Check that all numbers from 0 to n-1 are present in the permutation.
//     debug_assert!(
//         (0..n).all(|i| permutation.contains(&i)),
//         "permutation must consist of subsequent numbers"
//     );
//
//     if n <= 1 {
//         0
//     } else {
//         let mut rank = 0;
//         let f: usize = (0..n - 1).product();
//
//         todo!()
//     }
// }

/// Iterate over all permutations of the position sequence with length `num_pins`.
/// Length of the iterator is, of course, `factorial(num_pins)`.
pub fn all_position_sequences(num_pins: usize) -> impl Iterator<Item=PositionSequence> {
    let base: Vec<_> = (0..num_pins).collect();
    base.into_iter()
        .permutations(num_pins)
        .map(|s| PositionSequence::new(s))
}

pub fn group_index(position_sequence: &PositionSequence) -> usize {
    // Use a thread local scratch buffer to avoid allocations.
    thread_local! {
            static SCRATCH: RefCell<Vec<usize>> = RefCell::new(vec![]);
        }

    SCRATCH.with(|scratch| {
        let mut scratch = scratch.borrow_mut();
        scratch.clear();
        scratch.extend(0..position_sequence.len());
        modified_position_sequence(&mut scratch, position_sequence.sequence());
        group_index_from_modified_position_sequence(&scratch)
    })
}


/// Compute the index of a group. This is basically the rank of the permutation.
fn group_index_from_modified_position_sequence(modified_position_sequence: &[usize]) -> usize {
    let n = modified_position_sequence.len();

    modified_position_sequence.into_iter()
        .enumerate()
        .skip(1) // p_0 is always 0.
        .map(|(j, &p_j)| factorial_ratio(n, j + 1) * p_j)
        .sum()
}

/// a! / b!
fn factorial_ratio(a: usize, b: usize) -> usize {
    debug_assert!(b <= a);
    (b + 1..a + 1).product()
}


#[test]
fn test_factorial_ratio() {
    assert_eq!(factorial_ratio(0, 0), 1);
    assert_eq!(factorial_ratio(2, 1), 2);
    assert_eq!(factorial_ratio(4, 2), 4 * 3);
}

pub fn factorial(n: usize) -> usize {
    (1..n+1).product()
}

#[test]
fn test_factorial() {
    assert_eq!(factorial(0), 1);
    assert_eq!(factorial(1), 1);
    assert_eq!(factorial(3), 6);
}

/// Compute the 'modified position sequence' without allocations.
/// The 'modified position sequence' is used to compute the index of a group (index of the permutation).
fn modified_position_sequence(modified_position_sequence: &mut [usize], position_sequence: &[usize]) {
    debug_assert_eq!(modified_position_sequence.len(), position_sequence.len());

    let n = position_sequence.len();

    debug_assert!((0..n).all(|i| position_sequence.contains(&i)), "invalid position sequence");

    let modified_sequence = (0..n)
        .map(|i| {
            let s_i = position_sequence[i];
            position_sequence.iter()
                .take(i)
                .filter(|&&s_j| s_j < s_i)
                .count()
        });

    // Write modified sequence into output buffer.
    modified_position_sequence.iter_mut()
        .zip(modified_sequence)
        .for_each(|(dst, v)| *dst = v);
}


#[test]
fn test_modified_position_sequence() {
    let modified_sequence = |position_sequence: &[usize]| -> Vec<usize> {
        let mut m = vec![0; position_sequence.len()];
        modified_position_sequence(&mut m, &position_sequence);
        m
    };

    assert_eq!(modified_sequence(&[2, 0, 3, 1]), vec![0, 0, 2, 1]);
    assert_eq!(modified_sequence(&[0, 1, 2, 3]), vec![0, 1, 2, 3]);
    assert_eq!(modified_sequence(&[3, 2, 1, 0]), vec![0, 0, 0, 0]);
    assert_eq!(modified_sequence(&[3, 2, 0, 1]), vec![0, 0, 0, 1]);
}

#[test]
fn test_group_index() {
    // Test that permutations map 1-to-1 to group indices.

    let group_index = |position_sequence: &[usize]| -> usize {
        let mut m = vec![0; position_sequence.len()];
        modified_position_sequence(&mut m, &position_sequence);

        group_index_from_modified_position_sequence(&m)
    };

    let n = 4;
    let permutations = (0..n).rev().into_iter()
        .permutations(n);

    let permutation_indices: Vec<_> = permutations
        .map(|perm| group_index(&perm))
        .collect();

    {
        let f: usize = (1..n + 1).product();
        assert_eq!(permutation_indices.len(), f);
        assert!((0..f).all(|i| permutation_indices.contains(&i)));
    }

    // assert_eq!(group_index(&[3, 2, 1, 0]), 0);
    // assert_eq!(group_index(&[3, 2, 0, 1]), 1);
    // assert_eq!(group_index(&[3, 1, 2, 0]), 2);
    // assert_eq!(group_index(&[3, 1, 0, 2]), 3);
    // assert_eq!(group_index(&[0, 1, 2, 3]), 23);
}