Find intersection between curves in julia

Hello,
I have two curves with the values:

col_y = [8.02511788606e7, 9.98481296252e7, 2.328332483329e8, 6.148181496526e8, 9.453171952624e8, 1.2331292887706e9, 1.4394296796998e9, 1.460066157955e9, 1.3570403824485e9, 1.4269447174002e9]

vin_y = [1.434628836798e8, 1.318375623746e8, 1.108524843855e8, 6.91176196241e7, 4.19137382458e7, 1.48676444645e7, 1.00659726053e7, 7.8734512643e6]

Is it possible to find the intersection x between these two curves? what is the approach in julia?
Thank you

Do you have a model for what kind of curves these are (they look sigmoid, do you know what process generated them)? And is the x-scale the same or do they for instance span between the same limits (which is not the same thing since the number of elements is different)? Linear x-scale?

The reason you need to know is the answer will vary wildly depending on which way you interpret the question. A vector of values isn’t a function, in fact there is an uncountable infinity of functions which fit exactly.

The way I would do this is

  • Find an analytic function f_col, f_vin that fits each of the sample vectors (for instance using LsqFit.jl)
  • Form the function g(x) = f_col(x) - f_vin(x)
  • Find the zero(s) of g(x) (for instance using Roots.jl)
1 Like

You have points, not curves. The names suggest they are y-values, but it’s not clear what the corresponding x-values are for each.

Plotting these by assuming x=1:10 for the first and x=1:8 for the second gives:

image

So there is an intersection just after x=2, but this is based on a linear interpolation between points which Plots does automatically. As Gustaphe says you might have some other function in mind to interpolate between the points, which will result in a different point of intersection.

1 Like

For this problem, Dierckx is just fine:

using Dierckx
n = length(vin_y)
spl = Spline1D(1:n, col_y[1:n] - vin_y)
x0 = roots(spl)   # 2.4181
y0 = Spline1D(1:n, vin_y)(x0)  # 1.2583e8

2 Likes

Thanks for the answers. I made a model, but it does not exactly fit the data:

Thus, I wanted to check the time difference of intersection between the actual data and the model…
I’ll check the suggestions.

Here’s my way. I have an experimental version of LsqFit installed, so there will be a slight difference in syntax, but I just realized I’m in a bit of a hurry so I won’t fix that.

using Plots, LsqFit, Roots
default(;fontfamily="Computer Modern", label=nothing, seriestype=:scatter) # Better plot defaults
x_1 = 1:10
x_2 = 1:8
y_1 = [8.02511788606e7, 9.98481296252e7, 2.328332483329e8, 6.148181496526e8, 9.453171952624e8, 1.2331292887706e9, 1.4394296796998e9, 1.460066157955e9, 1.3570403824485e9, 1.4269447174002e9]
y_2 = [1.434628836798e8, 1.318375623746e8, 1.108524843855e8, 6.91176196241e7, 4.19137382458e7, 1.48676444645e7, 1.00659726053e7, 7.8734512643e6]
model(x, p) = 1e9*p[1]/(1+exp(-(x-p[2])/p[3])) # logistic function
f_1 = curve_fit(model, x_1, y_1, [1.0, 1.0, 1.0])
f_2 = curve_fit(model, x_2, y_2, [1e-1, 1.0, -1.0])
x_long = range(0.5,10.5;length=200)
pl_1 = plot(x_1, y_1)
plot!(pl_1, x_2, y_2)
plot!(pl_1, x_long, model.(x_long,Ref(f_1.param)); st=:line)
plot!(pl_1, x_long, model.(x_long,Ref(f_2.param)); st=:line)
g(x) = model(x, f_1.param) - model(x, f_2.param)
pl_2 = plot(g, x_long, st=:line)
plot(pl_1, pl_2; layout=(2,1)
find_zero(g, 2.0) # returns ~ 2.27

image

If you want the intersection between the actual data series, the answer is there is none. None of your sample points measure the same amount of bacteria of the two species.

1 Like

Thank you!
the Dierckx solution worked out of the box. By using the actual x values (sorry, I forgot to add them originally)

x = [0.0, 6.465599999999999, 12.544800000000002, 24.0, 29.107200000000002,  45.772800000000004, 53.392799999999994, 68.0616]

I got the exact intersection:

The more canonical solution by gustaphe is giving me troubles. The definition of the model is very similar to the model I have set:

model = true_logistic(x, p) = p[1] ./ (1 .+ exp.(-p[2].*(x.-p[3])))

that is a logistic regression but I set 1 instead of 1e9. When I define g and find_zero:

julia> g(x) = model(col_x[1:8], P_col) - model(col_x[1:8], P_vin)
g (generic function with 1 method)

julia> find_zero(g, col_x[2])
ERROR: Roots.ConvergenceFailed("Stopped")
Stacktrace:
 [1] find_zero(fs::Function, x0::Float64, M::Secant, N::AlefeldPotraShi; p::Nothing, tracks::Roots.NullTracks, verbose::Bool, kwargs::Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
   @ Roots ~/.julia/packages/Roots/H7pXT/src/order0.jl:57
 [2] #find_zero#52
   @ ~/.julia/packages/Roots/H7pXT/src/order0.jl:37 [inlined]
 [3] find_zero
   @ ~/.julia/packages/Roots/H7pXT/src/order0.jl:35 [inlined]
 [4] #find_zero#16
   @ ~/.julia/packages/Roots/H7pXT/src/find_zero.jl:704 [inlined]
 [5] find_zero(f::Function, x0::Float64)
   @ Roots ~/.julia/packages/Roots/H7pXT/src/find_zero.jl:704
 [6] top-level scope
   @ none:1

julia> find_zero(g, 2.0)
ERROR: Roots.ConvergenceFailed("Stopped")
Stacktrace:
 [1] find_zero(fs::Function, x0::Float64, M::Secant, N::AlefeldPotraShi; p::Nothing, tracks::Roots.NullTracks, verbose::Bool, kwargs::Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
   @ Roots ~/.julia/packages/Roots/H7pXT/src/order0.jl:57
 [2] #find_zero#52
   @ ~/.julia/packages/Roots/H7pXT/src/order0.jl:37 [inlined]
 [3] find_zero
   @ ~/.julia/packages/Roots/H7pXT/src/order0.jl:35 [inlined]
 [4] #find_zero#16
   @ ~/.julia/packages/Roots/H7pXT/src/find_zero.jl:704 [inlined]
 [5] find_zero(f::Function, x0::Float64)
   @ Roots ~/.julia/packages/Roots/H7pXT/src/find_zero.jl:704
 [6] top-level scope
   @ none:1

What am I getting wrong?

Likely you have a bad starting point. With only 1 intersection you should be able to use a bracketing method: find_zero(g, 1, 50)

The reason I multiplied by 1e9 is that LsqFit (and any fitting method) is better at finding similarly sized parameters, so I scaled up the model so we can expect p.≈1.0.

Your problem is in the definition of g. It needs to depend on x, not col_x.

1 Like

Ah, so simple. It worked!

julia> g(x) = model(x, P_col) - model(x, P_vin)
g (generic function with 1 method)

julia> find_zero(g, col_x[2])
10.018996507513226

Thank you

Glad it worked out!

Can I ask what function you’re fitting to in the figure? Because the logistic function seems like a better fit, and from the little I know of microbiology it’s not an unreasonable model, so it seems a little strange to only use the better fit to verify the worse.

Hi, it is the true_logistic; I see that does not fit the data completely but it is the same model that I used for another group of bacteria (where it was more fitting). For consistency, I am sticking to this…

Ah okay, I missed the log axis. The fit is decent :slight_smile: