add jacobian option MinAugMatrixBased to Fold/Hopf continuation

src/codim2/MinAugFold.jl

@@ -48,31 +48,8 @@ end
     res = 𝐅(x[begin:end-1], x[end], params)
     return vcat(res[1], res[2])
 end
-
 ###################################################################################################
-# Struct to invert the jacobian of the fold MA problem.
-struct FoldLinearSolverMinAug <: AbstractLinearSolver; end
-
-function foldMALinearSolver(x, p::𝒯, 𝐅::FoldProblemMinimallyAugmented, par,
-                            rhsu, rhsp;
-                            debugArray = nothing) where 𝒯
-    ################################################################################################
-    # debugArray is used as a temp to be filled with values used for debugging. If debugArray = nothing,
-    # then no debugging mode is entered. If it is AbstractArray, then it is populated
-    ################################################################################################
-    # Recall that the functional we want to solve is [F(x,p), σ(x,p)] where σ(x,p) is computed in the 
-    # function above. The Jacobian Jfold of the vector field is expressed at (x, p).
-    # We solve Jfold⋅res = rhs := [rhsu, rhsp]
-    # The Jacobian expression of the Fold problem is
-    #           ┌         ┐
-    #  Jfold =  │  J  dpF │
-    #           │ σx   σp │
-    #           └         ┘
-    # where σx := ∂_xσ and σp := ∂_pσ
-    # We recall the expression of
-    #  σx = -< w, d2F(x,p)[v, x2]>
-    # where (w, σ2) is solution of J'w + b σ2 = 0 with <a, w> = 1
-    ########################## Extraction of function names ########################################
+function _get_bordered_terms(𝐅::FoldProblemMinimallyAugmented, x, p::𝒯, par) where 𝒯
     a = 𝐅.a
     b = 𝐅.b
 
@@ -113,6 +90,48 @@ function foldMALinearSolver(x, p::𝒯, 𝐅::FoldProblemMinimallyAugmented, par
     rmul!(dJvdp, 𝒯(1/(2ϵ3)))
     σₚ = -dot(w, dJvdp)
 
+    return (;J_at_xp, JAd_at_xp, dₚF, σₚ, δ, ϵ2, v, w, par0, dJvdp, itv, itw)
+end
+###################################################################################################
+function jacobian(pdpb::FoldMAProblem{Tprob, MinAugMatrixBased}, X, par) where {Tprob}
+    𝐅 = pdpb.prob
+    x = @view X[begin:end-1]
+    p = X[end]
+
+    @unpack J_at_xp, JAd_at_xp, dₚF, σₚ, ϵ2, v, w, par0 = _get_bordered_terms(𝐅, x, p, par)
+
+    u1 = apply_jacobian(𝐅.prob_vf, x + ϵ2 * v, par0, w, true)
+    u2 = apply(JAd_at_xp, w) # TODO ON CONNAIT u2!!
+    σₓ = minus(u2, u1); rmul!(σₓ, 1 / ϵ2)
+
+    [_get_matrix(J_at_xp) dₚF ; σₓ' σₚ]
+end
+###################################################################################################
+# Struct to invert the jacobian of the fold MA problem.
+struct FoldLinearSolverMinAug <: AbstractLinearSolver; end
+
+function foldMALinearSolver(x, p::𝒯, 𝐅::FoldProblemMinimallyAugmented, par,
+                            rhsu, rhsp;
+                            debugArray = nothing) where 𝒯
+    ################################################################################################
+    # debugArray is used as a temp to be filled with values used for debugging. If debugArray = nothing,
+    # then no debugging mode is entered. If it is AbstractArray, then it is populated
+    ################################################################################################
+    # Recall that the functional we want to solve is [F(x,p), σ(x,p)] where σ(x,p) is computed in the 
+    # function above. The Jacobian Jfold of the vector field is expressed at (x, p).
+    # We solve Jfold⋅res = rhs := [rhsu, rhsp]
+    # The Jacobian expression of the Fold problem is
+    #           ┌         ┐
+    #  Jfold =  │  J  dpF │
+    #           │ σx   σp │
+    #           └         ┘
+    # where σx := ∂_xσ and σp := ∂_pσ
+    # We recall the expression of
+    #  σx = -< w, d2F(x,p)[v, x2]>
+    # where (w, σ2) is solution of J'w + b σ2 = 0 with <a, w> = 1
+
+    @unpack J_at_xp, JAd_at_xp, dₚF, σₚ, δ, ϵ2, v, w, par0, itv, itw = _get_bordered_terms(𝐅, x, p, par)
+
     if 𝐅.usehessian == false || has_hessian(𝐅) == false
         # We invert the jacobian of the Fold problem when the Hessian of x -> F(x, p) is not known analytically.
         # apply Jacobian adjoint
@@ -362,6 +381,10 @@ function continuation_fold(prob, alg::AbstractContinuationAlgorithm,
         foldpointguess = vcat(foldpointguess.u, foldpointguess.p)
         prob_f = FoldMAProblem(𝐅, FiniteDifferences(), foldpointguess, par, lens2, prob.plotSolution, prob.recordFromSolution)
         opt_fold_cont = @set options_cont.newton_options.linsolver = options_cont.newton_options.linsolver
+    elseif jacobian_ma == :MinAugMatrixBased
+        foldpointguess = vcat(foldpointguess.u, foldpointguess.p)
+        prob_f = FoldMAProblem(𝐅, MinAugMatrixBased(), foldpointguess, par, lens2, prob.plotSolution, prob.recordFromSolution)
+        opt_fold_cont = @set options_cont.newton_options.linsolver = options_cont.newton_options.linsolver
     else
         prob_f = FoldMAProblem(𝐅, nothing, foldpointguess, par, lens2, prob.plotSolution, prob.recordFromSolution)
         opt_fold_cont = @set options_cont.newton_options.linsolver = FoldLinearSolverMinAug()

src/codim2/MinAugHopf.jl

@@ -12,7 +12,7 @@ end
 ####################################################################################################
 @inline getvec(x, ::HopfProblemMinimallyAugmented) = get_vec_bls(x, 2)
 @inline getp(x, ::HopfProblemMinimallyAugmented) = get_par_bls(x, 2)
-
+###################################################################################################
 # this function encodes the functional
 function (𝐇::HopfProblemMinimallyAugmented)(x, p::𝒯, ω::𝒯, params) where 𝒯
     # These are the equations of the minimally augmented (MA) formulation of the 
@@ -49,34 +49,8 @@ end
     res = 𝐇(x[begin:end-2], x[end-1], x[end], params)
     return vcat(res[1], res[2], res[3])
 end
-################################################################################
-# Struct to invert the jacobian of the Hopf MA problem.
-struct HopfLinearSolverMinAug <: AbstractLinearSolver; end
-
-"""
-This function solves the linear problem associated with a linearization of the minimally augmented formulation of the Hopf bifurcation point. The keyword `debugArray` is used to debug the routine by returning several key quantities.
-"""
-function hopfMALinearSolver(x, p::𝒯, ω::𝒯, 𝐇::HopfProblemMinimallyAugmented, par,
-                             duu, dup, duω;
-                            debugArray = nothing) where 𝒯
-    ################################################################################################
-    # debugArray is used as a temp to be filled with values used for debugging. If debugArray = nothing, then no debugging mode is entered. If it is AbstractVector, then it is populated
-    ################################################################################################
-    # N = length(du) - 2
-    # The Jacobian J of the vector field is expressed at (x, p)
-    # the jacobian expression Jhopf of the hopf problem is
-    #           ┌             ┐
-    #  Jhopf =  │  J  dpF   0 │
-    #           │ σx   σp  σω │
-    #           └             ┘
-    ########## Resolution of the bordered linear system ########
-    # J * dX      + dpF * dp           = du => dX = x1 - dp * x2
-    # The second equation
-    #    <σx, dX> +  σp * dp + σω * dω = du[end-1:end]
-    # thus becomes
-    #   (σp - <σx, x2>) * dp + σω * dω = du[end-1:end] - <σx, x1>
-    # This 2 x 2 system is then solved to get (dp, dω)
-    ############### Extraction of function names #################
+###################################################################################################
+function _get_bordered_terms(𝐇::HopfProblemMinimallyAugmented, x, p::𝒯, ω::𝒯, par) where 𝒯
     a = 𝐇.a
     b = 𝐇.b
 
@@ -115,13 +89,67 @@ function hopfMALinearSolver(x, p::𝒯, ω::𝒯, 𝐇::HopfProblemMinimallyAugm
     # σω = dot(w, Complex{T}(0, 1) * v)
     σω = Complex{𝒯}(0, 1) * dot(w, v)
 
+    return (;J_at_xp, JAd_at_xp, dₚF, σₚ, δ, ϵ2, v, w, par0, dJvdp, itv, itw, σω)
+end
+###################################################################################################
+function jacobian(pdpb::HopfMAProblem{Tprob, MinAugMatrixBased}, X, par) where {Tprob}
+    𝐇 = pdpb.prob
+    x = @view X[begin:end-2]
+    p = X[end-1]
+    ω = X[end]
+    𝒯 = eltype(p)
+
+    @unpack J_at_xp, JAd_at_xp, dₚF, σₚ, ϵ2, v, w, par0, σω = _get_bordered_terms(𝐇, x, p, ω, par)
+
+    cw = conj(w)
+    vr = real(v); vi = imag(v)
+    u1r = apply_jacobian(𝐇.prob_vf, x + ϵ2 * vr, par0, cw, true)
+    u1i = apply_jacobian(𝐇.prob_vf, x + ϵ2 * vi, par0, cw, true)
+    u2 = apply(JAd_at_xp,  cw)
+    σxv2r = @. -(u1r - u2) / ϵ2
+    σxv2i = @. -(u1i - u2) / ϵ2
+    σₓ = @. σxv2r + Complex{𝒯}(0, 1) * σxv2i
+
+    Jhopf = hcat(_get_matrix(J_at_xp), dₚF, zero(dₚF))
+    Jhopf = vcat(Jhopf, vcat(real(σₓ), real(σₚ), real(σω))')
+    Jhopf = vcat(Jhopf, vcat(imag(σₓ), imag(σₚ), imag(σω))')
+end
+################################################################################
+# Struct to invert the jacobian of the Hopf MA problem.
+struct HopfLinearSolverMinAug <: AbstractLinearSolver; end
+
+"""
+This function solves the linear problem associated with a linearization of the minimally augmented formulation of the Hopf bifurcation point. The keyword `debugArray` is used to debug the routine by returning several key quantities.
+"""
+function hopfMALinearSolver(x, p::𝒯, ω::𝒯, 𝐇::HopfProblemMinimallyAugmented, par,
+                            duu, dup, duω;
+                            debugArray = nothing) where 𝒯
+    ################################################################################################
+    # debugArray is used as a temp to be filled with values used for debugging. If debugArray = nothing, then no debugging mode is entered. If it is AbstractVector, then it is populated
+    ################################################################################################
+    # N = length(du) - 2
+    # The Jacobian J of the vector field is expressed at (x, p)
+    # the jacobian expression Jhopf of the hopf problem is
+    #           ┌             ┐
+    #  Jhopf =  │  J  dpF   0 │
+    #           │ σx   σp  σω │
+    #           └             ┘
+    ########## Resolution of the bordered linear system ########
+    # J * dX      + dpF * dp           = du => dX = x1 - dp * x2
+    # The second equation
+    #    <σx, dX> +  σp * dp + σω * dω = du[end-1:end]
+    # thus becomes
+    #   (σp - <σx, x2>) * dp + σω * dω = du[end-1:end] - <σx, x1>
+    # This 2 x 2 system is then solved to get (dp, dω)
+
+    @unpack J_at_xp, JAd_at_xp, dₚF, σₚ, δ, ϵ2, v, w, par0, itv, itw, σω = _get_bordered_terms(𝐇, x, p, ω, par)
+
     # we solve J⋅x1 = duu and J⋅x2 = dₚF
     x1, x2, cv, (it1, it2) = 𝐇.linsolver(J_at_xp, duu, dₚF)
     ~cv && @debug "Linear solver for J did not converge"
 
     # the case of ∂_xσ is a bit more involved
     # we first need to compute the value of ∂_xσ written σx
-    # σx = zeros(Complex{T}, length(x))
     σx = similar(x, Complex{𝒯})
 
     if 𝐇.usehessian == false || has_hessian(𝐇) == false
@@ -366,6 +394,10 @@ function continuation_hopf(prob_vf, alg::AbstractContinuationAlgorithm,
         hopfpointguess = vcat(hopfpointguess.u, hopfpointguess.p)
         prob_h = HopfMAProblem(𝐇, FiniteDifferences(), hopfpointguess, par, lens2, prob_vf.plotSolution, prob_vf.recordFromSolution)
         opt_hopf_cont = @set options_cont.newton_options.linsolver = options_cont.newton_options.linsolver
+    elseif jacobian_ma == :MinAugMatrixBased
+        hopfpointguess = vcat(hopfpointguess.u, hopfpointguess.p)
+        prob_h = HopfMAProblem(𝐇, MinAugMatrixBased(), hopfpointguess, par, lens2, prob_vf.plotSolution, prob_vf.recordFromSolution)
+        opt_hopf_cont = @set options_cont.newton_options.linsolver = options_cont.newton_options.linsolver
     else
         prob_h = HopfMAProblem(𝐇, nothing, hopfpointguess, par, lens2, prob_vf.plotSolution, prob_vf.recordFromSolution)
         opt_hopf_cont = @set options_cont.newton_options.linsolver = HopfLinearSolverMinAug()
@@ -628,7 +660,7 @@ function (eig::HopfEig)(Jma, nev; kwargs...)
 end
 
 @views function (eig::HopfEig)(Jma::AbstractMatrix, nev; kwargs...)
-    eigenelts = eig.eigsolver(Jma[1:end-2,1:end-2], nev; kwargs...)
+    eigenelts = eig.eigsolver(Jma[1:end-2, 1:end-2], nev; kwargs...)
     return eigenelts
 end
 

test/testHopfMA.jl

@@ -66,7 +66,7 @@ end
 
 Jbru_ana(x, p) = ForwardDiff.jacobian(z->Fbru(z,p),x)
 
-n = 100
+n = 10
 par_bru = (α = 2., β = 5.45, D1 = 0.008, D2 = 0.004, l = 0.3)
 sol0 = vcat(par_bru.α * ones(n), par_bru.β/par_bru.α * ones(n))
 prob = BifurcationKit.BifurcationProblem(Fbru, sol0, par_bru, (@optic _.l);
@@ -117,6 +117,11 @@ rhs = rand(length(hopfpt))
 jac_hopf_fd = Jac_hopf_fdMA(Bd2Vec(hopfpt), par_bru)
 sol_fd = jac_hopf_fd \ rhs
 
+# test against analytical jacobian
+_hopf_ma_problem = BK.HopfMAProblem(hopfvariable, BK. MinAugMatrixBased(), Bd2Vec(hopfpt), par_bru, (@optic _.β), prob.plotSolution, prob.recordFromSolution)
+J_ana = BK.jacobian(_hopf_ma_problem, Bd2Vec(hopfpt), par_bru)
+@test norminf(J_ana - jac_hopf_fd) < 1e-3
+
 # create a linear solver
 hopfls = BK.HopfLinearSolverMinAug()
 tmpVecforσ = zeros(ComplexF64, 2+2n)
@@ -145,6 +150,7 @@ pbgopfperso = BK.BifurcationProblem((u, p) -> hopfvariable(u, p),
                 J = (x, p) -> Jac_hopf_MA(x, p, hopfvariable),)
 outhopf = BK.solve(pbgopfperso, Newton(), NewtonPar(verbose = false, linsolver = BK.HopfLinearSolverMinAug()))
 @test BK.converged(outhopf)
+
 # version with analytical Hessian = 2 P(du2) P(du1) QU + 2 PU P(du1) Q(du2) + 2PU P(du2) Q(du1)
 function d2F(x, p1, du1, du2)
     n = div(length(x),2)
@@ -165,18 +171,20 @@ outhopf = newton(br_d2f, 1)
 @test BK.converged(outhopf)
 
 br_hopf = continuation(br, ind_hopf, (@optic _.β),
-            ContinuationPar(dsmin = 0.001, dsmax = 0.05, ds= 0.01, p_max = 6.5, p_min = 0.0, a = 2., max_steps = 3, newton_options = NewtonPar(verbose = false)), jacobian_ma = :minaug)
+            ContinuationPar(dsmin = 0.001, dsmax = 0.05, ds= 0.01, p_max = 6.5, p_min = 0.0, a = 2., max_steps = 3, newton_options = NewtonPar(verbose = false)), 
+            jacobian_ma = :minaug)
 
-br_hopf = continuation(br_d2f, ind_hopf, (@optic _.β), ContinuationPar(dsmin = 0.001, dsmax = 0.05, ds= 0.01, p_max = 6.5, p_min = 0.0, a = 2., max_steps = 3, newton_options = NewtonPar(verbose = false)), jacobian_ma = :minaug)
+br_hopf = continuation(br_d2f, ind_hopf, (@optic _.β), 
+            ContinuationPar(dsmin = 0.001, dsmax = 0.05, ds= 0.01, p_max = 6.5, p_min = 0.0, a = 2., max_steps = 3, newton_options = NewtonPar(verbose = false)), 
+            jacobian_ma = :minaug)
 ####################################################################################################
 ind_hopf = 1
 hopfpt = BK.HopfPoint(br, ind_hopf)
 
 l_hopf = hopfpt.p[1]
-ωH       = hopfpt.p[2] |> abs
+ωH     = hopfpt.p[2] |> abs
 M = 20
 
-
 orbitguess = zeros(2n, M)
 phase = []; scalphase = []
 vec_hopf = geteigenvector(opt_newton.eigsolver, br.eig[br.specialpoint[ind_hopf].idx][2], br.specialpoint[ind_hopf].ind_ev-1)
@@ -192,7 +200,7 @@ l_hopf, Th, orbitguess2, hopfpt, vec_hopf = BK.guess_from_hopf(br, ind_hopf, opt
 
 prob = BifurcationKit.BifurcationProblem(Fbru, sol0, par_bru, (@optic _.l);
         J = Jbru_sp,
-        record_from_solution = (x, p) -> norminf(x))
+        record_from_solution = (x, p; k...) -> norminf(x))
 
 poTrap = PeriodicOrbitTrapProblem(prob, real.(vec_hopf), hopfpt.u, M, 2n    )
 

test/testJacobianFoldDeflation.jl

@@ -86,6 +86,12 @@ Jac_fold_fdMA(u0) = ForwardDiff.jacobian( u -> foldpbVec(u, par_chan), u0)
 J_fold_fd = Jac_fold_fdMA(Bd2Vec(foldpt))
 res_fd =  J_fold_fd \ rhs
 
+# test against analytical jacobian
+_fold_ma_problem = BK.FoldMAProblem(foldpb, BK. MinAugMatrixBased(), Bd2Vec(foldpt), par_chan, (@optic _.β), prob.plotSolution, prob.recordFromSolution)
+J_ana = BK.jacobian(_fold_ma_problem, Bd2Vec(foldpt), par_chan)
+@test norminf(J_ana - J_fold_fd) < 1e-5
+
+###
 Jac_fold_MA(u0, p, pb::FoldProblemMinimallyAugmented) = (return (x=u0, params=p, prob = pb))
 jacFoldSolver = BK.FoldLinearSolverMinAug()
 debugTmpForσ = zeros(n+1,n+1) # temporary array for debugging σ