Hi,
I implemented the stochastic universal sampling algorithm. One variant works with f32 and one with f64 and for some reason, I can provoke an index out of bounds panicks with the f32 variant, which is almost identical with the f64 variant.

pub fn random_choice_f32<'a, T>(samples: &'a [T], weights: &[f32], n: usize) -> Vec<&'a T> {
if weights.len() == 0 || n == 0 {
return Vec::new();
}
let sum: f32 = weights.iter().fold(0.0, |acc, &i| acc + i);
let spoke_gap: f32 = sum / n as f32;
// next_f32() ∈ [0.0, 1.0)
let spin = thread_rng().next_f32() * spoke_gap;
let mut i: usize = 0;
let mut accumulated_weights = weights[0];
let mut choices: Vec<&T> = Vec::with_capacity(n);
let mut current_spoke: f32 = spin;
println!("sum: {}, gap: {}", sum, spoke_gap);
for j in 0..n {
while accumulated_weights < current_spoke {
println!("accumulated_weights: {}, current_spoke: {}, spoke number: {}, i: {}", accumulated_weights, current_spoke, j, i);
i += 1;
accumulated_weights += weights[i];
}
choices.push(&samples[i]);
current_spoke += spoke_gap;
}
choices
}

The f64 variant is absolute identical, except it uses f64.
Following test provokes an index out of bounds panick:

#[test]
fn test_random_choice_f32() {
let capacity: usize = 500;
let mut samples: Vec<usize> = Vec::with_capacity(capacity);
let mut weights: Vec<f32> = Vec::with_capacity(capacity);
for i in 0..capacity {
samples.push(i);
weights.push(i as f32);
}
let number_choices = 10000;
let choices = RandomChoice::random_choice_f32(&samples, &weights, number_choices);
assert!(choices.len() == number_choices);
let mut weight_counter = BTreeMap::new();
for choice in choices {
let counter = weight_counter.entry(choice).or_insert(0);
*counter += 1;
}
let mut last_value: usize = 0;
for (_, value) in &weight_counter {
assert!((last_value as i32 - (*value) as i32).abs() <= 2);
last_value = *value;
}
}

The f64 test never panicked and the pseudo code you find in wiki, looks almost identical.
Anyone any idea? You can find the code on github.

I'm assuming that the reduced precision of f32 is causing some index you are calculating to round up to one-past-the-array-end, and the f64 version would have the same issue but it's triggered 2 ** (52 - 24) times less often.

Actually I am just comparing floats with the less operator and incrementing the index counter with i += 1 which has the type usize.
If you uncomment the println!(), you'll see the output with cargo test -- --nocapture.

That won't fix the crash in all cases though, because current_spoke == sum might be true. You could change it to while current_spoke < sum, which might work, but would still have the n-1 problem.

What I suggest instead is to treat samples which are just-past-the-end as if they were just-before-the-end, by modifying the inner loop condition:

while i < weights.len() && accumulated_weights < current_spoke {
println!("accumulated_weights: {}, current_spoke: {}, spoke number: {}, i: {}", accumulated_weights, current_spoke, j, i);
accumulated_weights += weights[i];
i += 1;
}

Side note: if your population has more than a few million individuals, doing this with all f32 variables will have enough roundoff error to substantially influence the result. If your population is that large, it'd probably be better to make all the local variables f64 and do weights[i] as f64 in all three places.

Thanks. I had an old version where I used the modulo operator, so if the current spoke is greater than the sum, it gets smaller with the modulo operator, but the pseudo code of wikipedia has no such check and in my opinion the accumulated spoke can't get bigger than sum, because of this:

let spoke_gap: f32 = sum / n as f32;
// next_f32() ∈ [0.0, 1.0)
let spin = thread_rng().next_f32() * spoke_gap;

So spin is smaller than spoke_gap and that's why current_spoke should always be smaller than the sum.
If the last spoke has the position sum - spoke_gap and you add spin, it is still less then the sum, because spin is less than spoke_gap, right?
So, there must be some implementation error.

The pseudo code on Wikipedia assumes infinitely accurate numbers. Your numbers are f32, which round off at 24 bits after the binary point. They can round in either direction, and if the random number is very close to one, the rounded result of current_spoke at the end of the loop could very well be greater than spin.