Merge pull request #37 from tlycken/teh/fasterinterp

Fix performance problems with interpolation

REQUIRE

@@ -1,3 +1,4 @@
 julia 0.4-
 
 WoodburyMatrices
+Ratios

src/Interpolations.jl

@@ -19,9 +19,9 @@ export
     # b-splines/b-splines.jl
     # extrapolation/extrapolation.jl
 
-using WoodburyMatrices
+using WoodburyMatrices, Ratios
 
-import Base: size, getindex
+import Base: convert, size, getindex, promote_rule
 
 abstract InterpolationType
 abstract GridType

src/b-splines/b-splines.jl

@@ -10,7 +10,7 @@ abstract Degree{N}
 immutable BSpline{D<:Degree} <: InterpolationType end
 BSpline{D<:Degree}(::Type{D}) = BSpline{D}
 
-type BSplineInterpolation{T,N,TCoefs,IT<:BSpline,GT<:GridType,pad} <: AbstractInterpolation{T,N,IT,GT}
+immutable BSplineInterpolation{T,N,TCoefs,IT<:BSpline,GT<:GridType,pad} <: AbstractInterpolation{T,N,IT,GT}
     coefs::Array{TCoefs,N}
 end
 function BSplineInterpolation{N,TCoefs,TWeights<:Real,IT<:BSpline,GT<:GridType,pad}(::Type{TWeights}, A::AbstractArray{TCoefs,N}, ::Type{IT}, ::Type{GT}, ::Val{pad})
@@ -26,7 +26,11 @@ function BSplineInterpolation{N,TCoefs,TWeights<:Real,IT<:BSpline,GT<:GridType,p
     BSplineInterpolation{T,N,TCoefs,IT,GT,pad}(A)
 end
 
-size{T,N,TCoefs,IT,GT,pad}(itp::BSplineInterpolation{T,N,TCoefs,IT,GT,pad}, d) = (d <= N ? size(itp.coefs, d) - 2*pad : 1)::Int
+@generated function size{T,N,TCoefs,IT,GT,pad}(itp::BSplineInterpolation{T,N,TCoefs,IT,GT,pad}, d)
+    quote
+        d <= $N ? size(itp.coefs, d) - 2*$pad : 1
+    end
+end
 
 function interpolate{TWeights,IT<:BSpline,GT<:GridType}(::Type{TWeights}, A, ::Type{IT}, ::Type{GT})
     Apad, Pad = prefilter(TWeights, A, IT, GT)

src/b-splines/indexing.jl

@@ -3,7 +3,9 @@ using Base.Cartesian
 import Base.getindex
 
 function getindex_impl{T,N,TCoefs,IT<:BSpline,GT<:GridType,Pad}(itp::Type{BSplineInterpolation{T,N,TCoefs,IT,GT,Pad}})
+    meta = Expr(:meta, :inline)
     quote
+        $meta
         @nexprs $N d->(x_d = xs[d])
 
         # Calculate the indices of all coefficients that will be used
@@ -49,6 +51,10 @@ end
     getindex_impl(itp)
 end
 
+@generated function getindex{T,N}(itp::BSplineInterpolation{T,N}, index::CartesianIndex{N})
+    :(getindex(itp, $(Base.IteratorsMD.cartindex_exprs((index,), (:index,))...)))
+end
+
 offsetsym(off, d) = off == -1 ? symbol(string("ixm_", d)) :
                     off ==  0 ? symbol(string("ix_", d)) :
                     off ==  1 ? symbol(string("ixp_", d)) :

src/b-splines/quadratic.jl

@@ -33,19 +33,19 @@ function coefficients{Q<:Quadratic}(::Type{BSpline{Q}}, N, d)
     symm, sym, symp =  symbol(string("cm_",d)), symbol(string("c_",d)), symbol(string("cp_",d))
     symfx = symbol(string("fx_",d))
     quote
-        $symm = 1//2 * ($symfx - 1//2)^2
-        $sym  = 3//4 - $symfx^2
-        $symp = 1//2 * ($symfx + 1//2)^2
+        $symm = sqr($symfx - SimpleRatio(1,2))/2
+        $sym  = SimpleRatio(3,4) - sqr($symfx)
+        $symp = sqr($symfx + SimpleRatio(1,2))/2
     end
 end
 
 function gradient_coefficients{Q<:Quadratic}(::Type{Q}, N, d)
     symm, sym, symp =  symbol(string("cm_",d)), symbol(string("c_",d)), symbol(string("cp_",d))
     symfx = symbol(string("fx_",d))
     quote
-        $symm = $symfx - 1//2
+        $symm = $symfx - SimpleRatio(1,2)
         $sym = -2 * $symfx
-        $symp = $symfx + 1//2
+        $symp = $symfx + SimpleRatio(1,2)
     end
 end
 
@@ -67,8 +67,8 @@ padding{BC<:BoundaryCondition}(::Type{BSpline{Quadratic{BC}}}) = Val{1}()
 padding(::Type{BSpline{Quadratic{Periodic}}}) = Val{0}()
 
 function inner_system_diags{T,Q<:Quadratic}(::Type{T}, n::Int, ::Type{Q})
-    du = fill(convert(T, 1//8), n-1)
-    d = fill(convert(T, 3//4), n)
+    du = fill(convert(T, SimpleRatio(1,8)), n-1)
+    d = fill(convert(T, SimpleRatio(3,4)), n)
     dl = copy(du)
     (dl,d,du)
 end
@@ -147,7 +147,9 @@ function prefiltering_system{T,TCoefs,GT<:GridType}(::Type{T}, ::Type{TCoefs}, n
     colspec[1,n] = colspec[2,1] = 1
     valspec = zeros(T,2,2)
     # [1,n]            [n,1]
-    valspec[1,1] = valspec[2,2] = 1//8
+    valspec[1,1] = valspec[2,2] = SimpleRatio(1,8)
 
     Woodbury(lufact!(Tridiagonal(dl, d, du), Val{false}), rowspec, valspec, colspec), zeros(TCoefs, n)
 end
+
+sqr(x) = x*x