Skip to content

outlier_aware_minimum can crash if values[0] and values[not_nan-1] are both NaN #31

Description

@olekspickle

I believe there is a bug in outlier_aware_minimum function:

fn outlier_aware_minimum<F: Float>(mut values: Vec<(F, F)>) -> Option<F> {
    let num = values.len();
    let mut not_nan = num;
    // move NaN to the end
    for i in 0..num {
        let i2 = not_nan - 1;
        if values[i].0.is_nan() {
            values.swap(i, i2);
            not_nan -= 1;  // switches with END, then decrements
        }
        if i == i2 { break; }
    }
...
values.sort_unstable_by(|a, b| a.0.partial_cmp(&b.0).unwrap());

outlier_aware_minimum NaN filter for-loop swaps NaN to end but advances i afterward. When both values[i] and values[not_nan-1] are NaN, the swap exchanges NaN, leaving NaN at position i while not_nan decrements. i moves past it, NaN survives, partial_cmp(NaN) panics.

Would be nicer to just:

// effectively pushes NaNs to the end
values.sort_unstable_by(|a, b| {
    a.0.partial_cmp(&b.0).unwrap_or(Ordering::Greater)
});

I think there is another place this can blow up:
wynn_extrapolate produces NaN (line 123): d2.recip() - d1.recip() overflows when Richardson differences are near f64::MIN_POSITIVE. The clamp if err1 < tiny uses strict <, so values equal to tiny pass through, and cancellation of two ~4.5e307 reciprocals yields NaN.
You might want to account for that:

 let abserr = err1
            + err2
            + if converged {
                tol2 * F::from(10.).unwrap()
            } else {
                (result - e2[i]).abs()
            };
        if result.is_nan() || abserr.is_nan() {
            continue;
        }
        derivatives.push((result, abserr));
    }
    if derivatives.is_empty() {
        None
    } else {
        Some(derivatives)
    }

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Type

    No type

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions