first implementation of condensation of parameters

simplify TMModel

examples/COModel.jl

@@ -58,7 +58,7 @@ args_po = (    record_from_solution = (x, p) -> begin
     plot_solution = plotSolution,
     normC = norminf)
 
-opts_po_cont = ContinuationPar(opts_br, dsmax = 1.95, ds= 2e-2, dsmin = 1e-6, p_max = 5., p_min=-5.,
+opts_po_cont = ContinuationPar(opts_br, dsmax = 1.9, ds= 2e-2, dsmin = 1e-6, p_max = 5., p_min=-5.,
 max_steps = 300, detect_bifurcation = 0, plot_every_step = 10)
 
 # @set! opts_po_cont.newton_options.verbose = false
@@ -130,9 +130,4 @@ hp_codim2 = continuation((@set br.alg.tangent = Bordered()), 2, (@lens _.k), Con
 
 BK.plot(sn_codim2, vars=(:q2, :x), branchlabel = "Fold", plotcirclesbif = true)
 plot!(hp_codim2, vars=(:q2, :x), branchlabel = "Hopf",plotcirclesbif = true)
-plot!(br,xlims=(0.6,1.5))
-
-BK.plot(sn_codim2, vars=(:k, :q2), branchlabel = "Fold")
-plot!(hp_codim2, vars=(:k, :q2), branchlabel = "Hopf",)
-
-BK.plot(hp_codim2, vars=(:q2, :x), branchlabel = "Hopf")
+plot!(br,xlims=(0.6,1.5))

examples/TMModel.jl

@@ -11,21 +11,21 @@ function TMvf!(dz, z, p, t = 0)
     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[3] = (U0 - u) / τF +  U0 * (1 - u) * E
     dz
 end
 
-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
+par_tm = (α = 1.5, τ = 0.013, J = 3.07, E0 = -2.0, τD = 0.200, U0 = 0.3, τF = 1.5, τS = 0.007)
 z0 = [0.238616, 0.982747, 0.367876 ]
 prob = BifurcationProblem(TMvf!, z0, par_tm, (@lens _.E0); record_from_solution = (x, p) -> (E = x[1], x = x[2], u = x[3]),)
 
-opts_br = ContinuationPar(p_min = -10.0, p_max = -0.9, ds = 0.04, dsmax = 0.1, n_inversion = 8, detect_bifurcation = 3, max_bisection_steps = 25, nev = 3)
+opts_br = ContinuationPar(p_min = -10.0, p_max = -0.9, dsmax = 0.1, n_inversion = 8, nev = 3)
 br = continuation(prob, PALC(tangent = Bordered()), opts_br; plot = true, normC = norminf)
 
 BK.plot(br, plotfold=false)
 ####################################################################################################
 # continuation parameters
-opts_po_cont = ContinuationPar(dsmax = 0.1, ds= 0.0005, dsmin = 1e-4, p_max = 0., p_min=-2., max_steps = 20, newton_options = NewtonPar(tol = 1e-10,  max_iterations = 9), nev = 3, tol_stability = 1e-6, detect_bifurcation = 2, plot_every_step = 20)
+opts_po_cont = ContinuationPar(opts_br, dsmin = 1e-4, ds = 1e-4, max_steps = 80, tol_stability = 1e-6, detect_bifurcation = 2, plot_every_step = 20)
 
 # arguments for periodic orbits
 function plotSolution(x, p; k...)
@@ -48,17 +48,17 @@ args_po = (	record_from_solution = (x, p) -> begin
     )
 
 br_potrap = continuation(br, 4, opts_po_cont,
-    PeriodicOrbitTrapProblem(M = 250, jacobian = :Dense, update_section_every_step = 1);
+    PeriodicOrbitTrapProblem(M = 150);
     verbosity = 2, plot = true,
     args_po...,
-    callback_newton = BK.cbMaxNorm(10.),
+    callback_newton = BK.cbMaxNorm(1.),
     )
 
 plot(br, br_potrap, markersize = 3)
 plot!(br_potrap.param, br_potrap.min, label = "")
 ####################################################################################################
 # continuation parameters
-opts_po_cont = ContinuationPar(dsmax = 0.1, ds= 0.0001, dsmin = 1e-4, p_max = 0., p_min=-5., max_steps = 110, newton_options = NewtonPar(verbose = false, tol = 1e-10,  max_iterations = 10), nev = 3, tol_stability = 1e-5, detect_bifurcation = 2, plot_every_step = 40, save_sol_every_step=1)
+opts_po_cont = ContinuationPar(opts_br, ds= 0.0001, dsmin = 1e-4, max_steps = 10, tol_stability = 1e-5, detect_bifurcation = 2, plot_every_step = 10)
 
 br_pocoll = @time continuation(
     br, 4, opts_po_cont,
@@ -81,14 +81,14 @@ using DifferentialEquations
 # this is the ODEProblem used with `DiffEqBase.solve`
 prob_ode = ODEProblem(TMvf!, copy(z0), (0., 1000.), par_tm; abstol = 1e-11, reltol = 1e-9)
 
-opts_po_cont = ContinuationPar(dsmax = 0.09, ds= -0.0001, dsmin = 1e-4, p_max = 0., p_min=-5., max_steps = 120, newton_options = NewtonPar(tol = 1e-6, max_iterations = 7), nev = 3, tol_stability = 1e-8, detect_bifurcation = 2, plot_every_step = 10, save_sol_every_step=1)
+opts_po_cont = ContinuationPar(opts_br, ds= -0.0001, dsmin = 1e-4, max_steps = 120, newton_options = NewtonPar(tol = 1e-6, max_iterations = 7), tol_stability = 1e-8, detect_bifurcation = 2, plot_every_step = 10, save_sol_every_step=1)
 
 br_posh = @time continuation(
     br, 4,
     # arguments for continuation
     opts_po_cont,
     # this is where we tell that we want Standard Shooting
-    ShootingProblem(15, prob_ode, Rodas4P(), parallel = true, jacobian = BK.AutoDiffDense(),);
+    ShootingProblem(15, prob_ode, Rodas4P(), parallel = true,);
     linear_algo = MatrixBLS(),
     verbosity = 2,
     plot = true,
@@ -99,22 +99,19 @@ plot(br_posh, br, markersize=3)
 ####################################################################################################
 # idem with Poincaré shooting
 
-opts_po_cont = ContinuationPar(dsmax = 0.02, ds= 0.001, p_max = 0., p_min=-5., max_steps = 50, newton_options = NewtonPar(tol = 1e-9, max_iterations=15), nev = 3, tol_stability = 1e-6, detect_bifurcation = 2, plot_every_step = 5)
+opts_po_cont = ContinuationPar(opts_br, dsmax = 0.02, ds= 0.0001, max_steps = 50, newton_options = NewtonPar(tol = 1e-9, max_iterations=15), tol_stability = 1e-6, detect_bifurcation = 2, plot_every_step = 5)
 
 br_popsh = @time continuation(
     br, 4,
     # arguments for continuation
     opts_po_cont,
     # this is where we tell that we want Poincaré Shooting
-    PoincareShootingProblem(5, prob_ode, Rodas4P(); parallel = true, jacobian = BK.AutoDiffDenseAnalytical());
-    alg = PALC(tangent = Bordered()),
-    ampfactor = 1.0, δp = 0.005,
+    PoincareShootingProblem(3, prob_ode, Rodas4P(); parallel = true);
     # usedeflation = true,
     linear_algo = MatrixBLS(),
     verbosity = 2, plot = true,
     args_po...,
-    callback_newton = BK.cbMaxNorm(1e1),
-    record_from_solution = (x, p) -> (return (max = getmaximum(p.prob, x, @set par_tm.E0 = p.p), period = getperiod(p.prob, x, @set par_tm.E0 = p.p))),
+    callback_newton = BK.cbMaxNorm(1e0),
     normC = norminf)
 
 plot(br, br_popsh, markersize=3)

src/periodicorbit/PeriodicOrbitCollocation.jl

@@ -1164,16 +1164,6 @@ $(SIGNATURES)
 This function extracts the indices of the blocks composing the matrix J which is a M x M Block matrix where each block N x N has the same sparsity.
 """
 function get_blocks(coll::PeriodicOrbitOCollProblem, Jac::SparseMatrixCSC)
-    # N, m, Ntst = size(coll)
-    # M = div(size(Jac,1)-1, N)
-    # I, J, K = findnz(Jac)
-    # out = [Vector{Int}() for i in 1:M+1, j in 1:M+1];
-    # for k in eachindex(I)
-    #     m, l = div(I[k]-1, N), div(J[k]-1, N)
-    #     push!(out[1+m, 1+l], k)
-    # end
-    # res = [length(m) for m in out]
-    # out
     N, m, Ntst = size(coll)
     blocks = N * ones(Int64, 1 + m * Ntst + 1); blocks[end] = 1
     n_blocks = length(blocks)
@@ -1185,3 +1175,81 @@ function get_blocks(coll::PeriodicOrbitOCollProblem, Jac::SparseMatrixCSC)
     end
     out
 end
+####################################################################################################
+@views function condensation_of_parameters(coll::PeriodicOrbitOCollProblem, J, rhs0)
+    N, m, Ntst = size(coll)
+    nbcoll = N * m
+    nⱼ = size(J, 1)
+    𝒯 = eltype(coll)
+
+    P = Matrix{𝒯}(LinearAlgebra.I(nⱼ))
+    blockⱼ = zeros(𝒯, nbcoll, nbcoll)
+    rg = 1:nbcoll
+    for _ in 1:Ntst
+        F = lu(J[rg, rg .+ N])
+        P[rg, rg] .= (F.P \ F.L)
+        rg = rg .+ nbcoll
+    end
+
+    Fₚ = lu(P)
+    Jcond = Fₚ \ J
+    rhs = Fₚ \ rhs0
+
+    Aᵢ = Matrix{𝒯}(undef, N, N)
+    Bᵢ = Matrix{𝒯}(undef, N, N)
+
+    r1 = 1:N
+    r2 = N*(m-1)+1:(m*N)
+    rN = 1:N
+
+    # solving for the external variables
+    In = I(N)
+    Jext = zeros(Ntst*N+N+1, Ntst*N+N+1)
+    Jext[end-N:end-1,end-N:end-1] .= In
+    Jext[end-N:end-1,1:N] .= -In
+    Jext[end, end] = J[end,end]
+    rhs_ext = 𝒯[]
+
+    # we solve for the external unknowns
+    for _ in 1:Ntst
+        Aᵢ .= Jcond[r2, r1]
+        Bᵢ .= Jcond[r2, r1 .+ nbcoll]
+
+        Jext[rN, rN] .= Aᵢ
+        Jext[rN, rN .+ N] .= Bᵢ
+        Jext[end, rN] .= Jcond[end, r1]
+        Jext[end, rN .+ N] .= Jcond[end, r1 .+ nbcoll]
+
+        append!(rhs_ext, copy(rhs[r2]))
+        r1 = r1 .+ nbcoll
+        r2 = r2 .+ nbcoll
+        rN = rN .+ N
+    end
+    append!(rhs_ext, rhs[r1])
+    append!(rhs_ext, [rhs[end]])
+
+    sol_ext = Jext \ rhs_ext
+
+    # we solver for the internal unknowns
+    sol_cop = copy(sol_ext[1:N])
+    r2 = N+1:(m)*N
+    r1 = 1:(m-1)*N
+    rsol = 1:(m-1)*N
+    rN_left = 1:N
+    rN = 1:N
+    for iₜ in 1:Ntst
+        Jtemp = UpperTriangular(Jcond[r1, r2])
+        left_part = Jcond[r1, rN_left]
+        right_part = Jcond[r1, r2[end]+1:r2[end]+N]
+        rhs_tmp = rhs[rsol] - left_part * sol_ext[rN] - right_part * sol_ext[rN .+ N]
+        sol_tmp = Jtemp \ rhs_tmp
+        append!(sol_cop, sol_tmp, sol_ext[rN .+ N])
+        r1 = r1 .+ nbcoll
+        r2 = r2 .+ nbcoll
+        rN_left = rN_left .+ nbcoll
+        rsol = rsol .+ nbcoll
+        rN = rN .+ N
+    end
+    push!(sol_cop, sol_ext[end])
+    sol_cop
+end

test/stuartLandauCollocation.jl

@@ -265,6 +265,21 @@ _indx = BifurcationKit.get_blocks(prob_col, Jco2);
 @time BifurcationKit.jacobian_poocoll_sparse_indx!(prob_col, Jco2, _ci, par_sl, _indx);
 @test norminf(Jco - Jco2) < 1e-14
 ####################################################################################################
+# condensation of parameters
+Ntst = 100
+m = 4
+prob_col = PeriodicOrbitOCollProblem(Ntst, m; prob_vf = prob_ana, N = N, ϕ = rand(N*( 1 + m * Ntst)), xπ = rand(N*( 1 + m * Ntst)))
+_ci = generate_solution(prob_col, t->cos(t) .* ones(N), 2pi);
+Jcofd = ForwardDiff.jacobian(z->prob_col(z, par_sl), _ci);
+Jco = @time BK.analytical_jacobian(prob_col, _ci, par_sl)
+
+Jco[1:end-1,end] .= 0
+Jco[end,1:end-1] .= 0; Jco[end, end] = 1
+
+_rhs = rand(size(Jco, 1))
+sol_bs = @time Jco \ _rhs;
+sol_cop = @time BK.condensation_of_parameters(prob_col, Jco, _rhs);
+####################################################################################################
 # test Hopf aBS
 let
     for jacPO in (BK.AutoDiffDense(), BK.DenseAnalytical(), BK.FullSparse())