From 0116b4fe6605d6b8c84f14f94d341f851e3d215e Mon Sep 17 00:00:00 2001 From: romain veltz Date: Sat, 21 Mar 2026 17:01:26 +0100 Subject: [PATCH 1/2] src: add AutoSparseDI option for jacobian based on DI --- Project.toml | 10 +- ext/SparseDIExt.jl | 62 ++++ src/periodicorbit/PeriodicOrbitCollocation.jl | 291 +++++++++++++++++- test/test_DI.jl | 46 +++ 4 files changed, 404 insertions(+), 5 deletions(-) create mode 100644 ext/SparseDIExt.jl create mode 100644 test/test_DI.jl diff --git a/Project.toml b/Project.toml index 014444b..c581bef 100644 --- a/Project.toml +++ b/Project.toml @@ -14,8 +14,16 @@ RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd" SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" +[weakdeps] +DifferentiationInterface = "a0c0ee7d-e4b9-4e03-894e-1c5f64a51d63" +SparseConnectivityTracer = "9f842d2f-2579-4b1d-911e-f412cf18a3f5" +SparseMatrixColorings = "0a514795-09f3-496d-8182-132a7b665d35" + +[extensions] +SparseDIExt = ["DifferentiationInterface", "SparseConnectivityTracer", "SparseMatrixColorings"] + [compat] -BifurcationKit = "^0.5.5" +BifurcationKit = "^0.5.6" DocStringExtensions = "^0.9" ForwardDiff = "^0.10, ^0.11" NonlinearEigenproblems = "^1.0.1" diff --git a/ext/SparseDIExt.jl b/ext/SparseDIExt.jl new file mode 100644 index 0000000..d6d1a1a --- /dev/null +++ b/ext/SparseDIExt.jl @@ -0,0 +1,62 @@ +module SparseDIExt + +using SparseConnectivityTracer, SparseMatrixColorings +import DifferentiationInterface as DI +import DDEBifurcationKit as DDEBK +import BifurcationKit as BK + +coll_residual_for_di!(result, u, coll, pars) = BK.residual!(coll, result, u, pars) + +function DDEBK.BK._generate_jacobian(coll::BK.PeriodicOrbitOCollProblem{Tprob}, + ::DDEBK.AutoSparseDI, + orbitguess, + pars; + k...) where {Tprob <: DDEBK.AbstractDDEBifurcationProblem} + backend = DI.AutoForwardDiff() + sparse_forward_backend = DI.AutoSparse( + backend; + # we use the following instead of TracerSparsityDetector + # because of searchsortedfirst in POSolution + # sparsity_detector = TracerLocalSparsityDetector(), + sparsity_detector = DI.DenseSparsityDetector(backend, atol = 1e-8), + coloring_algorithm = GreedyColoringAlgorithm(), + ) + out = copy(orbitguess) + + jac_prep_sparse_nonallocating = DI.prepare_jacobian(coll_residual_for_di!, + out, + sparse_forward_backend, + orbitguess, + DI.Constant(coll), + DI.Constant(pars), + # strict = Val(false) + ) + + jac_buffer = similar(sparsity_pattern(jac_prep_sparse_nonallocating), eltype(ones(length(out)))) + L1 = copy(jac_buffer) + return (DDEBK.AutoSparseDI(), + L1, + out, + jac_prep_sparse_nonallocating, + sparse_forward_backend) +end + +function BK.jacobian(coll::BK.PeriodicOrbitOCollProblem, + J::Tuple{DDEBK.AutoSparseDI, T1, T2, T3, T4}, + x, + p) where {T1, T2, T3, T4} + L1 = J[2] + out = J[3] + jac_prep_sparse_nonallocating = J[4] + sparse_forward_backend = J[5] + DI.jacobian!(coll_residual_for_di!, + out, + L1, + jac_prep_sparse_nonallocating, + sparse_forward_backend, + x, + DI.Constant(p)) + return L1 +end + +end # module \ No newline at end of file diff --git a/src/periodicorbit/PeriodicOrbitCollocation.jl b/src/periodicorbit/PeriodicOrbitCollocation.jl index 8cb3f5a..e48bb3c 100644 --- a/src/periodicorbit/PeriodicOrbitCollocation.jl +++ b/src/periodicorbit/PeriodicOrbitCollocation.jl @@ -1,3 +1,7 @@ +# TODO the bottleneck is computhing the jacobian. +# J: 0.011025 seconds (58.15 k allocations: 68.364 MiB) +# L: 0.001172 seconds (16 allocations: 3.752 MiB) + # TODO use getter from BK _get_gauss_nodes(coll) = coll.mesh_cache.gauss_nodes @@ -67,7 +71,8 @@ end end # udj = VectorOfArray([interp(mod(τ * period - d, period)) for d in _delays]) for (ind, d) in enumerate(_delays) - udj.u[ind] .= interp(τ * period - d) + # udj.u[ind] .= interp(τ * period - d) + udj.u[ind] .= BK.__interpolate_posolution(coll, τ - d/period, u, 1) end __po_coll_bc!(coll, outc[:, rg[l]], ∂gj[:, l], gj[:, l], udj, pars, period * dτj, outc[:, end]) end @@ -91,7 +96,16 @@ function BK.jacobian(coll::PeriodicOrbitOCollProblem{Tprob}, return analytical_jacobian_dde_cst(coll, x, p) end -# analytical jacobian for constant DDE +""" +using DifferentiationInterface to automatically derive the sparse jacobian. +""" +struct AutoSparseDI <: BK.AbstractJacobianSparseMatrix end + +function BK._generate_jacobian(coll::PeriodicOrbitOCollProblem{Tprob}, ::AutoSparseDI, orbitguess, pars; k...) where {Tprob <: AbstractDDEBifurcationProblem} + error("You need to import `DifferentiationInterface, SparseConnectivityTracer, SparseMatrixColorings` in order to use this jacobian") +end +######################################################################################## +# analytical jacobians for constant DDE for (fname, floquet) in ((:analytical_jacobian_dde_cst, false), (:analytical_jacobian_dde_cst_floquetgev, true), (:analytical_jacobian_dde_cst_floquetcoll, true), @@ -168,8 +182,8 @@ for (fname, floquet) in ((:analytical_jacobian_dde_cst, false), (ρD * ∂L[l2, l] - α * L[l2, l] * ρI) * In for (idelay, d) in enumerate(delays_v) # find interval where t-τ/period belongs - t0 = τ * period - d - τd = mod(t0, period) / period + t0 = τ - d/period + τd = mod(t0, 1) / 1 index_t = searchsortedfirst(mesh, τd) - 1 @assert 1 <= index_t <= Ntst "We have index_t = $index_t, which is out of bounds for mesh of size $(length(mesh)) and τd = $τd. Please open an issue on the website of BifurcationKit.jl" @@ -181,6 +195,9 @@ for (fname, floquet) in ((:analytical_jacobian_dde_cst, false), Jd[idelay][_rgX, rgNy_delay .+ (l2-1)*n] .+= -α .* JacDDE.Jd[idelay] .* β elseif ($(fname == :analytical_jacobian_dde_cst_floquetcoll) && t0 < 0) rgNy_delay = UnitRange(1, n) .+ ((m * n) * (index_t - 1)) + # fullmesh = coll.mesh_cache.full_mesh + # index_tau = searchsortedlast(fullmesh, τd)- 2 + # rgNy_delay = UnitRange(1, n) .+ index_tau * n Jd[_rgX, rgNy_delay .+ (l2-1)*n] .+= -α .* JacDDE.Jd[idelay] .* β else # case analytical_jacobian_dde_cst J[_rgX, rgNy_delay .+ (l2-1)*n] .+= -α .* JacDDE.Jd[idelay] .* β @@ -211,3 +228,269 @@ for (fname, floquet) in ((:analytical_jacobian_dde_cst, false), end # begin end # for-loop end +######################################################################################## +struct ExtendedSolution{Tpb, Tx, T} + extended_coll::Tpb + xc::Tx # AbstractMatrix + interval::Tuple{T, T} # (-τ/D, 1]) + initial_n::Int # number of unknowns initial +end + +function (sol::ExtendedSolution)(t) + if ~(sol.interval[1] <= t <= sol.interval[2]) + error("You passed t=$t and \n $(sol.interval)") + end + # extended problem + extended_coll = sol.extended_coll + n, m, Ntst = size(extended_coll) + xc = sol.xc + mesh = BK.getmesh(extended_coll) + index_t = searchsortedfirst(mesh, t) - 1 + if index_t <= 0 + return xc[:, 1] + elseif index_t > Ntst + return xc[:, end] + end + @assert mesh[index_t] <= t <= mesh[index_t+1] "Please open an issue on the website of BifurcationKit.jl" + σ = BK.σj(t, mesh, index_t) + # @assert -1 <= σ <= 1 "Strange value of $σ" + σs = BK.get_mesh_coll(extended_coll) + out = zeros(typeof(t), n) + rg = (1:m+1) .+ (index_t - 1) * m + for l in 1:m+1 + out .+= xc[:, rg[l]] .* BK.lagrange(l, σ, σs) + end + out +end + +@views function extended_sol(coll::BK.PeriodicOrbitOCollProblem, + periodic_sol::AbstractVector, + pars) + result = copy(periodic_sol) + periodic_solc = BK.get_time_slices(coll, periodic_sol) + period = BK.getperiod(coll, periodic_sol, nothing) + ratio = maximum(delays(coll.prob_vf, nothing, pars)) / period + + # from the mesh 0 = τ₁ < ... < τₙₜₛₜ₊₁ = 1, we build the new one: + # τ₋ₙ < ratio < ... < τₙₜₛₜ-1 < τ₁ < ... < τₙₜₛₜ₊₁ = 1 + mesh = BK.getmesh(coll) # τᵢ + times = BK.get_times(coll) # tᵢ + + extended_mesh = copy(mesh) + extended_times = copy(times) + extended_solution = copy(periodic_solc) + + for i = reverse(eachindex(mesh)) + new_t = mesh[i] - 1 + if new_t < 0 + pushfirst!(extended_mesh, new_t) + end + if new_t < -ratio + break + end + end + + for i = reverse(eachindex(times)) + new_t = times[i] - 1 + if new_t < 0 + pushfirst!(extended_times, new_t) + extended_solution = hcat(periodic_solc[:, i], extended_solution) + end + if new_t < -ratio + break + end + end + # @error "" ratio + # return extended_mesh + + extended_coll = deepcopy(coll) + @reset extended_coll.mesh_cache.τs = extended_mesh + @reset extended_coll.mesh_cache.Ntst = length(extended_mesh) - 1 + @reset extended_coll.mesh_cache.full_mesh = extended_times + + return ExtendedSolution(extended_coll, extended_solution, (-ratio, 1.0), length(coll)) +end + +# Continuation and Bifurcation Analysis of Delay Differential Equations page 10 + +@views function _residual_for_extended_meshv0!(coll::PeriodicOrbitOCollProblem{Tprob}, + outc::AbstractMatrix{𝒯}, + uc::AbstractMatrix{𝒯}, + period, + (L, ∂L), + pars, + u, # uc is a view of u[1:end-1] + ) where {Tprob <: AbstractDDEBifurcationProblem, 𝒯} + n, m, Ntst = size(coll) + # we want slices at fixed times, hence gj[:, j] is the fastest + # temporaries to reduce allocations + gj = BK.get_tmp(coll.cache.gj, u) # zeros(𝒯, n, m) + ∂gj = BK.get_tmp(coll.cache.∂gj, u) # zeros(𝒯, n, m) + uj = zeros(𝒯, n, m+1) + + # get P.O. interpolation which allows to get result(t) + # interp = BK.POSolution(coll, u, pars) + interp = extended_sol(coll, u, pars) + VF = coll.prob_vf + _delays = delays(VF, gj[:, 1], pars) + + # get the mesh of the collocation problem + mesh = BK.getmesh(coll) + σs = _get_gauss_nodes(coll) + udj = VectorOfArray([copy(uj[:, 1]) for _ in _delays]) + + # range for locating time slices + rg = UnitRange(1, m+1) + eq = 1 + for j in 1:Ntst + uj .= uc[:, rg] + LA.mul!(gj, uj, L) + LA.mul!(∂gj, uj, ∂L) + + # get the delayed states + dτj = (mesh[j+1] - mesh[j]) / 2 + + # compute the collocation residual + if mesh[j]>=0 + for l in 1:m + τ = BK.τj(σs[l], mesh, j) + # udj = VectorOfArray([interp(mod(τ * period - d, period)) for d in _delays]) + for (ind, d) in enumerate(_delays) + udj.u[ind] .= interp(τ - d/period) + # udj.u[ind] .= BK.__interpolate_posolution(coll, τ - d/period, u, 1) + end + __po_coll_bc!(coll, outc[:, rg[l]], ∂gj[:, l], gj[:, l], udj, pars, period * dτj, outc[:, end]) + end + end + rg = rg .+ m + eq += 1 + end + # add the periodicity condition + index = interp.initial_n + # @. outc[:, end] = uc[:, end] - uc[:, end-index+1] + return outc +end + +function _residual_for_extended_mesh(coll_ext::PeriodicOrbitOCollProblem, + interp, + u, + pars, + ) + uc = BK.get_time_slices(coll_ext, u) + period = BK.getperiod(coll_ext, u, nothing) + outc = 0*(uc) .+ 0 + out = vec(outc) + _residual_for_extended_mesh!(coll_ext, interp, outc, uc, period, BK.get_Ls(coll_ext.mesh_cache), pars, u) + return out +end + +@views function _residual_for_extended_mesh!(coll::PeriodicOrbitOCollProblem, + interp, + outc::AbstractMatrix{𝒯}, + uc::AbstractMatrix{𝒯}, + period, + (L, ∂L), + pars, + u, # uc is a view of u[1:end-1] + ) where {𝒯} + n, m, Ntst = size(coll) + # we want slices at fixed times, hence gj[:, j] is the fastest + # temporaries to reduce allocations + gj = BK.get_tmp(coll.cache.gj, u) # zeros(𝒯, n, m) + ∂gj = BK.get_tmp(coll.cache.∂gj, u) # zeros(𝒯, n, m) + uj = zeros(𝒯, n, m+1) + + # get P.O. interpolation which allows to get result(t) + # interp = BK.POSolution(coll, u, pars) + # interp = extended_sol(coll, u, pars) + VF = coll.prob_vf + _delays = delays(VF, gj[:, 1], pars) + + # get the mesh of the collocation problem + mesh = BK.getmesh(coll) + σs = _get_gauss_nodes(coll) + udj = VectorOfArray([copy(uj[:, 1]) for _ in _delays]) + + # range for locating time slices + rg = UnitRange(1, m+1) + eq = 1 + for j in 1:Ntst + uj .= uc[:, rg] + LA.mul!(gj, uj, L) + LA.mul!(∂gj, uj, ∂L) + + # get the delayed states + dτj = (mesh[j+1] - mesh[j]) / 2 + + # compute the collocation residual + if mesh[j]>=0 + for l in 1:m + τ = BK.τj(σs[l], mesh, j) + # udj = VectorOfArray([interp(mod(τ * period - d, period)) for d in _delays]) + for (ind, d) in enumerate(_delays) + udj.u[ind] .= interp(τ - d/period) + # udj.u[ind] .= BK.__interpolate_posolution(coll, τ - d/period, u, 1) + end + __po_coll_bc!(coll, outc[:, rg[l]], ∂gj[:, l], gj[:, l], udj, pars, period * dτj, outc[:, end]) + end + end + rg = rg .+ m + eq += 1 + end + # add the periodicity condition + index = interp.initial_n + # @. outc[:, end] = uc[:, end] - uc[:, end-index+1] +end + +function jacobian_extended_mesh(coll::PeriodicOrbitOCollProblem, + periodic_sol::AbstractVector, + pars) + period = BK.getperiod(coll, periodic_sol, pars) + interp = extended_sol(coll, periodic_sol, pars) + coll_ext = interp.extended_coll + extended_uc = interp.xc + extended_u = vcat(vec(extended_uc), period) + extended_outc = zero(extended_uc) .+ 0 + index = interp.initial_n + @error "" index size(extended_uc) + + # return _residual_for_extended_mesh(coll_ext, interp, extended_u, pars) + + # _residual_for_extended_mesh!(coll_ext, interp, outc, uc, period, BK.get_Ls(coll.mesh_cache), pars, u) + # return outc + + n, m, Ntst = size(coll) + + function residual(extended_u0) + interp0 = ExtendedSolution(coll_ext, extended_u0[1:end-1]', interp.interval, interp.initial_n) + _residual_for_extended_mesh(coll_ext, interp0, extended_u0, pars)[1:end-1-n][end-index+1:end] + end + + J = ForwardDiff.jacobian(residual, extended_u) + ncol = size(J,2)-size(J,1) + B = J[end-index+3:end, end-index:end-1] + A = J[3:index, 1:ncol+1] + + @error "" residual(extended_u) size(A) size(B) + return J, A, B + + + Mₜ = -B\A + Nₜ, N = size(Mₜ) + @error "" size(Mₜ) + if N <= Nₜ + @error "1" N <= Nₜ + It = Itilde(N, Nₜ) + M = It * Mₜ # same as Mₜ[end-N+1:end, :] + else + @error "2" N <= Nₜ + It = Itilde(N - Nₜ, N) + M = vcat(It, Mₜ) + end + + vals = LA.eigvals(M) + logvals = log.(complex.(vals)) + I = sortperm(logvals, by = real, rev = true) + # floquet exponents + σ = logvals[I] #.* 40.10727283620028 +end \ No newline at end of file diff --git a/test/test_DI.jl b/test/test_DI.jl new file mode 100644 index 0000000..5a97491 --- /dev/null +++ b/test/test_DI.jl @@ -0,0 +1,46 @@ +using DDEBifurcationKit, LinearAlgebra, SparseArrays +using BifurcationKit +const BK = BifurcationKit +const DDEBK = DDEBifurcationKit + +function Hutchinson(u, ud, p) + (;a,d,Δ) = p + d .* (Δ*u) .- a .* ud.u[1] .* (1 .+ u) +end + +delaysF(par) = [1.] + +# discretisation +Nx = 100; Lx = pi/2; +X = -Lx .+ 2Lx/Nx*(0:Nx-1) |> collect +h = 2Lx/Nx +Δ = spdiagm(0 => -2ones(Nx), 1 => ones(Nx-1), -1 => ones(Nx-1) ) / h^2; Δ[1,1]=Δ[end,end]=-1/h^2 + +pars = (a = 0.5, d = 1, τ = 1.0, Δ = Δ, N = Nx) +x0 = zeros(Nx) + +################################################################################ +# case where we specify the jacobian +function JacHutchinson(u, p) + (;a,d,Δ) = p + # we compute the jacobian at the steady state + J0 = d * Δ .- a .* Diagonal(u) + J1 = -a .* Diagonal(1 .+ u) + return J0, [J1] +end + +prob = ConstantDDEBifProblem(Hutchinson, delaysF, x0, pars, (@optic _.a); J = JacHutchinson) +optn = NewtonPar(verbose = true, eigsolver = DDE_DefaultEig()) +opts = ContinuationPar(p_max = 10., p_min = 0., newton_options = optn, ds = 0.01, detect_bifurcation = 3, nev = 5, dsmax = 0.2, n_inversion = 4) +br = continuation(prob, PALC(), opts; verbosity = 1, plot = false, normC = norminf) +################################################################################ +# periodic orbits +using DifferentiationInterface, SparseConnectivityTracer, SparseMatrixColorings + +br_pocoll = @time continuation( + br, 3, ContinuationPar(br.contparams; detect_bifurcation = 0, max_steps = 3, newton_options = NewtonPar(eigsolver = DDE_DefaultEig(), verbose = true), plot_every_step = 1), + PeriodicOrbitOCollProblem(20, 4; jacobian = DDEBK.AutoSparseDI()); + verbosity = 2, + # plot = true, + normC = norminf, + ) From 65e26f296890856ea79489cd5edbda9212a1dbb0 Mon Sep 17 00:00:00 2001 From: romain veltz Date: Sat, 21 Mar 2026 17:11:23 +0100 Subject: [PATCH 2/2] src: remove dev code from collocation --- src/periodicorbit/PeriodicOrbitCollocation.jl | 273 +----------------- 1 file changed, 1 insertion(+), 272 deletions(-) diff --git a/src/periodicorbit/PeriodicOrbitCollocation.jl b/src/periodicorbit/PeriodicOrbitCollocation.jl index e48bb3c..83a0bf3 100644 --- a/src/periodicorbit/PeriodicOrbitCollocation.jl +++ b/src/periodicorbit/PeriodicOrbitCollocation.jl @@ -1,8 +1,3 @@ -# TODO the bottleneck is computhing the jacobian. -# J: 0.011025 seconds (58.15 k allocations: 68.364 MiB) -# L: 0.001172 seconds (16 allocations: 3.752 MiB) - -# TODO use getter from BK _get_gauss_nodes(coll) = coll.mesh_cache.gauss_nodes @views function BK.residual!(coll::PeriodicOrbitOCollProblem{Tprob}, @@ -227,270 +222,4 @@ for (fname, floquet) in ((:analytical_jacobian_dde_cst, false), end # function end end # begin -end # for-loop end -######################################################################################## -struct ExtendedSolution{Tpb, Tx, T} - extended_coll::Tpb - xc::Tx # AbstractMatrix - interval::Tuple{T, T} # (-τ/D, 1]) - initial_n::Int # number of unknowns initial -end - -function (sol::ExtendedSolution)(t) - if ~(sol.interval[1] <= t <= sol.interval[2]) - error("You passed t=$t and \n $(sol.interval)") - end - # extended problem - extended_coll = sol.extended_coll - n, m, Ntst = size(extended_coll) - xc = sol.xc - mesh = BK.getmesh(extended_coll) - index_t = searchsortedfirst(mesh, t) - 1 - if index_t <= 0 - return xc[:, 1] - elseif index_t > Ntst - return xc[:, end] - end - @assert mesh[index_t] <= t <= mesh[index_t+1] "Please open an issue on the website of BifurcationKit.jl" - σ = BK.σj(t, mesh, index_t) - # @assert -1 <= σ <= 1 "Strange value of $σ" - σs = BK.get_mesh_coll(extended_coll) - out = zeros(typeof(t), n) - rg = (1:m+1) .+ (index_t - 1) * m - for l in 1:m+1 - out .+= xc[:, rg[l]] .* BK.lagrange(l, σ, σs) - end - out -end - -@views function extended_sol(coll::BK.PeriodicOrbitOCollProblem, - periodic_sol::AbstractVector, - pars) - result = copy(periodic_sol) - periodic_solc = BK.get_time_slices(coll, periodic_sol) - period = BK.getperiod(coll, periodic_sol, nothing) - ratio = maximum(delays(coll.prob_vf, nothing, pars)) / period - - # from the mesh 0 = τ₁ < ... < τₙₜₛₜ₊₁ = 1, we build the new one: - # τ₋ₙ < ratio < ... < τₙₜₛₜ-1 < τ₁ < ... < τₙₜₛₜ₊₁ = 1 - mesh = BK.getmesh(coll) # τᵢ - times = BK.get_times(coll) # tᵢ - - extended_mesh = copy(mesh) - extended_times = copy(times) - extended_solution = copy(periodic_solc) - - for i = reverse(eachindex(mesh)) - new_t = mesh[i] - 1 - if new_t < 0 - pushfirst!(extended_mesh, new_t) - end - if new_t < -ratio - break - end - end - - for i = reverse(eachindex(times)) - new_t = times[i] - 1 - if new_t < 0 - pushfirst!(extended_times, new_t) - extended_solution = hcat(periodic_solc[:, i], extended_solution) - end - if new_t < -ratio - break - end - end - # @error "" ratio - # return extended_mesh - - extended_coll = deepcopy(coll) - @reset extended_coll.mesh_cache.τs = extended_mesh - @reset extended_coll.mesh_cache.Ntst = length(extended_mesh) - 1 - @reset extended_coll.mesh_cache.full_mesh = extended_times - - return ExtendedSolution(extended_coll, extended_solution, (-ratio, 1.0), length(coll)) -end - -# Continuation and Bifurcation Analysis of Delay Differential Equations page 10 - -@views function _residual_for_extended_meshv0!(coll::PeriodicOrbitOCollProblem{Tprob}, - outc::AbstractMatrix{𝒯}, - uc::AbstractMatrix{𝒯}, - period, - (L, ∂L), - pars, - u, # uc is a view of u[1:end-1] - ) where {Tprob <: AbstractDDEBifurcationProblem, 𝒯} - n, m, Ntst = size(coll) - # we want slices at fixed times, hence gj[:, j] is the fastest - # temporaries to reduce allocations - gj = BK.get_tmp(coll.cache.gj, u) # zeros(𝒯, n, m) - ∂gj = BK.get_tmp(coll.cache.∂gj, u) # zeros(𝒯, n, m) - uj = zeros(𝒯, n, m+1) - - # get P.O. interpolation which allows to get result(t) - # interp = BK.POSolution(coll, u, pars) - interp = extended_sol(coll, u, pars) - VF = coll.prob_vf - _delays = delays(VF, gj[:, 1], pars) - - # get the mesh of the collocation problem - mesh = BK.getmesh(coll) - σs = _get_gauss_nodes(coll) - udj = VectorOfArray([copy(uj[:, 1]) for _ in _delays]) - - # range for locating time slices - rg = UnitRange(1, m+1) - eq = 1 - for j in 1:Ntst - uj .= uc[:, rg] - LA.mul!(gj, uj, L) - LA.mul!(∂gj, uj, ∂L) - - # get the delayed states - dτj = (mesh[j+1] - mesh[j]) / 2 - - # compute the collocation residual - if mesh[j]>=0 - for l in 1:m - τ = BK.τj(σs[l], mesh, j) - # udj = VectorOfArray([interp(mod(τ * period - d, period)) for d in _delays]) - for (ind, d) in enumerate(_delays) - udj.u[ind] .= interp(τ - d/period) - # udj.u[ind] .= BK.__interpolate_posolution(coll, τ - d/period, u, 1) - end - __po_coll_bc!(coll, outc[:, rg[l]], ∂gj[:, l], gj[:, l], udj, pars, period * dτj, outc[:, end]) - end - end - rg = rg .+ m - eq += 1 - end - # add the periodicity condition - index = interp.initial_n - # @. outc[:, end] = uc[:, end] - uc[:, end-index+1] - return outc -end - -function _residual_for_extended_mesh(coll_ext::PeriodicOrbitOCollProblem, - interp, - u, - pars, - ) - uc = BK.get_time_slices(coll_ext, u) - period = BK.getperiod(coll_ext, u, nothing) - outc = 0*(uc) .+ 0 - out = vec(outc) - _residual_for_extended_mesh!(coll_ext, interp, outc, uc, period, BK.get_Ls(coll_ext.mesh_cache), pars, u) - return out -end - -@views function _residual_for_extended_mesh!(coll::PeriodicOrbitOCollProblem, - interp, - outc::AbstractMatrix{𝒯}, - uc::AbstractMatrix{𝒯}, - period, - (L, ∂L), - pars, - u, # uc is a view of u[1:end-1] - ) where {𝒯} - n, m, Ntst = size(coll) - # we want slices at fixed times, hence gj[:, j] is the fastest - # temporaries to reduce allocations - gj = BK.get_tmp(coll.cache.gj, u) # zeros(𝒯, n, m) - ∂gj = BK.get_tmp(coll.cache.∂gj, u) # zeros(𝒯, n, m) - uj = zeros(𝒯, n, m+1) - - # get P.O. interpolation which allows to get result(t) - # interp = BK.POSolution(coll, u, pars) - # interp = extended_sol(coll, u, pars) - VF = coll.prob_vf - _delays = delays(VF, gj[:, 1], pars) - - # get the mesh of the collocation problem - mesh = BK.getmesh(coll) - σs = _get_gauss_nodes(coll) - udj = VectorOfArray([copy(uj[:, 1]) for _ in _delays]) - - # range for locating time slices - rg = UnitRange(1, m+1) - eq = 1 - for j in 1:Ntst - uj .= uc[:, rg] - LA.mul!(gj, uj, L) - LA.mul!(∂gj, uj, ∂L) - - # get the delayed states - dτj = (mesh[j+1] - mesh[j]) / 2 - - # compute the collocation residual - if mesh[j]>=0 - for l in 1:m - τ = BK.τj(σs[l], mesh, j) - # udj = VectorOfArray([interp(mod(τ * period - d, period)) for d in _delays]) - for (ind, d) in enumerate(_delays) - udj.u[ind] .= interp(τ - d/period) - # udj.u[ind] .= BK.__interpolate_posolution(coll, τ - d/period, u, 1) - end - __po_coll_bc!(coll, outc[:, rg[l]], ∂gj[:, l], gj[:, l], udj, pars, period * dτj, outc[:, end]) - end - end - rg = rg .+ m - eq += 1 - end - # add the periodicity condition - index = interp.initial_n - # @. outc[:, end] = uc[:, end] - uc[:, end-index+1] -end - -function jacobian_extended_mesh(coll::PeriodicOrbitOCollProblem, - periodic_sol::AbstractVector, - pars) - period = BK.getperiod(coll, periodic_sol, pars) - interp = extended_sol(coll, periodic_sol, pars) - coll_ext = interp.extended_coll - extended_uc = interp.xc - extended_u = vcat(vec(extended_uc), period) - extended_outc = zero(extended_uc) .+ 0 - index = interp.initial_n - @error "" index size(extended_uc) - - # return _residual_for_extended_mesh(coll_ext, interp, extended_u, pars) - - # _residual_for_extended_mesh!(coll_ext, interp, outc, uc, period, BK.get_Ls(coll.mesh_cache), pars, u) - # return outc - - n, m, Ntst = size(coll) - - function residual(extended_u0) - interp0 = ExtendedSolution(coll_ext, extended_u0[1:end-1]', interp.interval, interp.initial_n) - _residual_for_extended_mesh(coll_ext, interp0, extended_u0, pars)[1:end-1-n][end-index+1:end] - end - - J = ForwardDiff.jacobian(residual, extended_u) - ncol = size(J,2)-size(J,1) - B = J[end-index+3:end, end-index:end-1] - A = J[3:index, 1:ncol+1] - - @error "" residual(extended_u) size(A) size(B) - return J, A, B - - - Mₜ = -B\A - Nₜ, N = size(Mₜ) - @error "" size(Mₜ) - if N <= Nₜ - @error "1" N <= Nₜ - It = Itilde(N, Nₜ) - M = It * Mₜ # same as Mₜ[end-N+1:end, :] - else - @error "2" N <= Nₜ - It = Itilde(N - Nₜ, N) - M = vcat(It, Mₜ) - end - - vals = LA.eigvals(M) - logvals = log.(complex.(vals)) - I = sortperm(logvals, by = real, rev = true) - # floquet exponents - σ = logvals[I] #.* 40.10727283620028 -end \ No newline at end of file +end # for-loop end \ No newline at end of file