How to find inflection point of curve?

Luigi_Marongiu · November 28, 2020, 9:06am

Hello,
I have the function

y = ϑ - ((1/ω) * log(x/ε)) - 1/μ

Is there a simple way in Julia to find the flection point of the function?
if I set these paramaters:

mu    = 0.72        # maximum growth rate resistant (E. coli)
omega = 0.05        # outflow
tmax  = 4000.0      # time span 0-tmax
epsilon = 1           # concentration equivalent to 1 phage
Tp = 2.392731385751533
ε, ω, ϑ, μ = [epsilon, omega, Tp, mu, tmax]
Y = []
    X = 0.0:0.1:τ
    for x in X
        y = ϑ - ((1/ω) * log(x/ε)) - 1/μ
        push!(Y, y)
    end

I get the figure:

Is it possible to find the inflection point of the curve, which should be at around t=500?
Or it is all mathematics and I need to solve the point by mathematical analysis?
Thank you

dlfivefifty · November 28, 2020, 9:19am

ApproxFun.jl Can do this

rafael.guerra · November 28, 2020, 9:20am

You may need to find the point with minimum radius of curvature.

jonathanBieler · November 28, 2020, 9:49am

You can use automatic differentiation :

using ForwardDiff, Roots

f(x) = sin(x)
g = x->ForwardDiff.derivative(f,x)
find_zero(x->ForwardDiff.derivative(g,x), 3)

(you still need to check that the second derivative changes signs around the zero)

That said the second derivative of your function is 1/(ω*x^2) and thus has not inflection point.

You can also use broadcasting (or an array comprehension) instead of a for loop to evaluate your function :

f(x, ϑ, ω, ε, μ) = ϑ - ((1/ω) * log(x/ε)) - 1/μ
x = 0.0:0.1:τ
y = f.(x, ϑ, ω, ε, μ)

rafael.guerra · November 28, 2020, 11:22am

The maximum curvature (or minimum radius of curvature) occurs at time t = 14.14:

using ForwardDiff, Roots, Plots
function f(t)
    ε, ω, ϑ, μ = [1., 0.05 , 2.392731385751533,  0.72 ]
    return ϑ - (1/ω) * log(t/ε) - 1/μ
end
τ = 4000.0
t = 10 .^(-3:0.01:log10(τ))
y = f.(t)
plot(t,y)

df = t->ForwardDiff.derivative(f,t)
d2f = t->ForwardDiff.derivative(df,t)
rho(t) = (abs(1 + df(t)^2))^1.5 / abs(d2f(t))  # radius of curvature
drho = t->ForwardDiff.derivative(rho,t)        # ... and its derivative
plot(t,rho, xaxis=:log, yaxis=:log)
savefig("max_curvature.png")

julia> fzero(drho, 500)
14.142135623730953

max_curvature

To find something around t=500, the biologists definition of “inflection point” may be required.

Luigi_Marongiu · November 28, 2020, 1:59pm

Thank you all. The latter looks very interesting but:

julia> f_Tf(x, ϑ, ω, ε, μ) = ϑ - ((1/ω) * log(x/ε)) - 1/μ
f_Tf (generic function with 1 method)
julia> ε, ω, ϑ, μ, τ = [epsilon, omega, Tp, nu, tmax]
5-element Array{Float64,1}:
    1.0
    0.05
    2.392731385751533
    0.72
 4000.0
julia> x = 0.0:0.1:τ
0.0:0.1:4000.0
julia> y = f_Tf.(x, ϑ, ω, ε, μ)
40001-element Array{Float64,1}:
julia> df = t->ForwardDiff.derivative(f_Tf,t)
#19 (generic function with 1 method)
julia> d2f = t->ForwardDiff.derivative(df,t)
#21 (generic function with 1 method)
julia> rho(t) = (abs(1 + df(t)^2))^1.5 / abs(d2f(t))  # radius of curvature
rho (generic function with 1 method)
julia> drho = t->ForwardDiff.derivative(rho,t) 
#23 (generic function with 1 method)

julia> fzero(drho, 500)
ERROR: MethodError: no method matching f_Tf(::ForwardDiff.Dual{ForwardDiff.Tag{typeof(f_Tf),ForwardDiff.Dual{ForwardDiff.Tag{typeof(rho),Float64},Float64,1}},ForwardDiff.Dual{ForwardDiff.Tag{typeof(rho),Float64},Float64,1},1})
Closest candidates are:
  f_Tf(::Any, ::Any, ::Any, ::Any, ::Any) at none:1
Stacktrace:
 [1] derivative at /home/gigiux/.julia/packages/ForwardDiff/sdToQ/src/derivative.jl:14 [inlined]
 [2] (::var"#19#20")(::ForwardDiff.Dual{ForwardDiff.Tag{typeof(rho),Float64},Float64,1}) at ./none:1
 [3] rho(::ForwardDiff.Dual{ForwardDiff.Tag{typeof(rho),Float64},Float64,1}) at ./none:1
 [4] derivative at /home/gigiux/.julia/packages/ForwardDiff/sdToQ/src/derivative.jl:14 [inlined]
 [5] (::var"#23#24")(::Float64) at ./none:1
 [6] (::Roots.FnWrapper)(::Float64) at /home/gigiux/.julia/packages/Roots/Z5Mal/src/find_zero.jl:262
 [7] (::Roots.DerivativeFree{Roots.FnWrapper})(::Float64) at /home/gigiux/.julia/packages/Roots/Z5Mal/src/find_zero.jl:263 (repeats 2 times)
 [8] init_state(::Roots.Secant, ::Roots.DerivativeFree{Roots.DerivativeFree{Roots.FnWrapper}}, ::Tuple{Float64,Float64}) at /home/gigiux/.julia/packages/Roots/Z5Mal/src/derivative_free.jl:64
 [9] init_state(::Roots.Secant, ::Roots.DerivativeFree{Roots.DerivativeFree{Roots.FnWrapper}}, ::Float64) at /home/gigiux/.julia/packages/Roots/Z5Mal/src/derivative_free.jl:59
 [10] find_zero(::Roots.DerivativeFree{Roots.FnWrapper}, ::Float64, ::Roots.Secant, ::Roots.AlefeldPotraShi; tracks::Roots.NullTracks, verbose::Bool, kwargs::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}) at /home/gigiux/.julia/packages/Roots/Z5Mal/src/find_zero.jl:654
 [11] #find_zero#33 at /home/gigiux/.julia/packages/Roots/Z5Mal/src/derivative_free.jl:123 [inlined]
 [12] find_zero at /home/gigiux/.julia/packages/Roots/Z5Mal/src/derivative_free.jl:120 [inlined]
 [13] derivative_free(::Function, ::Float64; order::Int64, kwargs::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}) at /home/gigiux/.julia/packages/Roots/Z5Mal/src/fzero.jl:155
 [14] derivative_free(::Function, ::Float64) at /home/gigiux/.julia/packages/Roots/Z5Mal/src/fzero.jl:134
 [15] #fzero#48 at /home/gigiux/.julia/packages/Roots/Z5Mal/src/fzero.jl:65 [inlined]
 [16] fzero(::Function, ::Int64) at /home/gigiux/.julia/packages/Roots/Z5Mal/src/fzero.jl:63
 [17] top-level scope at none:1

Why the error?
Is 500 based on my initial eye-ball assumption?
Is the definition of rho universal (based on the curvature radius) or is specific to the equation f_Tf?

rafael.guerra · November 28, 2020, 2:10pm

If you click on the link I have provided in the earlier post, you will see the wikipedia page with the definition of the radius of curvature. It is universal for a real function (R–> R).

rafael.guerra · November 28, 2020, 4:20pm

@Luigi_Marongiu, the doc of Roots.jl says that it is a package for finding zeros of continuous scalar functions of a single real variable. If you follow the single variable syntax in my solution, there will be no such error.
For the Roots workflow, you could define a single variable auxiliary function:

f(t) = f_Tf(t, ϑ, ω, ε, μ) 
df = t->ForwardDiff.derivative(f,t)
# etc...

Geoffrey · November 28, 2020, 7:37pm

Hi,

I think what cause that the point of maximum curvature don’t fit on our intuition is that your function value is not bounded (it can give value from [-inf, inf]) so the plot we see and the impression that there is an elbow at 500 is biased by the fact that we only see a portion of the curve.

Another solution would be to find the best linear split, I mean the point that minimize the RMSE of the linear regression on the left and on the right of this point:

function f(t)
    ε, ω, ϑ, μ = (1.0, 0.05 , 2.392731385751533,  0.72)
    return ϑ - ((1 / ω) * log(t / ε)) - 1 / μ
end

# function to construct a Vandermonde matrix from a vector
function Vand(x, d)
    V = ones(length(x))
    for i in 1:d
        V = hcat(V, x .^ i)
    end
    return V
end

# function that given a test point compute the RMSE of the linear regression on the left and on the rigth of this point
function findchangepts(xi, f, xmin = eps(), xmax = 4000, n = 100)
    # take n sample of the function on the left
    l = range(xmin, xi - eps(); length = n)
    # perform the linear regression
    β = Vand(l, 1) \ f.(l)
    # compute the RMSE (using the n points)
    RMSE = sqrt(sum(t -> t * t, f.(l) .- Vand(l, 1) * β) / n)
    # do the same on the rigth
    r = range(xi + eps(), xmax; length = n)
    β = Vand(r, 1) \ f.(r)
    RMSE += sqrt(sum(t -> t * t, f.(r) .- Vand(r, 1) * β) / n)
    return RMSE
end

# look at the value for different change point
xi = range(100, 3000; length = 100)
rmse = [findchangepts(i, f) for i in xi]
plot(xi, rmse; legend = false)

Then we look for the minimum of the RMSE:

# use Optim to find the minimum of RMSE
using Optim
res = Optim.optimize(i -> findchangepts(i, f), 100, 3000)
k = Optim.minimizer(res) # 827.4879390290434

plot(x, y; legend = false)
vline!([k])
l = range(eps(), k - eps(); length = 100)
β1 = Vand(l, 1) \ f.(l)
r = range(k + eps(), τ; length = 100)
β2 = Vand(r, 1) \ f.(r)
plot!(l, Vand(l, 1) * β1)
plot!(r, Vand(r, 1) * β2)

changepts

But as you can see this point is quiet far from 500… Actually the more sample you will add on the right the more the split point will be drag on the right, so this point is totally dependent on xmax. To illustrate this you get a split point at circa 517 if you use xmax = 2500.

rafael.guerra · November 28, 2020, 10:28pm

Knowing the context/physical meaning of the curve, should help selecting some sensible criteria for the departure point.
With an arbitrary choice of setting the cut-off when the rate of change of the curve is 10 times above the minimum rate recorded at the end of the observation period, one gets t = 398:

dfmax = maximum(df.(t))     # maximum because all slopes are negative
dfcut(t) = df(t) - 10*dfmax  # cutoff at curve change rate = 10x minimum in time range
fzero(dfcut, 1)

julia> fzero(dfcut, 1)
398.1

Luigi_Marongiu · November 29, 2020, 9:02am

Thank you. The context is that of defining the optimal time Tf to administer a shot of phages in phage therapy, as defined by Payne and Jansen 2001. ϑ = Tp in the paper and corresponds to the time when the bacterial target reaches the optimal concentration. There are no limits indicated, but I guess it will go from 0 (start of the model) to a t max (end of the process), which in my case I set at 2500 (and extended to 4000 in the graph to see how it would change).

This time the definition of fzero worked, giving the value of 14.14 as reported. When using the values 517 and 398, I get models like these:

The use of 14.14 gives the same type of graphs. Even if I use 10 or 100 fold phage input, the graph does not change.
In theory, the red line should go to zero after the injection of the phage at t=398 or t=517
Given the intricacy of the subject, is there a good manual that explains how to do numeric analysis in Julia?

Tamas_Papp · November 29, 2020, 9:59am

Please note that numerical analysis is a very broad field; you may need to be more specific.

Incidentally, in the above graphs I don’t see the red lines crossing 0 at all. Are you sure that the problem is well-defined?

rafael.guerra · November 29, 2020, 10:04am

Concerning documentation on numerical analysis in Julia, the following course by Jan Verschelde from UIC seems to be a very good resource.
As you are simulating dynamic processes, you should look also at the detailed documentation of DifferentialEquations.jl.

Eben60 · November 29, 2020, 4:44pm

@Luigi_Marongiu

I couldn’t find the definition of “flection point”. Do you mean “inflection point”, i.e. the point where the at which the curvature changes sign? In this case it looks like your curve doesn’t have one.
Your plots (multiple colors) at different injection times look the same. If they are indeed what you get, then it’s probably due to some bug.
If I were you, I’d try to simulate the whole dynamic process using DifferentialEquations.jl as @rafael.guerra suggested, define your objective function, which could be some concentration of bacteria at a definite time point, or some integral over a time span, or whatever, and then try to vary the injection point and look for the optimum. It looks like your problem is 1D (time being the independent variable) and involves only a few dependent variables, so it’s probably not excessively compute intensive.

Luigi_Marongiu · November 29, 2020, 5:09pm

Yes “flection point” was a typo, sorry. It was correct the second time. The problem is indeed that the red line DOES NOT reach 0 (if I extend the y-axis to 1, the result does not change). The idea is that if the injection point is correct, the red lines should go to zero and stay there. Here there is no difference between 517, 398 or 14. The graphs are generated with DifferentialEquations.jl. I’ll see for the video course. Thanks

Eben60 · November 29, 2020, 5:18pm

And with 1000000000 fold phage input?

rafael.guerra · November 29, 2020, 6:33pm

@Luigi_Marongiu, the constants in your equation: y = ϑ - (1/ω) * log(x/ε) - 1/μ, do not have units defined. You should check the consistency of your units. Maybe some conversion factors are missing?

Luigi_Marongiu · November 29, 2020, 8:18pm

Same thing:

The whole model is:

function multiInfect!(du, u, p, t)
    μ, ν, κ, φ, ψ, ω, β, Β, η, Η = p
    #=
    du[1] = resident susceptible
    du[2] = resident infected
    du[3] = resident phage
    du[4] = invader susceptible
    du[5] = invader infected
    du[6] = invader phage
    =#
    sumBact = u[1]+u[2]+u[4]+u[5]
    du[1] = (μ * u[1]) * (1 - sumBact/κ) - (φ * u[1] * u[3]) - (ω * u[1])
    du[2] = (φ * u[1] * u[3]) - (η * u[2]) - (ω * u[2])
    du[3] = (β * η * u[2]) - (φ * u[1] * u[3]) - (ω * u[3])

    du[4] = (ν * u[4]) * (1 - sumBact/κ) - (ψ * u[4] * u[6]) - (ω * u[4])
    du[5] = (ψ * u[4] * u[6]) - (Η * u[5]) - (ω * u[5])
    du[6] = (Β * Η * u[5]) - (ψ * u[4] * u[6]) - (ω * u[6])
end
#=
    From: Payne J Theor Biol 2001;208:37
    μ       = growth rate (replication coef), a
    φ       = adsorption rate (transmission coef), b
    η       = lysis rate, k
    β       = burst size, L
    ω       = outflow (decay rate phage), m
    σ       = initial bacterial concentratoin, S0
    τ       = maximum running time
    ϑ       = Tp
=#
function ϝ_Tp(p) # proliferation-onset time
    ω, η, μ, β, φ, σ = p
    x = (1/μ) * log( ((η-μ)*ω) / (φ*η*β*σ) )
    return x
end
function ϝ_Vf(p)
    ε, ω, ϑ, μ = p
    # solution by looping
    Y = []
    X = 0.0:0.1:2*ϑ
    for x in X
        y = ε * exp( ω * (ϑ-x) + (ω/μ) * (exp( -μ * (ϑ-x) ) - 1) )
        push!(Y, y)
    end
    return minimum(Y)
end
ϝ_Tf(x, ϑ, ω, ε, μ) = ϑ - ((1/ω) * log(x/ε)) - 1/μ

# parameters
X     = 2500        # tmax for plot
mu    = 0.47        # maximum growth rate susceptible (B. longum)
nu    = 0.72        # maximum growth rate resistant (E. coli)
kappa = 2.2*10^7    # maximum population density
phi   = 10.0^-9     # adsorption rate
omega = 0.05        # outflow
eta   = 1.0         # lyse rate
beta  = 50.0        # burst size
tau   = 3.62        # latency time
tmax  = 4000.0      # time span 0-tmax
s0    = 50000.0     # initial susceptible population
i0    = 0.0         # initial infected population
v0    = 80.0        # initial phage population
psi   = 10.0^-9       # adsorption rate alternative
Eta   = 1.0           # lyse rate alternative
Beta  = 50.0          # burst size alternative
r_s0    = 50000.0     # initial resident susceptible population
r_i0    = 0.0         # initial resident infected population
r_v0    = 1000.0      # initial resident phage population
i_s0    = 0.0         # initial invader susceptible population
i_i0    = 0.0         # initial invader infected population
i_v0    = 0.0         # initial invader phage population
t_in    = 1000        # time of inoculum
i_in    = 500         # amount of invader
epsilon = 1           # concentration equivalent to 1 phage

# calculate critical parameters
Tp = ϝ_Tp([omega, Eta, nu, Beta, psi, r_s0])
Vf = ϝ_Vf([epsilon, omega, Tp, nu])
#=
 Tf (time of injection) from this post
=#
 
# implement
Tf = 398
tspan = (0.0, tmax)
u0 = [r_s0, r_i0, r_v0, i_s0, i_i0, i_v0]
parms = [mu, nu, kappa, phi, psi, omega, beta, Beta, eta, Eta]
condition_1(u, t, integrator) = t-t_in              # time invader inoculum
condition_2(u, t, integrator) = t-(t_in+Tf)         # time phage activation
affect_1!(integrator) = integrator.u[4] += i_in     # amount invader inoculum
affect_2!(integrator) = integrator.u[6] += 1000000000       # amount phage
cb1 = ContinuousCallback(condition_1, affect_1!)
cb2 = ContinuousCallback(condition_2, affect_2!)
delay = CallbackSet(cb1, cb2)
prob = ODEProblem(multiInfect!, u0, tspan, parms)
soln = solve(prob, AutoVern7(Rodas5()), callback=delay, tstops=[1000])

Eben60 · November 29, 2020, 8:53pm

As I understand it, physically (or, in your case rather “biologically”), very high phage dose at an arbitrary time point should eliminate all bacteria. If that doesn’t happen in your solution, then most probably due to some silly bug somewhere.

Luigi_Marongiu · November 30, 2020, 7:03am

Yes, it is worrisome. To give more context, The graph reported in the original paper for the definition of Tf gives this profile:

To note that the time t=20 is the end of the simulation (tmax) that in my case is 2500.
and the formula for Tf (ϑ - ((1/ω) * log(x/ε)) - 1/μ) [with ϑ = Tp = (1/μ) * log( ((η-μ)ω) / (φηβσ) )] gives a function that is similar to the hyperbolic or negative exponential:

The manuals report a/2b or a/e as important points, which should be similar to the inflection point I was looking for. The problem is to convert Tf into the hyperbolic or exponential format to find a and b…
Or perhaps it is just that the model reported in the paper is fundamentally different from the one I am studying.