add AutoSwitch continuation method

correct COP linear bordered solver

src/BifurcationKit.jl

@@ -49,6 +49,7 @@ module BifurcationKit
     include("continuation/Palc.jl")
     include("continuation/Multiple.jl")
     include("continuation/MoorePenrose.jl")
+    include("continuation/AutoSwitch.jl")
     include("DeflatedContinuation.jl")
 
     # wip
@@ -132,7 +133,7 @@ module BifurcationKit
     export DeflationOperator, DeflatedProblem
 
     # predictors for continuation
-    export Natural, PALC, Multiple, Secant, Bordered, DefCont, Polynomial, MoorePenrose, MoorePenroseLS
+    export Natural, PALC, Multiple, Secant, Bordered, DefCont, Polynomial, MoorePenrose, MoorePenroseLS, AutoSwitch
 
     # newton methods
     export NewtonPar, newton, newton_palc, newton_hopf, NonLinearSolution

src/continuation/AutoSwitch.jl

@@ -0,0 +1,61 @@
+struct AutoSwitch{Talg, T} <: AbstractContinuationAlgorithm
+    alg::Talg
+    tol_param::T
+end
+
+function AutoSwitch(;alg = PALC(), tol_param = 0.5)
+    return AutoSwitch(alg, tol_param)
+end
+
+Base.empty!(alg::AutoSwitch) = empty!(alg.alg)
+getθ(alg::AutoSwitch) = getθ(alg.alg)
+getdot(alg::AutoSwitch) = getdot(alg.alg)
+getlinsolver(alg::AutoSwitch) = getlinsolver(alg.alg)
+internal_adaptation!(alg::AutoSwitch, onoroff::Bool) = internal_adaptation!(alg.alg, onoroff)
+
+function update(alg::AutoSwitch, contParams::ContinuationPar, linear_algo) 
+    alg2 = update(alg.alg, contParams, linear_algo)
+    @set alg.alg = alg2
+end
+
+function getpredictor!(state::AbstractContinuationState,
+                        iter::AbstractContinuationIterable,
+                        alg_switch::AutoSwitch,
+                        nrm = false)
+    alg = alg_switch.alg
+    # we compute the tangent
+    # if the state has not converged, we dot not update the tangent
+    # state.z has been updated only if converged(state) == true
+    if converged(state)
+        @debug "Update tangent AutoSwitch"
+        gettangent!(state, iter, alg.tangent, getdot(alg))
+    end
+    # then update the predictor state.z_pred
+    addtangent!(state, nrm)
+end
+
+function initialize!(state::AbstractContinuationState,
+                    iter::AbstractContinuationIterable,
+                    alg::AutoSwitch,
+                    nrm = false)
+    initialize!(state, iter, alg.alg, nrm)
+end
+
+function corrector!(state::AbstractContinuationState,
+                    it::AbstractContinuationIterable,
+                    alg::AutoSwitch; 
+                    kwargs...)
+    τ = state.τ
+    λ = τ.p
+    θ = getθ(it)
+    dotθ = getdot(alg.alg)
+    @debug "" (1-θ)*abs(λ) dotθ(τ, θ)
+    if (1-θ)*abs(λ) > alg.tol_param
+        @debug "NATU" λ
+        corrector!(state, it, Natural(); kwargs...)
+    else
+        @debug "PALC" λ
+        corrector!(state, it, alg.alg; kwargs...)
+    end
+    return true
+end

src/periodicorbit/cop.jl

@@ -601,7 +601,8 @@ function (ls::COPBLS)(_Jc, dR,
                     R::AbstractVecOrMat, n::T,
                     ξu::T = T(1), ξp::T = T(1);
                     shift::Ts = nothing, 
-                    Mass::Tm = LinearAlgebra.I, 
+                    Mass::Tm = LinearAlgebra.I,
+                    dotp = nothing,
                     applyξu! = nothing)  where {T <: Number, Ts, Tm}
     Jc = _get_matrix(_Jc) # to handle FloquetWrapper
     if isnothing(shift)

test/simple_continuation.jl

@@ -1,6 +1,7 @@
 #using Revise
 # using Plots
 using Test
+using ForwardDiff
 using BifurcationKit, LinearAlgebra, SparseArrays
 const BK = BifurcationKit
 
@@ -100,6 +101,18 @@ BK.isinplace(prob)
 BK._getvectortype(prob)
 show(prob)
 
+br0 = @time continuation(prob, PALC(), opts)
+plot(br0)
+br0 = @time continuation(prob, 
+                BK.AutoSwitch(
+                    alg = PALC(tangent = Bordered()), 
+                    tol_param = 0.45,
+                    ),
+                ContinuationPar(opts0, max_steps = 47, detect_fold = false, newton_options = NewtonPar(max_iterations = 5), dsmax = 0.051); 
+                verbosity = 1)
+plot(br0, marker = :d)
+# plot!(br0.param, )
+
 br0 = @time continuation(prob, PALC(), opts; callback_newton = BK.cbMaxNormAndΔp(10,10)) #(6.20 k allocations: 409.469 KiB)
 BK._getfirstusertype(br0)
 BK.propertynames(br0)