Plots weird problem plotting a piecewise function at discontinuity points

using Plots;
g_k =0.0:0.01:10.0;
A1=1.0; A2=2.0; α=0.33;
κ2 =3.0;
κ̂2 = κ2 / (1.0 - ((A1/A2)^(1/α))); #3.41842
function f3p(k)
    if (k < κ̂2)
        return α*A1*(k^(α - 1.0))
    elseif (k > κ̂2)
        return α*A2*(( max(0.0, k - κ2) )^(α- 1.0))
plot(g_k, f3p)


The weird looking vertical line at the discontinuity at 3.41842 should not be there.
How can I remove it?

My only hack is to do it manually, but that’s not a good solution for more complicated problems

plot!(0.0:0.01:κ̂2, f3p, lab="f'(k)", c=1, lw=1)
plot!(κ̂2:0.01:10, f3p, lab="f'(k)", c=1, lw=1)

Well I don’t think there is an easy solution for exactly what you propose. plot essentially just takes 2 arrays of x and y coordinates and draws lines between them. How should it know which line to omit?

Some other ideas:

  1. Instead of using two plot calls you can mark the jump with the value NaN. IIRC NaNs break the line. So your picewise function could return NaN exactly at the jump (or a small neighbourhood) and if your are plotting enoigh points such that one of them is inside that window, then it should work out.
  2. Instead of plot/plot! you could use scatter/scatter! which just draws dots at the given location. If you play a bit with markersize and number of points you probably can get something reasonable as well.
  3. Maybe there is a package for piecewise functions that can help you?

Scatter() looks way too weird.
Here is a more minimal/parsimonious MWE

using Plots;
grid= 0.0:0.0001:1.0
g(x)= (x<0.5) ? (x) : (x-0.5)
plot(grid, g)
scatter!(grid, g)

I’m hoping for some kind of options to ignore discontinuities.

Yes you’ll need to play with some setting to make scatter look better. Did you try e.g. with color=:blue, markersize=0.25 or something like that? I am on mobile and can’t play with these settings myself.

Edit: You probably want to disable the border (black line) around the markers. You can try markerstrokewidth=0 for that.

It would probably make sense to use the information about the intervals that you have, by incorporating them into a type for which you can make a recipe which then makes gaps at the correct places. If you just have a blackbox function outputting numbers, you cannot really know whether two adjacent points are separated by a singularity or not, could always be continuous on a very small scale I guess?

As @jules mention, since this is just some blackbox you would need to provide it yourself with your own recipe, or break it into NaNs. scatter is probably okay too. Here’s one other approach you could take, identifying singularities by large slopes:

using Plots
function fix_singularities(y, x; slope_tol=100)
    _x = similar(x, 0)
    _y = similar(y, 0)
    for i in firstindex(y):(lastindex(y)-1)
        yᵢ, yᵢ₊₁ = y[i], y[i+1]
        xᵢ, xᵢ₊₁ = x[i], x[i+1]
        slope = (yᵢ₊₁ - yᵢ) / (xᵢ₊₁ - xᵢ)
        if abs(slope) < slope_tol
            push!(_x, x[i])
            push!(_y, y[i])
            push!(_x, x[i], (x[i] + x[i+1]) / 2)
            push!(_y, y[i], NaN)
    return _y, _x
g(x) = (x < 0.5) ? (x) : (x - 0.5)
x = 0.0:0.0001:1.0
y = g.(x)
_y, _x = fix_singularities(y, x)
plot(_x, _y) 

Maybe this is robust enough for your use case, if you are willing to recompute the x and y like this (you can of course define your own method for providing a function and a grid still).

1 Like

I was wrong to call it a singularity.
I meant discontinuity.

And I STRONGLY believe a good plotting package should have a fixdiscontinuity option.
Just b/c others don’t do it, doesn’t mean we shouldn’t try.

Random thought: Maybe Autodiff could identify true discontinuities? Just computing a local gradient does of course not contain this information. I am not too familiar with how all of these different packages work, but maybe one of the approaches could detect when there is a switch to a different code branch or just know due to code analysis when a jump occurs.

1 Like

Mathematica plots discontinuous functions correctly

Yes but this quite different from plotting values from a blackbox. In principle you could define a type in Julia that representa a piecewise function and define a plotting recipe that handles the discontinuity.

1 Like

I had no idea my original function was discontinuous or where the discontinuity points were

Some simple approach based on the derivative and on dense sampling (try also the tan() case):

using Statistics, Plots; gr(dpi=600)

A1 = 1.0; A2 = 2.0; α = 0.33
κ2 = 3.0
κ̂2 = κ2 / (1.0 - ((A1/A2)^(1/α)));  # 3.41842

f3p(k) = (k < κ̂2) ?  α*A1*(k^(α - 1.0)) : α*A2*((max(0.0, k - κ2))^(α- 1.0))
# f3p(k) = tan(k)               # try this too!

N = 1000                        # the higher/denser the sampling the better
k = range(0, 10, N)
g_k = f3p.(k)

δk = [diff(g_k); 0]
dp = abs.(δk[2:end-1])
ϵ = N * median(dp)              # discontinuity threshold
ix = [false; dp .> ϵ; false]
g_k[ix] .= NaN

plot(k, g_k)

The said functions in Mathematica is not the same thing as in Plots.jl. Mathematica knows it is defined piecewise (because of its symbolic nature) and does not assume continuity whereas, in Plots.jl, it is always numeric and it almost never knows it’s symbolic structure unless a recipe with a custom type defines piecewise function structure. Hence, like Matlab and several other plotting libraries, discontinuities in the lines of plot are inferred by NaN values.