From 4244ee744413efa4cb6be2ae6f32dd90580c00b8 Mon Sep 17 00:00:00 2001
From: Mauro Werder <mauro3@runbox.com>
Date: Fri, 22 Dec 2023 14:47:45 +0100
Subject: [PATCH] Added derivative output to piecewise_linear

---
 src/misc.jl | 34 +++++++++++++++++++++-------------
 1 file changed, 21 insertions(+), 13 deletions(-)

diff --git a/src/misc.jl b/src/misc.jl
index 5fae37c..c4bc1c7 100644
--- a/src/misc.jl
+++ b/src/misc.jl
@@ -235,11 +235,13 @@ end
 
 
 """
-    piecewiselinear(xs::AbstractVector, ys::AbstractVector)
+    piecewiselinear(xs::AbstractVector, ys::AbstractVector; derivative=false)
 
 Return a piecewise linear function `f(x) -> y`` through points (xs,ys).
 Extrapolation with x<xs[1]->ys[1] and for x>xs[end] -> ys[end] .
 
+If `derivative=true` it returns the derivative (slope) function instead.
+
 Notes:
 - xs needs monotonic, i.e. to be sorted in ascending or decending order
 - when length(xs)==1, then it returns a constant function.
@@ -254,26 +256,32 @@ function piecewiselinear(xs::AbstractVector, ys::AbstractVector)
     end
     rats = diff(ys)./diff(xs)
     if issorted(xs)
-        return function (xq)
+        return function (xq; derivative=false)
             # xq<xs[1] && error("cannot extrapolate")
             # xq>xs[end] && error("cannot extrapolate")
             # xq==xs[end] && return ys[end]
-            xq<=xs[1] && return ys[1]
-            xq>=xs[end] && return ys[end]
+            xq<xs[1] && return derivative ? ys[1]*0 : ys[1]
+            xq>=xs[end] && return derivative ? ys[end]*0 : ys[end]
             ii = searchsortedlast(xs, xq)
-            out = ys[ii] - (xs[ii]-xq)*rats[ii]
-            return out
+            if derivative
+                return rats[ii]
+            else
+                return ys[ii] - (xs[ii]-xq)*rats[ii]
+            end
         end
     elseif issorted(xs, Base.Order.Reverse)
-        return function (xq)
+        return function (xq; derivative=false)
             # xq>xs[1] && error("cannot extrapolate")
             # xq<xs[end] && error("cannot extrapolate")
             # xq==xs[end] && return ys[end]
-            xq>=xs[1] && return ys[1]
-            xq<=xs[end] && return ys[end]
+            xq>xs[1] && return derivative ? ys[1]*0 : ys[1]
+            xq<=xs[end] && return derivative ? ys[end]*0 : ys[end]
             ii = searchsortedlast(xs, xq, Base.Order.Reverse)
-            out = ys[ii] - (xs[ii]-xq)*rats[ii]
-            return out
+            if derivative
+                return rats[ii]
+            else
+                return ys[ii] - (xs[ii]-xq)*rats[ii]
+            end
         end
     else
         error("xs must be sorted")
@@ -300,7 +308,7 @@ julia> pf, span = parameterized_curve(x,y);
 julia> xy = hcat(pf.(r)...); plot(xy[1,:], xy[2,:])
 ````
 """
-function parameterized_curve(x,y)
+function parameterized_curve(x,y; derivative=false)
     # TODO: make dimension agnostic sometime
     pathlength = zeros(typeof(x[1]), length(x)-1) # assume all same datatype
     for c in (x,y)
@@ -309,5 +317,5 @@ function parameterized_curve(x,y)
     pathlength = [0, cumsum(sqrt.(pathlength))...]
     fx = piecewiselinear(pathlength, x)
     fy = piecewiselinear(pathlength, y)
-    return pl -> [fx(pl), fy(pl)], (0, pathlength[end])
+    return pl -> [fx(pl; derivative), fy(pl; derivative)], (0, pathlength[end])
 end