VectorOfArray is not an AbstractArray anymore

add Utils.jl for dispatch on BilinearMap
correct tests

Project.toml

@@ -22,6 +22,7 @@ ForwardDiff = "^0.10"
 NonlinearEigenproblems = "^1.0.1"
 Parameters = "0.12.3"
 RecursiveArrayTools = "^2.3, ^2.4, ^2.8, ^2.9, ^3"
+Setfield = "0.5.0, 0.7.0, 0.8.0, ^1.1"
 SciMLBase = "^1, ^2"
 julia = "1.9"
 

examples/neuron.jl

@@ -111,12 +111,12 @@ plot(br2, br_pocoll, br_pitchfork);plot!(br_pocoll, vars = (:param,:min))
 
 # plot the periodic orbit
 plot(layout = 2)
-    for ii = 1:10:110
-        solpo = BK.get_periodic_orbit(br_pocoll.γ.prob.prob, br_pocoll.sol[ii].x, 1)
-        plot!(solpo.t ./ solpo.t[end], solpo.u[1,:], label = "", subplot = 1)
-    end
-    xlabel!("t / period", subplot = 1)
-    plot!(br_pocoll, vars = (:param, :period), subplot = 2, xlims=(2.2,2.4))
+for ii = 1:10:110
+    solpo = BK.get_periodic_orbit(br_pocoll.γ.prob.prob, br_pocoll.sol[ii].x, 1)
+    plot!(solpo.t ./ solpo.t[end], solpo.u[1,:], label = "", subplot = 1)
+end
+xlabel!("t / period", subplot = 1)
+plot!(br_pocoll, vars = (:param, :period), subplot = 2, xlims=(2.2,2.4))
 
 plot(br2, br_pocoll, br_pitchfork);plot!(br_pocoll, vars = (:param,:min))
 ############

src/DDEBifurcationKit.jl

@@ -7,6 +7,7 @@ module DDEBifurcationKit
 
 
 	include("Problems.jl")
+	include("Utils.jl")
 	include("NormalForms.jl")
 	include("EigSolver.jl")
 	include("codim2/codim2.jl")
@@ -18,7 +19,4 @@ module DDEBifurcationKit
 
 	export ConstantDDEBifProblem, SDDDEBifProblem
 	export DDE_DefaultEig
-
-
-
 end

src/NormalForms.jl

@@ -45,7 +45,8 @@ function hopf_normal_form(prob::ConstantDDEBifProblem,
     R3 = BK.TrilinearMap((dx1, dx2, dx3) -> BK.d3F(prob, x0c, parbif, dx1, dx2, dx3) ./6 )
 
     # −LΨ001 = R01
-    R01 = (BK.residual(prob, x0, set(parbif, lens, p + δ)) .- BK.residual(prob, x0, set(parbif, lens, p - δ))) ./ (2δ)
+    R01 = (BK.residual(prob, x0, set(parbif, lens, p + δ)) .- 
+           BK.residual(prob, x0, set(parbif, lens, p - δ))) ./ (2δ)
     Ψ001 = Complex.(expθ(L, ls(Δ0, -R01)[1], 0))
 
     # (2iω−L)Ψ200 = R20(ζ,ζ)
@@ -117,14 +118,14 @@ function hopf_normal_form(prob::SDDDEBifProblem,
 end
 
 function BK.hopf_normal_form(prob::AbstractDDEBifurcationProblem,
-                    br::BK.AbstractBranchResult, 
-                    ind_hopf::Int;
-                    nev = length(BK.eigenvalsfrombif(br, id_bif)),
-                    verbose::Bool = false,
-                    lens = BK.getlens(br),
-                    Teigvec = BK._getvectortype(br),
-                    scaleζ = norm,
-                    detailed = false)
+                             br::BK.AbstractBranchResult, 
+                             ind_hopf::Int;
+                             nev = length(BK.eigenvalsfrombif(br, id_bif)),
+                             verbose::Bool = false,
+                             lens = BK.getlens(br),
+                             Teigvec = BK._getvectortype(br),
+                             scaleζ = norm,
+                             detailed = false)
     # the kwargs detailed is only here to allow to extend BK.hopf_normal_form
     @assert br.specialpoint[ind_hopf].type == :hopf "The provided index does not refer to a Hopf Point"
     verbose && println("#"^53*"\n--> Hopf Normal form computation")

src/Problems.jl

@@ -26,7 +26,7 @@ $(TYPEDFIELDS)
 
 """
 struct ConstantDDEBifProblem{Tvf, Tdf, Tu, Td, Tp, Tl <: Lens, Tplot, Trec, Tgets, Tδ} <: AbstractDDEBifurcationProblem
-    "Vector field, typically a [`BifFunction`](@ref)"
+    "Vector field, typically a [`BifFunction`](@ref). For more information, please look at the website https://bifurcationkit.github.io/DDEBifurcationKit.jl/dev/BifProblem"
     VF::Tvf
     "function delays. It takes the parameters and return the non-zero delays in an `AsbtractVector` form. Example: `delays = par -> [1.]`"
     delays::Tdf
@@ -85,7 +85,7 @@ function ConstantDDEBifProblem(F, delayF, u0, parms, lens = (@lens _);
                 δ = convert(eltype(u0), 1e-8)
                 )
 
-    Fc = (xd, p) -> F(xd[1], xd[2:end], p)
+    Fc = (xd, p) -> F(xd.u[1], xd.u[2:end], p)
     # J = isnothing(J) ? (x,p) -> ForwardDiff.jacobian(z -> F(z, p), x) : J
     dF = isnothing(dF) ? (x,p,dx) -> ForwardDiff.derivative(t -> Fc(x .+ t .* dx, p), 0.) : dF
     d1Fad(x, p, dx1) = ForwardDiff.derivative(t -> Fc(x .+ t .* dx1, p), 0.)
@@ -134,7 +134,7 @@ function jacobian_forwarddiff(prob::ConstantDDEBifProblem, x, p)
     xd = VectorOfArray([copy(x) for _ in eachindex(prob.delays0)])
     J0 = ForwardDiff.jacobian(z -> prob.VF.F(z, xd, p), x)
 
-    Jd = [ ForwardDiff.jacobian(z -> prob.VF.F(x, (@set xd[ii] = z), p), x) for ii in eachindex(prob.delays0)]
+    Jd = [ ForwardDiff.jacobian(z -> prob.VF.F(x, (@set xd.u[ii] = z), p), x) for ii in eachindex(prob.delays0)]
     return J0, Jd
 end
 
@@ -249,7 +249,7 @@ $(TYPEDFIELDS)
 
 """
 struct SDDDEBifProblem{Tvf, Tdf, Tu, Td, Tp, Tl <: Lens, Tplot, Trec, Tgets, Tδ} <: AbstractDDEBifurcationProblem
-    "Vector field, typically a [`BifFunction`](@ref)"
+    "Vector field, typically a [`BifFunction`](@ref). For more information, please look at the website https://bifurcationkit.github.io/DDEBifurcationKit.jl/dev/BifProblem"
     VF::Tvf
     "function delays. It takes the parameters and the state and return the non-zero delays in an `AsbtractVector` form. Example: `delays = (par, u) -> [1. + u[1]^2]`"
     delays::Tdf
@@ -346,7 +346,7 @@ function BK.jacobian(prob::SDDDEBifProblem, x, p)
     xd = VectorOfArray([x for _ in eachindex(prob.delays0)])
     J0 = ForwardDiff.jacobian(z -> prob.VF.F(z, xd, p), x)
 
-    Jd = [ ForwardDiff.jacobian(z -> prob.VF.F(x, (@set xd[ii] = z), p), x) for ii in eachindex(prob.delays0)]
+    Jd = [ ForwardDiff.jacobian(z -> prob.VF.F(x, (@set xd.u[ii] = z), p), x) for ii in eachindex(prob.delays0)]
     return JacobianDDE(prob, J0 + sum(Jd), J0, Jd, prob.delays(x, p))
 end
 

src/Utils.jl

@@ -0,0 +1,4 @@
+# we need to add this dispatch as AbstractVectorOfArray is not a AbstractArray
+(b::BK.TrilinearMap)(dx1::T, dx2::T, dx3::T) where {T <: RecursiveArrayTools.AbstractVectorOfArray{<: Real}} = b.tl(dx1, dx2, dx3)
+
+(b::BK.BilinearMap)(dx1::T, dx2::T) where {T <: RecursiveArrayTools.AbstractVectorOfArray{<: Real}} = b.bl(dx1, dx2)

test/codim2.jl

@@ -1,12 +1,12 @@
 # using Revise
 using Test, DDEBifurcationKit
-using Parameters, Setfield
+using Setfield
 using BifurcationKit
 const BK = BifurcationKit
 
 g(z) = (tanh(z − 1) + tanh(1))*cosh(1)^2
 function neuron2VF(x, xd, p)
-    @unpack a,b,c,d = p
+    (;a,b,c,d) = p
     [
         -x[1] - a * g(b*xd[1][1]) + c * g(d*xd[2][2]),
         -x[2] - a * g(b*xd[1][2]) + c * g(d*xd[2][1])
@@ -83,7 +83,7 @@ brfold = continuation(br2, 3, (@lens _.a),
 @test brfold.specialpoint[3].type == :zh
 ################################################################################
 function humpriesVF(x, xd, p)
-    @unpack κ1,κ2,γ,a1,a2,c = p
+    (;κ1,κ2,γ,a1,a2,c) = p
     [
         -γ * x[1] - κ1 * xd[1][1] - κ2 * xd[2][1]
     ]

test/normalform.jl

@@ -1,11 +1,11 @@
 # using Revise
 
-using DDEBifurcationKit, Parameters, Setfield, BifurcationKit
+using DDEBifurcationKit, Setfield, BifurcationKit
 using Test
 const BK = BifurcationKit
 
 function wrightVF(x, xd, p)
-    @unpack a = p
+    (;a) = p
     y = xd[1][1]
     [
         -a * y * (1 + x[1])
@@ -14,7 +14,7 @@ end
 
 delaysF(par) = [1.0]
 
-pars = (a=0.1,b=0.)
+pars = (a=0.1, b=0.)
 x0 = [0.]
 
 prob = ConstantDDEBifProblem(wrightVF, delaysF, x0, pars, (@lens _.a))
@@ -58,7 +58,7 @@ DDEBifurcationKit.newton_hopf(br, 2)
 ##########################################################################################
 # SDDE, test dummy Hopf normal form
 function humpriesVF(x, xd, p)
-    @unpack κ1,κ2,γ,a1,a2,c = p
+    (;κ1,κ2,γ,a1,a2,c) = p
     [
        -γ * x[1] - κ1 * xd[1][1] - κ2 * xd[2][1]
     ]

test/pocoll-sdde.jl

@@ -1,11 +1,11 @@
 # using Revise, Plots
-using DDEBifurcationKit, Parameters, Setfield, LinearAlgebra
+using DDEBifurcationKit, Setfield, LinearAlgebra
 using BifurcationKit
 const BK = BifurcationKit
 const DDEBK = DDEBifurcationKit
 
 function humpriesVF(x, xd, p)
-   @unpack κ1,κ2,γ,a1,a2,c = p
+   (;κ1,κ2,γ,a1,a2,c) = p
    [
       -γ * x[1] - κ1 * xd[1][1] - κ2 * xd[2][1]
    ]

test/pocoll.jl

@@ -1,11 +1,11 @@
 using Test, DDEBifurcationKit
-using Parameters, Setfield, LinearAlgebra
+using Setfield, LinearAlgebra
 using BifurcationKit
 const BK = BifurcationKit
 
 g(z) = (tanh(z − 1) + tanh(1))*cosh(1)^2
 function neuron2VF(x, xd, p)
-    @unpack a,b,c,d = p
+    (;a,b,c,d) = p
     [
         -x[1] - a * g(b*xd[1][1]) + c * g(d*xd[2][2]),
         -x[2] - a * g(b*xd[1][2]) + c * g(d*xd[2][1])