Proof of concept Recipes for Plots

Project.toml

@@ -11,13 +11,15 @@ Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca"
 Requires = "ae029012-a4dd-5104-9daa-d747884805df"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
 
 [compat]
 Arpack = "0.3, 0.4, 0.5"
 PlotUtils = "1.0"
 Colors = "0.12"
 ProgressMeter = "1"
 Requires = "1.0"
+RecipesBase = "1.1.1"
 julia = "1"
 
 [extras]

src/Pseudospectra.jl

@@ -116,6 +116,8 @@ include("transients.jl")
 include("plotter.jl")
 include("zooming.jl")
 
+include("plots/PSAPlots.jl")
+
 """
     new_matrix(A::AbstractMatrix, opts::Dict{Symbol,Any}=()) -> ps_data
 

src/plots/PSAPlots.jl

@@ -0,0 +1,26 @@
+module PSAPlots
+    using RecipesBase
+    using Pseudospectra
+    using LinearAlgebra
+    using PlotUtils
+    using Printf
+
+    import Pseudospectra: vec2ax, psa_compute
+    import Pseudospectra: psmode_inv_lanczos, oneeigcond
+
+    include("utils.jl")
+    include("defaults.jl")
+    #include("colormap.jl")
+
+    include("recipes/spectralportrait.jl")
+    include("recipes/plotmode.jl")
+    include("recipes/mtxpowersplot.jl")
+    include("recipes/mtxexpsplot.jl")
+    include("recipes/surfplot.jl")
+
+    export spectralportrait, spectralportrait!
+    export plotmode, plotmode!
+    export mtxpowersplot, mtxpowersplot!
+    export mtxexpsplot, mtxexpsplot!
+    export surfplot, surfplot!
+end

src/plots/colormap.jl

@@ -0,0 +1,88 @@
+
+"""
+    The colormap from EigTool. This provides clearly
+    distinguishable colors for a line-contour plot.
+"""
+function etcgrad()
+    cm = [
+        0   1.00000000000000   0.10000000000000
+        0   0.92500000000000   0.17500000000000
+        0   0.85000000000000   0.25000000000000
+        0   0.77500000000000   0.32500000000000
+        0   0.70000000000000   0.40000000000000
+        0   0.62500000000000   0.47500000000000
+        0   0.55000000000000   0.55000000000000
+        0   0.47500000000000   0.62500000000000
+        0   0.40000000000000   0.70000000000000
+        0   0.32500000000000   0.77500000000000
+        0   0.25000000000000   0.85000000000000
+        0   0.17500000000000   0.92500000000000
+        0   0.10000000000000   1.00000000000000
+        0.06153846153846   0.09230769230769   1.00000000000000
+        0.12307692307692   0.08461538461538   1.00000000000000
+        0.18461538461538   0.07692307692308   1.00000000000000
+        0.24615384615385   0.06923076923077   1.00000000000000
+        0.30769230769231   0.06153846153846   1.00000000000000
+        0.36923076923077   0.05384615384615   1.00000000000000
+        0.43076923076923   0.04615384615385   1.00000000000000
+        0.49230769230769   0.03846153846154   1.00000000000000
+        0.55384615384615   0.03076923076923   1.00000000000000
+        0.61538461538462   0.02307692307692   1.00000000000000
+        0.67692307692308   0.01538461538462   1.00000000000000
+        0.73846153846154   0.00769230769231   1.00000000000000
+        0.80000000000000                  0   1.00000000000000
+        0.80769230769231                  0   0.92307692307692
+        0.81538461538462                  0   0.84615384615385
+        0.82307692307692                  0   0.76923076923077
+        0.83076923076923                  0   0.69230769230769
+        0.83846153846154                  0   0.61538461538462
+        0.84615384615385                  0   0.53846153846154
+        0.85384615384615                  0   0.46153846153846
+        0.86153846153846                  0   0.38461538461538
+        0.86923076923077                  0   0.30769230769231
+        0.87692307692308                  0   0.23076923076923
+        0.88461538461538                  0   0.15384615384615
+        0.89230769230769                  0   0.07692307692308
+        0.90000000000000                  0                  0
+        0.86923076923077                  0                  0
+        0.83846153846154                  0                  0
+        0.80769230769231                  0                  0
+        0.77692307692308                  0                  0
+        0.74615384615385                  0                  0
+        0.71538461538462                  0                  0
+        0.68461538461538                  0                  0
+        0.65384615384615                  0                  0
+        0.62307692307692                  0                  0
+        0.59230769230769                  0                  0
+        0.56153846153846                  0                  0
+        0.53076923076923                  0                  0
+        0.50000000000000                  0                  0
+        0.54166666666667   0.05000000000000                  0
+        0.58333333333333   0.10000000000000                  0
+        0.62500000000000   0.15000000000000                  0
+        0.66666666666667   0.20000000000000                  0
+        0.70833333333333   0.25000000000000                  0
+        0.75000000000000   0.30000000000000                  0
+        0.79166666666667   0.35000000000000                  0
+        0.83333333333333   0.40000000000000                  0
+        0.87500000000000   0.45000000000000                  0
+        0.91666666666667   0.50000000000000                  0
+        0.95833333333333   0.55000000000000                  0
+        1.00000000000000   0.60000000000000                  0
+    ]
+    nc = size(cm,1)
+    cc = Vector{RGBA{Float64}}()
+    vv = zeros(nc)
+    for i in 1:nc
+        col = RGBA(RGB(cm[i,:]...))
+        push!(cc,col)
+        vv[i]=(i-1)/(nc-1)
+    end
+    if :ContinuousColorGradient in names(Plots.PlotUtils, all=true)
+        return PlotUtils.ContinuousColorGradient(cc,vv)
+    else
+        return PlotUtils.ColorGradient(cc,vv)
+    end
+end
+
+const eigtool_cgrad = etcgrad()

src/plots/defaults.jl

@@ -0,0 +1,8 @@
+
+_basic_psa_opts(zoom,ps_dict) = Dict{Symbol,Any}(
+    #:levels=>expandlevels(zoom.levels),
+    :recompute_levels=>zoom.autolev,
+    :proj_lev=>zoom.proj_lev,
+    :scale_equal=>zoom.scale_equal,
+    :real_matrix=>ps_dict[:Aisreal],
+    :verbosity=>0)

src/plots/recipes/mtxexpsplot.jl

@@ -0,0 +1,40 @@
+@userplot MtxExpsPlot
+
+
+RecipesBase.@recipe function f(p::MtxExpsPlot; dt=0.1, nmax=50, lbmethod=:none)
+    if length(p.args) != 1 || !(typeof(p.args[1]) <: PSAStruct)
+        error("mtxexpsplot must be given a PSAStruct. Got $(typeof(p.args))")
+    end
+    ps_data::PSAStruct = p.args[1] # NB hint the type to the compiler
+    A = ps_data.input_matrix
+
+    ts, transient = get_mtxexpnorm(A, dt, nmax)
+    alp = maximum(real.(ps_data.ps_dict[:ews]))
+
+    #layout := @layout [a; b]
+    #link := :x
+    seriestype := :path
+    label := ""
+
+    @series begin
+        #subplot := 1
+        label := "Matrix norm"
+        ts, transient
+    end
+    @series begin
+        #subplot := 1
+        label := "Spectral Radius"
+        ts, exp.(alp*ts)
+    end
+    #=
+    @series begin
+        subplot := 2
+        yscale := :log10
+        ts, transient
+    end
+    @series begin
+        subplot := 2
+        ts, exp.(alp*ts)
+    end
+    =#
+end

src/plots/recipes/mtxpowersplot.jl

@@ -0,0 +1,75 @@
+@userplot MtxPowersPlot
+
+"""
+    mtxpowersplot(A, nmax)
+
+Plot the matrix norm of `A^k` for `k` up to `nmax`
+
+# Examples
+```jldoctest
+julia> opts = Dict{Symbol,Any}(:npts => 40,:ax => [-110,110,-110,110])
+julia> A = diagm(1 => 100.0*ones(4)) + diagm(-2 => 1e-4*ones(3))
+5×5 Matrix{Float64}:
+ 0.0     100.0       0.0       0.0    0.0
+ 0.0       0.0     100.0       0.0    0.0
+ 0.0001    0.0       0.0     100.0    0.0
+ 0.0       0.0001    0.0       0.0  100.0
+ 0.0       0.0       0.0001    0.0    0.0
+ julia> ps_data = new_matrix(A,opts)
+ [...]
+ julia> driver!(ps_data, opts)
+ julia> mtxpowersplot(ps_data, nmax=20)
+```
+"""
+function mtxpowersplot end
+
+RecipesBase.@recipe function f(p::MtxPowersPlot;
+    nmax::Integer = 50,
+    lbmethod = :none,
+    lbdk = 0.25
+    )
+    if length(p.args) != 1 || !(typeof(p.args[1]) <: PSAStruct)
+        error("mtxpowersplot must be given a PSAStruct. Got $(typeof(p.args))")
+    end
+    ps_data::PSAStruct = p.args[1]
+    A = ps_data.input_matrix
+
+    powers, transient = get_mtxpowersnorm(A, nmax)
+
+    alp = maximum(abs.(ps_data.ps_dict[:ews]))
+
+    # NB. We need Plots to use @layout
+    #layout := @layout [a; b]
+    #link := :x # link the x-axis of both subplots
+
+    # normal plot
+    @series begin
+        seriestype := :path
+        subplot := 1
+        label := "Matrix norm"
+        powers, transient
+    end
+    @series begin
+        seriestype := :path
+        subplot := 1
+        label := "Spectral Radius"
+        powers, alp.^powers
+    end
+    # log plot
+    #=
+    @series begin
+        seriestype := :path
+        subplot := 2
+        label := ""
+        yscale := :log10
+        powers, transient
+    end
+    @series begin
+        seriestype := :path
+        subplot := 2
+        label := ""
+        yscale := :log10
+        powers, alp.^powers
+    end
+    =#
+end

src/plots/recipes/plotmode.jl

@@ -0,0 +1,71 @@
+@userplot PlotMode
+
+function get_mode_title(pseudo, approx, z, s)
+    prefix = pseudo ? "Pseudomode: " : "Eigenmode: "
+    rel = approx ? "\\approx" : "="
+    formatted = @sprintf("%15.2e", s)
+    infix = pseudo ? ("\$\\epsilon = " * formatted * "\$") : ("\$\\kappa(\\lambda) " * rel * formatted * "\$")
+    λ_str = @sprintf("%12g%+12gi", real(z), imag(z))
+    return prefix * infix * "\n  \$\\lambda = " * λ_str * "\$"
+end
+
+"""
+    plotmode(A, z, U=I; ...)
+
+Plot the (pseudo)eigenmode of `A` associated with (pseudo)eigenvalue `z`.
+
+# Arguments
+* `A` an arbitrary square matrix
+* `z` the real or complex (pseudo)eigenvalue whose eigenmode is to be plotted.
+* `U` optional Schur factor. The eigenmode is premultiplied by `U` before plotting.
+# Keyword arguments
+* `pseudo::Bool = false` Whether to plot the eigenmode or pseudomode.
+* `approx::Bool = false`
+* `verbosity::Int = 0` Verbosity parameter passed when computing the eigenmode.
+
+# Examples
+```jldoctest
+julia> A = Pseudospectra.grcar(40)
+julia> plotmode(A, 0.5+2.0im)
+```
+"""
+function plotmode end
+
+RecipesBase.@recipe function f(p::PlotMode, U = I; pseudo=false, approx=false, verbosity=0)
+    if length(p.args) != 2
+        error("plotmode must be given 2 arguments")
+    end
+    A = p.args[1]
+    z = p.args[2]
+
+    SS, q = get_mode(A, z, pseudo, verbosity)
+
+    q = U * q
+
+    #showlog = sum(abs.(diff(q))) * (1/length(q)) >= 10*eps()
+
+    # NB. The title is broken
+    #title_str = get_mode_title(pseudo, approx, z, SS[end])
+    #title := title_str
+
+    @series begin
+        seriestype := :path
+        linecolor := (pseudo ? :magenta : :cyan)
+        label := "Real part"
+        real.(q)
+    end
+
+    seriestype := :path
+    linecolor := :black
+    linestyle := :dash
+
+    @series begin
+        label := "Modulus"
+        abs.(q)
+    end
+
+    @series begin
+        label := ""
+        -abs.(q)
+    end
+end