Skip to content

1D spline smoothing with knot locations specified #96

@Divisible8737

Description

@Divisible8737

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

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Type

    No type

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions