From 9b0bbee9bdccf631acd74c5371ac5d9146fdfe9c Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Daniel=20Gonz=C3=A1lez?= <daniel.glez.arribas@gmail.com>
Date: Sat, 13 Jul 2024 11:57:06 +0200
Subject: [PATCH 01/11] Add accelerated quasilinear index lookup for
 LinearInterpolation

---
 src/interpolation_caches.jl  |  3 ++-
 src/interpolation_methods.jl |  9 ++++++++-
 src/interpolation_utils.jl   | 13 +++++++++++++
 3 files changed, 23 insertions(+), 2 deletions(-)

diff --git a/src/interpolation_caches.jl b/src/interpolation_caches.jl
index c7274471..049d867a 100644
--- a/src/interpolation_caches.jl
+++ b/src/interpolation_caches.jl
@@ -22,9 +22,10 @@ struct LinearInterpolation{uType, tType, IType, pType, T} <: AbstractInterpolati
     extrapolate::Bool
     idx_prev::Base.RefValue{Int}
     safetycopy::Bool
+    seems_linear::Bool
     function LinearInterpolation(u, t, I, p, extrapolate, safetycopy)
         new{typeof(u), typeof(t), typeof(I), typeof(p.slope), eltype(u)}(
-            u, t, I, p, extrapolate, Ref(1), safetycopy)
+            u, t, I, p, extrapolate, Ref(1), safetycopy, looks_quasilinear(t))
     end
 end
 
diff --git a/src/interpolation_methods.jl b/src/interpolation_methods.jl
index 5d09ceff..39bbcdb4 100644
--- a/src/interpolation_methods.jl
+++ b/src/interpolation_methods.jl
@@ -1,7 +1,14 @@
 function _interpolate(A, t)
     ((t < A.t[1] || t > A.t[end]) && !A.extrapolate) &&
         throw(ExtrapolationError())
-    val, idx_prev = _interpolate(A, t, A.idx_prev[])
+    idx_guess = if hasfield(typeof(A), :seems_linear) && A.seems_linear
+        f = (t - first(A.t)) / (last(A.t) - first(A.t))
+        i_0, i_f = firstindex(A.t), lastindex(A.t)
+        round(typeof(firstindex(A.t)), f * (i_f - i_0) + i_0)
+    else
+        A.idx_prev[]
+    end
+    val, idx_prev = _interpolate(A, t, idx_guess)
     A.idx_prev[] = idx_prev
     return val
 end
diff --git a/src/interpolation_utils.jl b/src/interpolation_utils.jl
index 466248b1..881d913c 100644
--- a/src/interpolation_utils.jl
+++ b/src/interpolation_utils.jl
@@ -106,6 +106,19 @@ function munge_data(U::StridedMatrix, t::AbstractVector, safetycopy::Bool)
     return readonly_wrap(U), readonly_wrap(t)
 end
 
+"""
+Determine if the abscissae are distributed in almost linear fashion, using as a cheap
+criterion to determine whether the middle element is close to the average of the first and
+last elements of the abscissae array, with a relative threshold of 10% by default
+"""
+function looks_quasilinear(t; threshold = 0.1)
+    length(t) <= 2 && return true
+    i_0, i_f = firstindex(t), lastindex(t)
+    t_0, t_f = first(t), last(t)
+    middle_index = round(typeof(i_0), 0.5 * (i_0 + i_f))
+    abs(0.5 * (t_0 + t_f) - t[middle_index]) < threshold * abs(t_0 - t_f)
+end
+
 # Don't nest ReadOnlyArrays
 readonly_wrap(a::AbstractArray) = ReadOnlyArray(a)
 readonly_wrap(a::ReadOnlyArray) = a

From 12ebad8ca72efc68c520f83d047cc22157a697b0 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Daniel=20Gonz=C3=A1lez?= <daniel.glez.arribas@gmail.com>
Date: Sat, 13 Jul 2024 12:49:50 +0200
Subject: [PATCH 02/11] Add user-facing mechanism to disable/enable linear-like
 guessing

---
 src/interpolation_caches.jl | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/src/interpolation_caches.jl b/src/interpolation_caches.jl
index 049d867a..c12f1719 100644
--- a/src/interpolation_caches.jl
+++ b/src/interpolation_caches.jl
@@ -23,7 +23,9 @@ struct LinearInterpolation{uType, tType, IType, pType, T} <: AbstractInterpolati
     idx_prev::Base.RefValue{Int}
     safetycopy::Bool
     seems_linear::Bool
-    function LinearInterpolation(u, t, I, p, extrapolate, safetycopy)
+    function LinearInterpolation(
+            u, t, I, p, extrapolate, safetycopy; guess_linear_t = nothing)
+        guess_linear_t = isnothing(guess_linear_t) ? looks_quasilinear(t) : guess_linear_t
         new{typeof(u), typeof(t), typeof(I), typeof(p.slope), eltype(u)}(
             u, t, I, p, extrapolate, Ref(1), safetycopy, looks_quasilinear(t))
     end

From 985e4545faa91658dd33028b9893f0de74ad3757 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Daniel=20Gonz=C3=A1lez?= <daniel.glez.arribas@gmail.com>
Date: Sun, 14 Jul 2024 18:26:16 +0200
Subject: [PATCH 03/11] Tighten quasi-regularity criterion

---
 src/interpolation_caches.jl |  7 ++++---
 src/interpolation_utils.jl  | 21 ++++++++++++---------
 2 files changed, 16 insertions(+), 12 deletions(-)

diff --git a/src/interpolation_caches.jl b/src/interpolation_caches.jl
index c12f1719..9694cd9a 100644
--- a/src/interpolation_caches.jl
+++ b/src/interpolation_caches.jl
@@ -24,10 +24,11 @@ struct LinearInterpolation{uType, tType, IType, pType, T} <: AbstractInterpolati
     safetycopy::Bool
     seems_linear::Bool
     function LinearInterpolation(
-            u, t, I, p, extrapolate, safetycopy; guess_linear_t = nothing)
-        guess_linear_t = isnothing(guess_linear_t) ? looks_quasilinear(t) : guess_linear_t
+            u, t, I, p, extrapolate, safetycopy; quasiregular = 1e-2)
+        quasiregular = quasiregular isa Bool ? quasiregular :
+                       looks_quasiregular(t; threshold = quasiregular)
         new{typeof(u), typeof(t), typeof(I), typeof(p.slope), eltype(u)}(
-            u, t, I, p, extrapolate, Ref(1), safetycopy, looks_quasilinear(t))
+            u, t, I, p, extrapolate, Ref(1), safetycopy, quasiregular)
     end
 end
 
diff --git a/src/interpolation_utils.jl b/src/interpolation_utils.jl
index 881d913c..fa46c214 100644
--- a/src/interpolation_utils.jl
+++ b/src/interpolation_utils.jl
@@ -106,17 +106,20 @@ function munge_data(U::StridedMatrix, t::AbstractVector, safetycopy::Bool)
     return readonly_wrap(U), readonly_wrap(t)
 end
 
-"""
-Determine if the abscissae are distributed in almost linear fashion, using as a cheap
-criterion to determine whether the middle element is close to the average of the first and
-last elements of the abscissae array, with a relative threshold of 10% by default
-"""
-function looks_quasilinear(t; threshold = 0.1)
+# """
+# Determine if the abscissae are regularly distributed, taking the standard \
+# deviation of the difference between the array of abscissae with respect to \
+# the straight line linking its first and last elements
+# """
+function looks_quasiregular(t; threshold = 1e-2)
     length(t) <= 2 && return true
-    i_0, i_f = firstindex(t), lastindex(t)
     t_0, t_f = first(t), last(t)
-    middle_index = round(typeof(i_0), 0.5 * (i_0 + i_f))
-    abs(0.5 * (t_0 + t_f) - t[middle_index]) < threshold * abs(t_0 - t_f)
+    t_span = t_f - t_0
+    tspan_over_N = t_span * length(t)^(-1)
+    norm_var = sum(
+        (t_i - t_0 - i * tspan_over_N)^2 for (i, t_i) in enumerate(t)
+    ) / (length(t) * t_span^2)
+    norm_var < threshold^2
 end
 
 # Don't nest ReadOnlyArrays

From 2ab91ed7e830f67b44e4d002a0b1277b2eba293a Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Daniel=20Gonz=C3=A1lez?= <daniel.glez.arribas@gmail.com>
Date: Mon, 29 Jul 2024 20:31:33 +0200
Subject: [PATCH 04/11] [WIP] Wording, toggling settings for linear lookup

---
 src/interpolation_caches.jl  | 32 +++++++++++++++++++++++---------
 src/interpolation_methods.jl |  3 ++-
 src/interpolation_utils.jl   |  5 ++++-
 3 files changed, 29 insertions(+), 11 deletions(-)

diff --git a/src/interpolation_caches.jl b/src/interpolation_caches.jl
index 9694cd9a..1f57808c 100644
--- a/src/interpolation_caches.jl
+++ b/src/interpolation_caches.jl
@@ -13,6 +13,10 @@ Extrapolation extends the last linear polynomial on each side.
 
   - `extrapolate`: boolean value to allow extrapolation. Defaults to `false`.
   - `safetycopy`: boolean value to make a copy of `u` and `t`. Defaults to `true`.
+  - `assume_linear_t`: boolean value to specify a faster index lookup behaviour for
+    evenly-distributed abscissae. Alternatively, a numerical threshold may be specified
+    for a test based on the normalized standard deviation of the difference with respect
+    to the straight line. Defaults to 1e-2
 """
 struct LinearInterpolation{uType, tType, IType, pType, T} <: AbstractInterpolation{T}
     u::uType
@@ -22,22 +26,22 @@ struct LinearInterpolation{uType, tType, IType, pType, T} <: AbstractInterpolati
     extrapolate::Bool
     idx_prev::Base.RefValue{Int}
     safetycopy::Bool
-    seems_linear::Bool
+    use_linear_lookup::Bool
     function LinearInterpolation(
-            u, t, I, p, extrapolate, safetycopy; quasiregular = 1e-2)
-        quasiregular = quasiregular isa Bool ? quasiregular :
-                       looks_quasiregular(t; threshold = quasiregular)
+            u, t, I, p, extrapolate, safetycopy, assume_linear_t)
+        linear_flag = seems_linear(assume_linear_t, t)
         new{typeof(u), typeof(t), typeof(I), typeof(p.slope), eltype(u)}(
-            u, t, I, p, extrapolate, Ref(1), safetycopy, quasiregular)
+            u, t, I, p, extrapolate, Ref(1), safetycopy, linear_flag)
     end
 end
 
-function LinearInterpolation(u, t; extrapolate = false, safetycopy = true)
+function LinearInterpolation(
+        u, t; extrapolate = false, safetycopy = true, assume_linear_t = 1e-2)
     u, t = munge_data(u, t, safetycopy)
     p = LinearParameterCache(u, t)
     A = LinearInterpolation(u, t, nothing, p, extrapolate, safetycopy)
     I = cumulative_integral(A)
-    LinearInterpolation(u, t, I, p, extrapolate, safetycopy)
+    LinearInterpolation(u, t, I, p, extrapolate, safetycopy, assume_linear_t)
 end
 
 """
@@ -66,11 +70,14 @@ struct QuadraticInterpolation{uType, tType, IType, pType, T} <: AbstractInterpol
     extrapolate::Bool
     idx_prev::Base.RefValue{Int}
     safetycopy::Bool
-    function QuadraticInterpolation(u, t, I, p, mode, extrapolate, safetycopy)
+    use_linear_lookup::Bool
+    function QuadraticInterpolation(
+            u, t, I, p, mode, extrapolate, safetycopy, assume_linear_t)
         mode ∈ (:Forward, :Backward) ||
             error("mode should be :Forward or :Backward for QuadraticInterpolation")
+        linear_flag = seems_linear(assume_linear_t, t)
         new{typeof(u), typeof(t), typeof(I), typeof(p.l₀), eltype(u)}(
-            u, t, I, p, mode, extrapolate, Ref(1), safetycopy)
+            u, t, I, p, mode, extrapolate, Ref(1), safetycopy, linear_flag)
     end
 end
 
@@ -233,6 +240,7 @@ struct ConstantInterpolation{uType, tType, IType, T} <: AbstractInterpolation{T}
     extrapolate::Bool
     idx_prev::Base.RefValue{Int}
     safetycopy::Bool
+    use_linear_lookup::Bool
     function ConstantInterpolation(u, t, I, dir, extrapolate, safetycopy)
         new{typeof(u), typeof(t), typeof(I), eltype(u)}(
             u, t, I, nothing, dir, extrapolate, Ref(1), safetycopy)
@@ -274,6 +282,7 @@ struct QuadraticSpline{uType, tType, IType, pType, tAType, dType, zType, T} <:
     extrapolate::Bool
     idx_prev::Base.RefValue{Int}
     safetycopy::Bool
+    use_linear_lookup::Bool
     function QuadraticSpline(u, t, I, p, tA, d, z, extrapolate, safetycopy)
         new{typeof(u), typeof(t), typeof(I), typeof(p.σ), typeof(tA),
             typeof(d), typeof(z), eltype(u)}(u,
@@ -358,6 +367,7 @@ struct CubicSpline{uType, tType, IType, pType, hType, zType, T} <: AbstractInter
     extrapolate::Bool
     idx_prev::Base.RefValue{Int}
     safetycopy::Bool
+    use_linear_lookup::Bool
     function CubicSpline(u, t, I, p, h, z, extrapolate, safetycopy)
         new{typeof(u), typeof(t), typeof(I), typeof(p.c₁), typeof(h), typeof(z), eltype(u)}(
             u,
@@ -454,6 +464,7 @@ struct BSplineInterpolation{uType, tType, pType, kType, cType, NType, T} <:
     extrapolate::Bool
     idx_prev::Base.RefValue{Int}
     safetycopy::Bool
+    use_linear_lookup::Bool
     function BSplineInterpolation(u,
             t,
             d,
@@ -586,6 +597,7 @@ struct BSplineApprox{uType, tType, pType, kType, cType, NType, T} <:
     extrapolate::Bool
     idx_prev::Base.RefValue{Int}
     safetycopy::Bool
+    use_linear_lookup::Bool
     function BSplineApprox(u,
             t,
             d,
@@ -730,6 +742,7 @@ struct CubicHermiteSpline{uType, tType, IType, duType, pType, T} <: AbstractInte
     extrapolate::Bool
     idx_prev::Base.RefValue{Int}
     safetycopy::Bool
+    use_linear_lookup::Bool
     function CubicHermiteSpline(du, u, t, I, p, extrapolate, safetycopy)
         new{typeof(u), typeof(t), typeof(I), typeof(du), typeof(p.c₁), eltype(u)}(
             du, u, t, I, p, extrapolate, Ref(1), safetycopy)
@@ -773,6 +786,7 @@ struct QuinticHermiteSpline{uType, tType, IType, duType, dduType, pType, T} <:
     extrapolate::Bool
     idx_prev::Base.RefValue{Int}
     safetycopy::Bool
+    use_linear_lookup::Bool
     function QuinticHermiteSpline(ddu, du, u, t, I, p, extrapolate, safetycopy)
         new{typeof(u), typeof(t), typeof(I), typeof(du),
             typeof(ddu), typeof(p.c₁), eltype(u)}(
diff --git a/src/interpolation_methods.jl b/src/interpolation_methods.jl
index 39bbcdb4..d8870568 100644
--- a/src/interpolation_methods.jl
+++ b/src/interpolation_methods.jl
@@ -1,7 +1,8 @@
 function _interpolate(A, t)
     ((t < A.t[1] || t > A.t[end]) && !A.extrapolate) &&
         throw(ExtrapolationError())
-    idx_guess = if hasfield(typeof(A), :seems_linear) && A.seems_linear
+    idx_guess = if hasfield(typeof(A), :use_linear_lookup) &&
+                   A.use_linear_lookup
         f = (t - first(A.t)) / (last(A.t) - first(A.t))
         i_0, i_f = firstindex(A.t), lastindex(A.t)
         round(typeof(firstindex(A.t)), f * (i_f - i_0) + i_0)
diff --git a/src/interpolation_utils.jl b/src/interpolation_utils.jl
index fa46c214..9e9b47a5 100644
--- a/src/interpolation_utils.jl
+++ b/src/interpolation_utils.jl
@@ -106,12 +106,15 @@ function munge_data(U::StridedMatrix, t::AbstractVector, safetycopy::Bool)
     return readonly_wrap(U), readonly_wrap(t)
 end
 
+seems_linear(assume_linear_t::Bool, _) = assume_linear_t
+seems_linear(assume_qr::Number, t) = looks_linear(t; threshold = assume_qr)
+
 # """
 # Determine if the abscissae are regularly distributed, taking the standard \
 # deviation of the difference between the array of abscissae with respect to \
 # the straight line linking its first and last elements
 # """
-function looks_quasiregular(t; threshold = 1e-2)
+function looks_linear(t; threshold = 1e-2)
     length(t) <= 2 && return true
     t_0, t_f = first(t), last(t)
     t_span = t_f - t_0

From c1e7d3f9769d84c18b2050220dd797eb87a69075 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Daniel=20Gonz=C3=A1lez?= <daniel.glez.arribas@gmail.com>
Date: Mon, 29 Jul 2024 22:28:37 +0200
Subject: [PATCH 05/11] Implement review suggestions + use in other caches

---
 ext/DataInterpolationsChainRulesCoreExt.jl |   2 +-
 src/derivatives.jl                         |  16 +-
 src/integral_inverses.jl                   |   8 +-
 src/integrals.jl                           |   4 +-
 src/interpolation_caches.jl                | 208 +++++++++++++++------
 src/interpolation_methods.jl               |  41 ++--
 src/interpolation_utils.jl                 |  12 +-
 7 files changed, 191 insertions(+), 100 deletions(-)

diff --git a/ext/DataInterpolationsChainRulesCoreExt.jl b/ext/DataInterpolationsChainRulesCoreExt.jl
index 9c33b09c..312e67e3 100644
--- a/ext/DataInterpolationsChainRulesCoreExt.jl
+++ b/ext/DataInterpolationsChainRulesCoreExt.jl
@@ -52,7 +52,7 @@ end
 
 function u_tangent(A::LinearInterpolation, t, Δ)
     out = zero(A.u)
-    idx = get_idx(A.t, t, A.idx_prev[])
+    idx = get_idx(A, t, A.idx_prev[])
     t_factor = (t - A.t[idx]) / (A.t[idx + 1] - A.t[idx])
     out[idx] = Δ * (one(eltype(out)) - t_factor)
     out[idx + 1] = Δ * t_factor
diff --git a/src/derivatives.jl b/src/derivatives.jl
index 75872095..e7cb97fd 100644
--- a/src/derivatives.jl
+++ b/src/derivatives.jl
@@ -18,7 +18,7 @@ function derivative(A, t, order = 1)
 end
 
 function _derivative(A::LinearInterpolation, t::Number, iguess)
-    idx = get_idx(A.t, t, iguess; idx_shift = -1, ub_shift = -1, side = :first)
+    idx = get_idx(A, t, iguess; idx_shift = -1, ub_shift = -1, side = :first)
     slope = get_parameters(A, idx)
     slope, idx
 end
@@ -108,7 +108,7 @@ function _derivative(A::LagrangeInterpolation{<:AbstractMatrix}, t::Number, idx)
 end
 
 function _derivative(A::AkimaInterpolation{<:AbstractVector}, t::Number, iguess)
-    idx = get_idx(A.t, t, iguess; idx_shift = -1, side = :first)
+    idx = get_idx(A, t, iguess; idx_shift = -1, side = :first)
     j = min(idx, length(A.c))  # for smooth derivative at A.t[end]
     wj = t - A.t[idx]
     (@evalpoly wj A.b[idx] 2A.c[j] 3A.d[j]), idx
@@ -130,14 +130,14 @@ end
 
 # QuadraticSpline Interpolation
 function _derivative(A::QuadraticSpline{<:AbstractVector}, t::Number, iguess)
-    idx = get_idx(A.t, t, iguess; lb = 2, ub_shift = 0, side = :first)
+    idx = get_idx(A, t, iguess; lb = 2, ub_shift = 0, side = :first)
     σ = get_parameters(A, idx - 1)
     A.z[idx - 1] + 2σ * (t - A.t[idx - 1]), idx
 end
 
 # CubicSpline Interpolation
 function _derivative(A::CubicSpline{<:AbstractVector}, t::Number, iguess)
-    idx = get_idx(A.t, t, iguess)
+    idx = get_idx(A, t, iguess)
     Δt₁ = t - A.t[idx]
     Δt₂ = A.t[idx + 1] - t
     dI = (-A.z[idx] * Δt₂^2 + A.z[idx + 1] * Δt₁^2) / (2A.h[idx + 1])
@@ -151,7 +151,7 @@ function _derivative(A::BSplineInterpolation{<:AbstractVector{<:Number}}, t::Num
     # change t into param [0 1]
     t < A.t[1] && return zero(A.u[1]), 1
     t > A.t[end] && return zero(A.u[end]), lastindex(t)
-    idx = get_idx(A.t, t, iguess)
+    idx = get_idx(A, t, iguess)
     n = length(A.t)
     scale = (A.p[idx + 1] - A.p[idx]) / (A.t[idx + 1] - A.t[idx])
     t_ = A.p[idx] + (t - A.t[idx]) * scale
@@ -173,7 +173,7 @@ function _derivative(A::BSplineApprox{<:AbstractVector{<:Number}}, t::Number, ig
     # change t into param [0 1]
     t < A.t[1] && return zero(A.u[1]), 1
     t > A.t[end] && return zero(A.u[end]), lastindex(t)
-    idx = get_idx(A.t, t, iguess)
+    idx = get_idx(A, t, iguess)
     scale = (A.p[idx + 1] - A.p[idx]) / (A.t[idx + 1] - A.t[idx])
     t_ = A.p[idx] + (t - A.t[idx]) * scale
     N = t isa ForwardDiff.Dual ? zeros(eltype(t), A.h) : A.N
@@ -192,7 +192,7 @@ end
 # Cubic Hermite Spline
 function _derivative(
         A::CubicHermiteSpline{<:AbstractVector{<:Number}}, t::Number, iguess)
-    idx = get_idx(A.t, t, iguess)
+    idx = get_idx(A, t, iguess)
     Δt₀ = t - A.t[idx]
     Δt₁ = t - A.t[idx + 1]
     out = A.du[idx]
@@ -204,7 +204,7 @@ end
 # Quintic Hermite Spline
 function _derivative(
         A::QuinticHermiteSpline{<:AbstractVector{<:Number}}, t::Number, iguess)
-    idx = get_idx(A.t, t, iguess)
+    idx = get_idx(A, t, iguess)
     Δt₀ = t - A.t[idx]
     Δt₁ = t - A.t[idx + 1]
     out = A.du[idx] + A.ddu[idx] * Δt₀
diff --git a/src/integral_inverses.jl b/src/integral_inverses.jl
index 31d0853c..aebee9d5 100644
--- a/src/integral_inverses.jl
+++ b/src/integral_inverses.jl
@@ -59,7 +59,7 @@ end
 
 function _interpolate(
         A::LinearInterpolationIntInv{<:AbstractVector{<:Number}}, t::Number, iguess)
-    idx = get_idx(A.t, t, iguess)
+    idx = get_idx(A, t, iguess)
     Δt = t - A.t[idx]
     x = A.itp.u[idx]
     slope = get_parameters(A.itp, idx)
@@ -104,13 +104,13 @@ end
 
 function _interpolate(
         A::ConstantInterpolationIntInv{<:AbstractVector{<:Number}}, t::Number, iguess)
-    idx = get_idx(A.t, t, iguess; ub_shift = 0)
+    idx = get_idx(A, t, iguess; ub_shift = 0)
     if A.itp.dir === :left
         # :left means that value to the left is used for interpolation
-        idx_ = get_idx(A.t, t, idx; lb = 1, ub_shift = 0)
+        idx_ = get_idx(A, t, idx; lb = 1, ub_shift = 0)
     else
         # :right means that value to the right is used for interpolation
-        idx_ = get_idx(A.t, t, idx; side = :first, lb = 1, ub_shift = 0)
+        idx_ = get_idx(A, t, idx; side = :first, lb = 1, ub_shift = 0)
     end
     A.u[idx] + (t - A.t[idx]) / A.itp.u[idx_], idx
 end
diff --git a/src/integrals.jl b/src/integrals.jl
index 03ea26c5..3ae893f7 100644
--- a/src/integrals.jl
+++ b/src/integrals.jl
@@ -8,9 +8,9 @@ function integral(A::AbstractInterpolation, t1::Number, t2::Number)
     ((t2 < A.t[1] || t2 > A.t[end]) && !A.extrapolate) && throw(ExtrapolationError())
     !hasfield(typeof(A), :I) && throw(IntegralNotFoundError())
     # the index less than or equal to t1
-    idx1 = get_idx(A.t, t1, 0)
+    idx1 = get_idx(A, t1, 0)
     # the index less than t2
-    idx2 = get_idx(A.t, t2, 0; idx_shift = -1, side = :first)
+    idx2 = get_idx(A, t2, 0; idx_shift = -1, side = :first)
 
     if A.cache_parameters
         total = A.I[idx2] - A.I[idx1]
diff --git a/src/interpolation_caches.jl b/src/interpolation_caches.jl
index aa452589..c4017460 100644
--- a/src/interpolation_caches.jl
+++ b/src/interpolation_caches.jl
@@ -29,9 +29,9 @@ struct LinearInterpolation{uType, tType, IType, pType, T} <: AbstractInterpolati
     cache_parameters::Bool
     use_linear_lookup::Bool
     function LinearInterpolation(u, t, I, p, extrapolate, cache_parameters, assume_linear_t)
-        linear_flag = seems_linear(assume_linear_t, t)
+        linear_lookup = seems_linear(assume_linear_t, t)
         new{typeof(u), typeof(t), typeof(I), typeof(p.slope), eltype(u)}(
-            u, t, I, p, extrapolate, Ref(1), cache_parameters, linear_flag)
+            u, t, I, p, extrapolate, Ref(1), cache_parameters, linear_lookup)
     end
 end
 
@@ -39,7 +39,8 @@ function LinearInterpolation(
         u, t; extrapolate = false, cache_parameters = false, assume_linear_t = 1e-2)
     u, t = munge_data(u, t)
     p = LinearParameterCache(u, t, cache_parameters)
-    A = LinearInterpolation(u, t, nothing, p, extrapolate, cache_parameters)
+    A = LinearInterpolation(
+        u, t, nothing, p, extrapolate, cache_parameters, assume_linear_t)
     I = cumulative_integral(A, cache_parameters)
     LinearInterpolation(u, t, I, p, extrapolate, cache_parameters, assume_linear_t)
 end
@@ -60,6 +61,10 @@ Extrapolation extends the last quadratic polynomial on each side.
 
   - `extrapolate`: boolean value to allow extrapolation. Defaults to `false`.
   - `cache_parameters`: precompute parameters at initialization for faster interpolation computations. Note: if activated, `u` and `t` should not be modified. Defaults to `false`.
+  - `assume_linear_t`: boolean value to specify a faster index lookup behaviour for
+    evenly-distributed abscissae. Alternatively, a numerical threshold may be specified
+    for a test based on the normalized standard deviation of the difference with respect
+    to the straight line. Defaults to 1e-2
 """
 struct QuadraticInterpolation{uType, tType, IType, pType, T} <: AbstractInterpolation{T}
     u::uType
@@ -75,22 +80,25 @@ struct QuadraticInterpolation{uType, tType, IType, pType, T} <: AbstractInterpol
             u, t, I, p, mode, extrapolate, cache_parameters, assume_linear_t)
         mode ∈ (:Forward, :Backward) ||
             error("mode should be :Forward or :Backward for QuadraticInterpolation")
-        linear_flag = seems_linear(assume_linear_t, t)
+        linear_lookup = seems_linear(assume_linear_t, t)
         new{typeof(u), typeof(t), typeof(I), typeof(p.l₀), eltype(u)}(
-            u, t, I, p, mode, extrapolate, Ref(1), cache_parameters, linear_flag)
+            u, t, I, p, mode, extrapolate, Ref(1), cache_parameters, linear_lookup)
     end
 end
 
-function QuadraticInterpolation(u, t, mode; extrapolate = false, cache_parameters = false)
+function QuadraticInterpolation(
+        u, t, mode; extrapolate = false, cache_parameters = false, assume_linear_t = 1e-2)
     u, t = munge_data(u, t)
+    linear_lookup = seems_linear(assume_linear_t, t)
     p = QuadraticParameterCache(u, t, cache_parameters)
-    A = QuadraticInterpolation(u, t, nothing, p, mode, extrapolate, cache_parameters)
+    A = QuadraticInterpolation(
+        u, t, nothing, p, mode, extrapolate, cache_parameters, linear_lookup)
     I = cumulative_integral(A, cache_parameters)
-    QuadraticInterpolation(u, t, I, p, mode, extrapolate, cache_parameters)
+    QuadraticInterpolation(u, t, I, p, mode, extrapolate, cache_parameters, linear_lookup)
 end
 
-function QuadraticInterpolation(u, t; extrapolate = false, cache_parameters = false)
-    QuadraticInterpolation(u, t, :Forward; extrapolate, cache_parameters)
+function QuadraticInterpolation(u, t; kwargs...)
+    QuadraticInterpolation(u, t, :Forward; kwargs...)
 end
 
 """
@@ -156,6 +164,10 @@ Extrapolation extends the last cubic polynomial on each side.
 
   - `extrapolate`: boolean value to allow extrapolation. Defaults to `false`.
   - `cache_parameters`: precompute parameters at initialization for faster interpolation computations. Note: if activated, `u` and `t` should not be modified. Defaults to `false`.
+  - `assume_linear_t`: boolean value to specify a faster index lookup behaviour for
+    evenly-distributed abscissae. Alternatively, a numerical threshold may be specified
+    for a test based on the normalized standard deviation of the difference with respect
+    to the straight line. Defaults to 1e-2
 """
 struct AkimaInterpolation{uType, tType, IType, bType, cType, dType, T} <:
        AbstractInterpolation{T}
@@ -168,7 +180,10 @@ struct AkimaInterpolation{uType, tType, IType, bType, cType, dType, T} <:
     extrapolate::Bool
     idx_prev::Base.RefValue{Int}
     cache_parameters::Bool
-    function AkimaInterpolation(u, t, I, b, c, d, extrapolate, cache_parameters)
+    use_linear_lookup::Bool
+    function AkimaInterpolation(
+            u, t, I, b, c, d, extrapolate, cache_parameters, assume_linear_t)
+        linear_lookup = seems_linear(assume_linear_t, t)
         new{typeof(u), typeof(t), typeof(I), typeof(b), typeof(c),
             typeof(d), eltype(u)}(u,
             t,
@@ -178,13 +193,16 @@ struct AkimaInterpolation{uType, tType, IType, bType, cType, dType, T} <:
             d,
             extrapolate,
             Ref(1),
-            cache_parameters
+            cache_parameters,
+            linear_lookup
         )
     end
 end
 
-function AkimaInterpolation(u, t; extrapolate = false, cache_parameters = false)
+function AkimaInterpolation(
+        u, t; extrapolate = false, cache_parameters = false, assume_linear_t = 1e-2)
     u, t = munge_data(u, t)
+    linear_lookup = seems_linear(assume_linear_t, t)
     n = length(t)
     dt = diff(t)
     m = Array{eltype(u)}(undef, n + 3)
@@ -205,9 +223,10 @@ function AkimaInterpolation(u, t; extrapolate = false, cache_parameters = false)
     c = (3.0 .* m[3:(end - 2)] .- 2.0 .* b[1:(end - 1)] .- b[2:end]) ./ dt
     d = (b[1:(end - 1)] .+ b[2:end] .- 2.0 .* m[3:(end - 2)]) ./ dt .^ 2
 
-    A = AkimaInterpolation(u, t, nothing, b, c, d, extrapolate, cache_parameters)
+    A = AkimaInterpolation(
+        u, t, nothing, b, c, d, extrapolate, cache_parameters, linear_lookup)
     I = cumulative_integral(A, cache_parameters)
-    AkimaInterpolation(u, t, I, b, c, d, extrapolate, cache_parameters)
+    AkimaInterpolation(u, t, I, b, c, d, extrapolate, cache_parameters, linear_lookup)
 end
 
 """
@@ -227,6 +246,10 @@ Extrapolation extends the last constant polynomial at the end points on each sid
   - `dir`: indicates which value should be used for interpolation (`:left` or `:right`).
   - `extrapolate`: boolean value to allow extrapolation. Defaults to `false`.
   - `cache_parameters`: precompute parameters at initialization for faster interpolation computations. Note: if activated, `u` and `t` should not be modified. Defaults to `false`.
+  - `assume_linear_t`: boolean value to specify a faster index lookup behaviour for
+    evenly-distributed abscissae. Alternatively, a numerical threshold may be specified
+    for a test based on the normalized standard deviation of the difference with respect
+    to the straight line. Defaults to 1e-2
 """
 struct ConstantInterpolation{uType, tType, IType, T} <: AbstractInterpolation{T}
     u::uType
@@ -238,18 +261,22 @@ struct ConstantInterpolation{uType, tType, IType, T} <: AbstractInterpolation{T}
     idx_prev::Base.RefValue{Int}
     cache_parameters::Bool
     use_linear_lookup::Bool
-    function ConstantInterpolation(u, t, I, dir, extrapolate, cache_parameters)
+    function ConstantInterpolation(
+            u, t, I, dir, extrapolate, cache_parameters, assume_linear_t)
+        linear_lookup = seems_linear(assume_linear_t, t)
         new{typeof(u), typeof(t), typeof(I), eltype(u)}(
-            u, t, I, nothing, dir, extrapolate, Ref(1), cache_parameters)
+            u, t, I, nothing, dir, extrapolate, Ref(1), cache_parameters, linear_lookup)
     end
 end
 
 function ConstantInterpolation(
-        u, t; dir = :left, extrapolate = false, cache_parameters = false)
+        u, t; dir = :left, extrapolate = false,
+        cache_parameters = false, assume_linear_t = 1e-2)
     u, t = munge_data(u, t)
-    A = ConstantInterpolation(u, t, nothing, dir, extrapolate, cache_parameters)
+    A = ConstantInterpolation(
+        u, t, nothing, dir, extrapolate, cache_parameters, assume_linear_t)
     I = cumulative_integral(A, cache_parameters)
-    ConstantInterpolation(u, t, I, dir, extrapolate, cache_parameters)
+    ConstantInterpolation(u, t, I, dir, extrapolate, cache_parameters, assume_linear_t)
 end
 
 """
@@ -267,6 +294,10 @@ Extrapolation extends the last quadratic polynomial on each side.
 
   - `extrapolate`: boolean value to allow extrapolation. Defaults to `false`.
   - `cache_parameters`: precompute parameters at initialization for faster interpolation computations. Note: if activated, `u` and `t` should not be modified. Defaults to `false`.
+  - `assume_linear_t`: boolean value to specify a faster index lookup behaviour for
+    evenly-distributed abscissae. Alternatively, a numerical threshold may be specified
+    for a test based on the normalized standard deviation of the difference with respect
+    to the straight line. Defaults to 1e-2
 """
 struct QuadraticSpline{uType, tType, IType, pType, tAType, dType, zType, T} <:
        AbstractInterpolation{T}
@@ -281,7 +312,9 @@ struct QuadraticSpline{uType, tType, IType, pType, tAType, dType, zType, T} <:
     idx_prev::Base.RefValue{Int}
     cache_parameters::Bool
     use_linear_lookup::Bool
-    function QuadraticSpline(u, t, I, p, tA, d, z, extrapolate, cache_parameters)
+    function QuadraticSpline(
+            u, t, I, p, tA, d, z, extrapolate, cache_parameters, assume_linear_t)
+        linear_lookup = seems_linear(assume_linear_t, t)
         new{typeof(u), typeof(t), typeof(I), typeof(p.σ), typeof(tA),
             typeof(d), typeof(z), eltype(u)}(u,
             t,
@@ -292,15 +325,18 @@ struct QuadraticSpline{uType, tType, IType, pType, tAType, dType, zType, T} <:
             z,
             extrapolate,
             Ref(1),
-            cache_parameters
+            cache_parameters,
+            linear_lookup
         )
     end
 end
 
 function QuadraticSpline(
         u::uType, t; extrapolate = false,
-        cache_parameters = false) where {uType <: AbstractVector{<:Number}}
+        cache_parameters = false, assume_linear_t = 1e-2) where {uType <:
+                                                                 AbstractVector{<:Number}}
     u, t = munge_data(u, t)
+    linear_lookup = seems_linear(assume_linear_t, t)
     s = length(t)
     dl = ones(eltype(t), s - 1)
     d_tmp = ones(eltype(t), s)
@@ -314,15 +350,18 @@ function QuadraticSpline(
     z = tA \ d
 
     p = QuadraticSplineParameterCache(z, t, cache_parameters)
-    A = QuadraticSpline(u, t, nothing, p, tA, d, z, extrapolate, cache_parameters)
+    A = QuadraticSpline(
+        u, t, nothing, p, tA, d, z, extrapolate, cache_parameters, linear_lookup)
     I = cumulative_integral(A, cache_parameters)
-    QuadraticSpline(u, t, I, p, tA, d, z, extrapolate, cache_parameters)
+    QuadraticSpline(u, t, I, p, tA, d, z, extrapolate, cache_parameters, linear_lookup)
 end
 
 function QuadraticSpline(
-        u::uType, t; extrapolate = false, cache_parameters = false) where {uType <:
-                                                                           AbstractVector}
+        u::uType, t; extrapolate = false, cache_parameters = false,
+        assume_linear_t = 1e-2) where {uType <:
+                                       AbstractVector}
     u, t = munge_data(u, t)
+    linear_lookup = seems_linear(assume_linear_t, t)
     s = length(t)
     dl = ones(eltype(t), s - 1)
     d_tmp = ones(eltype(t), s)
@@ -337,9 +376,10 @@ function QuadraticSpline(
     z = [z_s for z_s in eachslice(z_, dims = ndims(z_))]
 
     p = QuadraticSplineParameterCache(z, t, cache_parameters)
-    A = QuadraticSpline(u, t, nothing, p, tA, d, z, extrapolate, cache_parameters)
+    A = QuadraticSpline(
+        u, t, nothing, p, tA, d, z, extrapolate, cache_parameters, linear_lookup)
     I = cumulative_integral(A, cache_parameters)
-    QuadraticSpline(u, t, I, p, tA, d, z, extrapolate, cache_parameters)
+    QuadraticSpline(u, t, I, p, tA, d, z, extrapolate, cache_parameters, linear_lookup)
 end
 
 """
@@ -357,6 +397,10 @@ Second derivative on both ends are zero, which are also called "natural" boundar
 
   - `extrapolate`: boolean value to allow extrapolation. Defaults to `false`.
   - `cache_parameters`: precompute parameters at initialization for faster interpolation computations. Note: if activated, `u` and `t` should not be modified. Defaults to `false`.
+  - `assume_linear_t`: boolean value to specify a faster index lookup behaviour for
+    evenly-distributed abscissae. Alternatively, a numerical threshold may be specified
+    for a test based on the normalized standard deviation of the difference with respect
+    to the straight line. Defaults to 1e-2
 """
 struct CubicSpline{uType, tType, IType, pType, hType, zType, T} <: AbstractInterpolation{T}
     u::uType
@@ -369,7 +413,8 @@ struct CubicSpline{uType, tType, IType, pType, hType, zType, T} <: AbstractInter
     idx_prev::Base.RefValue{Int}
     cache_parameters::Bool
     use_linear_lookup::Bool
-    function CubicSpline(u, t, I, p, h, z, extrapolate, cache_parameters)
+    function CubicSpline(u, t, I, p, h, z, extrapolate, cache_parameters, assume_linear_t)
+        linear_lookup = seems_linear(assume_linear_t, t)
         new{typeof(u), typeof(t), typeof(I), typeof(p.c₁), typeof(h), typeof(z), eltype(u)}(
             u,
             t,
@@ -379,15 +424,17 @@ struct CubicSpline{uType, tType, IType, pType, hType, zType, T} <: AbstractInter
             z,
             extrapolate,
             Ref(1),
-            cache_parameters
+            cache_parameters,
+            linear_lookup
         )
     end
 end
 
 function CubicSpline(u::uType,
         t;
-        extrapolate = false, cache_parameters = false) where {uType <:
-                                                              AbstractVector{<:Number}}
+        extrapolate = false, cache_parameters = false,
+        assume_linear_t = 1e-2) where {uType <:
+                                       AbstractVector{<:Number}}
     u, t = munge_data(u, t)
     n = length(t) - 1
     h = vcat(0, map(k -> t[k + 1] - t[k], 1:(length(t) - 1)), 0)
@@ -405,16 +452,18 @@ function CubicSpline(u::uType,
              6(u[i + 1] - u[i]) / h[i + 1] - 6(u[i] - u[i - 1]) / h[i],
         1:(n + 1))
     z = tA \ d
-
+    linear_lookup = seems_linear(assume_linear_t, t)
     p = CubicSplineParameterCache(u, h, z, cache_parameters)
-    A = CubicSpline(u, t, nothing, p, h[1:(n + 1)], z, extrapolate, cache_parameters)
+    A = CubicSpline(
+        u, t, nothing, p, h[1:(n + 1)], z, extrapolate, cache_parameters, linear_lookup)
     I = cumulative_integral(A, cache_parameters)
-    CubicSpline(u, t, I, p, h[1:(n + 1)], z, extrapolate, cache_parameters)
+    CubicSpline(u, t, I, p, h[1:(n + 1)], z, extrapolate, cache_parameters, linear_lookup)
 end
 
 function CubicSpline(
-        u::uType, t; extrapolate = false, cache_parameters = false) where {uType <:
-                                                                           AbstractVector}
+        u::uType, t; extrapolate = false, cache_parameters = false,
+        assume_linear_t = 1e-2) where {uType <:
+                                       AbstractVector}
     u, t = munge_data(u, t)
     n = length(t) - 1
     h = vcat(0, map(k -> t[k + 1] - t[k], 1:(length(t) - 1)), 0)
@@ -431,9 +480,10 @@ function CubicSpline(
     z = [z_s for z_s in eachslice(z_, dims = ndims(z_))]
 
     p = CubicSplineParameterCache(u, h, z, cache_parameters)
-    A = CubicSpline(u, t, nothing, p, h[1:(n + 1)], z, extrapolate, cache_parameters)
+    A = CubicSpline(
+        u, t, nothing, p, h[1:(n + 1)], z, extrapolate, cache_parameters, assume_linear_t)
     I = cumulative_integral(A, cache_parameters)
-    CubicSpline(u, t, I, p, h[1:(n + 1)], z, extrapolate, cache_parameters)
+    CubicSpline(u, t, I, p, h[1:(n + 1)], z, extrapolate, cache_parameters, assume_linear_t)
 end
 
 """
@@ -453,6 +503,10 @@ Extrapolation is a constant polynomial of the end points on each side.
 ## Keyword Arguments
 
   - `extrapolate`: boolean value to allow extrapolation. Defaults to `false`.
+  - `assume_linear_t`: boolean value to specify a faster index lookup behaviour for
+    evenly-distributed abscissae. Alternatively, a numerical threshold may be specified
+    for a test based on the normalized standard deviation of the difference with respect
+    to the straight line. Defaults to 1e-2
 """
 struct BSplineInterpolation{uType, tType, pType, kType, cType, NType, T} <:
        AbstractInterpolation{T}
@@ -477,7 +531,9 @@ struct BSplineInterpolation{uType, tType, pType, kType, cType, NType, T} <:
             N,
             pVecType,
             knotVecType,
-            extrapolate)
+            extrapolate,
+            assume_linear_t)
+        linear_lookup = seems_linear(assume_linear_t, t)
         new{typeof(u), typeof(t), typeof(p), typeof(k), typeof(c), typeof(N), eltype(u)}(u,
             t,
             d,
@@ -488,13 +544,14 @@ struct BSplineInterpolation{uType, tType, pType, kType, cType, NType, T} <:
             pVecType,
             knotVecType,
             extrapolate,
-            Ref(1)
+            Ref(1),
+            linear_lookup
         )
     end
 end
 
 function BSplineInterpolation(
-        u, t, d, pVecType, knotVecType; extrapolate = false)
+        u, t, d, pVecType, knotVecType; extrapolate = false, assume_linear_t = 1e-2)
     u, t = munge_data(u, t)
     n = length(t)
     n < d + 1 && error("BSplineInterpolation needs at least d + 1, i.e. $(d+1) points.")
@@ -559,7 +616,7 @@ function BSplineInterpolation(
     c = vec(N \ u[:, :])
     N = zeros(eltype(t), n)
     BSplineInterpolation(
-        u, t, d, p, k, c, N, pVecType, knotVecType, extrapolate)
+        u, t, d, p, k, c, N, pVecType, knotVecType, extrapolate, assume_linear_t)
 end
 
 """
@@ -581,6 +638,10 @@ Extrapolation is a constant polynomial of the end points on each side.
 ## Keyword Arguments
 
   - `extrapolate`: boolean value to allow extrapolation. Defaults to `false`.
+  - `assume_linear_t`: boolean value to specify a faster index lookup behaviour for
+    evenly-distributed abscissae. Alternatively, a numerical threshold may be specified
+    for a test based on the normalized standard deviation of the difference with respect
+    to the straight line. Defaults to 1e-2
 """
 struct BSplineApprox{uType, tType, pType, kType, cType, NType, T} <:
        AbstractInterpolation{T}
@@ -607,8 +668,10 @@ struct BSplineApprox{uType, tType, pType, kType, cType, NType, T} <:
             N,
             pVecType,
             knotVecType,
-            extrapolate
+            extrapolate,
+            assume_linear_t
     )
+        linear_lookup = seems_linear(assume_linear_t, t)
         new{typeof(u), typeof(t), typeof(p), typeof(k), typeof(c), typeof(N), eltype(u)}(u,
             t,
             d,
@@ -620,13 +683,14 @@ struct BSplineApprox{uType, tType, pType, kType, cType, NType, T} <:
             pVecType,
             knotVecType,
             extrapolate,
-            Ref(1)
+            Ref(1),
+            linear_lookup
         )
     end
 end
 
 function BSplineApprox(
-        u, t, d, h, pVecType, knotVecType; extrapolate = false)
+        u, t, d, h, pVecType, knotVecType; extrapolate = false, assume_linear_t = 1e-2)
     u, t = munge_data(u, t)
     n = length(t)
     h < d + 1 && error("BSplineApprox needs at least d + 1, i.e. $(d+1) control points.")
@@ -712,7 +776,7 @@ function BSplineApprox(
     c[2:(end - 1)] .= vec(P)
     N = zeros(eltype(t), h)
     BSplineApprox(
-        u, t, d, h, p, k, c, N, pVecType, knotVecType, extrapolate)
+        u, t, d, h, p, k, c, N, pVecType, knotVecType, extrapolate, assume_linear_t)
 end
 
 """
@@ -730,6 +794,10 @@ It is a Cubic Hermite interpolation, which is a piece-wise third degree polynomi
 
   - `extrapolate`: boolean value to allow extrapolation. Defaults to `false`.
   - `cache_parameters`: precompute parameters at initialization for faster interpolation computations. Note: if activated, `u` and `t` should not be modified. Defaults to `false`.
+  - `assume_linear_t`: boolean value to specify a faster index lookup behaviour for
+    evenly-distributed abscissae. Alternatively, a numerical threshold may be specified
+    for a test based on the normalized standard deviation of the difference with respect
+    to the straight line. Defaults to 1e-2
 """
 struct CubicHermiteSpline{uType, tType, IType, duType, pType, T} <: AbstractInterpolation{T}
     du::duType
@@ -741,19 +809,24 @@ struct CubicHermiteSpline{uType, tType, IType, duType, pType, T} <: AbstractInte
     idx_prev::Base.RefValue{Int}
     cache_parameters::Bool
     use_linear_lookup::Bool
-    function CubicHermiteSpline(du, u, t, I, p, extrapolate, cache_parameters)
+    function CubicHermiteSpline(
+            du, u, t, I, p, extrapolate, cache_parameters, assume_linear_t)
+        linear_lookup = seems_linear(assume_linear_t, t)
         new{typeof(u), typeof(t), typeof(I), typeof(du), typeof(p.c₁), eltype(u)}(
-            du, u, t, I, p, extrapolate, Ref(1), cache_parameters)
+            du, u, t, I, p, extrapolate, Ref(1), cache_parameters, linear_lookup)
     end
 end
 
-function CubicHermiteSpline(du, u, t; extrapolate = false, cache_parameters = false)
+function CubicHermiteSpline(
+        du, u, t; extrapolate = false, cache_parameters = false, assume_linear_t = 1e-2)
     @assert length(u)==length(du) "Length of `u` is not equal to length of `du`."
     u, t = munge_data(u, t)
+    linear_lookup = seems_linear(assume_linear_t, t)
     p = CubicHermiteParameterCache(du, u, t, cache_parameters)
-    A = CubicHermiteSpline(du, u, t, nothing, p, extrapolate, cache_parameters)
+    A = CubicHermiteSpline(
+        du, u, t, nothing, p, extrapolate, cache_parameters, linear_lookup)
     I = cumulative_integral(A, cache_parameters)
-    CubicHermiteSpline(du, u, t, I, p, extrapolate, cache_parameters)
+    CubicHermiteSpline(du, u, t, I, p, extrapolate, cache_parameters, linear_lookup)
 end
 
 """
@@ -772,11 +845,16 @@ section 3.4 for more details.
 
   - `extrapolate`: boolean value to allow extrapolation. Defaults to `false`.
   - `cache_parameters`: precompute parameters at initialization for faster interpolation computations. Note: if activated, `u` and `t` should not be modified. Defaults to `false`.
+  - `assume_linear_t`: boolean value to specify a faster index lookup behaviour for
+    evenly-distributed abscissae. Alternatively, a numerical threshold may be specified
+    for a test based on the normalized standard deviation of the difference with respect
+    to the straight line. Defaults to 1e-2
 """
-function PCHIPInterpolation(u, t; extrapolate = false, cache_parameters = false)
+function PCHIPInterpolation(
+        u, t; extrapolate = false, cache_parameters = false, assume_linear_t = 1e-2)
     u, t = munge_data(u, t)
     du = du_PCHIP(u, t)
-    CubicHermiteSpline(du, u, t; extrapolate, cache_parameters)
+    CubicHermiteSpline(du, u, t; extrapolate, cache_parameters, assume_linear_t)
 end
 
 """
@@ -795,6 +873,10 @@ It is a Quintic Hermite interpolation, which is a piece-wise fifth degree polyno
 
   - `extrapolate`: boolean value to allow extrapolation. Defaults to `false`.
   - `cache_parameters`: precompute parameters at initialization for faster interpolation computations. Note: if activated, `u` and `t` should not be modified. Defaults to `false`.
+  - `assume_linear_t`: boolean value to specify a faster index lookup behaviour for
+    evenly-distributed abscissae. Alternatively, a numerical threshold may be specified
+    for a test based on the normalized standard deviation of the difference with respect
+    to the straight line. Defaults to 1e-2
 """
 struct QuinticHermiteSpline{uType, tType, IType, duType, dduType, pType, T} <:
        AbstractInterpolation{T}
@@ -808,18 +890,24 @@ struct QuinticHermiteSpline{uType, tType, IType, duType, dduType, pType, T} <:
     idx_prev::Base.RefValue{Int}
     cache_parameters::Bool
     use_linear_lookup::Bool
-    function QuinticHermiteSpline(ddu, du, u, t, I, p, extrapolate, cache_parameters)
+    function QuinticHermiteSpline(
+            ddu, du, u, t, I, p, extrapolate, cache_parameters, assume_linear_t)
+        linear_lookup = seems_linear(assume_linear_t, t)
         new{typeof(u), typeof(t), typeof(I), typeof(du),
             typeof(ddu), typeof(p.c₁), eltype(u)}(
-            ddu, du, u, t, I, p, extrapolate, Ref(1), cache_parameters)
+            ddu, du, u, t, I, p, extrapolate, Ref(1), cache_parameters, linear_lookup)
     end
 end
 
-function QuinticHermiteSpline(ddu, du, u, t; extrapolate = false, cache_parameters = false)
+function QuinticHermiteSpline(ddu, du, u, t; extrapolate = false,
+        cache_parameters = false, assume_linear_t = 1e-2)
     @assert length(u)==length(du)==length(ddu) "Length of `u` is not equal to length of `du` or `ddu`."
     u, t = munge_data(u, t)
+    linear_lookup = seems_linear(assume_linear_t, t)
     p = QuinticHermiteParameterCache(ddu, du, u, t, cache_parameters)
-    A = QuinticHermiteSpline(ddu, du, u, t, nothing, p, extrapolate, cache_parameters)
+    A = QuinticHermiteSpline(
+        ddu, du, u, t, nothing, p, extrapolate, cache_parameters, linear_lookup)
     I = cumulative_integral(A, cache_parameters)
-    QuinticHermiteSpline(ddu, du, u, t, I, p, extrapolate, cache_parameters)
+    QuinticHermiteSpline(
+        ddu, du, u, t, I, p, extrapolate, cache_parameters, linear_lookup)
 end
diff --git a/src/interpolation_methods.jl b/src/interpolation_methods.jl
index 99f3874d..69a5cf94 100644
--- a/src/interpolation_methods.jl
+++ b/src/interpolation_methods.jl
@@ -1,14 +1,7 @@
 function _interpolate(A, t)
     ((t < A.t[1] || t > A.t[end]) && !A.extrapolate) &&
         throw(ExtrapolationError())
-    idx_guess = if hasfield(typeof(A), :use_linear_lookup) &&
-                   A.use_linear_lookup
-        f = (t - first(A.t)) / (last(A.t) - first(A.t))
-        i_0, i_f = firstindex(A.t), lastindex(A.t)
-        round(typeof(firstindex(A.t)), f * (i_f - i_0) + i_0)
-    else
-        A.idx_prev[]
-    end
+    idx_guess = A.idx_prev[]
     val, idx_prev = _interpolate(A, t, idx_guess)
     A.idx_prev[] = idx_prev
     return val
@@ -23,7 +16,7 @@ function _interpolate(A::LinearInterpolation{<:AbstractVector}, t::Number, igues
         u1 = u2 = one(eltype(A.u))
         slope = t * get_parameters(A, idx)
     else
-        idx = get_idx(A.t, t, iguess)
+        idx = get_idx(A, t, iguess)
         t1, t2 = A.t[idx], A.t[idx + 1]
         u1, u2 = A.u[idx], A.u[idx + 1]
         slope = get_parameters(A, idx)
@@ -44,7 +37,7 @@ function _interpolate(A::LinearInterpolation{<:AbstractVector}, t::Number, igues
 end
 
 function _interpolate(A::LinearInterpolation{<:AbstractMatrix}, t::Number, iguess)
-    idx = get_idx(A.t, t, iguess)
+    idx = get_idx(A, t, iguess)
     Δt = t - A.t[idx]
     slope = get_parameters(A, idx)
     return A.u[:, idx] + slope * Δt, idx
@@ -53,7 +46,7 @@ end
 # Quadratic Interpolation
 _quad_interp_indices(A, t) = _quad_interp_indices(A, t, firstindex(A.t) - 1)
 function _quad_interp_indices(A::QuadraticInterpolation, t::Number, iguess)
-    idx = get_idx(A.t, t, iguess; idx_shift = A.mode == :Backward ? -1 : 0, ub_shift = -2)
+    idx = get_idx(A, t, iguess; idx_shift = A.mode == :Backward ? -1 : 0, ub_shift = -2)
     idx, idx + 1, idx + 2
 end
 
@@ -68,7 +61,7 @@ end
 
 # Lagrange Interpolation
 function _interpolate(A::LagrangeInterpolation{<:AbstractVector}, t::Number, iguess)
-    idx = get_idx(A.t, t, iguess)
+    idx = get_idx(A, t, iguess)
     findRequiredIdxs!(A, t, idx)
     if A.t[A.idxs[1]] == t
         return A.u[A.idxs[1]], idx
@@ -97,7 +90,7 @@ function _interpolate(A::LagrangeInterpolation{<:AbstractVector}, t::Number, igu
 end
 
 function _interpolate(A::LagrangeInterpolation{<:AbstractMatrix}, t::Number, iguess)
-    idx = get_idx(A.t, t, iguess)
+    idx = get_idx(A, t, iguess)
     findRequiredIdxs!(A, t, idx)
     if A.t[A.idxs[1]] == t
         return A.u[:, A.idxs[1]], idx
@@ -126,7 +119,7 @@ function _interpolate(A::LagrangeInterpolation{<:AbstractMatrix}, t::Number, igu
 end
 
 function _interpolate(A::AkimaInterpolation{<:AbstractVector}, t::Number, iguess)
-    idx = get_idx(A.t, t, iguess)
+    idx = get_idx(A, t, iguess)
     wj = t - A.t[idx]
     (@evalpoly wj A.u[idx] A.b[idx] A.c[idx] A.d[idx]), idx
 end
@@ -135,10 +128,10 @@ end
 function _interpolate(A::ConstantInterpolation{<:AbstractVector}, t::Number, iguess)
     if A.dir === :left
         # :left means that value to the left is used for interpolation
-        idx = get_idx(A.t, t, iguess; lb = 1, ub_shift = 0)
+        idx = get_idx(A, t, iguess; lb = 1, ub_shift = 0)
     else
         # :right means that value to the right is used for interpolation
-        idx = get_idx(A.t, t, iguess; side = :first, lb = 1, ub_shift = 0)
+        idx = get_idx(A, t, iguess; side = :first, lb = 1, ub_shift = 0)
     end
     A.u[idx], idx
 end
@@ -146,17 +139,17 @@ end
 function _interpolate(A::ConstantInterpolation{<:AbstractMatrix}, t::Number, iguess)
     if A.dir === :left
         # :left means that value to the left is used for interpolation
-        idx = get_idx(A.t, t, iguess; lb = 1, ub_shift = 0)
+        idx = get_idx(A, t, iguess; lb = 1, ub_shift = 0)
     else
         # :right means that value to the right is used for interpolation
-        idx = get_idx(A.t, t, iguess; side = :first, lb = 1, ub_shift = 0)
+        idx = get_idx(A, t, iguess; side = :first, lb = 1, ub_shift = 0)
     end
     A.u[:, idx], idx
 end
 
 # QuadraticSpline Interpolation
 function _interpolate(A::QuadraticSpline{<:AbstractVector}, t::Number, iguess)
-    idx = get_idx(A.t, t, iguess)
+    idx = get_idx(A, t, iguess)
     Cᵢ = A.u[idx]
     Δt = t - A.t[idx]
     σ = get_parameters(A, idx)
@@ -165,7 +158,7 @@ end
 
 # CubicSpline Interpolation
 function _interpolate(A::CubicSpline{<:AbstractVector}, t::Number, iguess)
-    idx = get_idx(A.t, t, iguess)
+    idx = get_idx(A, t, iguess)
     Δt₁ = t - A.t[idx]
     Δt₂ = A.t[idx + 1] - t
     I = (A.z[idx] * Δt₂^3 + A.z[idx + 1] * Δt₁^3) / (6A.h[idx + 1])
@@ -182,7 +175,7 @@ function _interpolate(A::BSplineInterpolation{<:AbstractVector{<:Number}},
     t < A.t[1] && return A.u[1], 1
     t > A.t[end] && return A.u[end], lastindex(t)
     # change t into param [0 1]
-    idx = get_idx(A.t, t, iguess)
+    idx = get_idx(A, t, iguess)
     t = A.p[idx] + (t - A.t[idx]) / (A.t[idx + 1] - A.t[idx]) * (A.p[idx + 1] - A.p[idx])
     n = length(A.t)
     N = t isa ForwardDiff.Dual ? zeros(eltype(t), n) : A.N
@@ -199,7 +192,7 @@ function _interpolate(A::BSplineApprox{<:AbstractVector{<:Number}}, t::Number, i
     t < A.t[1] && return A.u[1], 1
     t > A.t[end] && return A.u[end], lastindex(t)
     # change t into param [0 1]
-    idx = get_idx(A.t, t, iguess)
+    idx = get_idx(A, t, iguess)
     t = A.p[idx] + (t - A.t[idx]) / (A.t[idx + 1] - A.t[idx]) * (A.p[idx + 1] - A.p[idx])
     N = t isa ForwardDiff.Dual ? zeros(eltype(t), A.h) : A.N
     nonzero_coefficient_idxs = spline_coefficients!(N, A.d, A.k, t)
@@ -213,7 +206,7 @@ end
 # Cubic Hermite Spline
 function _interpolate(
         A::CubicHermiteSpline{<:AbstractVector{<:Number}}, t::Number, iguess)
-    idx = get_idx(A.t, t, iguess)
+    idx = get_idx(A, t, iguess)
     Δt₀ = t - A.t[idx]
     Δt₁ = t - A.t[idx + 1]
     out = A.u[idx] + Δt₀ * A.du[idx]
@@ -225,7 +218,7 @@ end
 # Quintic Hermite Spline
 function _interpolate(
         A::QuinticHermiteSpline{<:AbstractVector{<:Number}}, t::Number, iguess)
-    idx = get_idx(A.t, t, iguess)
+    idx = get_idx(A, t, iguess)
     Δt₀ = t - A.t[idx]
     Δt₁ = t - A.t[idx + 1]
     out = A.u[idx] + Δt₀ * (A.du[idx] + A.ddu[idx] * Δt₀ / 2)
diff --git a/src/interpolation_utils.jl b/src/interpolation_utils.jl
index d29959dc..1ae3ca28 100644
--- a/src/interpolation_utils.jl
+++ b/src/interpolation_utils.jl
@@ -113,7 +113,17 @@ function looks_linear(t; threshold = 1e-2)
     norm_var < threshold^2
 end
 
-function get_idx(tvec, t, iguess; lb = 1, ub_shift = -1, idx_shift = 0, side = :last)
+function get_idx(A::AbstractInterpolation, t, iguess; lb = 1,
+        ub_shift = -1, idx_shift = 0, side = :last)
+    iguess = if hasfield(typeof(A), :use_linear_lookup) &&
+                A.use_linear_lookup
+        f = (t - first(A.t)) / (last(A.t) - first(A.t))
+        i_0, i_f = firstindex(A.t), lastindex(A.t)
+        round(typeof(firstindex(A.t)), f * (i_f - i_0) + i_0)
+    else
+        iguess
+    end
+    tvec = A.t
     ub = length(tvec) + ub_shift
     return if side == :last
         clamp(searchsortedlastcorrelated(tvec, t, iguess) + idx_shift, lb, ub)

From 80ef84c1ea08c7b50f3b3e285709eeaa81b04cec Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Daniel=20Gonz=C3=A1lez?= <daniel.glez.arribas@gmail.com>
Date: Mon, 29 Jul 2024 22:40:22 +0200
Subject: [PATCH 06/11] Rename argument for consistency

---
 src/interpolation_utils.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/interpolation_utils.jl b/src/interpolation_utils.jl
index 1ae3ca28..9cbb69d6 100644
--- a/src/interpolation_utils.jl
+++ b/src/interpolation_utils.jl
@@ -95,7 +95,7 @@ function munge_data(U::StridedMatrix, t::AbstractVector)
 end
 
 seems_linear(assume_linear_t::Bool, _) = assume_linear_t
-seems_linear(assume_qr::Number, t) = looks_linear(t; threshold = assume_qr)
+seems_linear(assume_linear_t::Number, t) = looks_linear(t; threshold = assume_linear_t)
 
 # """
 # Determine if the abscissae are regularly distributed, taking the standard \

From 283775123b73e2616a430f3bb70496112a58299d Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Daniel=20Gonz=C3=A1lez?= <daniel.glez.arribas@gmail.com>
Date: Tue, 30 Jul 2024 13:46:18 +0200
Subject: [PATCH 07/11] Incorporate Sathvik's suggestions

---
 Project.toml                |  4 +++-
 docs/src/manual.md          |  1 +
 src/interpolation_caches.jl | 42 ++++++++++++++++++-------------------
 src/interpolation_utils.jl  | 14 ++++++-------
 test/interpolation_tests.jl | 18 ++++++++++++++++
 5 files changed, 50 insertions(+), 29 deletions(-)

diff --git a/Project.toml b/Project.toml
index 17fdd131..ee75d23d 100644
--- a/Project.toml
+++ b/Project.toml
@@ -25,6 +25,7 @@ DataInterpolationsSymbolicsExt = "Symbolics"
 
 [compat]
 Aqua = "0.8"
+BenchmarkTools = "1"
 ChainRulesCore = "1.24"
 FindFirstFunctions = "1.1"
 FiniteDifferences = "0.12.31"
@@ -45,6 +46,7 @@ julia = "1.10"
 
 [extras]
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
+BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
 FiniteDifferences = "26cc04aa-876d-5657-8c51-4c34ba976000"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
@@ -58,4 +60,4 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
 
 [targets]
-test = ["Aqua", "SafeTestsets", "ChainRulesCore", "Optim", "RegularizationTools", "Test", "StableRNGs", "FiniteDifferences", "QuadGK", "ForwardDiff", "Symbolics", "Zygote"]
+test = ["Aqua", "BenchmarkTools", "SafeTestsets", "ChainRulesCore", "Optim", "RegularizationTools", "Test", "StableRNGs", "FiniteDifferences", "QuadGK", "ForwardDiff", "Symbolics", "Zygote"]
diff --git a/docs/src/manual.md b/docs/src/manual.md
index 077cf29f..cf1cbc15 100644
--- a/docs/src/manual.md
+++ b/docs/src/manual.md
@@ -13,4 +13,5 @@ BSplineApprox
 CubicHermiteSpline
 PCHIPInterpolation
 QuinticHermiteSpline
+looks_linear
 ```
diff --git a/src/interpolation_caches.jl b/src/interpolation_caches.jl
index c4017460..82a53a13 100644
--- a/src/interpolation_caches.jl
+++ b/src/interpolation_caches.jl
@@ -17,7 +17,7 @@ Extrapolation extends the last linear polynomial on each side.
   - `assume_linear_t`: boolean value to specify a faster index lookup behaviour for
     evenly-distributed abscissae. Alternatively, a numerical threshold may be specified
     for a test based on the normalized standard deviation of the difference with respect
-    to the straight line. Defaults to 1e-2
+    to the straight line (see [`looks_linear`](@ref)). Defaults to 1e-2.
 """
 struct LinearInterpolation{uType, tType, IType, pType, T} <: AbstractInterpolation{T}
     u::uType
@@ -27,7 +27,7 @@ struct LinearInterpolation{uType, tType, IType, pType, T} <: AbstractInterpolati
     extrapolate::Bool
     idx_prev::Base.RefValue{Int}
     cache_parameters::Bool
-    use_linear_lookup::Bool
+    linear_lookup::Bool
     function LinearInterpolation(u, t, I, p, extrapolate, cache_parameters, assume_linear_t)
         linear_lookup = seems_linear(assume_linear_t, t)
         new{typeof(u), typeof(t), typeof(I), typeof(p.slope), eltype(u)}(
@@ -64,7 +64,7 @@ Extrapolation extends the last quadratic polynomial on each side.
   - `assume_linear_t`: boolean value to specify a faster index lookup behaviour for
     evenly-distributed abscissae. Alternatively, a numerical threshold may be specified
     for a test based on the normalized standard deviation of the difference with respect
-    to the straight line. Defaults to 1e-2
+    to the straight line. Defaults to 1e-2.
 """
 struct QuadraticInterpolation{uType, tType, IType, pType, T} <: AbstractInterpolation{T}
     u::uType
@@ -75,7 +75,7 @@ struct QuadraticInterpolation{uType, tType, IType, pType, T} <: AbstractInterpol
     extrapolate::Bool
     idx_prev::Base.RefValue{Int}
     cache_parameters::Bool
-    use_linear_lookup::Bool
+    linear_lookup::Bool
     function QuadraticInterpolation(
             u, t, I, p, mode, extrapolate, cache_parameters, assume_linear_t)
         mode ∈ (:Forward, :Backward) ||
@@ -167,7 +167,7 @@ Extrapolation extends the last cubic polynomial on each side.
   - `assume_linear_t`: boolean value to specify a faster index lookup behaviour for
     evenly-distributed abscissae. Alternatively, a numerical threshold may be specified
     for a test based on the normalized standard deviation of the difference with respect
-    to the straight line. Defaults to 1e-2
+    to the straight line. Defaults to 1e-2.
 """
 struct AkimaInterpolation{uType, tType, IType, bType, cType, dType, T} <:
        AbstractInterpolation{T}
@@ -180,7 +180,7 @@ struct AkimaInterpolation{uType, tType, IType, bType, cType, dType, T} <:
     extrapolate::Bool
     idx_prev::Base.RefValue{Int}
     cache_parameters::Bool
-    use_linear_lookup::Bool
+    linear_lookup::Bool
     function AkimaInterpolation(
             u, t, I, b, c, d, extrapolate, cache_parameters, assume_linear_t)
         linear_lookup = seems_linear(assume_linear_t, t)
@@ -249,7 +249,7 @@ Extrapolation extends the last constant polynomial at the end points on each sid
   - `assume_linear_t`: boolean value to specify a faster index lookup behaviour for
     evenly-distributed abscissae. Alternatively, a numerical threshold may be specified
     for a test based on the normalized standard deviation of the difference with respect
-    to the straight line. Defaults to 1e-2
+    to the straight line. Defaults to 1e-2.
 """
 struct ConstantInterpolation{uType, tType, IType, T} <: AbstractInterpolation{T}
     u::uType
@@ -260,7 +260,7 @@ struct ConstantInterpolation{uType, tType, IType, T} <: AbstractInterpolation{T}
     extrapolate::Bool
     idx_prev::Base.RefValue{Int}
     cache_parameters::Bool
-    use_linear_lookup::Bool
+    linear_lookup::Bool
     function ConstantInterpolation(
             u, t, I, dir, extrapolate, cache_parameters, assume_linear_t)
         linear_lookup = seems_linear(assume_linear_t, t)
@@ -297,7 +297,7 @@ Extrapolation extends the last quadratic polynomial on each side.
   - `assume_linear_t`: boolean value to specify a faster index lookup behaviour for
     evenly-distributed abscissae. Alternatively, a numerical threshold may be specified
     for a test based on the normalized standard deviation of the difference with respect
-    to the straight line. Defaults to 1e-2
+    to the straight line. Defaults to 1e-2.
 """
 struct QuadraticSpline{uType, tType, IType, pType, tAType, dType, zType, T} <:
        AbstractInterpolation{T}
@@ -311,7 +311,7 @@ struct QuadraticSpline{uType, tType, IType, pType, tAType, dType, zType, T} <:
     extrapolate::Bool
     idx_prev::Base.RefValue{Int}
     cache_parameters::Bool
-    use_linear_lookup::Bool
+    linear_lookup::Bool
     function QuadraticSpline(
             u, t, I, p, tA, d, z, extrapolate, cache_parameters, assume_linear_t)
         linear_lookup = seems_linear(assume_linear_t, t)
@@ -400,7 +400,7 @@ Second derivative on both ends are zero, which are also called "natural" boundar
   - `assume_linear_t`: boolean value to specify a faster index lookup behaviour for
     evenly-distributed abscissae. Alternatively, a numerical threshold may be specified
     for a test based on the normalized standard deviation of the difference with respect
-    to the straight line. Defaults to 1e-2
+    to the straight line. Defaults to 1e-2.
 """
 struct CubicSpline{uType, tType, IType, pType, hType, zType, T} <: AbstractInterpolation{T}
     u::uType
@@ -412,7 +412,7 @@ struct CubicSpline{uType, tType, IType, pType, hType, zType, T} <: AbstractInter
     extrapolate::Bool
     idx_prev::Base.RefValue{Int}
     cache_parameters::Bool
-    use_linear_lookup::Bool
+    linear_lookup::Bool
     function CubicSpline(u, t, I, p, h, z, extrapolate, cache_parameters, assume_linear_t)
         linear_lookup = seems_linear(assume_linear_t, t)
         new{typeof(u), typeof(t), typeof(I), typeof(p.c₁), typeof(h), typeof(z), eltype(u)}(
@@ -506,7 +506,7 @@ Extrapolation is a constant polynomial of the end points on each side.
   - `assume_linear_t`: boolean value to specify a faster index lookup behaviour for
     evenly-distributed abscissae. Alternatively, a numerical threshold may be specified
     for a test based on the normalized standard deviation of the difference with respect
-    to the straight line. Defaults to 1e-2
+    to the straight line. Defaults to 1e-2.
 """
 struct BSplineInterpolation{uType, tType, pType, kType, cType, NType, T} <:
        AbstractInterpolation{T}
@@ -521,7 +521,7 @@ struct BSplineInterpolation{uType, tType, pType, kType, cType, NType, T} <:
     knotVecType::Symbol
     extrapolate::Bool
     idx_prev::Base.RefValue{Int}
-    use_linear_lookup::Bool
+    linear_lookup::Bool
     function BSplineInterpolation(u,
             t,
             d,
@@ -641,7 +641,7 @@ Extrapolation is a constant polynomial of the end points on each side.
   - `assume_linear_t`: boolean value to specify a faster index lookup behaviour for
     evenly-distributed abscissae. Alternatively, a numerical threshold may be specified
     for a test based on the normalized standard deviation of the difference with respect
-    to the straight line. Defaults to 1e-2
+    to the straight line. Defaults to 1e-2.
 """
 struct BSplineApprox{uType, tType, pType, kType, cType, NType, T} <:
        AbstractInterpolation{T}
@@ -657,7 +657,7 @@ struct BSplineApprox{uType, tType, pType, kType, cType, NType, T} <:
     knotVecType::Symbol
     extrapolate::Bool
     idx_prev::Base.RefValue{Int}
-    use_linear_lookup::Bool
+    linear_lookup::Bool
     function BSplineApprox(u,
             t,
             d,
@@ -797,7 +797,7 @@ It is a Cubic Hermite interpolation, which is a piece-wise third degree polynomi
   - `assume_linear_t`: boolean value to specify a faster index lookup behaviour for
     evenly-distributed abscissae. Alternatively, a numerical threshold may be specified
     for a test based on the normalized standard deviation of the difference with respect
-    to the straight line. Defaults to 1e-2
+    to the straight line. Defaults to 1e-2.
 """
 struct CubicHermiteSpline{uType, tType, IType, duType, pType, T} <: AbstractInterpolation{T}
     du::duType
@@ -808,7 +808,7 @@ struct CubicHermiteSpline{uType, tType, IType, duType, pType, T} <: AbstractInte
     extrapolate::Bool
     idx_prev::Base.RefValue{Int}
     cache_parameters::Bool
-    use_linear_lookup::Bool
+    linear_lookup::Bool
     function CubicHermiteSpline(
             du, u, t, I, p, extrapolate, cache_parameters, assume_linear_t)
         linear_lookup = seems_linear(assume_linear_t, t)
@@ -848,7 +848,7 @@ section 3.4 for more details.
   - `assume_linear_t`: boolean value to specify a faster index lookup behaviour for
     evenly-distributed abscissae. Alternatively, a numerical threshold may be specified
     for a test based on the normalized standard deviation of the difference with respect
-    to the straight line. Defaults to 1e-2
+    to the straight line. Defaults to 1e-2.
 """
 function PCHIPInterpolation(
         u, t; extrapolate = false, cache_parameters = false, assume_linear_t = 1e-2)
@@ -876,7 +876,7 @@ It is a Quintic Hermite interpolation, which is a piece-wise fifth degree polyno
   - `assume_linear_t`: boolean value to specify a faster index lookup behaviour for
     evenly-distributed abscissae. Alternatively, a numerical threshold may be specified
     for a test based on the normalized standard deviation of the difference with respect
-    to the straight line. Defaults to 1e-2
+    to the straight line. Defaults to 1e-2.
 """
 struct QuinticHermiteSpline{uType, tType, IType, duType, dduType, pType, T} <:
        AbstractInterpolation{T}
@@ -889,7 +889,7 @@ struct QuinticHermiteSpline{uType, tType, IType, duType, dduType, pType, T} <:
     extrapolate::Bool
     idx_prev::Base.RefValue{Int}
     cache_parameters::Bool
-    use_linear_lookup::Bool
+    linear_lookup::Bool
     function QuinticHermiteSpline(
             ddu, du, u, t, I, p, extrapolate, cache_parameters, assume_linear_t)
         linear_lookup = seems_linear(assume_linear_t, t)
diff --git a/src/interpolation_utils.jl b/src/interpolation_utils.jl
index 9cbb69d6..3ee0cc5a 100644
--- a/src/interpolation_utils.jl
+++ b/src/interpolation_utils.jl
@@ -97,11 +97,11 @@ end
 seems_linear(assume_linear_t::Bool, _) = assume_linear_t
 seems_linear(assume_linear_t::Number, t) = looks_linear(t; threshold = assume_linear_t)
 
-# """
-# Determine if the abscissae are regularly distributed, taking the standard \
-# deviation of the difference between the array of abscissae with respect to \
-# the straight line linking its first and last elements
-# """
+"""
+Determine if the abscissae `t` are regularly distributed, taking the standard deviation of \
+the difference between the array of abscissae with respect to the straight line linking \
+its first and last elements, normalized by the range of `t`.
+"""
 function looks_linear(t; threshold = 1e-2)
     length(t) <= 2 && return true
     t_0, t_f = first(t), last(t)
@@ -115,8 +115,8 @@ end
 
 function get_idx(A::AbstractInterpolation, t, iguess; lb = 1,
         ub_shift = -1, idx_shift = 0, side = :last)
-    iguess = if hasfield(typeof(A), :use_linear_lookup) &&
-                A.use_linear_lookup
+    iguess = if hasfield(typeof(A), :linear_lookup) &&
+                A.linear_lookup
         f = (t - first(A.t)) / (last(A.t) - first(A.t))
         i_0, i_f = firstindex(A.t), lastindex(A.t)
         round(typeof(firstindex(A.t)), f * (i_f - i_0) + i_0)
diff --git a/test/interpolation_tests.jl b/test/interpolation_tests.jl
index 53fff25e..601047f6 100644
--- a/test/interpolation_tests.jl
+++ b/test/interpolation_tests.jl
@@ -2,6 +2,7 @@ using DataInterpolations
 using FindFirstFunctions: searchsortedfirstcorrelated
 using StableRNGs
 using Optim, ForwardDiff
+using BenchmarkTools
 
 function test_interpolation_type(T)
     @test T <: DataInterpolations.AbstractInterpolation
@@ -800,3 +801,20 @@ f_cubic_spline = c -> square(CubicSpline, c)
 @test ForwardDiff.derivative(f_quadratic_spline, 4.0) ≈ 8.0
 @test ForwardDiff.derivative(f_cubic_spline, 2.0) ≈ 4.0
 @test ForwardDiff.derivative(f_cubic_spline, 4.0) ≈ 8.0
+
+@testset "Linear lookup" begin
+    for N in (100, 1_000, 10_000, 100_000, 1_000_000) # Interpolant size
+        seed = 1234
+        t = collect(LinRange(0, 1, N)) # collect to avoid fast LinRange dispatch
+        y = rand(N)
+        A = LinearInterpolation(y, t)
+        A_fallback = LinearInterpolation(copy(y), copy(t); assume_linear_t = false)
+        @test A.linear_lookup
+        @test !(A_fallback.linear_lookup)
+
+        n_samples = 1_000
+        b_linear = @benchmark $A(rand()) samples=n_samples
+        b_fallback = @benchmark $A_fallback(rand()) samples=n_samples
+        @test mean(b_linear.times) < mean(b_fallback.times)
+    end
+end
\ No newline at end of file

From 4380c664fe9d9c29d33340a6b6462516bd1e1b47 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Daniel=20Gonz=C3=A1lez?= <daniel.glez.arribas@gmail.com>
Date: Tue, 30 Jul 2024 19:29:22 +0200
Subject: [PATCH 08/11] Update all docstrings to reference looks_linear

---
 src/interpolation_caches.jl | 20 ++++++++++----------
 1 file changed, 10 insertions(+), 10 deletions(-)

diff --git a/src/interpolation_caches.jl b/src/interpolation_caches.jl
index 82a53a13..71b2b042 100644
--- a/src/interpolation_caches.jl
+++ b/src/interpolation_caches.jl
@@ -64,7 +64,7 @@ Extrapolation extends the last quadratic polynomial on each side.
   - `assume_linear_t`: boolean value to specify a faster index lookup behaviour for
     evenly-distributed abscissae. Alternatively, a numerical threshold may be specified
     for a test based on the normalized standard deviation of the difference with respect
-    to the straight line. Defaults to 1e-2.
+    to the straight line (see [`looks_linear`](@ref)). Defaults to 1e-2.
 """
 struct QuadraticInterpolation{uType, tType, IType, pType, T} <: AbstractInterpolation{T}
     u::uType
@@ -167,7 +167,7 @@ Extrapolation extends the last cubic polynomial on each side.
   - `assume_linear_t`: boolean value to specify a faster index lookup behaviour for
     evenly-distributed abscissae. Alternatively, a numerical threshold may be specified
     for a test based on the normalized standard deviation of the difference with respect
-    to the straight line. Defaults to 1e-2.
+    to the straight line (see [`looks_linear`](@ref)). Defaults to 1e-2.
 """
 struct AkimaInterpolation{uType, tType, IType, bType, cType, dType, T} <:
        AbstractInterpolation{T}
@@ -249,7 +249,7 @@ Extrapolation extends the last constant polynomial at the end points on each sid
   - `assume_linear_t`: boolean value to specify a faster index lookup behaviour for
     evenly-distributed abscissae. Alternatively, a numerical threshold may be specified
     for a test based on the normalized standard deviation of the difference with respect
-    to the straight line. Defaults to 1e-2.
+    to the straight line (see [`looks_linear`](@ref)). Defaults to 1e-2.
 """
 struct ConstantInterpolation{uType, tType, IType, T} <: AbstractInterpolation{T}
     u::uType
@@ -297,7 +297,7 @@ Extrapolation extends the last quadratic polynomial on each side.
   - `assume_linear_t`: boolean value to specify a faster index lookup behaviour for
     evenly-distributed abscissae. Alternatively, a numerical threshold may be specified
     for a test based on the normalized standard deviation of the difference with respect
-    to the straight line. Defaults to 1e-2.
+    to the straight line (see [`looks_linear`](@ref)). Defaults to 1e-2.
 """
 struct QuadraticSpline{uType, tType, IType, pType, tAType, dType, zType, T} <:
        AbstractInterpolation{T}
@@ -400,7 +400,7 @@ Second derivative on both ends are zero, which are also called "natural" boundar
   - `assume_linear_t`: boolean value to specify a faster index lookup behaviour for
     evenly-distributed abscissae. Alternatively, a numerical threshold may be specified
     for a test based on the normalized standard deviation of the difference with respect
-    to the straight line. Defaults to 1e-2.
+    to the straight line (see [`looks_linear`](@ref)). Defaults to 1e-2.
 """
 struct CubicSpline{uType, tType, IType, pType, hType, zType, T} <: AbstractInterpolation{T}
     u::uType
@@ -506,7 +506,7 @@ Extrapolation is a constant polynomial of the end points on each side.
   - `assume_linear_t`: boolean value to specify a faster index lookup behaviour for
     evenly-distributed abscissae. Alternatively, a numerical threshold may be specified
     for a test based on the normalized standard deviation of the difference with respect
-    to the straight line. Defaults to 1e-2.
+    to the straight line (see [`looks_linear`](@ref)). Defaults to 1e-2.
 """
 struct BSplineInterpolation{uType, tType, pType, kType, cType, NType, T} <:
        AbstractInterpolation{T}
@@ -641,7 +641,7 @@ Extrapolation is a constant polynomial of the end points on each side.
   - `assume_linear_t`: boolean value to specify a faster index lookup behaviour for
     evenly-distributed abscissae. Alternatively, a numerical threshold may be specified
     for a test based on the normalized standard deviation of the difference with respect
-    to the straight line. Defaults to 1e-2.
+    to the straight line (see [`looks_linear`](@ref)). Defaults to 1e-2.
 """
 struct BSplineApprox{uType, tType, pType, kType, cType, NType, T} <:
        AbstractInterpolation{T}
@@ -797,7 +797,7 @@ It is a Cubic Hermite interpolation, which is a piece-wise third degree polynomi
   - `assume_linear_t`: boolean value to specify a faster index lookup behaviour for
     evenly-distributed abscissae. Alternatively, a numerical threshold may be specified
     for a test based on the normalized standard deviation of the difference with respect
-    to the straight line. Defaults to 1e-2.
+    to the straight line (see [`looks_linear`](@ref)). Defaults to 1e-2.
 """
 struct CubicHermiteSpline{uType, tType, IType, duType, pType, T} <: AbstractInterpolation{T}
     du::duType
@@ -848,7 +848,7 @@ section 3.4 for more details.
   - `assume_linear_t`: boolean value to specify a faster index lookup behaviour for
     evenly-distributed abscissae. Alternatively, a numerical threshold may be specified
     for a test based on the normalized standard deviation of the difference with respect
-    to the straight line. Defaults to 1e-2.
+    to the straight line (see [`looks_linear`](@ref)). Defaults to 1e-2.
 """
 function PCHIPInterpolation(
         u, t; extrapolate = false, cache_parameters = false, assume_linear_t = 1e-2)
@@ -876,7 +876,7 @@ It is a Quintic Hermite interpolation, which is a piece-wise fifth degree polyno
   - `assume_linear_t`: boolean value to specify a faster index lookup behaviour for
     evenly-distributed abscissae. Alternatively, a numerical threshold may be specified
     for a test based on the normalized standard deviation of the difference with respect
-    to the straight line. Defaults to 1e-2.
+    to the straight line (see [`looks_linear`](@ref)). Defaults to 1e-2.
 """
 struct QuinticHermiteSpline{uType, tType, IType, duType, dduType, pType, T} <:
        AbstractInterpolation{T}

From 179f1dbe9f10af5a88ff13d3c46f94237ad13b65 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Daniel=20Gonz=C3=A1lez?= <daniel.glez.arribas@gmail.com>
Date: Wed, 31 Jul 2024 15:39:07 +0200
Subject: [PATCH 09/11] Docstring changes for looks_linear

---
 docs/src/manual.md         | 7 ++++++-
 src/interpolation_utils.jl | 8 +++++---
 2 files changed, 11 insertions(+), 4 deletions(-)

diff --git a/docs/src/manual.md b/docs/src/manual.md
index cf1cbc15..08e4c32f 100644
--- a/docs/src/manual.md
+++ b/docs/src/manual.md
@@ -13,5 +13,10 @@ BSplineApprox
 CubicHermiteSpline
 PCHIPInterpolation
 QuinticHermiteSpline
-looks_linear
 ```
+
+# Utility Functions
+
+```@docs
+DataInterpolations.looks_linear
+```
\ No newline at end of file
diff --git a/src/interpolation_utils.jl b/src/interpolation_utils.jl
index 3ee0cc5a..6150fe04 100644
--- a/src/interpolation_utils.jl
+++ b/src/interpolation_utils.jl
@@ -98,9 +98,11 @@ seems_linear(assume_linear_t::Bool, _) = assume_linear_t
 seems_linear(assume_linear_t::Number, t) = looks_linear(t; threshold = assume_linear_t)
 
 """
-Determine if the abscissae `t` are regularly distributed, taking the standard deviation of \
-the difference between the array of abscissae with respect to the straight line linking \
-its first and last elements, normalized by the range of `t`.
+Determine if the abscissae `t` are regularly distributed, taking the standard deviation of
+the difference between the array of abscissae with respect to the straight line linking
+its first and last elements, normalized by the range of `t`. If this standard deviation is
+below the given `threshold` (default value: 1e-2), the vector looks linear. Internal
+function or structure - interface may change.
 """
 function looks_linear(t; threshold = 1e-2)
     length(t) <= 2 && return true

From 9e9876a78c2592181189b021fa53f4c1f158a8c4 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Daniel=20Gonz=C3=A1lez?= <daniel.glez.arribas@gmail.com>
Date: Wed, 31 Jul 2024 15:44:31 +0200
Subject: [PATCH 10/11] Formatting

---
 src/interpolation_utils.jl  | 6 ++++--
 test/interpolation_tests.jl | 2 +-
 2 files changed, 5 insertions(+), 3 deletions(-)

diff --git a/src/interpolation_utils.jl b/src/interpolation_utils.jl
index 6150fe04..ccd3db63 100644
--- a/src/interpolation_utils.jl
+++ b/src/interpolation_utils.jl
@@ -98,11 +98,13 @@ seems_linear(assume_linear_t::Bool, _) = assume_linear_t
 seems_linear(assume_linear_t::Number, t) = looks_linear(t; threshold = assume_linear_t)
 
 """
+    looks_linear(t; threshold = 1e-2)
+
 Determine if the abscissae `t` are regularly distributed, taking the standard deviation of
 the difference between the array of abscissae with respect to the straight line linking
 its first and last elements, normalized by the range of `t`. If this standard deviation is
-below the given `threshold` (default value: 1e-2), the vector looks linear. Internal
-function or structure - interface may change.
+below the given `threshold`, the vector looks linear (return true). Internal function -
+interface may change.
 """
 function looks_linear(t; threshold = 1e-2)
     length(t) <= 2 && return true
diff --git a/test/interpolation_tests.jl b/test/interpolation_tests.jl
index 601047f6..1d4e2e6b 100644
--- a/test/interpolation_tests.jl
+++ b/test/interpolation_tests.jl
@@ -817,4 +817,4 @@ f_cubic_spline = c -> square(CubicSpline, c)
         b_fallback = @benchmark $A_fallback(rand()) samples=n_samples
         @test mean(b_linear.times) < mean(b_fallback.times)
     end
-end
\ No newline at end of file
+end

From 72008673a2e9bb7ad7b1a230224b9c084a4fbee4 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Daniel=20Gonz=C3=A1lez?= <daniel.glez.arribas@gmail.com>
Date: Wed, 31 Jul 2024 16:03:09 +0200
Subject: [PATCH 11/11] Add newline at end of file

---
 docs/src/manual.md | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/docs/src/manual.md b/docs/src/manual.md
index 08e4c32f..3818a9e3 100644
--- a/docs/src/manual.md
+++ b/docs/src/manual.md
@@ -19,4 +19,4 @@ QuinticHermiteSpline
 
 ```@docs
 DataInterpolations.looks_linear
-```
\ No newline at end of file
+```