Duals in the boundary values: avoid infinite recursion #61
Add this suggestion to a batch that can be applied as a single commit.
This suggestion is invalid because no changes were made to the code.
Suggestions cannot be applied while the pull request is closed.
Suggestions cannot be applied while viewing a subset of changes.
Only one suggestion per line can be applied in a batch.
Add this suggestion to a batch that can be applied as a single commit.
Applying suggestions on deleted lines is not supported.
You must change the existing code in this line in order to create a valid suggestion.
Outdated suggestions cannot be applied.
This suggestion has been applied or marked resolved.
Suggestions cannot be applied from pending reviews.
Suggestions cannot be applied on multi-line comments.
Suggestions cannot be applied while the pull request is queued to merge.
Suggestion cannot be applied right now. Please check back later.
The previous behavior led to https://discourse.julialang.org/t/error-with-forwarddiff-in-turing/88633/16 . Dual numbers would just hit an infinite recursion because
float(::Dual)::Dual !<: AbstractFloat. This should handle more cases where float(T) !<: AbstractFloat too, like TrackedReal.Though I know that this probably isn't the best way to handle it. In a higher level interface like Integrals.jl, I think it would be better to detect such differentiation by boundary terms and specialize it. It's actually kind of funny and would be a good homework problem. If you put a dual number into the boundary with dual part one, i.e. , so then you integrate (f + df/dt eps) dt and by linearity of the integral it's just fundamental theorem of calculus approximated by GK! However, the argument above assumes that the dual part generated in the rules is exactly 1 (or rather, matches the partials of the boundary term so that it weighs it correctly). The algorithm doesn't exactly work like that, but it ends up being mathematically the same. What happens naturally in the code is that the rule values are multiplied by segment sizes as weights
f(a + (1-x[2i-1])*s)and it is there that thexis "imparted with" the correct dualness to then give the integral of the derivative.Of course, algorithmically it should probably just see this is the case, integrate f, and then slap f into the dual part. I'm not sure if QuadGK.jl should be mucking with dual numbers and autodiff overloads though, so I'll just do that at a higher level. But it's something to think about here.
The reason to not do a complete overload is because for differentiating with respect to other quantities like parameters, you can fuse the primal and derivative calls rather than having them as two separate integrals, so it's definitely faster to differentiate the algorithm here. Except for handling the boundary.
I don't know, it's pretty cool haha. Such a beautifully convoluted way to approximate f(x).