Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Faster _changebasis_sim #242

Merged
merged 4 commits into from
Aug 14, 2022
Merged

Faster _changebasis_sim #242

merged 4 commits into from
Aug 14, 2022

Conversation

hyrodium
Copy link
Owner

@hyrodium hyrodium commented Aug 1, 2022

using BasicBSpline
using BenchmarkTools

p = 3
P1 = BSplineSpace{p}(KnotVector(1:8))
P2 = BSplineSpace{p}(KnotVector(2,3,3,4,5,6,7,8))
changebasis(P1,P2)
@benchmark BasicBSpline._changebasis_sim(P1,P2)

On the current master

BenchmarkTools.Trial: 10000 samples with 1 evaluation.
 Range (min … max):  20.158 μs …  6.812 ms  ┊ GC (min … max): 0.00% … 98.70%
 Time  (median):     22.092 μs              ┊ GC (median):    0.00%
 Time  (mean ± σ):   23.798 μs ± 95.967 μs  ┊ GC (mean ± σ):  5.65% ±  1.40%

          ▂▅▆█▇▅▅▄▃                                            
  ▁▁▁▂▃▄▇████████████▇▆▆▄▅▄▄▃▃▃▃▂▂▂▂▂▂▂▂▂▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁ ▃
  20.2 μs         Histogram: frequency by time        28.3 μs <

 Memory estimate: 8.38 KiB, allocs estimate: 173.

With this PR

BenchmarkTools.Trial: 10000 samples with 1 evaluation.
 Range (min … max):  14.698 μs … 85.201 μs  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     15.739 μs              ┊ GC (median):    0.00%
 Time  (mean ± σ):   16.729 μs ±  2.978 μs  ┊ GC (mean ± σ):  0.00% ± 0.00%

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

 Memory estimate: 7.53 KiB, allocs estimate: 167.

@hyrodium
Copy link
Owner Author

hyrodium commented Aug 1, 2022

Getting faster with bsplinebasisall.

julia> using BasicBSpline

julia> using BenchmarkTools

julia> p = 3
3

julia> P1 = BSplineSpace{p}(KnotVector(1:8))
BSplineSpace{3, Int64}(KnotVector([1, 2, 3, 4, 5, 6, 7, 8]))

julia> P2 = BSplineSpace{p}(KnotVector(2,3,3,4,5,6,7,8))
BSplineSpace{3, Int64}(KnotVector([2, 3, 3, 4, 5, 6, 7, 8]))

julia> changebasis(P1,P2)
4×4 Matrix{Float64}:
 0.666667     0.0  -4.93432e-17   0.0
 0.333333     1.0  -2.22045e-16  -5.55112e-17
 3.70074e-17  0.0   1.0          -2.22045e-16
 0.0          0.0   0.0           1.0

julia> @benchmark BasicBSpline._changebasis_sim(P1,P2)
BenchmarkTools.Trial: 10000 samples with 1 evaluation.
 Range (min  max):   9.919 μs   5.473 ms  ┊ GC (min  max): 0.00%  99.17%
 Time  (median):     10.450 μs              ┊ GC (median):    0.00%
 Time  (mean ± σ):   11.455 μs ± 54.631 μs  ┊ GC (mean ± σ):  4.74% ±  0.99%

   ▁▇█▅▁                                                       
  ▂█████▇▅▄▃▂▂▂▃▃▃▂▂▃▂▂▂▂▂▂▂▂▂▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁ ▂
  9.92 μs         Histogram: frequency by time        16.2 μs <

 Memory estimate: 6.78 KiB, allocs estimate: 143.

@hyrodium
Copy link
Owner Author

hyrodium commented Aug 1, 2022

23.798 μs → 1.319 μs 🎉 18x faster 🎉 🎉 🎉

julia> @benchmark BasicBSpline._changebasis_sim(P1,P2)
BenchmarkTools.Trial: 10000 samples with 10 evaluations.
 Range (min  max):  988.800 ns  376.767 μs  ┊ GC (min  max): 0.00%  99.21%
 Time  (median):       1.097 μs               ┊ GC (median):    0.00%
 Time  (mean ± σ):     1.319 μs ±   6.258 μs  ┊ GC (mean ± σ):  8.18% ±  1.72%

   ▆▇██▇▆▅▅▃▂▂▂▂▂▂▂▂▂▂▃▃▃▂▂▂▂▁▂▂▁▁▁▁                 ▁          ▂
  ███████████████████████████████████▇█▇█▇█▇▇▆▇▇▆▇▇▇▇███▇▇▇▄▆▅▄ █
  989 ns        Histogram: log(frequency) by time       2.26 μs <

 Memory estimate: 1.97 KiB, allocs estimate: 25.

@hyrodium hyrodium changed the title Faster _changebasis_sim with MMatrix Faster _changebasis_sim Aug 1, 2022
@hyrodium
Copy link
Owner Author

hyrodium commented Aug 1, 2022

The test failures are related to the following issues around StaticArrays.jl.

0×0 matrix division

julia> using StaticArrays

julia> n = SMatrix{0,0,Float64}()
0×0 SMatrix{0, 0, Float64, 0} with indices SOneTo(0)×SOneTo(0)

julia> n \ n
0×0 SMatrix{0, 0, Float64, 0} with indices SOneTo(0)×SOneTo(0)

julia> n / n
ERROR: MethodError: one(::Type{Union{}}) is ambiguous. Candidates:
  one(::Union{Type{T}, T}) where T<:AbstractString in Base at strings/basic.jl:262
  one(::Union{Type{P}, P}) where P<:Dates.Period in Dates at /home/hyrodium/julia/julia-1.7.3/share/julia/stdlib/v1.7/Dates/src/periods.jl:54
  one(::Type{SM}) where SM<:(Union{LinearAlgebra.Adjoint{T, <:Union{StaticArray{Tuple{var"#s2"}, T, 1} where var"#s2", StaticArray{Tuple{var"#s3", var"#s4"}, T, 2} where {var"#s3", var"#s4"}}}, LinearAlgebra.Diagonal{T, <:StaticArray{Tuple{var"#s13"}, T, 1} where var"#s13"}, LinearAlgebra.Hermitian{T, <:StaticArray{Tuple{var"#s10", var"#s11"}, T, 2} where {var"#s10", var"#s11"}}, LinearAlgebra.LowerTriangular{T, <:StaticArray{Tuple{var"#s18", var"#s19"}, T, 2} where {var"#s18", var"#s19"}}, LinearAlgebra.Symmetric{T, <:StaticArray{Tuple{var"#s7", var"#s8"}, T, 2} where {var"#s7", var"#s8"}}, LinearAlgebra.Transpose{T, <:Union{StaticArray{Tuple{var"#s2"}, T, 1} where var"#s2", StaticArray{Tuple{var"#s3", var"#s4"}, T, 2} where {var"#s3", var"#s4"}}}, LinearAlgebra.UnitLowerTriangular{T, <:StaticArray{Tuple{var"#s24", var"#s25"}, T, 2} where {var"#s24", var"#s25"}}, LinearAlgebra.UnitUpperTriangular{T, <:StaticArray{Tuple{var"#s21", var"#s22"}, T, 2} where {var"#s21", var"#s22"}}, LinearAlgebra.UpperTriangular{T, <:StaticArray{Tuple{var"#s15", var"#s16"}, T, 2} where {var"#s15", var"#s16"}}, StaticArray{Tuple{var"#s1", var"#s3"}, T, 2} where {var"#s1", var"#s3"}} where T) in StaticArrays at /home/hyrodium/.julia/dev/StaticArrays/src/linalg.jl:108
  one(::Type{<:AbstractIrrational}) in Base at irrationals.jl:154
  one(::Type{T}) where T<:Number in Base at number.jl:334
Possible fix, define
  one(::Type{Union{}})
Stacktrace:
  [1] __lu(A::SMatrix{0, 0, Union{}, 0}, #unused#::Val{true})
    @ StaticArrays ~/.julia/dev/StaticArrays/src/lu.jl:112
  [2] macro expansion
    @ ~/.julia/dev/StaticArrays/src/lu.jl:85 [inlined]
  [3] _lu(A::SMatrix{0, 0, Union{}, 0}, pivot::Val{true}, check::Bool)
    @ StaticArrays ~/.julia/dev/StaticArrays/src/lu.jl:77
  [4] lu(A::SMatrix{0, 0, Union{}, 0}, pivot::Val{true}; check::Bool)
    @ StaticArrays ~/.julia/dev/StaticArrays/src/lu.jl:50
  [5] lu(A::SMatrix{0, 0, Union{}, 0}; check::Bool)
    @ StaticArrays ~/.julia/dev/StaticArrays/src/lu.jl:54
  [6] lu(A::SMatrix{0, 0, Union{}, 0})
    @ StaticArrays ~/.julia/dev/StaticArrays/src/lu.jl:54
  [7] macro expansion
    @ ~/.julia/dev/StaticArrays/src/solve.jl:55 [inlined]
  [8] _solve
    @ ~/.julia/dev/StaticArrays/src/solve.jl:45 [inlined]
  [9] \
    @ ~/.julia/dev/StaticArrays/src/solve.jl:1 [inlined]
 [10] /(A::SMatrix{0, 0, Float64, 0}, B::SMatrix{0, 0, Float64, 0})
    @ LinearAlgebra ~/julia/julia-1.7.3/share/julia/stdlib/v1.7/LinearAlgebra/src/generic.jl:1152
 [11] top-level scope
    @ REPL[4]:1

This is already reported in JuliaArrays/StaticArrays.jl#1066.

Matrix division with SMatrix{n, n, BigFloat}

julia> m = @SMatrix rand(BigFloat,3,3)
3×3 SMatrix{3, 3, BigFloat, 9} with indices SOneTo(3)×SOneTo(3):
 0.0533287  0.646633  0.713034
 0.203678   0.225494  0.889687
 0.560856   0.266583  0.283933

julia> m * m
3×3 SMatrix{3, 3, BigFloat, 9} with indices SOneTo(3)×SOneTo(3):
 0.534458  0.370379  0.815781
 0.555776  0.419728  0.598461
 0.243453  0.498473  0.717703

julia> m * inv(m)
3×3 SMatrix{3, 3, BigFloat, 9} with indices SOneTo(3)×SOneTo(3):
  1.0           0.0          4.31808e-78
  0.0           1.0          0.0
 -3.23856e-78  -4.31808e-78  1.0

julia> m / m
ERROR: setindex!() with non-isbitstype eltype is not supported by StaticArrays. Consider using SizedArray.
Stacktrace:
  [1] error(s::String)
    @ Base ./error.jl:33
  [2] setindex!
    @ ~/.julia/dev/StaticArrays/src/MArray.jl:39 [inlined]
  [3] macro expansion
    @ ~/.julia/dev/StaticArrays/src/indexing.jl:66 [inlined]
  [4] _setindex!_scalar
    @ ~/.julia/dev/StaticArrays/src/indexing.jl:46 [inlined]
  [5] setindex!
    @ ~/.julia/dev/StaticArrays/src/indexing.jl:42 [inlined]
  [6] macro expansion
    @ ~/.julia/dev/StaticArrays/src/indexing.jl:357 [inlined]
  [7] _setindex!
    @ ~/.julia/dev/StaticArrays/src/indexing.jl:320 [inlined]
  [8] setindex!
    @ ~/.julia/dev/StaticArrays/src/indexing.jl:274 [inlined]
  [9] _solve
    @ ~/.julia/dev/StaticArrays/src/solve.jl:36 [inlined]
 [10] \
    @ ~/.julia/dev/StaticArrays/src/solve.jl:1 [inlined]
 [11] /(A::SMatrix{3, 3, BigFloat, 9}, B::SMatrix{3, 3, BigFloat, 9})
    @ LinearAlgebra ~/julia/julia-1.7.3/share/julia/stdlib/v1.7/LinearAlgebra/src/generic.jl:1152
 [12] top-level scope
    @ REPL[19]:1

I'm not sure this should be fixed in StaticArrays.jl. (The error message can be updated, I guess.)

@hyrodium
Copy link
Owner Author

hyrodium commented Aug 4, 2022

I'm not sure this should be fixed in StaticArrays.jl. (The error message can be updated, I guess.)

Maybe, I'll fix this in a few weeks. (x-ref: JuliaArrays/StaticArrays.jl#1069)

@codecov
Copy link

codecov bot commented Aug 14, 2022

Codecov Report

Merging #242 (acd18a0) into master (fae9396) will increase coverage by 0.01%.
The diff coverage is 100.00%.

@@            Coverage Diff             @@
##           master     #242      +/-   ##
==========================================
+ Coverage   98.68%   98.70%   +0.01%     
==========================================
  Files          15       15              
  Lines        1372     1388      +16     
==========================================
+ Hits         1354     1370      +16     
  Misses         18       18              
Impacted Files Coverage Δ
src/_ChangeBasis.jl 100.00% <100.00%> (ø)

Help us with your feedback. Take ten seconds to tell us how you rate us. Have a feature suggestion? Share it here.

@hyrodium hyrodium merged commit 9d0049d into master Aug 14, 2022
@hyrodium hyrodium deleted the feature/faster_changebasis_sim branch August 22, 2022 11:58
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

1 participant