text:
module WikiText

export
    WikiText2, WikiText103, WikiText2Raw, WikiText103Raw,
    trainfile, validfile, testfile,
    isunk

using DataDeps

abstract type WikiTextCorpus end

abstract type WikiTextWordCorpus <: WikiTextCorpus end
abstract type WikiTextCharCorpus <: WikiTextCorpus end

"""
    WikiText2v1

WikiText2v1 corpus for word-level language modeling.
See https://blog.einstein.ai/the-wikitext-long-term-dependency-language-modeling-dataset/.
"""
struct WikiText2v1 <: WikiTextWordCorpus end

"""
    WikiText103v1

WikiText103v1 corpus for word-level language modeling.
See https://blog.einstein.ai/the-wikitext-long-term-dependency-language-modeling-dataset/.
"""
struct WikiText103v1 <: WikiTextWordCorpus end

"""
    WikiText2RawV1

WikiText2RawV1 corpus for character-level language modeling.
See https://blog.einstein.ai/the-wikitext-long-term-dependency-language-modeling-dataset/.
"""
struct WikiText2RawV1 <: WikiTextCharCorpus end

"""
    WikiText103RawV1

WikiText103RawV1 corpus for character-level language modeling.
See https://blog.einstein.ai/the-wikitext-long-term-dependency-language-modeling-dataset/.
"""
struct WikiText103RawV1 <: WikiTextCharCorpus end

const EOS = "<eos>"

const UNK = "<unk>"
isunk(token) = token == unk

const WikiText2 = WikiText2v1
const WikiText103 = WikiText103v1
const WikiText2Raw = WikiText2RawV1
const WikiText103Raw = WikiText103RawV1

# API:
trainfile(corpus) = filename(corpus, :train)
validfile(corpus) = filename(corpus, :valid)
testfile(corpus)  = filename(corpus, :test)

function filename(corpus::WikiTextCorpus, set)
    @assert set in [:train, :valid, :test]
    filename = "wiki.$set.$(suffix(corpus))"
    return joinpath(corpusdir(corpus), filename)
end
filename(corpus::Type{<:WikiTextCorpus}, set) = filename(corpus(), set)

suffix(corpus::WikiTextWordCorpus)   = "tokens"
suffix(corpus::WikiTextCharCorpus)   = "raw"

corpusdir(corpups::WikiText2v1)      = datadep"WikiText-2-v1"
corpusdir(corpups::WikiText103v1)    = datadep"WikiText-103-v1"
corpusdir(corpups::WikiText2RawV1)   = datadep"WikiText-2-raw-v1"
corpusdir(corpups::WikiText103RawV1) = datadep"WikiText-103-raw-v1"

function __init__()
    moveup(x) = mv(x, joinpath("..", x))
    register(DataDep("WikiText-2-v1",
                     """
                     Dataset: WikiText-2 word level language modeling dataset
                     Author: Stephen Merity
                     License: CC-SA 3.0
                     Website: https://einstein.ai/research/the-wikitext-long-term-dependency-language-modeling-dataset
                     Size: 12MB (unzipped)

                     Wikitext-2 word level long term dependency language modeling dataset.

                     """,
                     "https://s3.amazonaws.com/research.metamind.io/wikitext/wikitext-2-v1.zip",
                     "92675f1d63015c1c8b51f1656a52d5bdbc33aafa60cc47a218a66e7ee817488c",
                     post_fetch_method = function (zip)
                                             unpack(zip)
                                             pth(x) = joinpath("..", x)
                                             cd("wikitext-2") do
                                                 moveup("wiki.train.tokens")
                                                 moveup("wiki.test.tokens")
                                                 moveup("wiki.valid.tokens")
                                             end
                                             rm("wikitext-2")
                                         end))

    register(DataDep("WikiText-103-v1",
                     """
                     Dataset: WikiText-103 word level language modeling dataset
                     Author: Stephen Merity
                     License: CC-SA 3.0
                     Website: https://einstein.ai/research/the-wikitext-long-term-dependency-language-modeling-dataset
                     Size: 516 MB

                     Wikitext-103 word level long term dependency language modeling dataset.

                     """,
                     "https://s3.amazonaws.com/research.metamind.io/wikitext/wikitext-103-v1.zip",
                     "242ba0f20b329cfdf1ccc61e9e9e5b59becf189db7f7a81cd2a0e2fc31539590",
                     post_fetch_method = function (zip)
                         unpack(zip)
                         cd("wikitext-103") do
                             moveup("wiki.train.tokens")
                             moveup("wiki.valid.tokens")
                             moveup("wiki.test.tokens")
                         end
                     end))

    register(DataDep("WikiText-2-raw-v1",
                     """
                     Dataset: WikiText-2 raw character level language modeling dataset
                     Author: Stephen Merity
                     License: CC-SA 3.0
                     Website: https://einstein.ai/research/the-wikitext-long-term-dependency-language-modeling-dataset
                     Size: 12 MB (unzipped)

                     Wikitext-2 raw character level long term dependency language modeling dataset.

                     """,
                     "https://s3.amazonaws.com/research.metamind.io/wikitext/wikitext-2-raw-v1.zip",
                     "ef7edb566e3e2b2d31b29c1fdb0c89a4cc683597484c3dc2517919c615435a11",
                     post_fetch_method = function (zip)
                         unpack(zip)
                         cd("wikitext-2-raw") do
                             moveup("wiki.train.raw")
                             moveup("wiki.valid.raw")
                             moveup("wiki.test.raw")
                         end
                         rm("wikitext-2-raw")
                     end))

    register(DataDep("WikiText-103-raw-v1",
                     """
                     Dataset: WikiText-2 word level language modeling dataset
                     Author: Stephen Merity
                     License: CC-SA 3.0
                     Website: https://einstein.ai/research/the-wikitext-long-term-dependency-language-modeling-dataset
                     Size: 518MB (unzipped)

                     Wikitext-2 word level long term dependency language modeling dataset.

                     """,
                     "https://s3.amazonaws.com/research.metamind.io/wikitext/wikitext-103-raw-v1.zip",
                     "91c00ae287f0d699e18605c84afc9e45c192bc6b7797ff8837e5474655a33794",
                     post_fetch_method = function (zip)
                         unpack(zip)
                         cd("wikitext-103-raw") do
                             moveup("wiki.train.raw")
                             moveup("wiki.valid.raw")
                             moveup("wiki.test.raw")
                         end
                         rm("wikitext-103-raw")
                     end))
end

end # module


text:
function YakProblems(Rot=UnitQuaternion{Float64,CayleyMap}; scenario=:barrellroll, use_rot=Rot<:UnitQuaternion,
        costfun=:Quadratic, normcon=false)
    model = Dynamics.YakPlane(Rot, use_rot=use_rot)
    rsize = size(model)[1] - 9

    opts = SolverOptions(
        cost_tolerance_intermediate = 1e-1,
        penalty_scaling = 1000.,
        penalty_initial = 0.01
    )

    # Discretization
    tf = 1.25
    N = 101
    dt = tf/(N-1)
    # N = Int(1.25/dt+1)
    # dt = 0.025

    if scenario == :barrellroll
        ey = @SVector [0,1,0.]

        # Initial and final condition
        p0 = MRP(0.997156, 0., 0.075366) # initial orientation
        pf = MRP(0., -0.0366076, 0.) # final orientation (upside down)

        costfun == :QuatLQR ? sq = 0 : sq = 1
        if Rot <: UnitQuaternion
            rm_quat = @SVector [1,2,3,4,5,6,8,9,10,11,12,13]
        else
            rm_quat = @SVector [1,2,3,4,5,6,7,8,9,10,11,12]
        end

        x0 = Dynamics.build_state(model, [-3,0,1.5], p0, [5,0,0], [0,0,0])
        utrim  = @SVector  [41.6666, 106, 74.6519, 106]
        xf = Dynamics.build_state(model, [3,0,1.5], pf, [5,0,0], [0,0,0])

        # Objective
        Qf_diag = Dynamics.fill_state(model, 100, 500*sq, 100, 100.)
        Q_diag = Dynamics.build_state(model, [0.1,0.1,0.1],
            (@SVector fill(0.1*sq,rsize)), fill(0.1,3), fill(0.1,3))
        Qf = Diagonal(Qf_diag)
        Q = Diagonal(Q_diag)
        R = Diagonal(@SVector fill(1e-3,4))
        if costfun == :Quadratic
            costfun = LQRCost(Q, R, xf, utrim)
            costterm = LQRCost(Qf, R, xf, utrim)
        elseif costfun == :QuatLQR
            costfun = QuatLQRCost(Q, R, xf, utrim; w=0.1)
            costterm = QuatLQRCost(Qf, R, xf, utrim; w=200.0)
        elseif costfun == :ErrorQuad
            Q = Diagonal(Q_diag[rm_quat])
            Qf = Diagonal(Qf_diag[rm_quat])
            costfun = ErrorQuadratic(model, Q, R, xf, utrim)
            costterm = ErrorQuadratic(model, Qf, R, xf, utrim)
        end
        obj = Objective(costfun, costterm, N)

        # Constraints
        conSet = ConstraintSet(n,m,N)
        goal = GoalConstraint(xf)
        add_constraint!(conSet,goal,N:N)

    else
        throw(ArgumentError("$scenario isn't a known scenario"))
    end

    # Initialization
    U0 = [copy(utrim) for k = 1:N-1]

    # Build problem
    prob = Problem(model, obj, xf, tf, x0=x0, constraints=conSet)
    initial_controls!(prob, U0)
    prob, opts
end

ext: jl


text:
using PaddedViews

#TODO use parallel processing to speed up the unrolling process

#TODO build a reroll function to extract W from K

### unroll conv1d
@doc raw"""
    unroll(cLayer::Conv1D, AiS::Tuple, param::Symbol=:W)

unroll the `param` of `Conv1D` into 2D matrix

# Arguments

- `cLayer` := the layer of the paramters to unroll

- `AiS` := the `padded` input to determinde the size and shape of the output of `unroll`

- `param` := `Conv1D` parameter to be `unroll`ed

# Return

- `K` := 2D `Matrix` of the `param`
"""
function unroll(cLayer::Conv1D, AiS::Tuple, param::Symbol=:W)
    prevLayer = cLayer.prevLayer
    n_Hi, ci, m = AiS
    f_H = cLayer.f
    s_H = cLayer.s
    c = cLayer.channels
    n_H = (n_Hi - f_H) ÷ s_H + 1
    W = eval(:($cLayer.$param))
    T = eltype(W)
    B = cLayer.B
    K = nothing
    for ch=1:size(W)[end]
        k2d = nothing
        w3d = W[:,:,ch]
        k = nothing
        for i=range(1, step=s_H, length=n_H)
            w2d = PaddedView(0, w3d, (n_Hi,ci), (i,1))
            for w=1:size(w2d)[2]
                    if w==1
                        k2d = w2d[:,w]
                    else
                        k2d = vcat(k2d, w2d[:,w])
                    end #if h==1
            end #for w=1:size(w3d)[3]

            if i==1
                k = k2d
            else
                k = hcat(k, k2d)
            end #if i==1
        end #for i=1:n_W

        if ch == 1
            K = k
        else
            K = hcat(K, k)
        end #if ch==1
    end #for ch=1:size(W)[end]

    return transpose(K)

end #function unroll(cLayer::Conv2D)

###unroll conv2d
@doc raw"""
    unroll(cLayer::Conv2D, AiS::Tuple, param::Symbol=:W)

unroll the `param` of `Conv1D` into 2D matrix

# Arguments

- `cLayer` := the layer of the paramters to unroll

- `AiS` := the `padded` input to determinde the size and shape of the output of `unroll`

- `param` := `Conv1D` parameter to be `unroll`ed

# Return

- `K` := 2D `Matrix` of the `param`
"""
function unroll(cLayer::Conv2D, AiS::Tuple, param::Symbol=:W)
    prevLayer = cLayer.prevLayer
    n_Hi, n_Wi, ci, m = AiS
    HWi = n_Hi*n_Wi
    f_H, f_W = cLayer.f
    fHW = f_H * f_W
    s_H, s_W = cLayer.s
    c = cLayer.channels
    n_H = (n_Hi - f_H) ÷ s_H + 1
    n_W = (n_Wi - f_W) ÷ s_W + 1
    W = eval(:($cLayer.$param))
    T = eltype(W)
    B = cLayer.B
    K = nothing
    for ch=1:size(W)[end]
        w3da = W[:,:,:,ch]
        k2d = nothing
        k3d = nothing
        k4d = nothing
        k = nothing
        for i=range(1, step=s_H, length=n_H)
            w3db = PaddedView(0, w3da, (n_Hi, f_W,ci), (i,1,1))
            for j=range(1, step=s_W, length=n_W)
                w3d = PaddedView(0, w3db, (n_Hi, n_Wi,ci), (1,j,1))
                for w=1:size(w3d)[3]
                    for h=1:size(w3d)[1]
                        if h==1
                            k2d = w3d[h,:,w]
                        else
                            k2d = vcat(k2d, w3d[h,:,w])
                        end #if h==1
                    end #for h=1:size(w)[1]
                    if w==1
                        k3d = k2d
                    else
                        k3d = vcat(k3d, k2d)
                    end #if w==1
                end #for w=1:size(w3d)[3]

                if j==1
                    k4d = k3d
                else
                    k4d = hcat(k4d, k3d)
                end
            end #for j=1:n_H

            if i==1
                k = k4d
            else
                k = hcat(k, k4d)
            end #if i==1
        end #for i=1:n_W

        if ch==1
            K = k
        else
            K = hcat(K, k)
        end #if ch==1
    end #for ch=1:size(W)[end]

    return transpose(K)

end #function unroll(cLayer::Conv2D)

###unroll conv3d

@doc raw"""
    unroll(cLayer::Conv3D, AiS::Tuple, param::Symbol=:W)

unroll the `param` of `Conv3D` into 2D matrix

# Arguments

- `cLayer` := the layer of the paramters to unroll

- `AiS` := the `padded` input to determinde the size and shape of the output of `unroll`

- `param` := `Conv1D` parameter to be `unroll`ed

# Return

- `K` := 2D `Matrix` of the `param`
"""
function unroll(cLayer::Conv3D, AiS::Tuple, param::Symbol=:W)
    prevLayer = cLayer.prevLayer
    n_Hi, n_Wi, n_Di, ci, m = AiS
    HWDi = n_Hi*n_Wi*n_Di
    f_H, f_W, f_D = cLayer.f
    fHW = f_H * f_W * f_D
    s_H, s_W, s_D = cLayer.s
    c = cLayer.channels
    n_H = (n_Hi - f_H) ÷ s_H + 1
    n_W = (n_Wi - f_W) ÷ s_W + 1
    n_D = (n_Di - f_D) ÷ s_D + 1
    W = eval(:($cLayer.$param))
    T = eltype(W)
    B = cLayer.B
    K = nothing
    for ch=1:size(W)[end]
        w4da = W[:,:,:,:,ch]
        k2d = nothing
        k3d = nothing
        k4d = nothing
        k5d = nothing
        k = nothing
        k1 = nothing
        for l=range(1, step=s_D, length=n_D)
            w4db = PaddedView(0, w4da, (f_H, f_W, n_Di, ci), (1,1,l,1))
            for i=range(1, step=s_H, length=n_H)
                w4dc = PaddedView(0, w4db, (n_Hi, f_W, n_Di, ci), (i,1,1,1))
                for j=range(1, step=s_W, length=n_W)
                    w4d = PaddedView(0, w4dc, (n_Hi, n_Wi, n_Di, ci), (1,j,1,1))
                    for w=1:size(w4d)[4]
                        for d=1:size(w4d)[3]
                            for h=1:size(w4d)[1]
                                if h==1
                                    k2d = w4d[h,:,d,w]
                                else
                                    k2d = vcat(k2d, w4d[h,:,d,w])
                                end #if h==1
                            end #for h=1:size(w)[1]
                            if d==1
                                k3d = deepcopy(k2d)
                            else
                                k3d = vcat(k3d, k2d)
                            end #if w==1
                        end #for w=1:size(w3d)[3]
                        if w==1
                            k4d = deepcopy(k3d)
                        else
                            k4d = vcat(k4d, k3d)
                        end #if d==1
                    end #for d=1:size(w4d)[4]
                    if j==1
                        k5d = deepcopy(k4d)
                    else
                        k5d = hcat(k5d, k4d)
                    end
                end #for j=1:n_H

                if i==1
                    k = deepcopy(k5d)
                else
                    k = hcat(k, k5d)
                end #if i==1
            end #for i=1:n_W
            if l==1
                k1 = deepcopy(k)
            else
                k1 = hcat(k1, k)
            end #if l==1
        end #for l=range(1, step=s_D, length=n_D)
        if ch==1
            K = k1
        else
            K = hcat(K, k1)
        end #if ch==1
    end #for ch=1:size(W)[end]

    return transpose(K)

end #function unroll(cLayer::Conv2D)

export unroll

function unrollcol(cLayer::Conv2D, AiS::Tuple, param::Symbol=:W)
    prevLayer = cLayer.prevLayer
    n_Hi, n_Wi, ci, m = AiS
    HWi = n_Hi*n_Wi
    f_H, f_W = cLayer.f
    fHW = f_H * f_W
    s_H, s_W = cLayer.s
    c = cLayer.channels
    n_H = (n_Hi - f_H) ÷ s_H + 1
    n_W = (n_Wi - f_W) ÷ s_W + 1
    W = eval(:($cLayer.$param))
    T = eltype(W)
    B = cLayer.B
    K = nothing
    for ch=1:size(W)[end]
        w3da = W[:,:,:,ch]
        k2d = nothing
        k3d = nothing
        k4d = nothing
        k = nothing
        for i=range(1, step=s_W, length=n_W)
            w3db = PaddedView(0, w3da, (f_H, n_Wi,ci), (1,i,1))
            for j=range(1, step=s_H, length=n_H)
                w3d = PaddedView(0, w3db, (n_Hi, n_Wi,ci), (j,1,1))
                for cj=1:size(w3d)[3]
                    for w=1:size(w3d)[2]
                        if w==1
                            k2d = w3d[:,w,cj]
                        else
                            k2d = vcat(k2d, w3d[:,w,cj])
                        end #if h==1
                    end #for h=1:size(w)[1]
                    if cj==1
                        k3d = k2d
                        k2d = nothing
                    else
                        k3d = vcat(k3d, k2d)
                        k2d = nothing
                    end #if w==1
                end #for w=1:size(w3d)[3]

                if j==1
                    k4d = k3d
                    k3d = nothing
                else
                    k4d = hcat(k4d, k3d)
                    k3d = nothing
                end
            end #for j=1:n_H

            if i==1
                k = k4d
                k4d = nothing
            else
                k = hcat(k, k4d)
                k4d = nothing
            end #if i==1
        end #for i=1:n_W

        if ch==1
            K = k
            k = nothing
            Base.GC.gc()
        else
            K = hcat(K, k)
            k2d = nothing
            Base.GC.gc()
        end #if ch==1
    end #for ch=1:size(W)[end]

    return transpose(K)

end #function unroll(cLayer::Conv2D)


text:
using LightGraphs, SimpleWeightedGraphs

function bfs_depth(input_data, start_position, go_condition, stop_condition=(x->false))
  status = zeros(Int32, size(input_data)...)
  queue::Array{NamedTuple{(:position, :depth),Tuple{CartesianIndex, Int32}},1} = []
  if typeof(start_position) <: Tuple
    start_position = CartesianIndex(start_position)
  end
  push!(queue, (position=start_position, depth=0))
  while length(queue) > 0
    (position, depth) = popfirst!(queue)
    #println(position)
    if status[position] != 0
      continue
    end
    status[position] = depth
    if ! stop_condition(input_data[position])
      for movement in [CartesianIndex(-1,0), CartesianIndex(1,0), CartesianIndex(0,-1), CartesianIndex(0,1)]
        new_pos = position + movement
        if go_condition(input_data[new_pos])
          push!(queue, (position=new_pos, depth=depth + 1))
        end
      end
    end
  end
  return status
end

function generate_distance_graph(distances)
  println("distances: $distances")
  tokens = keys(distances)
  my_tokens = Set()
  for token in tokens
    push!(my_tokens, token[1])
    push!(my_tokens, token[2])
  end
  #g = path_graph(maximum(my_tokens))
  mg = SimpleWeightedDiGraph(maximum(my_tokens))
  for token in tokens
    println("Adding: $(token[1]), $(token[2]), $(distances[token])")
    add_edge!(mg, token[1], token[2], distances[token])
  end
  return mg
end

function dijkstra(distances, start_token)
  mg = generate_distance_graph(distances)

  ds = dijkstra_shortest_paths(mg, start_token)
  println(ds.dists)

  return ds.dists

end

max_stars_repo_name: buotex/aoc

max_issues_repo_name: buotex/aoc

max_forks_repo_name: buotex/aoc


text:
#--------------------------------------------------------------------
# DNSS.jl
# Soham 03-2022
# Plotting routines for 1D, 2D and Arrays.
# TODO: Plotting still royally sucks in Julia. I wouldn't
# waste time with this.
#--------------------------------------------------------------------

using PyPlot, LaTeXStrings, LinearAlgebra

"""
    Plot fields in 1D space
"""
function PyPlot. plot(u::Field{S}, path) where {S}
    rcParams = PyPlot.PyDict(PyPlot.matplotlib."rcParams")
    rcParams["font.size"] = 8
    rcParams["text.usetex"] = true
    rcParams["axes.titlesize"] = 8.0
    rcParams["axes.labelsize"] = 10.0
    rcParams["xtick.labelsize"] = 8.0
    rcParams["ytick.labelsize"] = 8.0
    rcParams["legend.fontsize"] =  8.0
    rcParams["figure.figsize"] = (10, 10)
    rcParams["figure.dpi"] = 300
    rcParams["savefig.dpi"] = 300
    w, h = plt[:figaspect](10)
    figure(figsize=(w,h))

    x = Field(u.space, x->x)
    plot(x.value, u.value, "-")

    # xlabel(L"$v$")
    # ylabel(L"$a(u_0, v)$")
    tight_layout()
    savefig("$path")
    close()
end

"""
    Plot fields in 2D space
"""
function PyPlot. contourf(w::Field{ProductSpace{S1, S2}}, levels; wmin=nothing, wmax=nothing) where {S1, S2}
    u  = Field(w.space.S1, u->u)
    v  = Field(w.space.S2, v->v)
    contourf(v.value, u.value, w.value, levels; vmin=wmin, vmax=wmax, cmap="Blues")
end

"""
    Plot fields in 2D array space
"""
function PyPlot. contourf(AoF::T, levels::Int, path::String) where {T<:Array{Field, 2}}

    rcParams = PyPlot.PyDict(PyPlot.matplotlib."rcParams")
    rcParams["font.size"] = 8
    rcParams["text.usetex"] = true
    rcParams["axes.titlesize"] = 8.0
    rcParams["axes.labelsize"] = 10.0
    rcParams["xtick.labelsize"] = 8.0
    rcParams["ytick.labelsize"] = 8.0
    rcParams["legend.fontsize"] =  8.0
    rcParams["figure.figsize"] = (10, 10)
    rcParams["figure.dpi"] = 300
    rcParams["savefig.dpi"] = 300
    # w, h = plt[:figaspect](10)
    # figure(figsize=(w,h))

    globalmin = minimum(AoF)
    globalmax = maximum(AoF)
    globallevels = collect(range(globalmin, stop=globalmax, length=levels))
    for index in CartesianIndices(AoF)
        if eltype(AoF[index].value) <: Number
            contourf(AoF[index], globallevels, wmin=globalmin, wmax=globalmax)
        end
    end
    colorbar()
    xlabel(L"$v$")
    ylabel(L"$u$")
    tight_layout()
    savefig("$path")
    close("all")
end

"""
Plot p convergence
"""
function plotpconv(n_::Vector, l_::Vector, path::String)
    rcParams = PyPlot.PyDict(PyPlot.matplotlib."rcParams")
    rcParams["font.size"] = 8
    rcParams["text.usetex"] = true
    rcParams["axes.titlesize"] = 8.0
    rcParams["axes.labelsize"] = 10.0
    rcParams["xtick.labelsize"] = 8.0
    rcParams["ytick.labelsize"] = 8.0
    rcParams["legend.fontsize"] =  8.0
    rcParams["figure.figsize"] = (10, 10)
    rcParams["figure.dpi"] = 300
    rcParams["savefig.dpi"] = 300

    w, h = plt[:figaspect](10)
    figure(figsize=(w,h))

    semilogy(n_, l_, "o--", linewidth=0.5,  markersize=1.0)
    xlabel(L"$N_p$")
    xticks(collect(12:2:24))
    ylabel(L"$\|\mathcal{C}\|_2$")
    tight_layout()
    savefig("$path")
    close()
end

"""
Plot h convergence
"""
function plothconv(n_::Vector, l_::Vector, path::String)
    rcParams = PyPlot.PyDict(PyPlot.matplotlib."rcParams")
    rcParams["font.size"] = 8
    rcParams["text.usetex"] = true
    rcParams["axes.titlesize"] = 8.0
    rcParams["axes.labelsize"] = 10.0
    rcParams["xtick.labelsize"] = 8.0
    rcParams["ytick.labelsize"] = 8.0
    rcParams["legend.fontsize"] =  8.0
    rcParams["figure.figsize"] = (10, 10)
    rcParams["figure.dpi"] = 300
    rcParams["savefig.dpi"] = 300
    ll_ = similar(l_)
    w, h = plt[:figaspect](10)
    figure(figsize=(w,h))

    for index in 2:length(l_)
        ll_[index] = l_[index - 1] / l_[index]
    end

    plot(log.(2, n_[2:end]), ll_[2:end], "o--", linewidth=0.5, markersize=1.0)
    xlabel(L"$\log_2(N_h)$")
    ylabel(L"$\Delta_{N_h} / \Delta_{N_{h+1}} $")
    tight_layout()
    savefig("$path")
    close()
end

function plotmodes(u::Field{ProductSpace{ChebyshevGL{Tag1, N1, T},
                                         ChebyshevGL{Tag2, N2, T}}}, path::String) where {Tag1, Tag2, N1, N2, T}
    rcParams = PyPlot.PyDict(PyPlot.matplotlib."rcParams")
    rcParams["font.size"] = 8
    rcParams["text.usetex"] = true
    rcParams["axes.titlesize"] = 8.0
    rcParams["axes.labelsize"] = 10.0
    rcParams["xtick.labelsize"] = 8.0
    rcParams["ytick.labelsize"] = 8.0
    rcParams["legend.fontsize"] =  8.0
    rcParams["figure.figsize"] = (10, 10)
    rcParams["figure.dpi"] = 300
    rcParams["savefig.dpi"] = 300
    w, h = plt[:figaspect](10)
    figure(figsize=(w,h))

    # Compute the maximum along the secondary diagonals
    ull = basistransform(u)
    B   = rotr90(ull.value,3)
    (lu, lv) = size(B)
    @assert lu == lv
    ul = [maximum(abs.(B[diagind(B, ind)])) for ind in -lu+1:lv-1]

    semilogy(ul, "o--", linewidth=0.5, markersize=1.0)
    xlabel(L"$l_u + l_v$")
    ylabel(L"$\mathrm{max}(|c_{uv}|)$")
    tight_layout()
    savefig("$path")
    close()
end


text:
module Markers
# Start Of Frame markers, non-differential, Huffman coding
const SOF0 = 0xFFC0
const SOF1 = 0xFFC1
const SOF2 = 0xFFC2
const SOF3 = 0xFFC3

# Start Of Frame markers, differential, Huffman coding
const SOF5 = 0xFFC5
const SOF6 = 0xFFC6
const SOF7 = 0xFFC7

#Start Of Frame markers, non-differential, arithmetic coding
const JPG = 0xFFC8
const SOF9 = 0xFFC9
const SOF10 = 0xFFCA
const SOF11 = 0xFFCB

# Start Of Frame markers, differential, arithmetic coding
const SOF13 = 0xFFCD
const SOF14 = 0xFFCE
const SOF15 = 0xFFCF

# Huffman table specification
const DHT = 0xFFC4

# Arithmetic coding conditioning specification
const DAC = 0xFFCC

# Restart interval termination
const RST0 = 0xFFD0
const RST1 = 0xFFD1
const RST2 = 0xFFD2
const RST3 = 0xFFD3
const RST4 = 0xFFD4
const RST5 = 0xFFD5
const RST6 = 0xFFD6
const RST7 = 0xFFD7

# Other markers
const SOI = 0xFFD8
const EOI = 0xFFD9
const SOS = 0xFFDA
const DQT = 0xFFDB
const DNL = 0xFFDC
const DRI = 0xFFDD
const DHP = 0xFFDE
const EXP = 0xFFDF
const APP0 = 0xFFE0
const APP1 = 0xFFE1
const APP2 = 0xFFE2
const APP3 = 0xFFE3
const APP4 = 0xFFE4
const APP5 = 0xFFE5
const APP6 = 0xFFE6
const APP7 = 0xFFE7
const APP8 = 0xFFE8
const APP9 = 0xFFE9
const APP10 = 0xFFEA
const APP11 = 0xFFEB
const APP12 = 0xFFEC
const APP13 = 0xFFED
const APP14 = 0xFFEE
const APP15 = 0xFFEF
const JPG0 = 0xFFF0
const JPG1 = 0xFFF1
const JPG2 = 0xFFF2
const JPG3 = 0xFFF3
const JPG4 = 0xFFF4
const JPG5 = 0xFFF5
const JPG6 = 0xFFF6
const JPG7 = 0xFFF7
const JPG8 = 0xFFF8
const JPG9 = 0xFFF9
const JPG10 = 0xFFFA
const JPG11 = 0xFFFB
const JPG12 = 0xFFFC
const JPG13 = 0xFFFD
const COM = 0xFFFE

end # module Markers

const markers = Dict{UInt16,Symbol}(map(s -> (getfield(Markers, s), s),
                                        filter(s -> getfield(Markers, s) isa UInt16,
                                               names(Markers, all=true))))

struct Marker{symbol} end

function Marker(bits::UInt16)
    0xFF00 < bits < 0xFFFF || error("invalid marker: 0x", hex(bits))
    Marker{get(markers, bits, bits)}()
end

Base.sizeof(::Type{<:Marker}) = 2
Base.sizeof(::Marker) = 2

const SOI = Marker{:SOI}()
const EOI = Marker{:EOI}()

Base.show(io::IO, m::Marker) = print(io, typeof(m), "()")

Base.write(io::IO, m::Marker) = write(io, hton(reinterpret(UInt16, m)))

function Base.reinterpret(::Type{UInt16}, ::Marker{symbol}) where symbol
    symbol isa Symbol ? getfield(Markers, symbol) : symbol
end

==(x::Marker, y::Marker) = reinterpret(UInt16, x) === reinterpret(UInt16, y)


text:
using StatsBase, Combinatorics, Plots ; pyplot()

matchExists1(n) = 1 - prod([k/365 for k in 365:-1:365-n+1])
matchExists2(n) = 1- factorial(365,365-big(n))/365^big(n)

function bdEvent(n)
    birthdays = rand(1:365,n)
    dayCounts = counts(birthdays, 1:365)
    return maximum(dayCounts) > 1
end

probEst(n) = sum([bdEvent(n) for _ in 1:N])/N

xGrid = 1:50
analyticSolution1 = [matchExists1(n) for n in xGrid]
analyticSolution2 = [matchExists2(n) for n in xGrid]
println("Maximum error: $(maximum(abs.(analyticSolution1 - analyticSolution2)))")

N = 10^3
mcEstimates = [probEst(n) for n in xGrid]

plot(xGrid, analyticSolution1, c=:blue, label="Analytic solution")
scatter!(xGrid, mcEstimates, c=:red, ms=6, msw=0, shape=:xcross,
	label="MC estimate", xlims=(0,50), ylims=(0, 1),
	xlabel="Number of people in room",
	ylabel="Probability of birthday match",
	legend=:topleft)

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
@printf("\nTesting L-BFGS routines\n")

function test_lbfgs()

Q = [32.0000 13.1000 -4.9000 -3.0000  6.0000  2.2000  2.6000  3.4000 -1.9000 -7.5000;
 13.1000 18.3000 -5.3000 -9.5000  3.0000  2.1000  3.9000  3.0000 -3.6000 -4.4000;
 -4.9000 -5.3000  7.7000  2.1000 -0.4000 -3.4000 -0.8000 -3.0000  5.3000  5.5000;
 -3.0000 -9.5000  2.1000 20.1000  1.1000  0.8000  -12.4000 -2.5000  5.5000  2.1000;
  6.0000  3.0000 -0.4000  1.1000  3.8000  0.6000  0.5000  0.9000 -0.4000 -2.0000;
  2.2000  2.1000 -3.4000  0.8000  0.6000  7.8000  2.9000 -1.3000 -4.3000 -5.1000;
  2.6000  3.9000 -0.8000  -12.4000  0.5000  2.9000 14.5000  1.7000 -4.9000  1.2000;
  3.4000  3.0000 -3.0000 -2.5000  0.9000 -1.3000  1.7000  6.6000 -0.8000  2.7000;
 -1.9000 -3.6000  5.3000  5.5000 -0.4000 -4.3000 -4.9000 -0.8000  7.9000  5.7000;
 -7.5000 -4.4000  5.5000  2.1000 -2.0000 -5.1000  1.2000  2.7000  5.7000 16.1000]

q = [2.9000, 0.8000, 1.3000, -1.1000, -0.5000, -0.3000, 1.0000, -0.3000, 0.7000, -2.1000]

xs = [1.0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09;
      0.09 1.0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08;
      0.08 0.09 1.0 0.01 0.02 0.03 0.04 0.05 0.06 0.07;
      0.07 0.08 0.09 1.0 0.01 0.02 0.03 0.04 0.05 0.06;
      0.06 0.07 0.08 0.09 1.0 0.01 0.02 0.03 0.04 0.05]'

dirs_ref = [
  -3.476000000000000e+01    -6.861170733797231e-01    -1.621334774299757e-01    -2.008976150849174e-01    -2.317011191832649e-01;
  -1.367700000000000e+01    -1.661270665201917e+00     2.870743130038511e-01     2.237224648542354e-01     2.980080835636926e-02;
   2.961000000000000e+00     2.217225828759783e-01    -5.485761164147891e-01     4.811889625788801e-02    -1.267017945785352e-01;
   3.756000000000000e+00     5.615134140894827e-01     9.992734938824949e-02    -6.855884193567087e-01     4.328230970765587e-02;
  -5.618000000000001e+00    -1.922426760799171e-01    -1.332550298134261e-02    -2.729265954345345e-02    -2.437461022925742e-01;
  -1.571000000000000e+00    -8.961101045874649e-02     5.326252573648003e-02     3.651730112313705e-02     1.349716200511426e-02;
  -4.121000000000000e+00    -3.044802963260585e-01    -6.299408068289100e-02     6.325330777317102e-02    -7.155992987801297e-04;
  -3.709000000000000e+00    -1.996235459345302e-01     1.525398352758626e-02     2.871281112230844e-02    -3.513449694839536e-03;
   4.010000000000000e-01     1.267604425710271e-01    -7.776943954825602e-02    -1.285590864125103e-01    -5.603489763638488e-02;
   7.639999999999999e+00     3.360845247013288e-01    -2.335884953507600e-02    -3.204963735369062e-03     5.612114259243499e-02
]

dirs = zeros(10, 5) # matrix of directions (to be filled in)

mem = 3;
x = zeros(10)

H = LBFGS(x, mem)
dir  = zeros(10)
println(H)

HH = LBFGS(ArrayPartition(x, x), mem)
dirdir = ArrayPartition(zeros(10), zeros(10))
println(HH)

x_old = [];
grad_old = [];

for i = 1:5

    x = xs[:,i]
    grad = Q*x + q

    if i > 1
      xx           = ArrayPartition(x, copy(x))
      xx_old       = ArrayPartition(x_old, copy(x_old))
      gradgrad     = ArrayPartition(grad, copy(grad))
      gradgrad_old = ArrayPartition(grad_old, copy(grad_old))
        @time update!(H, x, x_old, grad, grad_old)
        @time update!(HH,xx, xx_old, gradgrad, gradgrad_old)
    end

    dir_ref = dirs_ref[:,i]

    gradm = -grad
    @time mul!(dir, H, gradm)
    @test norm(dir-dir_ref, Inf)/(1+norm(dir_ref, Inf)) <= 1e-15

    gradm2 = ArrayPartition(-grad,copy(-grad))
    @time mul!(dirdir, HH, gradm2)
    @test norm(dirdir.x[1]-dir_ref, Inf)/(1+norm(dir_ref, Inf)) <= 1e-15
    @test norm(dirdir.x[2]-dir_ref, Inf)/(1+norm(dir_ref, Inf)) <= 1e-15

    x_old = x;
    grad_old = grad;

end

#testing reset

@test ones(size(H,1)) != H*ones(size(H,1))
@test ArrayPartition(ones.(size(HH,1))) != HH*ArrayPartition(ones.(size(HH,1))...)

AbstractOperators.reset!(H)
AbstractOperators.reset!(HH)

@test ones(size(H,1)) == H*ones(size(H,1))
@test ArrayPartition(ones.(size(HH,1))) == HH*ArrayPartition(ones.(size(HH,1))...)

end

test_lbfgs()


text:
# solve.jl
# Solve ODEs in three-stage chemostat experiment
# Callbacks determine when to dilute and which species to refresh
# Uses CVODE_BDF solver from Sundials
# July 2020, Nadanai Laohakunakorn (nadanai.laohakunakorn@ed.ac.uk)

using Sundials

function marktime!(integrator)
end

function condition(u,t,integrator)
    t%INTERVAL_DIL
end

###############################################
######### 1. Resource dependent model #########
###############################################

function solvemodel(grads,u0,params,TMAX,INTERVAL_DIL,tsave,DIL_FRAC,NSPECIES,SWITCHTIMES,FLAG)
	tspan = (0.0,TMAX);
	prob = ODEProblem(grads,u0,tspan,params);

	# Callbacks
	saved_values = SavedValues(Float64, Array{Float64});
	savecb = SavingCallback((u,t,integrator)->integrator(t,Val{1}),  # Save derivatives as well
	    saved_values,saveat=tsave);
    if FLAG=="SR"
	   contcb = ContinuousCallback(condition,diluteSR!;save_positions=(false,false))
    elseif FLAG=="PC"
        contcb = ContinuousCallback(condition,dilutePC!;save_positions=(false,false))
    elseif FLAG=="NC"
        contcb = ContinuousCallback(condition,diluteNC!;save_positions=(false,false))
    end
	periodcb = PeriodicCallback(marktime!,INTERVAL_DIL;save_positions=(false,false)) # hack to mark time points for dilution
	cb = CallbackSet(periodcb,contcb,savecb)

	# Solve
	sol = solve(prob,CVODE_BDF(),abstol=1e-10,reltol=1e-10, callback=cb,saveat=tsave);
	return(sol,saved_values)
end

###############################################
######## 2. Resource independent model ########
###############################################

function solvemodel_RI(grads,u0,params,TMAX,INTERVAL_DIL,tsave,DIL_FRAC,NSPECIES,SWITCHTIMES,FLAG)
    tspan = (0.0,TMAX);
    prob = ODEProblem(grads,u0,tspan,params);

    # Callbacks
    saved_values = SavedValues(Float64, Array{Float64});
    savecb = SavingCallback((u,t,integrator)->integrator(t,Val{1}),  # Save derivatives as well
        saved_values,saveat=tsave);
    if FLAG=="SR"
       contcb = ContinuousCallback(condition,diluteSR_RI!;save_positions=(false,false))
    elseif FLAG=="PC"
        contcb = ContinuousCallback(condition,dilutePC_RI!;save_positions=(false,false))
    elseif FLAG=="NC"
        contcb = ContinuousCallback(condition,diluteNC_RI!;save_positions=(false,false))
    end
    periodcb = PeriodicCallback(marktime!,INTERVAL_DIL;save_positions=(false,false)) # hack to mark time points for dilution
    cb = CallbackSet(periodcb,contcb,savecb)

    # Solve
    sol = solve(prob,CVODE_BDF(),abstol=1e-10,reltol=1e-10, callback=cb,saveat=tsave);
    return(sol,saved_values)
end

###############################################
## 3. Resource independent model  no loading ##
###############################################
function solvemodel_RInl(grads,u0,params,TMAX,INTERVAL_DIL,tsave,DIL_FRAC,NSPECIES,SWITCHTIMES,FLAG)
    tspan = (0.0,TMAX);
    prob = ODEProblem(grads,u0,tspan,params);

    # Callbacks
    saved_values = SavedValues(Float64, Array{Float64});
    savecb = SavingCallback((u,t,integrator)->integrator(t,Val{1}),  # Save derivatives as well
        saved_values,saveat=tsave);
    if FLAG=="SR"
       contcb = ContinuousCallback(condition,diluteSR_RInl!;save_positions=(false,false))
    elseif FLAG=="PC"
        contcb = ContinuousCallback(condition,dilutePC_RInl!;save_positions=(false,false))
    elseif FLAG=="NC"
        contcb = ContinuousCallback(condition,diluteNC_RInl!;save_positions=(false,false))
    end
    periodcb = PeriodicCallback(marktime!,INTERVAL_DIL;save_positions=(false,false)) # hack to mark time points for dilution
    cb = CallbackSet(periodcb,contcb,savecb)

    # Solve
    sol = solve(prob,CVODE_BDF(),abstol=1e-10,reltol=1e-10, callback=cb,saveat=tsave);
    return(sol,saved_values)
end

###############################################
######## 4. Coupled model TL saturation #######
###############################################

function solvemodel_1RI(grads,u0,params,TMAX,INTERVAL_DIL,tsave,DIL_FRAC,NSPECIES,SWITCHTIMES,FLAG)
    tspan = (0.0,TMAX);
    prob = ODEProblem(grads,u0,tspan,params);

    # Callbacks
    saved_values = SavedValues(Float64, Array{Float64});
    savecb = SavingCallback((u,t,integrator)->integrator(t,Val{1}),  # Save derivatives as well
        saved_values,saveat=tsave);
    if FLAG=="SR"
       contcb = ContinuousCallback(condition,diluteSR_1RI!;save_positions=(false,false))
    elseif FLAG=="PC"
        contcb = ContinuousCallback(condition,dilutePC_1RI!;save_positions=(false,false))
    elseif FLAG=="NC"
        contcb = ContinuousCallback(condition,diluteNC_1RI!;save_positions=(false,false))
    end
    periodcb = PeriodicCallback(marktime!,INTERVAL_DIL;save_positions=(false,false)) # hack to mark time points for dilution
    cb = CallbackSet(periodcb,contcb,savecb)

    # Solve
    sol = solve(prob,CVODE_BDF(),abstol=1e-10,reltol=1e-10, callback=cb,saveat=tsave);
    return(sol,saved_values)
end

###############################################
####### 5. Coupled model TXTL saturation ######
###############################################

function solvemodel_1RIf(grads,u0,params,TMAX,INTERVAL_DIL,tsave,DIL_FRAC,NSPECIES,SWITCHTIMES,FLAG)
    tspan = (0.0,TMAX);
    prob = ODEProblem(grads,u0,tspan,params);

    # Callbacks
    saved_values = SavedValues(Float64, Array{Float64});
    savecb = SavingCallback((u,t,integrator)->integrator(t,Val{1}),  # Save derivatives as well
        saved_values,saveat=tsave);
    if FLAG=="SR"
       contcb = ContinuousCallback(condition,diluteSR_1RIf!;save_positions=(false,false))
    elseif FLAG=="PC"
        contcb = ContinuousCallback(condition,dilutePC_1RIf!;save_positions=(false,false))
    elseif FLAG=="NC"
        contcb = ContinuousCallback(condition,diluteNC_1RIf!;save_positions=(false,false))
    end
    periodcb = PeriodicCallback(marktime!,INTERVAL_DIL;save_positions=(false,false)) # hack to mark time points for dilution
    cb = CallbackSet(periodcb,contcb,savecb)

    # Solve
    sol = solve(prob,CVODE_BDF(),abstol=1e-10,reltol=1e-10, callback=cb,saveat=tsave);
    return(sol,saved_values)
end


text:
using Test, YaoBlocks, YaoArrayRegister
using YaoBlocks: BlockMap

function heisenberg(n::Int; periodic::Bool = true)
    Sx(i) = put(n, i => X)
    Sy(i) = put(n, i => Y)
    Sz(i) = put(n, i => Z)

    res = map(1:(periodic ? n : n - 1)) do i
        j = mod1(i, n)
        Sx(i) * Sx(j) + Sy(i) * Sy(j) + Sz(i) * Sz(j)
    end
    Add(res)
end

const hm = heisenberg(4)

@testset "constructor:time evolution" begin
    te = TimeEvolution(hm, 0.2)
    # copy
    cte = copy(te)
    @test cte == te
    @test cte !== te

    # dispatch
    dispatch!(cte, [2.0])
    @test cte != te
    @test cte.dt == 2.0
    @test setiparams!(cte, 0.5).dt == 0.5
    @test setiparams!(cte, :random).dt != 0.5
    @test setiparams!(cte, :zero).dt == 0.0
    @test setiparams(cte, 0.5).dt == 0.5
    @test setiparams(cte, :random).dt != 0.5
    @test setiparams(cte, :zero).dt == 0.0
end

@testset "test imaginary time evolution" begin
    tei = TimeEvolution(hm, 0.2im)
    r = rand_state(4)
    r1 = copy(r) |> tei
    @test exp(Matrix(mat(hm)) * 0.2) * r.state ≈ r1.state

    @test applymatrix(tei) ≈ mat(tei)
    @test applymatrix(adjoint(tei)) ≈ applymatrix(tei)'
    @test !isunitary(tei)
    @test !isunitary(tei |> mat)
end

@testset "test time evolution" begin
    te = TimeEvolution(hm, 0.2)

    r = rand_state(4)
    r1 = copy(r) |> te
    @test exp(Matrix(mat(hm)) * -0.2im) * r.state ≈ r1.state
    @test applymatrix(adjoint(te)) ≈ applymatrix(te)'

    @test applymatrix(te) ≈ mat(te)
    @test isunitary(te)
    @test isunitary(te |> mat)
    cte = copy(te)
    @test cte == te
end

@testset "block map" begin
    @test BlockMap(ComplexF64, X) isa BlockMap{ComplexF64,typeof(X)}

    @test BlockMap(ComplexF64, X).block === X

    st = rand(ComplexF64, 2)
    @test BlockMap(ComplexF64, X) * st ≈ mat(X) * st
    @test ishermitian(BlockMap(ComplexF64, X))
    @test size(BlockMap(ComplexF64, X)) == (2, 2)
end

@testset "qudits" begin
    reg = zero_state(1; nlevel=3)
    H = randn(ComplexF64, 3,3)
    reg2 = apply(reg, time_evolve(matblock(H + H'; nlevel=3), 0.5))
    @test nlevel(reg2) == 3
    @test isnormalized(reg2)
end


text:
using PyCall
using Match

pushfirst!(PyVector(pyimport("sys")["path"]), "", "./grammar/", "./src/grammar/")
@pyimport eqn

function parseEqn(input)
    obj = eqn.parse(input)
    eqnToJulia(obj)
end

eqnToJulia(x) = x
eqnToJulia(::Nothing) = nothing
eqnToJulia(str::AbstractString) = Meta.parse(str)

function eqnToJulia(n::Real)
    integer = trunc(Int, n)
    decimal = n - integer
    decimal == 0 ? integer : n
end

function eqnToJulia(obj::PyObject)
    @match obj[:__class__][:__name__] begin
        "EqnAssignment" => eqnAssignment(obj)
        "EqnExpression" => eqnExpression(obj)
        "EqnFunction"   => eqnFunction(obj)
        "EqnIntegral"   => eqnIntegral(obj)
        _ => Meta.ParseError("Unknown PyObject type.")
    end
end

function eqnExpression(obj)
    addBrackets(x) = x
    addBrackets(ex::Expr) = "($(string(ex)))"
    args = map(addBrackets ∘ eqnToJulia, obj[:args])
    exprString = join(args, " ")
    Meta.parse(exprString)
end

function eqnAssignment(obj)
    lvalue = eqnToJulia(obj[:lvalue])
    rvalue = eqnToJulia(obj[:rvalue])
    Assignment(lvalue, rvalue)
end

function eqnFunction(obj)
    name = Symbol(obj[:name])
    args = [eqnToJulia(x) for x in obj[:args]]
    Expr(:call, name, args...)
end

function eqnIntegral(obj)
    ex = eqnToJulia(obj[:intexpr])
    diff = eqnToJulia(obj[:diff])

    a, b = eqnToJulia(obj[:a]), eqnToJulia(obj[:b])
    bounds = isa(a, Nothing) ? nothing : Bounds(a, b)
    createIntegral(ex, diff, bounds)
end


text:
function process!(case, strs)
    fronts = findall(x -> x[1] != '*', strs)

    ends = findall(x -> x[end] != '*', strs)

    if length(fronts) > 2 || length(ends) > 2
        println("Case $(case): *")
        return
    end

    f = String[]

end

function main()
    T = parse(Int, readline(stdin))
    for case in 1:T
        N = parse(Int, readline(stdin))
        strs = String[readline(stdin) for _ in 1:N]

        process!(case, strs)
    end
end

isinteractive() || main()

max_stars_repo_name: AnAverageHuman/competitive

max_issues_repo_name: AnAverageHuman/competitive

max_forks_repo_name: AnAverageHuman/competitive


text:
"""
    ArgumentWrapper

`ArgumentWrapper`s standardize the interface between `AbstractBasis` and the
various arguments with which a basis may be evaluated.
"""
abstract type ArgumentWrapper end

"""
    wrap(x)

Wrap `x` in an `ArgumentWrapper`, if needed.

# Examples
```
wrap(x::Number) = x
wrap(x::AbstractMatrix) = MatrixWrapper(x)
```
"""
wrap(x) = x

"""
    val(x)

Obtain the value of the argument `x`.

# Examples
```
val(x::MatrixVectorWrapper) = x.matrix * x.vector
```
"""
val(x) = x

Base.:+(s::UniformScaling, a::ArgumentWrapper) = a + s
Base.:-(a::ArgumentWrapper, s::UniformScaling) = a + (-s)

Base.Broadcast.broadcastable(a::ArgumentWrapper) = Ref(a)

"""
    MatrixWrapper(M) <: ArgumentWrapper

Wrap the matrix-like object `M` as an argument to `AbstractBasis`.
"""
struct MatrixWrapper{M} <: ArgumentWrapper
    matrix::M
end

wrap(x::AbstractMatrix) = MatrixWrapper(x)

Base.one(a::MatrixWrapper{<:Matrix}) = Matrix{eltype(a.matrix)}(I,size(a.matrix))
Base.one(a::MatrixWrapper{<:Diagonal}) = Diagonal(one.(diag(a.matrix)))
Base.one(a::MatrixWrapper{SparseMatrixCSC{Tv,Ti}}) where {Tv,Ti} = SparseMatrixCSC{Tv,Ti}(I,size(a.matrix))
val(a::MatrixWrapper) = a.matrix

Base.:+(a::MatrixWrapper, s::UniformScaling) = MatrixWrapper(a.matrix + s)
Base.:*(a::MatrixWrapper, b) = a.matrix*b
Base.:\(a::MatrixWrapper, b) = a.matrix\b
Base.inv(a::MatrixWrapper) = inv(a.matrix)
Base.:^(a::MatrixWrapper, p) = a.matrix^p

"""
    MatrixVectorWrapper(M,v, mv_inv = \\) <: ArgumentWrapper

Wrap the matrix-vector pair `M,v` as an argument to `AbstractBasis`.
"""
struct MatrixVectorWrapper{M,V,I} <: ArgumentWrapper
    matrix::M
    vector::V
    mv_inv::I
end

MatrixVectorWrapper(M,v) = MatrixVectorWrapper(M,v,\)

wrap((M,v)::NTuple{2,Any}) = MatrixVectorWrapper(M,v)
Base.one(a::MatrixVectorWrapper) = a.vector
val(a::MatrixVectorWrapper) = a.matrix*a.vector

Base.:+(a::MatrixVectorWrapper, s::UniformScaling) = MatrixVectorWrapper(a.matrix + s, a.vector)
Base.:*(a::MatrixVectorWrapper, b) = a.matrix*b
Base.:\(a::MatrixVectorWrapper, b) = a.matrix\b
Base.inv(a::MatrixVectorWrapper) = a.mv_inv(a.matrix,a.vector)
function Base.:^(a::MatrixVectorWrapper, p::Integer)
    m,v = a.matrix, a.vector
    for i = 1:p
        v = m*v
    end
    for i = 1:-p
        v = m\v
    end
    return v
end


text:
using Dates
using Setfield

using TeaFiles
using TeaFiles.Header: AbstractSection, ContentDescriptionSection, Field, ItemSection,
    NameValue, NameValueSection, TeaFileMetadata, TimeSection, field_type, field_type_id,
    get_section, minimum_item_start, tea_epoch

"""
Quick-and-dirty field construction. Note that this doesn't properly deal with field
offsets, which -- if reading a struct -- are likely to have e.g. 8 byte alignment.
For these tests that doesn't matter though.
"""
function _make_fields(name_types::Vector{Pair{String,DataType}})::Vector{Field}
    offset = Int32(0)
    result = Field[]
    for (name, type) in name_types
        type_id = field_type_id(type)
        push!(result, Field(type_id, offset, name))
        # BEWARE: this is only suitable for these tests. If mapping a structure then one
        # should use fieldoffset(T, i) to compute the offset of field i.
        offset += sizeof(type)
    end
    return result
end
_make_fields(types::Vector{DataType}) = _make_fields(["" => type for type in types])

"""
Quick-and-dirty ItemSection construction. If we are mapping a structure of type T, then we
should instead be making use of sizeof(T).
"""
function TeaFiles.Header.ItemSection(
    item_name::AbstractString,
    fields::AbstractVector{Field}
)
    # BEWARE: this is only suitable for these tests. If we are mapping a structure then the
    #   size on disk could be larger than this due to padding.
    item_size = if isempty(fields) 0 else
        sum(fields) do field
            sizeof(field_type(field))
        end
    end
    return ItemSection(item_size, item_name, fields)
end

"""
Get a metadata section for the example used in the specification.
"""
function _get_example_metadata()::TeaFileMetadata
    return TeaFileMetadata(200, 0, [
        ItemSection(
            "Tick",
            _make_fields(["time" => Int64, "price" => Float64, "volume" => Int64])
        ),
        ContentDescriptionSection("ACME prices"),
        NameValueSection([NameValue("decimals", Int32(2))]),
        TimeSection(719162, 86400000, [0])
    ])
end

"""This is the structure the metadata expects."""
struct Tick
    time::Int64
    price::Float64
    volume::Int64
end

struct DateTimeTick
    time::DateTime
    price::Float64
    volume::Int64
end

function _replace_section!(
    metadata::TeaFileMetadata,
    section::AbstractSection
)::TeaFileMetadata
    # Filter out the old section, and replace it with the new one.
    replace!(x -> isa(x, typeof(section)) ? section : x, metadata.sections)
    return metadata
end

"""
This is the same as the example metadata, but we replace the epoch in the time section
with one that is Julia-compatible.
"""
function _get_example_datetime_metadata()
    julia_epoch_utc = DateTime(0, 1, 1) - Millisecond(Dates.DATETIMEEPOCH)
    julia_epoch = tea_epoch(julia_epoch_utc)

    metadata = _get_example_metadata()

    # Create a new time section instance with the epoch changed.
    time_section = @set get_section(metadata, TimeSection).epoch = julia_epoch

    # Change the name of the item in item section.
    item_section = @set get_section(metadata, ItemSection).item_name = "DateTimeTick"

    _replace_section!(metadata, time_section)
    _replace_section!(metadata, item_section)

    # We've changed the sections, so update the location of the item section.
    return @set metadata.item_start = minimum_item_start(metadata.sections)
end

"""
Convert a vector of some struct T to a vector of equivalent namedtuples.
"""
function _structs_to_namedtuples(x::AbstractVector{T}) where {T}
    out_type = NamedTuple{tuple(fieldnames(T)...), Tuple{fieldtypes(T)...}}
    # Not the fastest, but oh well.
    return reinterpret(out_type, x)
end


text:
using ReachabilityAnalysis
Dict([
    "opts" => Dict(:T=>15.0, :alg=>GLGM06(δ=0.01)),
    "X0"   => Hyperrectangle([2.5, 0.0, 0.0, 0.0, 0.0, 0.0],
                            [2.5, 0.0, 0.2, 0.1, 0.0, 0.0]),
    "info" => Dict(
        "name" => "crane",
        "dim" => 6,
        "linear" => true,
        "hybrid" => false
    )
])

max_stars_repo_name: mforets/ReachabilityLibrary

max_issues_repo_name: mforets/ReachabilityLibrary

max_forks_repo_name: mforets/ReachabilityLibrary


text:
#-------------------------------------------------------------------
#* EMSO Model Library (EML) Copyright (C) 2004 - 2007 ALSOC.
#*
#* This LIBRARY is free software; you can distribute it and/or modify
#* it under the therms of the ALSOC FREE LICENSE as available at
#* http://www.enq.ufrgs.br/alsoc.
#*
#* EMSO Copyright (C) 2004 - 2007 ALSOC, original code
#* from http://www.rps.eng.br Copyright (C) 2002-2004.
#* All rights reserved.
#*
#* EMSO is distributed under the therms of the ALSOC LICENSE as
#* available at http://www.enq.ufrgs.br/alsoc.
#*--------------------------------------------------------------------
#* Author: Gerson Balbueno Bicca
#* $Id$
#*------------------------------------------------------------------
type DoublePipe_NTU
	DoublePipe_NTU()=begin
		new(
			DoublePipe_Basic(),
			DanaSwitcher(Dict{Symbol,Any}(
				:Brief=>"Flow Direction",
				:Valid=>["counter","cocurrent"],
				:Default=>"cocurrent"
			)),
			positive(Dict{Symbol,Any}(
				:Brief=>"Effectiveness estimate",
				:Default=>0.5
			)),
			NTU_Basic(Dict{Symbol,Any}(
				:Brief=>"NTU Method of Calculation",
				:Symbol=>" "
			)),
			[
				:(Method.Eft1 = 1),
				:(Method.NTU*Method.Cmin = _base_1.Details.Ud*_base_1.Geometry.Pi*_base_1.Geometry.DoInner*_base_1.Geometry.Lpipe),
				:(Method.Cmin = min([Method.Ch,Method.Cc])),
				:(Method.Cmax = max([Method.Ch,Method.Cc])),
				:(Method.Cr = Method.Cmin/Method.Cmax),
				:(Method.Eft = 1-exp(-Method.NTU)),
				:(Method.Eft = (1-exp(-Method.NTU*(1+Method.Cr)))/(1+Method.Cr)),
				:(Method.Eft = Method.NTU/(1+Method.NTU)),
				:(Method.Eft = (1-exp(-Method.NTU*(1-Method.Cr)))/(1-Method.Cr*exp(-Method.NTU*(1-Method.Cr)))),
				:(Method.NTU = 1),
				:(Method.Cmin = min([Method.Ch,Method.Cc])),
				:(Method.Cmax = max([Method.Ch,Method.Cc])),
				:(Method.Cr = 1),
				:(Method.Eft = Eftestimated),
				:(_base_1.Details.Q = Method.Eft*Method.Cmin*(_base_1.InletOuter.T-_base_1.InletInner.T)),
				:(_base_1.Details.Q = _base_1.Qestimated),
				:(Method.Ch = _base_1.InletOuter.F*_base_1.Outer.Properties.Average.Cp),
				:(Method.Cc = _base_1.InletInner.F*_base_1.Inner.Properties.Average.Cp),
				:(_base_1.Details.Q = Method.Eft*Method.Cmin*(_base_1.InletInner.T-_base_1.InletOuter.T)),
				:(_base_1.Details.Q = _base_1.Qestimated),
				:(Method.Cc = _base_1.InletOuter.F*_base_1.Outer.Properties.Average.Cp),
				:(Method.Ch = _base_1.InletInner.F*_base_1.Inner.Properties.Average.Cp),
			],
			[
				"Effectiveness Correction","Number of Units Transference","Minimum Heat Capacity","Maximum Heat Capacity","Thermal Capacity Ratio","Effectiveness","Effectiveness in Cocurrent Flow","Effectiveness in Counter Flow","Effectiveness in Counter Flow","Number of Units Transference","Minimum Heat Capacity","Maximum Heat Capacity","Thermal Capacity Ratio","Effectiveness","Duty","Duty","Hot Stream Heat Capacity","Cold Stream Heat Capacity","Duty","Duty","Cold Stream Heat Capacity","Hot Stream Heat Capacity",
			],
			[:FlowDirection,:Eftestimated,],
			[:Method,]
		)
	end
	_base_1::DoublePipe_Basic
	FlowDirection::DanaSwitcher
	Eftestimated::positive
	Method::NTU_Basic
	equations::Array{Expr,1}
	equationNames::Array{String,1}
	parameters::Array{Symbol,1}
	variables::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export DoublePipe_NTU
function setEquationFlow(in::DoublePipe_NTU)
	addEquation(1)
	let switch=CalculationApproach
		if switch=="Full"
			addEquation(2)
			addEquation(3)
			addEquation(4)
			addEquation(5)
			if Method.Cr == 0
				addEquation(6)
			else
				let switch=FlowDirection
					if switch=="cocurrent"
						addEquation(7)
					elseif switch=="counter"
						if Method.Cr == 1
							addEquation(8)
						else
							addEquation(9)
						end
					end
				end
			end
		elseif switch=="Simplified"
			addEquation(10)
			addEquation(11)
			addEquation(12)
			addEquation(13)
			addEquation(14)
		end
	end
	let switch=HotSide
		if HotSide==InletInner.T > InletOuter.T
			set(switch,"inner")
		end
		if HotSide==InletInner.T < InletOuter.T
			set(switch,"outer")
		end
		if switch=="outer"
			let switch=CalculationApproach
				if switch=="Full"
					addEquation(15)
				elseif switch=="Simplified"
					addEquation(16)
				end
			end
			addEquation(17)
			addEquation(18)
		elseif switch=="inner"
			let switch=CalculationApproach
				if switch=="Full"
					addEquation(19)
				elseif switch=="Simplified"
					addEquation(20)
				end
			end
			addEquation(21)
			addEquation(22)
		end
	end
end
function atributes(in::DoublePipe_NTU,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Icon]="icon/DoublePipe"
	fields[:Pallete]=true
	fields[:Brief]="Double Pipe Heat Exchanger - NTU Method"
	fields[:Info]="Thermal analysis of double pipe heat exchanger using the NTU Method.

== Specify ==
* The Inlet Inner stream
* The Inlet Outer stream
== Setting Parameters ==
* Flow Direction:
** counter flow
** cocurrent flow (Default)
* Heat Transfer Correlations:
** Laminar flow
*** Hausen (Default)
*** Schlunder
** Transition flow
*** Gnielinski (Default)
*** Hausen
** Turbulent flow
*** Petukhov (Default)
*** Sieder Tate
* Geometrical Parameters:
** DoInner 	- Outside Diameter of Inner Pipe
** DiInner	- Inside Diameter of Inner Pipe
** DiOuter	- Inside Diameter of Outer pipe
** Lpipe		- Effective Tube Length of one segment of Pipe
** Kwall		- Tube Wall Material Thermal Conductivity
* Fouling:
**Rfi	-	Inside Fouling Resistance
**Rfo	- 	Outside Fouling Resistance
"
	drive!(fields,_)
	return fields
end
DoublePipe_NTU(_::Dict{Symbol,Any})=begin
	newModel=DoublePipe_NTU()
	newModel.attributes=atributes(newModel,_)
	newModel
end


text:
# This is a test for cmdstan issue #547 (thanks to chris fisher & bob)

using Mamba, Stan

ProjDir = dirname(@__FILE__)
cd(ProjDir) do

  simplecode = "
  parameters {
    real y;
  }
  model {
    y ~ normal(0, 1);
  }  "

  global stanmodel, rc, sim
  stanmodel = Stanmodel(Sample(save_warmup=true, thin=1), name="simple", model=simplecode);
  rc, sim, cnames = stan(stanmodel, [Dict("y" => 0.)], CmdStanDir=CMDSTAN_HOME);
  rc == 0 && describe(sim)

end

# The first sample for y (in e.g. tmp/simple_samples_1.csv, ~line 39) should be 0.


text:
using DifferentialEquations, StaticArrays
using HTTP, CSV
using Interpolations
using Unitful, UnitfulAstro

# Downlaod the  data for SLy equation of state
df = CSV.File(HTTP.get("http://www.ioffe.ru/astro/NSG/NSEOS/sly4.dat").body;
    header=false,skipto=7,delim=" ",ignorerepeated=true)

ρ = df.Column3
p = df.Column4

# Interpolate and then extrapolate the data
p_interp = interpolate((ρ, ), p, Gridded(Linear()))
ρ_interp = interpolate((p,), ρ, Gridded(Linear()))
ρ_exterp = extrapolate(ρ_interp, Line())

# Define the tov odes
function tov(u, params, r)
    m, p = u
    G, C = params
    ρ = ρ_exterp(p)
    dmdr = 4 * pi * r^2 * ρ
    dpdr = - (G * m * ρ / r^2) * (1 + (p/(C^2 * ρ))) * (1 + (4* pi * r^3 * p/(C^2 * m))) * (1 - (2 * G * m/(C^2 * r)))^(-1)
    @SVector [dmdr, dpdr]
end

condition(u, r, integrator) = u[2]
affect!(integrator) = terminate!(integrator)
cb = ContinuousCallback(condition, affect!)

ρc = 1e15
u0 = @SVector [1e-10, p_interp(ρc)]
rspan = [1e-10, 1e10]
G = 6.6743 * 1e-8
c = 3 * 1e10

problem = ODEProblem(tov, u0, rspan, [G, c])
sol = solve(problem, callback=cb)

Msun = uconvert(NoUnits, 1 * u"Msun"/u"g")
sol.t[end]/1e5
sol.u[end][2]
sol.u[end][1]/Msun

# We solve for different values of the central density and plot the mass and radius of the star
ρcs = 10 .^ LinRange(14.5, log10(ρ[end]), 50)

Masses = zeros(length(ρcs))
Radii = zeros(length(ρcs))

for (idx, ρc) in enumerate(ρcs)
    u0 = @SVector [1e-10, p_interp(ρc)]
    problem = ODEProblem(tov, u0, rspan, [G, c])
    sol = solve(problem, callback=cb)
    Radii[idx] = sol.t[end]/1e5
    Masses[idx] = sol.u[end][1]/Msun
end

using Plots, LaTeXStrings
scatter(Radii, Masses, xlabel=L"$r \textrm{ [km]}$", ylabel=L"$\textrm{Mass }[M_\odot]$", label=L"\textrm{SLy}")
savefig("mass_vs_radius_SLy.png")


text:
using DataDeps

register(DataDep(
  "acute-inflammations-1",
  "http://archive.ics.uci.edu/ml/datasets/Acute+Inflammations",
  "http://archive.ics.uci.edu/ml/machine-learning-databases/acute/diagnosis.data",
  "5d0b61049992b8327aa847f916f3389ec08f74b43dbb3fc5222ac07185d7c88f",
  post_fetch_method=(path -> begin
    # First 2 bytes are garbage
    data = LegacyStrings.utf16(read(path)[3:end])
    data = replace(data, "," => ".")
    write(path, data)

    UCIData.process_dataset(path,
        target_index=7,
        feature_indices=1:6,
        categoric_indices=2:6,
        delim='\t',
    )
  end),
))

max_stars_repo_name: JackDunnNZ/uci-data

max_issues_repo_name: JackDunnNZ/uci-data

max_forks_repo_name: JackDunnNZ/uci-data


text:
## Batch result figure: MLE accuracy and spectra
using EstimateSDE
using Serialization
using Pipe: @pipe
using Plots
pgfplotsx()
using DataFrames
using Statistics
using StatsPlots
using REPL.TerminalMenus
using Printf
using Random

## Settings
name = "experiment01"
Tkld = 0.2
nspec = 500 # Number of points in a power spectrum
npsd = 5 # Number of power spectra plotted
ssize = (300,300) # figure size

## Load settings and aggregated results
results_dir = settings_resultsdir(name)
S = deserialize(joinpath(results_dir,"settings.jld"))
R = deserialize(joinpath(results_dir,"aggregate.jld")).df
df_eval = deserialize(joinpath(results_dir,"model-evaluation.jld")).df
df_data = deserialize(joinpath(results_dir,"data.jld")).df
df_eval = @pipe(
  innerjoin(df_eval,df_data,on=:index) |>
  rename(_,:label=>:kld_label,:T=>:kld_T))

figure_dir = joinpath(results_dir,"figures")
if !isdir(figure_dir) mkpath(figure_dir) end

## Preprocessing for plotting
em = Dict(SdeRootsEngine=>"root",ElimQ=>"coef")
for dfp in [df_eval,R]
  dfp.procname = map(x->split(String(x),"_")[end],dfp.proclabel)
  dfp.param = [em[x] for x in dfp.engine]
end

## KLD for reference MLE
dfr = @pipe(
  filter(R) do x
    x.method=="MLE" && x.param=="root" &&
    x.procname∈["A","B","C"] &&
    x.kld_label=="regular" && x.kld_T==Tkld
  end |>
  select(_,[:procname,:kld,:kld_std])
)

klds = Dict(zip(dfr.procname,dfr.kld))

## Plot relative KLD vs process
plt_kld = plot(xlabel="process",ylabel="KLD",size=ssize)
plot!(seriestype=:bar,dfr.procname,dfr.kld,linecolor=:black,bar_width=0.5,yerror=dfr.kld_std,
  label=nothing)
plot!(legend=:none, size=(350,300),xlabel="process",ylabel="KLD")

## Power spectra for representative simulation runs
plts_psd = []
for procname in ["A","B","C"]
  # Select case
  dfr = @pipe(
    filter(df_eval) do x
      x.procname==procname &&
      x.method=="MLE" && x.param=="root" &&
      x.kld_label=="regular" && x.kld_T==Tkld
    end |>
    select(_,[:proclabel,:procname,:data,:N,:Tav,:est,:kld_T,:kld])
  )

  # Plot actual psd
  rr = first(dfr)
  proc = testprocess(rr.proclabel)
  duration = rr.N*rr.Tav
  t_kld = 0:rr.kld_T:duration
  f = (0:nspec-1)./nspec/2Tkld
  h = powerspectrum(proc,f)
  plt = plot(f,h,label="process $(rr.procname)",linewidth=2)
  plot!(yscale=:log10,xlabel="f",ylabel="psd",size=ssize,
      fg_legend=nothing,bg_legend=nothing,legend=:topright)
  push!(plts_psd,plt)

  # Plot `npsd` power spectra
  rng = MersenneTwister(234230)
  iir = rand(rng,1:nrow(dfr),npsd)
  for i in 1:npsd
    row_eval = dfr[iir[i],:]
    ĥ = powerspectrum(row_eval.est,f)
    lab = i==1 ? @sprintf("MLE D %0.0f",klds[procname]) : nothing
    plot!(f,ĥ,label=lab,color=:orange,linealpha=0.8)
  end
  fn = 1/2S.data.Tav
  vline!([fn],linecolor=:black,linestyle=:dash,label=nothing)
end

# Compbine into 2x2 figure
plt = plot(plt_kld,plts_psd...;layout=(2,2),size=2 .* ssize)
display(plt)

## Save
for ext in ["tex","pdf"]
  f = "mle-accuracy.$ext"
  savefig(plt,joinpath(figure_dir,f))
  @info "Saved $f"
end


text:
import DSP: filt, filtfilt, resample

export fir, removedc, removedc!, demon
export upconvert, downconvert, rrcosfir, rcosfir
export mseq, gmseq, circconv, goertzel, pll
export sfilt, sfiltfilt, sresample, mfilter
export istft, whiten

"""
$(SIGNATURES)
Designs a `n`-tap FIR filter with a passband from `f1` to `f2` using the
specified `method`. If frame rate `fs` is not specified, `f1` and `f2` are given
in normalized units (1.0 being Nyquist). If `f1` is 0, the designed filter is
a lowpass filter, and if `f2` is `nothing` then it is a highpass filter.

This method is a convenience wrapper around `DSP.digitalfilter`.

# Examples:
```julia-repl
julia> lpf = fir(127, 0, 10kHz; fs=44.1kHz)   # design a lowpass filter
127-element Array{Float64,1}:
  ⋮

julia> hpf = fir(127, 10kHz; fs=44.1kHz)      # design a highpass filter
127-element Array{Float64,1}:
  ⋮

julia> bpf = fir(127, 1kHz, 5kHz; fs=44.1kHz) # design a bandpass filter
127-element Array{Float64,1}:
  ⋮
```
"""
function fir(n, f1, f2=nothing; fs=2.0, method=FIRWindow(hanning(n)))
  fs = inHz(fs)
  if f1 == 0
    f = Lowpass(inHz(f2); fs=fs)
  elseif f2 == nothing || inHz(f2) == fs/2
    f = Highpass(inHz(f1); fs=fs)
  else
    f = Bandpass(inHz(f1), inHz(f2); fs=fs)
  end
  return digitalfilter(f, method)
end

"""
$(SIGNATURES)
DC removal filter. Parameter `α` controls the cutoff frequency. Implementation
based on Lyons 2011 (3rd ed) real-time DC removal filter in Fig. 13-62(d).

See also: [`removedc`](@ref)
"""
function removedc!(s; α=0.95)
  for k = 1:size(s,2)
    for j = 2:size(s,1)
      s[j,k] += α*s[j-1,k]
    end
    s[2:end,k] .+= -s[1:end-1,k]
  end
  s *= sqrt(α)
  return s
end

"""
$(SIGNATURES)
DC removal filter. Parameter `α` controls the cutoff frequency. Implementation
based on Lyons 2011 (3rd ed) real-time DC removal filter in Fig. 13-62(d).

See also: [`removedc!`](@ref)
"""
removedc(s; α=0.95) = removedc!(copy(s); α=α)

"""
$(SIGNATURES)
Estimates DEMON spectrum. The output is highpass filtered with a `cutoff`
frequency and downsampled. Supported downsampling methods are `:rms` (default),
`:mean` and `:fir`.
"""
function demon(x; fs=framerate(x), downsample=250, method=:rms, cutoff=1.0)
  local y
  fs /= downsample
  for k = 1:size(x,2)
    if downsample == 1
      y1 = abs.(hilbert(x[:,k]))
    elseif method == :rms
      y1 = sqrt.(mean.(Periodograms.arraysplit(abs2.(hilbert(x[:,k])), downsample, 0)))
    elseif method == :mean
      y1 = mean.(Periodograms.arraysplit(abs.(hilbert(x[:,k])), downsample, 0))
    elseif method == :fir
      aaf = fir(127, 0, 0.48fs; fs=fs)
      y1 = filtfilt(aaf, abs.(hilbert(x[:,k])))[1:downsample:end]
    else
      throw(ArgumentError("Unknown method"))
    end
    if k == 1
      y = zeros(length(y1), size(x,2))
    end
    y[:,k] .= y1
  end
  maxflen = length(y)÷4
  mod(maxflen, 2) == 0 && (maxflen += 1)
  hpf = fir(min(127, maxflen), cutoff; fs=fs)
  signal(filtfilt(hpf, y), fs)
end

function demon(x::AbstractVector{T}; fs=framerate(x), downsample=250, method=:rms, cutoff=1.0) where T
  y = @view samples(x)[:,1:1]
  z = demon(y, fs=fs, downsample=downsample, method=method, cutoff=cutoff)
  signal(dropdims(samples(z), dims=2), framerate(z))
end

"""
$(SIGNATURES)
Converts baseband signal with `sps` symbols per passband sample to a real
passband signal centered around carrier frequency `fc`.
"""
function upconvert(s::AbstractVector, sps, fc, pulseshape=rrcosfir(0.25, sps); fs=framerate(s))
  pad = cld(length(pulseshape), 2*sps) - 1
  s = vcat(zeros(pad), analytic(s), zeros(pad))
  s = signal(resample(s, sps, pulseshape), sps*fs)
  √2 * real.(s .* cis.(2π * inHz(fc) * domain(s)))
end

function upconvert(s::AbstractMatrix, sps, fc, pulseshape=rrcosfir(0.25, sps); fs=framerate(s))
  out = Any[]
  for j in 1:nchannels(s)
    push!(out, upconvert(s[:,j], sps, fc, pulseshape; fs=fs))
  end
  hcat(out...)
end

"""
$(SIGNATURES)
Converts passband signal centered around carrier frequency `fc` to baseband,
and downsamples it by a factor of `sps`. If the `pulseshape` is specified to
be `nothing`, downsampling is performed without filtering.
"""
function downconvert(s::AbstractVector, sps, fc, pulseshape=rrcosfir(0.25, sps); fs=framerate(s))
  s = signal(analytic(s), fs)
  s .*= cis.(-2π * inHz(fc) * domain(s))
  sps == 1 && return signal(s, fs)
  pulseshape == nothing && return signal(s[1:sps:end,:], fs/sps)
  signal(resample(s, 1//sps, pulseshape), fs/sps)
end

function downconvert(s::AbstractMatrix, sps, fc, pulseshape=rrcosfir(0.25, sps); fs=framerate(s))
  out = Any[]
  for j in 1:nchannels(s)
    push!(out, downconvert(s[:,j], sps, fc, pulseshape; fs=fs))
  end
  hcat(out...)
end

"""
$(SIGNATURES)
Root-raised cosine filter.
"""
function rrcosfir(β, sps, span = β < 0.68 ? 33-floor(Int, 44β) : 4)
  # default span based on http://www.commsys.isy.liu.se/TSKS04/lectures/3/MichaelZoltowski_SquareRootRaisedCosine.pdf
  delay = fld(span*sps, 2)
  t = collect(-delay:delay)/sps
  h = Array{Float64}(undef, size(t))
  for i ∈ 1:length(t)
    if t[i] == 0
      h[i] = (1 + β*(4/π - 1))/sps
    elseif abs(t[i]) == 1/(4β)
      h[i] = β/(√2*sps) * ((1+2/pi)*sin(π/(4β)) + (1-2/pi)*cos(π/(4β)))
    else
      h[i] = (sin(π*t[i]*(1-β)) + 4β*t[i]*cos(π*t[i]*(1+β))) / (π*t[i]*(1 - (4β*t[i])^2)) / sps
    end
  end
  h / √sum(h.^2)
end

"""
$(SIGNATURES)
Raised cosine filter.
"""
function rcosfir(β, sps, span = β < 0.68 ? 33-floor(Int, 44β) : 4)
  # default span based on http://www.commsys.isy.liu.se/TSKS04/lectures/3/MichaelZoltowski_SquareRootRaisedCosine.pdf
  # since the span is for rrcosfir, for rcosfir, it is very conservative
  delay = fld(span*sps, 2)
  t = collect(-delay:delay)/sps
  h = Array{Float64}(undef, size(t))
  for i ∈ 1:length(t)
    if abs(t[i]) == 1/(2β)
      h[i] = π/(4sps) * sinc(1/(2β))
    else
      h[i] = sinc(t[i]) * cos(π * β * t[i]) / (1-(2β * t[i])^2) / sps
    end
  end
  h / √sum(h.^2)
end

"""
$(SIGNATURES)
Generates an m-sequence of length `2^m-1` or tap specification `m`.

m-sequences are sequences of `+1/-1` values with near-perfect discrete periodic
auto-correlation properties. All non-zero lag periodic auto-correlations
are -1. The zero-lag autocorrelation is `2^m-1`, where `m` is the shift register
length.

This function currently supports shift register lengths between 2 and 30.

# Examples:
```julia-repl
julia> x = mseq(3)                  # generate regular m-sequence
7-element Array{Float64,1}:
  1.0
  1.0
  1.0
 -1.0
  1.0
 -1.0
 -1.0

julia> x = mseq((1,3))              # generate m-sequence with specification (1,3)
7-element Array{Float64,1}:
  1.0
  1.0
  1.0
 -1.0
  1.0
 -1.0
 -1.0
```
"""
function mseq(m, θ=π/2)
  knownspecs = Dict(  # known m-sequences are specified as base 1 taps
       2 => (1,2),          3 => (1,3),          4 => (1,4),          5 => (2,5),
       6 => (1,6),          7 => (1,7),          8 => (1,2,7,8),      9 => (4,9),
      10 => (3,10),        11 => (9,11),        12 => (6,8,11,12),   13 => (9,10,12,13),
      14 => (4,8,13,14),   15 => (14,15),       16 => (4,13,15,16),  17 => (14,17),
      18 => (11,18),       19 => (14,17,18,19), 20 => (17,20),       21 => (19,21),
      22 => (21,22),       23 => (18,23),       24 => (17,22,23,24), 25 => (22,25),
      26 => (20,24,25,26), 27 => (22,25,26,27), 28 => (25,28),       29 => (27,29),
      30 => (7,28,29,30)
  )
  if m ∈ keys(knownspecs)
    spec = collect(knownspecs[m])
  else
    spec = collect(m)
    m = maximum(spec)
  end
  n = 2^m - 1
  reg = ones(UInt8, m)
  x = zeros(Float64, n)
  for j ∈ 1:n
    b = ⊻(reg[spec]...)
    reg = circshift(reg, 1)
    x[j] = 2.0*reg[1] - 1.0
    reg[1] = b
  end
  return x
end

"""
$(SIGNATURES)
Generates an generalized m-sequence of length `2^m-1` or tap specification `m`.

Generalized m-sequences are related to m-sequences but have an additional parameter
`θ`. When `θ = π/2`, generalized m-sequences become normal m-sequences. When
`θ < π/2`, generalized m-sequences contain a DC-component that leads to an exalted
carrier after modulation. When `θ` is `atan(√(2^m-1))`, the m-sequence
is considered to be _period matched_. Period matched m-sequences are complex sequences
with perfect discrete periodic auto-correlation properties, i.e., all non-zero lag
periodic auto-correlations are zero. The zero-lag autocorrelation is `2^m-1`,
where `m` is the shift register length.

This function currently supports shift register lengths between 2 and 30.

# Examples:
```julia-repl
julia> x = gmseq(3)         # generate period matched m-sequence
7-element Array{Complex{Float64},1}:
 0.3535533905932738 + 0.9354143466934853im
 0.3535533905932738 + 0.9354143466934853im
 0.3535533905932738 + 0.9354143466934853im
 0.3535533905932738 - 0.9354143466934853im
 0.3535533905932738 + 0.9354143466934853im
 0.3535533905932738 - 0.9354143466934853im
 0.3535533905932738 - 0.9354143466934853im

julia> x = gmseq(3, π/4)    # generate m-sequence with exalted carrier
7-element Array{Complex{Float64},1}:
 0.7071067811865476 + 0.7071067811865475im
 0.7071067811865476 + 0.7071067811865475im
 0.7071067811865476 + 0.7071067811865475im
 0.7071067811865476 - 0.7071067811865475im
 0.7071067811865476 + 0.7071067811865475im
 0.7071067811865476 - 0.7071067811865475im
 0.7071067811865476 - 0.7071067811865475im
```
"""
function gmseq(m, θ=atan(√(2^maximum(m)-1)))
  x = mseq(m) .+ 0im
  cos(θ) .+ 1im * sin(θ) .* x
end

"""
$(SIGNATURES)
Computes the circular convolution of `x` and `y`. Both vectors must be the same
length.
"""
function circconv(x::AbstractVector, y::AbstractVector=x)
  if length(x) != length(y)
    throw(ArgumentError("x and y must be of equal length"))
  end
  n = length(x)
  z = similar(x)
  for j ∈ 1:n
    z[j] = circshift(x, j-1)' * y
  end
  return z
end

"""
$(SIGNATURES)
Detects frequency `f` in input signal using the Goertzel algorithm.

The detection metric returned by this function is the complex output
of the Goertzel filter at the end of the input block. Typically, you
would want to compare the magnitude of this output with a threshold to
detect a frequency.

When a block size `n` is specified, the Goertzel algorithm in applied to
blocks of data from the original time series.
"""
function goertzel(x::AbstractVector, f, n; fs=framerate(x))
  signal(map(x1 -> goertzel(x1, f; fs=fs), partition(x, n)), fs/n)
end

function goertzel(x::AbstractVector, f; fs=framerate(x))
  n = length(x)
  m = inHz(f)/(inHz(fs)/n)
  w1 = 0
  w2 = 0
  for j ∈ 1:n
    w0 = 2 * cos(2π * m/n) * w1 - w2 + x[j]
    w2 = w1
    w1 = w0
  end
  w0 = 2 * cos(2π * m/n) * w1 - w2
  w0 - cis(-2π * m/n) * w1
end

function goertzel(x::AbstractMatrix, f, n; fs=framerate(x))
  count = cld(size(x,1), n)
  out = Array{ComplexF64}(undef, (count, nchannels(x)))
  for j ∈ 1:nchannels(x)
    out[:,j] = goertzel(x[:,j], f, n; fs=fs)
  end
  signal(out, fs/n)
end

function goertzel(x::AbstractMatrix, f; fs=framerate(x))
  out = Array{ComplexF64}(undef, nchannels(x))
  for j ∈ 1:nchannels(x)
    out[j] = goertzel(x[:,j], f; fs=fs)
  end
  out
end

"""
$(SIGNATURES)
Phased-lock loop to track dominant carrier frequency in input signal.
"""
function pll(x::AbstractVecOrMat, bandwidth=1e-3; fs=framerate(x))
  β = √bandwidth
  n = nchannels(x)
  ϕ = zeros(1,n)
  ω = zeros(1,n)
  y = similar(x, ComplexF64)
  for j ∈ 1:nframes(x)
    y[j,:] = cis.(ϕ)
    Δϕ = angle.(x[j,:] .* conj.(y[j,:])) .* abs.(x[j,:])
    ω .+= bandwidth * Δϕ
    ϕ .+= β*Δϕ .+ ω
  end
  signal(y, fs)
end

"""
$(SIGNATURES)
Same as [`filt`](https://docs.juliadsp.org/stable/filters/#DSP.filt),
but retains sampling rate information.
"""
sfilt(f, x, args...) = signal(filt(f, samples(x), args...), framerate(x))
sfilt(b, a, x, args...) = signal(filt(b, a, samples(x), args...), framerate(x))

"""
$(SIGNATURES)
Same as [`filtfilt`](https://docs.juliadsp.org/stable/filters/#DSP.Filters.filtfilt),
but retains sampling rate information.
"""
sfiltfilt(coef, x) = signal(filtfilt(coef, samples(x)), framerate(x))

"""
$(SIGNATURES)
Same as [`resample`](https://docs.juliadsp.org/stable/filters/#DSP.Filters.resample),
but correctly handles sampling rate conversion.
"""
sresample(x, rate, args...) = signal(resample(samples(x), rate, args...), rate * framerate(x))

"""
$(SIGNATURES)
Matched filter looking for reference signal `r` in signal `s`.
"""
function mfilter(r, s)
  f = conj.(reverse(samples(r)))
  n = length(r) - 1
  sfilt(f, padded(s, (0, n)))[n+1:end]
end

"""
$(SIGNATURES)
Compute the inverse short time Fourier transform (ISTFT) of STFT coefficients `X` which is based
on segments with `nfft` samples with overlap of `noverlap` samples. Refer to `DSP.Periodograms.spectrogram`
for description of the parameters.

For perfect reconstruction, the parameters `nfft`, `noverlap` and `window` in `stft` and
`istft` have to be the same, and the windowing must obey the constraint of "nonzero overlap add" (NOLA).
Implementation based on Zhivomirov 2019 and `istft` in `scipy`.

# Examples:
```julia-repl
julia> x = randn(1024)
1024-element Array{Float64,1}:
 -0.7903319156212055
 -0.564789077302601
  0.8621044972211616
  0.9351928359709288
  ⋮
  2.6158861993992533
  1.2980813993011973
 -0.010592954871694647
julia> X = stft(x, 64, 0)
33×31 Array{Complex{Float64},2}:
  ⋮
julia> x̂ = istft(X; nfft=64, noverlap=0)
1024-element Array{Float64,1}:
 -0.7903319156212054
 -0.5647890773026012
  0.8621044972211612
  0.9351928359709288
  ⋮
  2.6158861993992537
  1.2980813993011973
 -0.010592954871694371
```
"""
function istft(X::AbstractMatrix{Complex{T}};
               nfft::Int,
               noverlap::Int,
               onesided::Bool=true,
               window::Union{Function,AbstractVector,Nothing}=nothing) where {T<:AbstractFloat}
  # H. Zhivomirov, TEM Journal, Vol. 8, No. 1, pp. 56-64, 2019.
  (window === nothing) && (window = rect)
  win, norm2 = Periodograms.compute_window(window, nfft)

  nstep = nfft - noverlap
  nseg = size(X, 2)
  outputlength = nfft + (nseg-1) * nstep

  iX = onesided ? irfft(X, nfft, 1) : ifft(X, 1)
  iX .*= win
  x = zeros(eltype(iX), outputlength)
  normw = zeros(T, outputlength)
  for i = 1:nseg
      x[1+(i-1)*nstep:nfft+(i-1)*nstep] += iX[:,i]
      normw[1+(i-1)*nstep:nfft+(i-1)*nstep] += win .^ 2
  end

  (sum(normw[nfft÷2:end-nfft÷2] .> 1e-10) != length(normw[nfft÷2:end-nfft÷2])) && (
      @warn "NOLA condition failed, STFT may not be invertible")
  x .*= nstep/norm2
end

"""
$(SIGNATURES)
Spectral whitening of input signal `x` in the frequency domain. The parameters `nfft`,
`noverlap` and `window` are required for the computation of STFT coefficients of `x`.
Refer to `DSP.Periodograms.spectrogram` for description of the parameters. `γ` is a
scaling or degree-of-flattening factor. The algorithm is based on Lee 1986.
"""
function whiten(x::AbstractVector;
                nfft::Int,
                noverlap::Int,
                window::Union{Function,AbstractVector,Nothing}=nothing,
                γ=1)
  # M. W. Lee, Open-File Report 86-108, 1986.
  xstft = stft(x, nfft, noverlap; window=window)
  mag = abs.(xstft)
  logmag = log.(mag .+ eps(eltype(mag)))
  logmag .-= γ * mean(logmag; dims=2)
  istft(exp.(logmag) .* exp.(im .* angle.(xstft)); nfft=nfft, noverlap=noverlap, window=window)
end


text:
"""
    Environment([config::Config], [name::String], [root::AbstractPath])

An environment stores information about a playground environment and
provides methods fro interacting with them.

NOTE: In the future we might want to support different types of environments
(e.g., DockerEnvironment).
"""
type Environment
    config::Config
    name::AbstractString
    root::AbstractPath
end

Environment() = Environment(Config())

function Environment(config::Config)
    root = envpath(config)
    name = envname(config, root)
    Environment(config, name, root)
end

Environment(name::String) = Environment(Config(), name)
Environment(dir::AbstractPath) = Environment(Config(), dir)
Environment(config::Config, name::String) = Environment(config, name, envpath(config, name))

function Environment(config::Config, dir::AbstractPath)
    root = envpath(config, dir)
    name = envname(config, root)
    Environment(config, name, root)
end

function Environment(config::Config, dir::AbstractPath, name::String)
    p = if isempty(dir) && isempty(name)
        envpath(config)
    elseif isempty(dir) && !isempty(name)
        envpath(config, name)
    else
        envpath(config, dir)
    end

    debug(logger, "Environment: Name=$name, Path=$p")
    return Environment(config, name, p)
end

function init(env::Environment)
    for f in (root, bin, log, pkg)
        p = f(env)

        if !exists(p)
            mkdir(p; recursive=true)
            debug(logger, "$p created.")
        end
    end

    # If the playground name is set and the root path isn't already in the
    # shared folder then create a symlink.
    share_path = abs(join(env.config.share, env.name))
    if !isempty(env.name) && env.root != share_path
        symlink(env.root, share_path, exist_ok=true, overwrite=true)
    end
end

name(env::Environment) = env.name
root(env::Environment) = env.root
log(env::Environment) = join(env.root, p"log")
bin(env::Environment) = join(env.root, p"bin")
pkg(env::Environment) = join(env.root, p"packages")
depot(env::Environment) = join(env.root, p"depot")
julia(env::Environment) = join(env.root, p"bin", p"julia")
shell(env::Environment) = env.config.default_shell
function history(env::Environment, lang::Symbol)
    if lang === :shell
        env.config.isolated_shell_history
    elseif lang === :julia
        env.config.isolated_julia_history
    else
        error(logger, "Unsupported history language $lang")
    end
end

"""
    getenvs(env::Environmnet) -> Vector{Pair}

Generates a set of environment variable changes to make from the
`Environment`.
"""
function getenvs(env::Environment)
    envs = Pair[]
    path = "$(bin(env)):" * ENV["PATH"]
    envname = name(env) == "" ? basename(root(env)) : name(env)
    push!(
        envs,
        "PATH" => path,
        "PLAYGROUND_ENV" => envname,
        "JULIA_PKGDIR" => pkg(env),
        "JULIA_DEPOT_PATH" => depot(env),     # For Pkg3 support
    )

    if !isempty(shell(env))
        push!(envs, "SHELL" => shell(env))
    end

    if history(env, :julia)
        jh = join(root(env), p".julia_history")
        push!(envs, "JULIA_HISTORY" => jh)
    end

    if history(env, :shell)
        sh = join(root(env), p".shell_history")
        push!(envs, "HISTFILE" => sh)
    end

    return envs
end

"""
    set!(env::Environment, keyvals::Pair...) -> Dict

Applies a the `keyvals` to `ENV` and returns the old
`ENV` settings to be used for restore the `ENV` state.
"""
function set!(env::Environment, keyvals::Pair...)
    old = Dict{AbstractString, Any}()

    for (key, val) in keyvals
        logenv(key, val)
        old[key] = get(ENV, key, nothing)
        ENV[key] = val
    end

    return old
end

"""
    restore!(old::Dict{AbstractString, Any})

Applies the `old` dict to `ENV` to restore the environment state.
"""
function restore!(old::Dict{AbstractString, Any})
    for (key, val) in old
        if val !== nothing
            logenv(key, val)
            ENV[key] = val
        else
            debug(logger, "Deleted $key")
            delete!(ENV, key)
        end
    end
end

"""
    logenv(key, val)

Internal utility method used for nicely logging environment variable
changes.
"""
function logenv(key, val)
    if key == "PATH"
        paths = join(unique(split(val, ":")), "\n\t")
        debug(logger, "Set $key = \n\t$paths")
    else
        debug(logger, "Set $key = $val")
    end
end

"""
    withenv(f::Function, env::Environment)

Works the same as `withenv(f, keyvals...)`, but is specific to running `f` within
a playground `Environment`.
"""
function Base.withenv(f::Function, env::Environment)
    old = set!(env, getenvs(env)...)
    try
        f()
    finally
        restore!(old)
    end
end


text:
module BodePlots

using InspectDR
import InspectDR: Plot2D
import Printf: @sprintf
using NumericIO
using Colors

#Colour constants
#-------------------------------------------------------------------------------
const black = RGB24(0, 0, 0); const white = RGB24(1, 1, 1)
const red = RGB24(1, 0, 0); const green = RGB24(0, 1, 0); const blue = RGB24(0, 0, 1)

#Convenience aliases
#-------------------------------------------------------------------------------
mag(x) = abs(x); phase(x) = rad2deg(angle(x))
SI(x) = formatted(x, :SI, ndigits=3)
const j = im

#==Math: 2nd order transfer functions
===============================================================================#

function roots(a,b,c) #Solves: ax² + bx + c = 0
	r = sqrt(b*b-4*a*c)
	result = (-b .+ [r, -r]) / (2a)
	return (result[1], result[2])
end

#2nd order transfer function:
xf_o2(f, G, f_z, f_p1, f_p2) = #(f_z, f_p1, f_p2: in Hz)
	G*(1 .+ j*(f./f_z))./( (1 .+ j*(f./f_p1)).*(1 .+ j*(f./f_p2)) )

#Compute unity gain frequency
function f_0(f, G, f_z, f_p1, f_p2)
	invsq(x) = (result = 1/x; return result*result)
	_z² = invsq(f_z); G² = G*G
	_p1² = invsq(f_p1); _p2² = invsq(f_p2)
	r = roots(_p1²*_p2², _p1²+_p2²-G²*_z², 1-G²)
	posroot = max(r...) #Select positive root
	return posroot < 0 ? NaN : sqrt(posroot)
end

#Compute 3dB frequency
f_3dB(f, G, f_z, f_p1, f_p2) =
	f_0(f, sqrt(2), f_z, f_p1, f_p2)

#==Main functions
===============================================================================#

#Create new Bode plot:
#-------------------------------------------------------------------------------
function new(::Type{InspectDR.Plot})
	plot = InspectDR.bodeplot()
		strip_mag, strip_phase = plot.strips

	#Define extents & scales:
	strip_mag.yext_full = InspectDR.PExtents1D(-10, 25)
	strip_phase.yext_full = InspectDR.PExtents1D(-180, 0)

	return plot
end

function new()
	w = 500; h = w/1.6 #WANTCONST: Select plot width/height

	mplot = InspectDR.Multiplot()
	mplot.layout[:ncolumns] = 1
	#Bode plot looks better with wider aspect ratio:
	mplot.layout[:halloc_plot] = w; mplot.layout[:valloc_plot] = h

	add(mplot, new(InspectDR.Plot))
	return mplot
end

#Modify previously-generated Bode plot (base algorithm):
#-------------------------------------------------------------------------------
function update(plot::Plot2D, f::Vector, y::Vector)
	ldata = line(color=blue, width=3, style=:solid) #WANTCONST
	strip_mag, strip_phase = plot.strips

	plot.data = [] #Clear old data

	#Reset extents:
	plot.xext = InspectDR.PExtents1D()
	strip_mag.yext = InspectDR.PExtents1D()
	strip_phase.yext = InspectDR.PExtents1D()

	wfrm = add(plot, f, mag.(y), strip=1)
		wfrm.line = ldata
	wfrm = add(plot, f, phase.(y), strip=2)
		wfrm.line = ldata

	return
end

#Update annotation on Bode plot (enabled=false to clear annotation):
#-------------------------------------------------------------------------------
function update_annotation(plot::Plot2D, f0, fBW, phase0, enabled=true)
	lmarker = line(style=:dash, width=2.5) #WANTCONST
	lmarker_light = line(style=:dash, width=2.5, color=RGB24(.4,.4,.4)) #WANTCONST

	plot.userannot = [] #Clear old markers/text annotation

	if !enabled; return; end

	#Add specialized annotation
	afont = plot.layout[:font_annotation]

	#Add vertical markers to both plots:
	isfinite(f0) ? add(plot, vmarker(f0, lmarker, strip=0)) : nothing
	isfinite(fBW) ? add(plot, vmarker(fBW, lmarker, strip=0)) : nothing

	#Add annotation to Magnitude plot:
	add(plot, atext("0dB", y=0, xoffset_rel=0.5, yoffset=3, font=afont, align=:bc, strip=1))
	add(plot, hmarker(0, lmarker_light, strip=1))

	if isfinite(fBW)
		fstr = "f3dB=" * SI(fBW) * "Hz"
		add(plot, atext(fstr, x=fBW, xoffset=-3, yoffset_rel=0.5, font=afont, angle=-90, align=:bc, strip=1))
	end

	#Hack if f0>0: TODO: Fix bug where annotating @ 0 on log plot causes glitch.
	if isfinite(f0) && f0 > 0
		fstr = "f0=" * SI(f0) * "Hz"
		add(plot, atext(fstr, x=f0, xoffset=-3, yoffset_rel=0.5, font=afont, angle=-90, align=:bc, strip=1))
	end

	#Add annotation to Phase plot:
	if isfinite(phase0)
		pmargin = 180+phase0
		fstr = "PM=" * @sprintf("%.1f°", pmargin)
		add(plot, atext(fstr, y=(phase0-180)/2, xoffset_rel=0.5, font=afont, align=:cc, strip=2))
		add(plot, hmarker(phase0, lmarker, strip=2))
	end

	return
end

#Modify previously-generated Bode plot:
#-------------------------------------------------------------------------------
#(fmin, fmax, f_p1: in Hz)
function update(plot::Plot2D, fmin, fmax, G, f_z, f_p1, f_p2, annot=true)
	npts = 100 #WANTCONST

	_logspace(start, stop, n) = 10 .^ range(start, stop=stop, length=n)
	f = _logspace(log10(fmin), log10(fmax), npts)
	y = xf_o2(f, G, f_z, f_p1, f_p2)
	update(plot, f, y)

	f0 = f_0(f, G, f_z, f_p1, f_p2) #Unity gain frequency
	fBW = f_3dB(f, G, f_z, f_p1, f_p2)
	phase0 = phase(xf_o2(f0, G, f_z, f_p1, f_p2)) #Unity gain phase
	update_annotation(plot, f0, fBW, phase0, annot)

	return
end

#Modify previously-generated Bode plot:
#-------------------------------------------------------------------------------
#(fmin, fmax, f_p1: in Hz)
function update(mplot::InspectDR.Multiplot, fmin, fmax, GdB, f_z, f_p1, f_p2, annot=true)
	plot = mplot.subplots[1] #WANTCONST
	G = 10.0^(GdB/20) #WANTCONST
	update(plot, fmin, fmax, G, f_z, f_p1, f_p2, annot)
	return mplot
end

#==Un-exported interface
===============================================================================
	new()
	update(::InspectDR.Multiplot, fmin, fmax, GdB, f_z, f_p1, f_p2, annot=true)
=#

end #module


text:
using Test,TensorKit,MPSKit,MPSKitModels

@testset "xxz" begin
    @testset "nonsym_xxz_ham" begin
        th = MPSKitModels.nonsym_xxz_ham()
        ts = InfiniteMPS([ℂ^3],[ℂ^48]);
        (ts,pars,_) = find_groundstate(ts,th,Vumps(maxiter=400,verbose=false));
        (energies,Bs) = excitations(th,QuasiparticleAnsatz(),Float64(pi),ts,pars);
        @test energies[1] ≈ 0.41047925 atol=1e-4
    end
end

max_stars_repo_name: maartenvd/MPSKitModels

max_issues_repo_name: maartenvd/MPSKitModels

max_forks_repo_name: maartenvd/MPSKitModels


text:
using Memento
using JSON
using HTTP
using URIs

# To use included Project.toml to set up environment before running script, uncomment thsese lines
# using Pkg
# Pkg.activate(@__DIR__)
# Pkg.instantiate()

# Set up logging
logger = Memento.config!("info"; fmt="[{date}]: {msg}")
formatter = DefaultFormatter("[{date}]: {msg}"; date_fmt_string="U dd, yyyy HH:MM:SS")

handler = DefaultHandler("gddnsjl.log")
handler.fmt = formatter
push!(logger, handler)

# Function to parse config file
function config()
    path = joinpath(dirname(@__DIR__), "config.json")
    config = JSON.parsefile(path)
    return config
end

function main()
    ip = HTTP.get("https://domains.google.com/checkip").body |> x -> String(x)
    c = config()
    url = URI(parse(URI, "https://$(c["username"]):$(c["password"])@domains.google.com/nic/update"), query = Dict{String, Any}("hostname" => c["hostname"], "myip" => ip))
    rp = HTTP.post(url)
    info(logger, "hostname: $(c["hostname"]) IP: $(ip) response: $(String(rp.body))")
end

main()

max_stars_repo_name: camilovelezr/googleddns

max_issues_repo_name: camilovelezr/googleddns

max_forks_repo_name: camilovelezr/googleddns


text:
@testset "_minvol" begin
    Y = [
        -1 -1 1 2
        1 -1 -1 2
    ]

    @testset "centered" begin
        u, R = MinimumVolumeEllipsoids._minvol(Y, 1e-10)

        @test u == [0.5, 0, 0, 0.5]
        @test R.L * R.L' ≈ [5/2 3/2; 3/2 5/2] atol = 1e-10
    end

    @testset "not centered" begin
        X = [Y; ones(1, 4)]

        u, R = MinimumVolumeEllipsoids._minvol(X, 1e-10)

        @test u ≈ [9 / 32, 4 / 32, 9 / 32, 10 / 32] atol = 1e-10
        @test X * Diagonal(u) * X' ≈ [31/16 13/16 1/2; 13/16 31/16 1/2; 1/2 1/2 1] atol =
            1e-10
    end
end


text:
using StaticArrays,JLD2

"""
The file library.gz is a gzipped ascii file.  In this file, each
line refers to a deal.  Here is a typical one (the first, in fact):

JT852.93.KQ7.J82 AQ97.JT654.T6.A5 43.AK8.A542.7643 K6.Q72.J983.KQT9:88887777A9A977778888

There are 88 characters per line.  The first 16 give West's hand, with
suits separated by periods.  Then there is a space.  The next 16
characters give North's hand, then East's, and then South's.  Then
there is a colon, followed by 20 results (in hexadecimal).  The first
four refer to NT: the number of tricks that can be taken by N/S with S
leading, then E, then N, then W.  The next four refer to S trumps,
then H trumps, then D trumps, and finally C trumps.
"""
function readdoubledummy(stream::IO)
    data = []
    for line in eachline(stream)
        @assert length(line) == 88
        (hands,results) = split(line,':')
        hands = split(hands, ' ')
        h = zeros(UInt64,4)
        for player in 1:4 # w,n,e,s
            suits = reverse(split(hands[player],'.'))
            for suit in 1:4 # c,d,h,s
                for char in collect(suits[suit])
                    value = findfirst(CARDCHAR,char)
                    @assert value > 0
                    card = makecard(suit,value)
                    h[player] |= (1 << (card - 1))
                end
            end
        end
        r = zeros(UInt8,4,5)
        for trump in 1:5 # c,d,h,s,n
            for leader in 1:4 # w,n,e,s
                r[leader,trump] = hex2int(results[5-leader+4*(5-trump)])
            end
        end
        push!(data, (SVector{4}(h),SMatrix{4,5}(r)))
    end
    return data
end

function writedoubledummy(stream::IO,data)
    for (h,r) in data
        for player in 1:4
            for suit in 4:-1:1
                for value in 13:-1:1
                    card = makecard(suit,value)
                    if (h[player] & (1 << (card-1)) != 0)
                        print(stream, CARDCHAR[value])
                    end
                end
                if suit > 1; print(stream, '.'); end
            end
            if player < 4; print(stream, ' '); end
        end
        print(stream, ':')
        for trump in 5:-1:1
            for leader in 4:-1:1
                print(stream, int2hex(r[leader,trump]))
            end
        end
        println(stream)
    end
end

writedoubledummy(data)=writedoubledummy(STDOUT,data)

#= this is buggy
"""
Reads library.dat from GIB.

The file library.dat is a binary file.

Each hand in the library is stored using 26 bytes of information.
These 26 bytes comprise 4 32-bit integers, followed by 5 16-bit
integers.  The 32-bit integers refer to the four suits, with bits 0-1
indicating who holds the ace (I take West to be 0, North to be 1, and
so on; it doesn't matter as long as it's clockwise), bits 2-3 giving
the holder of the king, and so on.  I take spades to be suit 0, but
that doesn't matter at all.

The 5 16-bit integers indicate how many tricks can be taken by each
player and in each denomination.  The first integer refers to NT, and
the subsequent ones to the suits.  Bits 0-3 give the number of tricks
that can be taken by NS with player 0 (West) on lead; bits 4-7 the
number that can be taken by NS with player 1 on lead, and so on.
"""
function readlibrary(file)
    stream = open(file)
    printrecord(readrecord(stream)...)
    close(file)
end

# WNES
# JT852.93.KQ7.J82 AQ97.JT654.T6.A5 43.AK8.A542.7643 K6.Q72.J983.KQT9:88887777A9A977778888

function readrecord(stream)
    cards = zeros(Int,4)
    hands = falses(52,4)
    @show a = read(stream, UInt32, 4)
    for suit in 1:4
        ai = a[5-suit]
        for value in 1:13
            card = makecard(value, suit)
            owner = 1 + (ai >> (2*value+4)) & 3
            hands[card,owner] = true
            cards[owner] += 1
        end
    end
    @show cards
    @assert cards == [13,13,13,13]
    b = read(stream, UInt16, 5)
    results = zeros(UInt8,4,5)
    return hands,results
end

function printrecord(hands, results)
    for player in 1:4
        for suit in 4:-1:1
            for value in 13:-1:1
                card = makecard(value,suit)
                if hands[card,player]
                    print(CARDCHAR[value])
                end
            end
            if suit > 1; print('.'); end
        end
        if player < 4; print(' '); end
    end
    println()
end

=#

max_stars_repo_name: JuliaTagBot/play

max_issues_repo_name: JuliaTagBot/play

max_forks_repo_name: JuliaTagBot/play


text:
# RUN: julia -e "import Brutus; Brutus.lit(:emit_lowered)" --startup-file=no %s 2>&1 | FileCheck %s

add(x, y) = x + y
emit(add, Float64, Float64)

# CHECK: Core.MethodMatch(Tuple{typeof(Main.Main.add), Float64, Float64}, svec(), add(x, y) in Main.Main at /{{.*}}/test/Codegen/lower/add_float64.jl:3, true)after translating to MLIR in JLIR dialect:module  {
# CHECK-NEXT:   func nested @"Tuple{typeof(Main.add), Float64, Float64}"(%arg0: !jlir<"typeof(Main.add)">, %arg1: !jlir.Float64, %arg2: !jlir.Float64) -> !jlir.Float64 attributes {llvm.emit_c_interface} {
# CHECK-NEXT:     "jlir.goto"()[^bb1] : () -> ()
# CHECK-NEXT:   ^bb1:  // pred: ^bb0
# CHECK-NEXT:     %0 = "jlir.constant"() {value = #jlir<"#<intrinsic #12 add_float>">} : () -> !jlir<"typeof(Core.IntrinsicFunction)">
# CHECK-NEXT:     %1 = "jlir.call"(%0, %arg1, %arg2) : (!jlir<"typeof(Core.IntrinsicFunction)">, !jlir.Float64, !jlir.Float64) -> !jlir.Float64
# CHECK-NEXT:     "jlir.return"(%1) : (!jlir.Float64) -> ()
# CHECK-NEXT:   }
# CHECK-NEXT: }

# CHECK: module  {
# CHECK-NEXT:   func nested @"Tuple{typeof(Main.add), Float64, Float64}"(%arg0: !jlir<"typeof(Main.add)">, %arg1: f64, %arg2: f64) -> f64 attributes {llvm.emit_c_interface} {
# CHECK-NEXT:     %0 = addf %arg1, %arg2 : f64
# CHECK-NEXT:     return %0 : f64
# CHECK-NEXT:   }
# CHECK-NEXT: }

# CHECK:   llvm.func @"Tuple{typeof(Main.add), Float64, Float64}"(%arg0: !llvm.ptr<struct<"struct_jl_value_type", opaque>>, %arg1: f64, %arg2: f64) -> f64 attributes {llvm.emit_c_interface, sym_visibility = "nested"} {
# CHECK-NEXT:     %0 = llvm.fadd %arg1, %arg2  : f64
# CHECK-NEXT:     llvm.return %0 : f64
# CHECK-NEXT:   }
# CHECK-NEXT:   llvm.func @"_mlir_ciface_Tuple{typeof(Main.add), Float64, Float64}"(%arg0: !llvm.ptr<struct<"struct_jl_value_type", opaque>>, %arg1: f64, %arg2: f64) -> f64 attributes {llvm.emit_c_interface, sym_visibility = "nested"} {
# CHECK-NEXT:     %0 = llvm.call @"Tuple{typeof(Main.add), Float64, Float64}"(%arg0, %arg1, %arg2) : (!llvm.ptr<struct<"struct_jl_value_type", opaque>>, f64, f64) -> f64
# CHECK-NEXT:     llvm.return %0 : f64
# CHECK-NEXT:   }
# CHECK-NEXT: }

max_stars_repo_name: stephenneuendorffer/brutus

max_issues_repo_name: stephenneuendorffer/brutus

max_forks_repo_name: stephenneuendorffer/brutus


text:
export ThreadedModel

"""
    ThreadedModel
Experimental model wrapper than enables multi-threading. This will likely
be depreciated in the future for moving threading to executors similar to FLoops.jl
"""
struct ThreadedModel{M} <: AbstractModel
    model::M
end
@inline imanalytic(::Type{ThreadedModel{M}}) where {M} = imanalytic(M)
@inline visanalytic(::Type{ThreadedModel{M}}) where {M} = visanalytic(M)
@inline isprimitive(::Type{ThreadedModel{M}}) where {M} = isprimitive(M)

@inline visibility_point(m::ThreadedModel, u, v) = visibility_point(m.model, u, v)
@inline intensity_point(m::ThreadedModel, x, y) = intensity_point(m.model, x, y)

function intensitymap(::IsAnalytic, s::ThreadedModel, fovx::Number, fovy::Number, nx::Int, ny::Int; executor= pulse=ComradeBase.DeltaPulse())
    T = typeof(intensity_point(s, 0.0, 0.0))
    img = IntensityMap(zeros(T, ny, nx), fovx, fovy, pulse)
    intensitymap!(IsAnalytic(), img, s)
    return img
end

function intensitymap!(::IsAnalytic, img::AbstractIntensityMap, s::ThreadedModel)
    x,y = imagepixels(img)
    dx, dy = pixelsizes(img)
    Threads.@threads for I in CartesianIndices(img)
        iy,ix = Tuple(I)
        img[I] = intensity_point(s, x[ix], y[iy])*dx*dy
    end
    return img
end