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

**URL:** <https://discourse.julialang.org/t/hi-i-have-a-function-lets-say-f-sin-0-10-now-i-want-to-break/54826>\
**Category:** General Usage\
**Created:** [February 8, 2021, 6:56am UTC](https://discourse.julialang.org/t/hi-i-have-a-function-lets-say-f-sin-0-10-now-i-want-to-break/54826 "2021-02-08T06:56:40Z")\
**Posts on this page:** 9\
**Page:** 2

<div class="post-metadata">

**Author:** ![PeterSimon](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/petersimon/32/25193_2.png) [@PeterSimon](https://discourse.julialang.org/u/PeterSimon)\
**Post date:** [February 8, 2021, 10:25pm UTC](https://discourse.julialang.org/t/hi-i-have-a-function-lets-say-f-sin-0-10-now-i-want-to-break/54826/21 "2021-02-08T22:25:54Z")

</div>

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.

```plaintext
"""
    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:

```plaintext
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.")

```

![piecewise](https://global.discourse-cdn.com/julialang/original/3X/6/f/6f981a377563659e0160e216071a45940929e6a1.png)

---

<div class="post-metadata">

**Author:** ![rafael.guerra](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rafael.guerra/32/216610_2.png) [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)\
**Post date:** [February 8, 2021, 11:16pm UTC](https://discourse.julialang.org/t/hi-i-have-a-function-lets-say-f-sin-0-10-now-i-want-to-break/54826/22 "2021-02-08T23:16:03Z")

</div>

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

```julia
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")

```

![Dierckx_linear_spline_free_knots2](https://global.discourse-cdn.com/julialang/original/3X/f/3/f3614dc06cab57aa91a59cfc7f408bf32471e100.png)

---

<div class="post-metadata">

**Author:** ![PeterSimon](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/petersimon/32/25193_2.png) [@PeterSimon](https://discourse.julialang.org/u/PeterSimon)\
**Post date:** [February 9, 2021, 1:34am UTC](https://discourse.julialang.org/t/hi-i-have-a-function-lets-say-f-sin-0-10-now-i-want-to-break/54826/23 "2021-02-09T01:34:13Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![maxkapur](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/maxkapur/32/21208_2.png) [@maxkapur](https://discourse.julialang.org/u/maxkapur)\
**Post date:** [February 9, 2021, 2:10am UTC](https://discourse.julialang.org/t/hi-i-have-a-function-lets-say-f-sin-0-10-now-i-want-to-break/54826/24 "2021-02-09T02:10:08Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![rafael.guerra](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rafael.guerra/32/216610_2.png) [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)\
**Post date:** [February 9, 2021, 6:12am UTC](https://discourse.julialang.org/t/hi-i-have-a-function-lets-say-f-sin-0-10-now-i-want-to-break/54826/25 "2021-02-09T06:12:28Z")

</div>

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

```julia
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

```

---

<div class="post-metadata">

**Author:** ![Karajan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/karajan/32/8545_2.png) [@Karajan](https://discourse.julialang.org/u/Karajan)\
**Post date:** [February 9, 2021, 12:59pm UTC](https://discourse.julialang.org/t/hi-i-have-a-function-lets-say-f-sin-0-10-now-i-want-to-break/54826/26 "2021-02-09T12:59:44Z")

</div>

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:

```julia
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.

---

<div class="post-metadata">

**Author:** ![Iulian.Cioarca](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/iulian.cioarca/32/30166_2.png) [@Iulian.Cioarca](https://discourse.julialang.org/u/Iulian.Cioarca)\
**Post date:** [February 9, 2021, 2:04pm UTC](https://discourse.julialang.org/t/hi-i-have-a-function-lets-say-f-sin-0-10-now-i-want-to-break/54826/27 "2021-02-09T14:04:06Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![rafael.guerra](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rafael.guerra/32/216610_2.png) [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)\
**Post date:** [February 9, 2021, 2:50pm UTC](https://discourse.julialang.org/t/hi-i-have-a-function-lets-say-f-sin-0-10-now-i-want-to-break/54826/28 "2021-02-09T14:50:40Z")

</div>

Just for the record, the suggested use of [DifferentialEquations.jl](https://github.com/SciML/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.

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

```

![DiffEqs_sampling](https://global.discourse-cdn.com/julialang/original/3X/1/c/1c26e343b456a26a8be443c3a97c0b4ac6b02ee0.png)

---

<div class="post-metadata">

**Author:** ![rafael.guerra](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rafael.guerra/32/216610_2.png) [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)\
**Post date:** [February 10, 2021, 8:34am UTC](https://discourse.julialang.org/t/hi-i-have-a-function-lets-say-f-sin-0-10-now-i-want-to-break/54826/29 "2021-02-10T08:34:45Z")

</div>

> [@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.

[Previous page](https://discourse.julialang.org/t/hi-i-have-a-function-lets-say-f-sin-0-10-now-i-want-to-break/54826.md?page=1)