bifurcationkit / BifurcationKit.jl

A Julia package to perform Bifurcation Analysis
https://bifurcationkit.github.io/BifurcationKitDocs.jl/stable
Other
306 stars 37 forks source link

UndefVarError: getTrajectory not define (when running tutorial) #59

Closed TorkelE closed 2 years ago

TorkelE commented 2 years ago

I am running the tutorial, trying to figure out how to plot periodic orbits in ODE bifurcation diagrams.

I am trying the tutorial in: https://rveltz.github.io/BifurcationKit.jl/dev/tutorials/ode/tutorialsODE/#Neural-mass-equation-(Hopf-aBS) however, the section on "Branch of periodic orbits with Trapezoid method" gives an error:

UndefVarError: getTrajectory not defined

Stacktrace:
  [1] (::var"#5#8")(x::Vector{Float64}, p::NamedTuple{(:prob, :p), Tuple{PeriodicOrbitTrapProblem{typeof(TMvf), typeof(dTMvf), Nothing, Nothing, Nothing, Vector{Float64}, DefaultLS, BifurcationKit.TimeMesh{Int64}, Nothing}, Float64}})
    @ Main ./In[5]:7
  [2] #808#809
    @ ~/.julia/packages/BifurcationKit/eoTLY/src/periodicorbit/PeriodicOrbitUtils.jl:81 [inlined]
  [3] #808
    @ ~/.julia/packages/BifurcationKit/eoTLY/src/periodicorbit/PeriodicOrbitUtils.jl:81 [inlined]
  [4] ContResult(it::ContIterable{PeriodicOrbitTrapProblem{typeof(TMvf), typeof(dTMvf), Nothing, Nothing, Nothing, Vector{Float64}, DefaultLS, BifurcationKit.TimeMesh{Int64}, Nothing}, BifurcationKit.var"#753#758"{PeriodicOrbitTrapProblem{typeof(TMvf), typeof(dTMvf), Nothing, Nothing, Nothing, Vector{Float64}, DefaultLS, BifurcationKit.TimeMesh{Int64}, Nothing}}, Vector{Float64}, NamedTuple{(:α, :τ, :J, :E0, :τD, :U0, :τF, :τS), NTuple{8, Float64}}, Setfield.PropertyLens{:E0}, Float64, BifurcationKit.FloquetWrapperLS{DefaultLS}, FloquetQaD{DefaultEig{typeof(abs)}}, SecantPred, BorderingBLS{BifurcationKit.FloquetWrapperLS{DefaultLS}, Float64}, BifurcationKit.var"#818#822"{BifurcationKit.var"#818#819#823"{PeriodicOrbitTrapProblem{typeof(TMvf), typeof(dTMvf), Nothing, Nothing, Nothing, Vector{Float64}, DefaultLS, BifurcationKit.TimeMesh{Int64}, Nothing}, var"#6#9"}}, BifurcationKit.var"#808#812"{BifurcationKit.var"#808#809#813"{PeriodicOrbitTrapProblem{typeof(TMvf), typeof(dTMvf), Nothing, Nothing, Nothing, Vector{Float64}, DefaultLS, BifurcationKit.TimeMesh{Int64}, Nothing}, var"#5#8"}}, typeof(norminf), BifurcationKit.DotTheta{BifurcationKit.var"#169#171", BifurcationKit.var"#170#172"}, BifurcationKit.var"#800#804"{BifurcationKit.var"#800#801#805"{PeriodicOrbitTrapProblem{typeof(TMvf), typeof(dTMvf), Nothing, Nothing, Nothing, Vector{Float64}, DefaultLS, BifurcationKit.TimeMesh{Int64}, Nothing}, Int64}}, typeof(BifurcationKit.cbDefault), Nothing, String}, state::ContState{BorderedArray{Vector{Float64}, Float64}, Float64, Vector{ComplexF64}, Matrix{ComplexF64}, Tuple{Nothing, Nothing}})
    @ BifurcationKit ~/.julia/packages/BifurcationKit/eoTLY/src/Continuation.jl:238
  [5] continuation(it::ContIterable{PeriodicOrbitTrapProblem{typeof(TMvf), typeof(dTMvf), Nothing, Nothing, Nothing, Vector{Float64}, DefaultLS, BifurcationKit.TimeMesh{Int64}, Nothing}, BifurcationKit.var"#753#758"{PeriodicOrbitTrapProblem{typeof(TMvf), typeof(dTMvf), Nothing, Nothing, Nothing, Vector{Float64}, DefaultLS, BifurcationKit.TimeMesh{Int64}, Nothing}}, Vector{Float64}, NamedTuple{(:α, :τ, :J, :E0, :τD, :U0, :τF, :τS), NTuple{8, Float64}}, Setfield.PropertyLens{:E0}, Float64, BifurcationKit.FloquetWrapperLS{DefaultLS}, FloquetQaD{DefaultEig{typeof(abs)}}, SecantPred, BorderingBLS{BifurcationKit.FloquetWrapperLS{DefaultLS}, Float64}, BifurcationKit.var"#818#822"{BifurcationKit.var"#818#819#823"{PeriodicOrbitTrapProblem{typeof(TMvf), typeof(dTMvf), Nothing, Nothing, Nothing, Vector{Float64}, DefaultLS, BifurcationKit.TimeMesh{Int64}, Nothing}, var"#6#9"}}, BifurcationKit.var"#808#812"{BifurcationKit.var"#808#809#813"{PeriodicOrbitTrapProblem{typeof(TMvf), typeof(dTMvf), Nothing, Nothing, Nothing, Vector{Float64}, DefaultLS, BifurcationKit.TimeMesh{Int64}, Nothing}, var"#5#8"}}, typeof(norminf), BifurcationKit.DotTheta{BifurcationKit.var"#169#171", BifurcationKit.var"#170#172"}, BifurcationKit.var"#800#804"{BifurcationKit.var"#800#801#805"{PeriodicOrbitTrapProblem{typeof(TMvf), typeof(dTMvf), Nothing, Nothing, Nothing, Vector{Float64}, DefaultLS, BifurcationKit.TimeMesh{Int64}, Nothing}, Int64}}, typeof(BifurcationKit.cbDefault), Nothing, String})
    @ BifurcationKit ~/.julia/packages/BifurcationKit/eoTLY/src/Continuation.jl:482
  [6] continuation(Fhandle::PeriodicOrbitTrapProblem{typeof(TMvf), typeof(dTMvf), Nothing, Nothing, Nothing, Vector{Float64}, DefaultLS, BifurcationKit.TimeMesh{Int64}, Nothing}, Jhandle::Function, x0::Vector{Float64}, par::NamedTuple{(:α, :τ, :J, :E0, :τD, :U0, :τF, :τS), NTuple{8, Float64}}, lens::Setfield.PropertyLens{:E0}, contParams::ContinuationPar{Float64, BifurcationKit.FloquetWrapperLS{DefaultLS}, FloquetQaD{DefaultEig{typeof(abs)}}}, linearAlgo::BorderingBLS{BifurcationKit.FloquetWrapperLS{DefaultLS}, Float64}; bothside::Bool, kwargs::Base.Pairs{Symbol, Any, NTuple{7, Symbol}, NamedTuple{(:recordFromSolution, :linearPO, :verbosity, :plot, :plotSolution, :normC, :finaliseSolution), Tuple{BifurcationKit.var"#808#812"{BifurcationKit.var"#808#809#813"{PeriodicOrbitTrapProblem{typeof(TMvf), typeof(dTMvf), Nothing, Nothing, Nothing, Vector{Float64}, DefaultLS, BifurcationKit.TimeMesh{Int64}, Nothing}, var"#5#8"}}, Symbol, Int64, Bool, BifurcationKit.var"#818#822"{BifurcationKit.var"#818#819#823"{PeriodicOrbitTrapProblem{typeof(TMvf), typeof(dTMvf), Nothing, Nothing, Nothing, Vector{Float64}, DefaultLS, BifurcationKit.TimeMesh{Int64}, Nothing}, var"#6#9"}}, typeof(norminf), BifurcationKit.var"#800#804"{BifurcationKit.var"#800#801#805"{PeriodicOrbitTrapProblem{typeof(TMvf), typeof(dTMvf), Nothing, Nothing, Nothing, Vector{Float64}, DefaultLS, BifurcationKit.TimeMesh{Int64}, Nothing}, Int64}}}}})
    @ BifurcationKit ~/.julia/packages/BifurcationKit/eoTLY/src/Continuation.jl:0
  [7] continuationPOTrap(prob::PeriodicOrbitTrapProblem{typeof(TMvf), typeof(dTMvf), Nothing, Nothing, Nothing, Vector{Float64}, DefaultLS, BifurcationKit.TimeMesh{Int64}, Nothing}, orbitguess::Vector{Float64}, par::NamedTuple{(:α, :τ, :J, :E0, :τD, :U0, :τF, :τS), NTuple{8, Float64}}, lens::Setfield.PropertyLens{:E0}, contParams::ContinuationPar{Float64, DefaultLS, DefaultEig{typeof(real)}}, linearAlgo::BorderingBLS{DefaultLS, Float64}; jacobianPO::Symbol, updateSectionEveryStep::Int64, kwargs::Base.Pairs{Symbol, Any, NTuple{6, Symbol}, NamedTuple{(:recordFromSolution, :linearPO, :verbosity, :plot, :plotSolution, :normC), Tuple{var"#5#8", Symbol, Int64, Bool, var"#6#9", typeof(norminf)}}})
    @ BifurcationKit ~/.julia/packages/BifurcationKit/eoTLY/src/periodicorbit/PeriodicOrbitTrapeze.jl:916
  [8] #continuation#763
    @ ~/.julia/packages/BifurcationKit/eoTLY/src/periodicorbit/PeriodicOrbitTrapeze.jl:978 [inlined]
  [9] continuation(F::typeof(TMvf), dF::typeof(dTMvf), d2F::BifurcationKit.var"#d2f#284"{BifurcationKit.var"#d1f#282"{typeof(TMvf)}}, d3F::BifurcationKit.var"#d3f#286"{BifurcationKit.var"#d2f#284"{BifurcationKit.var"#d1f#282"{typeof(TMvf)}}}, br::ContResult{NamedTuple{(:E, :x, :u, :param, :itnewton, :itlinear, :ds, :θ, :n_unstable, :n_imag, :stable, :step), Tuple{Float64, Float64, Float64, Float64, Int64, Int64, Float64, Float64, Int64, Int64, Bool, Int64}}, Vector{ComplexF64}, Matrix{ComplexF64}, BifurcationKit.SpecialPoint{Float64, NamedTuple{(:E, :x, :u), Tuple{Float64, Float64, Float64}}, Vector{Float64}}, Nothing, ContinuationPar{Float64, DefaultLS, DefaultEig{typeof(real)}}, Nothing, NamedTuple{(:α, :τ, :J, :E0, :τD, :U0, :τF, :τS), NTuple{8, Float64}}, Setfield.PropertyLens{:E0}}, ind_bif::Int64, _contParams::ContinuationPar{Float64, DefaultLS, DefaultEig{typeof(real)}}, prob::PeriodicOrbitTrapProblem{Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, DefaultLS, BifurcationKit.TimeMesh{Int64}, Nothing}; Jᵗ::Nothing, δ::Float64, δp::Nothing, ampfactor::Int64, usedeflation::Bool, nev::Int64, updateSectionEveryStep::Int64, kwargs::Base.Pairs{Symbol, Any, NTuple{6, Symbol}, NamedTuple{(:linearPO, :verbosity, :plot, :recordFromSolution, :plotSolution, :normC), Tuple{Symbol, Int64, Bool, var"#5#8", var"#6#9", typeof(norminf)}}})
    @ BifurcationKit ~/.julia/packages/BifurcationKit/eoTLY/src/periodicorbit/PeriodicOrbits.jl:410
 [10] top-level scope
    @ In[5]:21
 [11] eval
    @ ./boot.jl:373 [inlined]
 [12] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
    @ Base ./loading.jl:1196

I get this after copying (and running) all the code in the tutorial

using Revise, ForwardDiff, Parameters, Setfield, Plots, LinearAlgebra
using BifurcationKit
const BK = BifurcationKit

# sup norm
norminf(x) = norm(x, Inf)

# vector field
function TMvf!(dz, z, p, t)
    @unpack J, α, E0, τ, τD, τF, U0 = p
    E, x, u = z
    SS0 = J * u * x * E + E0
    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
end

# out of place method
TMvf(z, p) = TMvf!(similar(z), z, p, 0)

# we group the differentials together
dTMvf(z,p) = ForwardDiff.jacobian(x -> TMvf(x,p), z)
jet = BK.getJet(TMvf, dTMvf)

# parameter values
par_tm = (α = 1.5, τ = 0.013, J = 3.07, E0 = -2.0, τD = 0.200, U0 = 0.3, τF = 1.5, τS = 0.007)

# initial condition
z0 = [0.238616, 0.982747, 0.367876]

# continuation options
opts_br = ContinuationPar(pMin = -10.0, pMax = -0.9,
    # parameters to have a smooth result
    ds = 0.04, dsmax = 0.05,
    # this is to detect bifurcation points precisely with bisection
    detectBifurcation = 3,
    # Optional: bisection options for locating bifurcations
    nInversion = 8, maxBisectionSteps = 25, nev = 3)

# continuation of equilibria
br, = continuation(TMvf, dTMvf, z0, par_tm, (@lens _.E0), opts_br;
    recordFromSolution = (x, p) -> (E = x[1], x = x[2], u = x[3]),
    tangentAlgo = BorderedPred(),
    plot = true, normC = norminf)

scene = plot(br, plotfold=false, markersize=3, legend=:topleft)

br

# newton parameters
optn_po = NewtonPar(verbose = true, tol = 1e-8,  maxIter = 10)

# continuation parameters
opts_po_cont = ContinuationPar(dsmax = 0.1, ds= -0.0001, dsmin = 1e-4, pMax = 0., pMin=-5.,
    maxSteps = 110, newtonOptions = (@set optn_po.tol = 1e-7),
    nev = 3, precisionStability = 1e-8, detectBifurcation = 3, plotEveryStep = 10, saveSolEveryStep=1)

# arguments for periodic orbits
args_po = ( recordFromSolution = (x, p) -> (xtt = BK.getTrajectory(p.prob, x, p.p);
        return (max = maximum(xtt.u[1,:]),
                min = minimum(xtt.u[1,:]),
                period = x[end])),
    plotSolution = (x, p; k...) -> begin
        xtt = BK.getTrajectory(p.prob, x, p.p)
        plot!(xtt.t, xtt.u[1,:]; label = "E", k...)
        plot!(xtt.t, xtt.u[2,:]; label = "x", k...)
        plot!(xtt.t, xtt.u[3,:]; label = "u", k...)
        plot!(br,subplot=1, putbifptlegend = false)
        end,
    normC = norminf)

Mt = 200 # number of time sections
    br_potrap, utrap = continuation(jet...,
    # we want to branch form the 4th bif. point
    br, 4, opts_po_cont,
    # we want to use the Trapeze method to locate PO
    PeriodicOrbitTrapProblem(M = Mt);
    # this jacobian is specific to ODEs
    # it is computed using AD and
    # updated inplace
    linearPO = :Dense,
    # regular continuation options
    verbosity = 2,  plot = true,
    args_po...)

scene = plot(br, br_potrap, markersize = 3)
plot!(scene, br_potrap.param, br_potrap.min, label = "")

Indeed, just trying to access the getTrajectory function:

BK.getTrajectory

gives the error. I checked, and I have the latest version of BifurcationKit ([0f109fa4] BifurcationKit v0.1.11).

rveltz commented 2 years ago

strange. Did you try master?

Also, note that the docs are automatically generated.

TorkelE commented 2 years ago

No, I will attempt this.

TorkelE commented 2 years ago

I tried on Master but didn't work. However, having the master file I could do a search for getTrajectory and the change log revealed it had been renamed getPeriodicOrbit. Changing getTrajectory to getPeriodicOrbit solved the issue. Thanks :)

rveltz commented 2 years ago

Sorry I did not see that this error would not be caught by automated docs.

rveltz commented 2 years ago

In this link, there is getPeriodic

https://bifurcationkit.github.io/BifurcationKitDocs.jl/dev/tutorials/ode/tutorialsODE/#Neural-mass-equation-(Hopf-aBS)

it seems you are on old docs