新しくなったBasicBSpline.jl(v0.8.3)の紹介

本記事はJuliaLang Advent Calendar 2022の19日目の記事です!
遅くなってごめんなさい！

はじめに

一昨年にBasicBSpline.jlを作ったので宣伝です！を書きました。
当時はBasicBSpline.jlのv0.1.1をリリースしたばかりでしたが、現在のBasicBSpline.jlの最新バージョンはv0.8.3です。
主な変更点は以下です。

• ノット列の型KnotVectorがパラメトリック型になってFloat64以外にもRationalやBigFloatに対応できるようになった
• NURBSに対応した
• Plots.jlに対応した
• 標準のBSplineSpaceが高速かつ汎用的になった
• 等間隔ノット列の型(UniformKnotVector)を用意して、一部の演算が高速化された
• 関数・型の名称変更 (e.g. Knots→KnotVector, bsplineunity→domain)

以前の記事は古くなってしまったので、本記事では新しくなったBasicBSpline.jlでのコード例や実装時に気をつけたことなどを書いていこうと思います。

B-splineってなに？

区分多項式を便利に扱えるように、関数空間の基底を取りやすいように道具を整備したものがB-splineです。
これまでに、色々なB-splineの資料を作ってきました。

• B-spline入門（線形代数がすこし分かる人向け）
  • 動画で学びたい人はこちら。
  • 数理的なモチベーションが伝わりやすい…はず
• NURBS多様体による形状表現
  • 証明を詳細にフォローしたい人はこちら。
• BasicBSpline.jlのドキュメント
  • Juliaで実行しながらB-splineを勉強したい人はこちら。
  • 英語で書いてます。
• BasicBSpline.jlを作ったので宣伝です！
  • 前述の記事でBasicBSpline.jlの説明としては古いですが、B-splineのモチベーションの説明としては悪くないと思います。

私の作ったものとは違いますが、数日前に投稿された次の動画が視覚的に綺麗に可視化されていてかなり良かったです。

https://www.youtube.com/watch?v=jvPPXbo87ds

書籍ではGeometric Modeling with Splinesがオススメです。

BasicBSpline.jlのつかいかた

B-spline自体のモチベーションは既に上記資料で説明したので、以下ではコード例を紹介します。

形状表現

B-splineを使えば、ノット列と制御点と重みを使って形状を定義することができます。
たとえばトーラスは以下のようになります。

using BasicBSpline
using StaticArrays
using Plots
plotly()

# 半径の定義
R = 3
r = 1

# 制御点の定義
a0 = [
    SVector( 1, 0, 0),
    SVector( 1, 1, 0),
    SVector( 0, 1, 0),
    SVector(-1, 1, 0),
    SVector(-1, 0, 0),
    SVector(-1,-1, 0),
    SVector( 0,-1, 0),
    SVector( 1,-1, 0),
    SVector( 1, 0, 0)
]
a1 = (R+r)*a0
a5 = (R-r)*a0
a2 = [p+r*SVector(0,0,1) for p in a1]
a3 = [p+r*SVector(0,0,1) for p in R*a0]
a4 = [p+r*SVector(0,0,1) for p in a5]
a6 = [p-r*SVector(0,0,1) for p in a5]
a7 = [p-r*SVector(0,0,1) for p in R*a0]
a8 = [p-r*SVector(0,0,1) for p in a1]
a9 = a1
a = hcat(a1,a2,a3,a4,a5,a6,a7,a8,a9)

# 重みの定義
w = [1,1/√2,1,1/√2,1,1/√2,1,1/√2,1]

# B-spline空間の定義
k = KnotVector([0,0,0,1,1,2,2,3,3,4,4,4])
P = BSplineSpace{2}(k)

# B-spline多様体の定義
M = RationalBSplineManifold(a, w*w', P, P)
plot(M)

BasicBSpline.jlのドキュメントの説明もどうぞ！

内挿補間

区分多項式による補間にもB-splineを使うことができます。

using BasicBSpline
using StaticArrays
using Plots

function interpolate(xs::AbstractVector, fs::AbstractVector{T}) where T
    # Cubic open B-spline space
    p = 3
    k = KnotVector(xs) + KnotVector([xs[1],xs[end]]) * p
    P = BSplineSpace{p}(k)

    # dimensions
    m = length(xs)
    n = dim(P)

    # The interpolant function has a f''=0 property at bounds.
    ddP = BSplineDerivativeSpace{2}(P)
    dda = [bsplinebasis(ddP,j,xs[1]) for j in 1:n]
    ddb = [bsplinebasis(ddP,j,xs[m]) for j in 1:n]

    # Compute the interpolant function (1-dim B-spline manifold)
    M = [bsplinebasis(P,j,xs[i]) for i in 1:m, j in 1:n]
    M = vcat(dda', M, ddb')
    y = vcat(zero(T), fs, zero(T))
    return BSplineManifold(M\y, P)
end

# Example inputs
xs = [1, 2, 3, 4, 6, 7]
fs = [1.3, 1.5, 2, 2.1, 1.9, 1.3]
f = interpolate(xs,fs)

# Plot
scatter(xs, fs)
plot!(t->f(t))

BasicBSpline.jlのドキュメントの説明もどうぞ！

関数フィッティング

関数にフィッティングするためのfittingcontrolpoints関数が用意されています。
以下のコードはsinに対して-6..6区間でフィッティングする例です。

using BasicBSpline
using StaticArrays
using Plots

f(t) = SVector(t,sin(t)) # 正弦波のパラメータ表示
t0,t1 = -6,6             # 左端と右端

p = 3
k = KnotVector(t0:t1)+p*KnotVector([t0,t1])
P = BSplineSpace{p}(k)

a = fittingcontrolpoints(f,P)  # B-spline曲線の制御点の計算
M = BSplineManifold(a,P)
plot(sin,-10,10)
plot!(M)

なめらかな曲線をJuliaでSVG出力するもどうぞ！

BasicBSpline.jlの内部実装

以降の本記事では、BasicBSpline.jlの実装においてJuliaの言語機能がどのように便利だったかをお伝えします！

かなり雑多な内容になりますが、Juliaの便利さ・面白さを少しでもお伝えできれば嬉しいです。
B-splineの知識があればより楽しめますが、B-splineから独立したJuliaの一般の話はこのように「！枠」で囲んであるので斜め読みでもどうぞ。

ノット列

実軸上に乗った広義単調増大の有限列が、ノット列と呼ばれるものです。
B-splineは区分多項式を扱う道具であり、ノットに囲まれた各区間が一つの多項式に対応します。

抽象型と型変換

Juliaの型システムは木構造^[1]を成しており、根がAnyです。
BasicBSplineにおけるノット列を表す抽象型がAbstractKnotVectorで、以下の3つの部分型を持ちます。

julia> subtypes(AbstractKnotVector)
3-element Vector{Any}:
 EmptyKnotVector{T} where T<:Real
 KnotVector{T} where T<:Real
 UniformKnotVector{T} where T<:Real

これらのノット列は以下のようにインスタンスを定義できます。

julia> k1 = KnotVector([1, 4, 2.])  # コンストラクタが引数をsortする
KnotVector([1.0, 2.0, 4.0])

julia> k2 = KnotVector(1:5)
KnotVector([1, 2, 3, 4, 5])

julia> k3 = UniformKnotVector(1:5)  # 等間隔なノット列を扱うための型
UniformKnotVector(1:5)

julia> k4 = EmptyKnotVector()  # 空のノット列を扱うための型
EmptyKnotVector{Bool}()

julia> k5 = KnotVector(Float64[])
KnotVector(Float64[])

異なる型どうしでも等号比較できます。

julia> k1 == k2  # KnotVector{Float64}とKnotVector{Int}
false

julia> k2 == k3  # KnotVector{Int}とUniformKnotVector
true

julia> k4 == k5  # EmptyKnotVectorとKnotVector
true

型変換も実装されています。

julia> convert(KnotVector{Int}, k3)
KnotVector([1, 2, 3, 4, 5])

julia> KnotVector(k3)  # convertとほぼ同じ
KnotVector([1, 2, 3, 4, 5])

julia> KnotVector(k4)
KnotVector(Bool[])

julia> convert(UniformKnotVector{Int}, k2)  # k2が等間隔だとしてもUniformKnotVectorへは変換できない
ERROR: MethodError: no method matching (UniformKnotVector{Int64})(::KnotVector{Int64})

ある2つの型を同じ抽象型の部分型にするかどうかは

• 相互に変換できる型
• 同一の関数の異なるメソッドの引数として指定される型

などを目安に決めると良い。型の包含関係は、数学的な集合としての包含関係とは異なることに注意。

効率的なコンストラクタ

KnotVector{T}型はノット列を内部的にVector{T}として保持していますが、標準的なコンストラクタでは引数がコピー・ソートされるようになっています。

julia> v1 = [2,4,3]
3-element Vector{Int64}:
 2
 4
 3

julia> k1 = KnotVector(v1)
KnotVector([2, 3, 4])

julia> v1[1] = 3  # v1の要素を書き換えても
3

julia> k1  # k1は変わらない
KnotVector([2, 3, 4])

しかし以下ケースではデフォルトのKnotVectorのコンストラクタには不満があります。

• パフォーマンスのために配列のコピーを避けたい
• 引数の配列がソートされている場合にデフォルトのソートを避けたい

このようなケースに対応するため、unsafe_knotvector関数が用意されています。

julia> v2 = [2,4,3]
3-element Vector{Int64}:
 2
 4
 3

julia> k2 = BasicBSpline.unsafe_knotvector(Int, v2)  # sortしない
KnotVector([2, 4, 3])

julia> v2[1] = 5  # v2の値を書き換えると…
5

julia> k2
KnotVector([5, 4, 3])  # k2も変わる。メモリを共有しているため

KnotVectorのインスタンスを生成するにはKnotVectorのコンストラクタを経由する必要があるので、デフォルトのKnotVectorコンストラクタの処理を避けるのは困難なように思えます。
このような動作を実現するには、KnotVectorのstructの内側にglobalつきで関数を定義することで実装できます。

struct KnotVector{T} <: AbstractKnotVector{T}
    vector::Vector{T}
    global unsafe_knotvector(::Type{T}, v) where T = new{T}(v)
end

Base.Rationalに対しても同様のunsafe_rationalが用意されており、同じようなstructの内側にglobalが入るように実装されています。

julia> Rational(2,4)
1//2

julia> Base.unsafe_rational(2,4)
2//4

デフォルトのコンストラクタを避けれるような実装にするには、structの内側にglobalつきでコンストラクタを定義すれば良い。

加法

BasicBSpline.jlでは、以下のようにノット列(KnotVector)の足し算(+)ができます。^[2]

julia> k1 = KnotVector([1, 2, 3])
KnotVector([1, 2, 3])

julia> k2 = KnotVector([2, 4.5])
KnotVector([2.0, 4.5])

julia> k1+k2
KnotVector([1.0, 2.0, 2.0, 3.0, 4.5])

julia> typeof(k1), typeof(k2)
(KnotVector{Int64}, KnotVector{Float64})

この+の操作はBase.:+関数にKnotVector用のメソッドを追加して実装されています。
Juliaの多重ディスパッチの仕組みがここで便利です。

Numberの部分型でなくとも、加法っぽい可換な二項演算に対してはBase.:+にメソッドを追加して良い。^[3]

型の昇格 (type promotion)

型の異なる足し算については、type promotionの仕組みによって同じ型に揃えられてから足し算が実行されます。
上記の例では、k1がKnotVector{Float64}に変換されてから足し算が実行されます。

BasicBSplineでは以下のようにpromote_ruleにメソッドを追加することで明示的に型変換を定義しています。^[4]

function Base.promote_rule(::Type{KnotVector{T}}, ::Type{KnotVector{S}}) where {T,S}
    KnotVector{promote_type(T,S)}
end

Juliaの型の昇格は暗黙の型変換に近いですが、暗黙さが低いです。型変換が言語にハードコーディングされている訳ではなく、ユーザーが自由に設定することができます。

B-spline空間 (区分多項式空間)

ノット列を区分点に持つ$p$次区分多項式全体の空間を$p$次B-spline空間と呼びます。^[5]
この空間は線形空間で、空間の次元や空間どうしの包含関係などが定義できます。

多項式次数と型パラメータ

BasicBSplineでは、B-spline空間を表す抽象型としてAbstractBSplineSpace{p}が用意されています。
ここでpは多項式次数を表す型パラメータで、構造体のフィールドに含まれる実装にはなっていません。
このような型パラメータを使った理由は以下になります。

• ほとんどの場合で多項式次数は実行時に決定できる。よって型パラメータを使って問題ない。
• 後述の生成関数(@generated)を使うことで高速な実装が実現できる。

structのフィールドにするか型パラメータにするかは

• 型パラメータにした際に速度向上が望めるか
• 実行時に型安定性を望めるか

を基準に考えれば良い。

プロット

BasicSplineでは、B-spline空間をプロットすることができます。

using BasicBSpline
using Plots

k = KnotVector([0.0, 0.5, 1.5, 2.5, 3.5, 5.5, 8.0, 9.0, 9.5, 10.0])
P = BSplineSpace{3}(k)
plot(P, ylims=(0,1))

空間をプロットすると言うと分かりにくいですが、有限次元の関数空間なのでその基底関数すべてをプロットすれば空間の可視化になっています。

このようなプロットの実装は、RecipesBaseの@recipeマクロを使えば実装が可能です。
BasicBSpline.jlでの実装例はこちらを参照してください。

独自に定義した型の標準的なプロットを追加するには、RecipesBase.@recipeでプロットを定義できる。
RecipesBase.jlパッケージは軽量なので、依存関係の追加を躊躇わなくて良い。

生成関数を使った基底関数の実装

BasicBSplineには、基底関数の値を計算するための関数としてbsplinebasisとbsplinebasisallの2種類が用意されています。

• bsplinebasis(P, i, t): i番目の基底関数のtでの値。
  • Cox-de Boorの漸化式に基づく素朴な実装
  • 実験に使う際に便利
• bsplinebasisall(P::BSplineSpace{p}, j, t): j番目の区間で非ゼロになっている基底関数のtでの値たち。
  • 戻り値の型はSVector{p+1}
  • 効率的な実装のためにはbsplinebasisよりもbsplinebasisallの方が良い

これらの関数を実装するために、生成関数(@generated)が使われています。
たとえば、bsplinebasisall関数は以下のように実装されています。

@inline function bsplinebasisall(P::AbstractBSplineSpace{0,T},i::Integer,t::S) where {T, S<:Real}
    U = StaticArrays.arithmetic_closure(promote_type(T,S))
    SVector(one(U),)
end

@inline function bsplinebasisall(P::AbstractBSplineSpace{1}, i::Integer, t::Real)
    k = knotvector(P)
    B1 = (k[i+2]-t)/(k[i+2]-k[i+1])
    B2 = (t-k[i+1])/(k[i+2]-k[i+1])
    return SVector(B1, B2)
end

@generated function bsplinebasisall(P::AbstractBSplineSpace{p}, i::Integer, t::Real) where p
    bs = [Symbol(:b,i) for i in 1:p]
    Bs = [Symbol(:B,i) for i in 1:p+1]
    K1s = [:((k[i+$(p+j)]-t)/(k[i+$(p+j)]-k[i+$(j)])) for j in 1:p]
    K2s = [:((t-k[i+$(j)])/(k[i+$(p+j)]-k[i+$(j)])) for j in 1:p]
    b = Expr(:tuple, bs...)
    B = Expr(:tuple, Bs...)
    exs = [:($(Bs[j+1]) = ($(K1s[j+1])*$(bs[j+1]) + $(K2s[j])*$(bs[j]))) for j in 1:p-1]
    Expr(:block,
        :($(Expr(:meta, :inline))),
        :(k = knotvector(P)),
        :($b = bsplinebasisall(_lower(P),i+1,t)),
        :($(Bs[1]) = $(K1s[1])*$(bs[1])),
        exs...,
        :($(Bs[p+1]) = $(K2s[p])*$(bs[p])),
        :(return SVector($(B)))
    )
end

@generatedマクロによる関数の定義は読みにくいですが、この生成関数によるbsplinebasisallの定義は以下のような実装とほぼ等価です。

@inline function bsplinebasisall(P::AbstractBSplineSpace{2}, i::Integer, t::Real)
    k = knotvector(P)
    (b1, b2) = bsplinebasisall(BSplineSpace{1}(k), i + 1, t)
    B1 = ((k[i + 3] - t) / (k[i + 3] - k[i + 1])) * b1
    B2 = ((k[i + 4] - t) / (k[i + 4] - k[i + 2])) * b2 + ((t - k[i + 1]) / (k[i + 3] - k[i + 1])) * b1
    B3 = ((t - k[i + 2]) / (k[i + 4] - k[i + 2])) * b2
    return SVector((B1, B2, B3))
end

@inline function bsplinebasisall(P::AbstractBSplineSpace{3}, i::Integer, t::Real)
    k = knotvector(P)
    (b1, b2, b3) = bsplinebasisall(BSplineSpace{2}(k), i + 1, t)
    B1 = ((k[i + 4] - t) / (k[i + 4] - k[i + 1])) * b1
    B2 = ((k[i + 5] - t) / (k[i + 5] - k[i + 2])) * b2 + ((t - k[i + 1]) / (k[i + 4] - k[i + 1])) * b1
    B3 = ((k[i + 6] - t) / (k[i + 6] - k[i + 3])) * b3 + ((t - k[i + 2]) / (k[i + 5] - k[i + 2])) * b2
    B4 = ((t - k[i + 3]) / (k[i + 6] - k[i + 3])) * b3
    return SVector((B1, B2, B3, B4))
end

@inline function bsplinebasisall(P::AbstractBSplineSpace{4}, i::Integer, t::Real)
    k = knotvector(P)
    (b1, b2, b3, b4) = bsplinebasisall(BSplineSpace{3}(k), i + 1, t)
    B1 = ((k[i + 5] - t) / (k[i + 5] - k[i + 1])) * b1
    B2 = ((k[i + 6] - t) / (k[i + 6] - k[i + 2])) * b2 + ((t - k[i + 1]) / (k[i + 5] - k[i + 1])) * b1
    B3 = ((k[i + 7] - t) / (k[i + 7] - k[i + 3])) * b3 + ((t - k[i + 2]) / (k[i + 6] - k[i + 2])) * b2
    B4 = ((k[i + 8] - t) / (k[i + 8] - k[i + 4])) * b4 + ((t - k[i + 3]) / (k[i + 7] - k[i + 3])) * b3
    B5 = ((t - k[i + 4]) / (k[i + 8] - k[i + 4])) * b4
    return SVector((B1, B2, B3, B4, B5))
end

# AbstractBSplineSpace{p}のpに関して無限につづく…

このような加算個の定義をコードに書くのは大変なので、@generatedで生成関数を定義するのが便利です。

Juliaでは生成関数によって型パラメータごとに最適化されたコードを生成することができる。

微分の実装

BasicBSplineでは微分のための型BSplineDerivativeSpace{r}が用意されています。

using BasicBSpline
using Plots

k = KnotVector([0.0, 0.5, 1.5, 2.5, 3.5, 5.5, 8.0, 9.0, 9.5, 10.0])
P = BSplineSpace{3}(k)
dP = BSplineDerivativeSpace{1}(P)
ddP = BSplineDerivativeSpace{2}(P)
plot(P, label="P")
plot!(dP, label="dP")
plot!(ddP, label="ddP")

このdPのような空間は、「もとの関数空間のそれぞれの元を$r$階導関数全体の空間」として定義されています。
微分によって定義された空間dPなどに対してもbsplinebasisやbsplinebasisall関数が定義されています。^[6]

自動微分に関わらず、微分のためのインターフェースを用意するのは良いこと。
単純に微分計算したいユースケースに対応できるし、自動微分のルールの実装にも役立つ。

B-spline多様体

「B-spline基底関数の線型結合で表される写像の像」をB-spline多様体と言います。^[7]
BasicBSplineにおけるB-spline多様体を表す型がBSplineManifoldです。

型の実装

BSplineManifoldは以下のように実装されています。

struct BSplineManifold{Dim,Deg,C,S<:NTuple{Dim, AbstractBSplineSpace}} <: AbstractBSplineManifold{Dim,Deg}
    bsplinespaces::S
    controlpoints::Array{C,Dim}
    function BSplineManifold(a::Array{C,Dim},Ps::S) where {S<:NTuple{Dim, AbstractBSplineSpace},C} where Dim
        if size(a) != dim.(Ps)
            msg = "The size of control points array $(size(a)) and dimensions of B-spline spaces $(dim.(Ps)) must be equal."
            throw(DimensionMismatch(msg))
        end
        Deg = degree.(Ps)
        new{Dim,Deg,C,S}(Ps,a)
    end
end

型パラメータが少し複雑ですが、それぞれ以下の意味を持っています。

• Dim
  • 多様体の次元を表す。
  • B-spline空間の数と制御点の配列の次元がこの次元に相当する
• Deg
  • 各次元に対応する多項式次数をタプルとして保持している
  • 例えば、第1方向に$3$、第2方向に$2$の多項式次数を持っていればDeg = (1,2)。
  • StaticArraysではタプル(1,2)の代わりにTuple{1,2}が使われているが、Tuple型の使い方から離れすぎているので、タプル(1,2)の方が"正しい"はず… (StaticArrays.jl#807)
  • 型パラメータの数は固定なので、可変長のパラメータを入れるためにはタプルにするのが良い。
• C
  • 制御点の型を表す。
  • ベクトル空間の元のように扱える型(Number, AbstractVectorなど)であればOK。
• S
  • 各次元のBSplineSpaceを保持するための型。
  • Degと情報が被っているので少し冗長。
  • しかし、AbstractBSplineManifoldにDegが含まれているのでどちらも削ることはできない。

型パラメータの数は型によって固定されているので、可変長のパラメータを入れるためにはタプルにするのが良い。

型パラメータは省略可能なので、重要度の高いものから並べるのが良い

プロット

BSplineSpaceと同様に、BSplineManifoldにもプロットが定義されています。

using BasicBSpline
using StaticArrays
using Plots
plotly()

k1 = KnotVector([0.0, 0.5, 1.5, 2.5, 3.5, 5.5, 8.0, 9.0, 9.5, 10.0])
k2 = KnotVector([0, 0, 1, 1])
P1 = BSplineSpace{3}(k1)
P2 = BSplineSpace{1}(k2)
a = [-j*SVector(cos(i-1), sin(i-1), rand()) for i in 1:dim(P1), j in 1:dim(P2)]
M = BSplineManifold(a, P1, P2)
plot(M)

を実行すると以下のようなプロットが得られれます。

BSplineManifoldのために用意されているプロットは

• 2次元平面上の曲線(BSplineManifold{1,Deg,<:StaticVector{2,<:Real}})
• 3次元空間中の曲線(BSplineManifold{1,Deg,<:StaticVector{3,<:Real}})
• 3次元空間中の曲面(BSplineManifold{2,Deg,<:StaticVector{3,<:Real}})

の3つです。2次元平面上の曲面もプロットできると便利ですが、現在のPlots.jlにはそのうような曲面をプロットする仕組みが用意されていないので実装できていません。

制御点を含んだ曲面などの、複雑なプロットもRecipesBase.@recipeで定義できる。
ただし、プロットの種類はPlots.jlでサポートしているものに限られてしまう。

写像の計算

B-spline多様体には標準的なパラメータが備わっており、BSplineManifoldのインスタンスがこの写像になっています。

julia> using StaticArrays

julia> P = BSplineSpace{2}(KnotVector([0,0,0,1,1,1]))
BSplineSpace{2, Int64}(KnotVector([0, 0, 0, 1, 1, 1]))

julia> a = [SVector(1,0), SVector(1,1), SVector(0,1)]
3-element Vector{SVector{2, Int64}}:
 [1, 0]
 [1, 1]
 [0, 1]

julia> M = BSplineManifold(a, P);

julia> M(0.4)  # Mが関数のように振る舞う
2-element SVector{2, Float64} with indices SOneTo(2):
 0.84
 0.64

julia> M(1.2)  # 定義域の外に出るとエラー
ERROR: DomainError with 1.2:
The input 1.2 is out of range.

これはJuliaのfunction-like objectの仕組みで実装できます。^[8]

function-like objectを使えば、インスタンスを関数のように扱えて便利

ChainRulesへの対応

ChainRulesとは、Juliaで自動微分を計算させる際のルールを定義するための仕組みのことです。
このルールとは、「$\sin(x)$を微分して$\cos(x)$になる」といった微分規則のことで、これをパッケージごとに定義することで自動微分の計算精度や計算速度を向上させることができます。

実際、BasicBSpline.jlでChainRules.jlを定義することで60倍以上高速化できました。

ChainRulesCore.fruleとChainRulesCore.rruleへメソッドを追加するだけで、微分の規則を定義できる。

ChainRules.jlでは複素共役がrruleに現れることがよくある。(たとえばsinの微分の定義)
この複素共役は以下が理由だが、B-splineのような実関数のみを扱っている限りは気にしなくて良い。

• リバースモードの自動微分においては転置(あるいは随伴)で書く方が都合が良い
• 複素数体で線形代数をやるときは転置より随伴の方が性質が良いことが多い(要出典)
• Juliaの記法ではtransposeよりもadjoint(')の方が書きやすい
• 以上の理由でChainRulesのrruleでは転置ではなく随伴を使う
• sin(x::Number)の場合でも、rruleでは随伴に相当する複素共役を取る

Q&A

いつv1.0を出すんですか？

2023年には出せそうな雰囲気です。
予定としては以下の順で取り組む予定です。

• 現在検討している破壊的変更(後述)に関するwarningの追加
• 破壊的変更を加えてv0.9.0をリリース
• Julia Discourseにパッケージアナウンス
• 1ヶ月くらいフィードバックを待ってからv1.0.0をリリース

破壊的変更の予定① (抽象型の再設計)

現在はB-spline空間を表す抽象型としてAbstractBSplineSpaceが用意されていて、BSplineSpaceとUniformBSplineSpaceを部分型にもっています。
これらの具象型はそれぞれノット列のKnotVectorとUniformKnotVectorに対応していますが、この設計が間違っていました。

正しくは、AbstractBSplineSpaceの抽象型は不要で、BSplineSpace型のパラメータとしてT <: AbstractKnotVectorを持つべきでした。
この理由は以下です：

• ノット列の具象型とB-spline空間の具象型が1対1対応する
• ノット列の具象型を増やすたびにB-spline空間の具象型を新しく作るのは面倒

破壊的変更の予定② (fittingcontrolpointsの移行)

BasicBSplineの野望として、「JuliaにおけるB-splineのスタンダードになる」があります！

これを実現するためにはパッケージは軽量で、不要な依存先を減らす必要があります。
コード例で紹介した関数のフィッティングでは内部的にGauss求積を使っており、このためにFastGaussQuadrature.jlに依存してしまっています。
関数フィッティングは便利ですがB-splineに必須ではないので、fittingcontolpoints関数を別のBasicBSplineFitting.jlパッケージに移行することを検討しています。

他のパッケージよりも優れてるんですか？

B-splineを扱うJuliaパッケージには他にもあって、例えば以下です。

• BSplineKit.jl
• BSplines.jl

BasicBSplineパッケージが優れている所は以下にあると思っています。

• 数学的な合理性を重視している
• リファインメントなどのアルゴリズムの実装
• 軽量な依存パッケージ

しかし、他のパッケージを詳細に調査した訳ではないので少し怪しい部分もあります。
詳しい方がいればissue#161で助けてください！

脚注

1. 子を持つ型が抽象型(e.g. Real)で、子を持てない型が具象型(e.g. Float64)です。

2. この「ノット列を足す操作」は可換なので、+の記号を採用しました。文字列の結合は非可換なのでJuliaでは*が使われます。∪を避けたのは重複した要素を保持することを強調したかったためです。

3. この辺りの基準は少し曖昧です。元々の関数の意図を大きく超えたメソッドを追加することはType III piracyと呼ばれており、避けた方が良いとされます。

4. このpromotionの仕組みの便利なところは、型Aと型Bの昇格のルールを決めるだけで、AとBを使ったメソッドすべてに型変換が利用できることにあります。
   例えば、AとBが両方ともRealの部分型であれば、+(::A, ::B), -(::A, ::B), *(::A, ::B), /(::A, ::B)などのすべてに定義が使えます。

5. 独自用語です。広く普及している名称では無いことに注意してください。B-spline基底関数が基底になっています。

6. $\frac{dB_{(i,p,k)}(t)}{dt}$は微分した関数空間の基底になっていないですが、インターフェースの共通化のためにこのような実装になっています。

7. 厳密には多様体ではありません。B-spline基底関数の線型結合で表される写像が単射でないこともありますし、Jacobi行列のランクが落ちていることもあります。また、ここでのB-spline基底関数はテンソル積によって多次元に拡張されていることに注意してください。

8. 現在は3次元のB-spline多様体までしか対応できていません。生成関数で生成するコードが「次元と各次元の多項式の次数」に対応する必要があり、面倒なのが理由です。技術的には任意の次元まで対応できるはずです。