text:
import Colors
export lsimplot, stepplot, impulseplot, bodeplot, nyquistplot, sigmaplot, marginplot, setPlotScale, gangoffour, gangoffourplot, gangofseven, pzmap, pzmap!, nicholsplot

_PlotScale = "log10"
_PlotScaleFunc = :log10
_PlotScaleStr = ""

"""`setPlotScale(str)`

Set the default scale of magnitude in `bodeplot` and `sigmaplot`.
`str` should be either `"dB"` or `"log10"`."""
function setPlotScale(str::AbstractString)
    if str == "dB"
        plotSettings = (str, :identity, "(dB)")
    elseif str == "log10"
        plotSettings = (str, :log10, "")
    else
        error("Scale must be set to either \"dB\" or \"log10\"")
    end
    global _PlotScale, _PlotScaleFunc, _PlotScaleStr
    _PlotScale, _PlotScaleFunc, _PlotScaleStr = plotSettings
end

"""
Get atributes from xlims or ylims
default to extrema(wmag) if xlims/ylims not defined or empty
"""
function getlims(xylims, plotattributes, wmag)
    lims = get(plotattributes, xylims, extrema(wmag))
    if !isa(lims, Tuple{<:Number, <:Number}) # If x/ylims not supplied as empty
        lims = extrema(wmag)
    end
    if !isempty(get_serieslist(plotattributes))
        subplot = get(plotattributes, :subplot, 0)
        subplot == 0 && (return lims)
        se = seriesextrema(xylims, plotattributes, subplot)
        lims = extremareducer(lims, se)
    end
    lims
end

get_serieslist(plotattributes) = plotattributes[:plot_object].series_list
get_serieslist(plotattributes, subplot) = plotattributes[:plot_object].subplots[subplot].series_list

function seriesextrema(xylims, plotattributes, subplot)
    serieslist = get_serieslist(plotattributes, subplot)
    isempty(serieslist) && (return (Inf, -Inf))
    sym = xylims === :xlims ? :x : :y
    mapreduce(extremareducer, serieslist) do series
        extrema(series[sym])
    end
end
extremareducer(x,y) = (min(x[1],y[1]),max(x[2],y[2]))

function getLogTicks(x, minmax)
    minx, maxx =  minmax
    major_minor_limit = 6
    minor_text_limit  = 8
    min               = minx <= 0 ? minimum(x) : ceil(log10(minx))
    max               = floor(log10(maxx))
    major             = exp10.(min:max)
    if Plots.backend() ∉ [Plots.GRBackend(), Plots.PlotlyBackend()]
        majorText = [latexstring("\$10^{$(round(Int64,i))}\$") for i = min:max]
    else
        majorText = ["10^{$(round(Int64,i))}" for i = min:max]
    end
    if max - min < major_minor_limit
        minor     = [j*exp10(i) for i = (min-1):(max+1) for j = 2:9]
        if Plots.backend() ∉ [Plots.GRBackend(), Plots.PlotlyBackend()]
            minorText = [latexstring("\$$j\\cdot10^{$(round(Int64,i))}\$") for i = (min-1):(max+1) for j = 2:9]
        else
            minorText = ["$j*10^{$(round(Int64,i))}" for i = (min-1):(max+1) for j = 2:9]
        end

        ind       = findall(minx .<= minor .<= maxx)
        minor     = minor[ind]
        minorText = minorText[ind]
        if length(minor) > minor_text_limit
            minorText = [" " for t in minorText]#fill!(minorText, L" ")
        end
        perm = sortperm([major; minor])
        return [major; minor][perm], [majorText; minorText][perm]

    else
        return major, majorText
    end
end

@userplot Lsimplot

"""`fig = lsimplot(sys::LTISystem, u, t; x0=0, method)`

`lsimplot(LTISystem[sys1, sys2...], u, t; x0, method)`

Calculate the time response of the `LTISystem`(s) to input `u`. If `x0` is
not specified, a zero vector is used.

Continuous time systems are discretized before simulation. By default, the
method is chosen based on the smoothness of the input signal. Optionally, the
`method` parameter can be specified as either `:zoh` or `:foh`.
"""
lsimplot

@recipe function lsimplot(p::Lsimplot; method=nothing)
    if length(p.args) < 3
        error("Wrong number of arguments")
    end
    systems,u,t = p.args[1:3]

    if !isa(systems,AbstractArray)
        systems = [systems]
    end
    if method == nothing
        method = _issmooth(u) ? :foh : :zoh
    end
    if !_same_io_dims(systems...)
        error("All systems must have the same input/output dimensions")
    end
    ny, nu = size(systems[1])
    layout --> (ny,1)
    s2i(i,j) = LinearIndices((ny,1))[j,i]
    for (si,s) in enumerate(systems)
        s = systems[si]
        y = length(p.args) >= 4 ? lsim(s, u, t, x0=p.args[4], method=method)[1] : lsim(s, u, t, method=method)[1]
        seriestype := iscontinuous(s) ? :path : :steppost
        for i=1:ny
            ytext = (ny > 1) ? "Amplitude to: y($i)" : "Amplitude"
            @series begin
                xguide  --> "Time (s)"
                yguide  --> ytext
                title   --> "System Response"
                subplot --> s2i(1,i)
                label     --> "\$G_{$(si)}\$"
                t,  y[:, i]
            end
        end
    end
end

@userplot Stepplot
@userplot Impulseplot
"""
`stepeplot(sys[, Tf]; kwargs...)`` or `stepplot(sys[, t]; kwargs...)``
Plot step response of  `sys` until final time `Tf` or at time points in the vector `t`.
If not defined, suitable values are chosen based on `sys`.
See also [`step`](@ref)

`kwargs` is sent as argument to Plots.plot.
"""
stepplot

"""
    `impulseplot(sys[, Tf]; kwargs...)`` or `impulseplot(sys[, t]; kwargs...)``
Plot impulse response of `sys` until final time `Tf` or at time points in the vector `t`.
If not defined, suitable values are chosen based on `sys`.
See also [`impulse`](@ref)

`kwargs` is sent as argument to Plots.plot.
"""
impulseplot

for (func, title, typ) = ((step, "Step Response", Stepplot), (impulse, "Impulse Response", Impulseplot))
    funcname = Symbol(func,"plot")

    @recipe function f(p::typ)
        systems = p.args[1]
        if !isa(systems, AbstractArray)
            systems = [systems]
        end
        if !_same_io_dims(systems...)
            error("All systems must have the same input/output dimensions")
        end
        ny, nu = size(systems[1])
        layout --> (ny,nu)
        titles = fill("", 1, ny*nu)
        title --> titles
        s2i(i,j) = LinearIndices((ny,nu))[i,j]
        for (si,s) in enumerate(systems)
            y,t = func(s, p.args[2:end]...)
            for i=1:ny
                for j=1:nu
                    ydata = reshape(y[:, i, j], size(t, 1))
                    style = iscontinuous(s) ? :path : :steppost
                    ttext = (nu > 1 && i==1) ? title*" from: u($j) " : title
                    titles[s2i(i,j)] = ttext
                    ytext = (ny > 1 && j==1) ? "Amplitude to: y($i)" : "Amplitude"
                    @series begin
                        seriestype := style
                        xguide --> "Time (s)"
                        yguide --> ytext
                        subplot --> s2i(i,j)
                        label --> "\$G_{$(si)}\$"
                        t, ydata
                    end
                end
            end
        end
    end

end

"""
    _processfreqplot(plottype, system::LTISystem, [w])
    _processfreqplot(plottype, system::AbstractVector{<:LTISystem}, [w])

    Calculate default frequency vector and put system in array of not already array.
    `plottype` is one of `Val{:bode}, Val{:nyquist}, ...`
    for which `_default_freq_vector` is defined.
    Check that system dimensions are compatible.
"""
_processfreqplot(plottype, system::LTISystem, args...) =
    _processfreqplot(plottype, [system], args...)
# Catch when system is not vector, with and without frequency input

# Cantch correct form
function _processfreqplot(plottype, systems::AbstractVector{<:LTISystem},
            w = _default_freq_vector(systems, plottype))

    if !_same_io_dims(systems...)
        error("All systems must have the same input/output dimensions")
    end
    return systems, w
end

@userplot Bodeplot
## FREQUENCY PLOTS ##
"""`fig = bodeplot(sys, args...)`, `bodeplot(LTISystem[sys1, sys2...], args...; plotphase=true, kwargs...)`

Create a Bode plot of the `LTISystem`(s). A frequency vector `w` can be
optionally provided. To change the Magnitude scale see `setPlotScale(str)`

`kwargs` is sent as argument to Plots.plot."""
bodeplot

@recipe function bodeplot(p::Bodeplot; plotphase=true, ylimsphase=(), unwrap=true)
    systems, w = _processfreqplot(Val{:bode}(), p.args...)
    ny, nu = size(systems[1])
    s2i(i,j) = LinearIndices((nu,(plotphase ? 2 : 1)*ny))[j,i]
    layout --> ((plotphase ? 2 : 1)*ny,nu)
    nw = length(w)
    xticks --> getLogTicks(w, getlims(:xlims, plotattributes, w))

    for (si,s) = enumerate(systems)
        mag, phase = bode(s, w)[1:2]
        if _PlotScale == "dB" # Set by setPlotScale(str) globally
            mag = 20*log10.(mag)
        end

        xlab = plotphase ? "" : "Frequency (rad/s)"
        group_ind = 0
        for j=1:nu
            for i=1:ny
                group_ind += 1
                magdata = vec(mag[:, i, j])
                if all(magdata .== -Inf)
                    # 0 system, don't plot anything
                    continue
                end
                phasedata = vec(phase[:, i, j])
                @series begin
                    grid      --> true
                    yscale    --> _PlotScaleFunc
                    xscale    --> :log10
                    if _PlotScale != "dB"
                        yticks    --> getLogTicks(magdata, getlims(:ylims, plotattributes, magdata))
                    end
                    xguide    --> xlab
                    yguide    --> "Magnitude $_PlotScaleStr"
                    subplot   --> s2i((plotphase ? (2i-1) : i),j)
                    title     --> "Bode plot from: u($j)"
                    label     --> "\$G_{$(si)}\$"
                    group     --> group_ind
                    w, magdata
                end
                plotphase || continue

                @series begin
                    grid      --> true
                    xscale    --> :log10
                    ylims      := ylimsphase
                    yguide    --> "Phase (deg)"
                    subplot   --> s2i(2i,j)
                    xguide    --> "Frequency (rad/s)"
                    label     --> "\$G_{$(si)}\$"
                    group     --> group_ind
                    w, unwrap ? ControlSystems.unwrap(phasedata.*(pi/180)).*(180/pi) : phasedata
                end

            end
        end
    end
end

@recipe function f(::Type{Val{:bodemag}}, x, y, z)
    w = x
    magdata = y
    seriestype := :path
    primary := false
    @series begin
        grid   --> true
        yscale --> :log10
        xscale --> :log10
        yguide --> "Magnitude"
        xticks --> getLogTicks(w,  getlims(:xlims, plotattributes, w))
        yticks --> getLogTicks(magdata,  getlims(:ylims, plotattributes, magdata))
        x := w; y := magdata
        ()
    end
    x := []
    y := []
    ()
end
@recipe function f(::Type{Val{:bodephase}}, x, y, z)
    w = x
    phasedata = y
    seriestype := :path
    primary := false
    @series begin
        grid   --> true
        xscale --> :log10
        yguide --> "Phase (deg)"
        xguide --> "Frequency (rad/s)"
        xticks --> getLogTicks(w, getlims(:xlims, plotattributes, w))
        x := w; y := phasedata
        ()
    end
    x := []
    y := []
    ()
end

@userplot Nyquistplot
"""`fig = nyquistplot(sys; gaincircles=true, kwargs...)`, `nyquistplot(LTISystem[sys1, sys2...]; gaincircles=true, kwargs...)`

Create a Nyquist plot of the `LTISystem`(s). A frequency vector `w` can be
optionally provided.

`gaincircles` plots the circles corresponding to |S(iω)| = 1 and |T(iω)| = 1, where S and T are
the sensitivity and complementary sensitivity functions.

`kwargs` is sent as argument to plot."""
nyquistplot
@recipe function nyquistplot(p::Nyquistplot; gaincircles=true)
    systems, w = _processfreqplot(Val{:nyquist}(), p.args...)
    ny, nu = size(systems[1])
    nw = length(w)
    layout --> (ny,nu)
    s2i(i,j) = LinearIndices((ny,nu))[j,i]
    # Ensure that `axes` is always a matrix of handles
    for (si,s) = enumerate(systems)
        re_resp, im_resp = nyquist(s, w)[1:2]
        for j=1:nu
            for i=1:ny
                redata      = re_resp[:, i, j]
                imdata      = im_resp[:, i, j]
                @series begin
                    ylims   --> (min(max(-20,minimum(imdata)),-1), max(min(20,maximum(imdata)),1))
                    xlims   --> (min(max(-20,minimum(redata)),-1), max(min(20,maximum(redata)),1))
                    title --> "Nyquist plot from: u($j)"
                    yguide --> "To: y($i)"
                    subplot --> s2i(i,j)
                    label --> "\$G_{$(si)}\$"
                    hover --> [Printf.@sprintf("ω = %.3f", w) for w in w]
                    (redata, imdata)
                end
                # Plot rings
                if gaincircles && si == length(systems)
                    v = range(0,stop=2π,length=100)
                    S,C = sin.(v),cos.(v)
                    @series begin
                        primary := false
                        linestyle := :dash
                        linecolor := :black
                        seriestype := :path
                        markershape := :none
                        (C,S)
                    end
                    @series begin
                        primary := false
                        linestyle := :dash
                        linecolor := :black
                        seriestype := :path
                        markershape := :none
                        (C .-1,S)
                    end
                end

            end
        end
    end
end

@userplot Nicholsplot

"""
fig = `nicholsplot{T<:LTISystem}(systems::Vector{T}, w::AbstractVector; kwargs...)`

Create a Nichols plot of the `LTISystem`(s). A frequency vector `w` can be
optionally provided.

Keyword arguments:
```
text = true
Gains = [12, 6, 3, 1, 0.5, -0.5, -1, -3, -6, -10, -20, -40, -60]
pInc = 30
sat = 0.4
val = 0.85
fontsize = 10
```

`pInc` determines the increment in degrees between phase lines.

`sat` ∈ [0,1] determines the saturation of the gain lines

`val` ∈ [0,1] determines the brightness of the gain lines

Additional keyword arguments are sent to the function plotting the systems and can be
used to specify colors, line styles etc. using regular Plots.jl syntax

This function is based on code subject to the two-clause BSD licence
Copyright 2011 Will Robertson
Copyright 2011 Philipp Allgeuer

"""
nicholsplot
@recipe function nicholsplot(p::Nicholsplot;
    text     = true,
    Gains    = [12, 6, 3, 1, 0.5, -0.5, -1, -3, -6, -10, -20, -40, -60],
    pInc     = 30,
    sat      = 0.4,
    val      = 0.85,
    fontsize = 10)

    systems, w = _processfreqplot(Val{:nyquist}(), p.args...)
    ny, nu = size(systems[1])

    if isdiscrete(systems[1])
        w_nyquist = 2π/systems[1].Ts
        w = w[w.<= w_nyquist]
    end
    nw = length(w)

    # Gain circle functions
    angle(x)        = unwrap(atan.(imag.(x),real.(x)))
    RadM(m)         = @. abs(m/(m^2-1))
    CentreM(m)      = @. m^2/(1-m^2)
    Ny(mdb,t)       = @. CentreM(10^(mdb/20))+RadM(10^(mdb/20)).*(cosd(t)+im.*sind(t))
    Niϕ(mdb,t)      = @. rad2deg((angle(Ny(mdb,t))))
    Ni_Ga(mdb,t)    = @. 20 * log10(abs(Ny(mdb,t)))

    # Phase circle functions
    Radϕ(ϕ)         = @. 1 / (2 * abs(sind(ϕ)))
    Nyℜ(ϕ,t)        = @. -0.5+Radϕ(ϕ)*cosd(t+mod(ϕ,180)-90)
    Nyℑ(ϕ,t)        = @. 1 / (2 .* tand(ϕ))+Radϕ(ϕ)*sind(t+mod(ϕ,180)-90)
    Niϕϕ(ϕ,t)       = @. rad2deg((angle(Nyℜ(ϕ,t)+im*Nyℑ(ϕ,t))))+360*(round(ϕ/360,RoundToZero)+(0t<0))
    Ni_Gaϕ(ϕ,t)     = @. 20 * log10(abs(Nyℜ(ϕ,t)+im*Nyℑ(ϕ,t)))
    Ni_La(ϕ)        = @. 0.090*10^(ϕ/60)
    getColor(mdb)   = convert(Colors.RGB,Colors.HSV(360*((mdb-minimum(Gains))/(maximum(Gains)-minimum(Gains)))^1.5,sat,val))

    megaangles      = vcat(map(s -> 180/π*angle(vec(freqresp(s, w))), systems)...)
    filter!(x-> !isnan(x), megaangles)
    extremeangles = extrema(megaangles)
    extremeangles = floor(extremeangles[1]/180)*180, ceil(extremeangles[2]/180)*180
    PCyc            = Set{Int}(floor.(Int,megaangles/360)) |> collect |> sort
    # PCyc            = extremeangles[1]:pInc:extremeangles[2]

    # yticks := (Float64[],String[])
    yguide --> "Open-loop gain [dB]"
    xguide --> "Open-loop phase [deg]"

    #  Gain circles
    for k=Gains
        ϕVals   =Niϕ(k,-180:1:180)
        GVals   =Ni_Ga(k,-180:1:180)
        for l in PCyc
            @series begin
                linewidth := 1
                linecolor := getColor(k)
                grid --> false
                if text
                    offset  = (l+1)
                    TextX   = Niϕ(k,210) .+offset
                    TextY   = Ni_Ga(k,210)
                    annotations := (TextX,TextY,Plots.text("$(string(k)) dB",fontsize))
                end
                ϕVals .+ 360(l+1),GVals
            end
        end
    end

    #  Phase circles
    PCycbottom = (PCyc[1] < 0 ? PCyc[1]*360 : (PCyc[1]-1)*360)
    PCyctop = (PCyc[end] < 0 ? (PCyc[end]+1)*360 : (PCyc[end])*360)

    Phi=(PCycbottom):pInc:(PCyctop)
    T1 = 10.0 .^range(-4,stop=log10(180), length=300)
    T2 = [T1; 360 .- reverse(T1,dims=1)]

    for k=(Phi .+ 180)
        if abs(sind(k))<1e-3
            @series begin
                linewidth := 1
                linecolor := Colors.RGB(0.75*[1, 1, 1]...)
                [k,k],[-110,25]
            end
            if cosd(5)>0
                TextX=k
                TextY=1
            else
                TextX=k
                TextY=-46.5
            end
        else
            @series begin
                linewidth := 1
                linecolor := Colors.RGB(0.75*[1,1,1]...)
                Niϕϕ(k,T2),Ni_Gaϕ(k,T2)
            end
            Offset=k-180*floor(Int,k/180);
            if sign(sind(k))==1
                TextX=Niϕϕ(k,Ni_La(180-Offset))
                TextY=Ni_Gaϕ(k,Ni_La(180-Offset))
            else
                TextX=Niϕϕ(k,-Ni_La(Offset))+360
                TextY=Ni_Gaϕ(k,-Ni_La(Offset))
            end
        end
        TextX
        annotations := (TextX,TextY,Plots.text("$(string(k))°",fontsize))

        title --> "Nichols chart"
        grid --> false
        legend --> false
        xguide --> "Phase [deg]"
        yguide --> "Magnitude [dB]"

    end

    extremas = extrema(Gains)
    # colors = [:blue, :cyan, :green, :yellow, :orange, :red, :magenta]
    for (sysi,s) = enumerate(systems)
        ℜresp, ℑresp        = nyquist(s, w)[1:2]
        ℜdata               = dropdims(ℜresp, dims=(2,3))
        ℑdata               = dropdims(ℑresp, dims=(2,3))
        mag                 = 20*log10.(sqrt.(ℜdata.^2 + ℑdata.^2))
        angles              = 180/π*angle(im*ℑdata.+ℜdata)
        extremas = extrema([extremas..., extrema(mag)...])
        @series begin
            linewidth --> 2
            hover --> [Printf.@sprintf("ω = %.3f", w) for w in w]
            angles, mag
        end
    end
    ylims --> extremas
    xlims --> extrema(PCyc*360)
    nothing

end

@userplot Sigmaplot
"""`sigmaplot(sys, args...)`, `sigmaplot(LTISystem[sys1, sys2...], args...)`

Plot the singular values of the frequency response of the `LTISystem`(s). A
frequency vector `w` can be optionally provided.

`kwargs` is sent as argument to Plots.plot."""
sigmaplot
@recipe function sigmaplot(p::Sigmaplot)
    systems, w = _processfreqplot(Val{:sigma}(), p.args...)
    ny, nu = size(systems[1])
    nw = length(w)
    title --> "Sigma Plot"
    xguide --> "Frequency (rad/s)",
    yguide --> "Singular Values $_PlotScaleStr"
    for (si, s) in enumerate(systems)
        sv = sigma(s, w)[1]
        if _PlotScale == "dB"
            sv = 20*log10.(sv)
        end
        for i in 1:size(sv, 2)
            @series begin
                xscale --> :log10
                yscale --> _PlotScaleFunc
                seriescolor --> si
                w, sv[:, i]
            end
        end
    end
end

"""`fig = marginplot(sys::LTISystem [,w::AbstractVector];  kwargs...)`, `marginplot(sys::Vector{LTISystem}, w::AbstractVector;  kwargs...)`

Plot all the amplitude and phase margins of the system(s) `sys`.
A frequency vector `w` can be optionally provided.

`kwargs` is sent as argument to Plots.plot."""
function marginplot(systems::Union{AbstractVector{T},T}, args...; kwargs...) where T<:LTISystem
    systems, w = _processfreqplot(Val{:bode}(), systems, args...)
    ny, nu = size(systems[1])
    fig = bodeplot(systems, w; kwargs...)
    s2i(i,j) = LinearIndices((ny,2nu))[j,i]
    titles = Array{AbstractString}(undef, nu,ny,2,2)
    titles[:,:,1,1] .= "Gm: "
    titles[:,:,2,1] .= "Pm: "
    titles[:,:,1,2] .= "Wgm: "
    titles[:,:,2,2] .= "Wpm: "
    for (si, s) in enumerate(systems)
        for j=1:nu
            for i=1:ny
                wgm, gm, wpm, pm, fullPhase = sisomargin(s[i,j],w, full=true, allMargins=true)
                # Let's be reasonable, only plot 5 smallest gain margins
                if length(gm) > 5
                    @warn "Only showing smallest 5 out of $(length(gm)) gain margins"
                    idx = sortperm(gm)
                    wgm = wgm[idx[1:5]]
                    gm = gm[idx[1:5]]
                end
                # Let's be reasonable, only plot 5 smallest phase margins
                if length(pm) > 5
                    @warn "Only showing \"smallest\" 5 out of $(length(pm)) phase margins"
                    idx = sortperm(pm)
                    wgm = wpm[idx[1:5]]
                    gm = pm[idx[1:5]]
                end
                if _PlotScale == "dB"
                    mag = 20 .* log10.(1 ./ gm)
                    oneLine = 0
                else
                    mag = 1 ./ gm
                    oneLine = 1
                end
                for k=1:length(wgm)
                    #Plot gain margins
                    Plots.plot!(fig, [wgm[k];wgm[k]], [1;mag[k]]; lab="", subplot=s2i(2i-1,j), group=si)
                end
                #Plot gain line at 1
                Plots.plot!(fig, [w[1],w[end]], [oneLine,oneLine], l=:dash, c=:gray, lab="", subplot=s2i(2i-1,j))
                titles[j,i,1,1] *= "["*join([Printf.@sprintf("%2.2f",v) for v in gm],", ")*"] "
                titles[j,i,1,2] *= "["*join([Printf.@sprintf("%2.2f",v) for v in wgm],", ")*"] "
                for k=1:length(wpm)
                    #Plot the phase margins
                    Plots.plot!(fig, [wpm[k];wpm[k]],[fullPhase[k];fullPhase[k]-pm[k]]; lab="", subplot=s2i(2i,j))
                    #Plot the line at 360*k
                    Plots.plot!(fig, [w[1],w[end]],(fullPhase[k]-pm[k])*ones(2); l=:dash, c=:gray, lab="", subplot=s2i(2i,j))
                end
                titles[j,i,2,1] *=  "["*join([Printf.@sprintf("%2.2f",v) for v in pm],", ")*"] "
                titles[j,i,2,2] *=  "["*join([Printf.@sprintf("%2.2f",v) for v in wpm],", ")*"] "
            end
        end
    end
    for j = 1:nu
        for i = 1:ny
            Plots.title!(fig, titles[j,i,1,1]*" "*titles[j,i,1,2], subplot=s2i(2i-1,j))
            Plots.title!(fig, titles[j,i,2,1]*" "*titles[j,i,2,2], subplot=s2i(2i,j))
        end
    end
    return fig
end

# HELPERS:

function _same_io_dims(systems::LTISystem...)
    sizes = map(size, systems)
    return all(s -> s == sizes[1], sizes)
end

function _default_time_data(systems::Vector{T}) where T<:LTISystem
    sample_times = [_default_dt(i) for i in systems]
    tfinal = 100*maximum(sample_times)
    return sample_times, tfinal
end
_default_time_data(sys::LTISystem) = _default_time_data(LTISystem[sys])

@userplot Pzmap
"""`fig = pzmap(fig, system, args...; kwargs...)`

Create a pole-zero map of the `LTISystem`(s) in figure `fig`, `args` and `kwargs` will be sent to the `scatter` plot command."""
pzmap
@recipe function pzmap(p::Pzmap)
    systems = p.args[1]
    if systems[1].nu + systems[1].ny > 2
        @warn("pzmap currently only supports SISO systems. Only transfer function from u₁ to y₁ will be shown")
    end
    seriestype := :scatter
    framestyle --> :zerolines
    title --> "Pole-zero map"
    legend --> false
    for system in systems
        z,p,k = zpkdata(system)
        if !isempty(z[1])
            @series begin
                markershape --> :c
                markersize --> 15.
                markeralpha --> 0.5
                real(z[1]),imag(z[1])
            end
        end
        if !isempty(p[1])
            @series begin
                markershape --> :x
                markersize := 15.
                real(p[1]),imag(p[1])
            end
        end

        if isdiscrete(system)
            v = range(0,stop=2π,length=100)
            S,C = sin.(v),cos.(v)
            @series begin
                linestyle --> :dash
                linecolor := :black
                grid --> true
                C,S
            end
        end
    end
end
pzmap(sys::LTISystem; kwargs...) = pzmap([sys]; kwargs...)
pzmap!(sys::LTISystem; kwargs...) = pzmap!([sys]; kwargs...)

"""`fig = gangoffourplot(P::LTISystem, C::LTISystem)`, `gangoffourplot(P::Union{Vector, LTISystem}, C::Vector; plotphase=false)`

Gang-of-Four plot.

`kwargs` is sent as argument to Plots.plot."""
function gangoffourplot(P::Union{Vector, LTISystem}, C::Vector, args...; plotphase=false, kwargs...)
    # Array of (S,D,N,T)
    if P isa LTISystem # Don't broadcast over scalar (with size?)
        P = [P]
    end
    sys = gangoffour.(P,C)
    fig = bodeplot([[sys[i][1] sys[i][2]; sys[i][3] sys[i][4]] for i = 1:length(C)], args..., plotphase=plotphase; kwargs...)
    titles = fill("", 1, plotphase ? 8 : 4)
    # Empty titles on phase
    titleIdx = plotphase ? [1,2,5,6] : [1,2,3,4]
    titles[titleIdx] = ["S = 1/(1+PC)", "D = P/(1+PC)", "N = C/(1+PC)", "T = PC/(1+PC)"]
    Plots.plot!(fig, title = titles)
    return fig
end

function gangoffourplot(P::LTISystem,C::LTISystem, args...; plotphase=false, kwargs...)
    gangoffourplot(P,[C], args...; plotphase=plotphase, kwargs...)
end


text:
#md # [![](https://mybinder.org/badge_logo.svg)](@__BINDER_ROOT_URL__/build_literate/Linear_Approximation.ipynb)

# # Linear Approximation
# The following code demonstrates the linear approximation of the generalized Bloch model and replicates Figs. 7 and 8 in the paper.

# For this analysis we need the following packages:
using DifferentialEquations
using BenchmarkTools
using LinearAlgebra
using MRIgeneralizedBloch
using Plots
plotlyjs(bg = RGBA(31/255,36/255,36/255,1.0), ticks=:native); #hide #!nb

# and we simulate a coupled spin system with the following parameters:
m₀ˢ = 0.1
m₀ᶠ = 1-m₀ˢ
R₁ = 1 # 1/s
R₂ᶠ = 1 / 50e-3 # 1/s
Rₓ = 70; # 1/s

# ## Linearized ``T_2^{s,l}``
# We demonstrate the linear approximation at the example of the Green's function corresponding to the [super-Lorentzian lineshape](http://dx.doi.org/10.1002/mrm.1910330404), which we interpolate to improve the perfomance:
G = interpolate_greens_function(greens_superlorentzian, 0, 1e-3 / 5e-6);

# The function [`precompute_R2sl`](@ref) returns another function, `R₂ˢˡ(Tʳᶠ, α, B1, T₂ˢ)`, that interpolates the linearized relaxation rate, as well as functions that describe its derivatives wrt. ``T_2^s`` and ``B_1``, respectively:
R₂ˢˡ, ∂R₂ˢˡ∂T₂ˢ, ∂R₂ˢˡ∂B₁ = precompute_R2sl(TRF_min=5, TRF_max=100, T2s_min=1, T2s_max=1, ω1_max=π/5, B1_max=1, greens=G);
# The derivatives are not used here and are just assigned for demonstration purposes.

# In order to replicate Fig. 7, we plot `R₂ˢˡ(Tʳᶠ, α, B₁, T₂ˢ)` for a varying ``α`` and ``T_\text{RF}/T_2^s``:
α = (0.01:.01:1) * π
TʳᶠoT₂ˢ = 5:100

TʳᶠoT₂ˢ_m = repeat(reshape(TʳᶠoT₂ˢ, 1, :), length(α), 1)
α_m = repeat(α, 1, size(TʳᶠoT₂ˢ_m, 2))

p = plot(xlabel="Tʳᶠ/T₂ˢ", ylabel="α/π", colorbar_title="T₂ˢˡ/T₂ˢ")
contour!(p, TʳᶠoT₂ˢ, α ./ π, 1 ./ R₂ˢˡ.(TʳᶠoT₂ˢ_m, α_m, 1, 1), fill = true)
#md Main.HTMLPlot(p) #hide

#src export
using Printf #src
io = open(expanduser(string("~/Documents/Paper/2021_MT_IDE/Figures/Linearized_T2s.txt")), "w") #src
write(io, "alphaopi TRFoT2s T2sloT2s \n") #src

for αi in α #src
    for TʳᶠoT₂ˢi in TʳᶠoT₂ˢ #src
        write(io, @sprintf("%1.3f %1.3f %1.3f \n", αi/π, TʳᶠoT₂ˢi, 1 / R₂ˢˡ(TʳᶠoT₂ˢi, αi, 1, 1))) #src
    end #src
    write(io, " \n") #src
end #src
close(io) #src

# ## Single RF Pulse
# To replicate Fig. 8a, we simulate and plot the dynamics of a coupled spin system during a single π-pulse, starting from thermal equilibrium.
Tʳᶠ = 100e-6 # s
T₂ˢ = 10e-6 # μs
m0_5D = [0,0,m₀ᶠ,m₀ˢ,1]
mfun(p, t; idxs=nothing) = typeof(idxs) <: Number ? 0 : m0_5D; # intialize history function, here with the ability to just call a single index

# The full generalized Bloch model is solved by
param = (π/Tʳᶠ, 1, 0, m₀ˢ, R₁, R₂ᶠ, Rₓ, R₁, T₂ˢ, G)
prob = DDEProblem(apply_hamiltonian_gbloch!, m0_5D, mfun, (0.0, Tʳᶠ), param)
sol_pi_full = solve(prob);

# and we evaluate the interpolated solution at the following time points:
t = (0:.01:1) * Tʳᶠ # s
Mpi_full = zeros(length(t),4)
for i in eachindex(t)
    Mpi_full[i,:] = sol_pi_full(t[i])[1:4]
end

# Further, we calculate the linear approximation, which is simulated in a 6D-space as it explicitly models ``x^s``:
m0_6D = [0,0,m₀ᶠ,0,m₀ˢ,1]

Mpi_appx = similar(Mpi_full)
for i in eachindex(t)
    H = exp(hamiltonian_linear(π/Tʳᶠ, 1, 0, t[i], m₀ˢ, R₁, Rₓ, R₁, R₂ᶠ, R₂ˢˡ(Tʳᶠ, π, 1, T₂ˢ)))
    Mpi_appx[i,:] = (H * m0_6D)[[1:3;5]]
end

# and plot the original generalized Bloch model and its linear approximation for comparison:
p = plot(xlabel="t [s]", ylabel="m/m₀")
plot!(p, t, Mpi_full[:,1] / m₀ᶠ, label="xᶠ/m₀ᶠ original model")
plot!(p, t, Mpi_appx[:,1] / m₀ᶠ, label="xᶠ/m₀ᶠ linear approximation")
plot!(p, t, Mpi_full[:,3] / m₀ᶠ, label="zᶠ/m₀ᶠ original model")
plot!(p, t, Mpi_appx[:,3] / m₀ᶠ, label="zᶠ/m₀ᶠ linear approximation")
plot!(p, t, Mpi_full[:,4] / m₀ˢ, label="zˢ/m₀ˢ original model")
plot!(p, t, Mpi_appx[:,4] / m₀ˢ, label="zˢ/m₀ˢ linear approximation")
#md Main.HTMLPlot(p) #hide

# We observe slight deviations of zˢ during the pulse, but a virtually perfect match at the end of the RF pulse.

# ## RF Pulses with Different Flip Angles
# To replicate Fig. 8b, we simulate the spin dynamics during multiple RF pulses with different flip angles α, each simulation starting from thermal equilibrium, and analyze the magnetization at the end of each pulse:
α = (.01:.01:1) * π

M_full = zeros(length(α), 4)
M_appx = similar(M_full)
for i in eachindex(α)
    param = (α[i]/Tʳᶠ, 1, 0, m₀ˢ, R₁, R₂ᶠ, Rₓ, R₁, T₂ˢ, G)
    prob = DDEProblem(apply_hamiltonian_gbloch!, m0_5D, mfun, (0.0, Tʳᶠ), param)
    M_full[i,:] = solve(prob)[end][1:4]

    u = exp(hamiltonian_linear(α[i]/Tʳᶠ, 1, 0, Tʳᶠ, m₀ˢ, R₁, R₂ᶠ, Rₓ, R₁, R₂ˢˡ(Tʳᶠ, α[i], 1, T₂ˢ))) * m0_6D
    M_appx[i,:] = u[[1:3;5]]
end

p = plot(xlabel="α/π", ylabel="m/m₀")
plot!(p, α/π, M_appx[:,1] / m₀ᶠ, label="xᶠ/m₀ˢ original model")
plot!(p, α/π, M_full[:,1] / m₀ᶠ, label="xᶠ/m₀ˢ linear approximation")
plot!(p, α/π, M_appx[:,3] / m₀ᶠ, label="zᶠ/m₀ˢ original model")
plot!(p, α/π, M_full[:,3] / m₀ᶠ, label="zᶠ/m₀ˢ linear approximation")
plot!(p, α/π, M_appx[:,4] / m₀ˢ, label="zˢ/m₀ˢ original model")
plot!(p, α/π, M_full[:,4] / m₀ˢ, label="zˢ/m₀ˢ linear approximation")
#md Main.HTMLPlot(p) #hide

# Visually, the linear approximation matches the full simulation well. The normalized root-mean-squared error of the linear approximation for ``x^f`` is
norm(M_appx[:,1] .- M_full[:,1]) / norm(M_full[:,1])
# for ``z^f`` it is
norm(M_appx[:,3] .- M_full[:,3]) / norm(M_full[:,3])
# and for ``z^s`` it is
norm(M_appx[:,4] .- M_full[:,4]) / norm(M_full[:,4])
# which confirms the good concordance.

# ## Benchmark
# We analyze the execution time for solving the full integro-differential equation:
param = (α[end]/Tʳᶠ, 1, 0, m₀ˢ, R₁, R₂ᶠ, Rₓ, R₁, T₂ˢ, G)
prob = DDEProblem(apply_hamiltonian_gbloch!, m0_5D, mfun, (0.0, Tʳᶠ), param)
@benchmark solve($prob)

# The `$` symbol *interpolates* the variable, which improves the accuracy of the timing measurement. We can compare this time to the time it takes to calculate the linear approximation, including the time it takes to evaluate the interpolated `R₂ˢˡ`:
@benchmark exp(hamiltonian_linear($(α[end]/Tʳᶠ), 1, 0, $Tʳᶠ, $m₀ˢ, $R₁, $R₂ᶠ, $Rₓ, $R₁, R₂ˢˡ($Tʳᶠ, $α[end], 1, $T₂ˢ))) * $m0_6D

# We can see that linear approximation is about 4 orders of magnitude faster compared to the full model.

#src export data for t vs. M figure
io = open(expanduser(string("~/Documents/Paper/2021_MT_IDE/Figures/Linearized_gBloch_M_during_Pulse.txt")), "w") #src
write(io, "t_us xf_full xf_appx zf_full zf_appx zs_full zs_appx \n") #src
for i = 1:3:length(t) #src
    write(io, @sprintf("%1.3f ", t[i]*1e6)) #src
    write(io, @sprintf("%1.3f ", Mpi_full[i,1] / m₀ᶠ)) #src
    write(io, @sprintf("%1.3f ", Mpi_appx[i,1] / m₀ᶠ)) #src
    write(io, @sprintf("%1.3f ", Mpi_full[i,3] / m₀ᶠ)) #src
    write(io, @sprintf("%1.3f ", Mpi_appx[i,3] / m₀ᶠ)) #src
    write(io, @sprintf("%1.3f ", Mpi_full[i,4] / m₀ˢ)) #src
    write(io, @sprintf("%1.3f ", Mpi_appx[i,4] / m₀ˢ)) #src
    write(io, "\n") #src
end #src
close(io) #src

#src export data for α vs. M figure
io = open(expanduser(string("~/Documents/Paper/2021_MT_IDE/Figures/Linearized_gBloch_vary_alpha.txt")), "w") #src
write(io, "alpha/pi xf_full xf_appx zf_full zf_appx zs_full zs_appx \n") #src
for i = 1:3:length(α) #src
    write(io, @sprintf("%1.3f ", α[i]/π)) #src
    write(io, @sprintf("%1.3f ", M_full[i,1] / m₀ᶠ)) #src
    write(io, @sprintf("%1.3f ", M_appx[i,1] / m₀ᶠ)) #src
    write(io, @sprintf("%1.3f ", M_full[i,3] / m₀ᶠ)) #src
    write(io, @sprintf("%1.3f ", M_appx[i,3] / m₀ᶠ)) #src
    write(io, @sprintf("%1.3f ", M_full[i,4] / m₀ˢ)) #src
    write(io, @sprintf("%1.3f ", M_appx[i,4] / m₀ˢ)) #src
    write(io, "\n") #src
end #src
close(io) #src

ext: jl


text:
### A Pluto.jl notebook ###
# v0.14.0

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ d84dbb08-b1a4-11eb-2d3b-0711fddd1347
html"""

Section 3.6

 Snowball Earth and hysteresis

Lecture Video

"""

# ╔═╡ a0b3813e-adab-11eb-2983-616cf2bb6f5e
begin
    import Pkg
    Pkg.activate(mktempdir())
    Pkg.add([
        Pkg.PackageSpec(name="DifferentialEquations", version="6"),
        Pkg.PackageSpec(name="Plots", version="1"),
        Pkg.PackageSpec(name="PlutoUI", version="0.7"),
    ])
    using DifferentialEquations, Plots, PlutoUI, LinearAlgebra
end

# ╔═╡ ef8d6690-720d-4772-a41f-b260d306b5b2
TableOfContents(title="📚 Table of Contents", indent=true, depth=4, aside=true)

# ╔═╡ 26e1879d-ab57-452a-a09f-49493e65b774
md"""
# Julia Concepts

- `sign(x)`

"""

# ╔═╡ fd12468f-de16-47cc-8210-9266ca9548c2
md"""
## The sign or signum(latin for "sign") function
"""

# ╔═╡ 30969341-9079-4732-bf55-d6bba2c2c16c
md"""
`sign(x)` returns 0 if x is 0, or ±1 for positive/negative x.
"""

# ╔═╡ 5bf3fe83-d096-4df1-8476-0a6500b01868
begin
	scatter(sign, -5:.1:5, legend=false, m=:c, ms=3, size=(600,300))
	title!("The sign function is discontinuous at 0")
	xlabel!("x")
	ylabel!("sign(x)")
end

# ╔═╡ af7f36d9-adca-48b8-95bb-ac620e6f1b4f
sign(Inf)

# ╔═╡ 790add0f-c83f-4824-92ae-53159ce58f64
sign(-Inf)

# ╔═╡ 21210cfa-0366-4019-86f7-158fdd5f21ad
md"""
# Mathematics:multiple equilibria.
## Using computation to explore hysteresis.
"""

# ╔═╡ 978e5fc0-ddd1-4e93-a243-a95d414123b9
md"""
The function  $f(y,a) = {\rm sign}(y) + a -y$ can be written

``f(y,a)= \left\{ \begin{array}{l} -1+a-y \ \ {\rm if} \ y<0  \ \
\textrm{(root at } a-1 \textrm{ if } a<1 )\\
\ \  \ 1+a-y \ \ {\rm if} \ y>0
\ \ \textrm{(root at } a+1 \textrm{ if } a>-1 )
\end{array} \right.``

(we will ignore y=0).  Notice that for -1<a<1 there are two roots.
"""

# ╔═╡ 6139554e-c6c9-4252-9d64-042074f68391
begin
	f(y,a) =   sign(y)  + a -  y
	f(y,a,t) = f(y,a) # just for the difeq solver
end

# ╔═╡ e115cbbc-9d49-4fa1-8701-fa48289a0916
md"""
The graph of `z=f(y,a)`  consists of two parallel half-planes. On the left below we intersect that graph with planes of constant `a`.  On the right, we have the intersection with `z=0`.
"""

# ╔═╡ bd65e7f9-ecf2-43ac-b5a8-99b03866a5c8
begin
	plot( y->f(y,a, 0), -9, -eps(), xlabel="y" , ylabel="sign(y)-y+a", lw=3,color=:blue)
	plot!(y->f(y,a, 0), eps() ,9, lw=3, color=:blue)
	ylims!(-5,10)
	xlims!(-9,9)
	hline!([0], ls=:dash, color=:pink, legend=false)

	if a< 1

	  annotate!( -2.2+a, -4, text(round(-1+a,digits=1), position=:right,11,:red))
	  scatter!([-1+a],[0],m=:circle,ms=10,color=:red)
	end
	if a>-1

	  annotate!( 1+a, -4, text(round(1+a,digits=1), position=:right,11,:blue))
		 scatter!([1+a],[0],m=:circle,ms=10,color=:blue)
	end
	vline!([ 1+a], ls=:dash, color=:gray)
	 vline!([-1+a], ls=:dash, color=:gray)
	vline!([0],ls=:dash, color=:pink)
	P1 = title!(" Plot of f(y,a=$a) vs y")

	begin
	plot(a->a-1, -5:1, lw=2, c=:blue)
		 annotate!( a, 3+a, "a=$a")
	plot!(a->a+1, -1:5, lw=2,  c=:blue, legend=false)
	xlabel!("a")
	ylabel!("y")
	title!("Solution to f(y,a)=0")

		if a<1
			annotate!(-4,  -.7+a,  text(round(-1+a,digits=1), position=:center,11,:red))
			hline!([-1+a],c=:gray,ls=:dash)
			scatter!([a],[-1+a],m=:circle,ms=10,color=:red)
		end
		if a> -1
						annotate!(-4,  1.3+a,  text(round(1+a,digits=1), position=:center,11,:blue))
			hline!([1+a],c=:gray,ls=:dash)
			scatter!([a],[1+a],m=:circle,ms=10,color=:blue)
		end
		ylims!(-5,5)
	P2 = annotate!(0,-4.7, text("Two roots if -1<a<1",6,:blue) )
		vline!([a],c=:gray,ls=:dash)
	plot( P1 ,P2,layout=(1,2))
end
end

# ╔═╡ 991f14bb-fc4a-4505-a3be-5ced2fb148b6

# ╔═╡ 61960e99-2932-4a8f-9e87-f64a7a043489
md"""
a = $(@bind a Slider(-5:.1:5, show_value=true, default=0) )

"""

# ╔═╡ 9acf9b8a-7fde-46c7-bf0c-c4dedc5a064b
md"""
The above are two section views of the surface z=f(y,a). Left: constant a.  Right: z =0.
"""

# ╔═╡ 9bc33f22-c065-4c1d-b06f-78c70415a111
# begin
# 	#s1=surface(-5:.1:-.1, -4:.1:4, f, legend=false )
# 	#s2=surface!(.00001:.1:4.1, -4:.1:4, f, legend=false )
# 	surface(-6:.05:4, -5:.05:4, f, legend=false, alpha=.5 )

# 	plot!(  [0,-6 ], [1, -5  ], [0,0], c=:red)
# 	plot!(  [0,5 ], [-1, 4  ], [0,0], c=:red)
# 	xlabel!("y")
# 	ylabel!("a")

# end

# ╔═╡ 6f61180d-6900-48fa-998d-36110e79d2dc
gr()

# ╔═╡ a09a0064-c6ab-4912-bc9f-96ab72b8bbca
sol = solve(  ODEProblem( f, y₀, (0,10.0), a ));

# ╔═╡ 5027e1f8-8c50-4538-949e-6c95c550016e
md"""
### Solution to y' = f(y,a)
with y(0)=y₀
"""

# ╔═╡ 465f637c-0555-498b-a881-a2f6e5714cbb
md"""
y₀ = $(@bind y₀ Slider(-6:.1:6, show_value=true, default=2.0) )
"""

# ╔═╡ a4686bca-90d6-4e02-961c-59f08fc37553
plotit(y₀, a)

# ╔═╡ f8ee2373-6af0-4d81-98fb-23bde10198ef
function plotit(y₀, a)

	sol = solve(  ODEProblem( f, y₀, (0,10.0), a ));
	p = plot( sol , legend=false, background_color_inside=:black , ylims=(-7,7), lw=3, c=:red)

	# plot direction field:
		xs = Float64[]
		ys = Float64[]

	    lrx = LinRange( xlims(p)..., 30)
		for x in  lrx
			for y in LinRange( ylims(p)..., 30)
				v = [1, f(y, a ,x) ]
				v ./=  (100* (lrx[2]-lrx[1]))

				push!(xs, x - v[1], x + v[1], NaN)
				push!(ys, y - v[2], y + v[2], NaN)

			end
		end

	if a<1 hline!( [-1+a],c=:white,ls=:dash); end
	if a>-1  hline!( [1+a],c=:white,ls=:dash);end
	    hline!( [0 ],c=:pink,ls=:dash)
		plot!(xs, ys, alpha=0.7, c=:yellow)
	    ylabel!("y")
	    annotate!(-.5,y₀,text("y₀",color=:red))
	if a>-1
	   annotate!(5,1+a,text(round(1+a,digits=3),color=:white))
	end
	if a<1
	   annotate!(5,-1+a-.4,text(round(-1+a,digits=3),color=:white))
	end
    title!("Solution to y'=f(y,a)")
return(p)
end

# ╔═╡ a94b5160-f4bf-4ddc-9ee6-581ea20c8b90
md"""
### Hysteresis: Increasing then decreasing ``a``
Let's increase a by .25 from -4 to 4 then decrease from -4 to 4.
Every time we change a, we let 10 units of time evolve, enough
to reach the equilibriumf for that a, and watch the y values.

We see that when -1<a<1 it's possible to be at the "negative" equilibrium
or the "positive" equilibrium, depending on how you got there.
"""

# ╔═╡ fd882095-6cc4-4927-967c-6b02d5b1ad95
let
    Δt = 10
    t = 0.0
	y₀ =-5
	a = -4
	sol = solve(  ODEProblem( f, y₀, (t,t+Δt), a ))
	p=plot(sol)
	annotate!(10,-5,text("a=",7))

	for i= 1:32
		a += .25
		t += Δt
		y₀ = sol(t)
		sol = solve(  ODEProblem( f, y₀, (t,t+Δt), a ))

		if -1≤a ≤1
		  p=plot!(sol,xlims=(0,t), legend=false, fillrange=-5,fillcolor=:gray,alpha=.5)
		else
		  p=plot!(sol,xlims=(0,t), legend=false)
		end

		if i%4==0
			annotate!(t,-5,text( round(a,digits=1), 7, :blue))
		end

	end
	for i= 1:32
		a -= .25
		t += Δt
		y₀ = sol(t)
		sol = solve(  ODEProblem( f, y₀, (t,t+Δt), a ))

		if -1≤a≤1
		  p=plot!(sol,xlims=(0,t), legend=false, fillrange=-5,fillcolor=:gray,alpha=.5)
		else
		  p=plot!(sol,xlims=(0,t), legend=false)
		end

		if i%4==0
			annotate!(t,-5,text( round(a,digits=1), 7, :red))
		end
	end
	ylabel!("y")

	as_svg(plot!())
end

# ╔═╡ 10693e53-3741-4388-b3b1-eba739ec01d0
md"""
 The dependence of the state of a system on its history, what we observe above,
is known as *hysteresis* (Greek (ὑστέρησις) for lagging behind).
"""

# ╔═╡ 11b3fb9e-5922-4350-9424-51fba33502d4
md"""
# Application to Snowball Earth, the ice-albedo feedback
(from Henri Drake's lecture)
##  Review of the last climate lecture.

Recall from [Lecture 20 (Part I)](https://www.youtube.com/watch?v=Gi4ZZVS2GLA&t=15s) that the the **zero-dimensional energy balance equation** is

\begin{gather}
\color{brown}{C \frac{dT}{dt}}
\; \color{black}{=} \; \color{orange}{\frac{(1 - α)S}{4}}
\; \color{black}{-} \; \color{blue}{(A - BT)}
\; \color{black}{+} \; \color{grey}{a \ln \left( \frac{[\text{CO}₂]}{[\text{CO}₂]_{\text{PI}}} \right)},
\end{gather}
"""

# ╔═╡ b71fca45-9687-4a51-8e1c-1f413e83e58d
html""""""

# ╔═╡ d993a2fc-2319-4f64-8a17-904a57593da2
md"""
Today, we will ignore changes in CO₂, so that

$\ln \left( \frac{ [\text{CO}₂]_{\text{PI}} }{[\text{CO}₂]_{\text{PI}}} \right) = \ln(1) = 0$

and the model simplifies to

\begin{gather}
\color{brown}{C \frac{dT}{dt}}
\; \color{black}{=} \; \color{orange}{\frac{(1 - α)S}{4}}
\; \color{black}{-} \; \color{blue}{(A - BT)}.
\end{gather}

The dynamics of this **Ordinary Differential Equation (ODE)** are quite simple because it is *linear*: we can rewrite it in the form

$\dot{T} = f(T(t))$ where $f(x) = \alpha x + \beta$ is a *linear* function of x. A linear ODE permits only one stable solution, $\dot{T} = 0$, which in Lecture 20 we found was Earth's pre-industrial temperature $T_{0} = 14$°C.

In this lecture, we show how a small modification that makes one term in our simple climate model non-linear completely changes its dynamics, allowing us to explain the existence of both "Snowball Earth" and the relatively warm pre-industrial climate that allowed humans to thrive.
"""

# ╔═╡ 4bdc0f0c-e696-4d87-b10c-8a0da9a0ee5b
md"""
## 1) Background: Snowball Earth

Geological evidence shows that the Neoproterozoic Era (550 to 1000 million years ago) is marked by two global glaciation events, in which Earth's surface was covered in ice and snow from the Equator to the poles (see review by [Pierrehumbert et al. 2011](https://www.annualreviews.org/doi/full/10.1146/annurev-earth-040809-152447)).
"""

# ╔═╡ 7b7b631e-2ba3-4ed3-bad0-ec6ecb70ad49
html"""

"""

# ╔═╡ 066743eb-c890-40b9-9f6b-9f79b7ebcbd2
md"""### 1.1) The ice-albedo feedback

In Lecture 20, we used a **constant** value $α =$ $(Model.hist.α) for Earth's planetary albedo, which is a reasonable thing to do for small climate variations relative to the present (such as the difference between the present-day and preindustrial climates). In the case of large variations, however, this approximation is not very reliable.

While oceans are dark and absorbant, $α_{ocean} \approx 0.05$, ice and snow are bright and reflective: $\alpha_{ice,\,snow} \approx 0.5$ to $0.9$. Thus, if much of the ocean's surface freezes over, we expect Earth's albedo to rise dramatically, causing more sunlight to be reflected to space, which in turn causes even more cooling and more of the ocean to freeze, etc. This *non-linear positive feedback effect* is referred to as the **ice-albedo feedback** (see illustration below).
"""

# ╔═╡ 70ec6ae9-601f-4862-96cb-f251d4b5a7fd
html"""

"""

# ╔═╡ 2bafd1a4-32a3-4787-807f-0a5132d66c28
md"""
We can represent the ice-albedo feedback crudely in our energy balance model by allowing the albedo to depend on temperature:

$\alpha(T) = \begin{cases}
\alpha_{i} & \mbox{if }\;\; T \leq -10\text{°C} &\text{(completely frozen)}\\
\alpha_{i} + (\alpha_{0}-\alpha_{i})\frac{T + 10}{20} & \mbox{if }\;\; -10\text{°C} \leq T \leq 10\text{°C} &\text{(partially frozen)}\\
\alpha_{0} &\mbox{if }\;\; T \geq 10\text{°C} &\text{(no ice)}
\end{cases}$
"""

# ╔═╡ 40b5e447-0cfb-4f35-8f95-6aa29793e5ad
begin
	T_example = -20.:1.:20.
	plot(size=(600, 400), ylims=(0.2, 0.6))
	plot!([-20, -10], [0.2, 0.2], fillrange=[0.6, 0.6], color=:lightblue, alpha=0.2, label=nothing)
	plot!([10, 20], [0.2, 0.2], fillrange=[0.6, 0.6], color=:red, alpha=0.12, label=nothing)
	plot!(T_example, α.(T_example), lw=3., label="α(T)", color=:black)
	plot!(ylabel="albedo α\n(planetary reflectivity)", xlabel="Temperature [°C]")
	annotate!(-15.5, 0.252, text("completely\nfrozen", 10, :darkblue))
	annotate!(15.5, 0.252, text("no ice", 10, :darkred))
	annotate!(-0.3, 0.252, text("partially frozen", 10, :darkgrey))

end

# ╔═╡ fca6c4ec-4d0c-4f97-b966-ce3a81a18710
md"""
### 1.2) Adding the ice-albedo feedback to our simple climate model

First, we program albedo as a function of temperature.
"""

# ╔═╡ b1fee17b-6522-4cf0-a614-5ff8aa8f8614
function α(T; α0=Model.α, αi=0.5, ΔT=10.)
	if T < -ΔT
		return αi
	elseif -ΔT <= T < ΔT
		return αi + (α0-αi)*(T+ΔT)/(2ΔT)
	elseif T >= ΔT
		return α0
	end
end

# ╔═╡ cfde8137-cfcd-46de-9c26-8abb64b6b3a9
md"""
To add this function into our energy balance model from [Lecture 20 (Part I)](https://www.youtube.com/watch?v=Gi4ZZVS2GLA&t=15s) (which we've copied into the cell below), all we have to do is overwrite the definition of the `timestep!` method to specify that the temperature-dependent albedo should be updated based on the current state:
"""

# ╔═╡ 4351b05f-f9bf-4046-9f95-a0a56b1e8cc9
module Model

const S = 1368; # solar insolation [W/m^2]  (energy per unit time per unit area)
const α = 0.3; # albedo, or planetary reflectivity [unitless]
const B = -1.3; # climate feedback parameter [W/m^2/°C],
const T0 = 14.; # preindustrial temperature [°C]

absorbed_solar_radiation(; α=α, S=S) = S*(1 - α)/4; # [W/m^2]
outgoing_thermal_radiation(T; A=A, B=B) = A - B*T;

const A = S*(1. - α)/4 + B*T0; # [W/m^2].

greenhouse_effect(CO2; a=a, CO2_PI=CO2_PI) = a*log(CO2/CO2_PI);

const a = 5.0; # CO2 forcing coefficient [W/m^2]
const CO2_PI = 280.; # preindustrial CO2 concentration [parts per million; ppm];
CO2_const(t) = CO2_PI; # constant CO2 concentrations

const C = 51.; # atmosphere and upper-ocean heat capacity [J/m^2/°C]

function timestep!(ebm)
	append!(ebm.T, ebm.T[end] + ebm.Δt*tendency(ebm));
	append!(ebm.t, ebm.t[end] + ebm.Δt);
end;

tendency(ebm) = (1. /ebm.C) * (
	+ absorbed_solar_radiation(α=ebm.α, S=ebm.S)
	- outgoing_thermal_radiation(ebm.T[end], A=ebm.A, B=ebm.B)
	+ greenhouse_effect(ebm.CO2(ebm.t[end]), a=ebm.a, CO2_PI=ebm.CO2_PI)
);

begin
	mutable struct EBM
		T::Array{Float64, 1}

		t::Array{Float64, 1}
		Δt::Float64

		CO2::Function

		C::Float64
		a::Float64
		A::Float64
		B::Float64
		CO2_PI::Float64

		α::Float64
		S::Float64
	end;

	# Make constant parameters optional kwargs
	EBM(T::Array{Float64, 1}, t::Array{Float64, 1}, Δt::Float64, CO2::Function;
		C=C, a=a, A=A, B=B, CO2_PI=CO2_PI, α=α, S=S) = (
		EBM(T, t, Δt, CO2, C, a, A, B, CO2_PI, α, S)
	);

	# Construct from float inputs for convenience
	EBM(T0::Float64, t0::Float64, Δt::Float64, CO2::Function;
		C=C, a=a, A=A, B=B, CO2_PI=CO2_PI, α=α, S=S) = (
		EBM([T0], [t0], Δt, CO2;
			C=C, a=a, A=A, B=B, CO2_PI=CO2_PI, α=α, S=S);
	);
end;

begin
	function run!(ebm::EBM, end_year::Real)
		while ebm.t[end] < end_year
			timestep!(ebm)
		end
	end;

	run!(ebm) = run!(ebm, 200.) # run for 200 years by default
end

CO2_hist(t) = CO2_PI * (1 .+ fractional_increase(t));
fractional_increase(t) = ((t .- 1850.)/220).^3;

begin
	hist = EBM(T0, 1850., 1., CO2_hist,
		C=C, B=B, A=A
	)
	run!(hist, 2020.)
end

begin
	CO2_RCP26(t) = CO2_PI * (1 .+ fractional_increase(t) .* min.(1., exp.(-((t .-1850.).-170)/100))) ;
	RCP26 = EBM(T0, 1850., 1., CO2_RCP26)
	run!(RCP26, 2100.)

	CO2_RCP85(t) = CO2_PI * (1 .+ fractional_increase(t) .* max.(1., exp.(((t .-1850.).-170)/100)));
	RCP85 = EBM(T0, 1850., 1., CO2_RCP85)
	run!(RCP85, 2100.)
end

end

# ╔═╡ e9942719-93cc-4203-8d37-8f91539104b1
function Model.timestep!(ebm)
	ebm.α = α(ebm.T[end]) # Added this line
	append!(ebm.T, ebm.T[end] + ebm.Δt*Model.tendency(ebm));
	append!(ebm.t, ebm.t[end] + ebm.Δt);
end

# ╔═╡ 0c1e3051-c491-4de7-a149-ce81c53f5841
md"""## 2) Multiple Equilibria
**OR: the existence of "alternate Earths"**

Human civilization flourished over the last several thousand years in part because Earth's global climate has been remarkably stable and forgiving. The preindustrial combination of natural greenhouse effect and incoming solar radiation yielded temperatures between the freezing and boiling points of water across most of the planet, allowing ecoystems based on liquid water to thrive.

The climate system, however, is rife with non-linear effects like the **ice-albedo effect**, which reveal just how fragile our habitable planet is and just how unique our stable pre-industrial climate was.

We learned in Lecture 20 that in response to temperature fluctuations, *net-negative feedbacks* act to restore Earth's temperature back towards a single equilibrium state in which absorbed solar radiation is balanced by outgoing thermal radiation. Here, we explore how *non-linear positive feedbacks* can temporarily result in a *net-positive feedback* and modify Earth's state space.
"""

# ╔═╡ d0d43dc3-4cda-4602-90c6-2d14a1e63871
md"""
### 2.1) Exploring the non-linear ice-albedo feedback

In [Lecture 20 (Part II)](https://www.youtube.com/watch?v=D3jpfeQCISU), we learned how introducing non-linear terms in *ordinary differential equations* can lead to complex state spaces that allow for multiple fixed points (e.g. for $\dot{x} = \mu + x^{2}$).

Let's explore how this plays out with the non-linear ice-albedo feedback by varying the initial condition $T_{0} \equiv T(t=0)$ and allowing the system to evolve for 200 years.
"""

# ╔═╡ 48431fc9-413f-4f72-8c6f-f48b42c29475
begin
	T0_slider = @bind T0_interact Slider(-60.:0.0025:30., default=24., show_value=true)
	md""" T₀ = $(T0_slider) °C"""
end

# ╔═╡ c719ff5f-421d-4d4d-87b7-34879ab188c5
begin
	ebm_interact = Model.EBM(Float64(T0_interact), 0., 1., Model.CO2_const)
	Model.run!(ebm_interact, 200)
end

# ╔═╡ 0b6bfab9-96e9-4d7a-a486-4582d7303244
begin
	p_interact = plot(ebm_interact.t, ebm_interact.T, label=nothing, lw=3)
	plot!([0.], [T0_interact], label=nothing, markersize=4, markershape=:circle)

	add_labels!(p_interact)
end |> as_svg

# ╔═╡ 9b6edae8-fbf0-4979-9536-c0f782ba70a7
md"We can get an overview of the behavior by drawing a set of these curves all on the same graph:"

# ╔═╡ 81d823a6-5c92-496a-96f1-a0f4762f1f05
begin
	p_equil = plot(xlabel="year", ylabel="temperature [°C]", legend=:bottomright, xlims=(0,205), ylims=(-60, 30.))

	plot!([0, 200], [-60, -60], fillrange=[-10., -10.], fillalpha=0.3, c=:lightblue, label=nothing)
	annotate!(120, -20, text("completely frozen", 10, :darkblue))

	plot!([0, 200], [10, 10], fillrange=[30., 30.], fillalpha=0.09, c=:red, lw=0., label=nothing)
	annotate!(120, 25, text("no ice", 10, :darkred))
	for T0_sample in (-60.:5.:30.)
		ebm = Model.EBM(T0_sample, 0., 1., Model.CO2_const)
		Model.run!(ebm, 200)

		plot!(p_equil, ebm.t, ebm.T, label=nothing)
	end

	T_un = 7.5472
	for δT in 1.e-2*[-2, -1., 0., 1., 2.]
		ebm_un = Model.EBM(T_un+δT, 0., 1., Model.CO2_const)
		Model.run!(ebm_un, 200)

		plot!(p_equil, ebm_un.t, ebm_un.T, label=nothing, linestyle=:dash)
	end

	plot!(p_equil, [200], [Model.T0], markershape=:circle, label="Our pre-industrial climate (stable ''warm'' branch)", color=:orange, markersize=8)
	plot!(p_equil, [200], [-38.3], markershape=:circle, label="Alternate universe pre-industrial climate (stable ''cold'' branch)", color=:aqua, markersize=8)
	plot!(p_equil, [200], [T_un], markershape=:diamond, label="Impossible alternate climate (unstable branch)", color=:lightgrey, markersize=8, markerstrokecolor=:white, alpha=1., markerstrokestyle=:dash)
	p_equil
end

# ╔═╡ 28924aef-9157-4490-afa5-7f232a5101f0
md"We see that for T₀ ⪆  $(round(T_un, digits=2)) °C, all of the curves seem to converge on the T = 14°C *equilibrium* (or *fixed point*) that we saw in Lecture 20. Curves that start below this value warm up and while curves that start above this value will cool down. For T₀ ⪅ $(round(T_un, digits=2)) °C, however, the temperatures converge on a much colder equilibrium around T = -40°C. This is the **Snowball Earth** equilibrium. These two states are referred to as **stable equilibria** because even if the state gets temporarily pushed slightly away from its equilibrium, it will eventually converge right back to its equilibrium.

So what happens is T₀ ≈ $(round(T_un, digits=2)) °C? For some exact temperature near there, there is indeed an equilibrim state: if you start with that temperature you will stay there forever. However, if the temperature starts off even *one one-hundredth of a degree* above or below this exact value, we see that temperatures eventually converge to one of the other two equilibria. Thus, we call this intermediate equilibrium an **unstable equilibrium**, because any *infinitesimal* push away will cause it to careen away towards another state.
"

# ╔═╡ 5eb888e1-ba95-4aa1-97c0-784dc9d9e6d5
md"""
### 2.2) Radiative stability analysis

We can understand why our model has two stable equilibria and one unstable equilibrium by applying concepts from dynamical systems theory.

Recall that, with fixed CO₂ concentrations, our energy balance model differential equation can be expressed as:

$C\frac{d\,T}{d\,t} = \text{ASR}(T) - \text{OTR}(T),$

where now the Absorbed Solar Radiation (ASR) is also temperature dependent because the albedo $α(T)$ is.

In particular, by plotting the right-hand-side tendency terms as a function of the state variable $T$, we can plot a stability diagram for our system that tells us whether the planet will warm ($C \frac{d\,T}{d\,t} > 0$) or cool ($C \frac{d\,T}{d\,t} < 0$).
"""

# ╔═╡ 18922c20-62dc-4524-8300-bbab4db828a9
begin
	T_samples = -60.:1.:30.
	OTR = Model.outgoing_thermal_radiation.(T_samples)
	ASR = [Model.absorbed_solar_radiation.(α=α(T_sample)) for T_sample in T_samples]
	imbalance = ASR .- OTR
end;

# ╔═╡ 22328073-bb46-4ec2-9fdf-17c0daff5741
begin
	p1_stability = plot(legend=:topleft, ylabel="energy flux [W/m²]", xlabel="temperature [°C]")
	plot!(p1_stability, T_samples, OTR, label="Outgoing Thermal Radiation", color=:blue, lw=2.)
	plot!(p1_stability, T_samples, ASR, label="Absorbed Solar Radiation", color=:orange, lw=2.)

	p2_stability = plot(ylims=(-50, 45), ylabel="energy flux [W/m²]", xlabel="temperature [°C]")
	plot!([-60., 30.], [-100, -100], fillrange=[0, 0], color=:blue, alpha=0.1, label=nothing)
	plot!([-60., 30.], [100, 100], fillrange=[0, 0], color=:red, alpha=0.1, label=nothing)
	annotate!(-58, -40, text("cooling", :left, :darkblue))
	annotate!(-58, 38, text("warming", :left, :darkred))
	plot!(p2_stability, T_samples, imbalance, label="Radiative Imbalance\n(ASR - OTR)", color=:black, lw=2.)
	plot!([7.542], [0], markershape=:diamond, markersize=6, color=:lightgrey, markerstrokecolor=:darkgrey, label=nothing)
	plot!([Model.T0], [0], markershape=:circle, markersize=6, color=:orange, markerstrokecolor=:black, label=nothing)
	plot!([-38.3], [0], markershape=:circle, markersize=6, color=:aqua, markerstrokecolor=:black, label=nothing)

	p_stability = plot(p1_stability, p2_stability, layout=(1,2), size=(680, 300))
end

# ╔═╡ 9c4ac33a-ebbd-47db-9057-91624b0a2497
md"
## 3) Transitioning to and from Snowball Earth

### 3.1) Turning up the Sun

Over the entire history of the Earth, the Sun is thought to have brightened by about 40%."

# ╔═╡ 2798028c-d971-45e4-9484-bdec7e8dc048
html"""

"""

# ╔═╡ 5fc03f0f-ae12-476f-93f8-6285ee7f5fc9
md"In the Neoproterozoic (~700 million years ago), the Sun was 93% as bright as it is today, such that the incoming solar radiation was $S =$ 1272 W/m², Earth's average temperature plunged to $T = -50$°C, and Earth's ice-covered surface had a high albedo (reflectivity) of $α_{i} = 0.5$.

### 3.2) Did the increasing brightness of the Sun melt the Snowball?
If we start out in the Neoproterozoic climate and all we do is increase solar insolation to today's value of $S =$ 1368 W/m², can we warm the planet up to the pre-industrial temperature of $T=14$°C?
"

# ╔═╡ eb0307f1-c828-41f9-bc49-2aedd5b588e1
begin
	md"""
	*Move the slider below to change the brightness of the Sun (solar insolation):* S = $(S) [W/m²]
	"""
end

# ╔═╡ ada47602-2cf6-4ef5-b1ba-0b1c6f81ad91
begin
	if extend_S
		solarSlider = @bind S Slider(Smin:2.:Smax, default=Sneo);
		md""" $(Smin) W/m² $(solarSlider) $(Smax) W/m²"""
	else
		solarSlider = @bind S Slider(Smin:2.:Smax_limited, default=Sneo);
		md""" $(Smin) W/m² $(solarSlider) $(Smax_limited) W/m²"""
	end
end

# ╔═╡ 5ebe385f-3542-4674-b419-fb54519aa945
begin
	S
	warming_mask = (1:size(Svec)[1]) .< size(Svec)./2
	p = plot(xlims=(Smin, Smax_limited), ylims=(-55, 75), title="Earth's solar insolation bifurcation diagram")
	# plot!([Model.S, Model.S], [-55, 75], color=:yellow, alpha=0.3, lw=8, label="Pre-industrial / present insolation")
	# if extend_S
	# 	plot!(p, xlims=(Smin, Smax))
	# 	if show_cold
	# 		plot!(Svec[warming_mask], Tvec[warming_mask], color=:blue,lw=3., alpha=0.5, label="cool branch")
	# 	end
	# 	if show_warm
	# 		plot!(Svec[.!warming_mask], Tvec[.!warming_mask], color="red", lw=3., alpha=0.5, label="warm branch")
	# 	end
	# 	if show_unstable
	# 		plot!(S_unstable, T_unstable, color=:darkgray, lw=3., alpha=0.4, ls=:dash, label="unstable branch")
	# 	end
	# end
	# plot!(legend=:topleft)
	# plot!(xlabel="solar insolation S [W/m²]", ylabel="Global temperature T [°C]")
	# plot!([Model.S], [Model.T0], markershape=:circle, label="Our preindustrial climate", color=:orange, markersize=8)
	# plot!([Model.S], [-38.3], markershape=:circle, label="Alternate preindustrial climate", color=:aqua, markersize=8)
	# plot!([Sneo], [Tneo], markershape=:circle, label="neoproterozoic (700 Mya)", color=:lightblue, markersize=8)
	plot_trajectory!(p, reverse(Straj), reverse(Ttraj), lw=9)

	plot!([Smin, Smax], [-60, -60], fillrange=[-10., -10.], fillalpha=0.3, c=:lightblue, label=nothing)
	annotate!(Smin+12, -19, text("completely\nfrozen", 10, :darkblue, :left))

	plot!([Smin, Smax], [10, 10], fillrange=[80., 80.], fillalpha=0.09, c=:red, lw=0., label=nothing)
	annotate!(Smin+12, 15, text("no ice", 10, :darkred, :left))
end |> as_svg

# ╔═╡ 72e822ef-0256-4612-8b6d-654ba6f6a83f
md"""
For insolations $S$ between $(Smin) W/m² and $(Smax_limited) W/m², temperatures always remain below -10°C and the planet remains completely frozen. What if we extend the upper limit on insolation so that the Sun becomes bright enough to start melting ice?
"""

# ╔═╡ 10815b22-b920-4cb7-a078-4d5ef457ebf2
md"""*Extend upper-limit of insolation* $(@bind extend_S CheckBox(default=false))"""

# ╔═╡ 984f0489-7430-43b9-8926-9d7c32e7da63
if extend_S
	md"""
	*"Cold" branch* $(@bind show_cold CheckBox(default=false))    |
	*"Warm" branch* $(@bind show_warm CheckBox(default=false))    |
	*Unstable branch* $(@bind show_unstable CheckBox(default=false))
	"""
else
	show_cold = true;
	nothing
end

# ╔═╡ bf0dd94d-e861-478a-9bb0-d34d17405fcc
if extend_S
	md"""
	##### Abrupt climate transitions

	In this model, temperature variations are fairly smooth unless temperatures rise above -10°C or fall below 10°C, in which case the *ice-albedo positive feedback* kicks in and causes an **abrupt climate transition**. While this is just a simple hypothetical model, these kinds of abrupt climate transitions show up all the time in the paleoclimate record and in more realistic climate model simulations.

	![](https://www.pnas.org/content/pnas/115/52/13288/F1.large.jpg)

	This simulation teaches us that we should not take the stability of our climate for granted and that pushing our present climate outside of its historical regime of stability could trigger catastrophic abrupt climate transitions.
	"""
end

# ╔═╡ 81088439-312a-4042-bcee-c021e5e3f6b7
md"""
### 3.3) If not the Sun, how did Snowball Earth melt?

The leading theory is that a slow but steady outgassing of CO₂ from volcanoes eventually caused a strong enough greenhouse gas effect to offset the cooling effect of the frozen surface's high albedo and raise temperatures above the melting point $-10$°C.
"""

# ╔═╡ 9d79d1ed-a538-4eda-bbe2-4b3858cc4e01
html"""

"""

# ╔═╡ 5b574c5a-ce4b-4804-aa6c-58e0f7459520
md"""
In **Homework 9**, you will extend the above model to include the effect of CO₂ and determine how much CO2 would need to be added to the snowball for it to melt.
"""

# ╔═╡ 66c35971-30fb-4743-b65e-d2e527e15be7
md"""### 4) Towards realistic climate modelling

In this simple model, the preindustrial climate of $T=14$°C is so warm that there is no ice anywhere on the planet. Indeed, the only two valid stable climates are one with *no ice* or one with *ice everywhere*.

So how did Earth's preindustrial climate, which was relatively stable for thousands of years, have substantial ice caps at the poles?

"""

# ╔═╡ b15da4ce-8fdf-4d75-9060-0cdf648c5355
md"""
**The "Aquaplanet", a simple three-dimensional ocean-world climate model**

An "Aquaplanet" is a three-dimensional global climate simulation of a hypothetical planet covered entirely by a single global Ocean. While this is of course very different from Earth, where 27% of the planet is covered in land, the "Aquaplanet" exhibits many of the same characteristics as Earth and is much more realistic than our zero-dimensional climate model above.

The video below shows that the Aquaplanet simulation exhibits a third equilibrium state, with a *mostly-liquid ocean but ice caps at the poles*, in addition to the two we found in our zero-dimensional model.

In **Homework 10** (2020 class), you will build a simple two-dimensional version of the aqua-planet and explore its stability.
"""

# ╔═╡ e3e0b838-19a9-4f9b-9026-0b546d3c416b
html"""

"""

# ╔═╡ c2db5c70-b854-48b5-ba20-b64e2c7cf680
md"""
## Appendix
"""

# ╔═╡ 8eee86b6-f8e4-4016-a203-8e5e2a8b44e4
md"#### Pluto Magic"

# ╔═╡ 764a1d24-cfc2-4e28-90bb-fba264dc8956
function ingredients(path::String)
	# this is from the Julia source code (evalfile in base/loading.jl)
	# but with the modification that it returns the module instead of the last object
	name = Symbol(basename(path))
	m = Module(name)
	Core.eval(m,
        Expr(:toplevel,
             :(eval(x) = $(Expr(:core, :eval))($name, x)),
             :(include(x) = $(Expr(:top, :include))($name, x)),
             :(include(mapexpr::Function, x) = $(Expr(:top, :include))(mapexpr, $name, x)),
             :(include($path))))
	m
end

# ╔═╡ 3e84f2e2-f2d3-41ce-903f-5596096c55ba
md"""
#### Bifurcation diagram helper functions
"""

# ╔═╡ 6a96738f-4864-48b8-b345-570c895633cd
function restart_ebm!(ebm)
	ebm.T = [ebm.T[end]]
	ebm.t = [ebm.t[1]]
end

# ╔═╡ 0b654bc2-9f01-42e5-a21d-e047ccbe43e2
begin
	ebm = Model.EBM(Tneo, 0., 5., Model.CO2_const)
	ebm.S = Sneo

	ntraj = 10
	Ttraj = repeat([NaN], ntraj)
	Straj = repeat([NaN], ntraj)

	md"**Data structures for storing trajectories of recent climates**"
end

# ╔═╡ 7b21aa83-f4fc-48d2-81ac-c120dd421d73
begin
	S
	restart_ebm!(ebm)
	ebm.S = S
	Model.run!(ebm, 500)

	insert!(Straj, 1, copy(S))
	pop!(Straj)

	insert!(Ttraj, 1, copy(ebm.T[end]))
	pop!(Ttraj)
end;

# ╔═╡ 80caaf52-062b-43c5-93c9-3bbdc2d61cec
function plot_trajectory!(p, x, y; lw=8)
	n = size(x,1)
	plot!(x, y, color="black", alpha=collect(1/n:1/n:1.), linewidth=collect(0.:lw/n:lw), label=nothing)
	plot!((x[end], y[end]), color="black", marker=:c, markersize=lw/2*1.2, label=nothing, markerstrokecolor=nothing, markerstrokewidth=0.)
	return p
end;

# ╔═╡ ff6b1d24-7870-4151-8928-213de7942c7d
begin
	Smin = 1200
	Smax = 1800
	Smax_limited = 1650
	Svec = Smin:1.:Smax
	Svec = vcat(Svec, reverse(Svec))
	Tvec = zeros(size(Svec))

	local T_restart = -100.
	for (i, S) = enumerate(Svec)
		ebm = Model.EBM(T_restart, 0., 5., Model.CO2_const);
		ebm.S = S
		Model.run!(ebm, 400.)
		T_restart = ebm.T[end]
		Tvec[i] = deepcopy(T_restart)
	end

	md"**Data structures for storing warm & cool branch climates**"
end

# ╔═╡ ff687478-30a2-4d6e-9a62-859c61db6a34
begin
	T_unstable_branch(S, A, B, αi, α0) = (
		(A/B-S/(4B)*(1-0.5(αi+α0))) /
		(1+S*(αi-α0)/(80B))
	)
	S_unstable = Smin:2:Smax
	T_unstable = T_unstable_branch.(S_unstable, Model.A, Model.B, 0.5, 0.3)
	T_unstable[.~(0.3 .< α.(T_unstable) .< 0.5)] .= NaN
	md"**Unstable branch solution**"
end

# ╔═╡ b2e4513a-76aa-4a8c-8ad4-490533789c74
function add_labels!(p)
	plot!(p, xlabel="year", ylabel="temperature [°C]", legend=:bottomright, xlims=(-5,205), ylims=(-60, 30.))

	plot!(p, [-5, 200], [-60, -60], fillrange=[-10., -10.], fillalpha=0.3, c=:lightblue, label=nothing)
	annotate!(120, -20, text("completely frozen", 10, :darkblue))

	plot!(p, [-5, 200], [10, 10], fillrange=[30., 30.], fillalpha=0.09, c=:red, lw=0., label=nothing)
	annotate!(p, 120, 25, text("no ice", 10, :darkred))
end

# ╔═╡ f5361e48-e613-4fed-8da3-5576feea776d
begin
	Sneo = Model.S*0.93
	Tneo = -48.
	md"**Initial conditions**"
end

# ╔═╡ Cell order:
# ╟─d84dbb08-b1a4-11eb-2d3b-0711fddd1347
# ╠═a0b3813e-adab-11eb-2983-616cf2bb6f5e
# ╠═ef8d6690-720d-4772-a41f-b260d306b5b2
# ╟─26e1879d-ab57-452a-a09f-49493e65b774
# ╟─fd12468f-de16-47cc-8210-9266ca9548c2
# ╟─30969341-9079-4732-bf55-d6bba2c2c16c
# ╟─5bf3fe83-d096-4df1-8476-0a6500b01868
# ╠═af7f36d9-adca-48b8-95bb-ac620e6f1b4f
# ╠═790add0f-c83f-4824-92ae-53159ce58f64
# ╟─21210cfa-0366-4019-86f7-158fdd5f21ad
# ╟─978e5fc0-ddd1-4e93-a243-a95d414123b9
# ╠═6139554e-c6c9-4252-9d64-042074f68391
# ╟─e115cbbc-9d49-4fa1-8701-fa48289a0916
# ╟─bd65e7f9-ecf2-43ac-b5a8-99b03866a5c8
# ╟─991f14bb-fc4a-4505-a3be-5ced2fb148b6
# ╟─61960e99-2932-4a8f-9e87-f64a7a043489
# ╟─9acf9b8a-7fde-46c7-bf0c-c4dedc5a064b
# ╟─9bc33f22-c065-4c1d-b06f-78c70415a111
# ╠═6f61180d-6900-48fa-998d-36110e79d2dc
# ╠═a09a0064-c6ab-4912-bc9f-96ab72b8bbca
# ╟─5027e1f8-8c50-4538-949e-6c95c550016e
# ╟─465f637c-0555-498b-a881-a2f6e5714cbb
# ╠═a4686bca-90d6-4e02-961c-59f08fc37553
# ╟─f8ee2373-6af0-4d81-98fb-23bde10198ef
# ╟─a94b5160-f4bf-4ddc-9ee6-581ea20c8b90
# ╟─fd882095-6cc4-4927-967c-6b02d5b1ad95
# ╟─10693e53-3741-4388-b3b1-eba739ec01d0
# ╟─11b3fb9e-5922-4350-9424-51fba33502d4
# ╟─b71fca45-9687-4a51-8e1c-1f413e83e58d
# ╟─d993a2fc-2319-4f64-8a17-904a57593da2
# ╟─4bdc0f0c-e696-4d87-b10c-8a0da9a0ee5b
# ╟─7b7b631e-2ba3-4ed3-bad0-ec6ecb70ad49
# ╟─066743eb-c890-40b9-9f6b-9f79b7ebcbd2
# ╟─70ec6ae9-601f-4862-96cb-f251d4b5a7fd
# ╟─2bafd1a4-32a3-4787-807f-0a5132d66c28
# ╠═40b5e447-0cfb-4f35-8f95-6aa29793e5ad
# ╟─fca6c4ec-4d0c-4f97-b966-ce3a81a18710
# ╠═b1fee17b-6522-4cf0-a614-5ff8aa8f8614
# ╟─cfde8137-cfcd-46de-9c26-8abb64b6b3a9
# ╟─4351b05f-f9bf-4046-9f95-a0a56b1e8cc9
# ╠═e9942719-93cc-4203-8d37-8f91539104b1
# ╟─0c1e3051-c491-4de7-a149-ce81c53f5841
# ╟─d0d43dc3-4cda-4602-90c6-2d14a1e63871
# ╟─48431fc9-413f-4f72-8c6f-f48b42c29475
# ╟─c719ff5f-421d-4d4d-87b7-34879ab188c5
# ╠═0b6bfab9-96e9-4d7a-a486-4582d7303244
# ╟─9b6edae8-fbf0-4979-9536-c0f782ba70a7
# ╠═81d823a6-5c92-496a-96f1-a0f4762f1f05
# ╟─28924aef-9157-4490-afa5-7f232a5101f0
# ╟─5eb888e1-ba95-4aa1-97c0-784dc9d9e6d5
# ╠═18922c20-62dc-4524-8300-bbab4db828a9
# ╠═22328073-bb46-4ec2-9fdf-17c0daff5741
# ╟─9c4ac33a-ebbd-47db-9057-91624b0a2497
# ╟─2798028c-d971-45e4-9484-bdec7e8dc048
# ╟─5fc03f0f-ae12-476f-93f8-6285ee7f5fc9
# ╠═eb0307f1-c828-41f9-bc49-2aedd5b588e1
# ╠═ada47602-2cf6-4ef5-b1ba-0b1c6f81ad91
# ╠═5ebe385f-3542-4674-b419-fb54519aa945
# ╟─72e822ef-0256-4612-8b6d-654ba6f6a83f
# ╟─10815b22-b920-4cb7-a078-4d5ef457ebf2
# ╠═984f0489-7430-43b9-8926-9d7c32e7da63
# ╠═bf0dd94d-e861-478a-9bb0-d34d17405fcc
# ╟─81088439-312a-4042-bcee-c021e5e3f6b7
# ╟─9d79d1ed-a538-4eda-bbe2-4b3858cc4e01
# ╟─5b574c5a-ce4b-4804-aa6c-58e0f7459520
# ╟─66c35971-30fb-4743-b65e-d2e527e15be7
# ╟─b15da4ce-8fdf-4d75-9060-0cdf648c5355
# ╟─e3e0b838-19a9-4f9b-9026-0b546d3c416b
# ╟─c2db5c70-b854-48b5-ba20-b64e2c7cf680
# ╟─8eee86b6-f8e4-4016-a203-8e5e2a8b44e4
# ╠═764a1d24-cfc2-4e28-90bb-fba264dc8956
# ╟─3e84f2e2-f2d3-41ce-903f-5596096c55ba
# ╠═6a96738f-4864-48b8-b345-570c895633cd
# ╠═0b654bc2-9f01-42e5-a21d-e047ccbe43e2
# ╠═7b21aa83-f4fc-48d2-81ac-c120dd421d73
# ╠═80caaf52-062b-43c5-93c9-3bbdc2d61cec
# ╠═ff6b1d24-7870-4151-8928-213de7942c7d
# ╠═ff687478-30a2-4d6e-9a62-859c61db6a34
# ╠═b2e4513a-76aa-4a8c-8ad4-490533789c74
# ╠═f5361e48-e613-4fed-8da3-5576feea776d


text:
@testset "Interface" begin
    @test NMF(5) == NMF(ncomps=5)
    @test isnothing(NMF().ncomps)
end

@testset "Decomposition" begin
    data = 4 .* randn(rng, 30, 101, 101) .+ 10
    angles = sort!(90randn(rng, 30)) |> normalize_par_angles

    # get sizes correct for ncomps
    for N in [1, 3, 5]
        A, w = decompose(NMF(ncomps=N), data, angles)
        @test size(A) == (N, 101 * 101)
        @test size(w) == (30, N)
    end
    @test_throws ErrorException decompose(NMF(40), data, angles)

    # default is to use whole cube
    S = reconstruct(NMF(), data, angles)
    @test size(S) == (30, 101, 101)
end

@testset "ADI Trivial" begin
    cube = ones(10, 100, 100)
    angles = zeros(10)

    reduced_5 = NMF(5)(cube, angles)
    reduced_10 = NMF(10)(cube, angles)

    @test size(reduced_5) == size(reduced_10) == (100, 100)
    @test all(x -> isapprox(x, 0, atol = 1e-9), reduced_5)
    @test all(x -> isapprox(x, 0, atol = 1e-9), reduced_10)
end

@testset "RDI Trivial" begin
    data = randn(rng, 30, 101, 101)
    angles = sort!(90rand(rng, 30)) |> normalize_par_angles

    S = data .- reconstruct(NMF(1), data, angles, zeros(10, 101, 101))
    @test S ≈ data rtol=2e-2
end


text:
using Pkg
Pkg.activate(".")

using PowerDynamics
using OrdinaryDiffEq
using CSV
using Plots
plotly()

##

dir = @__DIR__

## set per unit MV
const base_power = 1E6 # 1MW
const base_voltage = 20E3 # 20kV
const base_current = base_power / base_voltage # 50A
const base_admittance = base_power / base_voltage^2 # 0.0025Ω^-1
const ω = 2 * π * 50.0 # 314.1593rad/s
# per unit HV
const base_voltage_HV = 110E3 # 110kV
const base_admittance_HV = base_power / base_voltage_HV^2 # 8.264462809917356e-5
const base_current_HV = base_power / base_voltage_HV

## some additional functions that might develop into new PD.jl features
#include("PDpatches.jl")
# ADN control layer (needs ControlledPowerGrid from PDpatches)
include("$dir/control.jl")
# custom types
include("$dir/DGUnit.jl")
include("$dir/OLTC.jl")
# load actual data
include("$dir/cigre_static.jl")
# failure model
include("$dir/short_circuit.jl")
# static solution
include("$dir/power_flow.jl")
# load pre-defined test cases
include("$dir/mwe.jl")

## setup system

# define events
t_step = 1.0 # step in the power flow reference
t_fault = 3.0 # onset of node short circuit
t_duration = 0.15 # duration of node short circuit

S_total = complex(24.40048903027117, 9.39761634666317) * 1e6 / base_power
P_ref(t) = t > t_step ? 0.9 * real(S_total) : real(S_total)
Q_ref(t) = t > t_step ? 0.9 * imag(S_total) : imag(S_total)
#P_ref(t) = - real(S_total)
#Q_ref(t) = - imag(S_total)

# edge data
# ode = rhs(pg)
# ode.f.graph_data.e

## find operating point

cpg, power_flow =
    testbed_CIGRE(DGUnit; P_ref = t -> P_ref(0.0), Q_ref = t -> Q_ref(0.0))

## init from load flow sol

# set all voltages to slack value
icguess = initial_guess(cpg)
# only makes sense for the full model
icguess_pf = initial_guess(cpg, power_flow.V .* exp.(power_flow.A .* 1im))

op_pf = State(cpg, icguess_pf)
check_operationpoint(cpg, op_pf)

op_ss = find_steady_state(cpg, icguess)
check_operationpoint(cpg, op_ss)

op_ss_pf = find_steady_state(cpg, icguess_pf)
check_operationpoint(cpg, op_ss_pf)

op_v = State(cpg, find_valid_initial_condition(cpg, icguess_pf))
check_operationpoint(cpg, op_v)

begin
    n = length(cpg.nodes)
    #p = scatter(real.(op_pf[2:n, :u]), imag.(op_pf[2:n, :u]), label = "pf")
    p = scatter(
        real.(op_ss_pf[2:n, :u]),
        imag.(op_ss_pf[2:n, :u]),
        label = "ss-pf",
    )
    scatter!(p, real.(op_ss[2:n, :u]), imag.(op_ss[2:n, :u]), label = "ss")
end

## check dynamics

cpg, power_flow = testbed_CIGRE(DGUnit; P_ref = P_ref, Q_ref = Q_ref)

ode = rhs(cpg)
perturb = Perturbation(nsc_node, :θ, Dec(0.))
u0 = perturb(op_ss_pf)

prob = ODEProblem(ode, u0.vec, (0.0, 2.0))

_sol = solve(prob, Rodas5(), d_discontinuities = [t_step])
sol = PowerGridSolution(_sol, cpg)

plot(sol, 2:12, :v, label = false)
# hline(op_prestep[2:12, :v], c = :gray, label = "op")
# hline!(op_poststep[2:12, :v], c = :gray, label = "op")

## node short circuit

cpg, power_flow = testbed_CIGRE(DGUnit; P_ref = P_ref, Q_ref = Q_ref)

Ron = 4. # 2.5 Ω
Znsc = complex(Ron, 0.0) * base_admittance #  fixed by Emocosy
nsc_node = 9 #(8) in CIGRE - fixed by Emocosy

nsc = NodeShortCircuitAlt(;
    node_number = nsc_node,
    Y = (1 / Znsc),
    tspan_fault = (t_fault, t_fault + t_duration),
)

@time sol = simulate(nsc, cpg, op_ss_pf.vec, (0.0, 4.0))

plot(sol, 2:12, :v, label = false, xlabel="t [s]", ylabel="v [pu]")
plot(sol, 2:12, :θ, label = false, xlabel="t [s]", ylabel="PLL [rad]")
plot(sol, 2:12, :P_g, label = false, xlabel="t [s]", ylabel="Pg [pu]")
plot(sol, 2:12, :iqp, label = false, xlabel="t [s]", ylabel="Iq+ [pu]")

Δ = [
    cpg.controller(u, nothing, t) |> first
    for (t, u) in zip(sol.dqsol.t, sol.dqsol.u)
]

plot(sol.dqsol.t, first.(Δ), label="ΔP")
plot!(sol.dqsol.t, last.(Δ), label="ΔQ")

##
fcpg = nsc(cpg)

fault_op = find_operationpoint(fcpg, icguess_pf)

# TODO idea: solve PF for voltages
# TODO better: solve constraints first and then other variables, compar DE.jl

nsc_op = State(fcpg, find_valid_initial_condition(fcpg, op.vec))

function stepwise_init(cpg, x, p = nothing, t0 = 0.0)
    rpg = rhs(cpg)
    M = rpg.mass_matrix
    algebraic_vars = [all(iszero, c) for c in eachcol(M)]
    diff_vars = .!algebraic_vars
    algebraic_eqs = [all(iszero, r) for r in eachrow(M)]
    diff_eqs = .!algebraic_eqs

    out = copy(x)

    dview = @view out[diff_vars]
    # fix differential vars y and solve g(x, y)=0
    function solve_with_fixed_diff(du, u)
        u[diff_vars] .= dview
        rpg(du, u, p, t0)
        du .= algebraic_eqs .* du
    end

    res = nlsolve(
        solve_with_fixed_diff,
        x,
        autodiff = :forward,
        show_trace = false,
        method = :trust_region,
    )
    out[algebraic_vars] .= res.zero[algebraic_vars]
    println(res.x_converged, res.f_converged, res.residual_norm)
    return out
end

rand = State(fcpg, 10 * randn(systemsize(fcpg)))
ic = State(fcpg, stepwise_init(fcpg, op.vec .+ rand.vec))
ic2 = State(fcpg, relax_alg(fcpg, op.vec))

@show op[2:12, :v]
@show fault_op[2:12, :v]
@show nsc_op[2:12, :v]

@show op[2:12, :P_int]
@show nsc_op[2:12, :P_int]

timespan = (0.0, 4.0)
problem = ODEProblem{true}(rhs(cpg), op_ss_pf.vec, timespan)
integrator = init(
    problem,
    Rodas5(autodiff = false),
    #force_dtmin = true,
    #save_everystep = false,
    tstops = [t_fault, t_fault + t_duration],
)
# step to failure
step!(integrator, t_fault, true)
check_error(integrator)
pre_failure = State(cpg, copy(integrator.u))
@show pre_failure[2:12, :v]
@show pre_failure[2:12, :iqp]
integrator.f = rhs(nsc(cpg))
DiffEqBase.initialize_dae!(integrator, BrownFullBasicInit())
valid_post_failure = State(cpg, copy(integrator.u))
@show valid_post_failure[2:12, :v]
@show valid_post_failure[2:12, :iqp]
small_step = 0.00001
step!(integrator, small_step, true)
check_error(integrator)
end_of_failure = State(cpg, copy(integrator.u))
@show end_of_failure[2:12, :v]
@show end_of_failure[2:12, :iqp]

timespan = (t_fault, t_fault+t_duration)
x2 = find_valid_initial_condition(nsc(cpg), op.vec)

dx = similar(x2)
rpg(dx, x2, nothing, t_fault)
@show dx
findall(abs.(dx) .> 0.1)
findall(isnan.(dx))

problem = ODEProblem{true}(rhs(nsc(cpg)), x2, timespan)
sol = solve(
    problem,
    ABDF2(),
    #Rodas4(autodiff = true),
    #force_dtmin = true,
)

using NLsolve
using LinearAlgebra
using LineSearches

## plots

plot(sol, 2:12, :v, label = false)
hline!(fault_op[2:12, :v], c = :gray, label = "op")
vline!(
    [t_step, t_fault, t_fault + t_duration],
    label = "fault",
    c = :black,
    ls = :dot,
)

plot(sol, 2:12, :φ, label = false)
plot!(sol, 2:12, :θ, ls = :dot, c = :black, label = "PLL")
hline!(op[:, :φ], c = :gray, label = "op")
#vline!([t_fault, t_fault + t_duration], label = "fault", c = :black)
#ylims!(-0.01, 0.01)

plot(sol, 2:12, :P_g)
plot!(sol, 2:12, :Q_g)
hline!(op[2:12, :P_g], label = "op")
hline!(op[2:12, :Q_g], label = "op")
vline!([t_fault, t_fault + t_duration], label = "fault", c = :black, ls = :dot)

plot(sol, 2:12, :P_int)
plot!(sol, 2:12, :Q_int)
hline!(op[2:12, :P_int], label = "op")
hline!(op[2:2, :Q_int], label = "op")
vline!([t_fault, t_fault + t_duration], label = "fault", c = :black, ls = :dot)

Δ = [
    cpg.controller(u, nothing, t) |> first
    for (t, u) in zip(sol.dqsol.t, sol.dqsol.u)
]

plot(sol.dqsol.t, first.(Δ))
plot!(sol.dqsol.t, last.(Δ))

plot(sol, 2:12, :iqp)

inc = (last(sol.dqsol.prob.tspan) - first(sol.dqsol.prob.tspan)) / 100.0
t_array = range(
    first(sol.dqsol.prob.tspan);
    stop = last(sol.dqsol.prob.tspan),
    step = inc,
)

# currents
node = nsc_node
# current set point
I_r = [
    complex(sol(t, node, :P_g), -sol(t, node, :Q_g)) / sol(t, node, :v)
    for t in t_array
]
plot(t_array, abs.(I_r), label = "I_r_$node")

max_stars_repo_name: luap-pik/PowerDynIC

max_issues_repo_name: luap-pik/PowerDynIC

max_forks_repo_name: luap-pik/PowerDynIC


text:
@testset "test result csv" begin
	function test_csv_results(results, filename::AbstractString, filedir::AbstractString)
		file_tmp = joinpath(filedir, filename)
		file_data = DelimitedFiles.readdlm(file_tmp, ',')

		num_samples = size(results["inputs"]["pd"], 1)
		input_size = sum([size(x, 2) for x in values(results["inputs"])])
		output_size = sum([size(x, 2) for x in values(results["outputs"])])
		dual_size = sum([size(x, 2) for x in values(results["duals"]) if x isa Array]) - 1 # UAS not a var

		@test (size(file_data, 1) - 1) == num_samples
		@test size(file_data, 2) == (input_size + output_size + dual_size)
		rm(file_tmp)
		@test true
	end

    @testset "5-bus case file csv check" begin
		save_order = vcat(OPFLearn.DEFAULT_INPUTS, OPFLearn.DEFAULT_OUTPUTS, OPFLearn.DEFAULT_DUALS)
		data_dir = "data"
		net_file = "pglib_opf_case5_pjm.m"
		K = 10
		results = create_samples(net_file, K, net_path=data_dir)
        file_name = "TEMP_"*net_file*"_"*"dataset"
		save_results_csv(results, file_name, save_order=save_order, dir=data_dir)
		test_csv_results(results, file_name*".csv", data_dir)
	end

	@testset "5-bus case file csv check with specified vars" begin
		input_vars = ["pd", "qd"]
		output_vars = ["p_gen", "q_gen"]
		dual_vars = ["v_max"]
		save_order = vcat(input_vars, output_vars, dual_vars)
		data_dir = "data"
		net_file = "pglib_opf_case5_pjm.m"
		K = 10
		results = create_samples(net_file, K, net_path=data_dir, input_vars=input_vars,
							     output_vars=output_vars, dual_vars=dual_vars)
        file_name = "TEMP_"*net_file*"_"*"dataset"
		save_results_csv(results, file_name, save_order=save_order, dir=data_dir)
		test_csv_results(results, file_name*".csv", data_dir)
	end
end


text:
using Test, YaoBase, LuxurySparse
using BitBasis: invorder, reorder

@testset "batch normalize" begin
    s = rand(3, 4)
    batch_normalize!(s, 1)
    for i = 1:4
        @test sum(s[:, i]) ≈ 1
    end

    s = rand(3, 4)
    ss = batch_normalize(s, 1)
    for i = 1:4
        @test sum(s[:, i]) != 1
        @test sum(ss[:, i]) ≈ 1
    end
end

@testset "hilbertkron" begin
    A, B, C, D = [randn(2, 2) for i = 1:4]
    II = IMatrix(2)
    ⊗ = kron
    @test hilbertkron(4, [A, B], [3, 1]) ≈ II ⊗ A ⊗ II ⊗ B
    @test hilbertkron(4, [A ⊗ B, C], [3, 1]) ≈ A ⊗ B ⊗ II ⊗ C
    @test hilbertkron(4, [A ⊗ B], [1]) ≈ II ⊗ II ⊗ A ⊗ B
    @test hilbertkron(4, [A ⊗ B, C ⊗ D], [1, 3]) ≈ C ⊗ D ⊗ A ⊗ B

    U = randn(2, 2)
    U2 = randn(4, 4)
    m = U2 ⊗ II ⊗ U ⊗ II
    @test m == hilbertkron(5, [U, U2], [2, 4])
end

@testset "batched kron" begin
    A, B, C =
        rand(ComplexF64, 4, 4, 3), rand(ComplexF64, 4, 4, 3), rand(ComplexF64, 4, 4, 3)
    D = batched_kron(A, B, C)

    tD = zeros(ComplexF64, 64, 64, 3)
    for k = 1:3
        tD[:, :, k] = kron(A[:, :, k], B[:, :, k], C[:, :, k])
    end

    @test tD ≈ D

    B2 = reshape(transpose(reshape(permutedims(B, (3, 1, 2)), 3, 16)), 4, 4, 3)
    @test B2 isa Base.ReshapedArray
    @test Array(B2) ≈ B
    D2 = batched_kron(A, B2, C)
    @test tD ≈ D2
end

@testset "random matrices" begin
    mat = rand_unitary(8)
    @test isunitary(mat)
    mat = rand_hermitian(8)
    @test ishermitian(mat)

    @test ishermitian(sprand_hermitian(8, 0.5))
    @test isunitary(sprand_unitary(8, 0.5))
end

@testset "rotmat" begin
    theta = 0.5
    @test rot_mat(ComplexF64, Const.X, theta) ≈ ComplexF64[
        cos(theta / 2) -im*sin(theta / 2)
        -im*sin(theta / 2) cos(theta / 2)
    ]
    @test rot_mat(ComplexF64, Const.X, theta) |> eltype == ComplexF64
    @test rot_mat(ComplexF32, Const.X, theta) |> eltype == ComplexF32
end

using LuxurySparse: pmrand
@testset "reorder" begin
    ⊗ = kron
    PA = pmrand(2)
    PB = pmrand(2)
    PC = pmrand(2)
    @test reorder(PC ⊗ PB ⊗ PA, [3, 1, 2]) ≈ PB ⊗ PA ⊗ PC
    @test invorder(PC ⊗ PB ⊗ PA) ≈ PA ⊗ PB ⊗ PC
end

@testset "logdi" begin
    @test logdi(9, 3) == 2
    @test_throws ArgumentError logdi(9, 5)
end


text:
#    Problem : GROUP A
#    *********
#    The BRYBND problem by Oren.
#
#    Source: problem 31 in
#    J.J. More', B.S. Garbow and K.E. Hillstrom,
#    "Testing Unconstrained Optimization Software",
#    ACM Transactions on Mathematical Software, vol. 7(1), pp. 17-41, 1981.
#
#    See also Buckley#73 (p. 41) and Toint#18
#
#    BRYBND.SIF classification SUR2-AN-V-0
#
#    Number of variables is variable

f = (x) -> begin
	@warn "Does not agree with CUTEst.jl"
	mL = 5
	mU = 1
    fx = 0.0
	n = length(x)
    for i in firstindex(x):lastindex(x)-1
		bnd = 0.0
		0 > i-mL ? lo=0 : lo=i-mL
		n-1 < i+mU ? u=n-1 : u=i+mU

		for j in lo+1:i
        	bnd += x[j]*(1+x[j])
		end
        for j in i+2:u+1
			bnd += x[j]*(1+x[j])
		end

		group = (2.0 + 5.0x[i]^2)x[i] + 1 - bnd
		fx += group^2
    end
    return fx
end

g = (g, x) -> begin
	@warn "Does not agree with CUTEst.jl"
	mL = 5
	mU = 1
	n = length(x)
    for i in firstindex(x):lastindex(x)-1
		bnd = 0.0
		0 > i-mL ? lo=0 : lo=i-mL
		n-1 < i+mU ? u=n-1 : u=i+mU

		for j in lo+1:i
        	bnd += x[j]*(1+x[j])
		end
        for j in i+2:u+1
			bnd += x[j]*(1+x[j])
		end

		group = (2.0 + 5.0x[i]^2)x[i] + 1 - bnd
		g[i] += 2.0(2.0 +15.0x[i]^2)group

		for j in lo+1:i
            g[j] -= 2.0(1+ 2.0x[j])group
		end

        for j = i+2:u+1
            g[j] -= 2.0(1+ 2.0*x[j])group
		end
    end
    return g
end

fg! = (g, x) -> begin
    @warn "Does not agree with CUTEst.jl"
	mL = 5
	mU = 1
	n = length(x)
    for i in firstindex(x):lastindex(x)-1
		bnd = 0.0
		0 > i-mL ? lo=0 : lo=i-mL
		n-1 < i+mU ? u=n-1 : u=i+mU

		for j in lo+1:i
        	bnd += x[j]*(1+x[j])
		end
        for j in i+2:u+1
			bnd += x[j]*(1+x[j])
		end

		group = (2.0 + 5.0x[i]^2)x[i] + 1 - bnd
		fx += group^2
		g[i] += 2.0(2.0 +15.0x[i]^2)group

		for j in lo+1:i
            g[j] -= 2.0(1+ 2.0x[j])group
		end

        for j = i+2:u+1
            g[j] -= 2.0(1+ 2.0*x[j])group
		end
    end
    return fx, g
end

init = (n::Int=5000) -> begin
	n < 2 && @warn("BRYBND: number of variables must be ≥ 2")
	n = max(n, 2)

    return n, -1.0*ones(n)
end

@warn "TODO: Issue in CUTEst? See https://github.com/JuliaSmoothOptimizers/OptimizationProblems.jl/blob/main/src/brybnd.jl"
# TestSet["BRYBND"] = UncProgram("BRYBND", f, g!, fg!, init)


text:
# joint problem of avoiding one car and one pedestrian

# State type
struct PedCarMDPState
    crash::Bool
    ego::VehicleState
    ped::VehicleState
    car::VehicleState
    route::SVector{2, LaneTag}
end

# copy b to a
function Base.copyto!(a::PedCarMDPState, b::PedCarMDPState)
    a.crash = b.crash
    a.ego = b.ego
    a.ped = b.ped
    a.car = b.car
    a.route = b.route
end

function Base.hash(s::PedCarMDPState, h::UInt64 = zero(UInt64))
    return hash(s.crash, hash(s.ego, hash(s.ped, hash(s.car, hash(s.route, h)))))
end

function Base.:(==)(a::PedCarMDPState, b::PedCarMDPState)
    return a.crash == b.crash && a.ego == b.ego && a.ped == b.ped && a.car == b.car && a.route == b.route
end

# Action type
const PedCarMDPAction = UrbanAction

@with_kw mutable struct PedCarMDP <: MDP{PedCarMDPState, PedCarMDPAction}
    env::UrbanEnv = UrbanEnv(params=UrbanParams(nlanes_main=1,
                     crosswalk_pos =  [VecSE2(6, 0., pi/2), VecSE2(-6, 0., pi/2), VecSE2(0., -5., 0.)],
                     crosswalk_length =  [10.0, 10., 10.0],
                     crosswalk_width = [4.0, 4.0, 3.1],
                     stop_line = 22.0))
    ego_type::VehicleDef = VehicleDef()
    car_type::VehicleDef = VehicleDef()
    car_models::Dict{SVector{2, LaneTag}, DriverModel} = get_car_models(env, get_stop_model)
    ped_type::VehicleDef = VehicleDef(AgentClass.PEDESTRIAN, 1.0, 1.0)
    ped_model::DriverModel = ConstantPedestrian()
    max_acc::Float64 = 2.0
    pos_res::Float64 = 3.0
    vel_res::Float64 = 2.0
    vel_ped_res::Float64 = 1.0
    ego_start::Float64 = env.params.stop_line - ego_type.length/2
    ego_goal::LaneTag = LaneTag(2,1)
    off_grid::VecSE2 = VecSE2(UrbanEnv().params.x_min+VehicleDef().length/2, env.params.y_min+VehicleDef().width/2, 0)
    ΔT::Float64 = 0.5
    car_birth::Float64 = 0.3
    ped_birth::Float64 = 0.3
    collision_cost::Float64 = -1.0
    action_cost::Float64 = 0.
    goal_reward::Float64 = 1.
    γ::Float64 = 0.95
    _ped_grid::Dict{LaneTag, RectangleGrid{3}} = init_ped_grid(env, pos_res, vel_ped_res)
    _car_grid::Dict{LaneTag, RectangleGrid{2}} = init_car_grid(env, pos_res, vel_res)
    _l_grid::Dict{LaneTag, RectangleGrid{1}} = init_l_grid(env, pos_res)
    _v_grid::RectangleGrid{1} = init_v_grid(env, vel_res)
    _collision_checker::Dict{Tuple{Vararg{VehicleState, 3}}, Bool} = Dict{Tuple{Vararg{VehicleState, 3}}, Bool}()
    _car_proj::Dict{Tuple{Float64, Float64, LaneTag}, VehicleState} = car_projection(env, pos_res, vel_res)
    _ped_proj::Dict{Tuple{Float64, Float64, Float64, LaneTag}, VehicleState} = ped_projection(env, pos_res, vel_ped_res)
end

function POMDPs.reward(mdp::PedCarMDP, s::PedCarMDPState, a::PedCarMDPAction, sp::PedCarMDPState)
    r = mdp.action_cost
    ego = sp.ego
    if sp.crash
        r += mdp.collision_cost
    end
    if ego.posF.s >= get_end(mdp.env.roadway[mdp.ego_goal]) &&
       get_lane(mdp.env.roadway, ego).tag == mdp.ego_goal
        r += mdp.goal_reward
    end
    return r
end

function POMDPs.isterminal(mdp::PedCarMDP, s::PedCarMDPState)
    if s.crash
        return true
    elseif s.ego.posF.s >= get_end(mdp.env.roadway[mdp.ego_goal]) &&
       get_lane(mdp.env.roadway, s.ego).tag == mdp.ego_goal
       return true
   end
   return false
end

## Helpers

function POMDPs.discount(mdp::PedCarMDP)
    return mdp.γ
end

POMDPs.actions(mdp::PedCarMDP) = [PedCarMDPAction(-4.0), PedCarMDPAction(-2.0), PedCarMDPAction(0.0), PedCarMDPAction(2.0)]

function POMDPs.actionindex(mdp::PedCarMDP, action::PedCarMDPAction)
    if action.acc == -4.0
        return 1
    elseif action.acc == -2.0
        return 2
    elseif action.acc == 0.
        return 3
    else
        return 4
    end
end

function POMDPs.convert_s(::Type{V}, s::PedCarMDPState, mdp::PedCarMDP) where V<:AbstractArray
    n_routes = 4
    n_features = 4
    z = zeros(n_features*3 + n_routes + 1)
    z[1] = s.ego.posG.x / abs(mdp.env.params.x_min)
    z[2] = s.ego.posG.y / abs(mdp.env.params.y_min)
    z[3] = s.ego.posG.θ / π
    z[4] = s.ego.v / abs(mdp.env.params.speed_limit)
    z[5] = s.ped.posG.x / abs(mdp.env.params.x_min)
    z[6] = s.ped.posG.y / abs(mdp.env.params.y_min)
    z[7] = s.ped.posG.θ / π
    z[8] = s.ped.v / abs(mdp.env.params.speed_limit)
    z[9] = s.car.posG.x / abs(mdp.env.params.x_min)
    z[10] = s.car.posG.y / abs(mdp.env.params.y_min)
    z[11] = s.car.posG.θ / π
    z[12] = s.car.v / abs(mdp.env.params.speed_limit)
    # one hot encoding for the route
    routes = get_car_routes(mdp.env)
    for (i, r) in enumerate(routes)
        if r[1] == s.route[1] && r[end] == s.route[end]
            z[12+i] = 1.
        end
    end
    z[17] = float(s.crash)
    return z
end

function POMDPs.convert_s(::Type{PedCarMDPState}, z::V, mdp::PedCarMDP) where V<:AbstractArray{Float64}
    n_routes = 4
    n_features = 4
    @assert length(z) == n_features*3 + n_routes + 1
    ego_x = z[1]*abs(mdp.env.params.x_min)
    ego_y = z[2]*abs(mdp.env.params.y_min)
    ego_θ = z[3]*π
    ego_v = z[4]*abs(mdp.env.params.speed_limit)
    ego = VehicleState(VecSE2(ego_x, ego_y, ego_θ), mdp.env.roadway, ego_v)
    ped_x = z[5]*abs(mdp.env.params.x_min)
    ped_y = z[6]*abs(mdp.env.params.y_min)
    ped_θ = z[7]*π
    ped_v = z[8]*abs(mdp.env.params.speed_limit)
    ped = VehicleState(VecSE2(ped_x, ped_y, ped_θ), mdp.env.ped_roadway, ped_v)
    car_x = z[9]*abs(mdp.env.params.x_min)
    car_y = z[10]*abs(mdp.env.params.y_min)
    car_θ = z[11]*π
    car_v = z[12]*abs(mdp.env.params.speed_limit)
    car = VehicleState(VecSE2(car_x, car_y, car_θ), mdp.env.roadway, car_v)
    # one hot encoding for the route
    routes = get_car_routes(mdp.env)
    route = SVector{2, LaneTag}(LaneTag(0, 0), LaneTag(0, 0))
    for (i, r) in enumerate(routes)
        if z[12+i] == 1.
            route = SVector{2, LaneTag}(r[1], r[end])
        end
    end
    collision = Bool(z[17])
    return PedCarMDPState(collision, ego, ped, car, route)
end

function get_off_the_grid(mdp::PedCarMDP)
    return VehicleState(mdp.off_grid, Frenet(mdp.env.roadway[LaneTag(5,1)], 25.1, -26.5, pi/2), 0.)
    # return VehicleState(mdp.off_grid, mdp.env.roadway, 0.)
end

# create labels for model checking
function labeling(mdp::PedCarMDP)
    labels = Dict{PedCarMDPState, Vector{String}}()
    for s in ordered_states(mdp)
        if crash(mdp, s)
            labels[s] = ["crash"]
        elseif s.ego.posF.s >= get_end(mdp.env.roadway[mdp.ego_goal]) &&
            get_lane(mdp.env.roadway, s.ego).tag == mdp.ego_goal
            labels[s] = ["goal"]
        end
    end
    return labels
end

function crash(mdp::PedCarMDP, s::PedCarMDPState)
    return crash(mdp, s.ego, s.car, s.ped)
end

function crash(mdp::PedCarMDP, ego::VehicleState, car::VehicleState, ped::VehicleState)
    return collision_checker(ego, car, mdp.ego_type, mdp.car_type) || collision_checker(ego, ped, mdp.ego_type, mdp.ped_type)
end

max_stars_repo_name: MaximeBouton/AutomotivePOMDPs

max_issues_repo_name: MaximeBouton/AutomotivePOMDPs

max_forks_repo_name: MaximeBouton/AutomotivePOMDPs


text:
# Analyzing PowerModels' raw parser, it was concluded that b_to & b_fr data was present in DSHL section
function _handle_b(pwf_data::Dict, f_bus::Int, t_bus::Int, susceptance::Float64, circuit::Int, dict_dshl)
    b_fr, b_to = susceptance / 200, susceptance / 200
    if haskey(pwf_data, "DSHL") && haskey(dict_dshl, (f_bus, t_bus)) && haskey(dict_dshl[(f_bus, t_bus)], circuit)
        group = dict_dshl[(f_bus, t_bus)][circuit]
        b_fr += group["SHUNT FROM"] !== nothing && group["STATUS FROM"] == " L" ? group["SHUNT FROM"] / 100 : 0
        b_to += group["SHUNT TO"] !== nothing && group["STATUS TO"] == " L" ? group["SHUNT TO"] / 100 : 0
    end
    return b_fr, b_to
end

function _create_dict_dshl(data::Dict)
    dshl_dict = Dict{Tuple, Any}()

    for (k,v) in data
        sub_data = Dict{String, Any}()
        sub_data["SHUNT FROM"] = v["SHUNT FROM"]
        sub_data["SHUNT TO"] = v["SHUNT TO"]
        sub_data["STATUS FROM"] = v["STATUS FROM"]
        sub_data["STATUS TO"] = v["STATUS TO"]
        current_value = get(dshl_dict, (v["FROM BUS"], v["TO BUS"]), Dict{Int, Any}())
        current_value[v["CIRCUIT"]] = sub_data
        dshl_dict[(v["FROM BUS"], v["TO BUS"])] = current_value
    end
    return dshl_dict
end

function _pwf2pm_branch!(pm_data::Dict, pwf_data::Dict, branch::Dict; add_control_data::Bool=false)
    sub_data = Dict{String,Any}()

    sub_data["f_bus"] = pop!(branch, "FROM BUS")
    sub_data["t_bus"] = pop!(branch, "TO BUS")
    sub_data["br_r"] = pop!(branch, "RESISTANCE") / 100
    sub_data["br_x"] = pop!(branch, "REACTANCE") / 100

    dshl_dict = haskey(pwf_data, "DSHL") ? _create_dict_dshl(pwf_data["DSHL"]) : nothing
    b = _handle_b(pwf_data, sub_data["f_bus"], sub_data["t_bus"], branch["SHUNT SUSCEPTANCE"], branch["CIRCUIT"], dshl_dict)
    sub_data["g_fr"] = 0.0
    sub_data["b_fr"] = b[1]
    sub_data["g_to"] = 0.0
    sub_data["b_to"] = b[2]

    sub_data["tap"] = pop!(branch, "TAP")
    sub_data["shift"] = -pop!(branch, "LAG")
    sub_data["angmin"] = -60.0 # PowerModels.jl standard
    sub_data["angmax"] = 60.0 # PowerModels.jl standard
    sub_data["transformer"] = false

    if branch["STATUS"] == branch["OPENING FROM BUS"] == branch["OPENING TO BUS"] == 'L'
        sub_data["br_status"] = 1
    else
        sub_data["br_status"] = 0
    end

    if add_control_data
        sub_data["control_data"] = Dict{String,Any}()
        sub_data["control_data"]["tapmin"] = sub_data["tap"]
        sub_data["control_data"]["tapmax"] = sub_data["tap"]
        sub_data["control_data"]["circuit"] = branch["CIRCUIT"]

        # Transformer control_data fields are created and set to nothing
        sub_data["control_data"]["control_type"] = nothing
        sub_data["control_data"]["constraint_type"] = nothing
        sub_data["control_data"]["controlled_bus"] = nothing
        sub_data["control_data"]["vmsp"] = nothing
        sub_data["control_data"]["vmmin"] = nothing
        sub_data["control_data"]["vmmax"] = nothing
        sub_data["control_data"]["shift_control_variable"] = nothing
        sub_data["control_data"]["shiftmin"] = nothing
        sub_data["control_data"]["shiftmax"] = nothing
        sub_data["control_data"]["valsp"] = nothing
        sub_data["control_data"]["control"] = nothing
    end

    sub_data["source_id"] = ["branch", sub_data["f_bus"], sub_data["t_bus"], "01"]
    sub_data["index"] = length(pm_data["branch"]) + 1

    sub_data["rate_a"] = pop!(branch, "NORMAL CAPACITY")
    sub_data["rate_b"] = pop!(branch, "EMERGENCY CAPACITY")
    sub_data["rate_c"] = pop!(branch, "EQUIPAMENT CAPACITY")

    if sub_data["rate_a"] >= 9999
        delete!(sub_data, "rate_a")
    end
    if sub_data["rate_b"] >= 9999
        delete!(sub_data, "rate_b")
    end
    if sub_data["rate_c"] >= 9999
        delete!(sub_data, "rate_c")
    end

    idx = string(sub_data["index"])
    pm_data["branch"][idx] = sub_data

end

function _pwf2pm_DCSC_branch!(pm_data::Dict, pwf_data::Dict, branch::Dict; add_control_data::Bool=false)
    sub_data = Dict{String,Any}()

    sub_data["f_bus"] = pop!(branch, "FROM BUS")
    sub_data["t_bus"] = pop!(branch, "TO BUS")
    sub_data["br_r"] = 0
    sub_data["br_x"] = pop!(branch, "INITIAL VALUE") / 100

    sub_data["g_fr"] = 0.0
    sub_data["b_fr"] = 0.0
    sub_data["g_to"] = 0.0
    sub_data["b_to"] = 0.0

    sub_data["tap"] = 1.0
    sub_data["shift"] = 0
    sub_data["angmin"] = -60.0 # PowerModels.jl standard
    sub_data["angmax"] = 60.0 # PowerModels.jl standard
    sub_data["transformer"] = false

    if branch["STATUS"] == 'L'
        sub_data["br_status"] = 1
    else
        sub_data["br_status"] = 0
    end

    if add_control_data
        sub_data["control_data"] = Dict{String,Any}()
        sub_data["control_data"]["tapmin"] = 1.0
        sub_data["control_data"]["tapmax"] = 1.0
        sub_data["control_data"]["circuit"] = branch["CIRCUIT"]

        # Transformer control_data fields are created and set to nothing
        sub_data["control_data"]["control_type"] = nothing
        sub_data["control_data"]["constraint_type"] = nothing
        sub_data["control_data"]["controlled_bus"] = nothing
        sub_data["control_data"]["vmsp"] = nothing
        sub_data["control_data"]["vmmin"] = nothing
        sub_data["control_data"]["vmmax"] = nothing
        sub_data["control_data"]["shift_control_variable"] = nothing
        sub_data["control_data"]["shiftmin"] = nothing
        sub_data["control_data"]["shiftmax"] = nothing
        sub_data["control_data"]["valsp"] = nothing
        sub_data["control_data"]["control"] = nothing
    end

    sub_data["source_id"] = ["branch", sub_data["f_bus"], sub_data["t_bus"], "01"]
    sub_data["index"] = length(pm_data["branch"]) + 1

    sub_data["rate_a"] = pop!(branch, "NORMAL CAPACITY")
    sub_data["rate_b"] = pop!(branch, "EMERGENCY CAPACITY")
    sub_data["rate_c"] = pop!(branch, "EQUIPAMENT CAPACITY")

    if sub_data["rate_a"] >= 9999
        delete!(sub_data, "rate_a")
    end
    if sub_data["rate_b"] >= 9999
        delete!(sub_data, "rate_b")
    end
    if sub_data["rate_c"] >= 9999
        delete!(sub_data, "rate_c")
    end

    rep = findall(x -> x["f_bus"] == sub_data["f_bus"] && x["t_bus"] == sub_data["t_bus"] && x["circuit"] == sub_data["circuit"], pm_data["branch"])
    if length(rep) > 0
        @warn "Branch from $(sub_data["f_bus"]) to $(sub_data["t_bus"]) in circuit $(sub_data["circuit"]) is duplicated"
    end
    idx = string(sub_data["index"])
    pm_data["branch"][idx] = sub_data

end

function _pwf2pm_branch!(pm_data::Dict, pwf_data::Dict; add_control_data::Bool=false)
    pm_data["branch"] = Dict{String, Any}()
    if haskey(pwf_data, "DLIN")
        for (i,branch) in pwf_data["DLIN"]
            if !branch["TRANSFORMER"]
                _pwf2pm_branch!(pm_data, pwf_data, branch, add_control_data = add_control_data)
            end
        end
    end
    if haskey(pwf_data, "DCSC")
        for (i,csc) in pwf_data["DCSC"]
            _pwf2pm_DCSC_branch!(pm_data, pwf_data, csc, add_control_data = add_control_data)
        end
    end
end