Update refinement

[hyrodium]

No description provided.

[codecov]

Codecov Report

Attention: 26 lines in your changes are missing coverage. Please review.

Comparison is base (3899682) 93.79% compared to head (e041b4a) 94.90%.

Files	Patch %	Lines
src/_Refinement.jl	87.67%	26 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##             main     #357      +/-   ##
==========================================
+ Coverage   93.79%   94.90%   +1.10%     
==========================================
  Files          14       14              
  Lines        1660     1786     +126     
==========================================
+ Hits         1557     1695     +138     
+ Misses        103       91      -12     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

[hyrodium]

I thought the following function my_refinement_I may fix the type inference, but the performance was much worse.

julia> using BasicBSpline

julia> using BenchmarkTools

julia> using Test

julia> using GeometryBasics

julia> using SparseArrays

julia> function _a(value::T, index::Integer, dims::NTuple{Dim, Integer}) where {T, Dim}
           a = Array{T,Dim}(undef, dims)
           i = prod(index)
           for i in 1:(i-1)
               @inbounds a[i] = zero(value)
           end
           @inbounds a[i] = value
           return a
       end
_a (generic function with 1 method)

julia> function _isempty(R::NTuple{2, UnitRange{Int}}, J::CartesianIndex{2})
           return isempty(R[1][J[1]]) || isempty(R[2][J[2]])
       end
_isempty (generic function with 1 method)

julia> @generated function my_refinement_I(M::BSplineManifold{Dim}, P′::NTuple{Dim, BSplineSpace}) where {Dim}
           P = [Symbol(:P,i) for i in 1:Dim]
           P′ = [Symbol(:P,i,:(var"′")) for i in 1:Dim]
           A = [Symbol(:A,i) for i in 1:Dim]
           R = [Symbol(:R,i) for i in 1:Dim]
           i = [Symbol(:i,i) for i in 1:Dim]
           j = [Symbol(:j,i) for i in 1:Dim]
           n′ = [Symbol(:n,i,:(var"′")) for i in 1:Dim]
           Expr(
               :block,
               :($(Expr(:tuple, P...)) = bsplinespaces(M)),
               [:($(A[_i]) = changebasis_I($(P[_i]), $(P′[_i]))) for _i in 1:Dim]...,
               [:($(n′[_i]) = dim($(P′[_i]))) for _i in 1:Dim]...,
               [:($(R[_i]) = BasicBSpline._i_ranges_I($(A[_i]), $(P′[_i]))) for _i in 1:Dim]...,
               [:($(j[_i]) = findfirst(!isempty, $(R[_i]))) for _i in 1:Dim]...,
               (
                   ex = Expr(:tuple, [:($(R[_i])[$(j[_i])]) for _i in 1:Dim]...);
                   :(C = CartesianIndices($(ex)))
               ),
               (
                   ex = Expr(:ref, :a, [:(first($(R[_i])[$(j[_i])])) for _i in 1:Dim]...);
                   ex = Expr(:call, :*, [:($(A[_i])[first($(R[_i])[$(j[_i])]),$(j[_i])]) for _i in 1:Dim]..., ex);
                   :(value = $ex)
               ),
               (
                   ex = Expr(:ref, :a, [:($(i[_i])) for _i in 1:Dim]...);
                   Expr(:for, :(ci1 = C[2:end]),
                       Expr(
                           :block,
                           (
                               :($(Expr(:tuple, i...)) = ci1.I)
                           ),
                           (
                               ex = Expr(:call, :*, [:($(A[_i])[$(i[_i]),$(j[_i])]) for _i in 1:Dim]..., ex);
                               :(value += $(ex))
                           ),
                       )
                   ),
               ),
               :(j = $(Expr(:call, :*, j...))),
               :(a′ = _a(value, j, $(Expr(:tuple, n′...)))),
               (
                   ex = Expr(:tuple, [:(1:$(n′[_i])) for _i in 1:Dim]...);
                   Expr(:for, :(cj = CartesianIndices($ex)),
                       Expr(
                           :block,
                           :($(Expr(:tuple, j...)) = cj.I),
                           Expr(:if, :(isempty(R1[j1]) || isempty(R2[j2]) ),
                               (
                                   ex = Expr(:ref, :a′, [:($(j[_i])) for _i in 1:Dim]...);
                                   :($ex = zero(value))
                               ),
                               Expr(
                                   :block,
                                   (
                                       ex = Expr(:tuple, [:($(R[_i])[$(j[_i])]) for _i in 1:Dim]...);
                                       :(C = CartesianIndices($(ex)))
                                   ),
                                   (
                                       ex = Expr(:ref, :a, [:(first($(R[_i])[$(j[_i])])) for _i in 1:Dim]...);
                                       ex = (Expr(:call, :*, [:($(A[_i])[first($(R[_i])[$(j[_i])]),$(j[_i])]) for _i in 1:Dim]..., ex));
                                       :(value = $ex)
                                   ),
                                   (
                                       ex = Expr(:ref, :a, [:($(i[_i])) for _i in 1:Dim]...);
                                       Expr(:for, :(ci2 = C[2:end]),
                                           Expr(:block,
                                               :($(Expr(:tuple, i...)) = ci2.I),
                                               (
                                                   ex = Expr(:call, :*, [:($(A[_i])[$(i[_i]),$(j[_i])]) for _i in 1:Dim]..., ex);
                                                   :(value += $ex)
                                               )
                                           )
                                       ),
                                   ),
                                   (
                                       ex = Expr(:ref, :a′, [:($(j[_i])) for _i in 1:Dim]...);
                                       :($ex = value)
                                   )
                               )
                           )
                       )
                   )
               ),
               :(return BSplineManifold(a′, P′))
           )
       end
my_refinement_I (generic function with 1 method)

julia> P1 = BSplineSpace{2}(knotvector"1112121113")
BSplineSpace{2, Int64, KnotVector{Int64}}(KnotVector([1, 2, 3, 4, 4, 5, 6, 6, 7, 8, 9, 10, 10, 10]))

julia> P2 = BSplineSpace{1}(knotvector"11111 12")
BSplineSpace{1, Int64, KnotVector{Int64}}(KnotVector([1, 2, 3, 4, 5, 7, 8, 8]))

julia> P1′ = expandspace_I(P1, Val(1), KnotVector([3,3,3]))
BSplineSpace{3, Int64, KnotVector{Int64}}(KnotVector([1, 2, 3, 3, 3, 3, 3, 4, 4, 4, 5, 5, 6, 6, 6, 7, 7, 8, 8, 9, 9, 10, 10, 10, 10]))

julia> P2′ = expandspace_I(P2, Val(1), KnotVector([4,4,4]))
BSplineSpace{2, Int64, KnotVector{Int64}}(KnotVector([1, 2, 2, 3, 3, 4, 4, 4, 4, 4, 5, 5, 7, 7, 8, 8, 8]))

julia> n1 = dim(P1)
11

julia> n2 = dim(P2)
6

julia> a = rand(n1, n2);

julia> M = BSplineManifold(a, P1, P2);

julia> my_refinement_I(M, (P1′, P2′))
BSplineManifold{2, (3, 2), Float64, Int64, Tuple{BSplineSpace{3, Int64, KnotVector{Int64}}, BSplineSpace{2, Int64, KnotVector{Int64}}}}([0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.3595498013883678 0.4560889148520694 … 0.2575447846245697 0.24031646689580166; 0.24550313124722112 0.6203343192044695 … 0.12398973608093167 0.13267267832326124], (BSplineSpace{3, Int64, KnotVector{Int64}}(KnotVector([1, 2, 3, 3, 3, 3, 3, 4, 4, 4, 5, 5, 6, 6, 6, 7, 7, 8, 8, 9, 9, 10, 10, 10, 10])), BSplineSpace{2, Int64, KnotVector{Int64}}(KnotVector([1, 2, 2, 3, 3, 4, 4, 4, 4, 4, 5, 5, 7, 7, 8, 8, 8]))))

julia> refinement_I(M, (P1′, P2′)) == my_refinement_I(M, (P1′, P2′))
true

julia> @benchmark refinement_I(M, (P1′, P2′))
BenchmarkTools.Trial: 10000 samples with 1 evaluation.
 Range (min … max):  39.003 μs …  4.367 ms  ┊ GC (min … max): 0.00% … 91.28%
 Time  (median):     49.294 μs              ┊ GC (median):    0.00%
 Time  (mean ± σ):   52.387 μs ± 93.484 μs  ┊ GC (mean ± σ):  5.04% ±  2.84%

  ▁▆██▆▃▁▁▁▂▆█▇▆▄▃▃▂▂▂▂▂▁▁▁ ▁                                 ▂
  █████████████████████████▇█████▇▇█▇█▆██▇▇▇█▇█▆▇▆▇█▆▆▆▆▅▅▅▆▅ █
  39 μs        Histogram: log(frequency) by time      94.4 μs <

 Memory estimate: 47.36 KiB, allocs estimate: 374.

julia> @benchmark my_refinement_I(M, (P1′, P2′))
BenchmarkTools.Trial: 10000 samples with 1 evaluation.
 Range (min … max):  426.187 μs …   5.635 ms  ┊ GC (min … max): 0.00% … 79.62%
 Time  (median):     440.614 μs               ┊ GC (median):    0.00%
 Time  (mean ± σ):   486.117 μs ± 333.865 μs  ┊ GC (mean ± σ):  5.83% ±  7.62%

  ▅█▇▄▂▂▂▂▁ ▃▄▃▂▂▁                                              ▂
  ██████████████████▇▇▇▇▇▇▆▆▇▆▆▆▅▃▁▁▄▄▄▃▁▃▃▃▃▁▃▁▁▄▁▃▃▃▄▁▃▁▁▁▁▄▃ █
  426 μs        Histogram: log(frequency) by time        886 μs <

 Memory estimate: 346.58 KiB, allocs estimate: 10276.

[hyrodium]

It seems we can avoid @generated macro here. (cad4900)

[hyrodium]

Now the things that I would like to fix are finished.

TODOs for tomorrow:

• Add more test cases
• Check performance again
• Fix refinement for RationalBSplineManifold

x-ref: #354 (comment)