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

Hi! I have a function, let’s say f(ϕ) = sin(ϕ), ϕ = 0...10. Now I want to break that function down in a list of lines that approximates the real function well (maybe with a maximum or minimum line length, angle between lines, or number of lines). Is there a package that can help with that? Bonus points if it also works in polar coordinates.

Note that the original poster on Slack cannot see your response here on Discourse. Consider transcribing the appropriate answer back to Slack, or pinging the poster here on Discourse so they can follow this thread.
(Original message :slack:) (More Info)

A little more information based on the discussion on Slack:

I have a function that only works on lines (let’s say (0.0, 0.0) -> (0.1, 0.3)), that’s why I want to turn my sin (for example) into a piecewise linear function/list of lines.

I know that some combination of inputs (e.g. number of lines and maximum length of lines) might give trouble, we can ignore that for now.

I imagined that adaptive algorithms of integration/differentiation equations might do something close to what I’m looking for, but 1) I’m not sure about that being true, 2) I don’t know any package that does just that.

PS: Is there any way for me to change the title of the post? At the moment it’s quite … nondescript.

1 Like

Are you looking for something like this?

julia> function lines(f, r)
           step = Float64(r.step)
           map(x -> ((x, f(x)), (x+step, f(x+step))), r[1:end-1])
       end
lines (generic function with 2 methods)

julia> lines(sin, 0.0:0.1:1.0)
10-element Vector{Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}:
 ((0.0, 0.0), (0.1, 0.09983341664682815))
 ((0.1, 0.09983341664682815), (0.2, 0.19866933079506122))
 ((0.2, 0.19866933079506122), (0.30000000000000004, 0.2955202066613396))
 ((0.3, 0.29552020666133955), (0.4, 0.3894183423086505))
 ((0.4, 0.3894183423086505), (0.5, 0.479425538604203))
 ((0.5, 0.479425538604203), (0.6, 0.5646424733950354))
 ((0.6, 0.5646424733950354), (0.7, 0.644217687237691))
 ((0.7, 0.644217687237691), (0.7999999999999999, 0.7173560908995227))
 ((0.8, 0.7173560908995228), (0.9, 0.7833269096274834))
 ((0.9, 0.7833269096274834), (1.0, 0.8414709848078965))
1 Like

Thank you but this isn’t what I’m looking for.
If we stay with the example of sin, around 0 this would be a good approximation but around π/2 it might be to coarse. So I’m looking for an package/algorithm that picks the steps so that I can set “maximum or minimum line length, angle between lines, or number of lines” (or at least more than just number of lines). If I have the points, turning them into lines should be fine.

I see. In that case I can’t help you, unfortunately.

If you write the derivative f’(t) of your function and plug it in the DifferentialEquations.jl package, you can define an error tolerance and you will get as solution f(t), sampled at irregular times. Use the resulting array of irregular times to break down your original function.
PS: Yes, it sounds like using a sledgehammer to crack a nut.

2 Likes

Well, a sledgehammer indeed, but I prefer that over using my fingernails so I’ll give it a try if I don’t find anything better fitting. Thanks!

2 Likes

Another approach might be piecewise linear interpolation. I have never tried the PiecewiseLinearOpt.jl package but have a look at it, also the associated presentation by the author.

1 Like

It sounds like geometrically, you want to define a small circle around a point and then find the point where the function intersects the circle, then draw a circle around that point and find the next intersection and so on.

Fundamentally, this is a root finding problem in two dimensions. It should be possible to use Roots.jl or nlsolve, for example.

3 Likes

Hm, I’m not sure if I misunderstand the package but to me it sounds like optimization of a function that already is piecewise?

I guess that comes close to what I’m looking for but isn’t really that “adaptive”?

I’ve coded something crude up:

using Optim

f(x) = sin(x)
x = 0.0:0.01:3

function pwl(x, pnts)
    idx = length(pnts)
    for i in 2:length(pnts)
        if x <= pnts[i]
            idx = i
            break
        end
    end
    p1 = pnts[idx-1]
    p2 = pnts[idx]
    f1 = f(p1)
    f2 = f(p2)
    return f1 + (f2 - f1) * (x - p1) / (p2 - p1)
end

function cost(p)
    pnts = copy(p)
    push!(pnts, first(x))
    push!(pnts, last(x))
    sort!(pnts)

    pwlf = x -> pwl(x, pnts)
    
    sum(abs2, f.(x) .- pwlf.(x))
end

N = 19
x0 = range(0, 3, length = N)[2:end-1] |> collect
lower = zeros(N-2)
upper = ones(N-2) .* 3
result = optimize(cost, lower, upper, x0)

In the linear part of the sin there are very few points, around π/2 there are a lot.

Of course this won’t work in a number of cases (most of them, probably). This also only considers the number of points/lines, not their length or angle.

I guess with enough time I would be able to code something properly but of course I’d prefer to reuse something existing.

Isn’t that exactly what you want? Near the zeros (0, pi, 2pi, etc) the sine can be approximate for a large interval with a linear function. so you would expect one long piece there

Yes it is. I meant that the circle approach doesn’t seem to do this so I’m not sure it’s going to work for me.

no problem, I’m glad you found a solution that works for you

@Karajan, you are probably right. I thought that package would allow optimizing an initial set of N-breakpoints of the piecewise linear function. May be someone can confirm.

Well, I’m pretty certain this would break down with more complicated functions or more points so I’m not really happy with the Optim approach. But if there doesn’t exist anything else I’ll try my luck :smiley:

Dierckx linear splines in its free knots option, produces good results by tuning the admissible error tolerance parameter s. You can use the array xk with the knots along the x-axis.

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

Dierckx_linear_spline_free_knots

3 Likes

This look good, great find. I’ll have a play with it, thank you so much!

1 Like

Does this help?

1 Like

Looks also interesting! However, there don’t seem to be any docs so I’ll need to work through it to see if helps.

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