Skip to content

Commit 955d75e

Browse files
mkittimkitti
authored
Merge branch 'master' into regularize
2 parents 37b0ad9 + 13da9c4 commit 955d75e

File tree

2 files changed

+43
-36
lines changed

2 files changed

+43
-36
lines changed

‎src/monotonic/monotonic.jl

Lines changed: 23 additions & 19 deletions
Original file line numberDiff line numberDiff line change
@@ -107,15 +107,15 @@ end
107107
Monotonic interpolation up to third order represented by type, knots and
108108
coefficients.
109109
"""
110-
struct MonotonicInterpolation{T, TCoeffs, TEl, TInterpolationType<:MonotonicInterpolationType,
110+
struct MonotonicInterpolation{T, TCoeffs1, TCoeffs2, TCoeffs3, TEl, TInterpolationType<:MonotonicInterpolationType,
111111
TKnots<:AbstractVector{<:Number}, TACoeff <: AbstractArray{TEl,1}} <: AbstractInterpolation{T,1,DimSpec{TInterpolationType}}
112112

113113
it::TInterpolationType
114114
knots::TKnots
115115
A::TACoeff # constant parts of piecewise polynomials
116-
m::Vector{TCoeffs} # coefficients of linear parts of piecewise polynomials
117-
c::Vector{TCoeffs} # coefficients of quadratic parts of piecewise polynomials
118-
d::Vector{TCoeffs} # coefficients of cubic parts of piecewise polynomials
116+
m::Vector{TCoeffs1} # coefficients of linear parts of piecewise polynomials
117+
c::Vector{TCoeffs2} # coefficients of quadratic parts of piecewise polynomials
118+
d::Vector{TCoeffs3} # coefficients of cubic parts of piecewise polynomials
119119
end
120120

121121
function Base.:(==)(o1::MonotonicInterpolation, o2::MonotonicInterpolation)
@@ -134,10 +134,10 @@ itpflag(A::MonotonicInterpolation) = A.it
134134
coefficients(A::MonotonicInterpolation) = A.A
135135

136136
function MonotonicInterpolation(::Type{TWeights}, it::TInterpolationType, knots::TKnots, A::AbstractArray{TEl,1},
137-
m::Vector{TCoeffs}, c::Vector{TCoeffs}, d::Vector{TCoeffs}) where {TWeights, TCoeffs, TEl, TInterpolationType<:MonotonicInterpolationType, TKnots<:AbstractVector{<:Number}}
137+
m::Vector{TCoeffs1}, c::Vector{TCoeffs2}, d::Vector{TCoeffs3}) where {TWeights, TCoeffs1, TCoeffs2, TCoeffs3, TEl, TInterpolationType<:MonotonicInterpolationType, TKnots<:AbstractVector{<:Number}}
138138

139139
isconcretetype(TInterpolationType) || error("The spline type must be a leaf type (was $TInterpolationType)")
140-
isconcretetype(TCoeffs) || warn("For performance reasons, consider using an array of a concrete type (eltype(A) == $(eltype(A)))")
140+
isconcretetype(tcoef(A)) || warn("For performance reasons, consider using an array of a concrete type (eltype(A) == $(eltype(A)))")
141141

142142
check_monotonic(knots, A)
143143

@@ -148,22 +148,23 @@ function MonotonicInterpolation(::Type{TWeights}, it::TInterpolationType, knots:
148148
T = typeof(cZero * first(A))
149149
end
150150

151-
MonotonicInterpolation{T, TCoeffs, TEl, TInterpolationType, TKnots, typeof(A)}(it, knots, A, m, c, d)
151+
MonotonicInterpolation{T, TCoeffs1, TCoeffs2, TCoeffs3, TEl, TInterpolationType, TKnots, typeof(A)}(it, knots, A, m, c, d)
152152
end
153153

154-
function interpolate(::Type{TWeights}, ::Type{TCoeffs}, knots::TKnots,
155-
A::AbstractArray{TEl,1}, it::TInterpolationType) where {TWeights,TCoeffs,TEl,TKnots<:AbstractVector{<:Number},TInterpolationType<:MonotonicInterpolationType}
154+
function interpolate(::Type{TWeights}, ::Type{TCoeffs1},::Type{TCoeffs2},::Type{TCoeffs3}, knots::TKnots,
155+
A::AbstractArray{TEl,1}, it::TInterpolationType) where {TWeights,TCoeffs1,TCoeffs2,TCoeffs3,TEl,TKnots<:AbstractVector{<:Number},TInterpolationType<:MonotonicInterpolationType}
156156

157157
check_monotonic(knots, A)
158158

159159
# first we need to determine tangents (m)
160160
n = length(knots)
161-
m, Δ = calcTangents(TCoeffs, knots, A, it)
162-
c = Vector{TCoeffs}(undef, n-1)
163-
d = Vector{TCoeffs}(undef, n-1)
161+
m, Δ = calcTangents(TCoeffs1, knots, A, it)
162+
c = Vector{TCoeffs2}(undef, n-1)
163+
d = Vector{TCoeffs3}(undef, n-1)
164164
for k ∈ 1:n-1
165165
if TInterpolationType == LinearMonotonicInterpolation
166-
c[k] = d[k] = zero(TCoeffs)
166+
c[k] = zero(TCoeffs2)
167+
d[k] = zero(TCoeffs3)
167168
else
168169
xdiff = knots[k+1] - knots[k]
169170
c[k] = (3*Δ[k] - 2*m[k] - m[k+1]) / xdiff
@@ -177,7 +178,9 @@ end
177178
function interpolate(knots::AbstractVector{<:Number}, A::AbstractArray{TEl,1},
178179
it::TInterpolationType) where {TEl,TInterpolationType<:MonotonicInterpolationType}
179180

180-
interpolate(tweight(A), tcoef(A), knots, A, it)
181+
interpolate(tweight(A),typeof(oneunit(eltype(A)) / oneunit(eltype(knots))),
182+
typeof(oneunit(eltype(A)) / oneunit(eltype(knots))^2),
183+
typeof(oneunit(eltype(A)) / oneunit(eltype(knots))^3),knots,A,it)
181184
end
182185

183186
function (itp::MonotonicInterpolation)(x::Number)
@@ -282,6 +285,7 @@ function calcTangents(::Type{TCoeffs}, x::AbstractVector{<:Number},
282285
n = length(x)
283286
Δ = Vector{TCoeffs}(undef, n-1)
284287
m = Vector{TCoeffs}(undef, n)
288+
buff = Vector{typeof(first(Δ)^2)}(undef, n)
285289
for k in 1:n-1
286290
Δk = (y[k+1] - y[k]) / (x[k+1] - x[k])
287291
Δ[k] = Δk
@@ -290,10 +294,10 @@ function calcTangents(::Type{TCoeffs}, x::AbstractVector{<:Number},
290294
else
291295
# If any consecutive secant lines change sign (i.e. curve changes direction), initialize the tangent to zero.
292296
# This is needed to make the interpolation monotonic. Otherwise set tangent to the average of the secants.
293-
if Δk <= zero(Δk)
297+
if Δ[k-1] * Δk <= zero(Δk^2)
294298
m[k] = zero(TCoeffs)
295299
else
296-
m[k] = Δ[k-1] * (Δ[k-1] + Δk) / 2.0
300+
m[k] = (Δ[k-1] + Δk) / 2.0
297301
end
298302
end
299303
end
@@ -317,8 +321,8 @@ function calcTangents(::Type{TCoeffs}, x::AbstractVector{<:Number},
317321
end
318322
α = m[k] / Δk
319323
β = m[k+1] / Δk
320-
τ = 3.0 / sqrt(α^2 + β^2)
321-
if τ < 1.0 # if we're outside the circle with radius 3 then move onto the circle
324+
τ = 3.0 * oneunit(α) / sqrt(α^2 + β^2)
325+
if τ < 1.0 # if we're outside the circle with radius 3 then move onto the circle
322326
m[k] = τ * α * Δk
323327
m[k+1] = τ * β * Δk
324328
end
@@ -341,7 +345,7 @@ function calcTangents(::Type{TCoeffs}, x::AbstractVector{<:Number},
341345
Δ[k] = Δk
342346
if k == 1 # left endpoint
343347
m[k] = Δk
344-
elseif Δ[k-1] * Δk <= zero(TCoeffs)
348+
elseif Δ[k-1] * Δk <= zero(Δk^2)
345349
m[k] = zero(TCoeffs)
346350
else
347351
α = (1.0 + (x[k+1] - x[k]) / (x[k+1] - x[k-1])) / 3.0

‎test/monotonic/runtests.jl

Lines changed: 20 additions & 17 deletions
Original file line numberDiff line numberDiff line change
@@ -1,11 +1,14 @@
11

22
using Test
33
using Interpolations
4-
4+
using Unitful
55
@testset "Monotonic" begin
6-
x = [0.0, 0.2, 0.5, 0.6, 0.9, 1.0]
6+
xa = [[0.0, 0.2, 0.5, 0.6, 0.9, 1.0],
7+
[0.0, 0.2, 0.5, 0.6, 0.9, 1.0]u"s"]
78
y = [(true, [-3.0, 0.0, 5.0, 10.0, 18.0, 22.0]),
8-
(false, [10.0, 0.0, -5.0, 10.0, -8.0, -2.0])]
9+
(false, [10.0, 0.0, -5.0, 10.0, -8.0, -2.0]),
10+
(false, [10.0, 0.0, -5.0, 10.0, -8.0, -2.0]u"m"),
11+
(true, [-3.0, 0.0, 5.0, 10.0, 18.0, 22.0]u"m")]
912

1013
# second item indicates if interpolating function can overshoot
1114
# for non-monotonic data
@@ -19,12 +22,12 @@ using Interpolations
1922
(SteffenMonotonicInterpolation(), false)]
2023

2124
for (it, overshoot) in itypes
22-
for yi in 1:length(y)
25+
for yi in 1:length(y), x in xa
2326
monotonic, ys = y[yi]
2427
itp = interpolate(x, ys, it)
2528
for j in 1:6
2629
# checking values at nodes
27-
@test itp(x[j]) ≈ ys[j] rtol = 1.e-15 atol = 1.e-14
30+
@test itp(x[j]) ≈ ys[j] rtol = 1.e-15 atol = 1.e-14*unit(first(ys))
2831

2932
# checking overshoot for non-monotonic data
3033
# and monotonicity for monotonic data
@@ -34,25 +37,25 @@ using Interpolations
3437
if j < 6
3538
r = range(x[j], stop = x[j+1], length = 100)
3639
for k in 1:length(r)-1
37-
@test (ys[j+1] - ys[j]) * (itp(r[k+1]) - itp(r[k])) >= 0.0
40+
@test (ys[j+1] - ys[j]) * (itp(r[k+1]) - itp(r[k])) >= zero(first(ys)^2)
3841
end
3942
end
4043
end
4144
extFlatItp = extrapolate(itp, Flat())
42-
@test extFlatItp(x[1]-1) ≈ itp(x[1])
43-
@test extFlatItp(x[end]+1) ≈ itp(x[end])
45+
@test extFlatItp(x[1]-oneunit(x[1])) ≈ itp(x[1])
46+
@test extFlatItp(x[end]+oneunit(x[1])) ≈ itp(x[end])
4447
extThrowItp = extrapolate(itp, Throw())
45-
@test_throws BoundsError extThrowItp(x[1]-1)
46-
@test_throws BoundsError extThrowItp(x[end]+1)
48+
@test_throws BoundsError extThrowItp(x[1]-oneunit(x[1]))
49+
@test_throws BoundsError extThrowItp(x[end]+oneunit(x[1]))
4750
extLineItp = extrapolate(itp, Line())
48-
@test extLineItp(x[1]-1) ≈ itp(x[1]) - Interpolations.gradient1(itp, x[1])
49-
@test extLineItp(x[end]+1) ≈ itp(x[end]) + Interpolations.gradient1(itp, x[end])
51+
@test extLineItp(x[1]-oneunit(x[1])) ≈ itp(x[1]) - Interpolations.gradient1(itp, x[1])*oneunit(x[1])
52+
@test extLineItp(x[end]+oneunit(x[1])) ≈ itp(x[end]) + Interpolations.gradient1(itp, x[end])*oneunit(x[1])
5053
extReflectItp = extrapolate(itp, Reflect())
51-
@test extReflectItp(x[1]-0.1) ≈ itp(x[1]+0.1)
52-
@test extReflectItp(x[end]+0.1) ≈ itp(x[end]-0.1)
54+
@test extReflectItp(x[1]-0.1*oneunit(x[1])) ≈ itp(x[1]+0.1*oneunit(x[1]))
55+
@test extReflectItp(x[end]+0.1*oneunit(x[1])) ≈ itp(x[end]-0.1*oneunit(x[1]))
5356
extPeriodicItp = extrapolate(itp, Periodic())
54-
@test extPeriodicItp(x[1]-0.1) ≈ itp(x[end]-0.1)
55-
@test extPeriodicItp(x[end]+0.1) ≈ itp(x[1]+0.1)
57+
@test extPeriodicItp(x[1]-0.1*oneunit(x[1])) ≈ itp(x[end]-0.1*oneunit(x[1]))
58+
@test extPeriodicItp(x[end]+0.1*oneunit(x[1])) ≈ itp(x[1]+0.1*oneunit(x[1]))
5659
end
5760
end
5861

@@ -61,7 +64,7 @@ using Interpolations
6164
y2 = [1.0, 2.0, 3.0, 4.0]
6265
y3 = [-3.0 0.0 5.0 10.0 18.0 22.0]
6366

64-
for (it, overshoot) in itypes
67+
for (it, overshoot) in itypes, x in xa
6568
@test_throws ErrorException interpolate(xWrong, y2, it)
6669
@test_throws DimensionMismatch interpolate(x, y2, it)
6770
@test_throws MethodError interpolate(x, y3, it)

0 commit comments

Comments
 (0)