kbarbary / Dierckx.jl

Julia package for 1-d and 2-d splines
Other
157 stars 30 forks source link

1D spline smoothing with knot locations specified #96

Open Divisible8737 opened 5 months ago

Divisible8737 commented 5 months ago

Currently, the wrapper for Spline1D that allows for user-specified knot locations does not allow for smoothing: https://github.com/kbarbary/Dierckx.jl/blob/a5406e7e0e0901b91b7706cb42869e2ca08e1075/src/Dierckx.jl#L185

Instead, the smoothing parameter input is hardcoded to s=-1.0 here: https://github.com/kbarbary/Dierckx.jl/blob/a5406e7e0e0901b91b7706cb42869e2ca08e1075/src/Dierckx.jl#L234

and here: https://github.com/kbarbary/Dierckx.jl/blob/a5406e7e0e0901b91b7706cb42869e2ca08e1075/src/Dierckx.jl#L245

I propose the following nonbreaking modification:

function Spline1D(x::AbstractVector, y::AbstractVector,
    knots::AbstractVector;
    w::AbstractVector=ones(length(x)),
    k::Int=3, s::Real=-1.0, bc::AbstractString="nearest",
    periodic::Bool=false)
m = length(x)
length(y) == m || error("length of x and y must match")
length(w) == m || error("length of x and w must match")
m > k || error("k must be less than length(x)")
length(knots) <= m + k + 1 || error("length(knots) <= length(x) + k + 1 must hold")
first(x) < first(knots) || error("first(x) < first(knots) must hold")
last(x) > last(knots) || error("last(x) > last(knots) must hold")

# ensure inputs are of correct type
xin = convert(Vector{Float64}, x)
yin = convert(Vector{Float64}, y)
win = convert(Vector{Float64}, w)

# x knots
# (k+1) knots will be added on either end of interior knots.
n = length(knots) + 2(k + 1)
t = Vector{Float64}(undef, n)  # All knots
t[k+2:end-k-1] = knots

# outputs
c = Vector{Float64}(undef, n)
fp = Ref{Float64}(0)
ier = Ref{Int32}(0)

# workspace
lwrk = 0
if periodic
lwrk = m * (k + 1) + n*(8 + 5k)
else
lwrk = m * (k + 1) + n*(7 + 3k)
end
wrk = Vector{Float64}(undef, lwrk)
iwrk = Vector{Int32}(undef, n)

if !periodic
ccall((:curfit_, libddierckx), Nothing,
(Ref{Int32}, Ref{Int32},  # iopt, m
Ref{Float64}, Ref{Float64}, Ref{Float64},  # x, y, w
Ref{Float64}, Ref{Float64},  # xb, xe
Ref{Int32}, Ref{Float64},  # k, s
Ref{Int32}, Ref{Int32}, # nest, n
Ref{Float64}, Ref{Float64}, Ref{Float64},  # t, c, fp
Ref{Float64}, Ref{Int32}, Ref{Int32},  # wrk, lwrk, iwrk
Ref{Int32}),  # ier
-1, m, xin, yin, win, xin[1], xin[end], k, Float64(s),
n, n, t, c, fp, wrk, lwrk, iwrk, ier)
else
ccall((:percur_, libddierckx), Nothing,
(Ref{Int32}, Ref{Int32}, # iopt, m
Ref{Float64}, Ref{Float64}, Ref{Float64}, # x, y, w
Ref{Int32}, Ref{Float64}, # k, s
Ref{Int32}, Ref{Int32}, # nest, n
Ref{Float64}, Ref{Float64}, Ref{Float64}, # t, c, fp
Ref{Float64}, Ref{Int32}, Ref{Int32}, # wrk, lwrk, iwrk
Ref{Int32}), # ier
-1, m, xin, yin, win, k, Float64(s), n, n, t, c,
fp, wrk, lwrk, iwrk, ier)
end

ier[] <= 0 || error(_fit1d_messages[ier[]])
resize!(c, n - k - 1)

return Spline1D(t, c, k, _translate_bc(bc), fp[])
end