Skip to content

Smooth arclength interpolation second order derivatives #410

New issue
New issue
Closed
#411
Closed
Smooth arclength interpolation second order derivatives#410
#411
Labels
bugSomething isn't working
@yolhan83

Description

@yolhan83
yolhan83
opened on Mar 23, 2025 · edited by 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].

Activity

yolhan83
added
bugSomething isn't working
on Mar 23, 2025
yolhan83

yolhan83 commented on Mar 23, 2025

@yolhan83
yolhan83
on Mar 23, 2025
Author

@SouthEndMusic

SouthEndMusic

SouthEndMusic commented on Mar 23, 2025

@SouthEndMusic
SouthEndMusic
on Mar 23, 2025
Member

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

yolhan83 commented on Mar 23, 2025

@yolhan83
yolhan83
on Mar 23, 2025
Author

Thank you it works with double ForwardDiff now.

yolhan83
closed this as completedon Mar 23, 2025
SouthEndMusic
mentioned this on Mar 23, 2025
SouthEndMusic

SouthEndMusic commented on Mar 23, 2025

@SouthEndMusic
SouthEndMusic
on Mar 23, 2025
Member

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

yolhan83 commented on Mar 23, 2025

@yolhan83
yolhan83
on Mar 23, 2025
Author

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))
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

    bugSomething isn't working

    Type

    No type

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

      Development

      Participants

      @SouthEndMusic@yolhan83

      Issue actions
        Smooth arclength interpolation second order derivatives · Issue #410 · SciML/DataInterpolations.jl