Line data Source code
1 : """
2 : Defines interface for cubic Hermite Splines.
3 : """
4 : abstract type AbstractHermiteSpline{S<:Real, V} end
5 :
6 : """
7 : h::AbstractHermiteSpline
8 : t, x, dx = h[i]
9 :
10 : Get knotr value, value and derivative vector at knot index `i`.
11 :
12 : See also [`knot`](@ref), [`value`](@ref), [`derivative`](@ref)
13 : """
14 1 : Base.getindex(h::AbstractHermiteSpline, i) =
15 : error("not implemented") # keep only for documentation
16 :
17 2 : Base.lastindex(h::AbstractHermiteSpline) = length(h)
18 0 : Base.firstindex(h::AbstractHermiteSpline) = 1
19 :
20 : """
21 : eltype(h::AbstractHermiteSpline)
22 :
23 : Get type `Tuple` type of [`knot`](@ref), [`value`](@ref),
24 : [`derivative`](@ref). This is the type of `h[idx]`.
25 :
26 : See also [`knottype`](@ref), [`valuetype`](@ref)
27 : """
28 4 : Base.eltype(::Type{H}) where {S, V, H<:AbstractHermiteSpline{S, V}} =
29 : Tuple{S, V, V}
30 :
31 : """
32 : n = length(h::AbstractHermiteSpline)
33 :
34 : Get number of [`knots`](@ref)
35 : """
36 11819 : Base.length(h::AbstractHermiteSpline) = length(h.t)
37 :
38 : """
39 : for i in eachindex(h::AbstractHermiteSpline)
40 :
41 : Iterate over [`knots`](@ref) and control points.
42 : """
43 0 : Base.eachindex(h::AbstractHermiteSpline) = eachindex(h.t)
44 :
45 : """
46 : isempty(h::AbstractHermiteSpline)
47 :
48 : [`knot`](@ref) vector is empty or contains only one single knot.
49 :
50 : See also [`length`](@ref)
51 : """
52 34 : Base.isempty(h::AbstractHermiteSpline) = length(h) < 2
53 :
54 7 : Base.iterate(h::AbstractHermiteSpline, idx = 1) =
55 : idx > length(h) ? nothing : (h[idx], idx + 1)
56 :
57 : """
58 : knottype(h::AbstractHermiteSpline)
59 :
60 : Get type of [`knot`](@ref).
61 :
62 : See also [`eltype`](@ref), [`valuetype`](@ref)
63 : """
64 3 : knottype(h::AbstractHermiteSpline{S, V}) where {S, V} = S
65 :
66 : """
67 : valuetype(h::AbstractHermiteSpline)
68 :
69 : Get type of [`value`](@ref).
70 :
71 : See also [`eltype`](@ref), [`knottype`](@ref)
72 : """
73 4 : valuetype(h::AbstractHermiteSpline{S, V}) where {S, V} = V
74 :
75 :
76 : """
77 : t = knot(h::AbstractHermiteSpline, i)
78 :
79 : Get `i`-th knot.
80 :
81 : See also [`knots`](@ref)
82 : """
83 45078 : knot(h::AbstractHermiteSpline, i::Int) = first(h.t[i])
84 :
85 : """
86 : t0 = first(h::AbstractHermiteSpline)
87 :
88 : Get first [`knot`](@ref)
89 :
90 : See also [`last`](@ref), [`span`](@ref)
91 : """
92 22 : Base.first(h::AbstractHermiteSpline) = knot(h, 1)
93 :
94 : """
95 : t1 = last(h::AbstractHermiteSpline)
96 :
97 : Get last [`knot`](@ref)
98 :
99 : See also [`first`](@ref), [`span`](@ref)
100 : """
101 351 : Base.last(h::AbstractHermiteSpline) = knot(h, length(h))
102 :
103 : """
104 : t0, t1 = span(h::AbstractHermiteSpline)
105 :
106 : Get [`first`](@ref) and [`last`](@ref) [`knot`](@ref). They
107 : determine the parameter interval of `h`.
108 :
109 : !!! note
110 : The spline can certainly *extrapolated* beyond `t0` or `t1`.
111 : Note, however, that extrapolaiton is numerically unstable if the
112 : parameter is too far w.r.t. the length of the knot interval
113 : (evaluation of monomials away from ``[0, 1]``).
114 : """
115 8 : span(h::AbstractHermiteSpline) = (first(h), last(h))
116 :
117 : """
118 : isascending(h::AbstractHermiteSpline)
119 :
120 : Does `h` have [`knots`](@ref) in ascending order? This is
121 : determined by testing the first two knots.
122 :
123 : !!! note
124 : Requires *not* [`isempty`](@ref).
125 :
126 : See also [`ismonotone`](@ref)
127 : """
128 16054 : isascending(h::AbstractHermiteSpline) = (knot(h, 1) <= knot(h, 2))
129 :
130 9 : function ismonotone(h::AbstractHermiteSpline)
131 9 : isempty(h) && return true
132 :
133 7 : if isascending(h)
134 4 : all(knot(h, i-1) < knot(h, i) for i in 2:length(h))
135 : else
136 3 : all(knot(h, i-1) > knot(h, i) for i in 2:length(h))
137 : end
138 : end
139 :
140 343 : function ismonotone(h::AbstractHermiteSpline{S, V}, t::S) where {S, V}
141 343 : n = length(h)
142 343 : (n == 0) && return true
143 :
144 330 : s = last(h)
145 330 : (n == 1) && return (s != t) # isempty(h) == true
146 :
147 : # NOTE: isempty fails, because may insert first knot twice
148 : # NOTE: isascending cannot be evaluated for empty spline
149 :
150 317 : isascending(h) ? (s < t) : (s > t)
151 : end
152 :
153 12517 : _leq(h::AbstractHermiteSpline, s, t) = isascending(h) ? (s < t) : ( s > t)
154 :
155 : """
156 : ismonotone(h::AbstractHermiteSpline)
157 :
158 : Test if [`knot`](@ref) sequence is *strictly* monotone. This
159 : method iterates over all [`knots`](@ref) and is intended for
160 : (internal) checks.
161 :
162 : !!! note
163 : Returns `true` if [`isempty`](@ref)
164 :
165 : ismonotone(h::AbstractHermiteSpline, t)
166 :
167 : Test if `t` could be added to the knot sequence such that
168 : the sequence is still (strictly) monotone. This method is
169 : is called by [`push!`](@ref)
170 :
171 : !!! note
172 : Returns `true` if [`isempty`](@ref), returns `false` if
173 : `t == last(h)` (zero length knot interval).
174 :
175 : See also [`push!`](@ref)
176 : """
177 : ismonotone
178 :
179 : """
180 : push!(h::AbstractHermiteSpline, t, x, dx)
181 :
182 : Add [`knot`](@ref) `t` with [`value`](@ref) `x` and [`derivative`](@ref)
183 : `dx`.
184 :
185 : !!! note
186 : Requires (and checks) `ismonotone(h, t)`.
187 :
188 : See also [`ismonotone`](@ref).
189 : """
190 330 : function Base.push!(h::AbstractHermiteSpline{S, V},
191 : t::S, x::V, dx::V) where {S<:Real, V}
192 330 : @assert ismonotone(h, t)
193 330 : _push!(h, t, x, dx)
194 330 : h
195 : end
196 :
197 : """
198 : i = index(h::AbstractHermiteSpline, t[, lo=1])
199 :
200 : Find index of [`knot`](@ref) interval for parameter `t`.
201 :
202 : The optional argument `lo` defines the initial lower (left) bound for
203 : the binary serach, which must be a valid index!
204 :
205 : !!! note
206 : The method returns an invalid index (`0` or `length(h)`) if
207 : `t ∉ span(h)`. This method is intended for *internal use*.
208 :
209 : See also [`piece`](@ref), [`boundedindex`](@ref)
210 : """
211 4490 : index(h::AbstractHermiteSpline{S, V}, t::S, lo::Int = 1) where {S, V} =
212 : _searchsortedlast(h.t, t, lo, isascending(h))
213 :
214 : """
215 : i = boundedindex(h::AbstractHermiteSpline, t, [lo = 1])
216 :
217 : Find index of [`knot`](@ref) interval for parameter `t`.
218 : Returns the index of the first or last [`knot`](@ref) interval if
219 : `t < first(h)` or `t > last(h)`, respectvely, for *extrapolation*.
220 :
221 : !!! note
222 : The result is *undefined* if [`isempty`](@ref)
223 :
224 : See also [`index`](@ref)
225 : """
226 5325 : function boundedindex(h::AbstractHermiteSpline{S, V},
227 : t::S, idx::Int = 1) where {S, V}
228 5325 : lo = 1
229 :
230 5320 : if isindex(h, idx)
231 6015 : if _leq(h, knot(h, idx), t)
232 5380 : _leq(h, t, knot(h, idx + 1)) && return idx
233 3920 : lo = idx
234 : end
235 : end
236 :
237 4490 : i = index(h, t, lo)
238 4490 : i = max(1, i)
239 4490 : i = min(length(h) - 1, i)
240 4490 : i
241 : end
242 :
243 :
244 : """
245 : isindex(h::AbstractHermiteSpline, i, t)
246 :
247 : Is `i` a valid index for `h` such that `t` is within the `i`-th
248 : [`knot`](@ref) interval?
249 :
250 : !!! note
251 : Returns `false` if `t < first(s)` or `last(s) < t` (extrapolation).
252 :
253 : See also [`index`](@ref)
254 : """
255 561 : isindex(h::AbstractHermiteSpline{S, V}, i::Int, t::S) where {S, V} =
256 : isindex(h, i) && _leq(h, knot(h, i), t) && _leq(h, t, knot(h, i + 1))
257 :
258 : """
259 : isindex(h::AbstractHermiteSpline, i)
260 :
261 : Is `i` a valid index for a [`knot`](@ref) interval in `h`?
262 :
263 : See also [`index`](@ref)
264 : """
265 5881 : isindex(h::AbstractHermiteSpline, i::Int) = (1 <= i <= length(h) - 1)
266 :
267 11 : function piece(h::AbstractHermiteSpline, i::Int)
268 11 : t0, x0, dx0 = h[i]
269 11 : t1, x1, dx1 = h[i+1]
270 11 : HermitePiece(x0, dx0, x1, dx1, t0, t1)
271 : end
272 :
273 14 : piece(h::AbstractHermiteSpline{S, V}, t::S, i::Int = 1) where {S<:Real, V} =
274 : piece(h, boundedindex(h, t, i))
275 :
276 : """
277 : p = piece(h::AbstractHermiteSpline, i::Int)
278 : p = piece(h::AbstractHermiteSpline, t[, i::Int])
279 :
280 : Get [`HermitePiece`](@ref) for evaluation from index `i` or parameter `t`.
281 : No search if `t` is within `i`-th interval.
282 :
283 : !!! note
284 : - Requires *not* `isempty(h)`.
285 : - For *indices*, requires `1 <= i < length(h)`.
286 : - For *parameters*, computes [`index`](@ref) but returns the first or
287 : last piece if `t ∉ span(h)`.
288 : """
289 : piece
290 :
291 : """
292 : p = somepiece(h::AbstractHermiteSpline)
293 :
294 : Get some arbitrary, possibly invalid [`HermitePiece`](@ref) for `h`. The intended use
295 : is initialization of an undefine piece.
296 : """
297 6 : function somepiece(h::AbstractHermiteSpline{S, V}) where {S, V}
298 9 : isempty(h) || return piece(h, 1)
299 :
300 3 : v = zeros(V)
301 3 : HermitePiece(v, v, v, v, S(Inf), S(Inf))
302 : end
303 :
304 17 : function hermitespline(dim::Val{N},
305 : ::Type{S}, ::Type{H}) where {H<:AbstractHermiteSpline, N, S<:Real}
306 17 : H{S, _hs_vector(dim, S)}()
307 : end
308 :
309 1 : _hs_vector(::Val{1}, ::Type{S}) where {S<:Real} = S
310 16 : _hs_vector(::Val{N}, ::Type{S}) where {S<:Real, N} = SVector{N, S}
311 :
312 17 : function hermitespline(dim::Val{N},
313 : scalar::Type{S}, impl::Val{H}) where {H, N, S<:Real}
314 17 : hermitespline(dim, scalar, _hs_impl(impl))
315 : end
316 :
317 5 : _hs_impl(::Val{:t_x_dx}) = HermiteSpline_t_x_dx
318 21 : _hs_impl(::Val{:t_xdx}) = HermiteSpline_t_xdx
319 4 : _hs_impl(::Val{:txdx}) = HermiteSpline_txdx
320 :
321 20 : hermitespline(dim::Int, scalar::Type = Float64, impl::Symbol = :t_xdx) =
322 : hermitespline(Val(dim), scalar, Val(impl))
323 :
324 7 : hermitespline(::Type{SVector{N, S}},
325 : impl::Symbol = :t_xdx) where {N, S<:Real} =
326 : hermitespline(Val(N), S, Val(impl))
327 :
328 2 : hermitespline(::SVector{N, S},
329 : impl::Symbol = :t_xdx) where {N, S<:Real} =
330 : hermitespline(SVector{N, S}, impl)
331 :
332 3 : hermitespline(t::Vector{S}, x::Vector{V}, dx::Vector{V},
333 : impl::Symbol = :t_x_dx) where {S, V} =
334 : _hermitespline(Val(impl), t, x, dx)
335 :
336 11 : hermitespline(t::Vector{S}, xdx::Vector{Tuple{V, V}},
337 : impl::Symbol = :t_xdx) where {S, V} =
338 : _hermitespline(Val(impl), t, xdx)
339 :
340 3 : hermitespline(txdx::Vector{Tuple{S, V, V}},
341 : impl::Symbol = :txdx) where {S, V} =
342 : _hermitespline(Val(impl), txdx)
343 :
344 13 : _hermitespline(impl::Val{H}, args...) where {H} =
345 : _hermitespline(args..., _hs_impl(impl))
346 :
347 3 : _hermitespline(t::Vector{S}, x::Vector{V}, dx::Vector{V},
348 : ::Type{H}) where {H<:AbstractHermiteSpline, S, V} = H{S, V}(t, x, dx)
349 :
350 7 : _hermitespline(t::Vector{S}, xdx::Vector{Tuple{V, V}},
351 : ::Type{H}) where {H<:AbstractHermiteSpline, S, V} = H{S, V}(t, xdx)
352 :
353 3 : _hermitespline(txdx::Vector{Tuple{S, V, V}},
354 : ::Type{H}) where {H<:AbstractHermiteSpline, S, V} = H{S, V}(txdx)
355 :
356 : # TODO: still a dynamic dispatch on Val for hermitespline(t, x, dx, :t_xdx) ?
357 :
358 : """
359 : ### Empty spline
360 :
361 : s = hermitespline(n, type[, impl = :t_xdx])
362 : s = hermitespline(Val(n), type, Val(impl))
363 : s = hermitespline(SVector{n, type}[, impl = :t_xdx])
364 :
365 : Construct *empty* Hermite spline with scalar knot `type` and
366 : `n`-vectors (`SVector{n, type}`) for [`values`](@ref) and
367 : [`derivatives`](@ref). To "fill" data use, e.g., [`push!`](@ref) or
368 : [`append!`](@ref)
369 :
370 : The implementation defines the storage layout. Possible values are
371 : `:t_x_dx`, `:t_xdx` and `:txdx`, where speration by `"_"` indicates
372 : use of a `Vector` (e.g., for `:t_x_dx` one for [`knots`](@ref),
373 : [`values`](@ref) and [`derivatives`](@ref)).
374 :
375 : !!! note
376 : For `n==1`, i.e., scalar values and derivatives, values and derivatives
377 : are stored as scalars of `type` (not `SVector{1, type}`).
378 :
379 : ### Spline from data:
380 :
381 : s = hermitespline(t, x, dx[, impl = :t_x_dx])
382 : s = hermitespline(t, xdx[, impl = :t_xdx])
383 : s = hermitespline(txdx[, impl = :txdx])
384 :
385 : Construct Hermite spline from knots `t`, values `x` and derivatives `dx` or
386 : vectors of tuples.
387 :
388 : !!! warning "precondition"
389 : 1. The passed vectors must have *equal length*.
390 : 2. The passed vectors are for exclusive use by `hermitespline`,
391 : the constructor *doesn't copy*!
392 :
393 : See also [`AbstractHermiteSpline`](@ref)
394 : """
395 : hermitespline
396 :
397 3 : function Base.show(io::IO, h::AbstractHermiteSpline)
398 3 : if isempty(h)
399 2 : (length(h) == 0) && return print(io, "HermiteSpline(empty)")
400 1 : return print(io, "HermiteSpline(empty, single knot $(first(h)))")
401 : end
402 :
403 1 : print(io, "HermiteSpline($(length(h)) knots in $(span(h)))")
404 : end
405 :
406 :
407 :
408 : """
409 : t = knots(h::AbstractHermiteSpline)
410 :
411 : Get (possibly a generator of) knot sequence. (Use `collect` to create
412 : an `Array`.)
413 :
414 : See also [`knot`](@ref), [`values`](@ref), [`derivatives`](@ref)
415 : """
416 1 : knots(h::AbstractHermiteSpline) = (knot(h, i) for i in 1:length(h))
417 : # Note: May be specialized, generic implementation for documentation only!
418 :
419 : """
420 : value(h::AbstractHermiteSpline, i)
421 :
422 : Get value at [`knot`](@ref) index `i`.
423 :
424 : See also [`values`](@ref)
425 : """
426 210 : value(h::AbstractHermiteSpline, i::Int) = h[i][2]
427 :
428 : """
429 : t = values(h::AbstractHermiteSpline)
430 :
431 : Get (possibly a generator of) values per knot. (Use `collect` to
432 : create an `Array`.)
433 :
434 : See also [`value`](@ref), [`knots`](@ref), [`derivatives`](@ref)
435 : """
436 1 : Base.values(h::AbstractHermiteSpline) = (value(h, i) for i in 1:length(h))
437 :
438 : """
439 : derivative(h::AbstractHermiteSpline, i)
440 :
441 : Get derivative vector at [`knot`](@ref) index `i`.
442 :
443 : See also [`derivativess`](@ref)
444 : """
445 210 : derivative(h::AbstractHermiteSpline, i::Int) = h[i][3]
446 :
447 : """
448 : t = derivatives(h::AbstractHermiteSpline)
449 :
450 : Get (possibly a generator of) derivative vectors per knot. (Use `collect`
451 : to create an `Array`.)
452 :
453 : See also [`derivative`](@ref), [`knots`](@ref), [`valuess`](@ref)
454 : """
455 1 : derivatives(h::AbstractHermiteSpline) = (derivative(h, i) for i in 1:length(h))
456 :
457 : """
458 : similar(h::AbstractHermiteSpline)
459 :
460 : Get empty spline of `tyepof(h)`.
461 : """
462 3 : Base.similar(h::H) where {H<:AbstractHermiteSpline} = H()
463 :
464 : """
465 : hc = copy(h::AbstractHermiteSpline)
466 :
467 : Get (deep) copy of `h`.
468 : """
469 10 : Base.copy(h::H) where {H<:AbstractHermiteSpline} =
470 : H((copy(getproperty(h, name)) for name in propertynames(h))...)
471 :
472 : """
473 : copy(hdst::H, h::H) # where {H<:AbstractHermiteSpline}
474 :
475 : Deep copy of `h` into `hdst`
476 : """
477 3 : function Base.copy!(hdst::H, h::H) where {H<:AbstractHermiteSpline}
478 3 : for name in propertynames(h)
479 9 : copy!(getproperty(hdst, name), getproperty(h, name))
480 8 : end
481 3 : hdst
482 : end
483 :
484 : """
485 : reverse!(h::AbstractHermiteSpline)
486 :
487 : Reverse knot vector and coefficient vector(s).
488 : """
489 3 : function Base.reverse!(h::AbstractHermiteSpline)
490 3 : for name in propertynames(h)
491 9 : reverse!(getproperty(h, name))
492 8 : end
493 3 : h
494 : end
495 :
496 : """
497 : hr = reverse!(h::AbstractHermiteSpline)
498 :
499 : Copy `h` with reversed knot vector and coefficient vector(s).
500 : """
501 8 : Base.reverse(h::H) where {H<:AbstractHermiteSpline} =
502 : H((reverse(getproperty(h, name)) for name in propertynames(h))...)
503 :
504 :
505 : """
506 : append!(r::H, s::H) # H<:AbstractHermiteSpline
507 :
508 : Concatenate Hermite splines `r` and `s` by appending all [`knots`](@ref),
509 : [`values`](@ref) and [`derivative`](@ref) of `s` to `r`.
510 :
511 : !!! warning "precondition"
512 : The knot sequences of `s` and `t` must be compatible ([`ismonotone`](@ref)).
513 : If knots [`last`](@ref)`(s)` and [`first`](@ref)`(r)` are identical *and*
514 : their values and derivatives coincide, the first knot of `s` is ignored
515 : and will *not* be append.
516 :
517 : See also [`vcat`](@ref)
518 : """
519 12 : function Base.append!(r::H, s::H) where {S, V, H<:AbstractHermiteSpline{S, V}}
520 12 : (length(r) == 0) && return copy!(r, s)
521 12 : (length(s) == 0) && return r
522 :
523 12 : ifirst = 1
524 12 : n = length(r)
525 12 : if last(r) == first(s)
526 10 : (r[n] == s[1]) || error("double knot $(r[n]) (discontinuous)")
527 8 : ifirst = 2
528 : end
529 11 : (length(s) < ifirst) && return r
530 :
531 12 : (isascending(r) == isascending(s)) || error("mismatch in order of knots")
532 11 : ismonotone(r, knot(s, ifirst)) || error("overlapping knot sequence")
533 :
534 9 : for i in ifirst:length(s)
535 10 : push!(r, s[i]...)
536 11 : end
537 :
538 9 : r
539 : end
540 :
541 : """
542 : vcat(s::H, r::H, ...) # H<:AbstractHermiteSpline
543 :
544 : Construct a new Hermite spline by concatenating `r, s, ...`.
545 :
546 : !!! warning "precondition"
547 : Same requirements as for [`append!`](@ref)!
548 :
549 : See also [`append!`](@ref)
550 : """
551 7 : function Base.vcat(r::H, s::H, args...) where {S, V, H<:AbstractHermiteSpline{S, V}}
552 7 : r = copy(r)
553 :
554 7 : append!(r, s)
555 :
556 5 : for h::H in args
557 4 : append!(r, h)
558 3 : end
559 :
560 4 : r
561 : end
|
