Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Errors with BigFloat boundaries #60

Closed
araujoms opened this issue Feb 23, 2024 · 6 comments · Fixed by #61
Closed

Errors with BigFloat boundaries #60

araujoms opened this issue Feb 23, 2024 · 6 comments · Fixed by #61

Comments

@araujoms
Copy link
Contributor

The following function errors when I call it as erf(BigFloat). It should work, no? It doesn't complain if I call it with other funny types that are nevertheless isbits.

function erf(T)
    lb, ub = T(0), T(1)
    f(x) = exp(-x^2)
    sol = hquadrature(f, lb, ub)
    return sol[1] 
end
@stevengj
Copy link
Member

If you want to do 1d BigFloat quadrature, I would highly recommend using QuadGK instead, since QuadGK lets you increase the quadrature order (it is a good idea to increase the order roughly in proportion to the precision).

That being said, I agree that this should work (even if it is suboptimal). Maybe we can make use of JuliaArrays/StaticArrays.jl#799

@araujoms
Copy link
Contributor Author

The actual integral I need to solve is 2d, I only simplified it for the MWE.

Do you mean for example, in line 67 of HCubature.jl, declare ma = similar(a) instead of ma = MVector(a)?

@stevengj
Copy link
Member

stevengj commented Feb 23, 2024

Yes. similar is defined to return a mutable type, and apparently it should nowadays do the right thing with SVector:

julia> a = SA[1,2,3]
3-element SVector{3, Int64} with indices SOneTo(3):
 1
 2
 3

julia> ma = similar(a)
3-element MVector{3, Int64} with indices SOneTo(3):
 4537519120
  167772160
          0

julia> a = SA[big(1),2,3]
3-element SVector{3, BigInt} with indices SOneTo(3):
 1
 2
 3

julia> ma = similar(a)
3-element SizedVector{3, BigInt, Vector{BigInt}} with indices SOneTo(3):
 #undef
 #undef
 #undef

@stevengj
Copy link
Member

stevengj commented Feb 23, 2024

Sorry, not similar because we need to copy the contents. Maybe Base.copymutable(a)?

@stevengj
Copy link
Member

stevengj commented Feb 23, 2024

The annoyance here is that Base.copymutable(a), which calls similar and then copyto!, may be less efficient than MVector(a) for the cases where MVector does work.

It seems like StaticArrays.jl should ideally implement a specialized Base.copymutable method? JuliaArrays/StaticArrays.jl#1247

But for now Base.copymutable seems like the most reliable method. Would be good to benchmark it to check for regressions for Float64 integration — I doubt it is noticeable.

@araujoms
Copy link
Contributor Author

StaticArrays is above my pay grade.

But if this is the only place then I'm sure there will be no difference in performance. Looking around a bit, line 15 of genz-malik.jl will also certainly cause trouble. But there similar is enough, no copying is needed.

I can give it a try tomorrow; to see if I can fix it and benchmark it. Can you suggest a test case for benchmarking?

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging a pull request may close this issue.

2 participants