add inplace heurisitics to allow for passing inplace vector fields. U…

…nder the hood, out of place methods are still used.

clean the tests according to the additions of this commit

corrections of inplace problems constructor

improve Makie backend

examples/COModel.jl

@@ -14,7 +14,6 @@ function COm!(du, u, p, t = 0)
 	du[3] = q4 * (z - k * s)
 	du
 end
-COm(u,p) = COm!(similar(u),u,p,0)
 
 par_com = (q1 = 2.5, q2 = 2.0, q3 = 10., q4 = 0.0675, q5 = 1., q6 = 0.1, k = 0.4)
 

examples/SHpde_snaking.jl

@@ -33,7 +33,8 @@ sol0 = 1.1cos.(X) .* exp.(-0X.^2/(2*5^2))
 
 prob = BifurcationProblem(R_SH, sol0, parSH, (@lens _.λ); J = Jac_sp,
 	recordFromSolution = (x, p) -> (n2 = norm(x), nw = normweighted(x), s = sum(x), s2 = x[end ÷ 2], s4 = x[end ÷ 4], s5 = x[end ÷ 5]),
-	plotSolution = (x, p;kwargs...)->(plot!(X, x; ylabel="solution", label="", kwargs...)))
+	plotSolution = (x, p;kwargs...)->(plot!(X, x; ylabel="solution", label="", kwargs...))
+	)
 ####################################################################################################
 optnew = NewtonPar(verbose = false, tol = 1e-12)
 # allocations 357, 0.8ms

examples/TMModel.jl

@@ -13,16 +13,12 @@ function TMvf!(dz, z, p, t = 0)
 	dz
 end
 
-TMvf(z, p) = TMvf!(similar(z), z, p, 0)
-
 par_tm = (α = 1.5, τ = 0.013, J = 3.07, E0 = -2.0, τD = 0.200, U0 = 0.3, τF = 1.5, τS = 0.007) #2.87
 z0 = [0.238616, 0.982747, 0.367876 ]
-prob = BK.BifurcationProblem(TMvf, z0, par_tm, (@lens _.E0); recordFromSolution = (x, p) -> (E = x[1], x = x[2], u = x[3]),)
+prob = BifurcationProblem(TMvf!, z0, par_tm, (@lens _.E0); recordFromSolution = (x, p) -> (E = x[1], x = x[2], u = x[3]),)
 
 opts_br = ContinuationPar(pMin = -10.0, pMax = -0.9, ds = 0.04, dsmax = 0.125, nInversion = 8, detectBifurcation = 3, maxBisectionSteps = 25, nev = 3)
-opts_br = @set opts_br.newtonOptions.verbose = false
-br = continuation(prob, PALC(tangent=Bordered()), opts_br;
-	plot = true, normC = norminf)
+br = continuation(prob, PALC(tangent = Bordered()), opts_br; plot = true, normC = norminf)
 
 plot(br, plotfold=false, markersize=4, legend=:topleft)
 ####################################################################################################

examples/chan.jl

@@ -36,7 +36,7 @@ n = 101
 par = (α = 3.3, β = 0.01)
 sol0 = [(i-1)*(n-i)/n^2+0.1 for i=1:n]
 
-prob = BK.BifurcationProblem(F_chan, sol0, par, (@lens _.α); plotSolution = (x, p; kwargs...) -> (plot!(x;ylabel="solution",label="", kwargs...)))
+prob = BifurcationProblem(F_chan, sol0, par, (@lens _.α); plotSolution = (x, p; kwargs...) -> (plot!(x;ylabel="solution",label="", kwargs...)))
 
 optnewton = NewtonPar(tol = 1e-8, verbose = true)
 # ca fait dans les 63.59k Allocations

examples/codim2PO.jl

@@ -14,13 +14,13 @@ function Pop!(du, X, p, t = 0)
 	du[4] = 2pi * u + v - s * v
 	du
 end
-Pop(u,p) = Pop!(similar(u),u,p,0)
+# Pop(u,p) = Pop!(similar(u),u,p,0)
 
 par_pop = ( K = 1., r = 2π, a = 4π, b0 = 0.25, e = 1., d = 2π, ϵ = 0.2, )
 
 z0 = [0.1,0.1,1,0]
 
-prob = BK.BifurcationProblem(Pop, z0, par_pop, (@lens _.b0); recordFromSolution = (x, p) -> (x = x[1], y = x[2], u = x[3]))
+prob = BifurcationProblem(Pop!, z0, par_pop, (@lens _.b0); recordFromSolution = (x, p) -> (x = x[1], y = x[2], u = x[3]))
 
 opts_br = ContinuationPar(pMin = 0., pMax = 20.0, ds = 0.002, dsmax = 0.01, nInversion = 6, detectBifurcation = 3, maxBisectionSteps = 25, nev = 4, maxSteps = 20000)
 @set! opts_br.newtonOptions.verbose = true
@@ -236,7 +236,7 @@ getNormalForm(brpo_pd, 2, prm = true)
 opts_pocoll_pd = ContinuationPar(brpo_pd.contparams, detectBifurcation = 3, maxSteps = 40, pMin = 1.e-2, plotEveryStep = 1, dsmax = 1e-2, ds = 1e-3)
 @set! opts_pocoll_pd.newtonOptions.tol = 1e-10
 pd_po_coll2 = continuation(brpo_pd, 2, (@lens _.b0), opts_pocoll_pd;
-		verbosity = 3, plot = true,
+		verbosity = 3, plot = false,
 		detectCodim2Bifurcation = 2,
 		startWithEigen = false,
 		usehessian = false,
@@ -294,7 +294,7 @@ brpo_pd = continuation(prob2, ci, PALC(), ContinuationPar(opts_po_cont, dsmax =
 opts_pocoll_pd = ContinuationPar(brpo_pd.contparams, detectBifurcation = 3, maxSteps = 40, pMin = 1.e-2, plotEveryStep = 1, dsmax = 1e-2, ds = -1e-3)
 @set! opts_pocoll_pd.newtonOptions.tol = 1e-10
 pd_po_coll2 = continuation(brpo_pd, 2, (@lens _.b0), opts_pocoll_pd;
-		verbosity = 3, plot = true,
+		verbosity = 3, plot = false,
 		detectCodim2Bifurcation = 2,
 		startWithEigen = false,
 		usehessian = false,
@@ -336,7 +336,7 @@ pt = getNormalForm(brpo_pd_sh, 1)
 
 # codim 2 Fold
 opts_posh_fold = ContinuationPar(br_fold_sh.contparams, detectBifurcation = 3, maxSteps = 200, pMin = 0.01, pMax = 1.2)
-@error "ce foire la precision si cette tol est trop petite"
+@error "ça foire la precision si cette tol est trop petite"
 @set! opts_posh_fold.newtonOptions.tol = 1e-12
 fold_po_sh1 = continuation(br_fold_sh, 2, (@lens _.ϵ), opts_posh_fold;
 		verbosity = 2, plot = true,
@@ -361,7 +361,7 @@ fold_po_sh2 = continuation(br_fold_sh, 1, (@lens _.ϵ), opts_posh_fold;
 # codim 2 PD
 opts_posh_pd = ContinuationPar(brpo_pd_sh.contparams, detectBifurcation = 3, maxSteps = 40, pMin = -1.)
 @set! opts_posh_pd.newtonOptions.tol = 1e-12
-@error "ce foire la precision si cette tol est trop petite"
+@error "ça foire la precision si cette tol est trop petite"
 @set! opts_posh_pd.newtonOptions.verbose = true
 pd_po_sh = continuation(brpo_pd_sh, 1, (@lens _.b0), opts_posh_pd;
 		verbosity = 3, plot = true,
@@ -378,7 +378,7 @@ pd_po_sh = continuation(brpo_pd_sh, 1, (@lens _.b0), opts_posh_pd;
 plot(pd_po_sh)
 
 plot(fold_po_sh1, fold_po_sh2, branchlabel = ["FOLD", "FOLD"])
-	plot!(pd_po_sh, vars = (:ϵ, :b0), branchlabel = "PD")
+plot!(pd_po_sh, vars = (:ϵ, :b0), branchlabel = "PD")
 
 #####
 # find the NS case

src/Problems.jl

@@ -3,6 +3,7 @@ abstract type AbstractBifurcationProblem end
 # this type is based on the type BifFunction, see below
 # it provides all derivatives
 abstract type AbstractAllJetBifProblem <: AbstractBifurcationProblem end
+using SciMLBase: numargs
 
 getVectorType(::AbstractBifurcationProblem) = Nothing
 isInplace(::Union{AbstractBifurcationProblem, Nothing}) = false
@@ -12,6 +13,15 @@ isInplace(::Union{AbstractBifurcationProblem, Nothing}) = false
 # mesh for periodic orbit) and this approach seems efficient
 @inline getSolution(::AbstractBifurcationProblem, x) = x
 
+"""
+Determine if the vector field is of the form `f!(out,z,p)`.
+"""
+function isinplace(f)
+	m = minimum(numargs(f))
+	@assert 1<m<4 "You have too many/few arguments in your vector field F. It should be of the form `F(x,p)` or `F!(x,p)`."
+	return m == 3
+end
+
 """
 $(TYPEDEF)
 
@@ -124,6 +134,8 @@ $(TYPEDFIELDS)
 			end
 
 			getVectorType(::$op{Tvf, Tu, Tp, Tl, Tplot, Trec}) where {Tvf, Tu, Tp, Tl, Tplot, Trec} = Tu
+			plotSolution(prob::$op) = prob.plotSolution
+			recordFromSolution(prob::$op) = prob.recordFromSolution
 		end
 	else
         """
@@ -152,7 +164,7 @@ $(TYPEDFIELDS)
 	# forward getters
 	if op in (:BifurcationProblem, :ODEBifProblem, :PDEBifProblem)
 		@eval begin
-			function $op(F, u0, parms, lens = (@lens _);
+			function $op(_F, u0, parms, lens = (@lens _);
 							dF = nothing,
 							dFad = nothing,
 							J = nothing,
@@ -164,6 +176,12 @@ $(TYPEDFIELDS)
 							plotSolution = plotDefault,
 							delta = convert(eltype(u0), 1e-8),
 							inplace = false)
+				if inplace
+					F = _F
+				else
+					iip = isinplace(_F)
+					F = iip ? (x,p) -> _F(similar(x),x,p) : _F	
+				end
 				J = isnothing(J) ? (x,p) -> ForwardDiff.jacobian(z -> F(z, p), x) : J
 				dF = isnothing(dF) ? (x,p,dx) -> ForwardDiff.derivative(t -> F(x .+ t .* dx, p), zero(eltype(dx))) : dF
 				d1Fad(x,p,dx1) = ForwardDiff.derivative(t -> F(x .+ t .* dx1, p), zero(eltype(dx1)))
@@ -210,6 +228,7 @@ d3F(pb::AbstractAllJetBifProblem, x, p, dx1, dx2, dx3) = d3F(pb.VF, x, p, dx1, d
 d3Fc(pb::AbstractAllJetBifProblem, x, p, dx1, dx2, dx3) = d3Fc(pb.VF, x, p, dx1, dx2, dx3)
 hasHessian(pb::AbstractAllJetBifProblem) = hasHessian(pb.VF)
 hasAdjoint(pb::AbstractAllJetBifProblem) = hasAdjoint(pb.VF)
+hasAdjointMF(pb::AbstractAllJetBifProblem) = hasAdjointMF(pb.VF)
 getDelta(pb::AbstractAllJetBifProblem) = getDelta(pb.VF)
 
 for op in (:WrapPOTrap, :WrapPOSh, :WrapPOColl, :WrapTW)

src/codim2/NormalForms.jl

@@ -245,7 +245,6 @@ function bogdanovTakensNormalForm(prob_ma, L,
 	A22 = [[pAq(p1, q0, lens1), pAq(p0, q0, lens1)+pAq(p1, q1, lens1)] [pAq(p1, q0, lens2), pAq(p0, q0, lens2)+pAq(p1, q1, lens2)] ]
 
 	# solving the linear system of size n+2
-	# @infiltrate
 	c = 3dot(p0, H1100) - dot(p0, B(q1, q1))
 	H0010, K10, cv, it = bls(Val(:Block), L, J1s, (A12_1, A12_2), A22, q1, [dot(p1, B(q1, q1))/2, c])
 	@assert size(H0010) == size(x0)

src/periodicorbit/FlowDE.jl

@@ -1,4 +1,4 @@
-using SciMLBase: remake, solve, ODEProblem, EnsembleProblem, EnsembleThreads, DAEProblem, isinplace
+using SciMLBase: remake, solve, ODEProblem, EnsembleProblem, EnsembleThreads, DAEProblem, isinplace as isinplace_sciml
 
 struct FlowDE{Tprob, Talg, Tjac, TprobMono, TalgMono, Tkwde, Tcb, Tvjp, Tδ} <: AbstractFlow
 	"Store the ODEProblem associated to the flow of the Cauchy problem"
@@ -50,8 +50,8 @@ _getVectorField(prob::ODEProblem, x, p) = prob.f(x, p, prob.tspan[1])
 _getVectorField(prob::EnsembleProblem, x, p) = _getVectorField(prob.prob, x, p)
 _getVectorField(prob::EnsembleProblem, o, x, p) = _getVectorField(prob.prob, o, x, p)
 
-@inline _isinplace(pb::ODEProblem) = isinplace(pb)
-@inline _isinplace(pb::EnsembleProblem) = isinplace(pb.prob)
+@inline _isinplace(pb::ODEProblem) = isinplace_sciml(pb)
+@inline _isinplace(pb::EnsembleProblem) = isinplace_sciml(pb.prob)
 
 function vf(fl::FlowDE, x, p)
 	if _isinplace(fl.prob)

src/periodicorbit/PeriodicOrbitUtils.jl

@@ -111,6 +111,7 @@ end
 function modifyPORecord(probPO, kwargs, par, lens)
 	if :recordFromSolution in keys(kwargs)
 		_recordsol0 = get(kwargs, :recordFromSolution, nothing)
+		@assert ~isnothing(_recordsol0) "Please open an issue on the website."
 		return _recordsol = (x, p; k...) -> _recordsol0(x, (prob = probPO, p = p); k...)
 	else
 		if probPO isa AbstractPODiffProblem

src/periodicorbit/PeriodicOrbits.jl

@@ -135,7 +135,6 @@ FloquetWrapperLS(ls::FloquetWrapperLS) = ls
 
 # this is to use of MatrixBLS
 LinearAlgebra.hcat(shjac::FloquetWrapper, dR) = hcat(shjac.jacpb, dR)
-
 ####################################################################################################
 # different jacobian types which parametrize the way jacobians of PO are computed
 struct AutoDiffDense <: AbstractJacobianMatrix end

src/periodicorbit/ShootingDE.jl

@@ -2,7 +2,7 @@ using SciMLBase: ODEProblem, DAEProblem, EnsembleProblem, terminate!, solve, Vec
 const ODEType = Union{ODEProblem, DAEProblem}
 
 function getVectorField(prob::Union{ODEProblem, DAEProblem})
-	if isinplace(prob)
+	if isinplace_sciml(prob)
 		return (x, p) -> (out = similar(x); prob.f(out, x, p, prob.tspan[1]); return out)
 	else
 		return (x, p) -> prob.f(x, p, prob.tspan[1])

src/plotting/RecipesMakie.jl

@@ -14,6 +14,7 @@ function plot!(ax1, contres::AbstractBranchResult;
 				linewidthunstable = 1.0,
 				linewidthstable = 3.0linewidthunstable,
 				plotcirclesbif = true,
+				branchlabel = "",
 				applytoY = identity,
 				applytoX = identity)
 
@@ -26,7 +27,11 @@ function plot!(ax1, contres::AbstractBranchResult;
 	if hasstability(contres) && plotstability
 		linewidth = map(x -> isodd(x) ? linewidthstable : linewidthunstable, contres.stable)
 	end
-	lines!(ax1, contres, linewidth = linewidth)
+	if branchlabel == ""
+		lines!(ax1, contres, linewidth = linewidth)
+	else
+		lines!(ax1, contres, linewidth = linewidth, label = branchlabel)
+	end
 
 	# display bifurcation points
 	bifpt = filter(x -> (x.type != :none) && (x.type != :endpoint) && (plotfold || x.type != :fold) && (x.idx <= length(contres)-1), contres.specialpoint)
@@ -142,17 +147,19 @@ function plot(contres::AbstractBranchResult;
 	fig
 end
 
-function plot(contresV::AbstractBranchResult...; kP...)
-	if length(contresV) == 0; return ;end
+function plot(brs::AbstractBranchResult...; 
+		branchlabel = fill("", length(brs)),
+		kP...)
+	if length(brs) == 0; return ;end
 	fig = Figure(resolution = (1200, 700))
 	ax1 = fig[1, 1] = Axis(fig)
 
-	for (id, contres) in pairs(contresV)
+	for (id, contres) in pairs(brs)
 
 		ind1, ind2 = getPlotVars(contres, nothing)
 		xlab, ylab = getAxisLabels(ind1, ind2, contres)
 
-		plot!(ax1, contres; kP...)
+		plot!(ax1, contres; branchlabel = branchlabel[id], kP...)
 
 	end
 	GLMakie.axislegend(ax1, merge = true, unique = true)
@@ -234,4 +241,5 @@ function _plot_bifdiag_makie!(ax, bd::Vector{BifDiagNode}; code = (), level = (-
 	for b in bd
 		_plot_bifdiag_makie!(ax, b; code, level )
 	end
-end
+end
+####################################################################################################

src/plotting/RecipesPlots.jl

@@ -73,7 +73,7 @@ RecipesBase.@recipe function Plots(brs::AbstractBranchResult...;
 		vars = nothing,
 		plotstability = true,
 		plotspecialpoints = true,
-		branchlabel = fill("",length(brs)),
+		branchlabel = fill("", length(brs)),
 		linewidthunstable = 1.0,
 		linewidthstable = 2linewidthunstable,
 		applytoY = identity,

test/COModel.jl

@@ -3,8 +3,6 @@
 using Test, ForwardDiff, Parameters, LinearAlgebra
 using BifurcationKit, Test
 const BK = BifurcationKit
-
-norminf(x) = norm(x, Inf)
 ####################################################################################################
 function COm(u, p)
 	@unpack q1, q2, q3, q4, q5, q6, k = p
@@ -21,7 +19,7 @@ par_com = (q1 = 2.5, q2 = 1., q3 = 10., q4 = 0.0675, q5 = 1., q6 = 0.1, k = 0.4)
 
 z0 = [0.001137, 0.891483, 0.062345]
 
-prob = BK.BifurcationProblem(COm, z0, par_com, (@lens _.q2);
+prob = BifurcationProblem(COm, z0, par_com, (@lens _.q2);
 		recordFromSolution = (x, p) -> (x = x[1], y = x[2], s = x[3]))
 
 opts_br = ContinuationPar(pMin = 0.5, pMax = 2.3, ds = 0.002, dsmax = 0.01, nInversion = 6, detectBifurcation = 3, maxBisectionSteps = 25, nev = 3, maxSteps = 100)

test/codim2.jl

@@ -3,8 +3,6 @@ using BifurcationKit
 using Test, Parameters, Setfield, LinearAlgebra
 # using Plots
 const BK = BifurcationKit
-
-norminf(x) = norm(x, Inf)
 ####################################################################################################
 function COm(u, p)
 	@unpack q1,q2,q3,q4,q5,q6,k = p
@@ -20,7 +18,7 @@ end
 par_com = (q1 = 2.5, q2 = 2.0, q3 = 10., q4 = 0.0675, q5 = 1., q6 = 0.1, k = 0.4)
 z0 = [0.07,0.2,05]
 
-prob = BK.BifurcationProblem(COm, z0, par_com, (@lens _.q2); recordFromSolution = (x, p) -> (x = x[1], y = x[2], s = x[3]))
+prob = BifurcationProblem(COm, z0, par_com, (@lens _.q2); recordFromSolution = (x, p) -> (x = x[1], y = x[2], s = x[3]))
 
 opts_br = ContinuationPar(pMin = 0.6, pMax = 2.5, ds = 0.002, dsmax = 0.01, nInversion = 4, detectBifurcation = 3, maxBisectionSteps = 25, nev = 2, maxSteps = 20000)
 

test/codim2PO-OColl.jl

@@ -2,8 +2,6 @@
 using Test, ForwardDiff, Parameters, LinearAlgebra
 using BifurcationKit, Test
 const BK = BifurcationKit
-
-norminf(x) = norm(x, Inf)
 ###################################################################################################
 function Pop!(du, X, p, t = 0)
 	@unpack r,K,a,ϵ,b0,e,d = p
@@ -16,13 +14,12 @@ function Pop!(du, X, p, t = 0)
 	du[4] = 2pi * u + v - s * v
 	du
 end
-Pop(u,p) = Pop!(similar(u),u,p,0)
 
 par_pop = ( K = 1., r = 2π, a = 4π, b0 = 0.25, e = 1., d = 2π, ϵ = 0.2, )
 
 z0 = [0.1,0.1,1,0]
 
-prob = BK.BifurcationProblem(Pop, z0, par_pop, (@lens _.b0); recordFromSolution = (x, p) -> (x = x[1], y = x[2], u = x[3]))
+prob = BifurcationProblem(Pop!, z0, par_pop, (@lens _.b0); recordFromSolution = (x, p) -> (x = x[1], y = x[2], u = x[3]))
 
 opts_br = ContinuationPar(pMin = 0., pMax = 20.0, ds = 0.002, dsmax = 0.01, nInversion = 6, detectBifurcation = 3, maxBisectionSteps = 25, nev = 4, maxSteps = 20000)
 @set! opts_br.newtonOptions.verbose = false

test/codim2PO-shooting-mf.jl

@@ -4,8 +4,6 @@ using BifurcationKit, Test
 const BK = BifurcationKit
 
 using ComponentArrays # this is for SciMLSensitivity and adjoint of flow
-
-norminf(x) = norm(x, Inf)
 ###################################################################################################
 function Pop!(du, X, p, t = 0)
 	@unpack r,K,a,ϵ,b0,e,d, = p
@@ -24,7 +22,7 @@ par_pop = ComponentArray( K = 1., r = 2π, a = 4π, b0 = 0.25, e = 1., d = 2π,
 
 z0 = [0.1,0.1,1,0]
 
-prob = BK.BifurcationProblem(Pop, z0, par_pop, (@lens _.b0); recordFromSolution = (x, p) -> (x = x[1], y = x[2], u = x[3]))
+prob = BifurcationProblem(Pop, z0, par_pop, (@lens _.b0); recordFromSolution = (x, p) -> (x = x[1], y = x[2], u = x[3]))
 
 opts_br = ContinuationPar(pMin = 0., pMax = 20.0, ds = 0.002, dsmax = 0.01, nInversion = 6, detectBifurcation = 3, maxBisectionSteps = 25, nev = 4)
 @set! opts_br.newtonOptions.verbose = true

test/codim2PO-shooting.jl

@@ -2,8 +2,6 @@
 using Test, ForwardDiff, Parameters, LinearAlgebra
 using BifurcationKit, Test
 const BK = BifurcationKit
-
-norminf(x) = norm(x, Inf)
 ###################################################################################################
 function Pop!(du, X, p, t = 0)
 	@unpack r,K,a,ϵ,b0,e,d = p
@@ -16,13 +14,12 @@ function Pop!(du, X, p, t = 0)
 	du[4] = 2pi * u + v - s * v
 	du
 end
-Pop(u,p) = Pop!(similar(u),u,p,0)
 
 par_pop = ( K = 1., r = 2π, a = 4π, b0 = 0.25, e = 1., d = 2π, ϵ = 0.2, )
 
 z0 = [0.1,0.1,1,0]
 
-prob = BK.BifurcationProblem(Pop, z0, par_pop, (@lens _.b0); recordFromSolution = (x, p) -> (x = x[1], y = x[2], u = x[3]))
+prob = BifurcationProblem(Pop!, z0, par_pop, (@lens _.b0); recordFromSolution = (x, p) -> (x = x[1], y = x[2], u = x[3]))
 
 opts_br = ContinuationPar(pMin = 0., pMax = 20.0, ds = 0.002, dsmax = 0.01, nInversion = 6, detectBifurcation = 3, maxBisectionSteps = 25, nev = 4, maxSteps = 200)
 ################################################################################

test/event.jl

@@ -4,9 +4,7 @@ using Test
 using BifurcationKit, Setfield, ForwardDiff
 
 const BK = BifurcationKit
-
-
-####################################################################################################
+###################################################################################################
 _eve = BK.ContinuousEvent(2, (iter, state) -> (getp(state)+2, getx(state)[1]-1))
 _eved = BK.DiscreteEvent(1, (iter, state) -> getp(state)>-2)
 @test BK.hasCustomLabels(_eve) == false