Skip to content

Commit

Permalink
improve doc strings
Browse files Browse the repository at this point in the history
improve readme
  • Loading branch information
rveltz committed Jun 25, 2023
1 parent 5665d43 commit 365c0b5
Show file tree
Hide file tree
Showing 13 changed files with 79 additions and 51 deletions.
13 changes: 7 additions & 6 deletions README.md
Original file line number Diff line number Diff line change
Expand Up @@ -62,11 +62,16 @@ Most plugins are located in the organization [bifurcationkit](https://github.com
- [GridapBifurcationKit.jl](https://github.com/bifurcationkit/GridapBifurcationKit) bifurcation analysis of PDEs solved with the Finite Elements Method (FEM) using the package [Gridap.jl](https://github.com/gridap/Gridap.jl).
- [PeriodicSchurBifurcationKit.jl](https://github.com/bifurcationkit/PeriodicSchurBifurcationKit.jl) state of the art computation of Floquet coefficients, useful for computing the stability of periodic orbits.

## Examples of bifurcation diagrams
## Overview of capabilities

The list of capabilities is available [here](https://bifurcationkit.github.io/BifurcationKitDocs.jl/dev/capabilities/).

| ![](https://bifurcationkit.github.io/BifurcationKitDocs.jl/dev/tutorials/BDSH1d.png) | ![](https://bifurcationkit.github.io/BifurcationKitDocs.jl/dev/tutorials/mittlemannBD-1.png) |
## Examples of bifurcation diagrams (ODEs and PDEs)

| ![](https://github.com/bifurcationkit/BifurcationKitDocs.jl/blob/main/docs/src/tutorials/ode/nm-per.png?raw=true) | ![](https://bifurcationkit.github.io/BifurcationKitDocs.jl/dev/tutorials/ode/com-fig3.png?raw=true) |
|:-------------:|:-------------:|
| [simple ODE example](https://bifurcationkit.github.io/BifurcationKitDocs.jl/dev/tutorials/ode/tutorialsODE/#Neural-mass-equation-(Hopf-aBS)) | [Codimension 2 (ODE)](https://bifurcationkit.github.io/BifurcationKitDocs.jl/dev/tutorials/ode/tutorialCO/#CO-oxydation-(codim-2)) |
| ![](https://bifurcationkit.github.io/BifurcationKitDocs.jl/dev/tutorials/BDSH1d.png) | ![](https://bifurcationkit.github.io/BifurcationKitDocs.jl/dev/tutorials/mittlemannBD-1.png) |
| [Automatic Bif. Diagram in 1D Swift Hohenberg](https://bifurcationkit.github.io/BifurcationKitDocs.jl/dev/tutorials/Swift-Hohenberg1d/#d-Swift-Hohenberg-equation-(Automatic)) | [Automatic Bif. Diagram in 2D Bratu](https://bifurcationkit.github.io/BifurcationKitDocs.jl/dev/tutorials/mittelmannAuto/#Automatic-diagram-of-2d-Bratu–Gelfand-problem-(Intermediate)) |
| ![](https://bifurcationkit.github.io/BifurcationKitDocs.jl/dev/tutorials/sh2dbranches.png) | ![](https://bifurcationkit.github.io/BifurcationKitDocs.jl/dev/tutorials/bru-po-cont-3br.png) |
| [Snaking in 2D Swift Hohenberg](https://bifurcationkit.github.io/BifurcationKitDocs.jl/dev/tutorials/tutorials2/#d-Swift-Hohenberg-equation:-snaking,-Finite-Differences) | [Periodic orbits in 1D Brusselator](https://bifurcationkit.github.io/BifurcationKitDocs.jl/dev/tutorials/tutorials3/#d-Brusselator-(automatic))
Expand All @@ -75,7 +80,3 @@ Most plugins are located in the organization [bifurcationkit](https://github.com
| ![](https://bifurcationkit.github.io/BifurcationKitDocs.jl/dev/carrier.png) | ![](https://bifurcationkit.github.io/BifurcationKitDocs.jl/dev/tutorials/GPU-branch.png) |
| [Deflated Continuation in Carrier problem](https://bifurcationkit.github.io/BifurcationKitDocs.jl/dev/tutorials/tutorialCarrier/#Deflated-Continuation-in-the-Carrier-Problem) | [2D Swift Hohenberg on GPU](https://bifurcationkit.github.io/BifurcationKitDocs.jl/dev/tutorials/tutorials2b/#d-Swift-Hohenberg-equation-(non-local)-on-the-GPU,-periodic-BC-(Advanced)) |


## Main features

The list of capabilities is available [here](https://bifurcationkit.github.io/BifurcationKitDocs.jl/dev/capabilities/).
8 changes: 4 additions & 4 deletions examples/SHpde_snaking.jl
Original file line number Diff line number Diff line change
Expand Up @@ -37,9 +37,9 @@ prob = BifurcationProblem(R_SH, sol0, parSH, (@lens _.λ); J = Jac_sp,
plotSolution = (x, p;kwargs...)->(plot!(X, x; ylabel="solution", label="", kwargs...)))
####################################################################################################
optnew = NewtonPar(verbose = false, tol = 1e-12)
# allocations 357, 0.8ms
sol1 = @time newton(prob, optnew)
Plots.plot(X, sol1.u)
# allocations 357, 0.8ms
sol1 = @time newton(prob, optnew)
Plots.plot(X, sol1.u)

opts = ContinuationPar(dsmin = 0.0001, dsmax = 0.01, ds = 0.01,
newtonOptions = setproperties(optnew; maxIter = 30, tol = 1e-8), pMax = 1.,
Expand All @@ -63,7 +63,7 @@ kwargsC = (verbosity = 3,
callbackN = cb
)

brflat = @time continuation(prob, PALC(), opts; kwargsC...)
brflat = @time continuation(prob, PALC(), opts; kwargsC..., verbosity = 0)

plot(brflat, putspecialptlegend = false)
####################################################################################################
Expand Down
3 changes: 1 addition & 2 deletions examples/mittleman.jl
Original file line number Diff line number Diff line change
Expand Up @@ -158,7 +158,7 @@ code = ()
title!("#branches = $(size(getBranch(diagram, code)))")
# xlims!(0.01, 0.065, ylims=(2.5,6.5))

plot(getBranchesFromBP(diagram, 2); plotfold = false, legend = false, vars = (:param, :x))
plot(getBranchesFromBP(diagram, 2); plotfold = false, legend = false, vars = (:param, :n2))

getBranch(diagram, (2,1)) |> plot

Expand Down Expand Up @@ -212,7 +212,6 @@ res = BK.continuation(br, 2,

plot(res..., br ;plotfold= false)


δp = 0.005
deflationOp = DeflationOperator(2, 1.0, [zeros(2)])
success = [0]
Expand Down
2 changes: 1 addition & 1 deletion src/BifurcationKit.jl
Original file line number Diff line number Diff line change
Expand Up @@ -203,7 +203,7 @@ module BifurcationKit
export NewtonPar, newton, newtonDeflated, newtonPALC, newtonFold, newtonHopf, newtonBordered, NonLinearSolution

# continuation methods
export ContinuationPar, ContResult, GenericBifPoint, continuation, continuation!, continuationFold, continuationHopf, continuationPOTrap, continuationBordered, eigenvec, eigenvals, getSolx, getSolp, bifurcation_points
export ContinuationPar, ContResult, GenericBifPoint, continuation, continuation!, continuationFold, continuationHopf, continuationPOTrap, continuationBordered, eigenvec, eigenvals, getSolx, getSolp, bifurcation_points, SpecialPoint

# events
export ContinuousEvent, DiscreteEvent, PairOfEvents, SetOfEvents, SaveAtEvent, FoldDetectEvent, BifDetectEvent
Expand Down
20 changes: 15 additions & 5 deletions src/BifurcationPoints.jl
Original file line number Diff line number Diff line change
Expand Up @@ -8,24 +8,34 @@ istranscritical(bp::AbstractBranchPoint) = false
"""
$(TYPEDEF)
Structure to record specials point on the curve. There are two types of specials point that are recorded in this structure: bifurcation points and events (see https://bifurcationkit.github.io/BifurcationKitDocs.jl/dev/EventCallback/).
Structure to record special points on a curve. There are two types of special points that are recorded in this structure: bifurcation points and events (see https://bifurcationkit.github.io/BifurcationKitDocs.jl/dev/EventCallback/).
$(TYPEDFIELDS)
"""
@with_kw struct SpecialPoint{T, Tp, Tv, Tvτ} <: AbstractBifurcationPoint
"Description of the special point. In case of Events, this field records the user passed named to the event, or the default `:userD`, `:userC`. In case of bifurcation points, it can be the following:
"Description of the special points. In case of Events, this field records the user passed named to the event, or the default `:userD`, `:userC`. In case of bifurcation points, it can be one of the following:
- :bp Bifurcation point, simple eigenvalue crossing the imaginary axis
- :fold Fold point
- :hopf Hopf point
- :nd not documented bifurcation point. Detected by multiple eigenvalues crossing. Generally occurs in problems with symmetries or in cases where the continuation step size is too large and merge two different bifurcation points.
- :cusp Cusp point
- :gh Generalized Hopf point
- :gh Generalized Hopf point (also called Bautin point)
- :bt Bogdanov-Takens point
- :zh Zero-Hopf point
- :hh Hopf-Hopf point
- :ns Neimark-Sacker point
- :pd Period-doubling point
- :R1 Strong resonance 1:1 of periodic orbits
- :R2 Strong resonance 1:2 of periodic orbits
- :R3 Strong resonance 1:3 of periodic orbits
- :R4 Strong resonance 1:4 of periodic orbits
- :foldFlip Fold / Flip of periodic orbits
- :foldNS Fold / Neimark-Sacker of periodic orbits
- :pdNS Period-Doubling / Neimark-Sacker of periodic orbits
- :gpd Generalized Period-Doubling of periodic orbits
- :nsns Double Neimark-Sacker of periodic orbits
- :ch Chenciner bifurcation of periodic orbits
"
type::Symbol = :none

Expand Down Expand Up @@ -62,7 +72,7 @@ $(TYPEDFIELDS)
"Precision in the location of the special point"
precision::T = -1

"Interval containing the special point"
"Interval parameter containing the special point"
interval::Tuple{T, T} = (0, 0)
end

Expand Down Expand Up @@ -98,7 +108,7 @@ SpecialPoint(it::ContIterable, state::ContState, type::Symbol, status::Symbol, i
function _show(io::IO, bp::SpecialPoint, ii::Int, p::String = "p")
if bp.type == :none ; return; end
if bp.type == :endpoint
@printf(io, "- #%3i, \033[1m%5s\033[0m at %s ≈ %+4.8f, step = %3i\n", ii, "endpoint", p, bp.param, bp.step)
@printf(io, "- #%3i, \033[1m%5s\033[0m at %s ≈ %+4.8f, step = %3i\n", ii, "endpoint", p, bp.param, bp.step)
return
end
if bp.status == :converged
Expand Down
8 changes: 4 additions & 4 deletions src/Problems.jl
Original file line number Diff line number Diff line change
Expand Up @@ -17,14 +17,13 @@ $(TYPEDEF)
Structure to hold the vector field and its derivatives. It should rarely be called directly. Also, in essence, it is very close to `SciMLBase.ODEFunction`.
## Fields
$(TYPEDFIELDS)
## Methods
- `residual(pb::BifFunction, x, p)` calls `pb.F(x,p)`
- `jacobian(pb::BifFunction, x, p)` calls `pb.J(x, p)`
- `jacobian(pb::BifFunction, x, p)` calls `pb.J(x, p)`
- `dF(pb::BifFunction, x, p, dx)` calls `pb.dF(x,p,dx)`
- etc
"""
Expand Down Expand Up @@ -94,6 +93,7 @@ $(TYPEDFIELDS)
## Methods
- `reMake(pb; kwargs...)` modify a bifurcation problem
- `getu0(pb)` calls `pb.u0`
- `getParams(pb)` calls `pb.params`
- `getLens(pb)` calls `pb.lens`
Expand All @@ -104,8 +104,8 @@ $(TYPEDFIELDS)
- `isSymmetric(pb)` calls `isSymmetric(pb.prob)`
## Constructors
- `BifurcationProblem(F, u0, params, lens; J, Jᵗ, d2F, d3F, kwargs...)` and `kwargs` are the fields above. You can pass your own jacobian with `J` (see [`BifFunction`](@ref) for description of the jacobian function) and jacobian adjoint with `Jᵗ`. For example, this can be used to provide finite differences based jacobian.
- `BifurcationProblem(F, u0, params, lens)` all derivatives are computed using ForwardDiff.
- `BifurcationProblem(F, u0, params, lens; J, Jᵗ, d2F, d3F, kwargs...)` and `kwargs` are the fields above. You can pass your own jacobian with `J` (see [`BifFunction`](@ref) for description of the jacobian function) and jacobian adjoint with `Jᵗ`. For example, this can be used to provide finite differences based jacobian using `BifurcationKit.finiteDifferences`.
"""
struct $op{Tvf, Tu, Tp, Tl <: Lens, Tplot, Trec} <: AbstractAllJetBifProblem
Expand Down
35 changes: 25 additions & 10 deletions src/Results.jl
Original file line number Diff line number Diff line change
Expand Up @@ -23,12 +23,13 @@ $(TYPEDFIELDS)
# Associated methods
- `length(br)` number of the continuation steps
- `show(br)` display information about the branch
- `eigenvals(br, ind)` returns the eigenvalues for the ind-th continuation step
- `eigenvec(br, ind, indev)` returns the indev-th eigenvector for the ind-th continuation step
- `br[k+1]` gives information about the k-th step. A typical run yields the following
- `getNormalForm(br, ind)` compute the normal form of the ind-th points in `br.specialpoint`
- `getLens(br)` return the parameter axis used for the branch
- `getLenses(br)` return the parameter two axis used for the branch when 2 parameters continuation is used (Fold, Hopf, NS, PD)
- `br[k+1]` gives information about the k-th step. A typical run yields the following
```
julia> br[1]
(x = 0.0, param = 0.1, itnewton = 0, itlinear = 0, ds = -0.01, θ = 0.5, n_unstable = 2, n_imag = 2, stable = false, step = 0, eigenvals = ComplexF64[0.1 - 1.0im, 0.1 + 1.0im], eigenvecs = ComplexF64[0.7071067811865475 - 0.0im 0.7071067811865475 + 0.0im; 0.0 + 0.7071067811865475im 0.0 - 0.7071067811865475im])
Expand All @@ -53,14 +54,17 @@ julia> br.param
- `getParams(br)` Parameters passed to continuation and used in the equation `F(x, par) = 0`.
- `setParam(br, p0)` set the parameter value `p0` according to `::Lens` for the parameters of the problem `br.prob`
- `getLens(br)` get the lens used for the computation of the branch
- `continuation(br, ind)` performs automatic branch switching (aBS) from ind-th bifurcation point. Typically branching from equilibrium to equilibrium, or periodic orbit to periodic orbit.
- `continuation(br, ind, lens2)` performs two parameters `(getLens(br), lens2)` continuation of the ind-th bifurcation point.
- `continuation(br, ind, probPO::AbstractPeriodicOrbitProblem)` performs aBS from ind-th bifurcation point (which must be a Hopf bifurcation point) to branch of periodic orbits.
"""
@with_kw_noshow struct ContResult{Tkind <: AbstractContinuationKind, Tbr, Teigvals, Teigvec, Biftype, Tsol, Tparc, Tprob, Talg} <: AbstractResult{Tkind, Tprob}
"holds the low-dimensional information about the branch. More precisely, `branch[i+1]` contains the following information `(recordFromSolution(u, param), param, itnewton, itlinear, ds, θ, n_unstable, n_imag, stable, step)` for each continuation step `i`.\n
- `itnewton` number of Newton iterations
- `itlinear` total number of linear iterations during corrector
- `itlinear` total number of linear iterations during newton (corrector)
- `n_unstable` number of eigenvalues with positive real part for each continuation step (to detect stationary bifurcation)
- `n_imag` number of eigenvalues with positive real part and non zero imaginary part for each continuation step (to detect Hopf bifurcation).
- `stable` stability of the computed solution for each continuation step. Hence, `stable` should match `eig[step]` which corresponds to `branch[k]` for a given `k`.
- `n_imag` number of eigenvalues with positive real part and non zero imaginary part at current continuation step (useful to detect Hopf bifurcation).
- `stable` stability of the computed solution for each continuation step. Hence, `stable` should match `eig[step]` which corresponds to `branch[k]` for a given `k`.
- `step` continuation step (here equal `i`)"
branch::StructArray{Tbr}

Expand All @@ -70,16 +74,16 @@ julia> br.param
"Vector of solutions sampled along the branch. This is set by the argument `saveSolEveryStep::Int64` (default 0) in [`ContinuationPar`](@ref)."
sol::Tsol

"The parameters used for the call to `continuation` which produced this branch. Must be a ContinationPar"
"The parameters used for the call to `continuation` which produced this branch. Must be a [`ContinuationPar`](@ref)"
contparams::Tparc

"Type of solutions computed in this branch."
kind::Tkind = EquilibriumCont()

"Structure associated to the functional, useful for branch switching. For example, when computing periodic orbits, the functional `PeriodicOrbitTrapProblem`, `ShootingProblem`... will be saved here."
"Bifurcation problem used to compute the branch, useful for branch switching. For example, when computing periodic orbits, the functional `PeriodicOrbitTrapProblem`, `ShootingProblem`... will be saved here."
prob::Tprob = nothing

"A vector holding the set of detected bifurcation points. See [`SpecialPoint`](@ref) for a description of the fields."
"A vector holding the set of detected bifurcation points. See [`SpecialPoint`](@ref) for a list of special points."
specialpoint::Vector{Biftype}

"Continuation algorithm used for the computation of the branch"
Expand Down Expand Up @@ -115,11 +119,22 @@ getfirstusertype(br::AbstractBranchResult) = keys(br.branch[1])[1]
setParam(br::AbstractBranchResult, p0) = setParam(br.prob, p0)
Base.getindex(br::ContResult, k::Int) = (br.branch[k]..., eigenvals = haseigenvalues(br) ? br.eig[k].eigenvals : nothing, eigenvecs = haseigenvector(br) ? br.eig[k].eigenvecs : nothing)
Base.lastindex(br::ContResult) = length(br)

"""
$(SIGNATURES)
Return the solution for the ind-th point stored in br.sol
"""
@inline function getSolx(br::ContResult, ind::Int)
@assert hassolution(br) "You did not record the solution in the branch. Please set `saveSolEveryStep` in `ContinuationPar`"
return br.sol[ind].x
end

"""
$(SIGNATURES)
Return the parameter for the ind-th point stored in br.sol
"""
@inline function getSolp(br::ContResult, ind::Int)
@assert hassolution(br) "You did not record the solution in the branch. Please set `saveSolEveryStep` in `ContinuationPar`"
return br.sol[ind].p
Expand All @@ -143,9 +158,9 @@ Return the eigenvalues of the ind-th continuation step. `verbose` is used to tel
function eigenvals(br::AbstractBranchResult, ind::Int, verbose::Bool = false)
@assert br.eig[ind+1].step == ind "Error in indexing eigenvalues. Please open an issue on the website."
if verbose
println("--> For ", getLensSymbol(br), " = ", br.branch[ind].param)
println("--> There are ", br.branch[ind].n_unstable, " unstable eigenvalues")
println("--> Eigenvalues for continuation step ", br.eig[ind+1].step)
println("──> For ", getLensSymbol(br), " = ", br.branch[ind].param)
println("──> There are ", br.branch[ind].n_unstable, " unstable eigenvalues")
println("──> Eigenvalues for continuation step ", br.eig[ind+1].step)
end
~br.eig[ind+1].converged && @error "Eigen solver did not converged on the step!!"
br.eig[ind+1].eigenvals
Expand Down
8 changes: 4 additions & 4 deletions src/codim2/NormalForms.jl
Original file line number Diff line number Diff line change
Expand Up @@ -545,15 +545,15 @@ function bogdanovTakensNormalForm(_prob,

# Al-Hdaibat, B., W. Govaerts, Yu. A. Kuznetsov, and H. G. E. Meijer. “Initialization of Homoclinic Solutions near Bogdanov--Takens Points: Lindstedt--Poincaré Compared with Regular Perturbation Method.” SIAM Journal on Applied Dynamical Systems 15, no. 2 (January 2016): 952–80. https://doi.org/10.1137/15M1017491.
###########################
vext = real.(ζs[1])
vr = real.(ζs[1])
Lᵗ = hasAdjoint(prob_vf) ? jad(prob_vf, x0, parbif) : transpose(L)
_λ★, _ev★, _ = eigsolver(Lᵗ, nev)
Ivp = sortperm(_λ★, by = abs)
# in case the prob is HopfMA, we real it
zerov = real.(prob_ma.zero)
wext = real.(geteigenvector(eigsolver, _ev★, Ivp[1]))
q0, = bls(L, wext, vext, zero(𝒯), zerov, one(𝒯))
p1, = bls(Lᵗ, vext, wext, zero(𝒯), zerov, one(𝒯))
vl = real.(geteigenvector(eigsolver, _ev★, Ivp[1]))
q0, = bls(L, vl, vr, zero(𝒯), zerov, one(𝒯))
p1, = bls(Lᵗ, vr, vl, zero(𝒯), zerov, one(𝒯))
q1, = bls(L, p1, q0, zero(𝒯), q0, zero(𝒯))
p0, = bls(Lᵗ, q0, p1, zero(𝒯), p1, zero(𝒯))
# we want
Expand Down
Loading

0 comments on commit 365c0b5

Please sign in to comment.