Firwin

[SpookyYomo]

This PR aims to implement Firwin.

Blockers before merging:

• Decision on Err implementation in Library
• Special Functions
• Arg Struct pattern
• Documentation fixes #61

---

Original Message:

Currently opening as a draft PR to have eyes on the code before it grows even larger into eventually implementing firwin.

Questions I guess I need answers about specific to this repository to ensure the PR is decent:

1. What exactly in a function is allocating/otherwise for me to determine if the function name should be *_st or *_dyn?
2. File structure... As a first time contributor, it feels like it's just having many files and then being flattened by mod.rs to follow the scipy's submodule convention.... Individual files otherwise don't really make sense, and finding a function is best done either via grepping or through the [Source] button as created by cargo docs.
3. Doc's index.html: Maybe we should structure the Doc abit more clearly similar to Scipy's documentation where there's a clear header of Convolution, B-splines, Fitting, ...? I think this is achieveable with //!? Or I might be mistaken.
4. Unsafe code: I noticed some other functions using unsafe { F::from(-1).unwrap_unchecked) } or something similar. Is this the norm for "constants" which we can guarantee to be safe, and so I too should be doing this?

I will continue to force-push onto this branch with any suggestions such that each commit is consistent with its commit title and message.

[SpookyYomo]

a92506d I have added error types after looking through how the Rust ecosystem generally does errors. This might still not be ideal and needs more work.

We should leave panics instead to Python-facing side of the library.

[SpookyYomo]

a1783ca Implementing as a trait for different types was the only way I could see to ensure that the Float F was not pinned/polluted with i64: From<F> in type restrictions in everything that had to use i0.

Similar to scipy, we used Cephes, however, it was found that i0 of f32 as from the provided C library has an issue, so the corresponding test is ignored.

A macro could probably be done to implement all the functions for Vec and Array<F, D> after having them defined in the trait and manually done for f32/f64 (due to function names being a tad scattered).

[SpookyYomo]

4598cbd There's still some redundancy in argument specification of Firwin. But I'm not sure how to resolve it.

Nonetheless it behaves to expectation now.

[SpookyYomo]

Result

[SpookyYomo]

Result

[trueb2]

🎉 Awesome! This will take a little bit to review, but at first read through this looks like a great addition!

You have a ton of great comments.

Discussion

Your concerns over patterns are valid. I have avoided adding an idiomatic error type, but it is probably time to take a look at the bigger picture. Some of the idiomatic patterns are clashing with the Python interface such as Error vs throw and kwargs vs args Struct.

Error w/ and w/o alloc

The Error type looks useful What are your thoughts on making only certain members available w/o alloc enabled?

Docs and Readme

I think we are ready for a revamp of the README and the Docs index. We are starting to have a lot of features that could be easier to navigate and utilize. Some of the Rust interfaces that mirror python are also a little obtuse. They kinda need better docs or a high level wrapper interface.

Guidelines and Patterns

With so many new features added, it may be time for some better structure around patterns for the project to pursue:

• Error patterns
• Arg patterns
• Func naming patterns
• Documentation patterns
• File structure
• Unsafety

Maybe we can get docs and guidelines in place for v0.5.0.

For now, all of the alloc gated functions should be _dyn according to the current guidelines.

I am open to changes in the guidelines. For example, cargo could not detect if a missing function was feature gated when the project started, but now it would give a good compiler hint to enable alloc.

SciPy, Cephes, Special

I will look into putting together some prior arts and references that contributors and users could refer to. I don't really like the special-fun dependency and its BSD license. What are your thoughts on vendoring or copying https://github.com/stainless-steel/special into this project and extending for usage within sci-rs?

Unsafe Usage

I am 100% okay with unsafe if the situation requires. For example, a lot of the code in sci-rs runs on thumbv7em-none-eabihf, where there isn't good branch predictor + instruction pipeline performance. Allocating is expensive; however it is not hammered inside a loop. If we can reasonably elide an unwrap in a tight loop, I am for it. If the unwrap isn't in hot code, i.e. filter design, then maybe it should be converted into a sci_rs::Error and bubbled up.

File Size

Breaking up the big mod.rs files seems to speed up rust-analyzer. I think we might be doomed to ten thousand files and nested mods. If you have a solution or recommendation, I am for trying it out for a while.

Small file pros

• Faster rust-analyzer
• Better unit test isolation

Small file cons:

• Differs from scipy
• Harder to navigate (100s of mod.rs files)

[SpookyYomo]

The Error type looks useful What are your thoughts on making only certain members available w/o alloc enabled?

All returnable Errors should be limited to the scope enabled by the feature flag used during the compilation flag, at a minimum. Additional Errors being accessible through use sci_rs::error::* should be weighed against situations where it may prove useful. Otherwise, as a new Rustacean, I'm not familiar with best practices with regards to library level Errors.

I think we are ready for a revamp of the README and the Docs index. We are starting to have a lot of features that could be easier to navigate and utilize. Some of the Rust interfaces that mirror python are also a little obtuse. They kinda need better docs or a high level wrapper interface.

Parrallelising with this thought, I spent a bit of today mulling over providing if firwin_dyn should have a convenience macro that has Python styled args/kwargs firwin that then calls firwin_dyn under the hood, and returns firwin_dyn as is. This would also be comparably in line with the long-lasting goal as presented in the top level README.

Nonetheless, we should still try lean into Rust's compiler gurantees by having slightly more terse functions that are sci-rs specific that may not necessarily mirror scipy very well. For example, the last thing we would want is to take a Vector of 2-tuples for args/kwargs mirroring, which can be a nightmare where the 2nd argument is a mix of scalars, arrays, and hashmaps; we then let the macro mirror the Python styled variadic function, and perhaps even a sciprs/ folder to an scipy-expected interface from Python, which can also double as a demo on how to utilise said functions for sci-rs, and also is the one that does raise ValueError or as is similar to what scipy does for the PyO3-facing sciprs/.

For now, all of the alloc gated functions should be _dyn according to the current guidelines.

Just to double check, should the function get_window under the GetWindow trait instead by get_window_dyn then?

SciPy, Cephes, Special

I will look into putting together some prior arts and references that contributors and users could refer to. I don't really like the special-fun dependency and its BSD license. What are your thoughts on vendoring or copying https://github.com/stainless-steel/special into this project and extending for usage within sci-rs?

We would then have to wait for stainless-steel/special#18. I generally think the biggest struggle with regards to this is that I am only aware of the $$I_\alpha(z) = \sum_{m=0}^\infty \frac1{m!\Gamma(m+\alpha+1)}\left(\frac x2\right)^{2m+\alpha}$$ definition... which isn't computationally friendly (the Cephes implementation looks like black magic to me, I can't even fathom how Mathematica does it).

... it is not hammered inside a loop. If we can reasonably elide an unwrap in a tight loop, I am for it.

Would the following be considered as unwrap in a tight for loop? W::from(2).unwrap() in f901379#diff-b6a197656fff7d7ae5121775018d84d7cf2c05aa17d8adeb3aad10aefbbc9073R194

File Size

If you have a solution or recommendation, I am for trying it out for a while.

After having been immersed into rust through having implemented this for a while, I realise that it probably in fact not the wisest to mirror the file structure of scipy either, especially if I were to mirror the ufunc expectation (as I have had for i0) for several of the functions. Furthermore, there is nothing a little fd or rg cannot do to find the location of a function body, or even just a simple jump to definition, when utilizing the documented example.

The only sticking point I have to say is that I would probably expect
4598cbd#diff-9f66d7d9b9488cb41fa0ffe1a1ae1beac7ea115580e54d0005adcbf191a0752dR198
to have instead been use sci_rs::signal::filter::firwin, alot of our internal mod to segment sections of code could(?) be left as a vestige, but we should have use statement for the public interface be minimally identical to the scipy accessible, which would be import scipy.signal.filter.firwin.

[SpookyYomo]

What are your thoughts on https://www.netlib.org/slatec/src/besi.f from https://doi.org/10.1145/355719.355726? It seems to be public domain, but it certainly demands a Fortran to Rust translation.
There's also a few others as mentioned https://math.nist.gov/nesf/Sect3-3/Sect3-3.html.

Edit:
I've found https://github.com/Axect/puruspe, which has a f64 implementation for BesselI, but that means we would just taint the code base with generic restrictions of f64: From<F> everywhere for now.

[SpookyYomo]

As mentioned in #61, this PR is now split hopefully for easier understanding, and is rebased against the head of #61.

[matthew-hennefarth]

As an aside on the special math functions, I have a special-rs repository where I have functions for:

• Factorial, Double Factorial, K-Factorial
• Combinations and Permutations
• Bernoulli, Tangent, and Secant numbers (integers)
• Gamma Function (real and complex)
• Error Function (real)

A lot of them are translations either from cephos (the scipy backend) or boost (from C++).

What has prevented me from making more PR's here is that I was unsure how to deal with the function arguments quite honestly. Since, for example the gamma in scipy can take either real or complex values, one uses different algorithms to evaluate the gamma function. It is also a bit clunky to always have to pass a Complex{} object to the function call (ie promote a real -> complex number) and always return a complex value.

One can take inspiration from BLAS function signatures by using dgemm for f64 or sgemm for f32, and zgemm for Complex<f64>, but this then differs from the python interface. Alternatively, what I tried to do in my implementation was make a trait (such as Gamma trait (which cannot be used anymore because of a standard library Gamma trait soon apparently)) and then using macros to implement so that users can do

assert_eq!(4.0_f32.gamma(), 6.0); // Gamma(4) = 3!
assert!((0.0_f64).gamma().is_nan()); // Gamma(0) is undefined

it is then trivial to make a generic function gamma such that

pub fn gamma<T>(z: T) -> T 
where 
  z: Gamma
{
    z.gamma()
}

The possible downside of this approach is that Gamma is no longer generic over RealField for example but has to be implemented (albeit perhaps trivially).

[matthew-hennefarth]

Also, as an aside on all of the F:from().unwrap() for special numbers in the algorithms, one option is to make an interface that binds these special constants as static vars. For example, in my r_gamma implementation, I have special polynomial constants. Instead of littering .from().unwrap() everywhere, I generate a private RealGammaConst trait and implement it for the types I need (f32 and f64) (see here)

pub(crate) trait RealGammaConsts: Sized {
    const MIN_TO_USE_STIRLING: Self;
    const P: [Self; 7];
    const Q: [Self; 8];
    const INTERVAL: [Self; 2];
    const MIN_USE_SMALL: Self;
}


macro_rules! impl_realgammaconsts {
    ($($T: ty)*) => ($(
        impl RealGammaConsts for $T {
            const MIN_TO_USE_STIRLING: Self = 33.0;
            const P: [Self; 7] = [
                1.60119522476751861407E-4,
                1.19135147006586384913E-3,
                1.04213797561761569935E-2,
                4.76367800457137231464E-2,
                2.07448227648435975150E-1,
                4.94214826801497100753E-1,
                9.99999999999999996796E-1,
            ];
            const Q: [Self; 8] = [
                -2.31581873324120129819E-5,
                5.39605580493303397842E-4,
                -4.45641913851797240494E-3,
                1.18139785222060435552E-2,
                3.58236398605498653373E-2,
                -2.34591795718243348568E-1,
                7.14304917030273074085E-2,
                1.00000000000000000320E0,
            ];
            const INTERVAL: [Self; 2] = [2.0, 3.0];
            const MIN_USE_SMALL: Self = 1.0E-7;
        }
)*)
}

impl_realgammaconsts! {f32 f64}

Honestly, I do not know what the performance (or the maintenance) implications are, but is a possible solution where .from().unwrap() are used in for loops.

Again, this limits the ability to generically implement functions over any field type (and in particular, over types like BigInt or BigFloat.

[trueb2]

Lots of great discussion in this thread! It has been very thought provoking, so you can jump to the conclusion for a TLDR.

Thoughts on Direction

I have some overarching thoughts and would appreciate your feedback.

1. Features are most important in the short term.
   1. People want stuff to work for them.
   2. I don't want to break interfaces, but I expect it is inevitable.
2. Finalizing error handling must be resolved asap.
3. Deciding a special number representation scheme and primitive conversions is essential for further implementations.
   1. We will need to be careful and iterate.
4. I perceive that we only need to support floating point types the Rust compiler supports (f32, f64, and Complex<_>).
   1. f16 and f128 is on the horizon but who knows about the numerical stability of some of the math.
   2. BigInt, BigFloat, etc. are low on the priority totem pole.

Responses

The only sticking point I have to say is that I would probably expect
4598cbd#diff-9f66d7d9b9488cb41fa0ffe1a1ae1beac7ea115580e54d0005adcbf191a0752dR198
to have instead been use sci_rs::signal::filter::firwin, alot of our internal mod to segment sections of code could(?) be left as a vestige, but we should have use statement for the public interface be minimally identical to the scipy accessible, which would be import scipy.signal.filter.firwin.

I agree with you. I have been thinking about reorganizing imports or establishing a prelude pattern.

What are your thoughts on https://www.netlib.org/slatec/src/besi.f from https://doi.org/10.1145/355719.355726? It seems to be public domain, but it certainly demands a Fortran to Rust translation.
There's also a few others as mentioned https://math.nist.gov/nesf/Sect3-3/Sect3-3.html.

I'm all for translating or vendoring in some dependencies that are pretty relevant. I think we have experienced friction because of all the great design points under consideration here.

Just to double check, should the function get_window under the GetWindow trait instead by get_window_dyn then?

My commitment to _st and _dyn has been waning of late. It is a function coloring problem that is guarded by the compiler. Many new features I have needed lately involve explicit interface generics or allocation. As our feature set has been growing, even keeping a commitment to _st vs _dyn is waning as the cfg feature alloc attributes are a wart that just keeps spreading.

I would appreciate thoughts from others. I don't plan on getting rid of feature alloc, but I don't ever disable it anymore.

Would the following be considered as unwrap in a tight for loop? W::from(2).unwrap() in f901379#diff-b6a197656fff7d7ae5121775018d84d7cf2c05aa17d8adeb3aad10aefbbc9073R194

It could be a tight loop. In my usages, filter design is infrequent, but it could be frequent for some people. I would want to move the unwrap of 2 outside of the loop.

Also, as an aside on all of the F:from().unwrap() for special numbers in the algorithms, one option is to make an interface that binds these special constants as static vars.

I like this approach. consts in trait impls are a good use of Rust language features here.

Reflections

I recognize that some of the loose implementation thus far is creating some confusion. I feel this is the result of chasing features and performance whilst punting hard interface decisions. In a way, I am happy with this as we do have features and some critical interface design has fallen through as Rust support for const generics fizzled out in ways that were necessary for the original design.

I don't want this library to turn into a macro library (like a C all header library or C++ template library). Beyond the increase in compile times, I find those difficult to navigate, test, and contribute to ... However, macros for trait and function impls seem to be the best solution for both interface, constants, and conversion issues. Using goto definition and ending up at a (nested or proc) macro call site stinks, but we will have a lot of repeated code otherwise.

sci-rs has tried to be pythonic in its function interfaces. At times, it has followed algorithms written for C or Fortran. What it has never really tried to do is be Rusty. Maybe we should try ...

Another recent interface example

Recently, I needed the performance of filtering to be better. sci-rs is all about having a portable, performant implementation, but even the accelerated custom interface we have now wasn't enough. It used slices without a compile-time know length. I needed the compiler to spit out the best assembly — to unroll a loop and minimize branching that predicted poorly on Cortex-M4F. I haven't PR'd this into the repo because I am not sure what to make of this interface divergence and what it could mean for other interfaces like FIR filtering.

#[inline(always)]
pub fn sosfilt_fast32<const N: usize, const M: usize>(y: &[f32; N], sos: &mut [Sos<f32>; M], z: &mut [f32; N]) {
    _sosfilt32(y, sos, z);
}

This took an interface that was slightly dynamic and made it static. We could keep adding traits, but when should we stop adding traits and just copy paste (macro) an implementation? This interface change led to significant improvements but would live in the codebase as a static interface parallel to the dynamic slice length interface. This is reduces maintainability and compile time but unlocks optimal performance.

If we go looking at prior arts in other languages, many interfaces are happy to go terse+macro/codegen/template to get exactly the right version of the code. To achieve the desired behavior for each scipy feature, Rust's language features push us either towards:

1. static or dynamic dispatch over a trait per feature to let the compiler choose.
   1. a large divergence in interface
   2. easy to use for a Rustacean
   3. easiest to maintain
2. (proc)macro compile time code generation to get the exact implementation required.
   1. a large divergence in interface
   2. hard to maintain
   3. easiest to use
3. a heterogenous composition of shared traits, macro implementations of traits, and pub functions
   1. the hardest to maintain
   2. this is our current approach

Possible Interface Direction

Are we heading towards pythonic callsites via macro-style interface?

let a = gamma!(x);
let a = sosfilt!(f32, N, y, z, M, sos);

One can take inspiration from BLAS function signatures by using dgemm for f64 or sgemm for f32, and zgemm for Complex, but this then differs from the python interface.

I've had my embedded C hat on for a little bit, and there is a grace to the simplicity of some consts or macros. It can be quite easy for a user to navigate calling dgemm or sgemm. I like it for its usage simplicity but it would require some (proc)macros to get right. Compile times and maintainability may get out of hand. LLM agents may make this a tenable route.

Alternatively, what I tried to do in my implementation was make a trait (such as Gamma trait (which cannot be used anymore because of a standard library Gamma trait soon apparently)) and then using macros to implement so that users can do

I'm not a huge fan of macros to blanket trait implementations, but it is very common in rust libs. As you have in your example, there are problems including conflicts, compile time, and the orphan rule. nalgebra in particular suffers some serious compile time cost. There have been some serious issues with trait changes propagating through dependencies like our dependency on RealField.

I perceive the most efficient compile time and assembly to come from writing out the interface options. Using fewer Rust language features drastically improves the compile time as usage of traits and macros generates a lot of extra work for the compiler. This is the opposite of the pythonic interface direction.

Another Possible Interface Direction

Are we heading towards Rusty interface diverging from scipy? Here, each scipy function becomes a trait and one or more impls.

struct SosfiltTiled<F: SciReal, N: const usize, M: const usize> { ... }
impl<F: SciReal, N: const usize, M: const usize>  Sosfilt for SosfiltTiled <F, N, M> {
// do the specialized implementation
} 

struct SosfiltIterator<F: SciReal > { ... }
impl<F: SciReal> Sosfilt for SosfiltIterator<F> {
// do the specialized implementation
} 

let filter = SosfiltTiled { ... }
let a = filter.sosfilt();

This leans into Rust design, diverging from prior arts with C by relying on one name and diverging from python by removing the monadic control flow and naming. Each action requires capturing context, then using a sci-rs trait, allowing for optimal implementation but a huge number of trait implementations for the compiler to sort out.

It is a very flexible approach, but it will break nearly every interface we have. Heading in this direction means re-writing sci-rs

Crate structure

I think we are eventually going to be a workspace structured similarly to polars. (especially if we are going towards a Rusty interface)

• sci-rs: The umbrella crate that people would include in the Cargo.toml
• sci-rs-core: Required stuff like errors and primitive traits
• sci-rs-special: All the constants and glue to use them
• sci-rs-linalg: All the shared linear algebra stuff
• sci-rs-signal: All the scipy signal stuff
• sci-rs-stats: All the scipy stats stuff
• etc
• sciprs: The python package that people could pip install

Conclusion

I think we should start heading towards the Rusty interface that diverges from the pythonic style.

We can break up scipy functions into a struct and a trait. This allows defining constants and specialized implementations that unroll/tile/vectorize/etc. I think we need to establish our own floating point value trait like SciReal: RealField + RealGamma + .... A crate reorganization will be required to keep compile times down. The prelude and mod/use structure needs to update. The callsites for the users need to be more consistent across the project.

TLDR

My disposition is to take the following steps:

1. Merge this PR.
2. Update sci-rs error type and numerical abstractions
3. Begin crate reorganization, starting with sci-rs-special
   1. remove the BSD license dependency before the next release
4. Begin using struct + trait per feature to allow defining constants and specializations
   1. refactor existing implementations behind traits functions
5. Revamp docs around traits
   1. update contributing docs and CI
6. Release v0.5.0

Feedback is welcome! Thanks again!

[matthew-hennefarth]

Thoughts on Direction

1. Features are most important in the short term.
   i. People want stuff to work for them.
   ii. I don't want to break interfaces, but I expect it is inevitable.

I agree. Having a function which does the calculation and is tested, but maybe not having the best interface, is better than no function at all.

4. I perceive that we only need to support floating point types the Rust compiler supports (f32, f64, and Complex<_>).
   i. f16 and f128 is on the horizon but who knows about the numerical stability of some of the math.
   ii. BigInt, BigFloat, etc. are low on the priority totem pole.

That's fine. I think its good to have a scope of what numeric types are supported by default.

Another recent interface example

#[inline(always)]
pub fn sosfilt_fast32<const N: usize, const M: usize>(y: &[f32; N], sos: &mut [Sos<f32>; M], z: &mut [f32; N]) {
    _sosfilt32(y, sos, z);
}

This just seems wrong. While I am not an expert in Rust code, from a C++ perspective, this seems like it would create unnecessarily long compile times if sosfilt_fast32 is called on various different array lengths. In general, I don't think it is smart to make a function generic over an array length unless there is some very specific reason to. This is my opinion without any perspective on the performance consideration.

Another Possible Interface Direction

Are we heading towards pythonic callsites via macro-style interface?

I would avoid the macro-style interface to have variadic function parameters. I believe your conclusion on a Rusty interface makes the most sense.

Are we heading towards Rusty interface diverging from scipy? Here, each scipy function becomes a trait and one or more impls.

I think this makes the most sense in the Rust ecosystem. If I think about what a Rust end-user, one would like to make a function that generically accepts something that could call sosfilt or gamma or factorial on. From the python perspective, I don't believe there is any harm in having a small python-rust interface that essentially forwards python arguments to the corresponding rust function. This is already how scipy operates and interacts with Cephos for example.

My main issue is the necessity of having a struct associating each 'input' type to a function. But maybe this makes sense from the perspective of the signal functions. From the perspective of the special math functions, a lot of times just numbers (floats, ints, or complex) are the input type; however, floats (f32 and f63) will have 1 implementation, ints will have their own, and complex another implementation. So would you envision (taking the Gamma function as the example) the code structures like:

impl Gamma for f32 {
    fn gamma(self) -> Self {
        // some implementation
    }
}
impl Gamma for f64 {
    fn gamma(self) -> Self {
        // basically same implementation as f32?
    }
}

Though, I guess if there existed a SciReal (or just using RealField one could do:

fn r_gamma<F: RealField>(n: F) -> F {
  // Gamma implementation for real-valued
}

impl Gamma for f32 {
  fn gamma(self) -> Self {
    r_gamma(self)
  }
}

impl Gamma for f64 {
  fn gamma(self) -> Self {
    r_gamma(self)
  }
}

But then this is very similar to what I was already doing, except I was just using a macro to write the

impl Gamma for f32 {
  fn gamma(self) -> Self {
    r_gamma(self)
  }
}

impl Gamma for f64 {
  fn gamma(self) -> Self {
    r_gamma(self)
  }
}

part out. It is on this last point that I might be confused.

Conclusion

Overall, I think I generally agree with your viewpoints. I am not too stuck in my opinions and am happy to help where possible. I will defer to you and others about the general structure on the crate(s) for this project.

Also, should perhaps this discussion be moved then to an issue? Rather than on this PR?

[SpookyYomo]

Hi, please feel free to open as an issue if there is a need for more discussion. Clearly we need a rework on the relevant interfaces. Otherwise, I personally don't really have much more to add with regards to interfaces as mentioned by both. Personally, no objection to the proposed direction forward.

The code will still be in this branch and can easily be reworked to satisfy the relevant new interfaces if we aren't merging; or we can just take it in as is and fail-forward as Jacob has proposed. The only concern I have thus far is: Is having BSD code in the repo for some span suitable; even if not in releases? I'm not familiar with licensing (and hence why I just jumped the gun and imported Cephes to begin with, creating this debacle, my apologies.)

Interfaces aside, I found a book called Numerical Recipes by William T. Vetterling being conducted by my university that does have a domain "folding" similar to the original Bessel implementation (https://doi.org/10.1145/355719.355726). Perhaps I could cook that up in the interim? It's quite easily googleable.