catalystexp1.jl

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

using Catalyst
rn = @reaction_network begin
	(v0+v*(S*X)^n/((S*X)^n+(D*A)^n+K^n),d), ∅ ↔ X
	(τ*X,τ), ∅ ↔ A
end S D τ v0 v K n d

odefun = ODEFunction(convert(ODESystem,rn),jac=true)
F = (u,p) -> odefun(Array(u),p,0) # this is because MK does not like views
J = (u,p) -> odefun.jac(Array(u),p,0) # this is because MK does not like views

params = [1.,9.,0.001,0.01,2.,20.,3,0.05]

p_idx = 1			# The index of the bifurcation parameter.
p_span = (0.1,20.)   # The parameter range for the bifurcation diagram.
plot_var_idx = 1	 # The index of the variable we plot in the bifurcation diagram.

opts = ContinuationPar( dsmax = 0.95,		# 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 = 130,   # 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} determining to what extent bifurcation points are detected (0 means nothing is done, 3 both them and there localisation are detected).
						nInversion = 6, nev = 2,
						newtonOptions = NewtonPar(tol = 1e-9, verbose = false, maxIter = 5)) #Parameters to the newton solver (when finding fixed points) see BifurcationKit documentation.

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

@error 1==0 "Ca ne marche pas du tout sur les valeurs propres!!!"

br, = continuation(F, J, fill(0.,length(rn.states)), params_input, (@lens _[p_idx]) ,opts;			 # Gives our input.
	linearAlgo = MatrixBLS(),
	# tangentAlgo = BorderedPred(),
	verbosity = 3, plot=true,													# 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.
	)

####################################################################################################
# computation of the derivatives
using ForwardDiff
D(f, x, p, dx) = ForwardDiff.derivative(t->f(x .+ t .* dx, p), 0.)
d1F(x,p,dx1) = D((z, p0) -> F(z, p0), x, p, dx1)
d2F(x,p,dx1,dx2) = D((z, p0) -> d1F(z, p0, dx1), x, p, dx2)
d3F(x,p,dx1,dx2,dx3) = D((z, p0) -> d2F(z, p0, dx1, dx2), x, p, dx3)
jet  = (F, J, d2F, d3F)

hopfpt = computeNormalForm(jet..., br, 3)

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

# continuation parameters
opts_po_cont = ContinuationPar(dsmax = 0.3, ds= -0.0001, dsmin = 1e-4, pMax = p_span[2], pMin=p_span[1], maxSteps = 1000, newtonOptions = (@set optn_po.tol = 1e-8), nev = 2, precisionStability = 4e-2, detectBifurcation = 0, plotEveryStep = 20, saveSolEveryStep=1)

Mt = 200 # number of sections
	br_po, utrap = continuation(
	jet..., br, 3, opts_po_cont,
	PeriodicOrbitTrapProblem(M = Mt);
	# specify the parameters for the initial guess. Otherwise they are extracted from opts_po_cont
	ampfactor = 1., δp = 0.001,
	# update the section every 1 steps. Very important
	updateSectionEveryStep = 10,
	# tangent predictor
	tangentAlgo = BorderedPred(),
	# how the jacobian of the Periodic orbit problem is computed. Here an optimised version for ODEs
	linearPO = :Dense,
	verbosity = 3,	plot = true,
	printSolution = (x, p) -> (xtt=reshape(x[1:end-1],2,Mt); return (max = maximum(xtt[1,:]), min = minimum(xtt[1,:]), period = x[end])),
	# function to plot the solution along the branch
	plotSolution = (x, p; k...) -> begin
		xtt = BK.getTrajectory(p.prob, x, p.p)
		plot!(xtt.t, xtt.u[1,:]; label = "u1", markersize = 2, k...)
		plot!(xtt.t, xtt.u[2,:]; label = "u2", k...)
		plot!(br, br,subplot=1, putbifptlegend = false)
		end,
	# norm used for the computation
	normC = x->norm(x,Inf)
	)

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