optimize and test Bernstein

in place De Casteljau
combine De Casteljau's for _derivatives_1d
add test

src/Polynomials/BernsteinBases.jl

@@ -26,8 +26,6 @@ BernsteinBasis(args...) = UniformPolyBasis(Bernstein, args...)
 
 # 1D evaluation implementation
 
-# TODO Optimize with in-place De Casteljau
-
 """
     binoms(::Val{K})
 
@@ -36,22 +34,53 @@ Returns the tuple of binomials ( C₍ₖ₀₎, C₍ₖ₁₎, ..., C₍ₖₖ
 binoms(::Val{K}) where K = ntuple( i -> binomial(K,i-1), Val(K+1))
 
 
-# jth Bernstein poly of order K at x:
-# Bᵏⱼ(x) = binom(K,j) * x^j * (1-x)^(K-j)
-function _evaluate_1d!(::Type{Bernstein},::Val{K}, v::AbstractMatrix{T},x,d) where {K,T<:Number}
-  b = binoms(Val(K))
-  n = K + 1
-  @inbounds λ1 = x[d]
-  λ2 = one(T) - λ1
-
-  for i in 1:n
-    j = i-1
-    λ1_j = λ1^(j)
-    λ2_j = λ2^(K-j)
-    @inbounds v[d,i] = b[i] * λ1_j * λ2_j # order K
+function _evaluate_1d!(::Type{Bernstein},::Val{0},v::AbstractMatrix{T},x,d) where {T<:Number}
+  @inbounds v[d,1] = one(T)
+end
+
+@inline function _De_Casteljau_step_1D!(v,d,i,λ1,λ2)
+  # i = k+1
+
+  # vₖ <- xvₖ₋₁            # Bᵏₖ(x) = x*Bᵏ⁻¹ₖ₋₁(x)
+  v[d,i] = λ2*v[d,i-1]
+  # vⱼ <- xvⱼ₋₁ + (1-x)vⱼ  # Bᵏⱼ(x) = x*Bᵏ⁻¹ⱼ₋₁(x) + (1-x)*Bᵏ⁻¹ⱼ(x) for j = k-1, k-2, ..., 1
+  for l in i-1:-1:2
+    v[d,l] = λ2*v[d,l-1] + λ1*v[d,l]
   end
+  # v₀ <- (1-x)v₀          # Bᵏ₀(x) = (1-x)*Bᵏ⁻¹₀(x)
+  v[d,1] = λ1*v[d,1]
 end
 
+# jth Bernstein poly of order K at x:
+# Bᵏⱼ(x) = binom(K,j) * x^j * (1-x)^(K-j) = x*Bᵏ⁻¹ⱼ₋₁(x) + (1-x)*Bᵏ⁻¹ⱼ(x)
+function _evaluate_1d!(::Type{Bernstein},::Val{K},v::AbstractMatrix{T},x,d) where {K,T<:Number}
+  @inbounds begin
+    n = K + 1 # n > 1
+    λ2 = x[d]
+    λ1 = one(T) - λ2
+
+    # In place De Casteljau: init with B¹₀(x)=x and B¹₁(x)=1-x
+    v[d,1] = λ1
+    v[d,2] = λ2
+
+    for i in 3:n
+      _De_Casteljau_step_1D!(v,d,i,λ1,λ2)
+      ## vₖ <- xvₖ₋₁            # Bᵏₖ(x) = x*Bᵏ⁻¹ₖ₋₁(x)
+      #v[d,i] = λ2*v[d,i-1]
+      ## vⱼ <- xvⱼ₋₁ + (1-x)vⱼ  # Bᵏⱼ(x) = x*Bᵏ⁻¹ⱼ₋₁(x) + (1-x)*Bᵏ⁻¹ⱼ(x) for j = k-1, k-2, ..., 1
+      #for l in i-1:-1:2
+      #  v[d,l] = λ2*v[d,l-1] + λ1*v[d,l]
+      #end
+      ## v₀ <- (1-x)v₀          # Bᵏ₀(x) = (1-x)*Bᵏ⁻¹₀(x)
+      #v[d,1] = λ1*v[d,1]
+    end
+  end
+  # still optimisable for K > 2/3:
+  # - compute bj = binoms(k,j) at compile time (binoms(Val(K)) function)
+  # - compute vj = xʲ*(1-x)ᴷ⁻ʲ recursively in place like De Casteljau (saving half the redundant multiplications)
+  # - do it in a stack allocated cache (MVector, Bumber.jl)
+  # - @simd affect bj * vj in v[d,i] for all j
+end
 
 function _gradient_1d!(::Type{Bernstein},::Val{0},g::AbstractMatrix{T},x,d) where {T<:Number}
   @inbounds g[d,1] = zero(T)
@@ -66,19 +95,20 @@ end
 # (Bᵏⱼ)'(x) = K * ( Bᵏ⁻¹ⱼ₋₁(x) - Bᵏ⁻¹ⱼ(x) )
 #           = K * x^(j-1) * (1-x)^(K-j-1) * ((1-x)*binom(K-1,j-1) - x*binom(K-1,j))
 function _gradient_1d!(::Type{Bernstein},::Val{K}, g::AbstractMatrix{T},x,d) where {K,T<:Number}
-  b = binoms(Val(K-1))
-  n = K + 1
-
-  @inbounds λ1 = x[d]
-  λ2 = one(T) - λ1
-
-  @inbounds g[1] =-K * λ2^(K-1)
-  @inbounds g[n] = K * λ1^(K-1)
-  for i in 2:n-1
-    j = i-1
-    λ1_j = λ1^(j-1)
-    λ2_j = λ2^(K-j-1)
-    @inbounds g[d,i] = K * λ1_j * λ2_j *(λ2*b[i-1] - λ1*b[i]) # order K-1
+  @inbounds begin
+    n = K + 1 # n > 2
+
+    # De Casteljau for Bᵏ⁻¹ⱼ for j = k-1, k-2, ..., 1
+    _evaluate_1d!(Bernstein,Val(K-1),g,x,d)
+
+    # gₖ <- K*gₖ₋₁         # ∂ₓBᵏₖ(x) = K*Bᵏ⁻¹ₖ₋₁(x)
+    g[d,n] = K*g[d,n-1]
+    # gⱼ <- K(gⱼ₋₁ + gⱼ)   # ∂ₓBᵏⱼ(x) = K(Bᵏ⁻¹ⱼ₋₁(x) - Bᵏ⁻¹ⱼ(x)) for j = k-1, k-2, ..., 1
+    for l in n-1:-1:2
+      g[d,l] = K*(g[d,l-1] - g[d,l])
+    end
+    # g₀ <- K*g₀           # ∂ₓBᵏ₀(x) = -K*Bᵏ⁻¹₀(x)
+    g[d,1] = -K*g[d,1]
   end
 end
 
@@ -103,21 +133,96 @@ end
 #                  - 2x*(1-x)*binom(K-2,j-1) + (x)^2*binom(K-2,j)
 #              )
 function _hessian_1d!(::Type{Bernstein},::Val{K},h::AbstractMatrix{T},x,d) where {K,T<:Number}
-  b = binoms(Val(K-2))
-  n = K + 1
-  C = K*(K-1)
-
-  @inbounds λ1 = x[d]
-  λ2 = one(T) - λ1
-
-  @inbounds h[1] = C *    λ2^(K-2)
-  @inbounds h[2] = C * (-2λ2^(K-2) + (K-2)*λ2^(K-3)*λ1)
-  @inbounds h[n-1]=C * (-2λ1^(K-2) + (K-2)*λ1^(K-3)*λ2)
-  @inbounds h[n] = C *    λ1^(K-2)
-  for i in 3:n-2
-    j = i-1
-    λ1_j = λ1^(j-2)
-    λ2_j = λ2^(K-j-2)
-    @inbounds h[d,i] = C * λ1_j * λ2_j *(λ2*λ2*b[i-2] -2λ2*λ1*b[i-1] + λ1*λ1*b[i]) # order K-2
+  @inbounds begin
+    n = K + 1 # n > 3
+    KK = K*(K-1)
+
+    # De Casteljau for Bᵏ⁻²ⱼ for j = k-2, k-3, ..., 1
+    _evaluate_1d!(Bernstein,Val(K-2),h,x,d)
+
+    # hₖ   <- K(K-1)*hₖ₋₂
+    h[d,n] = KK*h[d,n-2]
+    # hₖ₋₁ <- K(K-1)*(-2*hₖ₋₁ + hₖ₋₂)
+    h[d,n-1] = KK*( h[d,n-3] -2*h[d,n-2] )
+
+    # hⱼ <- K(K-1)(hⱼ₋₂ -2hⱼ₋₁ + hⱼ)
+    for l in n-2:-1:3
+      h[d,l] = KK*( h[d,l-2] -2*h[d,l-1] + h[d,l] )
+    end
+
+    # h₁ <- K(K-1)*(-2h₀ + h₁)
+    h[d,2] = KK*( -2*h[d,1] + h[d,2] )
+    # h₀ <- K(K-1)*h₀
+    h[d,1] = KK*h[d,1]
+  end
+end
+
+function _derivatives_1d!(::Type{Bernstein},v::Val_01,t::NTuple{2},x,d)
+  @inline _evaluate_1d!(Bernstein, v, t[1], x, d)
+  @inline _gradient_1d!(Bernstein, v, t[2], x, d)
+end
+
+function _derivatives_1d!(::Type{Bernstein},::Val{K},t::NTuple{2},x,d) where K
+  @inbounds begin
+    n = K + 1 # n > 2
+    v, g = t
+
+    λ2 = x[d]
+    λ1 = one(eltype(v)) - λ2
+
+    # De Casteljau for Bᵏ⁻¹ⱼ for j = k-1, k-2, ..., 1
+    _evaluate_1d!(Bernstein,Val(K-1),v,x,d)
+
+    # Compute gradients as _gradient_1d!
+    g[d,n] = K*v[d,n-1]
+    @simd for l in n-1:-1:2
+      g[d,l] = K*(v[d,l-1] - v[d,l])
+    end
+    g[d,1] = -K*v[d,1]
+
+    # Last step of De Casteljau for _evaluate_1d!
+    _De_Casteljau_step_1D!(v,d,n,λ1,λ2)
+  end
+end
+
+function _derivatives_1d!(::Type{Bernstein},v::Val_012,t::NTuple{3},x,d)
+  @inline _evaluate_1d!(Bernstein, v, t[1], x, d)
+  @inline _gradient_1d!(Bernstein, v, t[2], x, d)
+  @inline _hessian_1d!( Bernstein, v, t[3], x, d)
+end
+
+function _derivatives_1d!(::Type{Bernstein},::Val{K},t::NTuple{3},x,d) where K
+  @inbounds begin
+    n = K + 1 # n > 3
+    v, g, h = t
+
+    KK = K*(K-1)
+    λ2 = x[d]
+    λ1 = one(eltype(v)) - λ2
+
+    # De Casteljau until Bᵏ⁻²ⱼ ∀j
+    _evaluate_1d!(Bernstein,Val(K-2),v,x,d)
+
+    # Compute hessians as _hessian_1d!
+    h[d,n] = KK*v[d,n-2]
+    h[d,n-1] = KK*( v[d,n-3] -2*v[d,n-2] )
+    @simd for l in n-2:-1:3
+      h[d,l] = KK*( v[d,l-2] -2*v[d,l-1] + v[d,l] )
+    end
+    h[d,2] = KK*( -2*v[d,1] + v[d,2] )
+    h[d,1] = KK*v[d,1]
+
+    # One step of De Casteljau to get Bᵏ⁻¹ⱼ ∀j
+    _De_Casteljau_step_1D!(v,d,n-1,λ1,λ2)
+
+    # Compute gradients as _gradient_1d!
+    g[d,n] = K*v[d,n-1]
+    @simd for l in n-1:-1:2
+      g[d,l] = K*(v[d,l-1] - v[d,l])
+    end
+    g[d,1] = -K*v[d,1]
+
+    # Last step of De Casteljau for _evaluate_1d!
+    _De_Casteljau_step_1D!(v,d,n,λ1,λ2)
   end
 end

src/Polynomials/PolynomialInterfaces.jl

@@ -204,7 +204,6 @@ function _gradient_nd!(
   @abstractmethod
 end
 
-
 """
     _hessian_nd!(b,xi,r,i,c,g,h,s)
 
@@ -312,6 +311,10 @@ function _hessian_1d!(::Type{<:Polynomial},::Val{K},h::AbstractMatrix{T},x,d) wh
   @abstractmethod
 end
 
+# Dispatch helpers for base cases
+const Val_01  = Union{Val{0},Val{1}}
+const Val_012 = Union{Val{0},Val{1},Val{2}}
+
 """
     _derivatives_1d!(PT::Type{<:Polynomial}, ::Val{K}, (c,g,...), x, d)
 
@@ -323,22 +326,21 @@ _gradient_1d!(PT, Val(K), g, x d)
 ```
 but with possible performance optimization.
 """
-function _derivatives_1d!(  ::Type{<:Polynomial},::Val{K},t::NTuple{N},x,d) where {K,N}
+function _derivatives_1d!(  ::Type{<:Polynomial},v::Val,t::NTuple{N},x,d) where N
   @abstractmethod
 end
 
-function _derivatives_1d!(PT::Type{<:Polynomial},::Val{K},t::NTuple{1},x,d) where K
-  _evaluate_1d!(PT, Val(K), t[1], x, d)
+function _derivatives_1d!(PT::Type{<:Polynomial},v::Val,t::NTuple{1},x,d)
+  @inline _evaluate_1d!(PT, v, t[1], x, d)
 end
 
-function _derivatives_1d!(PT::Type{<:Polynomial},::Val{K},t::NTuple{2},x,d) where K
-  _evaluate_1d!(PT, Val(K), t[1], x, d)
-  _gradient_1d!(PT, Val(K), t[2], x, d)
+function _derivatives_1d!(PT::Type{<:Polynomial},v::Val,t::NTuple{2},x,d)
+  @inline _evaluate_1d!(PT, v, t[1], x, d)
+  @inline _gradient_1d!(PT, v, t[2], x, d)
 end
 
-function _derivatives_1d!(PT::Type{<:Polynomial},::Val{K},t::NTuple{3},x,d) where K
-  _evaluate_1d!(PT, Val(K), t[1], x, d)
-  _gradient_1d!(PT, Val(K), t[2], x, d)
-  _hessian_1d!( PT, Val(K), t[3], x, d)
+function _derivatives_1d!(PT::Type{<:Polynomial},v::Val,t::NTuple{3},x,d)
+  @inline _evaluate_1d!(PT, v, t[1], x, d)
+  @inline _gradient_1d!(PT, v, t[2], x, d)
+  @inline _hessian_1d!( PT, v, t[3], x, d)
 end
-

test/PolynomialsTests/BernsteinBasesTests.jl

@@ -11,13 +11,30 @@ using ForwardDiff
 
 np = 3
 x = [Point(0.),Point(1.),Point(.4)]
-xi = x[1]
+x1 = x[1]
 
 # Only test 1D evaluations as tensor product structure is tested in monomial tests
 #
 V = Float64
-G = gradient_type(V,xi)
-H = gradient_type(G,xi)
+G = gradient_type(V,x1)
+H = gradient_type(G,x1)
+
+function test_internals(order,x,bx,∇bg,Hbx)
+  sz = (1,order+1)
+  for (i,xi) in enumerate(x)
+    v2 = zeros(sz)
+    Polynomials._evaluate_1d!(Bernstein,Val(order),v2,xi,1)
+    @test all( [ bxi[1]≈vxi[1] for (bxi,vxi) in zip(bx[i,:],v2[:,1]) ] )
+
+    g2 = zeros(sz)
+    Polynomials._gradient_1d!(Bernstein,Val(order),g2,xi,1)
+    @test all( [ bxi[1]≈vxi[1] for (bxi,vxi) in zip(∇bx[i,:],g2[:,1]) ] )
+
+    h2 = zeros(sz)
+    Polynomials._hessian_1d!(Bernstein,Val(order),h2,xi,1)
+    @test all( [ bxi[1]≈vxi[1] for (bxi,vxi) in zip(Hbx[i,:],h2[:,1]) ] )
+  end
+end
 
 
 # order 0 degenerated case
@@ -27,13 +44,12 @@ b = BernsteinBasis(Val(1),V,order)
 @test get_order(b) == 0
 @test get_orders(b) == (0,)
 
-v = V[1.0,]
-g = G[(0.0)]
-h = H[(0.0)]
+bx =   [ 1.; 1.; 1.;; ]
+∇bx = G[ 0.; 0.; 0.;; ]
+Hbx = H[ 0.; 0.; 0.;; ]
+
+test_internals(order,x,bx,∇bx,Hbx)
 
-bx = repeat(permutedims(v),np)
-∇bx = repeat(permutedims(g),np)
-Hbx = repeat(permutedims(h),np)
 test_field_array(b,x,bx,grad=∇bx,gradgrad=Hbx)
 test_field_array(b,x[1],bx[1,:],grad=∇bx[1,:],gradgrad=Hbx[1,:])
 
@@ -55,6 +71,8 @@ Hbx = H[  0. 0.
           0. 0.
           0. 0. ]
 
+test_internals(order,x,bx,∇bx,Hbx)
+
 test_field_array(b,x,bx,grad=∇bx,gradgrad=Hbx)
 test_field_array(b,x[1],bx[1,:],grad=∇bx[1,:],gradgrad=Hbx[1,:])
 
@@ -77,6 +95,8 @@ Hbx = H[  2. -4. 2.
           2. -4. 2.
           2. -4. 2. ]
 
+test_internals(order,x,bx,∇bx,Hbx)
+
 test_field_array(b,x,bx,≈, grad=∇bx,gradgrad=Hbx)
 test_field_array(b,x[1],bx[1,:],grad=∇bx[1,:],gradgrad=Hbx[1,:])
 
@@ -103,6 +123,8 @@ bx  = [      bernstein(order,n)( xi[1])  for xi in x,  n in 0:order]
 #    -6x+6 18x-12 -18x+6  6x
 Hbx = [ H(_H(bernstein(order,n))(xi[1])) for xi in x,  n in 0:order]
 
+test_internals(order,x,bx,∇bx,Hbx)
+
 test_field_array(b,x,bx,≈, grad=∇bx,gradgrad=Hbx)
 test_field_array(b,x[1],bx[1,:],grad=∇bx[1,:],gradgrad=Hbx[1,:])
 
@@ -116,6 +138,8 @@ bx  = [      bernstein(order,n)( xi[1])  for xi in x,  n in 0:order]
 ∇bx = [ G(_∇(bernstein(order,n))(xi[1])) for xi in x,  n in 0:order]
 Hbx = [ H(_H(bernstein(order,n))(xi[1])) for xi in x,  n in 0:order]
 
+test_internals(order,x,bx,∇bx,Hbx)
+
 test_field_array(b,x,bx,≈, grad=∇bx,gradgrad=Hbx)
 test_field_array(b,x[1],bx[1,:],grad=∇bx[1,:],gradgrad=Hbx[1,:])