To calculate the Kaiser window, I would like to implement the Modified Bessel Function of the first kind of order zero (I0). To make this work both with f32 and f64, I would like to use the num-traits crate.
This is my code so far:
use num_traits::{Float, NumAssignOps};
#[allow(non_snake_case)]
pub fn bessel_I0<Flt>(x: Flt) -> Flt
where
Flt: Float + NumAssignOps,
{
if !(x >= Flt::zero()) {
return Flt::nan();
}
let base = x * x / Flt::from(4).unwrap();
let mut power = Flt::one();
let mut coefficient = Flt::one();
let mut sum = Flt::one();
for j in 1.. {
power *= base;
if !power.is_finite() {
return Flt::infinity();
}
coefficient /= Flt::from(j * j).unwrap();
let addend = coefficient * power;
if addend.is_zero() {
break;
}
sum += addend;
}
sum
}
(Playground is currently buggy, so I'll provide a Playground link later)
I would like to know if my implementation is idiomatic or if it is bad style to implement it as I did. I particularly struggle a bit with the constants, such as 4, which I have to write like Flt::from(4).unwrap(), which uses num_traits::cast::NumCast::from.
Another thing that bugs me a bit is the num_traits::NumAssignOps trait. If I don't add it, I can't use the *= and += operators. But adding this bound makes all generic functions which use this function require that bound as well (which bloats up function signatures a bit). Why is NumAssignOps a separate trait? Does it have to do with certain implementations of floats which are immutable by design? I think they could somehow be wrapped to support +=. Afterall, Copy is a supertrait of num_traits::float::Float.
Also happy on any comments on numerical stability (if there are any).
The code seems fine except for the backwards check for x < 0 at the very beginning, as well as the the zero check of addend which could be moved up one line, directly checking whether coefficient is zero.
That's probably the only sensible way to create numeric constants of a generic type.
There's no such thing as an immutable type in Rust because mutability is a property of a binding, not a property of a type. I'm not an authority on num-traits, but it looks like the crate tries to follow the principle of composition/single responsibility by breaking up requirements into separate traits, instead of shoving everything in one giant trait.
Oh, that's indeed a mistake (don't even remember why I added it; I think because I looked at the graph and thought I0(x) is undefined for negative x, but that's wrong). Thanks for catching it.
Oh, right; at that point it's known that power is finite.
Here's the fixed function:
use num_traits::{Float, NumAssignOps};
#[allow(non_snake_case)]
pub fn bessel_I0<Flt>(x: Flt) -> Flt
where
Flt: Float + NumAssignOps,
{
let base = x * x / Flt::from(4).unwrap();
let mut power = Flt::one();
let mut coefficient = Flt::one();
let mut sum = Flt::one();
for j in 1.. {
power *= base;
if !power.is_finite() {
return Flt::infinity();
}
coefficient /= Flt::from(j * j).unwrap();
if coefficient.is_zero() {
break;
}
let addend = coefficient * power;
sum += addend;
}
sum
}
I found numeric_literals, but not sure if it's worth adding a dependency for a single constant (other than zero and one).
Hmmm, right. f32 and f64 don't come with "inner mutability" either, but you can reassign a variable of type f32 or f64. Sorry for my confusion.
But Float has a bunch of other supertraits as well:
My opinion is that generic functions working with primitive integer and floating-point types are more trouble than they're worth. It is much easier to split your function into two separate concrete functions bessel_I0_f32 and bessel_I0_f64. You don't need to bother with the absurd generic bounds that Rust requires, with by-reference and by-value implementations of operations, with possibly missing at the generic level functionality of concrete types (new methods on primitives are regularly added, are they also added to num-traits?). The precision of floating-point types may also affect your algorithms, and you can use the best possible algorithm in each case without trying to shoehorn it into a generic mold.
The implementation within your function really looks simple enough that it can be duplicated without any trouble. If you use later a more complex evaluation method, you can try extracting separate parts of it into private generic functions with more specific and simple trait bounds, or you can turn the body of the function into a macro (which is, admittedly, also a bit painful to work with).
By the way, I don't see a reason to use coefficient and power separately. You can directly compute the addend on the next step via the recurrence, and you will save a multiplication. I also don't see a good reason to check addend != 0. This is really only possible if x == 0, which is easier to handle separately at the start of the function. Either that, or you expect the loop run long enough to underflow the floating point type, which doesn't make any sense. It is way longer than required to compute the series up to any reasonable precision, and many of the terms in the series won't even affect the sum, since they will be smaller than the minimal representable difference from 1.
Also, if this is not just a learning problem and you actually need to compute bessel_I0 at many points, consider using a faster converging series than the naive Taylor expansion
You mean that certain numeric types should be used by-reference? (RefNum, NumRef) Does that also apply to working with num_traits::float::Float? I see that Copy is a supertrait of Float, though I would guess you always work by value? Or is that a false conclusion?
Sounds like you'd worry that 3rd party crates could easily be outdated in that matter?
I wonder how other crates that perform numeric operations deal with this? I looked at rustfft, and it provides an own trait FftNum for the API (that doesn't say anything about the internal implementation though, which I didn't look at).
I think providing two different functions for the API isn't so nice. But apparently rustfft didn't use num_traits::float::Float for some reason, but their own trait (that's only implemented for f32 and f64).
For such a small function, yes, but I will add more functions that will use this function. Should I duplicate them too? Eventually I'd end up with all code duplicated, I guess. That doesn't seem nice.
I thought this might help with numerical stability. Keeing the coefficient in a separate variable will ensure that coefficient will reach Flt::zero() quickly (in which case the iteration is stopped). I could consider putting coefficient and power into one variable (addend), which would then look like this:
use num_traits::Float;
#[allow(non_snake_case)]
pub fn bessel_I0<Flt: Float>(x: Flt) -> Flt {
let base = x * x / Flt::from(4).unwrap();
let mut addend = Flt::one();
let mut sum = Flt::one();
for j in 1.. {
addend = addend * base / Flt::from(j * j).unwrap();
if addend.is_zero() {
break;
}
if !addend.is_finite() {
return Flt::infinity();
}
sum = sum + addend;
}
sum
}
The problem here is that I didn't overlook how big j can practically grow. j * j needs to fit into an i32 to not cause overflow. Is it possible to prove that this algorithm is halting eventually (or before j * j overflows) for every input? Maybe this is more a theoretical problem and wouldn't cause any problem in practice.
But I now see you propose to provide a different termination condition anyway:
So I guess it's better to check if the result doesn't change anymore? Like this:
use num_traits::Float;
#[allow(non_snake_case)]
pub fn bessel_I0<Flt: Float>(x: Flt) -> Flt {
let base = x * x / Flt::from(4).unwrap();
let mut addend = Flt::one();
let mut sum = Flt::one();
for j in 1.. {
addend = addend * base / Flt::from(j * j).unwrap();
let old = sum;
sum = sum + addend;
if sum == old || !sum.is_finite() {
break;
}
}
sum
}
Regarding += and *=, I wonder if there is any downside in writing sum = sum + addend instead of sum += addend, beside a more verbose syntax. It would let me get rid of the NumAssignOps bound. Are there any floating point number implementations that wouldn't want to implement += or *= ?
I'm not an expert in numerical computation. Do you have any helpful links in that matter?
I mostly did this to learn how to use num-traits (and to see how well it works or if there is any trouble in writing generic numerical functions). I could also use an existing trait for my problem, but not sure if it's worth adding a big dependency.
One more question: What's your opinion on whether it's reasonable to violate the snake_case convention here by making the I in bessel_I0 capital?
In argmin I created a new trait ArgminFloat with all the trait bounds that I found reasonable for my use case and implemented this trait for all types which fulfill these trait bounds. I think this is similar to the approach in rustfft. You can find the implementation here. A simple float!(..) macro avoids having to write F::from(...).unwrap() everywhere (which still isn't great, but at least it leads to fewer distracting unwraps).
In linfa they did it in a similar way, but they implemented their Float trait on f32 and f64 explicitly and they added a convenient cast method.
I'm certainly not an authority on this topic, but I have found this to work pretty well.
So practically a "trait alias" (see also nightly feature) for all the traits you need. It seems to be common practice then, that each provider of floating point functions will provide their own trait, which specifies which functionality for floating point types is needed for a particular library.
I could create one for Float + NumAssignOps, I guess. But maybe desugaring x += y to x = x + y is not a big problem in my case. However, += (AddAssign) works by reference on the left hand side and might thus be faster for some implementations of Floats, I guess? That would be a good reason to provide my own trait.
I think there might be an advantage to implement the trait on f32 and f64 explicitly (and maybe even make it a sealed trait?). If you later need to add more supertraits to your own trait, this would break semver compatibility. See also what @afetisov said above:
I didn't consider semver compatibility in that matter.
Edit: Maybe a reasonable solution would be to simply document that new bounds might be added to whichever Float trait (alias) is being used, and that this doesn't break semver compatibility (i.e. you can only rely on f32 and f64 implementing it).
I see three and a half ways to go (which may be different for the API and internal/private functions):
Not use generics and only work on f32 and f64 directly
either by duplicating code manually
or through macros;
Using a trait alias (or trait with blanket implementation on stable Rust) that combines something like num_traits::float::Float with any other traits that are needed for a particular use case;
Using an own (possibly sealed?) trait that's directly implemented for f32 and f64.
Hmmmm… lots of possible trouble. Not sure what I should do.
Side question:
Does it also contain a (multivariate) Estimation of Distribution Algorithm as one possible strategy, and what do you think about those algorithms?
If I understand it right, the macro requires the type to be named F (source). I think that should be documented (currently doesn't seem to be documented). Maybe I'll create a private macro for my purpose. Something like:
use num_traits::Float;
macro_rules! flt {
($x:expr) => {
Flt::from($x).unwrap()
};
}
#[allow(non_snake_case)]
pub fn bessel_I0<Flt: Float>(x: Flt) -> Flt {
let base = x * x / flt!(4);
let mut addend = Flt::one();
let mut sum = Flt::one();
for j in 1.. {
addend = addend * base / flt!(j * j);
let old = sum;
sum = sum + addend;
if sum == old || !sum.is_finite() {
break;
}
}
sum
}
fn main() {
for x in [-1.0, 0.0, 0.5, 1.0, 2.0, 13.4, 500.0, 1000.0] {
println!("I_0({}) = {}", x, bessel_I0(x));
}
}
relative_kaiser_beta calculates the Kaiser window with an unknown scaling factor (which in some cases is okay because you perform some sort of normalization anyway).
kaiser_fn_with_beta returns a closure that returns values of a normalized (peak at 1.0) Kaiser window with a given beta. I used Currying here, so that I can write
kaiser_fn_with_beta(beta)(x)
but also store any intermediate values (that are the same for different x) as part of the closure:
let kaiser = kaiser_fn_with_beta(beta);
let y1 = kaiser(x1);
let y2 = kaiser(x2);
/* … */
I.e. this doesn't calculate scale more than once.
Is Currying a reasonable technique to use here? Or should I rather have a struct with a method?
Another problem: If I increase beta too much, I run into numerical problems: Playground. Again, I'm not experienced well-enough with numerics (or math in general) to come up with a way to work around it. Maybe I need to find a simplified representation of the scaled result (perhaps using a computer algebra software). In practice this doesn't matter, as beta should be reasonably small in most cases.
P.S.: Naïvely feeding this into WolframAlpha doesn't work.
P.S: Naïvely feeding this into WolframAlpha doesn't seem to work. I still get a result where a factor in the numerator grows quickly: β0, β2, β4, β6, …
I actually didn't consider sermver compatibility either. Implementing the trait explicitly of course solves that problem.
argmin currently does not contain such algorithms but they do sound very interesting and are certainly in the scope of the library. I'm happy to accept a PR if you (or anyone else) is interested in implementing
Not sure if implementing the trait explicitly solves the issue. Looking at your earlier example of linfa::Float, this is a public trait, and I believe it's non-sealed (though not sure if all its supertraits are really in public modules). I could use linfa with my floating point type for which I provide a linfa::Float implementation. Now let's assume std extends f32 or f64 with a new method. That method is added to linfa::Float as well. This would mean linfa would (theoretically) do a major version bump (see Rust reference). Adding an inherent item in contrast would be considered only potentially breaking (see Rust reference).
So the only solution would be to have a sealed trait or to specify that users of the crate may not rely on the trait being stable between new minor or patch versions (but only rely on being implemented by f32 and f64).
I mentioned these because they seem to be rarely used (though I'm no expert on that domain) but likely have some interesting properties, especially when functions behave a bit "weird". If I get again into programming something where I have an optimization problem, I'll take a look at argmin and consider adding these. Not sure if/when they are the best choice though.
