Small updates to Complete Polynomials #26
Conversation
Codecov Report@@ Coverage Diff @@
## master #26 +/- ##
==========================================
- Coverage 84.25% 81.31% -2.95%
==========================================
Files 9 9
Lines 972 1038 +66
==========================================
+ Hits 819 844 +25
- Misses 153 194 +41
Continue to review full report at Codecov.
|
|
Bad news... There's going to be more in this PR. My most recent commit adds functionality for computing derivative of basis functions for vectors of variables. Use is shown below: z = [1.0, 2.0, 3.0]
der_wrt_x = complete_polynomial_der(z, 3, 1)
der_wrt_y = complete_polynomial_der(z, 3, 2)
der_wrt_z = complete_polynomial_der(z, 3, 3)produces Notice computing the derivatives of basis functions is a little slower than computing the basis functions themselves. There may be a way to close this gap. |
|
Now there is a derivative function for matrices as well. This one seems to pay a higher performance penalty than the vector version. Would love for someone who knows more voodoo than me to have a look to see whether it can be sped up. |
src/BasisMatrices.jl
Outdated
| @@ -70,7 +70,8 @@ export BasisFamily, Cheb, Lin, Spline, Basis, | |||
| # functions | |||
| export nodes, get_coefs, funfitxy, funfitf, funeval, evalbase, | |||
| derivative_op, row_kron, evaluate, fit!, update_coefs!, | |||
| complete_polynomial, complete_polynomial!, n_complete | |||
| complete_polynomial, complete_polynomial!, n_complete, | |||
| complete_polynomial_der, complete_polynomial_der! | |||
There was a problem hiding this comment.
I was thinking that given that you've introduced a Derivative{D} type, that we could get away with only having the complete_polynomial and complete_polynomial! functions be exported. We could then infer if the user wanted a derivative based on whether or not a Derivative type was passed and/or a derivative keyword argument?
There was a problem hiding this comment.
Great point.
I don't think we can use a keyword argument for derivatives (it wouldn't dispatch on it) -- I think the following protocol seems sensible
# computes 3rd order complete polynomial for z
complete_polynomial(z, 3)
# computes derivatives with respect to variable 1 of 3rd order complete polynomial
complete_polynomial(z, 3, 1)
Thoughts?
| ) | ||
| out | ||
| end | ||
| end |
There was a problem hiding this comment.
I 100% agree that having methods that operate on a single point (passed as a Vector) is beneficial. We need to support this.
However, this seems like a lot of code duplication...
I wonder if the small difference in timing between the new method here for z::Vector and calling on z' has to do with copying the memory for doing the transpose? If so, maybe we could try defining the method on vector so it just calls the Matrix one with reshape(z, 1, length(z)) so there isn't a copy?
There was a problem hiding this comment.
I agree that it is kind of a lot of code duplication -- Some of it was a playground for me to try and learn some stuff. A bit surprisingly, complete_polynomial(reshape(z, 1, length(z)), 4) actually does worse than complete_polynomial(z', 4) does. It's also slower even if I just pass in the already reshaped z argument.
FWIW, I use a slightly different way of constructing the points for vectors than for matrices (slight might even be an overstatement of the difference), but it does seem to continue being faster even as the number of dimensions and degree of polynomials are ramped up (we're talking about microseconds here...).
Thoughts? If specializing turns out to keep it faster then my vote would be to keep it separate for now, but I strongly suspect there are ways to close the gaps (I just haven't figured out how).
src/complete.jl
Outdated
| function complete_polynomial{T}(z::Vector{T}, d::Int) | ||
| nvar = length(z) | ||
| out = Array(T, n_complete(nvar, d)) | ||
| complete_polynomial!(z, Degree{d}(), out)::Vector{T} |
There was a problem hiding this comment.
I realized that we can avoid the type assert by just returning out (implicitly of course 😉 ) after calling complete_polynomial! instead of relying on the return variable of complete_polynomial!.
There was a problem hiding this comment.
👏
Great. I'll fix this everywhere (with explicit returns).
There was a problem hiding this comment.
Hmmm. We still get intermediate type instability because it doesn't know what type gets returned by complete_polynomial!(z, Degree{d}(), out)... Is that a problem? It really seems like it should be able to infer that type.
There was a problem hiding this comment.
I don't think that the inability to infer the return type of complete_polynomial! should matter.
What I had in mind was this:
function complete_polynomial{T}(z::Vector{T}, d::Int)
nvar = length(z)
out = Array(T, n_complete(nvar, d))
complete_polynomial!(z, Degree{d}(), out)
out
endWe don't really care what the output of complete_polynomial! is because we never use it.
This adds some vector methods for generating basis functions of complete polynomials. Specializing on vectors gives a (small) speedup and is more convenient on some occasions.
Define the following
Timings on this branch are
Time on master is