Report index
StatProfilerHTML.jl report
Generated on Mon, 01 Apr 2024 21:01:18
File source code
Line Exclusive Inclusive Code
1 module Dierckx
2
3 using Dierckx_jll
4
5 export Spline1D,
6 Spline2D,
7 ParametricSpline,
8 evaluate,
9 derivative,
10 integrate,
11 roots,
12 evalgrid,
13 get_knots,
14 get_coeffs,
15 get_residual
16
17 import Base: show, ==
18
19 # ----------------------------------------------------------------------------
20 # 1-d splines
21
22 const _fit1d_messages = Dict(
23 2=>
24 """A theoretically impossible result was found during the iteration
25 process for finding a smoothing spline with fp = s: s too small.
26 There is an approximation returned but the corresponding weighted sum
27 of squared residuals does not satisfy the condition abs(fp-s)/s <
28 tol.""",
29 3=>
30 """The maximal number of iterations maxit (set to 20 by the program)
31 allowed for finding a smoothing spline with fp=s has been reached: s
32 too small. There is an approximation returned but the corresponding
33 weighted sum of squared residuals does not satisfy the condition
34 abs(fp-s)/s < tol.""",
35 10=>
36 """Error on entry, no approximation returned. The following conditions
37 must hold:
38 1<=k<=5
39 x[1] < x[2] < ... < x[end]
40 w[i] > 0.0 for all i
41
42 Additionally, if spline knots are given:
43 length(xknots) <= length(x) + k + 1
44 x[1] < xknots[1] < xknots[k+2] < ... < xknots[end] < x[end]
45 The schoenberg-whitney conditions: there must be a subset of data points
46 xx[j] such that t[j] < xx[j] < t[j+k+1] for j=1,2,...,n-k-1""")
47
48
49 const _eval1d_messages = Dict(
50 1=>
51 """Input point out of range""",
52 10=>
53 """Invalid input data. The following conditions must hold:
54 length(x) != 0 and xb <= x[1] <= x[2] <= ... x[end] <= xe""")
55
56 _translate_bc(bc::AbstractString) = (bc == "extrapolate" ? 0 :
57 bc == "zero" ? 1 :
58 bc == "error" ? 2 :
59 bc == "nearest" ? 3 :
60 error("unknown boundary condition: \"$bc\""))
61 _translate_bc(bc::Int) = (bc == 0 ? "extrapolate" :
62 bc == 1 ? "zero" :
63 bc == 2 ? "error" :
64 bc == 3 ? "nearest" : "")
65
66 mutable struct Spline1D
67 t::Vector{Float64}
68 c::Vector{Float64}
69 k::Int
70 bc::Int
71 fp::Float64
72 wrk::Vector{Float64}
73 end
74
75 # add a constructor that automatically creates the `work` array
76 Spline1D(t, c, k, bc, fp) = Spline1D(t, c, k, bc, fp, Vector{Float64}(undef, length(t)))
77
78 get_knots(spl::Spline1D) = spl.t[spl.k+1:end-spl.k]
79 get_coeffs(spl::Spline1D) = spl.c[1:end-spl.k+1]
80 get_residual(spl::Spline1D) = spl.fp
81
82 function reallycompact(a::Vector)
83 io = IOBuffer()
84 io_compact = IOContext(io, :compact => true)
85 if length(a) <= 5
86 show(io, a)
87 else
88 write(io, "[")
89 show(io_compact, a[1])
90 write(io, ",")
91 show(io_compact, a[2])
92 write(io, " \u2026 ")
93 show(io_compact, a[end-1])
94 write(io, ",")
95 show(io_compact, a[end])
96 write(io, "]")
97 write(io, " ($(length(a)) elements)")
98 end
99 seekstart(io)
100 return read(io, String)
101 end
102
103 function ==(s1::Spline1D, s2::Spline1D)
104 s1.t == s2.t && s1.c == s2.c && s1.k == s2.k && s1.bc == s2.bc && s1.fp == s2.fp
105 end
106
107 function show(io::IO, spl::Spline1D)
108 print(io, """Spline1D(knots=$(reallycompact(get_knots(spl))), k=$(spl.k), extrapolation=\"$(_translate_bc(spl.bc))\", residual=$(spl.fp))""")
109 end
110
111 function Spline1D(x::AbstractVector, y::AbstractVector;
112 w::AbstractVector=ones(length(x)),
113 k::Int=3, s::Real=0.0, bc::AbstractString="nearest",
114 periodic::Bool=false)
115 m = length(x)
116 length(y) == m || error("length of x and y must match")
117 length(w) == m || error("length of x and w must match")
118 m > k || error("k must be less than length(x)")
119 1 <= k <= 5 || error("1 <= k = $k <= 5 must hold")
120
121 # ensure inputs are of correct type
122 xin = convert(Vector{Float64}, x)
123 yin = convert(Vector{Float64}, y)
124 win = convert(Vector{Float64}, w)
125
126 nest = 0
127 if periodic
128 nest = max(m + 2k, 2k + 3)
129 else
130 nest = max(m + k + 1, 2k + 3)
131 end
132
133 # outputs
134 n = Ref{Int32}(0)
135 t = Vector{Float64}(undef, nest)
136 c = Vector{Float64}(undef, nest)
137 fp = Ref{Float64}(0)
138 ier = Ref{Int32}(0)
139
140 # workspace
141 lwrk = 0
142 if periodic
143 lwrk = m * (k + 1) + nest*(8 + 5k)
144 else
145 lwrk = m * (k + 1) + nest*(7 + 3k)
146 end
147 wrk = Vector{Float64}(undef, lwrk)
148 iwrk = Vector{Int32}(undef, nest)
149
150 if !periodic
151 ccall((:curfit_, libddierckx), Nothing,
152 (Ref{Int32}, Ref{Int32}, # iopt, m
153 Ref{Float64}, Ref{Float64}, Ref{Float64}, # x, y, w
154 Ref{Float64}, Ref{Float64}, # xb, xe
155 Ref{Int32}, Ref{Float64}, # k, s
156 Ref{Int32}, Ref{Int32}, # nest, n
157 Ref{Float64}, Ref{Float64}, Ref{Float64}, # t, c, fp
158 Ref{Float64}, Ref{Int32}, Ref{Int32}, # wrk, lwrk, iwrk
159 Ref{Int32}), # ier
160 0, m, xin, yin, win, xin[1], xin[end], k, Float64(s),
161 nest, n, t, c, fp, wrk, lwrk, iwrk, ier)
162 else
163 ccall((:percur_, libddierckx), Nothing,
164 (Ref{Int32}, Ref{Int32}, # iopt, m
165 Ref{Float64}, Ref{Float64}, Ref{Float64}, # x, y, w
166 Ref{Int32}, Ref{Float64}, # k, s
167 Ref{Int32}, Ref{Int32}, # nest, n
168 Ref{Float64}, Ref{Float64}, Ref{Float64}, # t, c, fp
169 Ref{Float64}, Ref{Int32}, Ref{Int32}, # wrk, lwrk, iwrk
170 Ref{Int32}), # ier
171 0, m, xin, yin, win, k, Float64(s), nest, n, t, c,
172 fp, wrk, lwrk, iwrk, ier)
173 end
174
175 ier[] <= 0 || error(_fit1d_messages[ier[]])
176
177 # resize output arrays
178 resize!(t, n[])
179 resize!(c, n[] - k - 1)
180
181 return Spline1D(t, c, k, _translate_bc(bc), fp[])
182 end
183
184 # version with user-supplied knots
185 function Spline1D(x::AbstractVector, y::AbstractVector,
186 knots::AbstractVector;
187 w::AbstractVector=ones(length(x)),
188 k::Int=3, bc::AbstractString="nearest",
189 periodic::Bool=false)
190 m = length(x)
191 length(y) == m || error("length of x and y must match")
192 length(w) == m || error("length of x and w must match")
193 m > k || error("k must be less than length(x)")
194 length(knots) <= m + k + 1 || error("length(knots) <= length(x) + k + 1 must hold")
195 first(x) < first(knots) || error("first(x) < first(knots) must hold")
196 last(x) > last(knots) || error("last(x) > last(knots) must hold")
197
198 # ensure inputs are of correct type
199 xin = convert(Vector{Float64}, x)
200 yin = convert(Vector{Float64}, y)
201 win = convert(Vector{Float64}, w)
202
203 # x knots
204 # (k+1) knots will be added on either end of interior knots.
205 n = length(knots) + 2(k + 1)
206 t = Vector{Float64}(undef, n) # All knots
207 t[k+2:end-k-1] = knots
208
209 # outputs
210 c = Vector{Float64}(undef, n)
211 fp = Ref{Float64}(0)
212 ier = Ref{Int32}(0)
213
214 # workspace
215 lwrk = 0
216 if periodic
217 lwrk = m * (k + 1) + n*(8 + 5k)
218 else
219 lwrk = m * (k + 1) + n*(7 + 3k)
220 end
221 wrk = Vector{Float64}(undef, lwrk)
222 iwrk = Vector{Int32}(undef, n)
223
224 if !periodic
225 ccall((:curfit_, libddierckx), Nothing,
226 (Ref{Int32}, Ref{Int32}, # iopt, m
227 Ref{Float64}, Ref{Float64}, Ref{Float64}, # x, y, w
228 Ref{Float64}, Ref{Float64}, # xb, xe
229 Ref{Int32}, Ref{Float64}, # k, s
230 Ref{Int32}, Ref{Int32}, # nest, n
231 Ref{Float64}, Ref{Float64}, Ref{Float64}, # t, c, fp
232 Ref{Float64}, Ref{Int32}, Ref{Int32}, # wrk, lwrk, iwrk
233 Ref{Int32}), # ier
234 -1, m, xin, yin, win, xin[1], xin[end], k, -1.0,
235 n, n, t, c, fp, wrk, lwrk, iwrk, ier)
236 else
237 ccall((:percur_, libddierckx), Nothing,
238 (Ref{Int32}, Ref{Int32}, # iopt, m
239 Ref{Float64}, Ref{Float64}, Ref{Float64}, # x, y, w
240 Ref{Int32}, Ref{Float64}, # k, s
241 Ref{Int32}, Ref{Int32}, # nest, n
242 Ref{Float64}, Ref{Float64}, Ref{Float64}, # t, c, fp
243 Ref{Float64}, Ref{Int32}, Ref{Int32}, # wrk, lwrk, iwrk
244 Ref{Int32}), # ier
245 -1, m, xin, yin, win, k, -1.0, n, n, t, c,
246 fp, wrk, lwrk, iwrk, ier)
247 end
248
249 ier[] <= 0 || error(_fit1d_messages[ier[]])
250 resize!(c, n - k - 1)
251
252 return Spline1D(t, c, k, _translate_bc(bc), fp[])
253 end
254
255
256 function _evaluate(t::Vector{Float64}, c::Vector{Float64}, k::Int,
257 x::Vector{Float64}, bc::Int)
258 bc in (0, 1, 2, 3) || error("bc = $bc not in (0, 1, 2, 3)")
259 m = length(x)
260 xin = convert(Vector{Float64}, x)
261 y = Vector{Float64}(undef, m)
262 ier = Ref{Int32}(0)
263 ccall((:splev_, libddierckx), Nothing,
264 (Ref{Float64}, Ref{Int32},
265 Ref{Float64}, Ref{Int32},
266 Ref{Float64}, Ref{Float64}, Ref{Int32},
267 Ref{Int32}, Ref{Int32}),
268 t, length(t), c, k, xin, y, m, bc, ier)
269
270 ier[] == 0 || error(_eval1d_messages[ier[]])
271 return y
272 end
273
274
6 (10 %) samples spent in _evaluate
6 (100 %) (ex.), 6 (100 %) (incl.) when called from evaluate line 296
function _evaluate(t::Vector{Float64}, c::Vector{Float64}, k::Int,
275 x::Real, bc::Int)
276 bc in (0, 1, 2, 3) || error("bc = $bc not in (0, 1, 2, 3)")
277 y = Ref{Float64}(0)
278 ier = Ref{Int32}(0)
279 6 (10 %) 6 (10 %) ccall((:splev_, libddierckx), Nothing,
280 (Ref{Float64}, Ref{Int32},
281 Ref{Float64}, Ref{Int32},
282 Ref{Float64}, Ref{Float64}, Ref{Int32},
283 Ref{Int32}, Ref{Int32}),
284 t, length(t), c, k, Float64(x), y, 1, bc, ier)
285
286 ier[] == 0 || error(_eval1d_messages[ier[]])
287 return y[]
288 end
289
290
291 evaluate(spline::Spline1D, x::AbstractVector) =
292 _evaluate(spline.t, spline.c, spline.k,
293 convert(Vector{Float64}, x), spline.bc)
294
295
296 1 (2 %) 7 (12 %)
7 (12 %) samples spent in evaluate
1 (14 %) (ex.), 7 (100 %) (incl.) when called from Spline1D line 1112
6 (100 %) samples spent calling _evaluate
evaluate(spline::Spline1D, x::Real) =
297 _evaluate(spline.t, spline.c, spline.k, x, spline.bc)
298
299
300 function _derivative(t::Vector{Float64}, c::Vector{Float64}, k::Int,
301 x::Vector{Float64}, nu::Int, bc::Int, wrk::Vector{Float64})
302 (1 <= nu <= k) || error("order of derivative must be positive and less than or equal to spline order")
303 m = length(x)
304 n = length(t)
305 y = Vector{Float64}(undef, m)
306 ier = Ref{Int32}(0)
307 ccall((:splder_, libddierckx), Nothing,
308 (Ref{Float64}, Ref{Int32}, # t, n
309 Ref{Float64}, Ref{Int32}, # c, k
310 Ref{Int32}, # nu
311 Ref{Float64}, Ref{Float64}, Ref{Int32}, # x, y, m
312 Ref{Int32}, Ref{Float64}, Ref{Int32}), # e, wrk, ier
313 t, n, c, k, nu, x, y, m, bc, wrk, ier)
314 ier[] == 0 || error(_eval1d_messages[ier[]])
315 return y
316 end
317
318 function _derivative(t::Vector{Float64}, c::Vector{Float64}, k::Int,
319 x::Real, nu::Int, bc::Int, wrk::Vector{Float64})
320 (1 <= nu <= k) || error("order of derivative must be positive and less than or equal to spline order")
321 n = length(t)
322 y = Ref{Float64}(0)
323 ier = Ref{Int32}(0)
324 ccall((:splder_, libddierckx), Nothing,
325 (Ref{Float64}, Ref{Int32}, # t, n
326 Ref{Float64}, Ref{Int32}, # c, k
327 Ref{Int32}, # nu
328 Ref{Float64}, Ref{Float64}, Ref{Int32}, # x, y, m
329 Ref{Int32}, Ref{Float64}, Ref{Int32}), # e, wrk, ier
330 t, n, c, k, nu, Float64(x), y, 1, bc, wrk, ier)
331 ier[] == 0 || error(_eval1d_messages[ier[]])
332 return y[]
333 end
334
335 # TODO: should the function name be evaluate, derivative, or grad?
336 # or should it be integrated with evaluate, above?
337 # (problem with that: derivative doesn't accept bc="nearest")
338 # TODO: should `nu` be `d`?
339 derivative(spline::Spline1D, x::AbstractVector, nu::Int=1) =
340 _derivative(spline.t, spline.c, spline.k,
341 convert(Vector{Float64}, x), nu, spline.bc, spline.wrk)
342
343
344 derivative(spline::Spline1D, x::Real, nu::Int=1) =
345 _derivative(spline.t, spline.c, spline.k, x, nu, spline.bc, spline.wrk)
346
347
348 # TODO: deprecate this?
349 derivative(spline::Spline1D, x; nu::Int=1) = derivative(spline, x, nu)
350
351
352 function _integrate(t::Vector{Float64}, c::Vector{Float64}, k::Int,
353 a::Real, b::Real, wrk::Vector{Float64})
354 n = length(t)
355 ccall((:splint_, libddierckx), Float64,
356 (Ref{Float64}, Ref{Int32},
357 Ref{Float64}, Ref{Int32},
358 Ref{Float64}, Ref{Float64},
359 Ref{Float64}),
360 t, n, c, k,
361 Float64(a), Float64(b), wrk)
362 end
363
364 integrate(spline::Spline1D, a::Real, b::Real) =
365 _integrate(spline.t, spline.c, spline.k, a, b, spline.wrk)
366
367 # TODO roots for parametric splines
368 # note: default maxn in scipy.interpolate is 3 * (length(spline.t) - 7)
369 function roots(spline::Spline1D; maxn::Integer=8)
370 if spline.k != 3
371 error("root finding only supported for cubic splines (k=3)")
372 end
373 n = length(spline.t)
374 zeros = Vector{Float64}(undef, maxn)
375 m = Vector{Int32}(undef, 1)
376 ier = Vector{Int32}(undef, 1)
377 ccall((:sproot_, libddierckx), Nothing,
378 (Ref{Float64}, Ref{Int32}, # t, n
379 Ref{Float64}, Ref{Float64}, # c, zeros
380 Ref{Int32}, Ref{Int32}, # mest, m
381 Ref{Int32}), # ier
382 spline.t, n, spline.c, zeros,
383 maxn, m, ier)
384
385 if ier[1] == 0
386 return zeros[1:m[1]]
387 elseif ier[1] == 1
388 @warn("number of zeros exceeded maxn; only first maxn zeros returned")
389 return zeros
390 elseif ier[1] == 10
391 error("Invalid input data.")
392 else
393 error("unknown error code in sproot: $(ier[1])")
394 end
395 end
396
397 # ----------------------------------------------------------------------------
398 # parametric splines
399
400 mutable struct ParametricSpline
401 t::Vector{Float64}
402 c::Matrix{Float64}
403 k::Int
404 bc::Int
405 fp::Float64
406 wrk::Vector{Float64}
407 end
408
409 ParametricSpline(t, c, k, bc, fp) =
410 ParametricSpline(t, c, k, bc, fp, Vector{Float64}(undef, length(t)))
411
412 get_knots(spl::ParametricSpline) = spl.t[spl.k+1:end-spl.k]
413 get_coeffs(spl::ParametricSpline) = spl.c[:, 1:end-spl.k+1]
414 get_residual(spl::ParametricSpline) = spl.fp
415
416 function show(io::IO, spl::ParametricSpline)
417 print(io, """ParametricSpline(knots=$(reallycompact(get_knots(spl))), k=$(spl.k), extrapolation=\"$(_translate_bc(spl.bc))\", residual=$(spl.fp))""")
418 end
419
420 function ==(s1::ParametricSpline, s2::ParametricSpline)
421 s1.t == s2.t && s1.c == s2.c && s1.k == s2.k && s1.bc == s2.bc && s1.fp == s2.fp
422 end
423
424
425 function ParametricSpline(x::AbstractMatrix;
426 w::AbstractVector=ones(size(x, 2)),
427 k::Int=3, s::Real=0.0, bc::AbstractString="nearest",
428 periodic::Bool=false)
429 return _ParametricSpline(nothing, x, nothing, w, k, s, bc, periodic)
430 end
431
432 # version with user-supplied u
433 function ParametricSpline(u::AbstractVector, x::AbstractMatrix;
434 ub::Real=u[1], ue::Real=u[end],
435 w::AbstractVector=ones(size(x, 2)),
436 k::Int=3, s::Real=0.0, bc::AbstractString="nearest",
437 periodic::Bool=false)
438 return _ParametricSpline(u, x, nothing, w, k, s, bc, periodic)
439 end
440
441 # version with user-supplied knots
442 function ParametricSpline(x::AbstractMatrix, knots::AbstractVector;
443 w::AbstractVector=ones(size(x, 2)),
444 k::Int=3, bc::AbstractString="nearest",
445 periodic::Bool=false)
446 return _ParametricSpline(nothing, x, knots, w, k, -1.0, bc, periodic)
447 end
448
449 # version with user-supplied u and knots
450 function ParametricSpline(u::AbstractVector, x::AbstractMatrix,
451 knots::AbstractVector;
452 w::AbstractVector=ones(size(x, 2)),
453 k::Int=3, bc::AbstractString="nearest",
454 periodic::Bool=false)
455 return _ParametricSpline(u, x, knots, w, k, -1.0, bc, periodic)
456 end
457
458 function _ParametricSpline(u::Union{AbstractVector, Nothing}, x::AbstractMatrix,
459 knots::Union{AbstractVector, Nothing},
460 w::AbstractVector, k::Int, s::Real,
461 bc::AbstractString, periodic::Bool)
462 idim, m = size(x)
463 if periodic
464 x[:, 1] == x[:, end] || error("for periodic splines x[:,1] and x[:,end] must match")
465 end
466
467 length(w) == m || error("number of data points and length of w must match")
468 0 < idim < 11 || error("number of dimension must be between 1 and 10")
469 m > k || error("number of data points must be greater than k")
470 1 <= k <= 5 || error("1 <= k = $k <= 5 must hold")
471
472 local ipar, uin, ub, ue
473 if u != nothing
474 all(u[1:end-1] .< u[2:end]) || error("u[i] must be strictly increasing")
475 ipar = 1
476 uin = convert(Vector{Float64}, u)
477 ub = uin[1]
478 ue = uin[end]
479 else
480 ipar = 0
481 uin = Vector{Float64}(undef, size(x, 2))
482 ub = 0.0
483 ue = 1.0
484 end
485
486 local iopt, nest::Int, t
487 if knots != nothing
488 length(knots) <= m + k + 1 || error("length(knots) <= length(x) + k + 1 must hold")
489 first(x) < first(knots) || error("first(x) < first(knots) must hold")
490 last(x) > last(knots) || error("last(x) > last(knots) must hold")
491 iopt = -1
492 nest = length(knots) + 2(k + 1)
493 t = Vector{Float64}(undef, nest)
494 t[k+2:end-k-1] = knots
495 else
496 iopt = 0
497 nest = m + 2k
498 if s == 0
499 nest = periodic ? m + 2k : m + k + 1
500 end
501 nest = max(nest, 2k + 3)
502 t = Vector{Float64}(undef, nest)
503 end
504
505 xin = convert(Matrix{Float64}, x)
506 win = convert(Vector{Float64}, w)
507 n = Ref{Int32}(nest)
508 c = Vector{Float64}(undef, idim*nest)
509 fp = Ref{Float64}(0)
510 ier = Ref{Int32}(0)
511
512 local lwrk::Int
513 if periodic
514 lwrk = m*(k + 1) + nest*(7 + idim + 5k)
515 else
516 lwrk = m*(k + 1) + nest*(6 + idim + 3k)
517 end
518 wrk = Vector{Float64}(undef, lwrk)
519 iwrk = Vector{Int32}(undef, nest)
520
521 if !periodic
522 ccall((:parcur_, libddierckx), Nothing,
523 (Ref{Int32}, Ref{Int32}, # iopt, ipar
524 Ref{Int32}, Ref{Int32}, # idim, m
525 Ref{Float64}, Ref{Int32}, # u, mx
526 Ref{Float64}, Ref{Float64}, # x, w
527 Ref{Float64}, Ref{Float64}, # ub, ue
528 Ref{Int32}, Ref{Float64}, # k, s
529 Ref{Int32}, Ref{Int32}, # nest, n
530 Ref{Float64}, Ref{Int32}, # t, nc
531 Ref{Float64}, Ref{Float64}, # c, fp
532 Ref{Float64}, Ref{Int32}, # wrk, lwrk
533 Ref{Int32}, Ref{Int32}), #iwrk, ier
534 iopt, ipar,
535 idim, m,
536 uin, length(x),
537 xin, win,
538 ub, ue,
539 k, s,
540 nest, n,
541 t, length(c),
542 c, fp,
543 wrk, lwrk,
544 iwrk, ier)
545 else
546 ccall((:clocur_, libddierckx), Nothing,
547 (Ref{Int32}, Ref{Int32}, # iopt, ipar
548 Ref{Int32}, Ref{Int32}, # idim, m
549 Ref{Float64}, Ref{Int32}, # u, mx
550 Ref{Float64}, Ref{Float64}, # x, w
551 Ref{Int32}, Ref{Float64}, # k, s
552 Ref{Int32}, Ref{Int32}, # nest, n
553 Ref{Float64}, Ref{Int32}, # t, nc
554 Ref{Float64}, Ref{Float64}, # c, fp
555 Ref{Float64}, Ref{Int32}, # wrk, lwrk
556 Ref{Int32}, Ref{Int32}), #iwrk, ier
557 iopt, ipar,
558 idim, m,
559 uin, length(x),
560 xin, win,
561 k, s,
562 nest, n,
563 t, length(c),
564 c, fp,
565 wrk, lwrk,
566 iwrk, ier)
567 end
568
569 ier[] <= 0 || error(_fit1d_messages[ier[]])
570
571 resize!(t, n[])
572 c = [c[n[]*(j-1) + i] for j=1:idim, i=1:n[]-k-1]
573
574 return ParametricSpline(t, c, k, _translate_bc(bc), fp[])
575 end
576
577 _evaluate(t::Vector{Float64}, c::Matrix{Float64}, k::Int,
578 x::Vector{Float64}, bc::Int) =
579 mapslices(v -> _evaluate(t, v, k, x, bc), c, dims=[2])
580
581 _evaluate(t::Vector{Float64}, c::Matrix{Float64}, k::Int,
582 x::Real, bc::Int) =
583 vec(mapslices(v -> _evaluate(t, v, k, x, bc), c, dims=[2]))
584
585 function evaluate(spline::ParametricSpline, x::AbstractVector)
586 xin = convert(Vector{Float64}, x)
587 _evaluate(spline.t, spline.c, spline.k, xin, spline.bc)
588 end
589
590 evaluate(spline::ParametricSpline, x::Real) =
591 _evaluate(spline.t, spline.c, spline.k, x, spline.bc)
592
593 _derivative(t::Vector{Float64}, c::Matrix{Float64}, k::Int,
594 x::Vector{Float64}, nu::Int, bc::Int, wrk::Vector{Float64}) =
595 mapslices(v -> _derivative(t, v, k, x, nu, bc, wrk), c, dims=[2])
596
597 _derivative(t::Vector{Float64}, c::Matrix{Float64}, k::Int,
598 x::Real, nu::Int, bc::Int, wrk::Vector{Float64}) =
599 vec(mapslices(v -> _derivative(t, v, k, x, nu, bc, wrk), c, dims=[2]))
600
601 derivative(spline::ParametricSpline, x::AbstractVector, nu::Int=1) =
602 _derivative(spline.t, spline.c, spline.k,
603 convert(Vector{Float64}, x), nu, spline.bc, spline.wrk)
604
605 # TODO: deprecate this
606 derivative(spline::ParametricSpline, x; nu::Int=1) = derivative(spline, x, nu)
607
608 derivative(spline::ParametricSpline, x::Real, nu::Int=1) =
609 _derivative(spline.t, spline.c, spline.k, x, nu, spline.bc, spline.wrk)
610
611 _integrate(t::Vector{Float64}, c::Matrix{Float64}, k::Int,
612 a::Real, b::Real, wrk::Vector{Float64}) =
613 vec(mapslices(v -> _integrate(t, v, k, a, b, wrk), c, dims=[2]))
614
615 integrate(spline::ParametricSpline, a::Real, b::Real) =
616 _integrate(spline.t, spline.c, spline.k, a, b, spline.wrk)
617
618 # ----------------------------------------------------------------------------
619 # 2-d splines
620
621 # NOTE REGARDING ARGUMENT ORDER: In the "grid" version of the Spline2D
622 # constructor and evaluators, the fortran functions expects z to have
623 # shape (my, mx), but we'd rather have x be the fast axis in z. So,
624 # in the ccall()s in these methods, all the x and y related inputs are
625 # swapped with regard to what the Fortran documentation says.
626
627 const _fit2d_messages = Dict(
628 -3=>
629
630 """The coefficients of the spline returned have been computed as the
631 minimal norm least-squares solution of a (numerically) rank deficient
632 system (deficiency=%i). If deficiency is large, the results may be
633 inaccurate. Deficiency may strongly depend on the value of eps.""",
634
635 1=>
636
637 """The required storage space exceeds the available storage space:
638 nxest or nyest too small, or s too small. Try increasing s.""",
639 # The weighted least-squares spline corresponds to the current set of knots.
640
641 2=>
642
643 """A theoretically impossible result was found during the iteration
644 process for finding a smoothing spline with fp = s: s too small or
645 badly chosen eps. Weighted sum of squared residuals does not satisfy
646 abs(fp-s)/s < tol.""",
647
648 3=>
649
650 """the maximal number of iterations maxit (set to 20 by the program)
651 allowed for finding a smoothing spline with fp=s has been reached: s
652 too small. Weighted sum of squared residuals does not satisfy
653 abs(fp-s)/s < tol.""",
654
655 4=>
656
657 """No more knots can be added because the number of b-spline
658 coefficients (nx-kx-1)*(ny-ky-1) already exceeds the number of data
659 points m: either s or m too small. The weighted least-squares spline
660 corresponds to the current set of knots.""",
661
662 5=>
663
664 """No more knots can be added because the additional knot would
665 (quasi) coincide with an old one: s too small or too large a weight to
666 an inaccurate data point. The weighted least-squares spline
667 corresponds to the current set of knots.""",
668
669 10=>
670
671 """Error on entry, no approximation returned. The following conditions
672 must hold:
673 xb<=x[i]<=xe, yb<=y[i]<=ye, w[i]>0, i=0..m-1
674 If iopt==-1, then
675 xb<tx[kx+1]<tx[kx+2]<...<tx[nx-kx-2]<xe
676 yb<ty[ky+1]<ty[ky+2]<...<ty[ny-ky-2]<ye""")
677
678 const _eval2d_message = (
679 """Invalid input data. Restrictions:
680 length(x) != 0, length(y) != 0
681 x[i-1] <= x[i] for i=2,...,length(x)
682 y[j-1] <= y[j] for j=2,...,length(y)
683 """)
684
685 mutable struct Spline2D
686 tx::Vector{Float64}
687 ty::Vector{Float64}
688 c::Vector{Float64}
689 kx::Int
690 ky::Int
691 fp::Float64
692 end
693
694 get_knots(spl::Spline2D) = (spl.tx[spl.kx+1:end-spl.kx],
695 spl.ty[spl.ky+1:end-spl.ky])
696 get_residual(spl::Spline2D) = spl.fp
697
698 function ==(s1::Spline2D, s2::Spline2D)
699 s1.tx == s2.tx && s1.ty == s2.ty && s1.c == s2.c && s1.kx == s2.kx && s1.ky == s2.ky && s1.fp == s2.fp
700 end
701
702 # Helper functions for calculating required size of work arrays in surfit.
703 # Note that x and y here are as in the Fortran documentation.
704 # These are translated from scipy/interpolate/src/fitpack.pyf
705 function calc_surfit_lwrk1(m, kx, ky, nxest, nyest)
706 u = nxest - kx - 1
707 v = nyest - ky - 1
708 km = max(kx, ky) + 1
709 ne = max(nxest, nyest)
710 bx = kx*v + ky + 1
711 by = ky*u + kx + 1
712 b1 = b2 = 0
713 if (bx<=by)
714 b1 = bx
715 b2 = bx + v - ky
716 else
717 b1 = by
718 b2 = by + u - kx
719 end
720 return u*v*(2 + b1 + b2) + 2*(u+v+km*(m+ne)+ne-kx-ky) + b2 + 1
721 end
722
723 function calc_surfit_lwrk2(m, kx, ky, nxest, nyest)
724 u = nxest - kx - 1
725 v = nyest - ky - 1
726 bx = kx * v + ky + 1
727 by = ky * u + kx + 1
728 b2 = (bx <= by ? bx + v - ky : by + u - kx)
729 return u * v * (b2 + 1) + b2
730 end
731
732 # Construct spline from unstructured data
733 function Spline2D(x::AbstractVector, y::AbstractVector, z::AbstractVector;
734 w::AbstractVector=ones(length(x)), kx::Int=3, ky::Int=3,
735 s::Real=0.0)
736
737 # array sizes
738 m = length(x)
739 (length(y) == length(z) == m) || error("lengths of x, y, z must match")
740 (length(w) == m) || error("length of w must match other inputs")
741
742 nxest = max(kx+1+ceil(Int,sqrt(m/2)), 2*(kx+1))
743 nyest = max(ky+1+ceil(Int,sqrt(m/2)), 2*(ky+1))
744 nmax = max(nxest, nyest)
745
746 eps = 1.0e-16
747
748 # bounds
749 xb = minimum(x)
750 xe = maximum(x)
751 yb = minimum(y)
752 ye = maximum(y)
753
754 # ensure arrays are of correct type
755 xin = convert(Vector{Float64}, x)
756 yin = convert(Vector{Float64}, y)
757 zin = convert(Vector{Float64}, z)
758 win = convert(Vector{Float64}, w)
759
760 # return values
761 nx = Ref{Int32}()
762 tx = Vector{Float64}(undef, nxest)
763 ny = Ref{Int32}()
764 ty = Vector{Float64}(undef, nyest)
765 c = Vector{Float64}(undef, (nxest-kx-1) * (nyest-ky-1))
766 fp = Ref{Float64}()
767 ier = Ref{Int32}()
768
769 # work arrays
770 # Note: in lwrk1 and lwrk2, x and y are swapped on purpose.
771 lwrk1 = calc_surfit_lwrk1(m, ky, kx, nyest, nxest)
772 lwrk2 = calc_surfit_lwrk2(m, ky, kx, nyest, nxest)
773 kwrk = m + (nxest - 2*kx - 1) * (nyest - 2*ky - 1)
774 wrk1 = Vector{Float64}(undef, lwrk1)
775 wrk2 = Vector{Float64}(undef, lwrk2)
776 iwrk = Vector{Int32}(undef, kwrk)
777
778 ccall((:surfit_, libddierckx), Nothing,
779 (Ref{Int32}, Ref{Int32}, # iopt, m
780 Ref{Float64}, Ref{Float64}, # y, x
781 Ref{Float64}, Ref{Float64}, # z, w
782 Ref{Float64}, Ref{Float64}, # yb, ye
783 Ref{Float64}, Ref{Float64}, # xb, xe
784 Ref{Int32}, Ref{Int32}, Ref{Float64}, # ky, kx, s
785 Ref{Int32}, Ref{Int32}, Ref{Int32}, # nyest, nxest, nmax
786 Ref{Float64}, # eps
787 Ref{Int32}, Ref{Float64}, # ny, tx
788 Ref{Int32}, Ref{Float64}, # nx, tx
789 Ref{Float64}, Ref{Float64}, # c, fp
790 Ref{Float64}, Ref{Int32}, # wrk1, lwrk1
791 Ref{Float64}, Ref{Int32}, # wrk2, lwrk2
792 Ref{Int32}, Ref{Int32}, Ref{Int32}), # iwrk, kwrk, ier
793 0, m, yin, xin, zin, win, yb, ye, xb, xe, ky, kx,
794 Float64(s), nyest, nxest, nmax, eps, ny, ty, nx, tx,
795 c, fp, wrk1, lwrk1, wrk2, lwrk2, iwrk, kwrk, ier)
796
797 while ier[] > 10
798 # lwrk2 is too small, i.e., there is not enough workspace
799 # for computing the minimal least-squares solution of a rank
800 # deficient system of linear equations. ier gives the
801 # requested value for lwrk2. Rerun with that value in "continue"
802 # mode, with iopt = 1.
803 lwrk2 = ier[]
804 resize!(wrk2, lwrk2)
805 ccall((:surfit_, libddierckx), Nothing,
806 (Ref{Int32}, Ref{Int32}, # iopt, m
807 Ref{Float64}, Ref{Float64}, # y, x
808 Ref{Float64}, Ref{Float64}, # z, w
809 Ref{Float64}, Ref{Float64}, # yb, ye
810 Ref{Float64}, Ref{Float64}, # xb, xe
811 Ref{Int32}, Ref{Int32}, Ref{Float64}, # ky, kx, s
812 Ref{Int32}, Ref{Int32}, Ref{Int32}, # nyest, nxest, nmax
813 Ref{Float64}, # eps
814 Ref{Int32}, Ref{Float64}, # ny, tx
815 Ref{Int32}, Ref{Float64}, # nx, tx
816 Ref{Float64}, Ref{Float64}, # c, fp
817 Ref{Float64}, Ref{Int32}, # wrk1, lwrk1
818 Ref{Float64}, Ref{Int32}, # wrk2, lwrk2
819 Ref{Int32}, Ref{Int32}, Ref{Int32}), # iwrk, kwrk, ier
820 1, m, yin, xin, zin, win, yb, ye, xb, xe, ky, kx, s,
821 nyest, nxest, nmax, eps, ny, ty, nx, tx, c, fp,
822 wrk1, lwrk1, wrk2, lwrk2, iwrk, kwrk, ier)
823 end
824
825 if (ier[] == 0 || ier[] == -1 || ier[] == -2)
826 # good values, pass.
827 elseif ier[] < -2
828 @warn("""
829 The coefficients of the spline returned have been
830 computed as the minimal norm least-squares solution of a
831 (numerically) rank deficient system. The rank is $(-ier[]).
832 The rank deficiency is $((nx[]-kx-1)*(ny[]-ky-1)+ier[]).
833 Especially if the rank deficiency is large the results may
834 be inaccurate.""")
835 # "The results could also seriously depend on the value of
836 # eps" (not in message because eps is currently not an input)
837 else
838 error(_fit2d_messages[ier[]])
839 end
840
841 # Resize output arrays to the size actually used.
842 resize!(tx, nx[])
843 resize!(ty, ny[])
844 resize!(c, (nx[] - kx - 1) * (ny[] - ky - 1))
845
846 return Spline2D(tx, ty, c, kx, ky, fp[])
847 end
848
849
850
851 # Construct spline from data on a grid.
852 function Spline2D(x::AbstractVector, y::AbstractVector, z::AbstractMatrix;
853 kx::Int=3, ky::Int=3, s::Real=0.0)
854 mx = length(x)
855 my = length(y)
856 @assert size(z, 1) == mx && size(z, 2) == my
857
858 mx > kx || error("length(x) must be greater than kx")
859 my > ky || error("length(y) must be greater than ky")
860
861 # Bounds
862 xb = x[1]
863 xe = x[end]
864 yb = y[1]
865 ye = y[end]
866 nxest = mx+kx+1
867 nyest = my+ky+1
868
869 # ensure arrays are of correct type
870 xin = convert(Vector{Float64}, x)
871 yin = convert(Vector{Float64}, y)
872 zin = convert(Matrix{Float64}, z)
873
874 # Return values
875 nx = Ref{Int32}()
876 tx = Vector{Float64}(undef, nxest)
877 ny = Ref{Int32}()
878 ty = Vector{Float64}(undef, nyest)
879 c = Vector{Float64}(undef, (nxest-kx-1) * (nyest-ky-1))
880 fp = Ref{Float64}()
881 ier = Ref{Int32}()
882
883 # Work arrays.
884 # Note that in lwrk, x and y are swapped with respect to the Fortran
885 # documentation. See "NOTE REGARDING ARGUMENT ORDER" above.
886 lwrk = (4 + nyest * (mx+2*ky+5) + nxest * (2*kx+5) +
887 my*(ky+1) + mx*(kx+1) + max(mx, nyest))
888 wrk = Vector{Float64}(undef, lwrk)
889 kwrk = 3 + mx + my + nxest + nyest
890 iwrk = Vector{Int32}(undef, kwrk)
891
892 ccall((:regrid_, libddierckx), Nothing,
893 (Ref{Int32}, # iopt
894 Ref{Int32}, Ref{Float64}, # my, y
895 Ref{Int32}, Ref{Float64}, # mx, x
896 Ref{Float64}, # z
897 Ref{Float64}, Ref{Float64}, # yb, ye
898 Ref{Float64}, Ref{Float64}, # xb, xe
899 Ref{Int32}, Ref{Int32}, Ref{Float64}, # ky, kx, s
900 Ref{Int32}, Ref{Int32}, # nyest, nxest
901 Ref{Int32}, Ref{Float64}, # ny, ty
902 Ref{Int32}, Ref{Float64}, # nx, tx
903 Ref{Float64}, Ref{Float64}, # c, fp
904 Ref{Float64}, Ref{Int32}, # wrk, lwrk
905 Ref{Int32}, Ref{Int32}, # iwrk, lwrk
906 Ref{Int32}), # ier
907 0f0, my, yin, mx, xin, zin, yb, ye, xb, xe, ky, kx,
908 Float64(s), nyest, nxest, ny, ty, nx, tx, c, fp,
909 wrk, lwrk, iwrk, kwrk, ier)
910
911 if !(ier[] == 0 || ier[] == -1 || ier[] == -2)
912 error(_fit2d_messages[ier[]])
913 end
914
915 # Resize output arrays to the size actually used.
916 resize!(tx, nx[])
917 resize!(ty, ny[])
918 resize!(c, (nx[] - kx - 1) * (ny[] - ky - 1))
919
920 return Spline2D(tx, ty, c, kx, ky, fp[])
921 end
922
923
924 # Evaluate spline at individual points
925 function evaluate(spline::Spline2D, x::AbstractVector, y::AbstractVector)
926 m = length(x)
927 @assert length(y) == m
928
929 xin = convert(Vector{Float64}, x)
930 yin = convert(Vector{Float64}, y)
931
932 ier = Ref{Int32}()
933 lwrk = spline.kx + spline.ky + 2
934 wrk = Vector{Float64}(undef, lwrk)
935 z = Vector{Float64}(undef, m)
936
937 ccall((:bispeu_, libddierckx), Nothing,
938 (Ref{Float64}, Ref{Int32}, # ty, ny
939 Ref{Float64}, Ref{Int32}, # tx, nx
940 Ref{Float64}, # c
941 Ref{Int32}, Ref{Int32}, # ky, kx
942 Ref{Float64}, Ref{Float64}, Ref{Float64}, # y, x, z
943 Ref{Int32}, # m
944 Ref{Float64}, Ref{Int32}, Ref{Int32}), # wrk, lwrk, ier
945 spline.ty, length(spline.ty),
946 spline.tx, length(spline.tx),
947 spline.c, spline.ky, spline.kx, yin, xin, z, m,
948 wrk, lwrk, ier)
949
950 ier[] == 0 || error(_eval2d_message)
951
952 return z
953 end
954
955 function evaluate(spline::Spline2D, x::Real, y::Real)
956 ier = Ref{Int32}()
957 lwrk = spline.kx + spline.ky + 2
958 wrk = Vector{Float64}(undef, lwrk)
959 z = Ref{Float64}()
960 ccall((:bispeu_, libddierckx), Nothing,
961 (Ref{Float64}, Ref{Int32}, # ty, ny
962 Ref{Float64}, Ref{Int32}, # tx, nx
963 Ref{Float64}, # c
964 Ref{Int32}, Ref{Int32}, # ky, kx
965 Ref{Float64}, Ref{Float64}, Ref{Float64}, # y, x, z
966 Ref{Int32}, # m
967 Ref{Float64}, Ref{Int32}, Ref{Int32}), # wrk, lwrk, ier
968 spline.ty, length(spline.ty),
969 spline.tx, length(spline.tx),
970 spline.c, spline.ky, spline.kx, y, x, z, 1,
971 wrk, lwrk, ier)
972
973 ier[] == 0 || error(_eval2d_message)
974 return z[]
975 end
976
977 # Evaluate spline on the grid spanned by the input arrays.
978 function evalgrid(spline::Spline2D, x::AbstractVector, y::AbstractVector)
979 mx = length(x)
980 my = length(y)
981
982 xin = convert(Vector{Float64}, x)
983 yin = convert(Vector{Float64}, y)
984
985 lwrk = mx*(spline.kx + 1) + my*(spline.ky + 1)
986 wrk = Vector{Float64}(undef, lwrk)
987 kwrk = mx + my
988 iwrk = Vector{Int32}(undef, kwrk)
989 ier = Ref{Int32}()
990 z = Matrix{Float64}(undef, mx, my)
991
992 ccall((:bispev_, libddierckx), Nothing,
993 (Ref{Float64}, Ref{Int32}, # ty, ny
994 Ref{Float64}, Ref{Int32}, # tx, nx
995 Ref{Float64}, # c
996 Ref{Int32}, Ref{Int32}, # ky, kx
997 Ref{Float64}, Ref{Int32}, # y, my
998 Ref{Float64}, Ref{Int32}, # x, mx
999 Ref{Float64}, # z
1000 Ref{Float64}, Ref{Int32}, # wrk, lwrk
1001 Ref{Int32}, Ref{Int32}, # iwrk, kwrk
1002 Ref{Int32}), # ier
1003 spline.ty, length(spline.ty),
1004 spline.tx, length(spline.tx),
1005 spline.c, spline.ky, spline.kx, yin, my, xin, mx, z,
1006 wrk, lwrk, iwrk, kwrk, ier)
1007
1008 ier[] == 0 || error(_eval2d_message)
1009
1010 return z
1011 end
1012
1013 function _derivative(tx::Vector{Float64}, ty::Vector{Float64}, c::Vector{Float64},
1014 kx::Int, ky::Int, nux::Int, nuy::Int, x::Vector{Float64}, y::Vector{Float64},
1015 wrk::Vector{Float64}, iwrk::Vector{Float64})
1016 (0 <= nux < kx) || error("order of derivative must be positive and less than the spline order")
1017 (0 <= nuy < ky) || error("order of derivative must be positive and less than the spline order")
1018 mx = length(x)
1019 my = length(y)
1020 nx = length(tx)
1021 ny = length(ty)
1022 lwrk = length(wrk)
1023 lwrkmin = mx * (kx + 1 - nux) + my * (ky + 1 - nuy) + (nx - kx - 1) * (ny - ky - 1)
1024 lwrk >= lwrkmin || error("Length of wrk must be at least $lwrkmin")
1025 kwrk = length(iwrk)
1026 kwrk >= mx + my || error("length of iwrk must be greater than or equal to length(x) + length(y) = $(mx + my)")
1027 z = Vector{Float64}(undef, mx * my)
1028 ier = Ref{Int32}(0)
1029 # the order of x and y are switched compared to the fortran implementation
1030 # here x refers to rows and y to columns
1031 ccall((:parder_, libddierckx), Nothing,
1032 (Ref{Float64}, Ref{Int32}, # ty, ny
1033 Ref{Float64}, Ref{Int32}, # tx, nx
1034 Ref{Float64}, Ref{Int32}, Ref{Int32}, # c, ky, kx
1035 Ref{Int32}, Ref{Int32}, # nuy, nux
1036 Ref{Float64}, Ref{Int32}, Ref{Float64}, Ref{Int32}, Ref{Float64}, # y, my, x, mx, z
1037 Ref{Float64}, Ref{Int32}, Ref{Float64}, Ref{Int32}, # wrk, lwrk, iwrk, kwrk
1038 Ref{Int32}), # ier
1039 ty, ny, tx, nx, c, ky, kx, nuy, nux, y, my, x, mx, z, wrk, lwrk, iwrk, kwrk, ier)
1040 ier[] == 0 || error(_eval2d_message)
1041 return reshape(z, mx, my)
1042 end
1043
1044 function derivative(spline::Spline2D, x::AbstractVector, y::AbstractVector, nux::Int = 1, nuy::Int = 1)
1045 mx = length(x)
1046 my = length(y)
1047 nx = length(spline.tx)
1048 ny = length(spline.ty)
1049
1050 kx = spline.kx
1051 ky = spline.ky
1052
1053 lwrkmin = mx * (kx + 1 - nux) + my * (ky + 1 - nuy) + (nx - kx - 1) * (ny - ky - 1)
1054 lwrk = lwrkmin
1055 wrk = Vector{Float64}(undef, lwrk)
1056 kwrk = mx + my
1057 iwrk = Vector{Float64}(undef, kwrk)
1058
1059 _derivative(spline.tx, spline.ty, spline.c,
1060 spline.kx, spline.ky, nux, nuy,
1061 convert(Vector{Float64}, x),
1062 convert(Vector{Float64}, y),
1063 wrk, iwrk)
1064 end
1065
1066 function derivative(spline::Spline2D, x::Real, y::AbstractVector, nux::Int = 1, nuy::Int = 1)
1067 z = derivative(spline, [Float64(x)], y, nux, nuy)
1068 vec(z)
1069 end
1070
1071 function derivative(spline::Spline2D, x::AbstractVector, y::Real, nux::Int = 1, nuy::Int = 1)
1072 z = derivative(spline, x, [Float64(y)], nux, nuy)
1073 vec(z)
1074 end
1075 function derivative(spline::Spline2D, x::Real, y::Real, nux::Int = 1, nuy::Int = 1)
1076 z = derivative(spline, [Float64(x)], [Float64(y)], nux, nuy)
1077 z[]
1078 end
1079
1080 # TODO: deprecate this?
1081 function derivative(spline::Spline2D, x, y; nux::Int = 1, nuy::Int = 1)
1082 derivative(spline, x, y, nux, nuy)
1083 end
1084
1085 # 2D integration
1086 function integrate(spline::Spline2D, xb::Real, xe::Real, yb::Real, ye::Real)
1087 nx = length(spline.tx)
1088 ny = length(spline.ty)
1089
1090 kx = spline.kx
1091 ky = spline.ky
1092
1093 wrk = Vector{Float64}(undef, nx + ny -kx - ky -2)
1094 ccall((:dblint_, libddierckx), Float64,
1095 (Ref{Float64}, Ref{Int32}, # tx, nx
1096 Ref{Float64}, Ref{Int32}, # ty, ny
1097 Ref{Float64}, Ref{Int32}, # c, kx
1098 Ref{Int32}, #ky
1099 Ref{Float64}, Ref{Float64}, # xb, xe
1100 Ref{Float64}, Ref{Float64}, # yb, ye
1101 Ref{Float64}), # wrk
1102 spline.tx, nx,
1103 spline.ty, ny,
1104 spline.c, spline.kx,
1105 spline.ky,
1106 xb, xe,
1107 yb, ye,
1108 wrk)
1109 end
1110
1111 # call synonyms for evaluate():
1112 7 (12 %)
7 (12 %) samples spent in Spline1D
3 (43 %) (incl.) when called from calc_aero_forces! line 303
3 (43 %) (incl.) when called from calc_aero_forces! line 301
1 (14 %) (incl.) when called from calc_aero_forces! line 302
7 (100 %) samples spent calling evaluate
(spl::Spline1D)(x::Real) = evaluate(spl, x)
1113 (spl::Spline1D)(x::AbstractVector) = evaluate(spl, x)
1114 (spl::ParametricSpline)(x::Real) = evaluate(spl, x)
1115 (spl::ParametricSpline)(x::AbstractVector) = evaluate(spl, x)
1116 (spl::Spline2D)(x::Real, y::Real) = evaluate(spl, x, y)
1117 (spl::Spline2D)(x::AbstractVector, y::AbstractVector) =
1118 evaluate(spl, x, y)
1119
1120 end # module