Smooth arclength interpolation second order derivatives

[yolhan83]

Hello, thank you for added this functionality, it seems like the second order derivative is not working for now

using DataInterpolations
t = 1:10
u = [rand(3) for _ in t]
interp = QuadraticInterpolation(u,t)
curve = SmoothArcLengthInterpolation(interp)
DataInterpolations.derivative(curve,1.0,1) #works
DataInterpolations.derivative(curve,1.0,2) #do not work

with the error,

1-element ExceptionStack:
LoadError: MethodError: no method matching Float64(::ForwardDiff.Dual{ForwardDiff.Tag{DataInterpolations.var"#56#57"{SmoothArcLengthInterpolation{Matrix{Float64}, Vector{Float64}, Nothing, Float64, Float64, QuadraticInterpolation{Vector{Vector{Float64}}, UnitRange{Int64}, Nothing, Vector{Vector{Float64}}, Vector{Float64}}, Float64}, FindFirstFunctions.Guesser{Vector{Float64}}}, Float64}, Float64, 1})
The type `Float64` exists, but no method is defined for this combination of argument types when trying to construct it.

Closest candidates are:
  (::Type{T})(::Real, ::RoundingMode) where T<:AbstractFloat
   @ Base rounding.jl:265
  (::Type{T})(::T) where T<:Number
   @ Core boot.jl:900
  Float64(::Float32)
   @ Base float.jl:341
  ...

Stacktrace:
  [1] convert(::Type{Float64}, x::ForwardDiff.Dual{ForwardDiff.Tag{DataInterpolations.var"#56#57"{SmoothArcLengthInterpolation{Matrix{Float64}, Vector{Float64}, Nothing, Float64, Float64, QuadraticInterpolation{Vector{Vector{Float64}}, UnitRange{Int64}, Nothing, Vector{Vector{Float64}}, Vector{Float64}}, Float64}, FindFirstFunctions.Guesser{Vector{Float64}}}, Float64}, Float64, 1})
    @ Base .\number.jl:7
  [2] setindex!(A::Vector{Float64}, x::ForwardDiff.Dual{ForwardDiff.Tag{DataInterpolations.var"#56#57"{SmoothArcLengthInterpolation{Matrix{Float64}, Vector{Float64}, Nothing, Float64, Float64, QuadraticInterpolation{Vector{Vector{Float64}}, UnitRange{Int64}, Nothing, Vector{Vector{Float64}}, Vector{Float64}}, Float64}, FindFirstFunctions.Guesser{Vector{Float64}}}, Float64}, Float64, 1}, i::Int64)
    @ Base .\array.jl:987
  [3] macro expansion
    @ .\broadcast.jl:973 [inlined]
  [4] macro expansion
    @ .\simdloop.jl:77 [inlined]
  [5] copyto!
    @ .\broadcast.jl:972 [inlined]
  [6] copyto!
    @ .\broadcast.jl:925 [inlined]
  [7] materialize!
    @ .\broadcast.jl:883 [inlined]
  [8] materialize!
    @ .\broadcast.jl:880 [inlined]
  [9] _derivative(A::SmoothArcLengthInterpolation{Matrix{Float64}, Vector{Float64}, Nothing, Float64, Float64, QuadraticInterpolation{Vector{Vector{Float64}}, UnitRange{Int64}, Nothing, Vector{Vector{Float64}}, Vector{Float64}}, Float64}, t::ForwardDiff.Dual{ForwardDiff.Tag{DataInterpolations.var"#56#57"{SmoothArcLengthInterpolation{Matrix{Float64}, Vector{Float64}, Nothing, Float64, Float64, QuadraticInterpolation{Vector{Vector{Float64}}, UnitRange{Int64}, Nothing, Vector{Vector{Float64}}, Vector{Float64}}, Float64}, FindFirstFunctions.Guesser{Vector{Float64}}}, Float64}, Float64, 1}, iguess::FindFirstFunctions.Guesser{Vector{Float64}})
    @ DataInterpolations C:\Users\yolha\.julia\packages\DataInterpolations\M6B1l\src\derivatives.jl:338
 [10] #56
    @ C:\Users\yolha\.julia\packages\DataInterpolations\M6B1l\src\derivatives.jl:11 [inlined]
 [11] derivative
    @ C:\Users\yolha\.julia\packages\ForwardDiff\UBbGT\src\derivative.jl:14 [inlined]
 [12] derivative(A::SmoothArcLengthInterpolation{Matrix{Float64}, Vector{Float64}, Nothing, Float64, Float64, QuadraticInterpolation{Vector{Vector{Float64}}, UnitRange{Int64}, Nothing, Vector{Vector{Float64}}, Vector{Float64}}, Float64}, t::Float64, order::Int64)
    @ DataInterpolations C:\Users\yolha\.julia\packages\DataInterpolations\M6B1l\src\derivatives.jl:9
 [13] top-level scope
    @ c:\Users\yolha\Desktop\juju_tests\Tests\main.jl:85
 [14] eval
    @ .\boot.jl:430 [inlined]
 [15] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
    @ Base .\loading.jl:2734
 [16] invokelatest(::Any, ::Any, ::Vararg{Any}; kwargs::@Kwargs{})
    @ Base .\essentials.jl:1055
 [17] invokelatest(::Any, ::Any, ::Vararg{Any})
    @ Base .\essentials.jl:1052
 [18] inlineeval(m::Module, code::String, code_line::Int64, code_column::Int64, file::String; softscope::Bool)
    @ VSCodeServer c:\Users\yolha\.vscode\extensions\julialang.language-julia-1.138.1\scripts\packages\VSCodeServer\src\eval.jl:271
 [19] (::VSCodeServer.var"#69#74"{Bool, Bool, Bool, Module, String, Int64, Int64, String, VSCodeServer.ReplRunCodeRequestParams})()
    @ VSCodeServer c:\Users\yolha\.vscode\extensions\julialang.language-julia-1.138.1\scripts\packages\VSCodeServer\src\eval.jl:181
 [20] withpath(f::VSCodeServer.var"#69#74"{Bool, Bool, Bool, Module, String, Int64, Int64, String, VSCodeServer.ReplRunCodeRequestParams}, path::String)
    @ VSCodeServer c:\Users\yolha\.vscode\extensions\julialang.language-julia-1.138.1\scripts\packages\VSCodeServer\src\repl.jl:276
 [21] (::VSCodeServer.var"#68#73"{Bool, Bool, Bool, Module, String, Int64, Int64, String, VSCodeServer.ReplRunCodeRequestParams})()
    @ VSCodeServer c:\Users\yolha\.vscode\extensions\julialang.language-julia-1.138.1\scripts\packages\VSCodeServer\src\eval.jl:179
 [22] hideprompt(f::VSCodeServer.var"#68#73"{Bool, Bool, Bool, Module, String, Int64, Int64, String, VSCodeServer.ReplRunCodeRequestParams})
    @ VSCodeServer c:\Users\yolha\.vscode\extensions\julialang.language-julia-1.138.1\scripts\packages\VSCodeServer\src\repl.jl:38
 [23] #67
    @ c:\Users\yolha\.vscode\extensions\julialang.language-julia-1.138.1\scripts\packages\VSCodeServer\src\eval.jl:150 [inlined]
 [24] with_logstate(f::VSCodeServer.var"#67#72"{Bool, Bool, Bool, Module, String, Int64, Int64, String, VSCodeServer.ReplRunCodeRequestParams}, logstate::Base.CoreLogging.LogState)
    @ Base.CoreLogging .\logging\logging.jl:522
 [25] with_logger
    @ .\logging\logging.jl:632 [inlined]
 [26] (::VSCodeServer.var"#66#71"{VSCodeServer.ReplRunCodeRequestParams})()
    @ VSCodeServer c:\Users\yolha\.vscode\extensions\julialang.language-julia-1.138.1\scripts\packages\VSCodeServer\src\eval.jl:263
 [27] #invokelatest#2
    @ .\essentials.jl:1055 [inlined]
 [28] invokelatest(::Any)
    @ Base .\essentials.jl:1052
 [29] (::VSCodeServer.var"#64#65")()
    @ VSCodeServer c:\Users\yolha\.vscode\extensions\julialang.language-julia-1.138.1\scripts\packages\VSCodeServer\src\eval.jl:34
in expression starting at c:\Users\yolha\Desktop\juju_tests\Tests\main.jl:85

same for ForwardDiff derivative

ForwardDiff.derivative(curve,1.0)

and FiniteDiff simply return [0,0,0].

[yolhan83]

@SouthEndMusic

[SouthEndMusic]

Adding in_place = false to the constructor should solve the ForwardDiff problem (I see I forgot to add that to the docstrings). The second derivate should be zero in the line segments.

[yolhan83]

Thank you it works with double ForwardDiff now.

[SouthEndMusic]

No problem, I added some clarifications in #411. I also made in_place = false the default, which is technically breaking but I think that doesn't cause to many problems in this case? :)

[yolhan83]

For your interest, this will work like a charm in my case

using DataInterpolations
using Plots
using ForwardDiff
using LinearAlgebra
using QuadGK

# generate some points (in real this could be the mass-centre of your artery)
N = 100 # Number of points
data = begin
    s = range(0,40,N)
    [[cos(0.2*i),sin(0.2*i),i] for i in s]
end

# Some parametrisation of the curve
x = range(0,1,N)
y = data
itp1 = CubicSpline(y,x)

# check the curve
Y = itp1(x)
plot(Tuple.(Y))
plot!(Tuple.(y))

# calculate the arc length
L,_ = quadgk(i->norm(ForwardDiff.derivative(itp1,i)),x)

# deduce the new parametrisation
nx = range(0,L,100)
itp = CubicSpline(Y,nx)
curve = SmoothArcLengthInterpolation(itp;m=100,in_place=false)

# check the new parametrisation
nnx = range(0,L,100)
Y = curve(nx)
plot(Tuple.(eachcol(Y)))
plot!(Tuple.(y))

[norm(ForwardDiff.derivative(curve,i)) for i in nx] ≈ ones(length(nx))