Skip to content

1D spline smoothing with knot locations specified #96

Description

@matthewferrara

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

Activity

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

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