Repository: hyrodium/BasicBSpline.jl
Pull request #156 (2022-01-30 06:10:14)

`innerproduct_R(::UniformBSplineSpace)` with Euler's triangle
This PR solves #154.

```julia
julia> using BasicBSpline, BenchmarkTools

julia> k = UniformKnotVector(1:34)
UniformKnotVector(1:34)

julia> P1 = UniformBSplineSpace{4}(k)
UniformBSplineSpace{4, Int64, UnitRange{Int64}}(UniformKnotVector(1:34))

julia> P2 = BSplineSpace(P1)
BSplineSpace{4, Int64}(KnotVector([1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34]))

julia> BasicBSpline.innerproduct_R(P1) ≈ BasicBSpline.innerproduct_R(P2)
true

julia> @benchmark BasicBSpline.innerproduct_R(P1)
BenchmarkTools.Trial: 10000 samples with 42 evaluations.
 Range (min … max):  909.119 ns … 37.527 μs  ┊ GC (min … max):  0.00% … 95.52%
 Time  (median):       1.089 μs              ┊ GC (median):     0.00%
 Time  (mean ± σ):     1.384 μs ±  2.335 μs  ┊ GC (mean ± σ):  12.88% ±  7.42%

     ▃▇█▇▆▄▃▄▄▅▃▂▂▁                        ▂▁▁▁                ▂
  ▅▅█████████████████▇▇▇▆▆▅▅▆▅▆▅▅▄▃▅▂▃▂▃▃▅███████▇▆▆▆▆▆▇▇▆▆▆▆▅ █
  909 ns        Histogram: log(frequency) by time      2.62 μs <

 Memory estimate: 6.78 KiB, allocs estimate: 2.

julia> @benchmark BasicBSpline.innerproduct_R(P2)
BenchmarkTools.Trial: 10000 samples with 1 evaluation.
 Range (min … max):  18.395 μs …  1.846 ms  ┊ GC (min … max): 0.00% … 98.57%
 Time  (median):     19.116 μs              ┊ GC (median):    0.00%
 Time  (mean ± σ):   19.668 μs ± 25.155 μs  ┊ GC (mean ± σ):  1.80% ±  1.39%

        ▁▃▅▆█▅▄▂                                               
  ▁▁▁▂▃▆████████▇▇▆▄▃▃▃▃▃▃▃▃▃▃▃▂▂▂▂▂▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁ ▂
  18.4 μs         Histogram: frequency by time        22.1 μs <

 Memory estimate: 8.28 KiB, allocs estimate: 12.
```


Patch:
diff --git a/src/_Fitting.jl b/src/_Fitting.jl
index 6cde57e2..6e7ba923 100644
--- a/src/_Fitting.jl
+++ b/src/_Fitting.jl
@@ -6,7 +6,7 @@ Calculate a matrix
 A_{ij}=\int_{I} B_{(i,p,k)}(t) B_{(j,p,k)}(t) dt
 ```
 """
-function innerproduct_I(P::BSplineSpace{p}) where p
+function innerproduct_I(P::AbstractBSplineSpace{p}) where p
     n = dim(P)
     A = zeros(n,n)
     k = knotvector(P)
@@ -38,7 +38,7 @@ Calculate a matrix
 A_{ij}=\int_{\mathbb{R}} B_{(i,p,k)}(t) B_{(j,p,k)}(t) dt
 ```
 """
-function innerproduct_R(P::BSplineSpace{p}) where p
+function innerproduct_R(P::AbstractBSplineSpace{p}) where p
     n = dim(P)
     A = innerproduct_I(P).data
     k = knotvector(P)
@@ -71,7 +71,7 @@ function innerproduct_R(P::BSplineSpace{p}) where p
     return Symmetric(A)
 end
 
-function innerproduct_I(func, Ps::Tuple{<:BSplineSpace{p₁}}) where {p₁}
+function innerproduct_I(func, Ps::Tuple{<:AbstractBSplineSpace{p₁}}) where {p₁}
     P₁, = Ps
     n₁ = dim(P₁)
     k₁ = knotvector(P₁)
@@ -98,7 +98,7 @@ function innerproduct_I(func, Ps::Tuple{<:BSplineSpace{p₁}}) where {p₁}
     return b
 end
 
-function innerproduct_I(func, Ps::Tuple{<:BSplineSpace{p₁},<:BSplineSpace{p₂}}) where {p₁,p₂}
+function innerproduct_I(func, Ps::Tuple{<:AbstractBSplineSpace{p₁},<:AbstractBSplineSpace{p₂}}) where {p₁,p₂}
     P₁,P₂ = Ps
     n₁,n₂ = dim.(Ps)
     k₁,k₂ = knotvector.(Ps)
@@ -137,7 +137,7 @@ function innerproduct_I(func, Ps::Tuple{<:BSplineSpace{p₁},<:BSplineSpace{p₂
     return b
 end
 
-function innerproduct_I(func, Ps::Tuple{<:BSplineSpace{p₁},<:BSplineSpace{p₂},<:BSplineSpace{p₃}}) where {p₁,p₂,p₃}
+function innerproduct_I(func, Ps::Tuple{<:AbstractBSplineSpace{p₁},<:AbstractBSplineSpace{p₂},<:AbstractBSplineSpace{p₃}}) where {p₁,p₂,p₃}
     P₁,P₂,P₃ = Ps
     n₁,n₂,n₃ = dim.(Ps)
     k₁,k₂,k₃ = knotvector.(Ps)
@@ -191,7 +191,7 @@ function innerproduct_I(func, Ps::Tuple{<:BSplineSpace{p₁},<:BSplineSpace{p₂
     return b
 end
 
-function innerproduct_R(func, Ps::Tuple{<:BSplineSpace{p₁}}) where {p₁}
+function innerproduct_R(func, Ps::Tuple{<:AbstractBSplineSpace{p₁}}) where {p₁}
     P₁, = Ps
     n₁ = dim(P₁)
     k₁ = knotvector(P₁)
@@ -214,7 +214,7 @@ function innerproduct_R(func, Ps::Tuple{<:BSplineSpace{p₁}}) where {p₁}
     return b
 end
 
-function innerproduct_R(func, Ps::Tuple{<:BSplineSpace{p₁},<:BSplineSpace{p₂}}) where {p₁,p₂}
+function innerproduct_R(func, Ps::Tuple{<:AbstractBSplineSpace{p₁},<:AbstractBSplineSpace{p₂}}) where {p₁,p₂}
     P₁,P₂ = Ps
     n₁,n₂ = dim.(Ps)
     k₁,k₂ = knotvector.(Ps)
@@ -246,7 +246,7 @@ function innerproduct_R(func, Ps::Tuple{<:BSplineSpace{p₁},<:BSplineSpace{p₂
     return b
 end
 
-function innerproduct_R(func, Ps::Tuple{<:BSplineSpace{p₁},<:BSplineSpace{p₂},<:BSplineSpace{p₃}}) where {p₁,p₂,p₃}
+function innerproduct_R(func, Ps::Tuple{<:AbstractBSplineSpace{p₁},<:AbstractBSplineSpace{p₂},<:AbstractBSplineSpace{p₃}}) where {p₁,p₂,p₃}
     P₁,P₂,P₃ = Ps
     n₁,n₂,n₃ = dim.(Ps)
     k₁,k₂,k₃ = knotvector.(Ps)
@@ -290,6 +290,38 @@ function innerproduct_R(func, Ps::Tuple{<:BSplineSpace{p₁},<:BSplineSpace{p₂
     return b
 end
 
+function innerproduct_R(P::UniformBSplineSpace{p,T,<:AbstractUnitRange}) where {p,T}
+    U = StaticArrays.arithmetic_closure(T)
+    n = dim(P)
+    A = zeros(U,n,n)
+    m = factorial(2p+1)
+
+    for q in 0:p
+        a = U(eulertriangle(2p+1,q))/m
+        for i in 1:n-(p-q)
+            A[i,i+(p-q)] = a
+        end
+    end
+    return Symmetric(A)
+end
+
+function innerproduct_R(P::UniformBSplineSpace{p,T}) where {p,T}
+    U = StaticArrays.arithmetic_closure(T)
+    d = step(_vec(knotvector(P)))
+    n = dim(P)
+    A = zeros(U,n,n)
+    m = factorial(2p+1)
+
+    for q in 0:p
+        a = U(eulertriangle(2p+1,q)*d)/m
+        for i in 1:n-(p-q)
+            A[i,i+(p-q)] = a
+        end
+    end
+    return Symmetric(A)
+end
+
+
 """
 * func: Real -> ℝ-vector space
 """
diff --git a/src/_util.jl b/src/_util.jl
index 0e325023..82b6fda4 100644
--- a/src/_util.jl
+++ b/src/_util.jl
@@ -5,3 +5,6 @@
     return U(a/b)
 end
 @inline _d(a,b) = _d(promote(a,b)...)
+
+# Euler's triangle number (http://oeis.org/wiki/Eulerian_numbers,_triangle_of)
+eulertriangle(n,k) = sum((-1)^j*binomial(n+1,j)*(k-j+1)^n for j in 0:k)


Test patch:
diff --git a/test/test_Fitting.jl b/test/test_Fitting.jl
index 42677318..aaeec8fa 100644
--- a/test/test_Fitting.jl
+++ b/test/test_Fitting.jl
@@ -174,4 +174,31 @@
         a_fit = fittingcontrolpoints(M, (P1′, P2′, P3′))
         @test norm(a_fit - a_ref) < 1.0e-6
     end
+
+    @testset "uniform" begin
+        k = UniformKnotVector(1:42)
+        for p in 0:3
+            P = UniformBSplineSpace{p}(k)
+            @test BasicBSpline.innerproduct_R(P) ≈ BasicBSpline.innerproduct_R(BSplineSpace(P))
+            @test BasicBSpline.innerproduct_R(P) isa Symmetric{Float64,Matrix{Float64}}
+            @test BasicBSpline.innerproduct_R(BSplineSpace(P)) isa Symmetric{Float64,Matrix{Float64}}
+        end
+
+        k = UniformKnotVector(1:2:42)
+        for p in 0:3
+            P = UniformBSplineSpace{p}(k)
+            @test BasicBSpline.innerproduct_R(P) ≈ BasicBSpline.innerproduct_R(BSplineSpace(P))
+            @test BasicBSpline.innerproduct_R(P) isa Symmetric{Float64,Matrix{Float64}}
+            @test BasicBSpline.innerproduct_R(BSplineSpace(P)) isa Symmetric{Float64,Matrix{Float64}}
+        end
+
+        k = UniformKnotVector(1:2:42//1)
+        for p in 0:3
+            P = UniformBSplineSpace{p}(k)
+            @test BasicBSpline.innerproduct_R(P) ≈ BasicBSpline.innerproduct_R(BSplineSpace(P))
+            @test BasicBSpline.innerproduct_R(P) isa Symmetric{Rational{Int},Matrix{Rational{Int}}}
+            @test BasicBSpline.innerproduct_I(P) isa Symmetric{Float64,Matrix{Float64}}  # TODO: This should be Rational
+            @test BasicBSpline.innerproduct_R(BSplineSpace(P)) isa Symmetric{Float64,Matrix{Float64}}
+        end
+    end
 end