I just started experimenting with the recent DynamicQuantities.jl package, which aims to provide the same functionality as Unitful.jl but with a statically-typed Quantity to improve performance. It looks like QuadGK currently doesn't support these DynamicQuantities Quantity types, but I'm not familiar enough with the codebase to quite understand why.

I'm interested in contributing by generating a PR for this, but would really appreciate a nudge in the right direction. It looks like the issue spawns from approximately here, where either a @generated cachedrule isn't being produced at compile-time or maybe the generic cachedrule below it isn't working.

Any ideas that might help me understand what I can do to help integrate this functionality?

Thanks,
Mike

Minimal reproduction (Julia v1.9.3):

julia> using DynamicQuantities  # [06fc5a27] DynamicQuantities v0.7.0

julia> using QuadGK  # [1fd47b50] QuadGK v2.8.2

julia> quadgk(t -> 5u"m/s"*t, 0u"s", 10u"s")
ERROR: MethodError: no method matching cachedrule(::Type{Quantity{Float64, Dimensions{DynamicQuantities.FixedRational{Int32, 25200}}}}, ::Int64)

Closest candidates are:
  cachedrule(::Union{Type{Complex{BigFloat}}, Type{BigFloat}}, ::Integer)
   @ QuadGK C:\Users\*\.julia\packages\QuadGK\XqIlh\src\gausskronrod.jl:253
  cachedrule(::Type{T}, ::Integer) where T<:Number
   @ QuadGK C:\Users\*\.julia\packages\QuadGK\XqIlh\src\gausskronrod.jl:264

I used to be able to do

julia> g(y) = quadgk(x -> x, y, 1.0)
g (generic function with 1 method)
julia> g(0.0)
(0.5, 0.0)
julia> g(1.0)
(0.0, 0.0)

Doing an integral where the upper and lower bound were identical gave the answer 0.0, as expected. After upgrading to v2.9.2, I get the error:

julia> g(1.0)
ERROR: DomainError with 1.0:
roundoff error detected near endpoint of the initial interval (1.0, 1.0)
Stacktrace:
  [1] throw_endpoint_error(x::Float64, a::Float64, b::Float64)
    @ QuadGK ~/.julia/packages/QuadGK/PWpum/src/adapt.jl:186
  [2] #check_endpoint_roundoff#48
    @ ~/.julia/packages/QuadGK/PWpum/src/adapt.jl:181 [inlined]
  [3] check_endpoint_roundoff
    @ ~/.julia/packages/QuadGK/PWpum/src/adapt.jl:179 [inlined]
  [4] #6
    @ ~/.julia/packages/QuadGK/PWpum/src/adapt.jl:15 [inlined]
  [5] ntuple
    @ ./ntuple.jl:48 [inlined]
  [6] do_quadgk(f::var"#5#6", s::Tuple{…}, n::Int64, atol::Nothing, rtol::Nothing, maxevals::Int64, nrm::typeof(LinearAlgebra.norm), segbuf::Nothing)
    @ QuadGK ~/.julia/packages/QuadGK/PWpum/src/adapt.jl:13
  [7] #51
    @ ~/.julia/packages/QuadGK/PWpum/src/adapt.jl:257 [inlined]
  [8] handle_infinities(workfunc::QuadGK.var"#51#52"{…}, f::var"#5#6", s::Tuple{…})
    @ QuadGK ~/.julia/packages/QuadGK/PWpum/src/adapt.jl:146
  [9] #quadgk#50
    @ QuadGK ~/.julia/packages/QuadGK/PWpum/src/adapt.jl:256 [inlined]
 [10] quadgk
    @ QuadGK ~/.julia/packages/QuadGK/PWpum/src/adapt.jl:254 [inlined]
 [11] g(y::Float64)
    @ Main ./REPL[2]:1
 [12] top-level scope
    @ REPL[3]:1

which seems to be caused by #91 .

Is there anyway to restore the functionality to accept identical lower and upper bounds? The above function is just an MWE, but it's not uncommon to have an integral inside a function where the lower bound is an argument and can vary to the upper bound.

As mentioned on discourse, the kronrod routine yields NaN for larger n, e.g.

kronrod(550,1e-11,40)

It looks like this is due to spurious underflow in the _kronrodjacobi function. It forms arrays s and t by calculations that involve repeatedly multiplying by the b matrix, corresponding to the squared off-diagonal elements of the Jacobi matrix. For the canonical interval (-1,1) with a constant weight function, the entries of the b matrix asymptote to 0.25, so this effectively divides by 4 over and over and eventually underflows.

The underflow is spurious because the final Kronrod–Jacobi matrix only depends on ratios of s and t entries, so we are free to rescale them by any overall constant factor that we want. This fact is actually mentioned in a somewhat obscure way by Laurie's paper:

Actually, as recommended by Reichel [23], the coefficients were scaled to the interval (−2, 2) before performing the computation, and the resulting points and weights (obtained by Gautschi’s routine gauss from [10]) were scaled back to (−1, 1). On a machine with a small exponent range, this subterfuge is necessary to avoid underflow of the modified moments.

Rescaling the interval to (-2,2) works because it multiplies the integrals by 2, and hence the squared Jacobi matrix entries get multiplied by 4, causing b to asymptote to 1 and removing the exponential shrinkage.

We could do the same trick, but it's not clear how to generalize this to arbitrary Jacobi matrices. Instead, it should suffice to normalize the s and t matrices (e.g. scaling by their infinity norm) after each iteration. This should only slow things down by a small constant factor (1–2% in my current experiments).

DataStructures 0.18 is out, can this package be revised to support it? I am having conflicts between this package and the latest version of DiffEqBase due to this dependency version.

Hi, I want to cite this package in a book and can't figure out the author field. So far I'm using

@misc{QuadGK,
title="{QuadGK.jl}",
year=2021,
note="Julia Package",
howpublished="https://github.com/JuliaMath/QuadGK.jl"
}

which works well with the publisher's LaTeX style. So .... author = ?
Thanks,

-- Tim

This seems wrong:

julia> quadgk_count(x -> 1, 0, 1, order=1)
(2.4444444444019076, 0.5555555555356156, 23333345)

This is possibly related to #30, but not sure. I am integrating a function that is basically the Gamma PDF (very spiked near zero, but tapers off as x -> infinity). The function does not return NaN for any value (I fix gpdf(0)=0, or else it would give NaN here).

quadgk(x->gpdf(x,l1,d1,s1),0,Inf)
ERROR: DomainError with 0.9999999999999964:
integrand produced NaN in the interval (0.9999999999999929, 1.0)
Stacktrace:
 [1] evalrule(::getfield(QuadGK, Symbol("##12#18")){getfield(Main, Symbol("##41#42")),Float64}, ::Float64, ::Float64, ::Array{Float64,1}, ::Array{Float64,1}, ::Array{Float64,1}, ::typeof(LinearAlgebra.norm)) at C:\Users\tbeason\.julia\packages\QuadGK\v6Zrs\src\evalrule.jl:41
 [2] adapt(::Function, ::Array{QuadGK.Segment{Float64,Float64,Float64},1}, ::Float64, ::Float64, ::Int64, ::Array{Float64,1}, ::Array{Float64,1}, ::Array{Float64,1}, ::Int64, ::Float64, ::Float64, ::Int64, ::typeof(LinearAlgebra.norm)) at C:\Users\tbeason\.julia\packages\QuadGK\v6Zrs\src\adapt.jl:39
 [3] do_quadgk(::getfield(QuadGK, Symbol("##12#18")){getfield(Main, Symbol("##41#42")),Float64}, ::Tuple{Float64,Float64}, ::Int64, ::Nothing, ::Nothing, ::Int64, ::typeof(LinearAlgebra.norm)) at C:\Users\tbeason\.julia\packages\QuadGK\v6Zrs\src\adapt.jl:28
 [4] handle_infinities(::getfield(Main, Symbol("##41#42")), ::Tuple{Float64,Float64}, ::Int64, ::Nothing, ::Nothing, ::Int64, ::Function) at C:\Users\tbeason\.julia\packages\QuadGK\v6Zrs\src\adapt.jl:83
 [5] #quadgk#21(::Nothing, ::Nothing, ::Int64, ::Int64, ::Function, ::typeof(quadgk), ::Function, ::Float64, ::Float64) at C:\Users\tbeason\.julia\packages\QuadGK\v6Zrs\src\adapt.jl:155
 [6] #quadgk#20 at C:\Users\tbeason\.julia\packages\QuadGK\v6Zrs\src\adapt.jl:155 [inlined]
 [7] quadgk(::Function, ::Int64, ::Float64) at C:\Users\tbeason\.julia\packages\QuadGK\v6Zrs\src\adapt.jl:151
 [8] top-level scope at REPL[109]:1

The error indicates that the NaN arises in the handle_infinities function, which I guess it where you transform the integration domain from (0,Inf) to something usable. If I use a large number instead of Inf, I do not get the error.

MWE

using SpecialFunctions, QuadGK
l1 = exp(-1.2604534253)
d1 = -0.916695708
s1 = 8.428294026

function gpdf(t,lambda,delta,sigma)
    t == 0 && return 0.0
    d2 = delta^(-2)
    constantpiece = lambda*abs(delta/sigma)/gamma(d2)
    LT = lambda*t
    varpiece = LT^(1/(delta*sigma)-1) * exp(-LT^(delta/sigma))
    val = constantpiece*varpiece
    return val
end

quadgk(x->gpdf(x,l1,d1,s1),0,Inf)
quadgk(x->gpdf(x,l1,d1,s1),0,1e30)

Hi,

Is there another way than using Sympy from Python?

If i do need to use it, as my integral is within a function and part of a larger numerical problem, should i then reconvert all my variables into Julia types (I mean versus symbolic - type Sym- that I had to convert them into to use integrate)?

Many thanks.

I think there's something wrong with the way quadgk handles this integral, the error estimate is almost zero but the answer is substantially inaccurate. Just replacing 0.499 with 0.4999 fixes it:

julia> check(f,a,b) = let (x1,e1) = QuadGK.quadgk(f,a,b); (x2,e2) = QuadGK.quadgk(f,big(a),big(b)); (Float64(abs(x1-x2)), e1); end
check (generic function with 1 method)

julia> check(x -> exp(abs(x-0.499)), 0, 1)
(1.0000000833402914e-6, 2.220446049250313e-16)

julia> check(x -> exp(abs(x-0.4999)), 0, 1)
(7.284892445276705e-17, 1.1102230246251565e-16)

It's an easy integral, so zero error estimate can be expected, but there could be something going wrong with the integral value it actually returns.

When I run the integral of f(x)="0.21702437*(0.6227007 + 1.4335294/(x + 1.6759317))/(x + 1.52924537 - 0.2316349/x)" from x=0 to x=1, by running this code in Python:

jl.eval("""

my_tuple = quadgk(x -> 0.21702437*(0.6227007 + 1.4335294/(x + 1.6759317))/(x + 1.52924537 - 0.2316349/x), 0, 1, rtol=0.01)
println(my_tuple)
""")

I consistently get a value of "0.17259607425281273" for my_tuple[1].
I would expect to get a value of infinity or undefined.

Interestingly, when I use PySR python package's eval_tree_array and I pass in a tree whose definition is:

x1 = Node(;feature=1)
tree = 0.21702437*(0.6227007 + 1.4335294/(x1 + 1.6759317))/(x1 + 1.52924537 - 0.2316349/x1)

tree = Node{Float32}(tree)

my_tuple = quadgk(x -> (eval_tree_array(tree, reshape([Float32(x)], 1, 1), options))[1][1], 0, 1, rtol=0.01)

I consistently get a value of "0.14588867513106551" for my_tuple[1].
Perhaps the difference between 0.17259 and 0.14588 is because I converted the tree to Float32 (instead of Float64).

If you open up Desmos.com and plot this function f(x) over x=0 to x=1, you will find that the function goes to negative infinity as x approaches 0.16 from the to-the-right direction, and the function goes to positive nfinity as x approaches 0.16 from the to-the-left direction.

And I'm pretty sure that in reality, the integral of this function is undefined. Although perhaps, the negative infinity cancels out the positive infinity, which means the value of the definite integral should be finite?

Thanks for all the help, in advance.

rulekey(::Type{BigFloat}, n) = (BigFloat, precision(BigFloat), n)
rulekey(T,n) = (T,n)

The function returns either a triplet or a pair, this would make it easier to reason about if it returned always the same structure and types

Hi, I have just tried to integrate a Laplace distribution as

using QuadGK, Distributions
pc = Laplace(0, 1)
quadgk(x -> pdf(pc, x), -1e8, 1e8, rtol=1e-8)
(0.0, 0.0)

I obtain zero, which is obviously wrong as it should be one. The manual says that I should add an explicit point, when the point is non-smooth, which I do

quadgk(x -> pdf(pc, x), -1e8, 0, 1e8, rtol=1e-8)
(0.0, 0.0)

But it does not help either. I think that the method just omit to sample around the zero sufficiently.

QuadGK is holding my DataStructures to v0.15.

What is the reason for these lines in the Project.toml?

[compat]
DataStructures = "0.11, 0.12, 0.13, 0.14, 0.15"

Hi there, I would like to ask a question regarding quadrature points and weights in the intervals [0, +Inf] and [-Inf, +Inf]. At the moment, it seems that I could not get quadrature points and weights in these intervals using the command gauss(N,a,b). However, it is possible to find quadgk(x -> exp(-x^2), -Inf, +Inf) or quadgk(x -> exp(-x^2), 0, +Inf). So, my question here is that will the command gauss(N,a,b) support the intervals [0, Inf] and [-Inf, +Inf] in a near future? Or I have to perform a change of variable to scale these intervals into finite intervals to get the quadrature points and weights?

JuliaError: Exception 'UndefVarError: integral not defined' occurred while calling julia code:

When I ran this code:
"
using QuadGK

try
integral, err = quadgk(x -> exp(1/(x-0.5)), 0, 1, rtol=1e-8)
catch
print("Error")
end
print(integral)
print(err)
"

I am generating a bunch of expression trees somewhat randomly, where each expression tree represents a particular function of x.
For example, one expression tree could be:
*
/
cos 5
|
x0
which represents the function "cos(x0)*5". I can call this function tree_as_func(x).

I need to run "quadgk(x -> tree_as_func(x), 0, 1, rtol=1e-8)" where "tree_as_func" is the function of x that the current expression tree represents.
I do not know if "tree_as_func(x)" blows up to infinity on the domain x=[0, 1) or if its integral over x=[0, 1) is infinity.

As a result, I need to be able to catch errors. But this try-catch block fails to catch the error.
How can I do this?

I ran

@code_warntype quadgk( x->x, 0., 1. )

today and it says, that the return type doesn't get inferred. Is it me or is it a bug?

Hi there,

I would like to ask you a question regarding quadrature rules for the weight function w(x)=1/(1+x^2). I hope this is fine to post my question here. My question is that I would like to compute the integral I = \int_{a}^{b} w(x)f(x)dx using quadrature rules, where w(x) is a weight function and f(x) can be approximated by a polynomial of degree 2N−1 or less on [a, b]. This is because it is costly to compute f(x). Therefore, I am looking for a quadrature rule for the weight function w(x)=1/(1+x^2). I had a look at https://en.wikipedia.org/wiki/Gaussian_quadrature. However, I could not find a quadrature rule for the weight function that I need. But I found that your package can compute the integral I for arbitrary weight functions. So, do you think that it is fine that I use the function "QuadGK.gauss" to obtain the quadrature points x[i] and weights w[i] and can compute the integral I from sum(w .* f.(x)). Or do you have any suggestions for me? Thank you very much in advance!

When you follow the deprecation advice for using this package on Julia 0.6 and beyond,

julia> quadgk()
ERROR: quadgk() has been moved to the package QuadGK.jl.
Run Pkg.add("QuadGK") to install QuadGK on Julia v0.6 and later, and then run `using QuadGK`.

It does not actually work due to a clash with the deprecation warning itself.

julia> using QuadGK
julia> quadgk()
WARNING: both QuadGK and Base export "quadgk"; uses of it in module Main must be qualified
ERROR: UndefVarError: quadgk not defined

I believe the above was working on MASTER sometime over Winter 2016 / Spring 2017, as I had added using QuadGK to all my codes to make the deprecation warning go away!

Instead you must import the function(s) you want explicitly,

julia> import QuadGK.quadgk
julia> quadgk
quadgk (generic function with 3 methods)

Or reference the QuadGK package every time you call the function,

julia> using QuadGK
julia> QuadGK.quadgk
quadgk (generic function with 3 methods)

See also discussion at:
JuliaLang/julia#19741

Can the QuadGK functions be forcibly exported into the namespace when the packaged is using'd ? That would seem to be the expected behaviour by the user.

Hi,

wanted to know if this is the expected behavior. When doing

function f!(b, x)
   b[1] = exp(-x^2)
   b[2] = exp(-x^2)
   return nothing 
end

buf = zeros(Float64, 2); quadgk!(f!, buf, 0.0, Inf)

I get ERROR: MethodError: objects of type QuadGK.InplaceIntegrand{typeof(f!),Array{Float64,1},Array{Float64,1}} are not callable. For finite upper bounds though everything seems to work. In my use case however, I need to compute integrals for which the upper bound is Inf, and I do not a priori know if and when I can truncate the integral for some large upper bound.

While this package allows evaluation of integrands which return arrays, it's not efficient for the integrand to be an allocating f(t). Instead, quadgk should support a non-allocating form f!(out,t) for this case, for this version of the function do minimal internal allocation by caching necessary arrays and doing everything in-place.

Reported by @sosock in PainterQubits/Unitful.jl#597 (comment):

julia> using Unitful, QuadGK

julia> quadgk(x->exp(-x/(1.0u"m")), 0.0u"m", Inf*u"m")
ERROR: DomainError with Inf m:
integrand produced NaN m in the interval (0.0 m, Inf m)
Stacktrace:
[...]

In the documentation for kronrod() I see references to both wg and gw. Also there are two entries for kronrod() on that page, the second near the bottom.

Following is my code to solve two-dimensional integral in Julia.

using Statistics,Distributions,Distributed,StatsFuns,Roots,QuadGK
p=3;n0=3;n1=1;n2=29;d=0.25;
inARL=200;n=n1+n2;w1=chisqinvcdf(p,1-(n0-n1)/n2);inpa1=chisqcdf(p,w1);

function inpa(h2)
function fvv(v)
return(nchisqcdf(p,(n1/n2)v,(n*h2)/n2)*chisqpdf(p,v));end;
return(1/(1-inpa1-quadgk(fvv,w1,Inf,rtol=1e-10)[1])-inARL);end;

h=fzero(inpa,0.25,50);
outpa1=nchisqcdf(p,n1*d^2,w1);

function fz1(z)
function fv1(v)
lam=(n1/n2)*((z+sqrt(n1)*d)^2+v)+(2*n*sqrt(n1)*d*z)/(n2)+((n*d)^2)/(n2);
return(nchisqcdf(p,lam,(n*h)/n2)*chisqpdf(p-1,v));end;
return(normpdf(z)*quadgk(fv1,0,Inf,rtol=1e-10)[1]);end;
int1=quadgk(fz1,-Inf,-sqrt(w1)-sqrt(n1)*d,rtol=1e-10)[1];# problem

I have a function if I have a function f(x) which is the integral another function (using QuadGK) and x is the upper limit of integration. If I want to take an automatic derivative of that function using ForwardDiff, it results in a stack overflow. Here are instructions to reproduce the error

using QuadGK
using ForwardDiff
D(f::Function, x::Number) = ForwardDiff.derivative(f, x)

f(x) = quadgk(y -> y^2, 0, x)
D(x -> quadgk(y -> y^2, 0, x), 1)

ERROR: StackOverflowError:
Stacktrace:
 [1] #quadgk#15(::Array{Any,1}, ::Function, ::Function, ::ForwardDiff.Dual{ForwardDiff.Tag{##3#5,Int64},Float64,1}, ::ForwardDiff.Dual{ForwardDiff.Tag{##3#5,Int64},Float64,1}) at /Users/mason/.julia/v0.6/QuadGK/src/QuadGK.jl:241 (repeats 65479 times)
 [2] derivative(::##3#5, ::Int64) at /Users/mason/.julia/v0.6/ForwardDiff/src/derivative.jl:14
 [3] D(::Function, ::Int64) at ./REPL[4]:1

quadgk is not compatible with AD (well, with ForwardDiff.jl) because of cachedrule. I am trying to solve an optimization problem using JuMP where some nonlinear constraints involve integrals. I can compute them using other methods, but it would be nice if QuadGK "just worked"

Here is the rather uninformative error

ERROR: StackOverflowError:
Stacktrace:
 [1] cachedrule(#unused#::Type{ForwardDiff.Dual{ForwardDiff.Tag{JuMP.var"#137#139"{typeof(integratesegmentgk)}, Float64}, Float64, 5}}, n::Int64) (repeats 79984 times)
   @ QuadGK C:\Users\beasont\.julia\packages\QuadGK\ENhXl\src\gausskronrod.jl:253

The reference to integratesegmentgk is my user-defined function in the constraint. It simply calls quadgk, but for completeness:

function integratesegmentgk(mL,mU,a,b,c)
	hfun(x) = exp(-a -b*x -c*x^2)
	res,tol = quadgk(hfun,mL,mU)
	return res
end

I am using FastGaussQuadrature.jl in the meantime.

@stevengj Do you have any input on what should be done here before this is tagged and registered?

Some problem with the DOCUMENTER_KEY ssh key according to the logs.

Initialized empty Git repository in /tmp/jl_B4CaCu/.git/
[email protected]: Permission denied (publickey).
fatal: Could not read from remote repository.

Please make sure you have the correct access rights
and the repository exists.
┌ Error: Git failed to fetch [email protected]:JuliaMath/QuadGK.jl.git
│ This can be caused by a DOCUMENTER_KEY variable that is not correctly set up.
│ Make sure that the environment variable is properly set up as a Base64-encoded string
│ of the SSH private key. You may need to re-generate the keys with DocumenterTools.
└ @ Documenter ~/.julia/packages/Documenter/6vUwN/src/Documenter.jl:552

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!

I think I might have discovered a recent regression in QuadGK when using the quadgk function with ranges specified in Unitful types. I tried to dig into the source by following the stack trace, and it seems like it might be related to the recent changes to handle_infinities, but I'm not familiar enough with this package or Julia's underlying type mechanics to nail down the issue at the moment.

Steps to Reproduce

I was able to reproduce the issue using the following minimal code snippet. This snippet should integrate a constant 1 Joule value over the domain [0,10] seconds, which should result in a value of 10 Joule-seconds.

using QuadGK, Unitful
f(t) = 1.0u"J"
quadgk(f, 0.0u"s", 10.0u"s")

Test Configuration 1 (v2.6.0, Error)

System information for this configuration:

versioninfo()
# Julia Version 1.8.5
# Commit 17cfb8e65e (2023-01-08 06:45 UTC)
# Platform Info:
#   OS: Windows (x86_64-w64-mingw32)
#   CPU: 32 × AMD Ryzen 9 3950X 16-Core Processor
#   WORD_SIZE: 64
#   LIBM: libopenlibm
#   LLVM: libLLVM-13.0.1 (ORCJIT, znver2)
#   Threads: 1 on 32 virtual cores

import Pkg
Pkg.status(["Unitful", "QuadGK"])
# Status `C:\Users\mikei\.julia\environments\v1.8\Project.toml`
#   [1fd47b50] QuadGK v2.6.0
#   [1986cc42] Unitful v1.12.2

Running the code snippet above produces an error

using QuadGK, Unitful

f(t) = 1.0u"J"
# f (generic function with 1 method)

quadgk(f, 0.0u"s", 10.0u"s")
#=
ERROR: MethodError: no method matching Quantity{Float64, 𝐓 , Unitful.FreeUnits{(s,), 𝐓 , nothing}}(::Int64)
The applicable method may be too new: running in world age 32472, while current world is 32475.
Closest candidates are:
  Quantity{T, D, U}(::Number) where {T, D, U} at C:\Users\mikei\.julia\packages\Unitful\fbiNW\src\types.jl:151 (method too new to be called from this world context.)
  (::Type{T})(::T) where T<:Number at boot.jl:772
  (::Type{T})(::AbstractChar) where T<:Union{AbstractChar, Number} at char.jl:50
  ...
Stacktrace:
  [1] convert(#unused#::Type{Quantity{Float64, 𝐓, Unitful.FreeUnits{(s,), 𝐓, nothing}}}  , x::Int64)
    @ Base .\number.jl:7
  [2] one(#unused#::Type{Quantity{Float64, 𝐓, Unitful.FreeUnits{(s,), 𝐓, nothing}}}  )
    @ Base .\number.jl:334
  [3] #s14#1
    @ C:\Users\mikei\.julia\packages\QuadGK\kf0xA\src\gausskronrod.jl:260 [inlined]
  [4] var"#s14#1"(T::Any, ::Any, #unused#::Any, n::Any)
    @ QuadGK .\none:0
  [5] (::Core.GeneratedFunctionStub)(::Any, ::Vararg{Any})
    @ Core .\boot.jl:582
  [6] do_quadgk(f::typeof(f), s::Tuple{Quantity{Float64, 𝐓, Unitful.FreeUnits{(s,), 𝐓, nothing}}, Quantity{Float64, 𝐓
, Unitful.FreeUnits{(s,), 𝐓, nothing}}} , n::Int64, atol::Nothing, rtol::Nothing, maxevals::Int64, nrm::typeof(LinearAlgebra.norm), segbuf::Nothing)
    @ QuadGK C:\Users\mikei\.julia\packages\QuadGK\kf0xA\src\adapt.jl:7
  [7] #28
    @ C:\Users\mikei\.julia\packages\QuadGK\kf0xA\src\adapt.jl:186 [inlined]
  [8] handle_infinities(workfunc::QuadGK.var"#28#29"{Nothing, Nothing, Int64, Int64, typeof(LinearAlgebra.norm), Nothing}, f::typeof(f), s::Tuple{Quantity{Float64, 𝐓, Unitful.FreeUnits{(s,), 𝐓, nothing}}, Quantity{Float64, 𝐓, Unitful.Fre
eUnits{(s,), 𝐓, nothing}}} )
    @ QuadGK C:\Users\mikei\.julia\packages\QuadGK\kf0xA\src\adapt.jl:115
  [9] quadgk(::Function, ::Quantity{Float64, 𝐓, Unitful.FreeUnits{(s,), 𝐓, nothing}}  , ::Vararg{Quantity{Float64, 𝐓, U
nitful.FreeUnits{(s,), 𝐓, nothing}}} ; atol::Nothing, rtol::Nothing, maxevals::Int64, order::Int64, norm::Function, segbuf::Nothing)
    @ QuadGK C:\Users\mikei\.julia\packages\QuadGK\kf0xA\src\adapt.jl:185
 [10] quadgk(::Function, ::Quantity{Float64, 𝐓, Unitful.FreeUnits{(s,), 𝐓, nothing}}  , ::Vararg{Quantity{Float64, 𝐓, U
nitful.FreeUnits{(s,), 𝐓, nothing}}} )
    @ QuadGK C:\Users\mikei\.julia\packages\QuadGK\kf0xA\src\adapt.jl:183
 [11] top-level scope
    @ REPL[3]:1
=#

Test Configuration 2 (v.2.5.0, Working)

For the next configuration, I simply pinned QuadGK to version 2.5.0 in Julia's Pkg REPL and re-ran the test. System information for this configuration:

versioninfo()
# Julia Version 1.8.5
# Commit 17cfb8e65e (2023-01-08 06:45 UTC)
# Platform Info:
#   OS: Windows (x86_64-w64-mingw32)
#   CPU: 32 × AMD Ryzen 9 3950X 16-Core Processor
#   WORD_SIZE: 64
#   LIBM: libopenlibm
#   LLVM: libLLVM-13.0.1 (ORCJIT, znver2)
#   Threads: 1 on 32 virtual cores

import Pkg
Pkg.status(["Unitful", "QuadGK"])
# ⌃ [1fd47b50] QuadGK v2.5.0 ⚲
#   [1986cc42] Unitful v1.12.2
# Info Packages marked with ⌃ have new versions available and may be upgradable.

In this configuration, the test runs as expected.

using QuadGK, Unitful

f(t) = 1.0u"J"
# f (generic function with 1 method)

quadgk(f, 0.0u"s", 10.0u"s")
# (10.0 J s, 0.0 J s)

The issue is that segbufs don't always work with infinite limits, such as when the limits have units.
This is a MWE

using QuadGK, Unitful

Tu = typeof(1.0u"m")
segbuf = QuadGK.alloc_segbuf(Tu, Tu, Tu)

quadgk(x -> 1/(1+oneunit(x)/x), 0.0u"m", 1.0u"m", segbuf=segbuf) # OK (0.3068528194400547 m, 1.9931833961095435e-11 m)

quadgk(x -> 1/(1+oneunit(x)/x), 0.0u"m", Inf*u"m", segbuf=segbuf) # broken

which gives the error

ERROR: DimensionError: m and 0.0 are not dimensionally compatible.
Stacktrace:
  [1] convert(#unused#::Type{Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}}, x::Float64)
    @ Unitful ~/.julia/packages/Unitful/jiUeO/src/conversion.jl:106
  [2] QuadGK.Segment{Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}, Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}, Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}}(a::Float64, b::Float64, I::Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}, E::Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}})
    @ QuadGK ~/.julia/dev/QuadGK/src/evalrule.jl:3
  [3] convert(#unused#::Type{QuadGK.Segment{Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}, Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}, Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}}}, s::QuadGK.Segment{Float64, Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}, Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}})
    @ QuadGK ~/.julia/dev/QuadGK/src/evalrule.jl:10
  [4] setindex!(A::Vector{QuadGK.Segment{Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}, Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}, Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}}}, x::QuadGK.Segment{Float64, Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}, Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}}, i1::Int64)
    @ Base ./array.jl:969
  [5] macro expansion
    @ ./broadcast.jl:974 [inlined]
  [6] macro expansion
    @ ./simdloop.jl:77 [inlined]
  [7] copyto!
    @ ./broadcast.jl:973 [inlined]
  [8] copyto!
    @ ./broadcast.jl:926 [inlined]
  [9] materialize!
    @ ./broadcast.jl:884 [inlined]
 [10] materialize!
    @ ./broadcast.jl:881 [inlined]
 [11] do_quadgk(f::QuadGK.var"#18#27"{var"#11#12", Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}, Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}}, s::Tuple{Float64, Float64}, n::Int64, atol::Nothing, rtol::Nothing, maxevals::Int64, nrm::typeof(LinearAlgebra.norm), segbuf::Vector{QuadGK.Segment{Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}, Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}, Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}}})
    @ QuadGK ~/.julia/dev/QuadGK/src/adapt.jl:39
 [12] #50
    @ ~/.julia/dev/QuadGK/src/adapt.jl:235 [inlined]
 [13] handle_infinities(workfunc::QuadGK.var"#50#51"{Nothing, Nothing, Int64, Int64, typeof(LinearAlgebra.norm), Vector{QuadGK.Segment{Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}, Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}, Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}}}}, f::var"#11#12", s::Tuple{Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}, Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}})
    @ QuadGK ~/.julia/dev/QuadGK/src/adapt.jl:126
 [14] quadgk(::Function, ::Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}, ::Vararg{Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}}; atol::Nothing, rtol::Nothing, maxevals::Int64, order::Int64, norm::Function, segbuf::Vector{QuadGK.Segment{Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}, Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}, Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}}})
    @ QuadGK ~/.julia/dev/QuadGK/src/adapt.jl:234
 [15] top-level scope
    @ REPL[11]:1

An issue was caught in SciML/DataInterpolations.jl#178 where tests started failing with [email protected]

I have made a MWE to illustrate the issue (using a simple LinearInterpolation from DataInterpolations):

u = collect(0:10)
t = collect(0:10)
A = LinearInterpolation(u,t)
qint, err = quadgk(A, 0.0, 10.0)

It is a straight line with slope 1. So the integral should be 50.0.

with [email protected]

julia> qint, err = quadgk(A, 0.0, 10.0)
(50.0, 0.0)

with [email protected]

julia> qint, err = quadgk(A, 0.0, 10.0)
(51.76911234248546, 4.274582013885939)

Dear QuadGK developers:
I'm attempting to calculate multiple 1D integrals in parallel using GPU via CUDA.jl. However currently QuadGK does not seem to be GPU-compatible (with CUDA.jl). A simple example below generates errors ( see the attached text file ). Is it possible, perhaps for the 1D integrals with limited options, to become compatible with CUDA.jl? Thanks.

using CUDA
using QuadGK

Z = Array([z for z in 1:1:4]);
Zcu = CuArray([z for z in 1:1:4]);

function f_int(z_ob) 

    iii(z) =  z_ob - z  
    
    return quadgk(z -> iii(z), 0, 1, rtol=1e-4)[1]
    
end

f_int.(Z)  # This works 

f_int.(Zcu) # This does not work

QUADgk error with CUDA.txt

Could you please tag a new release for quadgk!?

Two minimal working examples:

using QuadGK, Unitful
f(x) = x^2u"m"
quadgk(f, 0, 1)

and

using QuadGK, Unitful
f(x) = x^2
quadgk(f, 0u"m", 1u"m")

both crash Julia for me.

Seems to be directly related to #13. In hindsight, this probably could have been a comment in that thread. Should I remove this issue and move its contents there?

Ref: https://discourse.julialang.org/t/package-compatibility-caps/15301

This is probably related to #38, I am trying to integrate with quadgk a function from 1 to infinity which should be converging, and encounter this error (julia 1.8.5 and QuadGK v2.8.2).

using QuadGK
quadgk(x->1/x^1.1,1,Inf)
DomainError with 0.9999999999999964:
integrand produced NaN in the interval (0.9999999999999929, 1.0)

Stacktrace:
  [1] evalrule(f::QuadGK.var"#14#23"{var"#1#2", Float64, Float64}, a::Float64, b::Float64, x::Vector{Float64}, w::Vector{Float64}, gw::Vector{Float64}, nrm::typeof(LinearAlgebra.norm))
    @ QuadGK ~/.julia/packages/QuadGK/XqIlh/src/evalrule.jl:37
  [2] adapt(f::QuadGK.var"#14#23"{var"#1#2", Float64, Float64}, segs::Vector{QuadGK.Segment{Float64, Float64, Float64}}, I::Float64, E::Float64, numevals::Int64, x::Vector{Float64}, w::Vector{Float64}, gw::Vector{Float64}, n::Int64, atol::Float64, rtol::Float64, maxevals::Int64, nrm::typeof(LinearAlgebra.norm))
    @ QuadGK ~/.julia/packages/QuadGK/XqIlh/src/adapt.jl:48
  [3] do_quadgk(f::QuadGK.var"#14#23"{var"#1#2", Float64, Float64}, s::Tuple{Float64, Float64}, n::Int64, atol::Nothing, rtol::Nothing, maxevals::Int64, nrm::typeof(LinearAlgebra.norm), segbuf::Nothing)
    @ QuadGK ~/.julia/packages/QuadGK/XqIlh/src/adapt.jl:37
  [4] #46
    @ ~/.julia/packages/QuadGK/XqIlh/src/adapt.jl:219 [inlined]
  [5] handle_infinities(workfunc::QuadGK.var"#46#47"{Nothing, Nothing, Int64, Int64, typeof(LinearAlgebra.norm), Nothing}, f::var"#1#2", s::Tuple{Float64, Float64})
    @ QuadGK ~/.julia/packages/QuadGK/XqIlh/src/adapt.jl:110
  [6] quadgk(::Function, ::Float64, ::Vararg{Float64}; atol::Nothing, rtol::Nothing, maxevals::Int64, order::Int64, norm::Function, segbuf::Nothing)
    @ QuadGK ~/.julia/packages/QuadGK/XqIlh/src/adapt.jl:218
  [7] quadgk(::Function, ::Float64, ::Vararg{Float64})
    @ QuadGK ~/.julia/packages/QuadGK/XqIlh/src/adapt.jl:216
  [8] quadgk(::Function, ::Int64, ::Vararg{Any}; kws::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ QuadGK ~/.julia/packages/QuadGK/XqIlh/src/adapt.jl:213
  [9] quadgk(::Function, ::Int64, ::Vararg{Any})
    @ QuadGK ~/.julia/packages/QuadGK/XqIlh/src/adapt.jl:213
 [10] top-level scope
    @ In[2]:1

I think that the error is due to a call by quadgk of the integrand at the upper bound, which should not happen given the documentation of quadgk. https://docs.juliahub.com/QuadGK/hq5Ol/2.8.0/quadgk-examples/

This is confirmed by doing the change of variable done by quadgk manually (x=1+t/(1-t), dx = 1/(1-t)^2dt) as can be seen in the following code.

using QuadGK
function integrand(t)
    if t==1.0
        println(t)
    end

    return 1/(1+t/(1-t))^1.1 * 1/(1-t)^2
end
quadgk(integrand,0,1)
1.0
DomainError with 0.9999999999999964:
integrand produced NaN in the interval (0.9999999999999929, 1.0)

Stacktrace:
 [1] evalrule(f::typeof(integrand), a::Float64, b::Float64, x::Vector{Float64}, w::Vector{Float64}, gw::Vector{Float64}, nrm::typeof(LinearAlgebra.norm))
   @ QuadGK ~/.julia/packages/QuadGK/XqIlh/src/evalrule.jl:37
 [2] adapt(f::typeof(integrand), segs::Vector{QuadGK.Segment{Float64, Float64, Float64}}, I::Float64, E::Float64, numevals::Int64, x::Vector{Float64}, w::Vector{Float64}, gw::Vector{Float64}, n::Int64, atol::Float64, rtol::Float64, maxevals::Int64, nrm::typeof(LinearAlgebra.norm))
   @ QuadGK ~/.julia/packages/QuadGK/XqIlh/src/adapt.jl:48
 [3] do_quadgk(f::typeof(integrand), s::Tuple{Int64, Int64}, n::Int64, atol::Nothing, rtol::Nothing, maxevals::Int64, nrm::typeof(LinearAlgebra.norm), segbuf::Nothing)
   @ QuadGK ~/.julia/packages/QuadGK/XqIlh/src/adapt.jl:37
 [4] #46
   @ ~/.julia/packages/QuadGK/XqIlh/src/adapt.jl:219 [inlined]
 [5] handle_infinities(workfunc::QuadGK.var"#46#47"{Nothing, Nothing, Int64, Int64, typeof(LinearAlgebra.norm), Nothing}, f::typeof(integrand), s::Tuple{Int64, Int64})
   @ QuadGK ~/.julia/packages/QuadGK/XqIlh/src/adapt.jl:118
 [6] quadgk(::Function, ::Int64, ::Vararg{Int64}; atol::Nothing, rtol::Nothing, maxevals::Int64, order::Int64, norm::Function, segbuf::Nothing)
   @ QuadGK ~/.julia/packages/QuadGK/XqIlh/src/adapt.jl:218
 [7] quadgk(::Function, ::Int64, ::Vararg{Int64})
   @ QuadGK ~/.julia/packages/QuadGK/XqIlh/src/adapt.jl:216
 [8] top-level scope
   @ In[15]:1

Following is my code to solve two-dimensional integral in Julia.

using Statistics,Distributions,Distributed,StatsFuns,Roots,QuadGK
p=3;n0=3;n1=1;n2=4;d=0.25;
inARL=200;n=n1+n2;w1=chisqinvcdf(p,1-(n0-n1)/n2);inpa1=chisqcdf(p,w1);

function inpa(h2)
function fvv(v)
return(nchisqcdf(p,(n1/n2)v,(n*h2)/n2)*chisqpdf(p,v));end;
return(1/(1-inpa1-quadgk(fvv,w1,Inf,rtol=1e-10)[1])-inARL);end;

h=fzero(inpa,0.25,50);
outpa1=nchisqcdf(p,n1*d^2,w1);

function fz1(z)
function fv1(v)
lam=(n1/n2)*((z+sqrt(n1)*d)^2+v)+(2*n*sqrt(n1)*d*z)/(n2)+((n*d)^2)/(n2);
return(nchisqcdf(p,lam,(n*h)/n2)*chisqpdf(p-1,v));end;
return(normpdf(z)*quadgk(fv1,0,Inf,rtol=1e-10)[1]);end;
int1=quadgk(fz1,-Inf,-sqrt(w1)-sqrt(n1)*d,rtol=1e-10)[1];# problem

A simple syntax error when calling QuadGK returns an error giving the impression that QuadGK is not installed:

julia> using QuadGK

julia> g(x)=x^2
g (generic function with 1 method)

julia> quadgk(f,(0,1))
ERROR: quadgk(f, (0, 1)) has been moved to the package QuadGK.jl.
Run Pkg.add("QuadGK") to install QuadGK on Julia v0.6 and later, and then run `using QuadGK`.
Stacktrace:
 [1] #quadgk#771(::Array{Any,1}, ::Function, ::Function, ::Vararg{Any,N} where N) at ./deprecated.jl:698
 [2] quadgk(::Function, ::Tuple{Int64,Int64}) at ./deprecated.jl:698

julia> 

This can cause quite a bit confusion for students learning Julia (and even some of the more familiar!).

By the way, great package!

Since a few days I notice below issue (both Julia 0.6 and Julia 0.7-):

ERROR: QuadGK: fetch failed to get commit 064d2f0, please file an issue at https://github.com/JuliaMath/QuadGK.jl/issues
Stacktrace:
[1] fetch(::Base.LibGit2.GitRepo, ::String, ::String) at ./pkg/write.jl:25
[2] #3 at ./pkg/write.jl:53 [inlined]
[3] with(::Base.Pkg.Write.##3#4{String,String}, ::Base.LibGit2.GitRepo) at ./libgit2/types.jl:608
[4] resolve(::Dict{String,Base.Pkg.Types.VersionSet}, ::Dict{String,Dict{VersionNumber,Base.Pkg.Types.Available}}, ::Dict{String,Tuple{VersionNumber,Bool}}, ::Dict{String,Base.Pkg.Types.Fixed}, ::Dict{String,VersionNumber}, ::Set{String}) at ./pkg/entry.jl:542
[5] update(::String, ::Set{String}) at ./pkg/entry.jl:461
[6] (::Base.Pkg.Dir.##4#7{Array{Any,1},Base.Pkg.Entry.#update,Tuple{String,Set{String}}})() at ./pkg/dir.jl:36
[7] cd(::Base.Pkg.Dir.##4#7{Array{Any,1},Base.Pkg.Entry.#update,Tuple{String,Set{String}}}, ::String) at ./file.jl:70
[8] #cd#1(::Array{Any,1}, ::Function, ::Function, ::String, ::Vararg{Any,N} where N) at ./pkg/dir.jl:36
[9] update() at ./pkg/pkg.jl:228

Suggestion to make it more clear in the documentation that a function is required in QuadGK to evaluate the data and that
for purely numerical input (e.g., as recorded by instruments), other packages are recommended (e.g. Romberg.jl for equally spaced data or Trapz.jl for arbitrarily spaced data).

Right now the I and E types (returned from QuadGK) seem undefined in Segment:

# integration segment (a,b), estimated integral I, and estimated error E
struct Segment
    a::Number
    b::Number
    I
    E
end
isless(i::Segment, j::Segment) = isless(i.E, j.E)

It could be replaced by something like:

struct Segment{T<:Real}
    a::Number
    b::Number
    I::T
    E::T
end
isless(i::Segment, j::Segment) = isless(i.E, j.E)

I thought this was supported, but then this integral should give 2πi by the residue theorem:

julia> quadgk(z -> cos(z)/z, 1, im, -1, -im)
(0.50046301807924 + 4.71238898038469im, 4.027482497221434e-9)

and the result looks clearly wrong. 😢

Thank you for a great package!
I would like to use Zygote on quadgk but below is a self-contained example, where it does not work because quadgk calls setindex!. Can I do anything to work around this? I coded up a small example of differentiating quadgk of a constant function and it did work so it seems that there is hope for Zygote and QuadGK to play together nicely.

using Zygote, QuadGK

"P(0<X<x)"
function F(c, x)
    if x < c
        return 0.
    elseif x < 1.
        return exp(c)*(x-c)^2/2
    elseif x < 1. + c
        return exp(c)*(1-c)*(2x - c - 1.)/2.
    elseif x < 2.
        return exp(c)*(4x - x^2 - 2c - 2)/2.
    else
        return exp(c)*(1-c)
    end
end

"pdf of X"
function f(c, x)
    if x < c || x > 2.
        return 0.
    elseif x < 1.
        return exp(c)*(x-c)
    elseif x < 1. + c
        return exp(c)*(1-c)
    else
        return exp(c)*(2-x)
    end
end

function D(c)
    exp(c)*(1-c)
end

function W(c1, c2)
    val, acc = quadgk(t -> F(c1, t)*f(c2, t), c1, 2)
    val
end

function prob_of_win(c1, c2)
    ((1-D(c1))*D(c2) + W(c2, c1))/(D(c1) + D(c2) - D(c1)*D(c2))
end

gradient(prob_of_win, 0.25, 0.25) # ERROR: Mutating arrays is not supported -- called setindex!(Vector{QuadGK.Segment{Float64, Float64, Float64}}, ...)