psychopath/sub_crates/halton/build.rs

// Copyright (c) 2012 Leonhard Gruenschloss (leonhard@gruenschloss.org)
//
// Permission is hereby granted, free of charge, to any person obtaining a copy
// of this software and associated documentation files (the "Software"), to deal
// in the Software without restriction, including without limitation the rights to
// use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies
// of the Software, and to permit persons to whom the Software is furnished to do
// so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included in
// all copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
// SOFTWARE.
//
// Adapted from Python to Rust and to generate Rust instead of C by Nathan Vegdahl

// Generate Rust code for evaluating Halton points with Faure-permutations for different bases.

use std::{env, fs::File, io::Write, path::Path};

/// How many components to generate.
const NUM_DIMENSIONS: usize = 128;

fn main() {
    let out_dir = env::var("OUT_DIR").unwrap();
    let dest_path = Path::new(&out_dir).join("halton.rs");
    let mut f = File::create(&dest_path).unwrap();

    // Init prime number array.
    let primes = {
        let mut primes = Vec::new();
        let mut candidate = 1;
        for _ in 0..NUM_DIMENSIONS {
            loop {
                candidate += 1;
                if is_prime(candidate) {
                    primes.push(candidate);
                    break;
                }
            }
        }
        primes
    };

    // Init Faure permutations.
    let faure = {
        let mut faure: Vec<Vec<usize>> = Vec::new();
        for b in 0..(primes.last().unwrap() + 1) {
            let perm = get_faure_permutation(&faure, b);
            faure.push(perm);
        }
        faure
    };

    // Write the beginning bits of the file
    f.write_all(
        format!(
            r#"
// Copyright (c) 2012 Leonhard Gruenschloss (leonhard@gruenschloss.org)
//
// Permission is hereby granted, free of charge, to any person obtaining a copy
// of this software and associated documentation files (the "Software"), to deal
// in the Software without restriction, including without limitation the rights to
// use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies
// of the Software, and to permit persons to whom the Software is furnished to do
// so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included in
// all copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
// SOFTWARE.

// This file is automatically generated.

// Compute points of the Halton sequence with with Faure-permutations for different bases.

pub const MAX_DIMENSION: u32 = {};
"#,
            NUM_DIMENSIONS
        )
        .as_bytes(),
    )
    .unwrap();

    // Write the sampling function
    f.write_all(
        format!(
            r#"
pub fn sample(index: u32, dimension: u32) -> f32 {{
    let mut index = index;

    match dimension {{"#
        )
        .as_bytes(),
    )
    .unwrap();

    // Write the special-cased first dimension
    f.write_all(
        format!(
            r#"
        // Special case: radical inverse in base 2, with direct bit reversal.
        0 => {{
            index = (index << 16) | (index >> 16);
            index = ((index & 0x00ff00ff) << 8) | ((index & 0xff00ff00) >> 8);
            index = ((index & 0x0f0f0f0f) << 4) | ((index & 0xf0f0f0f0) >> 4);
            index = ((index & 0x33333333) << 2) | ((index & 0xcccccccc) >> 2);
            index = ((index & 0x55555555) << 1) | ((index & 0xaaaaaaaa) >> 1);
            return (index as f32) * (1.0 / ((1u64 << 32) as f32));
        }}"#,
        )
        .as_bytes(),
    )
    .unwrap();

    // The rest of the dimensions.
    for i in 1..NUM_DIMENSIONS {
        let base = primes[i];

        // Based on the permutation table size, we process multiple digits at once.
        let mut digits = 1;
        let mut pow_base = base;
        while pow_base * base <= 500 {
            // Maximum permutation table size.
            pow_base *= base;
            digits += 1;
        }

        let mut max_power = pow_base;
        let mut powers = Vec::new();
        while (max_power * pow_base) < (1 << 32) {
            // 32-bit unsigned precision
            powers.push(max_power);
            max_power *= pow_base;
        }

        // Build the permutation table.
        let perm = (0..pow_base)
            .map(|j| invert(&faure, base, j, digits))
            .collect::<Vec<_>>();
        let perm_string = {
            let mut perm_string = String::new();
            for i in perm.iter() {
                let s = format!("{}, ", i);
                perm_string.push_str(&s);
            }
            perm_string
        };

        let mut power = max_power / pow_base;
        f.write_all(
            format!(
                r#"

        {} => {{
            static PERM{}: [u16; {}] = [{}];"#,
                i,
                base,
                perm.len(),
                perm_string
            )
            .as_bytes(),
        )
        .unwrap();

        f.write_all(
            format!(
                r#"
            return unsafe {{(
                *PERM{}.get_unchecked((index % {}) as usize) as u32 * {}"#,
                base, pow_base, power
            )
            .as_bytes(),
        )
        .unwrap();

        // Advance to next set of digits.
        let mut div = 1;
        while power / pow_base > 1 {
            div *= pow_base;
            power /= pow_base;
            f.write_all(
                format!(
                    r#"
                + *PERM{}.get_unchecked(((index / {}) % {}) as usize) as u32 * {}"#,
                    base, div, pow_base, power
                )
                .as_bytes(),
            )
            .unwrap();
        }

        f.write_all(
            format!(
                r#"
                + *PERM{}.get_unchecked(((index / {}) % {}) as usize) as u32
            )}} as f32
            * (0.999999940395355224609375f32 / ({}u32 as f32)); // Results in [0,1).
        }}
        "#,
                base,
                div * pow_base,
                pow_base,
                max_power
            )
            .as_bytes(),
        )
        .unwrap();
    }

    f.write_all(
        format!(
            r#"
        _ => panic!("Halton sampling: exceeded max dimensions."),
    }}
}}
    "#
        )
        .as_bytes(),
    )
    .unwrap();
}

/// Check primality. Not optimized, since it's not performance-critical.
fn is_prime(p: usize) -> bool {
    for i in 2..p {
        if (p % i) == 0 {
            return false;
        }
    }
    return true;
}

/// Computes the Faure digit permutation for 0, ..., b - 1.
fn get_faure_permutation(faure: &Vec<Vec<usize>>, b: usize) -> Vec<usize> {
    if b < 2 {
        return vec![0];
    } else if b == 2 {
        return vec![0, 1];
    } else if (b & 1) != 0 {
        // odd
        let c = (b - 1) / 2;

        return (0..b)
            .map(|i| {
                if i == c {
                    return c;
                }

                let f: usize = faure[b - 1][i - ((i > c) as usize)];
                f + ((f >= c) as usize)
            })
            .collect();
    } else {
        // even
        let c = b / 2;

        return (0..b)
            .map(|i| {
                if i < c {
                    2 * faure[c][i]
                } else {
                    2 * faure[c][i - c] + 1
                }
            })
            .collect();
    }
}

/// Compute the radical inverse with Faure permutations.
fn invert(faure: &Vec<Vec<usize>>, base: usize, mut index: usize, digits: usize) -> usize {
    let mut result = 0;
    for _ in 0..digits {
        let remainder = index % base;
        index = index / base;
        result = result * base + faure[base][remainder];
    }
    return result;
}