Various suggestions #19
Conversation
| end # module IntegralsGaussTuranExt | ||
| function GaussTuranQuadrature.GaussTuranComputeRule(::Type{T}, n::Int, s::Int; kws...) where {T} | ||
| rhs = inv.(T.(1:(2 * (s + 1) * n + 1))) # 1 ./ collect(range(one(T), T(2 * (s + 1) * n + 1))) | ||
| ϕ = let v = Val(2s + 1) |
There was a problem hiding this comment.
Can you explain what is going on here? I understand what this code is supposed to do but not how it does it 😅 I'm not very familiar with the do and let keywords
There was a problem hiding this comment.
The first line converts integers to the appropriate float type and inverts them to compute the analytical integrals. The second line uses a let block to hopefully improve type-stability within the closure, since Val(2s+1) is computed at runtime. The do syntax creates an anonymous function used to evaluate each element of the tuple and in this case also compute derivatives via an iteration. Here I am using a let block since _f gets reassigned in the closure, which can lead to performance issues because of the language. I should still check to see if this is type stable, which would be the main goal as far as performance goes
There was a problem hiding this comment.
I may not be very consistent with my usage of let blocks. Usually I would add them if profiling reveals a type instability, but here I just added them pre-emptively
There was a problem hiding this comment.
Thanks for the explanations, I have a few more questions:
- Why do you need the
*1at the end of line 198? - If I understand correctly,
_fis the value of the 'current derivative' and__fthe value of the 'next derivative'. If so, can you add an explanation of that in the code?
| order = f.order === nothing ? Val(n_derivs(I)) : | ||
| f.order isa Integer ? Val(f.order) : | ||
| f.order | ||
| __f = order isa Val{1} ? x -> (_f(x),) : x -> value(derivatives(_f, x, 1, order)) |
There was a problem hiding this comment.
Can you also explain this part?
There was a problem hiding this comment.
Since TaylorDiff.jl expects the order to be a Val(n), the first line infers the order from the rule if the user doesn't provide an order and otherwise uses the user's order. The second line uses the TaylorDiff.jl API to construct the integrand with derivatives, although the no-derivative case requires special treatment.
There was a problem hiding this comment.
Thanks. Some questions:
- Is this nested/chained use of ternary operators common practice? Then it is just something for me to get used to. Otherwise it would be better to write it out with
if else. - Would it be an idea to store
f.orderas aVal?
| function GaussTuranRule(n::Integer, s::Integer; domain::Tuple{<:Number, <:Number} = (0, 1)) | ||
| dom = promote(domain...) | ||
| Δx = dom[2] - dom[1] | ||
| T = typeof(float(real(one(Δx)))) |
There was a problem hiding this comment.
Why this convoluted way of determining the type?
There was a problem hiding this comment.
TLDR; this is how QuadGK.jl does it.
The issue is that you want a real, unitless, floating-point type for optimization, but the domain could be specified with integer endpoints, or could be a complex segment, so this is the best way to get a type that works.
| return evalrule(f, W, X, a, b) | ||
| end | ||
|
|
||
| function evalrule(f, W, X, a, b) |
There was a problem hiding this comment.
Can you also explain this function?
There was a problem hiding this comment.
Your previous version of the function used out = zero(eltype(X)), which is wrong if the integrand returns some generic type supporting addition and scalar multiplication. This version of the function correctly initializes out but has to unwrap the first iteration of the loop to do so.
|
Regarding the accuracy of the tabulated rules: I hadn't looked into that yet, I just took a part of the parameter space that didn't take too long to compute, and the supplied optimization parameters are not very well informed. It would indeed be nice to only tabulate the rule if it passes an accuracy criterion (and maybe raise an error if the user tries to use such a rule which apparently is hard to derive).
and their results aren't very good as you said. It's also unclear whether they mean 120 digits in binary or in decimal. The optimization algorithm can also be performed with |
|
Thanks for your comments about the tabulated rules. I think it's totally fine to tabulate the rule only when the accuracy criterion is passed because these quadrature rules are rather niche and as long as you have enough rules for your needs then I think the package has already accomplished more than what some of the original literature presents. As for the accuracy of the endpoints, I can see why doing the calculation at a higher precision and then reducing precision as necessary for subsequent calculations is desirable, however I don't think this should be done automatically for the user, except for the tabulated rules where we can actually guarantee full precision. (Also, digits usually mean decimal places, while bits correspond to binary.) So, I think I will print a message whenever a quadrature rule is being computed and say that it is up to the user to verify its accuracy. This check could also be done a posteriori as a convenience for the user, but is most likely non-essential as long as we inform the user about which of the tabulated rules they can choose from. For reference, similar packages which obtain quadrature rules from nonlinear optimization will generally error if you request a rule that is not already computed (see https://github.com/JoshuaTetzner/GeneralizedGaussianQuadrature.jl/blob/main/src/quadrature/gqlog.jl) because this is by no means a solved problem. |
|
Thanks a lot for your answers to my questions, I learned a lot from it (: |
Hi,
This PR does the following:
The main issue is the accuracy of the monomial rules, which according to my test are not all accurate to within 3 digits of the working precision. A "fix" is to tabulate only the rules that pass the test or otherwise to find a way to make the optimization succeed.
Accuracy test results