Hi! I have a function, let's say `f(ϕ) = sin(ϕ), ϕ = 0...10`. Now I want to break

[PeterSimon]

Here’s another approach that also lets you specify the maximum absolute error between the piecewise linear interpolant and the original function. It uses the amazing ApproxFun package.

"""
    piecewise(fun, a, b, maxerr)

Return a vector of points on the interval `[a,b]` at which a piecewise linear interpolant
to the function `fun` will yield a maximum absolute error at most `maxerr`, using the `ApproxFun` package.
"""
function piecewise(fun, a, b, maxerr)
    dom = Domain(a..b)
    f = Fun(fun, dom)
    ya, yb = fun(a), fun(b)
    m = (yb-ya)/(b-a)
    flinear = Fun(x -> m*(x-a) + ya, dom) 
    extr = extrema(flinear-f)
    err = maximum(abs.(extr))
    if err ≤ maxerr
        [a,b]
    else
        l = piecewise(fun, a, (a+b)/2, maxerr)
        r = piecewise(fun, (a+b)/2, b, maxerr)
        vcat(l, r[2:end]) # eliminate duplicate midpoint
    end
end

Here is an example of its use:

using ApproxFun
using Plots
f(x) = x*sin(1/x)
fun = f
a, b = 0.03, 0.2
maxerr = 0.002
knots =  piecewise(fun, a, b, maxerr)
fknots = fun.(knots)
plot(range(a, b, length=500), fun, label = "Original Function")
scatter!(knots, fknots, label = "Interp. Knots")
plot!(knots, fknots, label = "Linear Interp.")

[rafael.guerra]

Fyi, a similar result with Dierckx linear splines with free knots:

using Dierckx, Plots; gr()
x = LinRange(0.03, 0.2, 200)
y = x .* sin.(1 ./ x)
spl1 = Spline1D(x, y; k=1, bc="nearest", s=3e-5)
xk = get_knots(spl1)
plot(x, y, label="x*sin(1/x)", legend =:bottom)
scatter!(xk, spl1(xk), ms=3, label="Dierckx linear spline free knots")

[PeterSimon]

Sorry, I didn’t read your previous post carefully enough to see that you are also specifying the maximum allowable error. I’ve edited my post to indicate that it is just another possible way of doing the same thing.

[maxkapur]

Just a mathematical idea here, but why insist that the knots lie on the curve itself? All you want to do is minimize the average error from the function. (And do you know which error function you are trying to minimize? E.g. the absolute distance, the squared distance, or something else?)

Moreover, if the function you are approximating is a sine wave, then choosing points on the original curve ensures that all of your interpolated values are underestimates (in absolute value) of the function. Depending on the application, this could be an undesirable form of bias, and you may prefer a model where there are some overestimates to “cancel out” the underestimates.

Suppose we are using the easy model where the knots are evenly spaced in x. You could start by defining your knots as [(x, f(x)) for x in 0:step:10], then “wiggle” the y coordinates locally (say, allowing them to vary by ±10%) in order to minimize total error. Notice that for the sine function, slightly increasing the y-values in the neighborhood of x = \pi/2 from the true value of \sin(y) reduces the average distance between the curve and the linear spline at all points.

Combining this approach with smarter selection of the x coordinates might yield a more robust model.

[rafael.guerra]

@maxkapur, please note that the Dierckx linear spline free knots do not lie on the curve.

julia> spl1(xk) .- f.(xk)
67-element Vector{Float64}:
  0.0003490975575791283
 -0.0006339828718615326
 -2.1984047666268134e-5
  0.00019359272018671866
  0.00020392286807298377
  ⋮
 -0.00026440715875873655
 -0.00025290289025928225
 -0.00026277647778455426
 -0.0002755482682172683
 -0.0002023539394988283

[Karajan]

Alright, I had some time to try out all the responses!

@PeterSimon

Here’s another approach that also lets you specify the maximum absolute error between the piecewise linear interpolant and the original function. It uses the amazing ApproxFun package.

This works quite well and a plus point is that I know what the algorithm does and what maxerr means – with the other options I’d have to dig a little. This one works on the function directly and doesn’t require prior “discretization”.

@rafael.guerra

Dierckx linear splines in its free knots option, produces good results by tuning the admissible error tolerance parameter s .

On the other hand I noticed that with the Dierckx option you need a list of points but the points don’t have to be equally spaced – that might come in handy. (I have thought about going at the problem the other way around: take very small steps first and then recombine elements according to some criterion.)

@Iulian.Cioarca

Just take a look at the demo in the example folder. If you have issues please ping me since I’m the author.

This also works well! It seems you have a nice little package there, have you thought about registering it?
For the record the code:

using AdaptiveSampling
x = 0:0.02:2π
y = sin.(4x) .+ 2
thr = 0.01
tc, yc, idx2save = AdaptiveSampling.anchor_point(x, y, thr)
scatter(x, y, label = "Original Function")
scatter!(tc, yc, ms = 2, label = "Interp. Knots")

A interesting thing here is that it only seems to pick existing points of the list of x and doesn’t create new ones. This might also be helpful if I want to restrict the result in this way.

I’ll have to try around a little, see how I can bend my application to these different approaches and which ones will work best in the end.

[Iulian.Cioarca]

Yes, it picks existing points from the original array. The original array doesn’t need to be sampled with a constant sampling rate. One drawback is that if you have a noisy signal some of the “stored” points might end up on the noise peaks and dips, but the condition that a local region variation is smaller than thr is always satisfied.( I usually add a wavelet noise reduction step before compression).
To be honest I didn’t think anyone would use the package. I only use it for lossy compression of large oscilloscope waveforms. But I will consider registering it:)

[rafael.guerra]

Just for the record, the suggested use of DifferentialEquations.jl for this problem does not seem to be a good idea. The variable time-step logic is different, as illustrated in simple example below for one of the solvers.

using DifferentialEquations, Plots
dfdt(u,p,t)  = cos(t)   # define the system 
u0 = 0; tspan = (0.0,π)  # initial conditions & timespan
prob = ODEProblem(dfdt, u0, tspan)  # define the problem
sol0 = solve(prob,Tsit5(),reltol=1e-8,abstol=1e-12)  # solve it
plot(sol0, label="f(t) = sin(t) from DiffEqs")
scatter!(sol0.t, sol0.u, label="@DiffEqs time steps", legend=:bottom)

[rafael.guerra]

@Karajan

(I have thought about going at the problem the other way around: take very small steps first and then recombine elements according to some criterion.)

This is exactly what is done in the linear splines with free knots code above, where small steps are taken first and then Dierckx.jl estimates the minimum number of knots along x-axis that are required to meet the accuracy that was set.