I think the project should be called "NIntegrate".
BTW, that is not a serious suggestion; it is just that Wolfram Language (aka Mathematica) has both `Integrate` and `NIntegrate` for symbolic and numeric integration, respectively.
> Does the function oscillate over the region of integration? If it does, then make sure that the step size is chosen to be smaller than the wave length of the function.
One way to think about is that these techniques work by integrating exactly the polynomial that interpolates the function where you're sampling it, so you need to resolve the features of the function to get good accuracy.
The Nyquist sampling theorem is of course proved by considering a Fourier transform, which is given by an integral, so the relation to integration in general should not be surprising.
It looks a bit sloppy to hardcode so many constants in a single file: `src/gauss_quadrature/legendre.rs`. Isn't it possible to generate them with the help of rust macros in the same way Julia uses metaprogramming?
Gaussian quadrature points are typically solved numerically. There's a good chance these ultimately came from a table.
Additionally, compile time floating-point evaluation is limited. When I looked around recently, I didn't see a rust equivalent of gcem; any kind of transcendental function evaluation (which finding Gaussian quadrature points absolutely would require) would not allow compile-time evaluation.
It's relatively straightforward to build transcendental functions out of the basic operations and the stdlib support will eventually get there, but rust's float story is still a work in progress. They're trying to do things more properly and document semantics better than C and C++ have.
I was under the impression that macros can execute arbitrary code, surely some FP would not be big problem. And if not macros then build.rs script certainly could do the trick.
build.rs can definitely execute arbitrary code, which means that a lot of places (including, IIRC crates.io) will forbid crates that rely on build.rs. I ended up refactoring my build.rs into a separate sub-application in finl_unicode that built data tables which are then checked into git and used pre-built. I include the sub-app in the source code so that anyone with access to the repo could continue development if I were to die tomorrow.
I ran into some issues with crates.io and my build.rs when I first released the crate, although it’s long enough ago, that I don’t remember exactly what the issue was. It might also have been that the build.rs script downloaded raw data files from unicode.org
crates.io doesn't care what your build.rs does because it doesn't try to compile your code, neither now or ever in the past. There would be no point in it trying to compile your code; there are lots of crates that are bindings to C libraries that crates.io's builders can't be expected to have, crates that target architectures which crates.io can't be expected to have builders for, etc.
You wouldn't have to recompile them every time. What if you didn't necessarily use macros but auto-generated it in a file that you keep separate from the other code at the bottom?
What I would do in these cases is to define the general computation function, but special-case it to return the hard-coded value for specific common inputs if it's being evaluated at compile time. Then add a test to verify both behaviors.
It seems like it is lacking the functionality R's integrate has for handling infinite boundaries, but I suppose you could implement that yourself on the outside.
For what it's worth,
use integrate::adaptive_quadrature::simpson::adaptive_simpson_method;
use statrs::distribution::{Continuous, Normal};
fn dnorm(x: f64) -> f64 {
Normal::new(0.0, 1.0).unwrap().pdf(x)
}
fn main() {
let result = adaptive_simpson_method(dnorm, -100.0, 100.0, 1e-2, 1e-8);
println!("Result: {:?}", result);
}
prints Result: Ok(1.000000000053865)
It does seem to be a usability hazard that the function being integrated is defined as a fn, rather than a Fn, as you can't pass closures that capture variables, requiring the weird dnorm definition
How many evaluations of the underlying function does it make? (Hoping someone will fire up their R interpreter and find out.)
Or, probably, dnorm is a probability distribution which includes a likeliness function, and a cumulative likeliness function, etc. I bet it doesn't work on arbitrary functions.
In the rectangle method, there is "let x = a + i * h + (h / F1::from(2)...)"
I didn't check, but I wonder if it is not better to do x = a + (i+0.5)*h... My reasoning is that if "i" is very big, then i * h can be much bigger than h/2 and thus h/2 may be eaten by numerical precision when added to i*h... And then you have to convince the optimizer to do it the way you want. But well, it's late...
I don't see any explicit SIMD in here. Is the rust compiler able to emit SIMD instructions automatically in cases like this? (I guess I could compile and disassemble to check... )
In my experience Rust is very good about using simd for loading and not great at using it automatically for math. This is from some experimentation at work and checking disassembly so ymmv
The Fn trait could be used, which prevents mutation, but allows a lot of useful closures. I should note, a motivated user could provide a junk function no matter what the type accepted is
I think the project should be called "NIntegrate".
BTW, that is not a serious suggestion; it is just that Wolfram Language (aka Mathematica) has both `Integrate` and `NIntegrate` for symbolic and numeric integration, respectively.
> Does the function oscillate over the region of integration? If it does, then make sure that the step size is chosen to be smaller than the wave length of the function.
Nyquist limit, but for numerical integration?
One way to think about is that these techniques work by integrating exactly the polynomial that interpolates the function where you're sampling it, so you need to resolve the features of the function to get good accuracy.
The Nyquist sampling theorem is of course proved by considering a Fourier transform, which is given by an integral, so the relation to integration in general should not be surprising.
Well in both cases it comes down to aliasing I think (high frequency wave presents as low frequency wave with too big step size)
It looks a bit sloppy to hardcode so many constants in a single file: `src/gauss_quadrature/legendre.rs`. Isn't it possible to generate them with the help of rust macros in the same way Julia uses metaprogramming?
Gaussian quadrature points are typically solved numerically. There's a good chance these ultimately came from a table.
Additionally, compile time floating-point evaluation is limited. When I looked around recently, I didn't see a rust equivalent of gcem; any kind of transcendental function evaluation (which finding Gaussian quadrature points absolutely would require) would not allow compile-time evaluation.
Support for float const fns was merged just a couple months ago and hasn't been officially announced yet.
Support for constant float operations was released in Rust 1.82! https://blog.rust-lang.org/2024/10/17/Rust-1.82.0.html
IIRC, that only supports elementary arithmetic operations. Useful but not general.
It's relatively straightforward to build transcendental functions out of the basic operations and the stdlib support will eventually get there, but rust's float story is still a work in progress. They're trying to do things more properly and document semantics better than C and C++ have.
I was under the impression that macros can execute arbitrary code, surely some FP would not be big problem. And if not macros then build.rs script certainly could do the trick.
build.rs can definitely execute arbitrary code, which means that a lot of places (including, IIRC crates.io) will forbid crates that rely on build.rs. I ended up refactoring my build.rs into a separate sub-application in finl_unicode that built data tables which are then checked into git and used pre-built. I include the sub-app in the source code so that anyone with access to the repo could continue development if I were to die tomorrow.
There are many crates with build.rs scripts on crates.io, because they host just the source code with the Cargo.{toml,lock}.
I ran into some issues with crates.io and my build.rs when I first released the crate, although it’s long enough ago, that I don’t remember exactly what the issue was. It might also have been that the build.rs script downloaded raw data files from unicode.org
crates.io doesn't care what your build.rs does because it doesn't try to compile your code, neither now or ever in the past. There would be no point in it trying to compile your code; there are lots of crates that are bindings to C libraries that crates.io's builders can't be expected to have, crates that target architectures which crates.io can't be expected to have builders for, etc.
this was mainly to use an Iteration free method in this paper: https://www.cfm.brown.edu/faculty/gk/APMA2560/Handouts/GL_qu...
this method is much faster and simpler.
Probably but that would slow down compilation a lot.
Exactly, it's not like the constants are gonna change.
You wouldn't have to recompile them every time. What if you didn't necessarily use macros but auto-generated it in a file that you keep separate from the other code at the bottom?
What I would do in these cases is to define the general computation function, but special-case it to return the hard-coded value for specific common inputs if it's being evaluated at compile time. Then add a test to verify both behaviors.
I was always amazed that R can do:
Can we do the same in this library?It seems like it is lacking the functionality R's integrate has for handling infinite boundaries, but I suppose you could implement that yourself on the outside.
For what it's worth,
prints Result: Ok(1.000000000053865)It does seem to be a usability hazard that the function being integrated is defined as a fn, rather than a Fn, as you can't pass closures that capture variables, requiring the weird dnorm definition
You will be completely blown away, then, from what Wolfram Language (aka Mathematica) can do. (When it comes to numerical integration.)
https://reference.wolfram.com/language/tutorial/NIntegrateOv...
for ]-inf, inf[ integrals, you can use Gauss Hermite method, just keep in mind to multiply your function with exp(x^2).
I got Result: 1.0000000183827922.How many evaluations of the underlying function does it make? (Hoping someone will fire up their R interpreter and find out.)
Or, probably, dnorm is a probability distribution which includes a likeliness function, and a cumulative likeliness function, etc. I bet it doesn't work on arbitrary functions.
R integrate is just a wrapper around quadpack. It works with arbitrary functions, but arguably dnorm is pretty well behaved.
In the rectangle method, there is "let x = a + i * h + (h / F1::from(2)...)"
I didn't check, but I wonder if it is not better to do x = a + (i+0.5)*h... My reasoning is that if "i" is very big, then i * h can be much bigger than h/2 and thus h/2 may be eaten by numerical precision when added to i*h... And then you have to convince the optimizer to do it the way you want. But well, it's late...
herbie recommends `fma(h, (i + 0.5), a)`, but also doesn't report any accuracy problems with the original either
yeah, fused mul-add is certainly better. Dunno how one epxresses that in rust though :-) Ahhh seems like there at least f64::mul_add() in stdlib :-)
I don't see any explicit SIMD in here. Is the rust compiler able to emit SIMD instructions automatically in cases like this? (I guess I could compile and disassemble to check... )
In my experience Rust is very good about using simd for loading and not great at using it automatically for math. This is from some experimentation at work and checking disassembly so ymmv
There are common library extensions for that.
This is a good joke, if anyone is wondering.
Fast, small, lightweight compared to what?
It's a minimal implementation (small, lightweight) of well-known and established algorithms (fast).
Those are adjectives not comparatives.
Thanks for showing this! It is very motivating to develop (and finish) my Raku numerical integration project.
Thanks! That’s awesome to hear—I’d love to see how your Raku numerical integration project turns out!
You can email me if you want to, I'll be happy to help.
Is there a technical reason to now allow closures as the integrand?
Mayve because they aren't guaranteed to be actual functions (in the mathematical sense) and could return random values
The Fn trait could be used, which prevents mutation, but allows a lot of useful closures. I should note, a motivated user could provide a junk function no matter what the type accepted is