dextorious / numericalintegration.jl Goto Github PK

Basic numerical integration routines for presampled data.

License: Other

Julia 100.00%

numericalintegration.jl's Introduction

NumericalIntegration

This is a simple package to provide functionality for numerically integrating presampled data (meaning you can't choose arbitrary nodes). If you have the ability to evaluate your integrand at arbitrary points, please consider using better tools for the job (such as the excellent FastGaussQuadrature.jl).

Do note that while the code is trivial, it has not been extensively tested and does not focus on numerical precision. Issues, suggestions and pull requests are welcome.

Example usage

# setup some data
x = collect(-π/2 : π/1000 : π/2)
y = sin.(x)

# integrate using the default Trapezoidal method
integrate(x, y)

# integrate using a specific method
integrate(x, y, SimpsonEven())

# compute cumulative integral
Y = cumul_integrate(x, y)

# compute cumulative integral for each column of an array
z = [sin.(x) cos.(x) exp.(x/pi)]
Z = cumul_integrate(x, z)

# compute cumulative integral for each line of an array
zp = permutedims(z)
ZP = cumul_integrate(x, zp, dims=1)

# Multidimensional integration
Y = [i*j for i=x,j=y]
integrate((x,y), Y)

The currently available methods are:

Trapezoidal (default)
TrapezoidalEven
TrapezoidalFast
TrapezoidalEvenFast
SimpsonEven
SimpsonEvenFast
RombergEven

Currently cumulative integrals and multidimensional integrals are restricted to using Trapezoidal methods.

All methods containing "Even" in the name assume evenly spaced data. All methods containing "Fast" omit basic correctness checks and focus on performance. Consequently, the fast methods will segfault or produce incorrect results if you supply incorrect data (vectors of different lengths, etc.). RombergEven needs a power of 2 + 1 points (so 9, 17, 33, 65, 129, 257, 513, 1025...) evenly spaced for it to work. Useful when control over accuracy is needed.

numericalintegration.jl's People

Contributors

Stargazers

Watchers

numericalintegration.jl's Issues

Update TagBot

Due to the fact that GitHub disables actions with cron jobs on repositories that don't have activity for 60days, TagBot and CompatHelper might get disabled (see https://github.com/dextorious/NumericalIntegration.jl/actions?query=workflow%3ATagBot). For TagBot there is a solution which gets triggered by comments (see https://discourse.julialang.org/t/ann-required-updates-to-tagbot-yml/49249). I've made a PR for this in #35

Because of this issue the latest tag was not registered (and will have to be manually created I think).

initial value for cumulative integration

When I use culum_integrate, it appears that it integrates 1 index beyond the array.

using NumericalIntegration
a = LinRange(0,1,100)
print(integrate(a,a), ' ')
print(last(cumul_integrate(a,a)))

returns:

0.5 0.5000510152025299

Is there a way to define the initial value of the integration so that it returns an integration for 0 to 1 with the first index being the initial value?

Should I contribute here?

I need to integrate presampled 1d functions and this package seems to do exactly that job 😄
I now consider contributing. However there is a zoo of integration packages and before I do any work, I wanted to make sure that this is not abandoned in favor of some other package?

Multidimensional example is not working

Hi,

The following example (from the README file) is not working:

julia> using NumericalIntegration

julia> x = range(-π, π, length=50)
-3.141592653589793:0.1282282715750936:3.141592653589793

julia> typeof(x)
StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}

julia> y = sin.(x);

julia> ω = [xᵢ * yᵢ for xᵢ in x, yᵢ in y];

julia> res = integrate((x, y), ω)

julia> typeof(res)
Nothing

The current code also treats StepRangeLen in a different way than Array, although it seems not to be working as well (but for a different reason?):

julia> x = collect(range(-π, π, length=50));

julia> typeof(x)
Array{Float64,1}

julia> y = sin.(x);

julia> ω = [xᵢ * yᵢ for xᵢ in x, yᵢ in y];

julia> res = integrate((x, y), ω)
ERROR: knot-vectors must be sorted in increasing order
Stacktrace:
 [1] error(::String) at ./error.jl:33
 [2] check_gridded at /home/schneider/.julia/packages/Interpolations/qHlUr/src/gridded/gridded.jl:51 [inlined]
 [3] check_gridded at /home/schneider/.julia/packages/Interpolations/qHlUr/src/gridded/gridded.jl:52 [inlined]
 [4] Interpolations.GriddedInterpolation(::Type{Float64}, ::Tuple{Array{Float64,1},Array{Float64,1}}, ::Array{Float64,2}, ::Interpolations.Gridded{Interpolations.Linear}) at /home/schneider/.julia/packages/Interpolations/qHlUr/src/gridded/gridded.jl:32
 [5] interpolate(::Type{Float64}, ::Type{Float64}, ::Tuple{Array{Float64,1},Array{Float64,1}}, ::Array{Float64,2}, ::Interpolations.Gridded{Interpolations.Linear}) at /home/schneider/.julia/packages/Interpolations/qHlUr/src/gridded/gridded.jl:66
 [6] interpolate at /home/schneider/.julia/packages/Interpolations/qHlUr/src/gridded/gridded.jl:83 [inlined]
 [7] #LinearInterpolation#52 at /home/schneider/.julia/packages/Interpolations/qHlUr/src/convenience-constructors.jl:16 [inlined]
 [8] LinearInterpolation at /home/schneider/.julia/packages/Interpolations/qHlUr/src/convenience-constructors.jl:16 [inlined]
 [9] integrate(::Tuple{Array{Float64,1},Array{Float64,1}}, ::Array{Float64,2}, ::TrapezoidalFast) at /home/schneider/.julia/packages/NumericalIntegration/JrtZK/src/NumericalIntegration.jl:210
 [10] integrate at /home/schneider/.julia/packages/NumericalIntegration/JrtZK/src/NumericalIntegration.jl:203 [inlined]
 [11] integrate(::Tuple{Array{Float64,1},Array{Float64,1}}, ::Array{Float64,2}) at /home/schneider/.julia/packages/NumericalIntegration/JrtZK/src/NumericalIntegration.jl:339
 [12] top-level scope at REPL[41]:1

I'm using the latest version of the package:

julia> versioninfo()
Julia Version 1.4.1
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: AMD Ryzen 5 1600 Six-Core Processor
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-8.0.1 (ORCJIT, znver1)
Environment:
  JULIA_NUM_THREADS = 12

(@v1.4) pkg> st NumericalIntegration
Status `~/.julia/environments/v1.4/Project.toml`
  [e7bfaba1] NumericalIntegration v0.3.2

Multidimensional integration

Thank you for your work on the package.

Multidimensional integration would be a good feature to have. Using the dims keyword argument as in the base function sum.

Best,
Maxandre

`dims` ordering

Hi,

Is there a reason dims is ordered differently than the multi-dimensional arrays? I have to do:

julia> x=1:3
julia> y=1:5
julia> Z = randn(3, 5)
julia> integrate(x, Z; dims=2) # instead of `dims=1`
julia> integrate(y, Z; dims=1) # instead of `dims=2`

Thanks!

Slower than Trapz.jl

Comparison against Trapz.jl

I believe NumericalIntegration.jl package is much better structured and complete, though the integration routine could use some speed improvement, indeed the simple code in Trapz.jl is 13-16 times faster than the TrapezoidalFast() mode in this package.

If i used your package incorrently then i apologize in advance.
EDIT: old benchmark Old link

Updated benchmark after bugfix in Trapz: Here

Adding Romberg Method - Question about tests

I've added Romberg's method for higher precision integration (with the inconvenient that you need 2^n + 1 evenly sampled points). It's my first project in Julia and I am having problems adding support for Static Arrays. What's the recommended procedure? Should i send a PR with the test failing? Should I add an exception to the test so it doesn't test Romberg? Or should I not send a PR until the SVector test works as well?

Support Julia 1.0

Currently, import NumericalIntegration works on Julia 0.6 and Julia 0.7, but fails on Julia 1.0.

Could we add support for using NumericalIntegration.jl on Julia 1.0?

cc: @dextorious

Could you register the latest version?

Hi!

Thanks for making this package. The GitHub page shows a v0.3.0 tag, but there is no registered version of this package with that version in the General registry.

TagBot trigger issue

This issue is used to trigger TagBot; feel free to unsubscribe.

If you haven't already, you should update your TagBot.yml to include issue comment triggers.
Please see this post on Discourse for instructions and more details.

If you'd like for me to do this for you, comment TagBot fix on this issue.
I'll open a PR within a few hours, please be patient!

InexactError caused by milliseconds in DateTime

I am trying to compute the numerical integral from a time series. This fails due to a type conversion error that I do not really understand.

julia> using Dates, NumericalIntegration

julia> t = now() .- Millisecond.(90:-10:0)
10-element Vector{DateTime}:
 2021-05-29T16:28:27.071
 2021-05-29T16:28:27.081
 2021-05-29T16:28:27.091
 2021-05-29T16:28:27.101
 2021-05-29T16:28:27.111
 2021-05-29T16:28:27.121
 2021-05-29T16:28:27.131
 2021-05-29T16:28:27.141
 2021-05-29T16:28:27.151
 2021-05-29T16:28:27.161

julia> y = rand(10)
10-element Vector{Float64}:
 0.9599798922023242
 0.6541501943939838
 0.37632284117351933
 0.7371585093938842
 0.3863796441086791
 0.033734632856622815
 0.7643338520694836
 0.15143170117912086
 0.9771867088054458
 0.1608486283477646

julia> integrate(t, y)
ERROR: InexactError: Int64(16.14130086596308)
Stacktrace:
 [1] Int64
   @ ./float.jl:723 [inlined]
 [2] convert
   @ ./number.jl:7 [inlined]
 [3] Millisecond
   @ /buildworker/worker/package_linuxarmv7l/build/usr/share/julia/stdlib/v1.6/Dates/src/type                                                                                                 s.jl:34 [inlined]
 [4] *
   @ /buildworker/worker/package_linuxarmv7l/build/usr/share/julia/stdlib/v1.6/Dates/src/peri                                                                                                 ods.jl:94 [inlined]
 [5] integrate(x::Vector{DateTime}, y::Vector{Float64}, #unused#::TrapezoidalFast)
   @ NumericalIntegration ~/.julia/packages/NumericalIntegration/p29gp/src/NumericalIntegrati                                                                                                 on.jl:83
 [6] integrate
   @ ~/.julia/packages/NumericalIntegration/p29gp/src/NumericalIntegration.jl:58 [inlined]
 [7] #integrate#13
   @ ~/.julia/packages/NumericalIntegration/p29gp/src/NumericalIntegration.jl:333 [inlined]
 [8] integrate(x::Vector{DateTime}, y::Vector{Float64})
   @ NumericalIntegration ~/.julia/packages/NumericalIntegration/p29gp/src/NumericalIntegrati                                                                                                 on.jl:333
 [9] top-level scope
   @ REPL[5]:1

I don't think this is really a bug in NumericalIntegration.jl, since you can easily produce the same fault directly without the package.

julia> [t[i] - t[i-1] for i in 2:length(t)] .* [(y[i] + y[i-1])/2 for i in 2:length(y)]
ERROR: InexactError: Int64(8.07065043298154)
Stacktrace:
  [1] Int64
    @ ./float.jl:723 [inlined]
  [2] convert
    @ ./number.jl:7 [inlined]
  [3] Millisecond
    @ /buildworker/worker/package_linuxarmv7l/build/usr/share/julia/stdlib/v1.6/Dates/src/types.jl:34 [inlined]
  [4] *
    @ /buildworker/worker/package_linuxarmv7l/build/usr/share/julia/stdlib/v1.6/Dates/src/periods.jl:94 [inlined]
  [5] _broadcast_getindex_evalf
    @ ./broadcast.jl:648 [inlined]
  [6] _broadcast_getindex
    @ ./broadcast.jl:621 [inlined]
  [7] getindex
    @ ./broadcast.jl:575 [inlined]
  [8] macro expansion
    @ ./broadcast.jl:984 [inlined]
  [9] macro expansion
    @ ./simdloop.jl:77 [inlined]
 [10] copyto!
    @ ./broadcast.jl:983 [inlined]
 [11] copyto!
    @ ./broadcast.jl:936 [inlined]
 [12] copy
    @ ./broadcast.jl:908 [inlined]
 [13] materialize(bc::Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{1}, Nothing, typeof(*), Tuple{Vector{Millisecond}, Vector{Float64}}})
    @ Base.Broadcast ./broadcast.jl:883
 [14] top-level scope
    @ REPL[6]:1

As a workaround, you can use Dates.value to coerce a usable differential.

julia> Dates.value.([t[i] - t[i-1] for i in 2:length(t)]) .* [(y[i] + y[i-1])/2 for i in 2:length(y)]
9-element Vector{Float64}:
 8.07065043298154
 5.152365177837516
 5.567406752837018
 5.6176907675128165
 2.1005713848265097
 3.9903424246305317
 4.578827766243022
 5.643092049922833
 5.690176685766052

julia> sum(Dates.value.([t[i] - t[i-1] for i in 2:length(t)]) .* [(y[i] + y[i-1])/2 for i in 2:length(y)])
46.41112344255784

julia> integrate(Dates.value.(t), y)
46.41112344255784

Problems with Romberg and Simpson integration

Dear developer, hi to everybody! I am new to Julia, I am physicist and I am trying to understand if it makes sense to use julia for my next scientific projects. I was trying to perform some numerical integration and I found your package, so I decided to make some simple tests.
k = 2^4
x = LinRange(0, π, k+1)
y = sin.(x)
println(integrate(x, y, Trapezoidal()))
println(integrate(x, y, SimpsonEven()))
println(integrate(x, y, RombergEven()))

obtaining
1.9935703437723393
2.222549924725271
2.000005549979671
In order to be sure I am doing things correctly, I decided to cross-check the results performing the same calculation in Python
k = 2**4
x = np.linspace(0, np.pi, k+1)
y = np.sin(x)
print(np.trapz(y,x))
print(integrate.simps(y,x))
print(integrate.romb(y=y, dx=x[1]-x[0], show=False))
obtaining
1.9935703437723393
2.0000165910479355
1.9999999945872902
As you can see, the trapz integrations looks fine but there are huge discrepancies in the Simpson and Romberg results. I looked into your issues and I found this Romber.jl, so I decided to use it finding
println(romberg(x, y))
(1.9999999945872906, 5.555392379896773e-6)
This result is almost the same of the scipy simpson method!
Furthermore, I found that if try
k = 2^5
x = LinRange(0, π, k+1)
y = sin.(x)
println(integrate(x, y, RombergEven()))
I obtain
1.9999999945872906
Which is exactly the result of Romberg.jl and scipy! it looks that the result of RombergEven with 2^(k+1)+1 points is the same of scipy and romberg with 2^k+1 points.
Am I using your package in the wrong way or there is something else?
Thank you for your time!

Support for `Interpolations.jl` v0.14

Many state-of-the-art Julia Frameworks rely on Interpolations.jl v0.14, but NumericalIntegration.jl forces v0.13.
This results in many libraries aren't usabel together with NumericalIntegrations.jl.

Would be nice to add compat entry like:
Interpolations = "0.13, 0.14"

Thanks and best regards!

cumul_integrate gives wrong result?

julia> x = collect(1.0:5.0)
5-element Vector{Float64}:
 1.0
 2.0
 3.0
 4.0
 5.0

julia> y = x.^2
5-element Vector{Float64}:
  1.0
  4.0
  9.0
 16.0
 25.0

julia> cumul_integrate(x, y)
5-element Vector{Float64}:
  2.5
  5.0
 11.5
 24.0
 44.5

The result is different from
https://www.mathworks.com/help/matlab/ref/cumtrapz.html

Should we expect an array like [2.5, 9.0, 21.5, 42.0] from cumul_integrate?

improved implementation of Romberg integration

I posted an implementation on discourse (https://discourse.julialang.org/t/romberg-jl-simple-romberg-integration-in-julia/39313/8?u=stevengj) using Richardson.jl that is much more general (is not restricted to power-of-two sizes), seems to be more accurate than your current code, and is vastly simpler. I'm not sure if you would be interested in adopting it here?

My implementation could be sped up by inlining the sum calls and a few other tweaks, if necessary