[HN Gopher] Show HN: Rust library for numerical integration of r...
___________________________________________________________________
Show HN: Rust library for numerical integration of real-valued
functions
Author : mtantaoui
Score : 99 points
Date : 2024-11-19 12:41 UTC (10 hours ago)
(HTM) web link (github.com)
(TXT) w3m dump (github.com)
| mtantaoui wrote:
| Integrate is a fast, small, lightweight Rust library for
| performing numerical integration of real-valued functions. It is
| designed to integrate functions, providing a simple and efficient
| way to approximate definite integrals using various numerical
| methods.
|
| Integrate supports a variety of numerical integration techniques:
| - Newton-Cotes methods: - Rectangle Rule. -
| Trapezoidal Rule. - Simpson's Rule. - Newton's 3/8
| Rule.
|
| - Gauss quadrature methods: - Gauss-Legendre.
| - Gauss-Laguerre. - Gauss-Hermite. - Gauss-Chebyshev
| First Kind. - Gauss-Chebyshev Second Kind.
|
| - Adaptive Methods: - Adaptive Simpson's method
|
| - Romberg's method.
| antononcube wrote:
| 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.
| zppln wrote:
| Fast, small, lightweight compared to what?
| akkad33 wrote:
| Those are adjectives not comparatives.
| legobmw99 wrote:
| Is there a technical reason to now allow closures as the
| integrand?
| n_plus_1_acc wrote:
| Mayve because they aren't guaranteed to be actual functions
| (in the mathematical sense) and could return random values
| legobmw99 wrote:
| 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
| JanisErdmanis wrote:
| 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?
| ok123456 wrote:
| 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.
| AlotOfReading wrote:
| Support for float const fns was merged just a couple months
| ago and hasn't been officially announced yet.
| ok123456 wrote:
| IIRC, that only supports elementary arithmetic operations.
| Useful but not general.
| AlotOfReading wrote:
| 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.
| BD103 wrote:
| Support for constant float operations was released in Rust
| 1.82! https://blog.rust-
| lang.org/2024/10/17/Rust-1.82.0.html
| zokier wrote:
| 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.
| dhosek wrote:
| 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.
| n_plus_1_acc wrote:
| There are many crates with build.rs scripts on crates.io,
| because they host just the source code with the
| Cargo.{toml,lock}.
| dhosek wrote:
| 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
| two_handfuls wrote:
| Probably but that would slow down compilation a lot.
| simlevesque wrote:
| Exactly, it's not like the constants are gonna change.
| blharr wrote:
| 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?
| mtantaoui wrote:
| 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.
| wjholden wrote:
| I was always amazed that R can do: >
| integrate(dnorm, -Inf, +Inf) 1 with absolute error <
| 9.4e-05
|
| Can we do the same in this library?
| Buttons840 wrote:
| 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.
| thrasibule wrote:
| R integrate is just a wrapper around quadpack. It works with
| arbitrary functions, but arguably dnorm is pretty well
| behaved.
| legobmw99 wrote:
| 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
| antononcube wrote:
| 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...
| mtantaoui wrote:
| for ]-inf, inf[ integrals, you can use Gauss Hermite method,
| just keep in mind to multiply your function with exp(x^2).
| use integrate::{
| gauss_quadrature::hermite::gauss_hermite_rule, };
| use statrs::distribution::{Continuous, Normal}; fn
| dnorm(x: f64) -> f64 { Normal::new(0.0,
| 1.0).unwrap().pdf(x)* x.powi(2).exp() } fn
| main() { let n: usize = 170; let result
| = gauss_hermite_rule(dnorm, n); println!("Result:
| {:?}", result); }
|
| I got Result: 1.0000000183827922.
| antononcube wrote:
| Thanks for showing this! It is very motivating to develop (and
| finish) my Raku numerical integration project.
| mtantaoui wrote:
| 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.
| cozzyd wrote:
| 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... )
| jvanderbot wrote:
| 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.
| wiz21c wrote:
| 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...
| zokier wrote:
| herbie recommends `fma(h, (i + 0.5), a)`, but also doesn't
| report any accuracy problems with the original either
| zackangelo wrote:
| > 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?
| im3w1l wrote:
| Well in both cases it comes down to aliasing I think (high
| frequency wave presents as low frequency wave with too big step
| size)
___________________________________________________________________
(page generated 2024-11-19 23:00 UTC)