bifurcationkit / bifurcationkit.jl Goto Github PK

julia> res = @time probMPsh(par_hopf)(initpo, initpo)
ERROR: Event repeated at the same time. Please report this error
Stacktrace:
 [1] apply_callback!(::OrdinaryDiffEq.ODEIntegrator{Tsit5,false,Array{Float64,1},Float64,NamedTuple{(:r, :μ, :ν, :c3, :c5),NTuple{5,Float64}},Float64,Float64,Float64,Array{Array{Float64,1},1},ODESolution{Float64,2,Array{Array{Float64,1},1},Nothing,Nothing,Array{Float64,1},Array{Array{Array{Float64,1},1},1},ODEProblem{SubArray{Float64,1,Array{Float64,2},Tuple{Base.Slice{Base.OneTo{Int64}},Int64},true},Tuple{Float64,Float64},false,NamedTuple{(:r, :μ, :ν, :c3, :c5),NTuple{5,Float64}},ODEFunction{false,typeof(Fsl),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},DiffEqBase.StandardODEProblem},Tsit5,OrdinaryDiffEq.InterpolationData{ODEFunction{false,typeof(Fsl),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Array{Array{Float64,1},1},Array{Float64,1},Array{Array{Array{Float64,1},1},1},OrdinaryDiffEq.Tsit5ConstantCache{Float64,Float64}},DiffEqBase.DEStats},ODEFunction{false,typeof(Fsl),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},OrdinaryDiffEq.Tsit5ConstantCache{Float64,Float64},OrdinaryDiffEq.DEOptions{Float64,Float64,Float64,Float64,typeof(DiffEqBase.ODE_DEFAULT_NORM),typeof(opnorm),CallbackSet{Tuple{VectorContinuousCallback{PseudoArcLengthContinuation.var"#pSection#599"{typeof(section!)},PseudoArcLengthContinuation.var"#affect!#600",Nothing,typeof(DiffEqBase.INITIALIZE_DEFAULT),Float64,Int64,Nothing}},Tuple{}},typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN),typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE),typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK),DataStructures.BinaryHeap{Float64,DataStructures.LessThan},DataStructures.BinaryHeap{Float64,DataStructures.LessThan},Nothing,Nothing,Int64,Array{Float64,1},Array{Float64,1},Array{Float64,1}},Array{Float64,1},Float64,DiffEqBase.CallbackCache{Array{Float64,N} where N,Array{Float64,N} where N}}, ::VectorContinuousCallback{PseudoArcLengthContinuation.var"#pSection#599"{typeof(section!)},PseudoArcLengthContinuation.var"#affect!#600",Nothing,typeof(DiffEqBase.INITIALIZE_DEFAULT),Float64,Int64,Nothing}, ::Float64, ::Float64, ::Int64) at /Users/rveltz/.julia/packages/DiffEqBase/YIwj5/src/callbacks.jl:667

using DifferentialEquations, DiffEqOperators, ForwardDiff
	using PseudoArcLengthContinuation, LinearAlgebra, Parameters, Setfield
	const PALC = PseudoArcLengthContinuation

norminf = x -> norm(x, Inf)
function Fsl!(f, u, p)
	@unpack r, μ, ν, c3, c5 = p
	u1 = u[1]
	u2 = u[2]

	ua = u1^2 + u2^2

	f[1] = r * u1 - ν * u2 - ua * (c3 * u1 - μ * u2) - c5 * ua^2 * u1
	f[2] = r * u2 + ν * u1 - ua * (c3 * u2 + μ * u1) - c5 * ua^2 * u2

	return f
end

function Fsl(x, p, t)
	f = similar(x)
	Fsl!(f, x, p)
end

Fode(f, x, p, t) = Fsl!(f, x, p)
Jsl(x, dx, p) = ForwardDiff.derivative(t -> Fsl(x .+ t .* dx, p), 0.)
# JacOde(J, x, p, t) = copyto!(J, Jsl(x, p))

par_sl = (r = 0.5, μ = 0., ν = 1.0, c3 = 1.0, c5 = 0.0,)
u0 = [.001, .001]

par_hopf = (@set par_sl.r = 0.1)
probsundials = ODEProblem(Fsl, u0, (0., 100.), par_hopf)

function section!(out, x)
	M = 4
	for ii=1:M
		th = 2pi*ii/M
		tmp = [-sin(th), cos(th)]
		out[ii] = dot(x, tmp)
	end
	out
end

# standard simple shooting
M = 4
probMPsh = p -> PALC.PoincareShootingProblem(u -> Fsl(u, p, 0.), p, probsundials, Tsit5(), M, section!; atol = 1e-10, rtol = 1e-9)

r = 1.04

initpo = Float64[]
for ii=1:M
	th = 2pi*ii/M
	push!(initpo, r*cos(th))
	push!(initpo, r*sin(th))

end

res = @time probMPsh(par_hopf)(initpo);
norminf(res)
# 0.005197 seconds (61.15 k allocations: 4.358 MiB)
# 0.4953547991423255

using DiffEqBiological,BifurcationKit
using Setfield,LinearAlgebra,Plots

system = @reaction_network begin
     k₁,		Y → 2X
     2k₂,		2X → X + Y
     k₃,		X + Y → Y
     k₄,		X → 0
end k₁ k₂ k₃ k₄

# get rate function and jacobian, for use in continuation
F! = oderhsfun(system)
function F(u::Vector{T},p::NamedTuple) where T<:Number
	f = similar(u)
	F!(f,u,[ p[key] for (key, value) ∈ paramsmap(system) ],Inf)
	return f
end

∂Fu! = jacfun(system)
function ∂Fu(u::Vector{T},p::NamedTuple) where T<:Number
	J = zeros(T,length(u),length(u))
	∂Fu!(J,u,[ p[key] for (key, value) ∈ paramsmap(system) ],Inf)
	return J
end

# define the parameter direction that you want to explore
pRange = range(0,12,length=200)
parameter_name = (@lens _.k₁)
p = (k₁=0.0, k₂=1.0, k₃=1.0, k₄=1.5)

# options for branch continuation along parameter direction
options = ContinuationPar( maxSteps = 5000,
	dsmax=10*step(pRange), dsmin=step(pRange), ds=step(pRange),
	pMin = minimum(pRange), pMax=maximum(pRange),
	newtonOptions = NewtonPar(tol = 1e-8),
	saveEigenvectors = false
)

# options for detecting multiple branches
deflation = DeflationOperator(2.0, dot, 1.0, [ zeros(length(species(system))) ])
termination_condition(x, f, J, residual, iteration, niter, options; kwargs...) = residual <1e3
root_search(x,p,id) = x .+ rand(length(x))

##################################### main continuation method
@assert( det(∂Fu(randn(length(species(system))),p)) ≠ 0, "species must be linearly independent")
branches, = continuation( F, ∂Fu, p, parameter_name, options, deflation; showplot=false,
	callbackN = termination_condition, perturbSolution = root_search )

#################################### display results
plot(size=(500,500),ylabel="steady states", xlabel="parameter $(first(typeof(parameter_name).parameters))")
plot!(branches,lw=3,label="")

# Fetch stuff
using PseudoArcLengthContinuation, SparseArrays, LinearAlgebra, Plots, Setfield, UnPack
const PALC = PseudoArcLengthContinuation
using DiffEqBiological # Required to generate my system

# Taken from the docs. Don't know why we are doing this.
normInf = x -> norm(x, Inf)

# Generates my functions
# A reactio network. This generates the stuff for a model system with a single variable and 5 parameters.
rn = @reaction_network begin
	v0+hill(X,v,K,n), 0 --> X
	d, X --> 0
end v0 v K n d

# The reaction network structures contain the right functions, however, they are designed for DifferentialEquations and requires modification.
# The original functions do not return the output, but rather modifes an input vector/matrix with the desired values. Also, the parameters works differently (I tried it the natural way, ran into problems, but I will worry about that when I got the main methodology working).
function F(x,p)
	@unpack v0,v,K,n,d = p
	output = copy(x)
	rn.f(output,x,[v0, v, K, n, d],0.)
	return output
end
function J(x,p)
	@unpack v0,v,K,n,d = p
	output = zeros(length(x),length(x))
	rn.jac(output,x,[v0, v, K, n, d],0.)
	return output
end

# Just taken from the docs. I modified pMin and pMax to 0.1 and 10.0. I sat pMin to 0.1 because that was my best guess at the initial value.
opts = PALC.ContinuationPar(dsmax = 0.1, dsmin = 1e-5, ds = -0.001, maxSteps = 130, pMin = 0.1, pMax = 10.0, saveSolEveryNsteps = 0, computeEigenValues = true, detectBifurcation = true, newtonOptions = NewtonPar(tol = 1e-8, verbose = true))

# Set the parameters. I sat v to 5. because I thought that might be the final value.
pars = (v0=0.1, v=5., K=0.5, n=3., d=1.2)

# Tries to define the iterator. The fixed point of the system for v=0.1 is 0.08333333333333334.
iter = PALC.PALCIterable(F, J, [0.08333333333333334], pars, (@lens _.v), opts; verbosity = 2)

# This is just taken directly from the example.
resp = Float64[]
resx = Float64[]
for state in iter
	push!(resx, getx(state)[1])
	push!(resp, getp(state))
end

#This plot
plot(resp, resx; label = "", xlabel = "v")

# Fetch packages
using ModelingToolkit, OrdinaryDiffEq, BifurcationKit, Plots, Setfield, LinearAlgebra

# Declare model.
@parameters t s d T v0 n
@variables X(t) A1(t) A2(t) A3(t)
D = Differential(t)

eqs = [D(X) ~ v0 + ((s*X)^n)/((s*X)^n + (d*A3)^n + 1) - X,
       D(A1) ~ (X-A1)/T,
       D(A2) ~ (A1-A2)/T,
       D(A3) ~ (A2-A3)/T]

@named sys = ODESystem(eqs);

# Prepares input 
params_input = [0.25,1.0,0.1,3,0.1]  # Order is "d,s,v0,n,T", check in "sys.ps". Should be fine.
p_span = (1.,4.);
#odefun = ODEFunction(convert(ODESystem,network),jac=true)
odefun = ODEFunction(sys,jac=true)
F = (u,p) -> odefun(u,p,0)      
J = (u,p) -> odefun.jac(u,p,0)
jet = BifurcationKit.getJet(F, J; matrixfree=false)

u0 = solve(ODEProblem(sys,rand(length(sys.states)),(0.,100.0),params_input),Rosenbrock23()).u[end];

# Sets options
dsmax = 1e-2
dsmin = 1e-5
ds = sqrt(dsmax*dsmin);
maxSteps = 10000
opts_br = ContinuationPar(pMin = p_span[1], pMax = p_span[2], dsmax = dsmax, dsmin = dsmin, ds=ds, maxSteps=maxSteps,
detectBifurcation = 3, nInversion = 6, maxBisectionSteps = 25,nev = 4);

# Calculates the diagram.
bif = bifurcationdiagram(jet..., u0, params_input, (@lens _[2]), 2,
    (x,p,level)->setproperties(opts_br);
    tangentAlgo = SecantPred(),
    recordFromSolution=(x, p) -> x[1], verbosity = 0, plot=false);

# Plots the diagram.
plot(bif)

using Catalyst, LinearAlgebra

rn = @reaction_network begin
    k1,  X2 → X1
    k2,  X1 → X2
    k3,  X3 → X5
    k4,  X5 → X4
    k5,  X4 → X6
    k6,  X4 → X3
    k7,  X6 → X5
    k8,  X1 + X5 → X2 + X3
    k9,  X1 + X6 → X2 + X4
end k1 k2 k3 k4 k5 k6 k7 k8 k9
Catalyst.reorder_states!(rn, [2,1,3,5,4,6])

Xs = [1.6563516576558746, 3.5232975392793924e8, 1.0, 0.0038843824790898043, 3.7708607296986027e-10, 0.08021161610521155]
Ks = [2.8384448207706775e-9, 2.276606366932523e-10, 1.0000091094346868, 29473.54319518398, 0.01816749288551057, 0.002345143593835089, 2.499174147678733e-5, 1.6010576150026472e9, 0.0005160737862267585]

N = netstoichmat(rn)
W = conservationlaws(N)

## Singular system
odesys = convert(ODESystem, rn; combinatoric_ratelaws=false)
odefun = ODEFunction(odesys; jac=true)

## Smooth system
# Conservative laws
states_new = states(rn)
Ts = W * Xs
states_new[2] = Ts[1] - states(rn)[1]
states_new[3] = Ts[2] - sum(states(rn)[4:6])

# Build the new system
odefun_new = odefun(states_new, parameters(rn), odesys.iv)
D = Differential(odesys.iv)
eqs_new = equations(odesys)

ind_new = [1,4,5,6]
eqs_new = [D(states(rn)[i]) ~ odefun_new[i] for i in ind_new]

@named odesys_smooth = ODESystem(
    eqs_new,
    odesys.iv,
    [states(rn)[i] for i in ind_new],
    parameters(rn))

odefun_smooth = ODEFunction(odesys_smooth; jac=true)
F_smooth(u, p) = odefun_smooth(u, p, 0.)
J_smooth(u, p) = odefun_smooth.jac(u, p, 0.)
## Sanity check
Xs_smooth = [Xs[i] for i in [1,4,5,6]]
F_smooth(Xs_smooth, Ks)
det(J_smooth(Xs_smooth, Ks))

p_idx = 9                     # The index of the bifurcation parameter.
p_span = (0.000, 1.)   # The parameter range for the bifurcation diagram.
plot_var_idx = 4          # The index of the variable we plot in the bifurcation diagram.
# I know about the option `bothside=true`, but I am not sure how to use it, so first I track the solution `Xs` along `k9` using [HomotopyContinuation.jl.](https://github.com/JuliaHomotopyContinuation/HomotopyContinuation.jl)
import HomotopyContinuation
const HC = HomotopyContinuation

# Translate vars and params to HC
hc_vars = HC.Variable.(Symbol.([states(rn)[i] for i in ind_new]))
hc_params = Vector{Any}(undef, length(Ks))
hc_param = HC.Variable(Symbol(Catalyst.parameters(rn)[p_idx]))

# Set all parameters except `k9`, for which we want to track.
hc_params .= Ks
hc_params[p_idx] = hc_param

# Build HC.system with one parameter
hc_eqs_smooth = F_smooth(hc_vars, hc_params)
hc_sys_smooth = HC.System(hc_eqs_smooth; variables=hc_vars, parameters=[hc_param])

# Track `Xs` along `k9`
x₀ = HC.solve(hc_sys_smooth, Xs_smooth; start_parameters=[Ks[p_idx]], target_parameters=[p_span[1]])
Xs_init_smooth = HC.real_solutions(x₀)[1]
# We obtain
Xs_init_smooth = [7.774317758843825, 0.00011544400492545825, 8.034530803010681e-11, 0.08392084801717949]

using BifurcationKit, Setfield, Plots
const BK = BifurcationKit

cont_opt =
    ContinuationPar(
        dsmax=1e-2, # Maximum arclength value of the pseudo-arc length continuation method.
        dsmin=1e-7, # Minimum arclength
        ds=1e-3, # Initial arclength
        theta=0.01,
        pMin=p_span[1],    # Minimum p-vale (if hit, the method stops).
        pMax=p_span[2],    # Maximum p-vale (if hit, the method stops).
        maxSteps=1000000,   # The maximum number of steps.
        # Value in {0,1,2,3} determining to what extent bifurcation points are detected
        # (0 nothing, 3 both them and their localization are detected).
        detectBifurcation=3,
        nInversion=6,
        # saveToFile        = true,
        # Parameters to the newton solver (when finding fixed points).
        newtonOptions=NewtonPar(tol=1e-9,
                                verbose=false,
                                maxIter=20,
                                )
    )

params_input = setindex!(copy(Ks), p_span[1], p_idx) # The input parameter values have to start at the first index of our parameter span.
## Sanity check
F_smooth(Xs_init_smooth, params_input)

branches, _ =
    continuation(F_smooth, # The vector field
                 J_smooth, # The Jacobian
                 Xs_init_smooth,
                 params_input,
                 (@lens _[p_idx]),
                 cont_opt;
                 plot=false, # We do not want to display, or plot, intermediary results.
                 recordFromSolution=(x, p) -> (x4 = x[2]), # 2 is the new index for X4
                 verbosity=0,
                 callbackN=(x, f, J, res, iteration, itlinear, options; kwargs...) -> res < 1e7,  # Parameter for the continuation method, see BifurcationKit documentation.
                 # filename = "branch_k9_X4"
                 )
#  1,  hopf at p ≈ +0.00051615 ∈ (+0.00051601, +0.00051615), |δp|=1e-07, [converged], δ = ( 2,  2), step =  35, eigenelements in eig[ 36], ind_ev =   2
#  2,  hopf at p ≈ +0.07592376 ∈ (+0.07592131, +0.07592376), |δp|=2e-06, [converged], δ = (-2, -2), step =  48, eigenelements in eig[ 49], ind_ev =   2

# Newton parameters for periodic orbits
newton_opt_po = NewtonPar(tol=1e-9,
                          maxIter=25)

# continuation parameters
cont_opt_po = ContinuationPar(
    dsmax=0.1,
    ds=-0.001,
    dsmin=1e-4,
    newtonOptions=(@set newton_opt_po.tol = 1e-8),
    precisionStability=1e-2,
    detectBifurcation=1,
    saveSolEveryStep=1,
)

# define the sup norm
norminf(x) = norm(x, Inf)
Mt = 101 # number of time sections
jet_smooth = BK.getJet(F_smooth, J_smooth; matrixfree=false)

branch_po, utrap =
    continuation(
        jet_smooth...,
        branches,
        1,  # Periodic orbits for the first Hopf point.
        cont_opt,
        PeriodicOrbitTrapProblem(M=Mt);
        # help branching from Hopf
        # usedeflation=true,
        usedeflation=false,
        # algorithm to solve linear associated with periodic orbit problem
        # linearPO=:FullSparseInplace,
        linearPO=:Dense,
        # tangent algorithm along the branch
        tangentAlgo=BorderedPred(),
        recordFromSolution = (x, p) -> (xtt=reshape(x[1:end-1],2,Mt);
                                        return (max = maximum(xtt[1,:]),
                                                min = minimum(xtt[1,:]),
                                                period = x[end])),
        # recordFromSolution=(x, p) -> (x4 = x[2]),
        callbackN=(x, f, J, res, iteration, itlinear, options; kwargs...) -> res < 1e7,  # Parameter for
        # finaliseSolution=(z, tau, step, contResult; prob=nothing, kwargs...) -> begin
        #     # limit the period
        #     z.u[end] < 100
        # end,
        normC=norminf
    )
# ERROR: MethodError: no method matching create_array(::Type{Base.ReshapedArray{

import SymbolicUtils
@inline function SymbolicUtils.Code.create_array(A::Type{<:Base.ReshapedArray{T,N,P,MI}}, S, nd::Val, d::Val, elems...) where {T,N,P,MI}
    SymbolicUtils.Code.create_array(P, S, nd, d, elems...)
end

## Refinement of Hopf points
outfold, hist, flag =
    newton(F_smooth,
           J_smooth,
           branches, 1; # Refine first point.
           issymmetric=false, Jᵗ=nothing, d2F=nothing,
           normN=norminf, options=branches.contparams.newtonOptions,
           # bdlinsolver = BorderingBLS(options.linsolver),
           startWithEigen=false
           )

# Fetch Packages
using DifferentialEquations, Catalyst
using BifurcationKit, Setfield, LinearAlgebra

# Declare Model
rn = @reaction_network begin
    kDeg,       (w,w2,w2v,v,w2v2,vP,σB,w2σB) ⟶ ∅
    kDeg,       vPp ⟶ phos
    (kBw,kDw),  2w ⟷ w2
    (kB1,kD1),  w2 + v ⟷ w2v
    (kB2,kD2),  w2v + v ⟷ w2v2
    kK1,        w2v ⟶ w2 + vP
    kK2,        w2v2 ⟶ w2v + vP
    (kB3,kD3),  w2 + σB ⟷ w2σB
    (kB4,kD4),  w2σB + v ⟷ w2v + σB
    (kB5,kD5),  vP + phos ⟷ vPp
    kP,         vPp ⟶ v + phos
    v0*((1+F*σB)/(K+σB)),     ∅ ⟶ σB
    λW*v0*((1+F*σB)/(K+σB)),  ∅ ⟶ w
    λV*v0*((1+F*σB)/(K+σB)),  ∅ ⟶ v
    (pStress,1), 0 <--> phos
end kBw kDw kD kB1 kB2 kB3 kB4 kB5 kD1 kD2 kD3 kD4 kD5 kK1 kK2 kP kDeg v0 F K λW λV pInit pStress;

# Set Parameter Values
kBw = 3600; kDw = 18; kD = 18;
kB1 = 3600; kB2 = 3600; kB3 = 3600; kB4 = 1800; kB5 = 3600;
kD1 = 18; kD2 = 18; kD3 = 18; kD4 = 1800; kD5 = 18;
kK1 = 36; kK2 = 12; #kK2 reduced to a third to achive arcitability. kK2 = 36 in the original model.
kP = 180; kDeg = 0.7;
v0 = 0.4; F = 30; K = 0.2;
λW = 4; λV = 4.5;

η = 0.05
pInit = 0.001
pStress = 0.4
parameters = [kBw, kDw, kD, kB1, kB2, kB3, kB4, kB5, kD1, kD2, kD3, kD4, kD5, kK1, kK2, kP, kDeg, v0, F, K, λW, λV, pInit, pStress, η];

# Find good u0 (know good initial guess should be needed, but had some other error when I just used a zero vector)
oprob = ODEProblem(rn,ones(length(rn.states)),(0.,50.),parameters)
sol = solve(oprob,Rosenbrock23());

# Prepare Bifurcation Comoutation
#rn = noise_modulation_model.system
odefun = ODEFunction(convert(ODESystem,rn),jac=true)
F = (u,p) -> odefun(u,p,0)      
J = (u,p) -> odefun.jac(u,p,0);
params = parameters;
p_idx = 15       
p_span = (12.,36.)
plot_var_idx = 7;

# Set Bifurcation Options
opts = ContinuationPar( dsmax = 0.05,        # Maximum arclength value of the pseudo-arc length continuation method.
                        dsmin = 1e-4,        # Minimum arclength value of the pseudo-arc length continuation method.
                        ds=0.001,            # Initial arclength value of the pseudo-arc length continuation method (should be positive).
                        maxSteps = 100000,   # The maximum number of steps.
                        pMin = p_span[1],    # Minimum p-vale (if hit, the method stops).
                        pMax = p_span[2],    # Maximum p-vale (if hit, the method stops).
                        detectBifurcation=3, # Value in {0,1,2,3} determening to what extent bofurcation points are detected (0 means nothing is done, 3 both them and there localisation are detected).
                        newtonOptions = NewtonPar(tol = 1e-9, verbose = false, maxIter = 15)); #Parameters to the newton solver (when finding fixed points) see BifurcationKit documentation.

DO = DeflationOperator( 2.0,    # Algorithm parameter required when using deflated continuation, see BifurcationKit documentation.
                        dot,    # Algorithm parameter required when using deflated continuation, see BifurcationKit documentation.
                        1.,     # Algorithm parameter required when using deflated continuation, see BifurcationKit documentation.
                        [sol.u[end]]); # Guess(es) of the fixed point for the initial parameter set. Do not need to be exact.

# Calculate the diagram
params_input = setindex!(copy(params),p_span[1],p_idx);
@time branches, = continuation(F, J, params_input, (@lens _[p_idx]) ,opts , DO,       # Gives our input.
    verbosity = 0, showplot=false,                                                    # We do not want to display, or plot, intermediary results.
    printSolution=(x, p) -> x[plot_var_idx],                                          # How we wish to print the output in the diagram. Here we simply want the value of the target varriable.
    perturbSolution = (x,p,id) -> (x  .+ 0.8 .* rand(length(x))),                     # Parameter for the continuation method, see BifurcationKit documentation.
    callbackN = (x, f, J, res, iteration, itlinear, options; kwargs...) -> res <1e7)  # Parameter for the continuation method, see BifurcationKit documentation.

using PseudoArcLengthContinuation, LinearAlgebra, Plots, PyPlot, Optim
const PALC = PseudoArcLengthContinuation
using ForwardDiff
using Parameters, Setfield

function f(x, p)
    [p.α*x[1] + x[1]^3 ]
end
function Jf(x, p)
    Jf = zeros(1, 1)
    Jf[1] = p.α+3x[1]^2
    Jf
end

Fmit = f
D(f, x, p, dx) = ForwardDiff.derivative(t->f(x .+ t .* dx, p), 0.)
d1Fmit(x,p,dx1) = D((z, p0) -> Fmit(z, p0), x, p, dx1)
d2Fmit(x,p,dx1,dx2) = D((z, p0) -> d1Fmit(z, p0, dx1), x, p, dx2)
d3Fmit(x,p,dx1,dx2,dx3) = D((z, p0) -> d2Fmit(z, p0, dx1, dx2), x, p, dx3)
jet = (f, d1Fmit, d2Fmit, d3Fmit)
jet = (f, Jf, d2Fmit, d3Fmit)

alpha_min = -10.
alpha_max = 10.
alpha_start = 2.

# options for Newton solver
opt_newton = PALC.NewtonPar(tol = 1e-8, verbose = true, maxIter = 20)

# options for continuation
opts_br = ContinuationPar(dsmin = 0.001, dsmax = 0.05, ds = -0.01, pMax = alpha_max, pMin = alpha_min,
	                  detectBifurcation = 2, nev = 30, plotEveryNsteps = 10, newtonOptions = opt_newton,
	                  maxSteps = 100, precisionStability = 1e-6, nInversion = 4, dsminBisection = 1e-7, maxBisectionSteps = 25)

sol_start, _, _ = newton( f, [1.], (α=.1,), opt_newton)
par_mit = (α=alpha_start,)

br, _ = @time PALC.continuation(
    Fmit, Jf,
    sol_start, par_mit, (@lens _.α), opts_br;
    # printSolution = (x, p) -> norm(x),
    # plotSolution = (x, p; kwargs...) -> plotsol!(x ; kwargs...),
    plot = false, verbosity = 3)

br2, _ = continuation(jet..., br, 1, setproperties(opts_br, ds = -0.01, maxSteps = 20), plot = false, printSolution = (x,p) -> x[1], verbosity=3)

Plots.plot([br,br2])

using PseudoArcLengthContinuation, LinearAlgebra, Plots
const PALC = PseudoArcLengthContinuation
using ForwardDiff
using Parameters, Setfield

function f(x, p)
    [p.α*x[1] + x[1]^3 ]
end
function Jf(x, p)
    Jf = zeros(1, 1)
    Jf[1] = p.α+3x[1]^2
    dx -> Jf*dx
end
Fmit = f

# eigensolver
struct MyStruct <: PALC.AbstractEigenSolver
end
function (m::MyStruct)(J, nev)
    return [(J(1))[1]], [1.], true, 1
end

function D(f, x, p, dx)
    ε=1e-7
    (f(x + ε*dx, p) - f(x - ε*dx, p))/(2ε)
end
d1Fmit(x,p,dx1) = D((z, p0) -> Fmit(z, p0), x, p, dx1)
d2Fmit(x,p,dx1,dx2) = D((z, p0) -> d1Fmit(z, p0, dx1), x, p, dx2)
d3Fmit(x,p,dx1,dx2,dx3) = D((z, p0) -> d2Fmit(z, p0, dx1, dx2), x, p, dx3)
jet = (f, d1Fmit, d2Fmit, d3Fmit)
jet = (f, Jf, d2Fmit, d3Fmit)

alpha_min = -10.
alpha_max = 10.
alpha_start = 2.

# options for Newton solver
opt_newton = PALC.NewtonPar(tol = 1e-8, verbose = true, maxIter = 20, eigsolver = MyStruct())
opt_newton = PALC.NewtonPar(tol = 1e-8, verbose = true, maxIter = 20, eigsolver = MyStruct(), linsolver=GMRESKrylovKit())

# options for continuation
opts_br = ContinuationPar(dsmin = 0.001, dsmax = 0.05, ds = -0.01, pMax = alpha_max, pMin = alpha_min,
	                  detectBifurcation = 2, nev = 30, plotEveryNsteps = 10, newtonOptions = opt_newton,
	                  maxSteps = 100, precisionStability = 1e-6, nInversion = 4, dsminBisection = 1e-7, maxBisectionSteps = 25)

sol_start, _, _ = newton( f, [1.], (α=.1,), opt_newton)
par_mit = (α=alpha_start,)

br, _ = @time PALC.continuation(
    Fmit, Jf,
    sol_start, par_mit, (@lens _.α), opts_br;
    # printSolution = (x, p) -> norm(x),
    # plotSolution = (x, p; kwargs...) -> plotsol!(x ; kwargs...),
    plot = false, verbosity = 3)

br2, _ = continuation(jet..., br, 1, setproperties(opts_br, ds = -0.01, maxSteps = 20), plot = false, printSolution = (x,p) -> x[1], verbosity=3, Jt=Jf)
br3, _ = continuation(jet..., br, 1, setproperties(opts_br, ds = 0.01, maxSteps = 20), plot = false, printSolution = (x,p) -> x[1], verbosity=3, Jt=Jf)

Plots.plot([br,br2,br3], marker=2)

function Fitz(u,p)
    x,y = u
    @unpack I,a,b,c = p
    du = similar(u)
    
    du[1] = x-x^3/3-y + I
    du[2] = a*(b*x-c*y)
    
    return du
end

sol0 = [0.1,0.1]

par_Fitz = (I=-2.0,a=0.08,b=1.,c=0.8)

ls = GMRESKrylovKit(dim = 100)
# defaultLS() not defined

optnewton = NewtonPar(tol = 1e-11, verbose = true, linsolver = ls)
out, _, _ = @time newton( Fitz, sol0, par_Fitz, optnewton)

optnew = NewtonPar(verbose = true, tol = 1e-10)
optcont = ContinuationPar(dsmin = 0.01, dsmax = 0.1, ds= 0.05, pMax = 5.1, pMin = -2.1,
    newtonOptions = setproperties(optnew; maxIter = 50, tol = 1e-10), 
	maxSteps = 300, plotEveryStep = 40, 
	detectBifurcation = 3, nInversion = 4, tolBisectionEigenvalue = 1e-12, dsminBisection = 1e-5)

br, _ = @time continuation(Fitz, out, par_Fitz, (@lens _.I),
		optcont; plot = true, verbosity = 0,
		printSolution = (x, p) -> x[1])


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

# index of the Hopf point in br.bifpoint
norminf = x -> norm(x, Inf)
ind_hopf = 1
hopfpoint, _, flag = @time newton(Fitz, dFbpSecBif,
	br, ind_hopf, par_Fitz, (@lens _.I); normN = norminf)
flag && printstyled(color=:red, "--> We found a Hopf Point at μ = ", hopfpoint.p[1],
    ", from μ = ", br.bifpoint[ind_hopf].param, "\n")

# automatic branch switching from Hopf point
opt_po = NewtonPar(tol = 1e-5, verbose = true, maxIter = 15)
opts_po_cont = ContinuationPar(dsmin = 0.01, dsmax = 0.1, ds = 0.1, pMax = 5.5,pMin=-2.1, maxSteps = 200,
    newtonOptions = opt_po, saveSolEveryStep = 2, plotEveryStep = 1, nev = 11, precisionStability = 1e-4,
	detectBifurcation = 3, dsminBisection = 1e-4, maxBisectionSteps = 10);

# number of time slices for the periodic orbit
M = 150
probFD = PeriodicOrbitTrapProblem(M = M)
br_po, = continuation(
	# arguments for branch switching from the first 
	# Hopf bifurcation point
	jet..., br, 1,
	# arguments for continuation
	opts_po_cont, probFD;
	# OPTIONAL parameters
	# we want to jump on the new branch at phopf + δp
	# ampfactor is a factor to increase the amplitude of the guess
	δp = 0.05, ampfactor = 1.1,
	# specific method for solving linear system
	# of Periodic orbits with trapeze method
	linearPO = :FullLU,
# regular options for continuation
	verbosity = 3,	plot = true, 
	printSolution = (x,p)->x[1], normC = norminf)

fieldnames(typeof(br_po))
plot(br)
plot!(br_po)

parameters = ContinuationPar() # inital hyperparameters

for i=1:10 # some optimisation loop
       parameters.ds *= i*parameters.ds # for example changing ds
end

using PseudoArcLengthContinuation, LinearAlgebra, Plots
const PALC = PseudoArcLengthContinuation
using ForwardDiff
using Parameters, Setfield

function f(x, p)
    [p.α*x[1] + x[1]^3 ]
end
function Jf(x, p)
    Jf = zeros(1, 1)
    Jf[1] = p.α+3x[1]^2
    dx -> Jf*dx
end

# eigensolver
struct MyStruct <: PALC.AbstractEigenSolver
end
function (m::MyStruct)(J, nev)
    return [(J(1))[1]], [1.], true, 1
end

Fmit = f
D(f, x, p, dx) = ForwardDiff.derivative(t->f(x .+ t .* dx, p), 0.)
function D(f, x, p, dx)
    ε=1e-10
    (f(x + ε*dx, p) - f(x, p))/ε
end
d1Fmit(x,p,dx1) = D((z, p0) -> Fmit(z, p0), x, p, dx1)
d2Fmit(x,p,dx1,dx2) = D((z, p0) -> d1Fmit(z, p0, dx1), x, p, dx2)
d3Fmit(x,p,dx1,dx2,dx3) = D((z, p0) -> d2Fmit(z, p0, dx1, dx2), x, p, dx3)
jet = (f, d1Fmit, d2Fmit, d3Fmit)
jet = (f, Jf, d2Fmit, d3Fmit)

alpha_min = -10.
alpha_max = 10.
alpha_start = 2.

# options for Newton solver
opt_newton = PALC.NewtonPar(tol = 1e-8, verbose = true, maxIter = 20, eigsolver = MyStruct())
# opt_newton = PALC.NewtonPar(tol = 1e-8, verbose = true, maxIter = 20, eigsolver = MyStruct(), linsolver=GMRESKrylovKit())

# options for continuation
opts_br = ContinuationPar(dsmin = 0.001, dsmax = 0.05, ds = -0.01, pMax = alpha_max, pMin = alpha_min,
	                  detectBifurcation = 2, nev = 30, plotEveryNsteps = 10, newtonOptions = opt_newton,
	                  maxSteps = 100, precisionStability = 1e-6, nInversion = 4, dsminBisection = 1e-7, maxBisectionSteps = 25)

sol_start, _, _ = newton( f, [1.], (α=.1,), opt_newton)
par_mit = (α=alpha_start,)

br, _ = @time PALC.continuation(
    Fmit, Jf,
    sol_start, par_mit, (@lens _.α), opts_br;
    # printSolution = (x, p) -> norm(x),
    # plotSolution = (x, p; kwargs...) -> plotsol!(x ; kwargs...),
    plot = false, verbosity = 3)

br2, _ = continuation(jet..., br, 1, setproperties(opts_br, ds = -0.01, maxSteps = 20), plot = false, printSolution = (x,p) -> x[1], verbosity=3, Jt=Jf)
br3, _ = continuation(jet..., br, 1, setproperties(opts_br, ds = 0.01, maxSteps = 20), plot = false, printSolution = (x,p) -> x[1], verbosity=3, Jt=Jf)

Plots.plot([br,br2,br3])

using Revise, Parameters, Setfield, Plots, LinearAlgebra
using DifferentialEquations
using BifurcationKit
const BK = BifurcationKit

# define the sup norm
norminf(x) = norm(x, Inf)

# function to record information from a solution
recordFromSolution(x, p) = (u1 = x[1], u2 = x[2])

####################################################################################################

function F!(f, x, p, t=0)
	@unpack a, b1, b2, d, λ, ω = p
	f[1] = x[2]
	f[2] = - d * x[2] - a * x[1] - b1 * (x[1]^2) - b2 * (x[1]^3) + λ * sin(ω * t)
    return f
end
F(x,p,t) = F!(similar(x),x,p,t)

# Defining the Jacobian
function dF(x,p,t = 0)
	@unpack a, b1, b2, d, λ, ω = p
    M = zeros(2,2)
    M[1,1] = 0
    M[1,2] = 1
    M[2,1] = -a -b1*2*x[1] -b2*3*(x[1]^2)
    M[2,2] = -d
    return M
 end

# parameters of the model
P = @with_kw (a=-1., b1=1., b2=0., d=0.1, λ=0.06, ω=0.81)
par_Helm = P(-1.0, 1.0, 0.0, 0.1, 0.06, 0.81)
par0 = (@set par_Helm.a = -1.15)

# initial condition guess
z0 = [1.; 0.001]
vf = ODEFunction(F; jac = dF)       #ODEfunction with the Jacobian
prob = ODEProblem(vf, z0, (0.0, 1000*2*π/par0.ω), par0)   #ODEProblem encoding
sol = solve(prob, Tsit5())

z0 = sol.u[end]
prob = ODEProblem(vf, z0, (0.0, 2*π/par0.ω), par0)   #ODEProblem encoding
sol = solve(prob, Tsit5())

# this encodes the functional for the Shooting problem
probSh = ShootingProblem(
	# we pass the ODEProblem encoding the flow and the time stepper
	prob, Tsit5(),
	# this is for the phase condition, you can pass your own section as well
	[vec(sol[end])],
	# enable threading
	parallel = true,
	# these are options passed to the ODE time stepper
	atol = 1e-10, rtol = 1e-8)

# Newton interation for low energy level
zp0 = vcat(vec(sol[end]),2*π/par_Helm.ω)

ls = GMRESIterativeSolvers(reltol = 1e-7, N = length(zp0), maxiter = 100)
newtonOpt = NewtonPar(verbose = true, tol = 1e-8, maxIter = 25, linsolver = ls)

# continuation options
opts_cont = ContinuationPar(dsmin = 0.001, dsmax = 1.0, ds= 0.01, pMin = 0.02, pMax = 0.1,
	maxSteps = 50, newtonOptions = newtonOpt, nev = 2,
	# options to detect bifurcations
	detectBifurcation = 3,   # this gives an error!
        detectFold = true, nInversion = 8, maxBisectionSteps = 25,
	precisionStability = 1e-8, saveSolEveryStep = 1)

br, = continuation(
    probSh,
    zp0,
    par0,
    (@lens _.λ),
    opts_cont,
    linearAlgo = MatrixFreeBLS(@set ls.N = ls.N + 1);
    verbosity = 3,
    plot = true,
    recordFromSolution = recordFromSolution,
	normC = norminf
)

# Fetch packages.
using Catalyst
using BifurcationKit, LinearAlgebra, Setfield, SparseArrays, Plots

# Create functions.
rn = @reaction_network begin
    p, 0 --> X
    d, X --> 0
end p d;
odefun = ODEFunction(convert(ODESystem,rn),jac=true)
F = (u,p) -> odefun(u,p,0)
J = (u,p) -> odefun.jac(u,p,0);

# Set initial parameter values.
p = (p=1.,d=0.1);
# Sets initial guess
x0 = 0.5;

# Set solving options
opts = BifurcationKit.ContinuationPar(dsmax = 0.051, dsmin = 1e-3, ds=0.001, maxSteps = 140, pMin = -3., saveSolEveryStep = 0, newtonOptions = NewtonPar(tol = 1e-8, verbose = false), saveEigenvectors = false);

# Set Deflation Options
DO = DeflationOperator(2.0, dot, .001, [[0.]]);

# Create bifurcation diagram
brdc, = continuation(F,J, x0, (@lens _),opts,DO; showplot=false, verbosity = 0);

# Plot the diagram
plot(brdc[1])

brdc, = continuation(F,J,p, (@lens _.p),opts,DO; showplot=false, verbosity = 0);

brdc, = continuation(F,J,x0, p, (@lens _.p),opts,DO; showplot=false, verbosity = 0);
brdc, = continuation(F,J,p, x0, (@lens _.p),opts,DO; showplot=false, verbosity = 0);

# Fetch packages.
using Catalyst
using BifurcationKit, LinearAlgebra, Setfield, SparseArrays, Plots

# Create functions.
rn2 = @reaction_network begin
    (v0+(S*X)^n/((S*X)^n+(D*A)^n+1),1), ∅ ↔ X
    (τ*X,τ), ∅ ↔ A
end v0 S n D τ;
odefun2 = ODEFunction(convert(ODESystem,rn2),jac=true)
F2 = (u,p) -> odefun2(u,p,0)
J2 = (u,p) -> odefun2.jac(u,p,0);

# Set initial parameter values.
p2 = (v0=0.1, S=1., n=3, D=2., τ=0.01);
# Sets initial guess
x02 = [0.5, 0.5];

# Set solving options
opts2 = BifurcationKit.ContinuationPar(dsmax = 0.051, dsmin = 1e-3, ds=0.001, maxSteps = 1000000, pMin = 0.1, pMax = 5., saveSolEveryStep = 0, newtonOptions = NewtonPar(tol = 1e-8, verbose = false), saveEigenvectors = false);

# Set Deflation Options
DO2 = DeflationOperator(2.0, dot, .001, [[0.,0.]]);

# Create bifurcation diagram
brdc2, = continuation(F2,J2,p2, (@lens _.S),opts2,DO2; showplot=false, verbosity = 0);

# Plot the diagram
plot(brdc2[1])

termination_condition(x, f, J, residual, iteration, niter, options; kwargs...) = residual <1e3
root_search(x,p,id) = x .+ rand(length(x))

# Fetch packages.
using BifurcationKit
using Catalyst
using OrdinaryDiffEq
using Plots
using Setfield

# Declare our model.
hillAR(x1,x2,v,K,n) = (v*x1^n)/(x1^n+x2^n+K^n)
ss_network = @reaction_network begin
    v0+hillAR(S*X,D*Y,v,K,n), ∅ → X
    v0+hillAR(S*Y,D*X,v,K,n), ∅ → Y
    1,(X,Y) → ∅
end S D v0 v n K;

# Generate system functions.
odefun = ODEFunction(convert(ODESystem,rn),jac=true)
F = (u,p) -> odefun(u,p,0)      
J = (u,p) -> odefun.jac(u,p,0);
jet = BifurcationKit.getJet(F, J; matrixfree=false)

# Select parameters
params_in = [1.0,0.5,0.1,1.0,2,0.5]; 
p_idx = 1; 
p_span = (0.1,2.5)
params = setindex!(copy(params_in),p_span[1],p_idx);

# Set options.
dsmax = params[p_idx]/500
dsmin = min(params[p_idx]/10000,1e-4)
ds = sqrt(dsmax*dsmin);
opts_br = ContinuationPar(pMin = p_span[1], pMax = p_span[2], dsmax = dsmax, dsmin = dsmin, ds=ds, maxSteps = 1000000,
    detectBifurcation = 3,
    nInversion = 6, maxBisectionSteps = 25,
    nev = 3);

# Compute diagram.
bif = bifurcationdiagram(jet..., u0, params, (@lens _[p_idx]), 3,
    (x,p,level)->setproperties(opts_br); 
    tangentAlgo = NaturalPred(),
    recordFromSolution=(x, p) -> x[v_idx], verbosity = 0, plot=false);

# Plot the diagram.
plot(bif)

plot(plot(bif.γ),plot(bif.child[1].γ),plot(bif.child[2].γ),plot(bif.child[3].γ),plot(bif.child[4].γ),size=(1200,600))

(x,p,level)->setproperties(opts_br;ds=-ds*(-1)^level)

F = (x, p) -> F(x, p)

using Revise, ForwardDiff, Parameters, Setfield, Plots, LinearAlgebra
using BifurcationKit
const BK = BifurcationKit

# sup norm
norminf(x) = norm(x, Inf)

# vector field
function TMvf!(dz, z, p, t)
	@unpack J, α, E0, τ, τD, τF, U0 = p
	E, x, u = z
	SS0 = J * u * x * E + E0
	SS1 = α * log(1 + exp(SS0 / α))
	dz[1] = (-E + SS1) / τ
	dz[2] =	(1.0 - x) / τD - u * x * E
	dz[3] = (U0 - u) / τF +  U0 * (1.0 - u) * E
	dz
end

# out of place method
TMvf(z, p) = TMvf!(similar(z), z, p, 0)

# we group the differentials together
dTMvf(z,p) = ForwardDiff.jacobian(x -> TMvf(x,p), z)
jet = BK.getJet(TMvf, dTMvf)

# parameter values
par_tm = (α = 1.5, τ = 0.013, J = 3.07, E0 = -2.0, τD = 0.200, U0 = 0.3, τF = 1.5, τS = 0.007)

# initial condition
z0 = [0.238616, 0.982747, 0.367876]

# continuation options
opts_br = ContinuationPar(pMin = -10.0, pMax = -0.9,
	# parameters to have a smooth result
	ds = 0.04, dsmax = 0.05,
	# this is to detect bifurcation points precisely with bisection
	detectBifurcation = 3,
	# Optional: bisection options for locating bifurcations
	nInversion = 8, maxBisectionSteps = 25, nev = 3)

# continuation of equilibria
br, = continuation(TMvf, dTMvf, z0, par_tm, (@lens _.E0), opts_br;
	recordFromSolution = (x, p) -> (E = x[1], x = x[2], u = x[3]),
	tangentAlgo = BorderedPred(),
	plot = true, normC = norminf)

scene = plot(br, plotfold=false, markersize=3, legend=:topleft)

br

# newton parameters
optn_po = NewtonPar(verbose = true, tol = 1e-8,  maxIter = 10)

# continuation parameters
opts_po_cont = ContinuationPar(dsmax = 0.1, ds= -0.0001, dsmin = 1e-4, pMax = 0., pMin=-5.,
	maxSteps = 110, newtonOptions = (@set optn_po.tol = 1e-7),
	nev = 3, precisionStability = 1e-8, detectBifurcation = 3, plotEveryStep = 10, saveSolEveryStep=1)

# arguments for periodic orbits
args_po = (	recordFromSolution = (x, p) -> (xtt = BK.getTrajectory(p.prob, x, p.p);
		return (max = maximum(xtt.u[1,:]),
				min = minimum(xtt.u[1,:]),
				period = x[end])),
	plotSolution = (x, p; k...) -> begin
		xtt = BK.getTrajectory(p.prob, x, p.p)
		plot!(xtt.t, xtt.u[1,:]; label = "E", k...)
		plot!(xtt.t, xtt.u[2,:]; label = "x", k...)
		plot!(xtt.t, xtt.u[3,:]; label = "u", k...)
		plot!(br,subplot=1, putbifptlegend = false)
		end,
	normC = norminf)

Mt = 200 # number of time sections
	br_potrap, utrap = continuation(jet...,
	# we want to branch form the 4th bif. point
	br, 4, opts_po_cont,
	# we want to use the Trapeze method to locate PO
	PeriodicOrbitTrapProblem(M = Mt);
	# this jacobian is specific to ODEs
	# it is computed using AD and
	# updated inplace
	linearPO = :Dense,
	# regular continuation options
	verbosity = 2,	plot = true,
	args_po...)

scene = plot(br, br_potrap, markersize = 3)
plot!(scene, br_potrap.param, br_potrap.min, label = "")

BK.getTrajectory

Cont.continuation( .... recordfunction = u -> [ u[1], u[3] ] )

# Fetch Packages
using DifferentialEquations, Catalyst
using BifurcationKit, Setfield, LinearAlgebra

# Declare Model
rn = @reaction_network begin
    kDeg,       (w,w2,w2v,v,w2v2,vP,σB,w2σB) ⟶ ∅
    kDeg,       vPp ⟶ phos
    (kBw,kDw),  2w ⟷ w2
    (kB1,kD1),  w2 + v ⟷ w2v
    (kB2,kD2),  w2v + v ⟷ w2v2
    kK1,        w2v ⟶ w2 + vP
    kK2,        w2v2 ⟶ w2v + vP
    (kB3,kD3),  w2 + σB ⟷ w2σB
    (kB4,kD4),  w2σB + v ⟷ w2v + σB
    (kB5,kD5),  vP + phos ⟷ vPp
    kP,         vPp ⟶ v + phos
    v0*((1+F*σB)/(K+σB)),     ∅ ⟶ σB
    λW*v0*((1+F*σB)/(K+σB)),  ∅ ⟶ w
    λV*v0*((1+F*σB)/(K+σB)),  ∅ ⟶ v
    
    (pStress,1), 0 <--> phos
    1,  vPp ⟶ vP 
end kBw kDw kD kB1 kB2 kB3 kB4 kB5 kD1 kD2 kD3 kD4 kD5 kK1 kK2 kP kDeg v0 F K λW λV pInit pStress;

# Set Parameter Values
kBw = 3600; kDw = 18; kD = 18;
kB1 = 3600; kB2 = 3600; kB3 = 3600; kB4 = 1800; kB5 = 3600;
kD1 = 18; kD2 = 18; kD3 = 18; kD4 = 1800; kD5 = 18;
kK1 = 36; kK2 = 12; #kK2 reduced to a third to achive arcitability. kK2 = 36 in the original model.
kP = 180; kDeg = 0.7;
v0 = 0.4; F = 30; K = 0.2;
λW = 4; λV = 4.5;

η = 0.05
pInit = 0.001
pStress = 0.4
parameters = [kBw, kDw, kD, kB1, kB2, kB3, kB4, kB5, kD1, kD2, kD3, kD4, kD5, kK1, kK2, kP, kDeg, v0, F, K, λW, λV, pInit, pStress, η];

# Find good u0 (know good initial guess should be needed, but had some other error when I just used a zero vector)
oprob = ODEProblem(rn,ones(length(rn.states)),(0.,50.),parameters)
sol = solve(oprob,Rosenbrock23());

# Prepare Bifurcation Comoutation
#rn = noise_modulation_model.system
odefun = ODEFunction(convert(ODESystem,rn),jac=true)
F = (u,p) -> odefun(u,p,0)      
J = (u,p) -> odefun.jac(u,p,0);
params = parameters;
p_idx = 15       
p_span = (12.,36.)
plot_var_idx = 7;

# Set Bifurcation Options
opts = ContinuationPar( dsmax = 0.05,        # Maximum arclength value of the pseudo-arc length continuation method.
                        dsmin = 1e-4,        # Minimum arclength value of the pseudo-arc length continuation method.
                        ds=0.001,            # Initial arclength value of the pseudo-arc length continuation method (should be positive).
                        maxSteps = 100000,   # The maximum number of steps.
                        pMin = p_span[1],    # Minimum p-vale (if hit, the method stops).
                        pMax = p_span[2],    # Maximum p-vale (if hit, the method stops).
                        detectBifurcation=3, # Value in {0,1,2,3} determening to what extent bofurcation points are detected (0 means nothing is done, 3 both them and there localisation are detected).
                        newtonOptions = NewtonPar(tol = 1e-9, verbose = false, maxIter = 15)); #Parameters to the newton solver (when finding fixed points) see BifurcationKit documentation.

DO = DeflationOperator( 2.0,    # Algorithm parameter required when using deflated continuation, see BifurcationKit documentation.
                        dot,    # Algorithm parameter required when using deflated continuation, see BifurcationKit documentation.
                        1.,     # Algorithm parameter required when using deflated continuation, see BifurcationKit documentation.
                        [sol.u[end]]); # Guess(es) of the fixed point for the initial parameter set. Do not need to be exact.

# Calculate the diagram
params_input = setindex!(copy(params),p_span[1],p_idx);
@time branches, = continuation(F, J, params_input, (@lens _[p_idx]) ,opts , DO,       # Gives our input.
           verbosity = 0, plot=false, 
           recordFromSolution=(x, p) -> x[plot_var_idx], 
           perturbSolution = (x,p,id) -> (x  .+ 0.8 .* rand()),
           callbackN = BifurcationKit.cbMaxNorm(1e7));  
``

# 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));)

# Plots the diagram.
plot(diagram)

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));)

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)

using BifurcationKit,Parameters,StaticArrays,ForwardDiff
using BifurcationKit: DeflationOperator
using LinearAlgebra: dot

########################################################### setup
function F(x,params) where T<:Number
	@unpack k = params
	f = first(x)*first(k)
	F = similar(x,typeof(f))

	F[1] = k[2] + x[1] - x[1]^(k[1]+1)/(k[1]+1)
	F[2] = x[1] - x[2]
	return F
end

J(x,p) = ForwardDiff.jacobian( x->F(x,p), x )
p = ( k=[1.0,1.0], )

using LinearAlgebra: lu
import LinearAlgebra: factorize
factorize(A::StaticMatrix) = lu(A) # default factorization for StaticArrays

########################################################### normal arrays
x = [-1.0,42.0]
deflated_root = [Inf,Inf]
deflation = DeflationOperator(1.0, dot, 1.0, [ deflated_root ] )

u, = newton( F, J, x, p, NewtonPar() )                      # success
u, = newton( F, J, x, p, NewtonPar(), deflation )      # success

########################################################### static immutable arrays
x = @SVector[-1.0,42.0]
deflated_root = @SVector[Inf,Inf]
deflation = DeflationOperator(1.0, dot, 1.0, [ deflated_root ] )

u, = newton( F, J, x, p, NewtonPar() )                      # success
u, = newton( F, J, x, p, NewtonPar(), deflation )      # no method matching ( ::DeflationProblem, ::SArray, ... )

########################################################### static mutable arrays
x = @MVector[-1.0,42.0]
deflated_root = @MVector[Inf,Inf]
deflation = DeflationOperator(1.0, dot, 1.0, [ deflated_root ] )

u, = newton( F, J, x, p, NewtonPar() )                      # success
u, = newton( F, J, x, p, NewtonPar(), deflation )      # success

########################################################### sized mutable arrays
x = SizedVector{2}([-1.0,42.0])
deflated_root = SizedVector{2}([Inf,Inf])
deflation = DeflationOperator(1.0, dot, 1.0, [ deflated_root ] )

u, = newton( F, J, x, p, NewtonPar() )                      # success
u, = newton( F, J, x, p, NewtonPar(), deflation )      # success

using Plots, BifurcationKit, Setfield, ForwardDiff, Parameters
const BK = BifurcationKit

function f(x,p)
	@unpack a,g,γ,r,d = p
	M,C = x
	T = 1 - M - C
	dM = @. a * M * C  - (g * M) / (M + T) + γ * M * T
	dC = @. r * T * C - d * C - a * M * C
	return [dM, dC]
end

jet = BK.getJet(f; matrixfree = false)

# Starting variables
start = [0.2, 0.6]
# Default values for parameters
parameters = (a = 0.1, g = 0.01, γ = 0.8, r = 1., d = 0.44)

# Options for Newton-Krylov solver and continuation
optnewton = NewtonPar(tol = 1e-12, maxIter = 100)
options = ContinuationPar(pMin = 0.01, pMax = 1., dsmin = 0.0001, dsmax = 0.05, ds = 0.01,
	maxSteps = 100, nev = 30, nInversion = 8, detectBifurcation = 3, newtonOptions = optnewton,
	dsminBisection = 1e-7, maxBisectionSteps = 100)

# Solving for equilibrium from starting point using Newton-Krylon solver
solution, = newton(f, start, parameters, optnewton)
solution[2]

br, = continuation(f, solution, parameters, (@lens _.g), options, 
	plot = true, 
	recordFromSolution = (x,p) -> (C = x[2], M = x[1])
)

br1, = continuation(jet..., br, 1, setproperties(options;ds = 0.0001, maxSteps = 30), recordFromSolution = (x,p) -> (C = x[2], M = x[1]), bothside = true, plot = true)
br2, = continuation(jet..., br1, 2, setproperties(options; maxSteps = 20, ds = 0.0001), recordFromSolution = (x,p) -> (C = x[2], M = x[1]), bothside = true,plot = true)

plot(br, br1; branchlabel = ["br", "br1"], legend = :topright, ylims = (0,0.8), vars = (:param, :C))
plot(br, br1; branchlabel = ["br", "br1"], legend = :topright, ylims = (0,1), vars = (:param, :M))