diff --git a/News.md b/News.md index 87152caa..627613f2 100644 --- a/News.md +++ b/News.md @@ -6,6 +6,9 @@ All notable changes to this project will be documented in this file (hopefully). ## Near future - Setfield.jl has been changed to Accessors.jl +## [0.3.7] +- remove `Requires.jl` and use extensions. This requires julia>=1.9 + ## [0.3.5] - add field `save_solution` to `BifurcationProblem`. Allows to save problem state along with the current solution. Useful for periodic orbits for example where we can save the phase condition and the mesh when the latter is adapted. diff --git a/Project.toml b/Project.toml index 35297a0d..e9bff7de 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "BifurcationKit" uuid = "0f109fa4-8a5d-4b75-95aa-f515264e7665" authors = ["Romain Veltz "] -version = "0.3.6" +version = "0.3.7" [deps] ArnoldiMethod = "ec485272-7323-5ecc-a04f-4719b315124d" @@ -22,7 +22,6 @@ Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01" RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd" Reexport = "189a3867-3050-52da-a836-e630ba90ab69" -Requires = "ae029012-a4dd-5104-9daa-d747884805df" SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462" Setfield = "efcf1570-3423-57d1-acb7-fd33fddbac46" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" @@ -42,14 +41,23 @@ LinearMaps = "2.7, ^3.0" OrdinaryDiffEq = "^6.33" Parameters = "0.12.3" PreallocationTools = "0.4.1" -RecipesBase = "1.0, 1.1" +RecipesBase = "1.0, 1.1, 1.3" RecursiveArrayTools = "^2.3, ^2.4, ^2.8, ^2.9, ^3" Reexport = "1.2.2" -Requires = "1.1.2" SciMLBase = "^1, ^2" Setfield = "0.5.0, 0.7.0, 0.8.0, ^1.1" StructArrays = "0.4, ^0.6" -julia = "1.5" +julia = "1.9" + +[weakdeps] +JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819" +Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" +GLMakie = "e9467ef8-e4e7-5192-8a1a-b1aee30e663a" + +[extensions] +JLD2Ext = "JLD2" +PlotsExt = "Plots" +GLMakieExt = "GLMakie" [extras] AbstractDifferentiation = "c29ec348-61ec-40c8-8164-b8c60e9d9f3d" diff --git a/examples/COModel.jl b/examples/COModel.jl index 6457d57a..058e6393 100644 --- a/examples/COModel.jl +++ b/examples/COModel.jl @@ -27,7 +27,7 @@ br = @time continuation(prob, PALC(), opts_br; normC = norminf, bothside = true) -BK.plot(br, plotfold=false, )#markersize=4, legend=:topright, ylims=(0,0.16)) +BK.plot(br)#markersize=4, legend=:topright, ylims=(0,0.16)) #################################################################################################### # periodic orbits function plotSolution(x, p; k...) @@ -42,7 +42,7 @@ end function plotSolution(ax, x, p; ax1 = nothing, k...) @info "plotsol Makie" xtt = BK.get_periodic_orbit(p.prob, x, p.p) - lines!(ax1, br; putspecialptlegend = false) + lines!(ax1, br) lines!(ax, xtt.t, xtt[1,:]; k...) lines!(ax, xtt.t, xtt[2,:]; k...) lines!(ax, xtt.t, xtt[3,:]; k...) @@ -99,7 +99,7 @@ brpo = @time continuation(br, 5, opts_po_cont, args_po... ) -BK.plot(br, brpo, putspecialptlegend = false, branchlabel = ["eq","max"]) +BK.plot(br, brpo, branchlabel = ["eq","max"]) xlims!((1.037, 1.055)) scatter!(br) plot!(brpo.param, brpo.min, label = "min", xlims = (1.037, 1.055)) diff --git a/ext/GLMakieExt/GLMakieExt.jl b/ext/GLMakieExt/GLMakieExt.jl new file mode 100644 index 00000000..1762a981 --- /dev/null +++ b/ext/GLMakieExt/GLMakieExt.jl @@ -0,0 +1,29 @@ +module GLMakieExt + using GLMakie, BifurcationKit + import BifurcationKit: plot, + plot!, + hasbranch, + plot_branch_cont, + plot_periodic_potrap, + plot_periodic_shooting!, + plot_periodic_shooting, + plot_eigenvals, + compute_eigenelements, + get_lens_symbol, + _getfirstusertype, + ContResult, + DCResult, + AbstractBranchResult, + get_plot_vars, + get_axis_labels, + _hasstability, + filter_bifurcations, + get_color, + colorbif, + get_plot_backend, + BK_Makie + + # TODO block precompilation + get_plot_backend() = BK_Makie() + include("plot.jl") +end diff --git a/src/plotting/RecipesMakie.jl b/ext/GLMakieExt/plot.jl similarity index 82% rename from src/plotting/RecipesMakie.jl rename to ext/GLMakieExt/plot.jl index 9bd77b3e..1921dfed 100644 --- a/src/plotting/RecipesMakie.jl +++ b/ext/GLMakieExt/plot.jl @@ -1,5 +1,5 @@ -# https://github.com/SciML/SciMLBase.jl/issues/427 -# https://github.com/VarLad/SciMLMakie.jl +using GLMakie: Point2f0 + function GLMakie.convert_arguments(::PointBased, contres::AbstractBranchResult, vars = nothing, applytoY = identity, applytoX = identity) ind1, ind2 = get_plot_vars(contres, vars) return ([Point2f0(i, j) for (i, j) in zip(map(applytoX, getproperty(contres.branch, ind1)), map(applytoY, getproperty(contres.branch, ind2)))],) @@ -35,7 +35,7 @@ function plot!(ax1, contres::AbstractBranchResult; end ax1.xlabel = xlab ax1.ylabel = ylab - + # display bifurcation points bifpt = filter(x -> (x.type != :none) && (x.type != :endpoint) && (plotfold || x.type != :fold) && (x.idx <= length(contres)-1), contres.specialpoint) if length(bifpt) >= 1 && plotspecialpoints #&& (ind1 == :param) @@ -43,64 +43,64 @@ function plot!(ax1, contres::AbstractBranchResult; bifpt = filterBifurcations(bifpt) end scatter!(ax1, - [applytoX(getproperty(contres[pt.idx], ind1)) for pt in bifpt], - [applytoY(getproperty(contres[pt.idx], ind2)) for pt in bifpt]; - marker = map(x -> (x.status == :guess) && (plotcirclesbif==false) ? :rect : :circle, bifpt), - markersize = 10, - color = map(x -> colorbif[x.type], bifpt), - ) + [applytoX(getproperty(contres[pt.idx], ind1)) for pt in bifpt], + [applytoY(getproperty(contres[pt.idx], ind2)) for pt in bifpt]; + marker = map(x -> (x.status == :guess) && (plotcirclesbif==false) ? :rect : :circle, bifpt), + markersize = 10, + color = map(x -> colorbif[x.type], bifpt), + ) end - + # add legend for bifurcation points if putspecialptlegend && length(bifpt) >= 1 bps = unique(x -> x.type, [pt for pt in bifpt if (pt.type != :none && (plotfold || pt.type != :fold))]) (length(bps) == 0) && return for pt in bps scatter!(ax1, - [applytoX(getproperty(contres[pt.idx], ind1))], - [applytoY(getproperty(contres[pt.idx], ind2))]; - color = get_color(pt.type), - markersize = 10, - label = "$(pt.type)") + [applytoX(getproperty(contres[pt.idx], ind1))], + [applytoY(getproperty(contres[pt.idx], ind2))]; + color = get_color(pt.type), + markersize = 10, + label = "$(pt.type)") end - # GLMakie.axislegend(ax1, merge = true, unique = true) + GLMakie.axislegend(ax1, merge = true, unique = true) end ax1 end function plot_branch_cont(contres::ContResult, - state, - iter, - plotuserfunction; - plotfold = false, - plotstability = true, - plotspecialpoints = true, - putspecialptlegend = true, - filterspecialpoints = false, - linewidthunstable = 1.0, - linewidthstable = 3.0linewidthunstable, - plotcirclesbif = true, - applytoY = identity, - applytoX = identity) + state, + iter, + plotuserfunction; + plotfold = false, + plotstability = true, + plotspecialpoints = true, + putspecialptlegend = true, + filterspecialpoints = false, + linewidthunstable = 1.0, + linewidthstable = 3.0linewidthunstable, + plotcirclesbif = true, + applytoY = identity, + applytoX = identity) sol = getsolution(state) if length(contres) == 0; return ; end - + # names for axis labels ind1, ind2 = get_plot_vars(contres, nothing) xlab, ylab = get_axis_labels(ind1, ind2, contres) - + # stability linewidth linewidth = linewidthunstable if _hasstability(contres) && plotstability linewidth = map(x -> isodd(x) ? linewidthstable : linewidthunstable, contres.stable) end - + fig = Figure(size = (1200, 700)) ax1 = fig[1:2, 1] = Axis(fig, xlabel = String(xlab), ylabel = String(ylab), tellheight = true) - + ax2 = fig[1, 2] = Axis(fig, xlabel = "step [$(state.step)]", ylabel = String(xlab)) lines!(ax2, contres.step, contres.param, linewidth = linewidth) - + if compute_eigenelements(iter) eigvals = contres.eig[end].eigenvals ax_ev = fig[3, 1:2] = Axis(fig, xlabel = "ℜ", ylabel = "ℑ") @@ -113,23 +113,23 @@ function plot_branch_cont(contres::ContResult, end lines!(ax_ev, [0, 0], [maxIm, minIm], color = :blue, linewidth = linewidthunstable) end - + # plot arrow to indicate the order of computation if length(contres) > 1 x = contres.branch[end].param y = getproperty(contres.branch,1)[end] u = contres.branch[end].param - contres.branch[end-1].param - v = getproperty(contres.branch,1)[end]-getproperty(contres.branch,1)[end-1] + v = getproperty(contres.branch,1)[end] - getproperty(contres.branch,1)[end-1] GLMakie.arrows!(ax1, [x], [y], [u], [v], color = :green, arrowsize = 20,) end - + plot!(ax1, contres; plotfold, plotstability, plotspecialpoints, putspecialptlegend, filterspecialpoints, linewidthunstable, linewidthstable, plotcirclesbif, applytoY, applytoX) - + if isnothing(plotuserfunction) == false ax_perso = fig[2, 2] = Axis(fig, tellheight = true) plotuserfunction(ax_perso, sol.u, sol.p; ax1 = ax1) end - + display(fig) fig end @@ -142,30 +142,25 @@ function plot(contres::AbstractBranchResult; kP...) fig = Figure() ax1 = fig[1, 1] = Axis(fig, xlabel = String(xlab), ylabel = String(ylab), tellheight = true) - - plot!(ax1, contres; kP...) + plot!(ax1, contres; kP...) display(fig) fig, ax1 end +plot(brdc::DCResult; kP...) = plot(brdc.branches...; kP...) + function plot(brs::AbstractBranchResult...; - branchlabel = fill("", length(brs)), - kP...) + branchlabel = ["$i" for i=1:length(brs)], + kP...) if length(brs) == 0; return ;end fig = Figure() ax1 = fig[1, 1] = Axis(fig) for (id, contres) in pairs(brs) - - ind1, ind2 = get_plot_vars(contres, nothing) - xlab, ylab = get_axis_labels(ind1, ind2, contres) - plot!(ax1, contres; branchlabel = branchlabel[id], kP...) - end - # GLMakie.axislegend(ax1, merge = true, unique = true) - + GLMakie.axislegend(ax1, merge = true, unique = true) display(fig) fig, ax1 end @@ -218,7 +213,7 @@ end function plot(bd::BifDiagNode; code = (), level = (-Inf, Inf), k...) if ~hasbranch(bd); return; end - fig = Figure(size = (1200, 700)) + fig = Figure() ax = fig[1, 1] = Axis(fig) _plot_bifdiag_makie!(ax, bd; code, level, k...) @@ -244,4 +239,4 @@ function _plot_bifdiag_makie!(ax, bd::Vector{BifDiagNode}; code = (), level = (- _plot_bifdiag_makie!(ax, b; code, level, k... ) end end -#################################################################################################### +#################################################################################################### \ No newline at end of file diff --git a/ext/JLD2Ext/JLD2Ext.jl b/ext/JLD2Ext/JLD2Ext.jl new file mode 100644 index 00000000..d3c15d54 --- /dev/null +++ b/ext/JLD2Ext/JLD2Ext.jl @@ -0,0 +1,5 @@ +module JLD2Ext + using JLD2, BifurcationKit + import BifurcationKit: save_to_file, AbstractContinuationIterable, ContResult + include("save.jl") +end diff --git a/ext/JLD2Ext/save.jl b/ext/JLD2Ext/save.jl new file mode 100644 index 00000000..766ab7e2 --- /dev/null +++ b/ext/JLD2Ext/save.jl @@ -0,0 +1,50 @@ +""" +Save solution / data in JLD2 file +- `filename` is for example "example.jld2" +- `sol` is the solution +- `p` is the parameter +- `i` is the index of the solution to be saved +""" +function save_to_file(iter::AbstractContinuationIterable, sol, p, i::Int64, br::ContResult) + if iter.contparams.save_to_file == false; return nothing; end + filename = iter.filename + # this allows to save two branches forward/backward in case + # bothside = true is passed to continuation + fd = iter.contparams.ds >=0 ? "fw" : "bw" + + # create a group in the JLD format + jldopen(filename*".jld2", "a+") do file + if haskey(file, "sol-$fd-$i") + delete!(file, "sol-$fd-$i") + end + mygroup = JLD2.Group(file, "sol-$fd-$i") + mygroup["sol"] = sol + mygroup["param"] = p + end + + jldopen(filename*"-branch.jld2", "a+") do file + if haskey(file, "branch"*fd) + delete!(file, "branch"*fd) + end + file["branch"*fd] = br + end +end + +# final save of branch, in case bothsided = true is used +function save_to_file(iter::AbstractContinuationIterable, br::ContResult) + if iter.contparams.save_to_file == false; return nothing; end + filename = iter.filename + + jldopen(filename*"-branch.jld2", "a+") do file + if haskey(file, "branchfw") + delete!(file, "branchfw") + end + if haskey(file, "branchbw") + delete!(file, "branchbw") + end + if haskey(file, "branch") + delete!(file, "branch") + end + file["branch"] = br + end +end \ No newline at end of file diff --git a/ext/PlotsExt/PlotsExt.jl b/ext/PlotsExt/PlotsExt.jl new file mode 100644 index 00000000..2aa15309 --- /dev/null +++ b/ext/PlotsExt/PlotsExt.jl @@ -0,0 +1,22 @@ +module PlotsExt + using Plots, BifurcationKit + import BifurcationKit: plot_branch_cont, + plot_periodic_potrap, + plot_periodic_shooting!, + plot_periodic_shooting, + plot_eigenvals, + compute_eigenelements, + get_lens_symbol, + _getfirstusertype, + ContResult, + DCResult, + AbstractBranchResult, + get_plot_vars, + get_axis_labels, + _hasstability, + filter_bifurcations, + get_color, + AbstractResult + + include("plot.jl") +end diff --git a/ext/PlotsExt/plot.jl b/ext/PlotsExt/plot.jl new file mode 100644 index 00000000..e9f60fba --- /dev/null +++ b/ext/PlotsExt/plot.jl @@ -0,0 +1,94 @@ +""" +Plot the branch of solutions during the continuation +""" +function plot_branch_cont(contres::ContResult, + state, + iter, + plotuserfunction) + sol = getsolution(state) + l = compute_eigenelements(iter) ? Plots.@layout([a{0.5w} [b; c]; e{0.2h}]) : Plots.@layout([a{0.5w} [b; c]]) + Plots.plot(layout = l ) + + Plots.plot!(contres ; + filterspecialpoints = true, + putspecialptlegend = false, + xlabel = get_lens_symbol(contres), + ylabel = _getfirstusertype(contres), + label = "", + plotfold = false, + subplot = 1) + + plotuserfunction(sol.u, sol.p; subplot = 3) + + # put arrow to indicate the order of computation + length(contres) > 1 && plot!([contres.branch[end-1:end].param], [getproperty(contres.branch,1)[end-1:end]], label = "", arrow = true, subplot = 1) + + if compute_eigenelements(iter) + eigvals = contres.eig[end].eigenvals + Plots.scatter!(real.(eigvals), imag.(eigvals), subplot=4, label = "", markerstrokewidth = 0, markersize = 3, color = :black, xlabel = "ℜ", ylabel = "ℑ") + # add stability boundary + maxIm = maximum(imag, eigvals) + minIm = minimum(imag, eigvals) + Plots.plot!([0, 0], [maxIm, minIm], subplot=4, label = "", color = :blue) + end + + Plots.plot!(contres; + vars = (:step, :param), + putspecialptlegend = false, + plotspecialpoints = false, + xlabel = "step [$(state.step)]", + ylabel = get_lens_symbol(contres), + label = "", + subplot = 2) |> display +end + +#################################################################################################### +# plotting function of the periodic orbits +function plot_periodic_potrap(x, M, Nx, Ny; ratio = 2, kwargs...) + @assert ratio > 0 "You need at least one component" + n = Nx*Ny + outpo = reshape(x[begin:end-1], ratio * n, M) + po = reshape(x[1:n,1], Nx, Ny) + rg = 2:6:M + for ii in rg + po = hcat(po, reshape(outpo[1:n,ii], Nx, Ny)) + end + heatmap!(po; color = :viridis, fill=true, xlabel = "space", ylabel = "space", kwargs...) + for ii in eachindex(rg) + Plots.plot!([ii*Ny, ii*Ny], [1, Nx]; color = :red, width = 3, label = "", kwargs...) + end +end + +function plot_periodic_potrap(outpof, n, M; ratio = 2) + @assert ratio > 0 "You need at least one component" + outpo = reshape(outpof[begin:end-1], ratio * n, M) + if ratio == 1 + Plots.heatmap(outpo[1:n,:]', ylabel="Time", color=:viridis) + else + Plots.plot( + Plots.heatmap(outpo[1:n,:]', ylabel="Time", color=:viridis), + Plots.heatmap(outpo[n+2:end,:]', color=:viridis)) + end +end +#################################################################################################### +function plot_periodic_shooting!(x, M; kwargs...) + N = div(length(x), M); plot!(x; label = "", kwargs...) + for ii in 1:M + Plots.plot!([ii*N, ii*N], [minimum(x), maximum(x)] ;color = :red, label = "", kwargs...) + end +end + +function plot_periodic_shooting(x, M; kwargs...) + Plots.plot() + plot_periodic_shooting!(x, M; kwargs...) +end +#################################################################################################### +function plot_eigenvals(br::AbstractResult, with_param = true; var = :param, k...) + p = getproperty(br.branch, var) + data = mapreduce(x -> x.eigenvals, hcat, br.eig) + if with_param + plot(p, real.(data'); k...) + else + plot(real.(data'); k...) + end +end \ No newline at end of file diff --git a/src/BifurcationKit.jl b/src/BifurcationKit.jl index 52894b7b..b7dfc01f 100644 --- a/src/BifurcationKit.jl +++ b/src/BifurcationKit.jl @@ -1,5 +1,5 @@ module BifurcationKit - using Printf, Dates, LinearMaps, BlockArrays, RecipesBase, StructArrays, Requires + using Printf, Dates, LinearMaps, BlockArrays, RecipesBase, StructArrays using Reexport @reexport using Setfield: @lens, @set, @set!, Lens import Setfield @@ -103,85 +103,12 @@ module BifurcationKit # plotting include("plotting/Utils.jl") + include("plotting/RecipesPlots.jl") # wrappers for SciML include("Diffeqwrap.jl") - using Requires - - function __init__() - # if Plots.jl is available, then we allow plotting of solutions - @require Plots="91a5bcdd-55d7-5caf-9e0b-520d859cae80" begin - using .Plots - include("plotting/RecipesPlots.jl") - get_plot_backend() = BK_Plots() - end - @require AbstractPlotting="537997a7-5e4e-5d89-9595-2241ea00577e" begin - using .AbstractPlotting: @recipe, layoutscene, Figure, Axis, lines! - include("plotting/RecipesMakie.jl") - end - - @require GLMakie="e9467ef8-e4e7-5192-8a1a-b1aee30e663a" begin - @info "Loading GLMakie code in BifurcationKit" - using .GLMakie: @recipe, Figure, Axis, lines!, PointBased, Point2f0, scatter! - include("plotting/RecipesMakie.jl") - get_plot_backend() = BK_Makie() - end - - @require JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819" begin - using .JLD2 - """ - Save solution / data in JLD2 file - - `filename` is for example "example.jld2" - - `sol` is the solution - - `p` is the parameter - - `i` is the index of the solution to be saved - """ - function save_to_file(iter::AbstractContinuationIterable, sol, p, i::Int64, br::ContResult) - if iter.contparams.save_to_file == false; return nothing; end - filename = iter.filename - # this allows to save two branches forward/backward in case - # bothside = true is passed to continuation - fd = iter.contparams.ds >=0 ? "fw" : "bw" - - # create a group in the JLD format - jldopen(filename*".jld2", "a+") do file - if haskey(file, "sol-$fd-$i") - delete!(file, "sol-$fd-$i") - end - mygroup = JLD2.Group(file, "sol-$fd-$i") - mygroup["sol"] = sol - mygroup["param"] = p - end - - jldopen(filename*"-branch.jld2", "a+") do file - if haskey(file, "branch"*fd) - delete!(file, "branch"*fd) - end - file["branch"*fd] = br - end - end - - # final save of branch, in case bothsided = true is used - function save_to_file(iter::AbstractContinuationIterable, br::ContResult) - if iter.contparams.save_to_file == false; return nothing; end - filename = iter.filename - - jldopen(filename*"-branch.jld2", "a+") do file - if haskey(file, "branchfw") - delete!(file, "branchfw") - end - if haskey(file, "branchbw") - delete!(file, "branchbw") - end - if haskey(file, "branch") - delete!(file, "branch") - end - file["branch"] = br - end - end - end - end + function save_to_file end # linear solvers export norminf diff --git a/src/DeflatedContinuation.jl b/src/DeflatedContinuation.jl index a39fb60b..77f333f3 100644 --- a/src/DeflatedContinuation.jl +++ b/src/DeflatedContinuation.jl @@ -72,6 +72,8 @@ end Base.lastindex(br::DCResult) = lastindex(br.branches) Base.getindex(br::DCResult, k::Int) = getindex(br.branches, k) Base.length(br::DCResult) = length(br.branches) +get_plot_vars(br::DCResult, vars) = get_plot_vars(br[1], vars) +_hasstability(br::DCResult) = _hasstability(br[1]) function Base.show(io::IO, brdc::DCResult; comment = "", prefix = " ") printstyled(io, "Deflated continuation result, # branches = $(length(brdc.branches))", "\n", color=:cyan, bold = true) @@ -164,12 +166,24 @@ function get_states_contResults(iter::DefContIterable, roots::Vector{Tvec}) wher end # plotting functions -function plot_DCont_branch(branches, nbrs::Int, nactive::Int, nstep::Int) - plot(branches..., label = "", title = "$nbrs branches, actives = $(nactive), step = $nstep") +function plot_DCont_branch(branches, + nbrs::Int, + nactive::Int, + nstep::Int) + _bcd_plot = get_plot_backend() + if _bcd_plot == BK_Plots() + plot(branches...; label = "", title = "$nbrs branches, actives = $(nactive), step = $nstep") + else + plot(branches...) + end for br in branches - length(br) > 1 && plot!([br.branch[end-1:end].param], [getproperty(br.branch,1)[end-1:end]], label = "", arrow = true, color = :red) + (length(br) > 1 && _bcd_plot == BK_Plots()) && plot!([br.branch[end-1:end].param], + [getproperty(br.branch,1)[end-1:end]], + label = "", arrow = true, color = :red) end - scatter!([br.branch[1].param for br in branches], [br.branch[1][1] for br in branches], marker = :cross, color=:green, label = "") |> display + _bcd_plot == BK_Plots() && scatter!([br.branch[1].param for br in branches], + [br.branch[1][1] for br in branches], + marker = :cross, color=:green, label = "") |> display end plotAllDCBranch(branches) = display(plot(branches..., label = "")) diff --git a/src/Utils.jl b/src/Utils.jl index 4b66ca80..f9da60af 100644 --- a/src/Utils.jl +++ b/src/Utils.jl @@ -6,6 +6,7 @@ abstract type AbstractPlotBackend end struct BK_NoPlot <: AbstractPlotBackend end struct BK_Plots <: AbstractPlotBackend end struct BK_Makie <: AbstractPlotBackend end + get_plot_backend() = BK_NoPlot() #################################################################################################### # functions for parameter handling diff --git a/src/plotting/RecipesPlots.jl b/src/plotting/RecipesPlots.jl index caddcbd3..fcc29bc9 100644 --- a/src/plotting/RecipesPlots.jl +++ b/src/plotting/RecipesPlots.jl @@ -169,78 +169,6 @@ RecipesBase.@recipe function Plots(brs::DCResult; end end -""" -Plot the branch of solutions during the continuation -""" -function plot_branch_cont(contres::ContResult, - state, - iter, - plotuserfunction) - sol = getsolution(state) - l = compute_eigenelements(iter) ? Plots.@layout([a{0.5w} [b; c]; e{0.2h}]) : Plots.@layout([a{0.5w} [b; c]]) - plot(layout = l ) - - plot!(contres ; filterspecialpoints = true, putspecialptlegend = false, - xlabel = get_lens_symbol(contres), - ylabel = _getfirstusertype(contres), - label = "", plotfold = false, subplot = 1) - - plotuserfunction(sol.u, sol.p; subplot = 3) - - # put arrow to indicate the order of computation - length(contres) > 1 && plot!([contres.branch[end-1:end].param], [getproperty(contres.branch,1)[end-1:end]], label = "", arrow = true, subplot = 1) - - if compute_eigenelements(iter) - eigvals = contres.eig[end].eigenvals - scatter!(real.(eigvals), imag.(eigvals), subplot=4, label = "", markerstrokewidth = 0, markersize = 3, color = :black, xlabel = "ℜ", ylabel = "ℑ") - # add stability boundary - maxIm = maximum(imag, eigvals) - minIm = minimum(imag, eigvals) - plot!([0, 0], [maxIm, minIm], subplot=4, label = "", color = :blue) - end - - plot!(contres; vars = (:step, :param), putspecialptlegend = false, plotspecialpoints = false, xlabel = "step [$(state.step)]", ylabel = get_lens_symbol(contres), label = "", subplot = 2) |> display - -end - -#################################################################################################### -# plotting function of the periodic orbits -function plot_periodic_potrap(x, M, Nx, Ny; ratio = 2, kwargs...) - @assert ratio > 0 "You need at least one component" - n = Nx*Ny - outpo = reshape(x[begin:end-1], ratio * n, M) - po = reshape(x[1:n,1], Nx, Ny) - rg = 2:6:M - for ii in rg - po = hcat(po, reshape(outpo[1:n,ii], Nx, Ny)) - end - heatmap!(po; color = :viridis, fill=true, xlabel = "space", ylabel = "space", kwargs...) - for ii in eachindex(rg) - plot!([ii*Ny, ii*Ny], [1, Nx]; color = :red, width = 3, label = "", kwargs...) - end -end - -function plot_periodic_potrap(outpof, n, M; ratio = 2) - @assert ratio > 0 "You need at least one component" - outpo = reshape(outpof[begin:end-1], ratio * n, M) - if ratio == 1 - heatmap(outpo[1:n,:]', ylabel="Time", color=:viridis) - else - plot(heatmap(outpo[1:n,:]', ylabel="Time", color=:viridis), - heatmap(outpo[n+2:end,:]', color=:viridis)) - end -end -#################################################################################################### -function plot_periodic_shooting!(x, M; kwargs...) - N = div(length(x), M); plot!(x; label = "", kwargs...) - for ii in 1:M - plot!([ii*N, ii*N], [minimum(x), maximum(x)] ;color = :red, label = "", kwargs...) - end -end - -function plot_periodic_shooting(x, M; kwargs...) - plot();plot_periodic_shooting!(x, M; kwargs...) -end #################################################################################################### RecipesBase.@recipe function Plots(sol::SolPeriodicOrbit; indx = nothing @@ -289,16 +217,6 @@ RecipesBase.@recipe function f(bd::Nothing) nothing end #################################################################################################### -function plot_eigenvals(br::AbstractResult, with_param = true; var = :param, k...) - p = getproperty(br.branch, var) - data = mapreduce(x -> x.eigenvals, hcat, br.eig) - if with_param - plot(p, real.(data'); k...) - else - plot(real.(data'); k...) - end -end -#################################################################################################### # plot recipe for codim 2 plot # TODO Use dispatch for this RecipesBase.@recipe function Plots(contres::AbstractResult{Tk, Tprob}; diff --git a/src/plotting/Utils.jl b/src/plotting/Utils.jl index c6e0bffd..911135fe 100644 --- a/src/plotting/Utils.jl +++ b/src/plotting/Utils.jl @@ -82,4 +82,12 @@ function filter_bifurcations(bifpt) 0 < ii <= length(bifpt) && push!(res, (type = bifpt[ii].type, idx = bifpt[ii].idx, param = bifpt[ii].param, printsol = bifpt[ii].printsol, status = bifpt[ii].status) ) return res[2:end] -end \ No newline at end of file +end + +function plot end +function plot! end +function plot_branch_cont end +function plot_periodic_potrap end +function plot_periodic_shooting! end +function plot_periodic_shooting end +function plot_eigenvals end \ No newline at end of file