bifurcationkit / BifurcationKit.jl

A Julia package to perform Bifurcation Analysis
https://bifurcationkit.github.io/BifurcationKitDocs.jl/stable
Other
301 stars 36 forks source link

Questions regarding, and oddities with, the bifurcationdiagram function #30

Closed TorkelE closed 3 years ago

TorkelE commented 3 years ago

So, I have been playing around with the bifurcationdiagram function. It's really cool! However, there are a few things I do not understand, and a few other things which seem weird, which I have not been able to understand from reading the docs.

For a starter, I run this code to generate my bifurcation diagram:

# Fetch the required packages.
using Revise, Plots
using BifurcationKit, Setfield, ForwardDiff
using Catalyst, OrdinaryDifEq
const BK = BifurcationKit

# The input (a system, a parameter set, a target parameter, and a range to varry it over).
system = @reaction_network begin
    v0+v*(S*σ)^n/((S*σ)^n+(D*A)^n+K^n), ∅ → σ
    d, σ → ∅
    (τ*σ,τ), ∅ ↔ A
end S D τ v v0 K n d;
pRange = (0.1,20.)
parameter_vals = [10.,9.,0.005,2.,0.01,20.,3,0.05]
p = :S
pRange = (0.1,20.)

# Extract required stuff from the input.
fun = ODEFunction(convert(ODESystem,system),jac=true)
parameters = copy(parameter_vals); parameters[1] = pRange[1]
p_idx = findfirst(Symbol.(system.ps) .== p)

# Differentiate the function.
D(f, x, p, dx) = ForwardDiff.derivative(t -> f(x .+ t .* dx, p), 0.)
FbpSecBif(u, p) = fun(u,p,0)
dFbpSecBif(x,p) =  ForwardDiff.jacobian( z-> FbpSecBif(z,p), x)
d1FbpSecBif(x,p,dx1) = D((z, p0) -> FbpSecBif(z, p0), x, p, dx1)
d2FbpSecBif(x,p,dx1,dx2) = D((z, p0) -> d1FbpSecBif(z, p0, dx1), x, p, dx2)
d3FbpSecBif(x,p,dx1,dx2,dx3) = D((z, p0) -> d2FbpSecBif(z, p0, dx1, dx2), x, p, dx3)
jet = (FbpSecBif, dFbpSecBif, d2FbpSecBif, d3FbpSecBif)

# Gets a starting fixed point, sets options.
x0 = solve(ODEProblem(system,[1.,1.],(0.,2000),parameters),Rosenbrock23()).u[end];
opt_newton = NewtonPar(tol = 1e-9, verbose = false, maxIter = 20)
opts_br = ContinuationPar(pMin=pRange[1], pMax=pRange[2], maxSteps = 50000,
            dsmin = (pRange[2]-pRange[1])*1e-9, dsmax = (pRange[2]-pRange[1])*1e-3, ds = (pRange[2]-pRange[1])*1e-9, 
            detectBifurcation = 3, nev = 2, newtonOptions = opt_newton, nInversion = 4, 
            tolBisectionEigenvalue = 1e-8, dsminBisection = 1e-9);

# Creates the bifurcation diagram.
diagram = bifurcationdiagram(jet..., x0, parameters, (@lens _[p_idx]), 2, (args...) -> setproperties(opts_br; pMin=pRange[1], pMax=pRange[2], ds = (pRange[2]-pRange[1])*1e-9, dsmax = (pRange[2]-pRange[1])*1e-3, nInversion = 8, detectBifurcation = 3,dsminBisection =1e-18, tolBisectionEigenvalue=1e-11, maxBisectionSteps=20, newtonOptions = (@set opt_newton.verbose=false));)

now this generates a set of error messages, which seems very worrying:

##################################################
---> Automatic computation of bifurcation diagram

────────────────────────────────────────────────────────────────────────────────
--> New branch level = 2, dim(Kernel) = 1, code = 0, from bp #1 at p = 9.063030706988087
────────────────────────────────────────────────────────────────────────────────
--> New branch level = 2, dim(Kernel) = 1, code = 0, from bp #2 at p = 4.228100778475947
────────────────────────────────────────────────────────────────────────────────
--> New branch level = 2, dim(Kernel) = 2, code = 0, from bp #3 at p = 10.755971961392246
#####################################################
--> Hopf Normal form computation
--> Hopf bifurcation point is: SuperCritical
┌ Error: It seems the point is a Saddle-Node bifurcation. The normal form is (a = 0.0018436142964569494, b1 = 0.015793706727377507, b2 = 0.3066201972182164, b3 = 0.41398034761966024).
└ @ BifurcationKit /home/torkelloman/.julia/packages/BifurcationKit/CO5k1/src/NormalForms.jl:203
┌ Error: It seems the point is a Saddle-Node bifurcation. The normal form is (a = -0.007753919955655971, b1 = -0.002178682035144704, b2 = 0.0016086450861372025, b3 = -0.00015770032398444782).
└ @ BifurcationKit /home/torkelloman/.julia/packages/BifurcationKit/CO5k1/src/NormalForms.jl:203
┌ Error: MethodError(BifurcationKit.var"#multicontinuation##kw"(), ((nev = 2, issymmetric = false, usedeflation = false, Teigvec = Array{Float64,1}, ζs = nothing, δp = nothing, verbosedeflation = false, scaleζ = LinearAlgebra.norm, maxIterDeflation = 50, perturb = identity, plotSolution = BifurcationKit.var"#344#349"{BifurcationKit.var"#344#345#350"{Int64,GenericBifPoint{Float64,Array{Float64,1}},NamedTuple{(:current, :maxlevel),Tuple{Int64,Int64}},BifurcationKit.var"#342#347"{BifurcationKit.var"#342#343#348"},String,BifDiagNode{ContResult{Float64,Array{Complex{Float64},1},Array{Complex{Float64},2},GenericBifPoint{Float64,Array{Float64,1}},Nothing,Nothing,Array{Float64,1},Setfield.IndexLens{Tuple{Int64}}},Array{BifDiagNode,1}}}}(BifurcationKit.var"#344#345#350"{Int64,GenericBifPoint{Float64,Array{Float64,1}},NamedTuple{(:current, :maxlevel),Tuple{Int64,Int64}},BifurcationKit.var"#342#347"{BifurcationKit.var"#342#343#348"},String,BifDiagNode{ContResult{Float64,Array{Complex{Float64},1},Array{Complex{Float64},2},GenericBifPoint{Float64,Array{Float64,1}},Nothing,Nothing,Array{Float64,1},Setfield.IndexLens{Tuple{Int64}}},Array{BifDiagNode,1}}}(3, GenericBifPoint{Float64,Array{Float64,1}}
│   type: Symbol hopf
│   idx: Int64 1413
│   param: Float64 10.755971961392246
│   norm: Float64 35.945080086698646
│   printsol: Float64 35.945080086698646
│   x: Array{Float64}((2,)) [25.417009879598147, 25.417009879598147]
│   tau: BorderedArray{Array{Float64,1},Float64}
│   ind_ev: Int64 2
│   step: Int64 1412
│   status: Symbol converged
│   δ: Tuple{Int64,Int64}
│   precision: Float64 2.4654062258377962e-6
│   interval: Tuple{Float64,Float64}
│ , (current = 1, maxlevel = 2), BifurcationKit.var"#342#347"{BifurcationKit.var"#342#343#348"}(BifurcationKit.var"#342#343#348"()), "0", Bifurcation diagram. Root branch (level 1) has 0 children and is such that:
│ Branch number of points: 1947
│ Branch of Equilibrium
│ Bifurcation points:
│  (ind_ev = index of the bifurcating eigenvalue e.g. `br.eig[idx].eigenvals[ind_ev]`)
│ - #  1,    bp at p ≈ +9.06303071 ∈ (+9.05981691, +9.06303071), |δp|=3e-03, [converged], δ = ( 1,  0), step = 358, eigenelements in eig[359], ind_ev =   1
│ - #  2,    bp at p ≈ +4.22810078 ∈ (+4.22810078, +4.22810078), |δp|=3e-12, [converged], δ = ( 1,  0), step = 568, eigenelements in eig[569], ind_ev =   2
│ - #  3,  hopf at p ≈ +10.75597196 ∈ (+10.75596950, +10.75597196), |δp|=2e-06, [converged], δ = (-2, -2), step = 1412, eigenelements in eig[1413], ind_ev =   2
│ Fold points:
│ - #  1, fold at p ≈ 9.06308297 ∈ (9.06308297, 9.06308297), |δp|=-1e+00, [    guess], δ = ( 0,  0), step = 354, eigenelements in eig[354], ind_ev =   0
│ - #  2, fold at p ≈ 4.22810078 ∈ (4.22810078, 4.22810078), |δp|=-1e+00, [    guess], δ = ( 0,  0), step = 565, eigenelements in eig[565], ind_ev =   0
│ ))), BifurcationKit.multicontinuation, FbpSecBif, dFbpSecBif, Branch number of points: 1947
│ Branch of Equilibrium
│ Bifurcation points:
│  (ind_ev = index of the bifurcating eigenvalue e.g. `br.eig[idx].eigenvals[ind_ev]`)
│ - #  1,    bp at p ≈ +9.06303071 ∈ (+9.05981691, +9.06303071), |δp|=3e-03, [converged], δ = ( 1,  0), step = 358, eigenelements in eig[359], ind_ev =   1
│ - #  2,    bp at p ≈ +4.22810078 ∈ (+4.22810078, +4.22810078), |δp|=3e-12, [converged], δ = ( 1,  0), step = 568, eigenelements in eig[569], ind_ev =   2
│ - #  3,  hopf at p ≈ +10.75597196 ∈ (+10.75596950, +10.75597196), |δp|=2e-06, [converged], δ = (-2, -2), step = 1412, eigenelements in eig[1413], ind_ev =   2
│ Fold points:
│ - #  1, fold at p ≈ 9.06308297 ∈ (9.06308297, 9.06308297), |δp|=-1e+00, [    guess], δ = ( 0,  0), step = 354, eigenelements in eig[354], ind_ev =   0
│ - #  2, fold at p ≈ 4.22810078 ∈ (4.22810078, 4.22810078), |δp|=-1e+00, [    guess], δ = ( 0,  0), step = 565, eigenelements in eig[565], ind_ev =   0
│ , BifurcationKit.HopfBifPoint{Array{Float64,1},Float64,Float64,Array{Float64,1},Setfield.IndexLens{Tuple{Int64}},Array{Complex{Float64},1},Array{Complex{Float64},1},NamedTuple{(:a, :b),Tuple{Complex{Float64},Complex{Float64}}}}([25.417009879598147, 25.417009879598147], 10.755971961392246, 0.015804494505662976, [10.755971961392246, 9.0, 0.005, 2.0, 0.01, 20.0, 3.0, 0.05], (@lens _[1]), Complex{Float64}[-0.9571165934996202 + 0.0im, -0.0842773205962544 + 0.2771735190829004im], Complex{Float64}[-0.5213874537927115 - 0.16204406999169177im, 0.009147349963915332 + 1.8102101966238962im], (a = -0.004816693942145319 + 2.6835568552196435e-5im, b = 5.225986319528161e-5 + 0.0001546547428260476im), :SuperCritical), ContinuationPar{Float64,DefaultLS,DefaultEig{typeof(real)}}
│   dsmin: Float64 1.99e-8
│   dsmax: Float64 0.019899999999999998
│   ds: Float64 1.99e-8
│   theta: Float64 0.5
│   doArcLengthScaling: Bool false
│   gGoal: Float64 0.5
│   gMax: Float64 0.8
│   thetaMin: Float64 0.001
│   a: Float64 0.5
│   tangentFactorExponent: Float64 1.5
│   pMin: Float64 0.1
│   pMax: Float64 20.0
│   maxSteps: Int64 50000
│   finDiffEps: Float64 1.0e-9
│   newtonOptions: NewtonPar{Float64,DefaultLS,DefaultEig{typeof(real)}}
│   η: Float64 150.0
│   saveToFile: Bool false
│   saveSolEveryStep: Int64 0
│   nev: Int64 2
│   saveEigEveryStep: Int64 1
│   saveEigenvectors: Bool true
│   plotEveryStep: Int64 10
│   precisionStability: Float64 1.0e-10
│   detectFold: Bool true
│   detectBifurcation: Int64 3
│   dsminBisection: Float64 1.0e-18
│   nInversion: Int64 8
│   maxBisectionSteps: Int64 20
│   tolBisectionEigenvalue: Float64 1.0e-11
│   detectLoop: Bool false
│ ), 0x0000000000006dde)
└ @ BifurcationKit /home/torkelloman/.julia/packages/BifurcationKit/CO5k1/src/bifdiagram/BifurcationDiagram.jl:115

However, I am still actually able to plot the diagram, and it looks correct:

# Plots the diagram.
plot(diagram)

bif_example

Now the first query is what all these error messages are? What went wrong, but in such a way that the plot still turned out fine?

Next, are there options for modifying the plot diagram? The unstable/stable regions are not very easy to distinguish due to the line thickness and if you change the line thickness (e.g. to make it thicker, and easier to see) the difference disappears. Typically I use a different colour, or dashed lines, for the unstable states. Is there an easy way of doing this? Alternatively, is there a way of extracting the separate branches, as divided by the bifurcation points (I only managed to get the full path)?

Also, the bifurcation diagram labels the bifurcation dots weirdly doesn't look correct, the dots themself are missing in the label (and the dots in the diagram, especially the Hopf one, are very small).

Next, a question. It seems to me that I am setting most of the options twice:

opt_newton = NewtonPar(tol = 1e-9, verbose = false, maxIter = 20)
opts_br = ContinuationPar(pMin=pRange[1], pMax=pRange[2], maxSteps = 50000,
            dsmin = (pRange[2]-pRange[1])*1e-9, dsmax = (pRange[2]-pRange[1])*1e-3, ds = (pRange[2]-pRange[1])*1e-9, 
            detectBifurcation = 3, nev = 2, newtonOptions = opt_newton, nInversion = 4, 
            tolBisectionEigenvalue = 1e-8, dsminBisection = 1e-9);

diagram = bifurcationdiagram( ..., (args...) -> setproperties(opts_br; pMin=pRange[1], pMax=pRange[2], ds = (pRange[2]-pRange[1])*1e-9, dsmax = (pRange[2]-pRange[1])*1e-3, nInversion = 8, detectBifurcation = 3,dsminBisection =1e-18, tolBisectionEigenvalue=1e-11, maxBisectionSteps=20, newtonOptions = (@set opt_newton.verbose=false));)

Would I always want to do this, or could I go for something simpler?

Finally, when I created diagrams using continuation I didn't have to define so many derivatives as here:

D(f, x, p, dx) = ForwardDiff.derivative(t -> f(x .+ t .* dx, p), 0.)
FbpSecBif(u, p) = fun(u,p,0)
dFbpSecBif(x,p) =  ForwardDiff.jacobian( z-> FbpSecBif(z,p), x)
d1FbpSecBif(x,p,dx1) = D((z, p0) -> FbpSecBif(z, p0), x, p, dx1)
d2FbpSecBif(x,p,dx1,dx2) = D((z, p0) -> d1FbpSecBif(z, p0, dx1), x, p, dx2)
d3FbpSecBif(x,p,dx1,dx2,dx3) = D((z, p0) -> d2FbpSecBif(z, p0, dx1, dx2), x, p, dx3)
jet = (FbpSecBif, dFbpSecBif, d2FbpSecBif, d3FbpSecBif)

rather, in those cases, things worked OK with only the Jacobian (which is autogenerated by the @reaction_network macro). Why are these required here? Would it be possible to do without them (and what would the negatives of that be?).

rveltz commented 3 years ago

Hi,

Thank you for trying hard! I appreciate it.

Before I answer all the comments: bifurcationdiagram will give you nothing for your diagram. Your diagram is too simple... bifurcationdiagram allows you to branch recursively at stationary bifurcation points (e.g. not Hopf points) and you dont have any. In this case, you are better off using deflated continuation with detection of bifurcation points (see other issue you posted on).

My error messages are much cleaner (I am on master):

────────────────────────────────────────────────────────────────────────────────
--> New branch level = 2, dim(Kernel) = 1, code = 0, from bp #1 at p = 9.063030706988087
- #  1,    bp at p ≈ +9.06303071 ∈ (+9.05981691, +9.06303071), |δp|=3e-03, [converged], δ = ( 1,  0), step = 358, eigenelements in eig[359], ind_ev =   1
┌ Error: It seems the point is a Saddle-Node bifurcation. The normal form is (a = 0.0018436142964569494, b1 = 0.015793706727377507, b2 = 0.3066201972182164, b3 = 0.41398034761966024).
└ @ BifurcationKit ~/work/prog_gd/julia/dev/dev1/BifurcationKit/src/NormalForms.jl:203
────────────────────────────────────────────────────────────────────────────────
--> New branch level = 2, dim(Kernel) = 1, code = 0, from bp #2 at p = 4.228100778475947
- #  2,    bp at p ≈ +4.22810078 ∈ (+4.22810078, +4.22810078), |δp|=3e-12, [converged], δ = ( 1,  0), step = 568, eigenelements in eig[569], ind_ev =   2
┌ Error: It seems the point is a Saddle-Node bifurcation. The normal form is (a = -0.007753919955655971, b1 = -0.002178682035144704, b2 = 0.0016086450861372025, b3 = -0.00015770032398444782).
└ @ BifurcationKit ~/work/prog_gd/julia/dev/dev1/BifurcationKit/src/NormalForms.jl:203
────────────────────────────────────────────────────────────────────────────────
--> New branch level = 2, dim(Kernel) = 2, code = 0, from bp #3 at p = 10.755971961392275
- #  3,  hopf at p ≈ +10.75597196 ∈ (+10.75596950, +10.75597196), |δp|=2e-06, [converged], δ = (-2, -2), step = 1412, eigenelements in eig[1413], ind_ev =   2

Basically it is trying to branch from the bifurcation points which are folds or Hopf. It only considers the fold points which are first labelled as regular bifurcation points until it computes the normal form and see they are folds points. It sends an error because sometimes, the normal form is close to a Pitchfork one and I want the user to be able to know. I should probably find another way to tell so that the user can turn it off.

So basically: nothing went wrong!

Plotting

Here is a very important point concerning plotting. When the bifurcation points are displayed with square, it means they are poorly located ( bisection was not used for those). For the size, you might prefer

plot(diagram, markersize = 3)

if you change the line thickness (e.g. to make it thicker, and easier to see) the difference disappears.

I can add an option to the plotting recipe for that. Super easy.

Typically I use a different colour, or dashed lines, for the unstable states.

I tried before using dashed for unstable but it is too unreliable. See here unless you have a suggestion. As for different colours, it becomes a mess with many branches.

Is there an easy way of doing this?

The plotting recipe is quite simple, see here

Alternatively, is there a way of extracting the separate branches, as divided by the bifurcation points (I only managed to get the full path)?

Well, I never needed it so I did not write it.

Also, the bifurcation diagram labels the bifurcation dots weirdly doesn't look correct, the dots themself are missing in the label (and the dots in the diagram, especially the Hopf one, are very small).

It is because the fold points overlay them. You can do

plot(diagram, markersize = 3, plotfold=false)

Next, a question. It seems to me that I am setting most of the options twice:

I'd say yes. Passing (args...) -> opts_br should be enough. I could even use dispatch to simply pass opts_br.

Why are these required here? Would it be possible to do without them (and what would the negatives of that be?).

They are used for the computation of normal forms. No automatic branching without them. You can pass a finitedifferences version in any cases. The most crutial thing is to know the bifurcation point location precisely. For this, you only need to set nInversion high enough. Then, if you are not super precise on d2F and d3F, it should be able to branch from simple BP (dimKernel = 1) pretty easily.

TorkelE commented 3 years ago

This is incredibly helpful. I see what you mean fir the more complicated bifurcation diagrams such as in https://rveltz.github.io/BifurcationKit.jl/dev/BifurcationDiagram/ yeah, I can see that that's a different level of stuff. I'll move over to deflated continuations, as you say, that should probably be enough for everything you encounter in biochemistry.

Is there a list of plot options, for when plotting these things, that I haven't found (with stuff like markersize)? The only real reason I would want to extract separate branches was that that would give me the option of customising the plotting of the branches.

rveltz commented 3 years ago

I would use deflated continuation (DC) when possible, it is quite powerful even considering its drawbacks:

  1. it can be slow because you "waste" time testing poor initial guess
  2. it cannot branch periodic orbit (yet)
  3. you get sometimes too many solutions for you to handle. I never showed the same as this example but computed with DC, it would scare the rabbit away :D

Nevertheless, for "small" systems <100, I'd say it is superior to regular continuation.

Is there a list of plot options, for when plotting these things, that I haven't found

Yes, see here. Then, the plotting options are recipes. Hence, you can do plot(br) but also plot!(br) or even scatter(br). Finally, you can pass all options of Plots.jl (linewidth,markersize...)

rveltz commented 3 years ago

I added a keyword argument to plot to tune the linewidth of the stable / unstable part, see here.

rveltz commented 3 years ago

Should we close this?

TorkelE commented 3 years ago

Got a bit lost over the holidays. Yes, we should.