text:
export objective, smooth_objective

# NOTE: RobustLoss are not always everywhere  smooth but "smooth-enough".
const SmoothLoss = Union{L2Loss, LogisticLoss, MultinomialLoss, RobustLoss}

"""
$SIGNATURES

Return the objective function (sum of loss + penalty) of a Generalized Linear Model.
"""
objective(glr::GLR, n) = glr.loss + glr.penalty * ifelse(glr.scale_penalty_with_samples, n, 1.)

"""
$SIGNATURES

Return a function computing the objective at a given point `θ`.
Note that the [`apply_X`](@ref) takes care of a potential intercept.
"""
objective(glr::GLR, X, y; c::Int=0) =
    θ -> objective(glr, size(X, 1))(y, apply_X(X, θ, c), view_θ(glr, θ))

"""
$SIGNATURES

Return a function computing the smooth part of the objective at a given
evaluation point `θ`.
"""
smooth_objective(glr::GLR, X, y; c::Int=0) =
    θ -> smooth_objective(glr, size(X, 1))(y, apply_X(X, θ, c), view_θ(glr, θ))

"""
$SIGNATURES

Return the smooth part of the objective function of a GLR.
"""
smooth_objective(glr::GLR{<:SmoothLoss,<:ENR}, n) = glr.loss + get_l2(glr.penalty) * ifelse(glr.scale_penalty_with_samples, n, 1.)

smooth_objective(::GLR) = @error "Case not implemented yet."

"""
$SIGNATURES

Return a model corresponding to the smooth part of the objective.
"""
get_smooth(glr::GLR) =
    GLR(glr.loss, get_l2(glr.penalty),
        glr.fit_intercept, glr.penalize_intercept, glr.scale_penalty_with_samples)

"""
$SIGNATURES

Helper function to compute the residuals.
"""
function get_residuals!(r, X, θ, y)
    apply_X!(r, X, θ)
    r .-= y
end


text:
using Documenter
using Xfoil, AirfoilPrep, GeometricTools, CCBlade, FLOWVLM, FLOWUnsteady, FLOWVPM
vlm = FLOWVLM
gt = GeometricTools
vpm = FLOWVPM
# using FLOWUnsteady

makedocs(
    sitename = "FLOWUnsteady",
    format = Documenter.HTML(),
    modules = [Xfoil, AirfoilPrep, GeometricTools, CCBlade, FLOWVLM, FLOWUnsteady, FLOWVPM, vlm, vpm],
    # modules = [FLOWUnsteady],
    pages = ["Home" => "index.md",
            "Tutorials" => Any[
                "tutorials/installation-instructions.md",
                "tutorials/getting-started.md",
                "tutorials/aeroacoustics.md",
                "tutorials/colab.md",
                ],
            "How to " => Any[
                "how-to-guide/paraview-visualization.md",
                "how-to-guide/define-complex-wings.md",
                "how-to-guide/define-complex-systems.md",
                "how-to-guide/define-a-rotor.md",
                "how-to-guide/define-kinematic-maneuvers.md",
                "how-to-guide/set-up-runtime-functions.md",
                "how-to-guide/select-a-solver.md",
                "how-to-guide/wake-stability.md",],
            "Reference" => Any[
                "reference/FLOWVLMfunctions.md",
                "reference/VPMfunctions.md",
                "reference/FLOWUnsteadyfunctions.md",],

            "Theory" => Any[
                "theory/FLOWVLMtheory.md",
                "theory/VPMtheory.md",
                "theory/FLOWUnsteadytheory.md",],
            ]
)

# Documenter can also automatically deploy documentation to gh-pages.
# See "Hosting Documentation" and deploydocs() in the Documenter manual
# for more information.
deploydocs(
    repo   = "github.com/byuflowlab/FLOWUnsteady.git",
)

ext: jl


text:
function rebuild_conv(x, shape)
    w = Array{Float32}(undef, shape...)
    filter_x, filter_y = shape[1:2] .+ 1
    @inbounds for (i, j, k, m) in Iterators.product([1:s for (i, s) in enumerate(shape)]...)
        w[filter_x - i, filter_y - j, k, m] = Float32(x[m][k][j][i])
    end
    w
end

function rebuild_dense(x, shape)
    w = Array{Float32}(undef, shape...)
    @inbounds for (i, j) in Iterators.product([1:s for s in shape]...)
        w[i, j] = Float32(x[i][j])
    end
    w
end

function loadweights!(model, weights)
    for (p, w) in zip(params(model), weights)
        if length(p) == 1
            continue
        end
        copyto!(p, w)
    end
end

function load_conv(conv, paths, pid)
    path = paths[pid]
    @load path p
    copyto!(conv.weight, rebuild_conv(p, size(conv.weight)))
    pid += 1
    pid
end

function load_bn(bn, paths, pid)
    path = paths[pid]
    @load path p
    copyto!(bn.γ, p)
    pid += 1
    path = paths[pid]
    @load path p
    copyto!(bn.β, p)
    pid += 1
    path = paths[pid]
    @load path p
    copyto!(bn.μ, p)
    pid += 1
    path = paths[pid]
    @load path p
    copyto!(bn.σ², p)
    pid += 1
    pid
end

function load_dense(d, paths, pid)
    path = paths[pid]
    @load path p
    copyto!(d.W, rebuild_dense(p, size(d.W)))
    pid += 1
    path = paths[pid]
    @load path p
    copyto!(d.b, p)
    pid += 1
    pid
end

function load_entry(entry, paths, pid)
    if entry isa Conv
        pid = load_conv(entry, paths, pid)
    elseif entry isa BatchNorm
        pid = load_bn(entry, paths, pid)
    elseif entry isa Dense
        pid = load_dense(entry, paths, pid)
    elseif entry isa Chain
        for e in entry
            pid = load_entry(e, paths, pid)
        end
    elseif entry isa SkipConnection
        pid = load_entry(entry.layers, paths, pid)
        pid = load_entry(entry.connection, paths, pid)
    elseif entry isa Shortcut
        pid = load_entry(entry.shortcut, paths, pid)
    end
    pid
end

function load()
    path = "resnet18-pretrained.bson"
    dir = raw"C:\Users\tonys\projects\julia\ResNet\resnet18-pretrained"
    weights_paths = readdir(dir, join=true)

    pid = 1
    model = resnet(18, 1000)
    # @info "Loading model"
    # pid = load_entry(model.entry, weights_paths, pid)
    # pid = load_entry(model.encoder, weights_paths, pid)
    # pid = load_entry(model.head, weights_paths, pid)
    # @info "Loaded"

    # parameters = params(model) .|> cpu
    # @save path parameters
    @load path parameters
    loadweights!(model, parameters)
    testmode!(model)

    images = [
        raw"C:\Users\tonys\Downloads\kazan.jpg",
        raw"C:\Users\tonys\Downloads\el.jpg",
        raw"C:\Users\tonys\Downloads\pug.jpg",
    ]
    for path in images
        image = FileIO.load(path)
        image = Float32.(Images.channelview(image))
        image = permutedims(image, (3, 2, 1))
        image = Flux.unsqueeze(image, 4)

        entry = image |> model.entry
        features = [model.encoder[1](entry)]
        for e in model.encoder[2:end]
            push!(features, e(features[end]))
        end
        y = Flux.softmax(model.head(features[end])[:, 1])

        # y = Flux.softmax(model(image)[:, 1])
        top5 = sortperm(y)[end - 4:end]
        println(path)
        println(top5 .- 1)
        println(y[top5])
        println("===========")
    end
    nothing
end


text:
# Compares EM variants on given network structure
# Example of commaind line call:
#     julia --threads=4 --color=yes em.jl assets/nltcs.spn assets/nltcs.train.csv assets/nltcs.valid.csv 100
#using DataFrames, CSV
using DelimitedFiles # much faster than above
using Random
using SPNetworks
import SPNetworks: NLL
import SPNetworks.ParameterLearning: SEM, initialize, converged, update

if length(ARGS) < 3
    println("Usage: julia --color=yes em.jl spn_filename train_filename validation_filename [batchsize]")
    exit()
end
spn_filename = normpath(ARGS[1])
train_filename = normpath(ARGS[2])
valid_filename = normpath(ARGS[3])
batchsize = 100 # default
if length(ARGS) >= 3
    batchsize = parse(Int, ARGS[4])
end
println("Maximum batchsize: $batchsize")
println("Threads: ", Threads.nthreads())

println("SPN: $spn_filename")
spn = SumProductNetwork(spn_filename)
@show summary(spn)
#SPNetworks.save(spn,normpath("$train_filename"))
nvars = length(scope(spn))
println("Training data: $train_filename")
# Load datasets
@time tdata = readdlm(train_filename, ',', Float64) .+ 1;
@show summary(tdata)
println("Validation data: $valid_filename")
@time vdata = readdlm(valid_filename, ',', Float64) .+ 1;
@show summary(vdata)
# initialize Stochastic EM learner
learner = SEM(spn)
Random.seed!(3)
#initialize(learner) # generate random weights for each sum node
avgnll = NLL(spn, tdata)
runnll = 0.0
println("It: $(learner.steps) \t train NLL: $avgnll \t held-out NLL: $(NLL(spn, vdata))")
indices = shuffle!(Vector(1:size(tdata,1)))
# Running Expectation Maximization
while !converged(learner) && learner.steps < 31
    global avgnll, runnll
    sid = rand(1:(length(indices)-batchsize))
    batch = view(tdata, indices[sid:(sid+batchsize-1)], :) # extract minibatch sample
    # if learner.steps % 2 == 0
        η = max(0.95^learner.steps, 0.3) # learning rate
        @time update(learner, batch, η)
        testnll = NLL(spn, vdata)
        batchnll = NLL(spn, batch)
        # running average NLL
        avgnll *= (learner.steps-1)/learner.steps # discards initial NLL
        avgnll += batchnll/learner.steps
        runnll = (1-η)*runnll + η*batchnll
        println("It: $(learner.steps) \t avg NLL: $avgnll \t mov NLL: $runnll \t batch NLL: $batchnll \t held-out NLL: $testnll \t Learning rate: $η")
    # else
    #     @time update(learner, batch, 0.9^learner.steps)
    #     println("It: $(learner.steps) \t NLL: $(learner.score)")
    # end
end


text:
#########
# Trial #
#########

mutable struct Trial
    params::Parameters
    times::Vector{Float64}
    gctimes::Vector{Float64}
    memory::Int
    allocs::Int
end

Trial(params::Parameters) = Trial(params, Float64[], Float64[], typemax(Int), typemax(Int))

function Base.:(==)(a::Trial, b::Trial)
    return a.params == b.params &&
           a.times == b.times &&
           a.gctimes == b.gctimes &&
           a.memory == b.memory &&
           a.allocs == b.allocs
end

Base.copy(t::Trial) = Trial(copy(t.params), copy(t.times), copy(t.gctimes), t.memory, t.allocs)

function Base.push!(t::Trial, time, gctime, memory, allocs)
    push!(t.times, time)
    push!(t.gctimes, gctime)
    memory < t.memory && (t.memory = memory)
    allocs < t.allocs && (t.allocs = allocs)
    return t
end

function Base.deleteat!(t::Trial, i)
    deleteat!(t.times, i)
    deleteat!(t.gctimes, i)
    return t
end

Base.length(t::Trial) = length(t.times)
Base.getindex(t::Trial, i::Number) = push!(Trial(t.params), t.times[i], t.gctimes[i], t.memory, t.allocs)
Base.getindex(t::Trial, i) = Trial(t.params, t.times[i], t.gctimes[i], t.memory, t.allocs)
Base.lastindex(t::Trial) = length(t)

function Base.sort!(t::Trial)
    inds = sortperm(t.times)
    t.times = t.times[inds]
    t.gctimes = t.gctimes[inds]
    return t
end

Base.sort(t::Trial) = sort!(copy(t))

Base.time(t::Trial) = time(minimum(t))
gctime(t::Trial) = gctime(minimum(t))
memory(t::Trial) = t.memory
allocs(t::Trial) = t.allocs
params(t::Trial) = t.params

# returns the index of the first outlier in `values`, if any outliers are detected.
# `values` is assumed to be sorted from least to greatest, and assumed to be right-skewed.
function skewcutoff(values)
    current_values = copy(values)
    while mean(current_values) > median(current_values)
        deleteat!(current_values, length(current_values))
    end
    return length(current_values) + 1
end

skewcutoff(t::Trial) = skewcutoff(t.times)

function rmskew!(t::Trial)
    sort!(t)
    i = skewcutoff(t)
    i <= length(t) && deleteat!(t, i:length(t))
    return t
end

function rmskew(t::Trial)
    st = sort(t)
    return st[1:(skewcutoff(st) - 1)]
end

trim(t::Trial, percentage = 0.1) = t[1:max(1, floor(Int, length(t) - (length(t) * percentage)))]

#################
# TrialEstimate #
#################

mutable struct TrialEstimate
    params::Parameters
    time::Float64
    gctime::Float64
    memory::Int
    allocs::Int
end

function TrialEstimate(trial::Trial, t, gct)
    return TrialEstimate(params(trial), t, gct, memory(trial), allocs(trial))
end

function Base.:(==)(a::TrialEstimate, b::TrialEstimate)
    return a.params == b.params &&
           a.time == b.time &&
           a.gctime == b.gctime &&
           a.memory == b.memory &&
           a.allocs == b.allocs
end

Base.copy(t::TrialEstimate) = TrialEstimate(copy(t.params), t.time, t.gctime, t.memory, t.allocs)

function Base.minimum(trial::Trial)
    i = argmin(trial.times)
    return TrialEstimate(trial, trial.times[i], trial.gctimes[i])
end

function Base.maximum(trial::Trial)
    i = argmax(trial.times)
    return TrialEstimate(trial, trial.times[i], trial.gctimes[i])
end

Statistics.median(trial::Trial) = TrialEstimate(trial, median(trial.times), median(trial.gctimes))
Statistics.mean(trial::Trial) = TrialEstimate(trial, mean(trial.times), mean(trial.gctimes))
Statistics.std(trial::Trial) = TrialEstimate(trial, std(trial.times), std(trial.gctimes))

Base.isless(a::TrialEstimate, b::TrialEstimate) = isless(time(a), time(b))

Base.time(t::TrialEstimate) = t.time
gctime(t::TrialEstimate) = t.gctime
memory(t::TrialEstimate) = t.memory
allocs(t::TrialEstimate) = t.allocs
params(t::TrialEstimate) = t.params

##############
# TrialRatio #
##############

mutable struct TrialRatio
    params::Parameters
    time::Float64
    gctime::Float64
    memory::Float64
    allocs::Float64
end

function Base.:(==)(a::TrialRatio, b::TrialRatio)
    return a.params == b.params &&
           a.time == b.time &&
           a.gctime == b.gctime &&
           a.memory == b.memory &&
           a.allocs == b.allocs
end

Base.copy(t::TrialRatio) = TrialRatio(copy(t.params), t.time, t.gctime, t.memory, t.allocs)

Base.time(t::TrialRatio) = t.time
gctime(t::TrialRatio) = t.gctime
memory(t::TrialRatio) = t.memory
allocs(t::TrialRatio) = t.allocs
params(t::TrialRatio) = t.params

function ratio(a::Real, b::Real)
    if a == b # so that ratio(0.0, 0.0) returns 1.0
        return one(Float64)
    end
    return Float64(a / b)
end

function ratio(a::TrialEstimate, b::TrialEstimate)
    ttol = max(params(a).time_tolerance, params(b).time_tolerance)
    mtol = max(params(a).memory_tolerance, params(b).memory_tolerance)
    p = Parameters(params(a); time_tolerance = ttol, memory_tolerance = mtol)
    return TrialRatio(p, ratio(time(a), time(b)), ratio(gctime(a), gctime(b)),
                      ratio(memory(a), memory(b)), ratio(allocs(a), allocs(b)))
end

gcratio(t::TrialEstimate) =  ratio(gctime(t), time(t))

##################
# TrialJudgement #
##################

struct TrialJudgement
    ratio::TrialRatio
    time::Symbol
    memory::Symbol
end

function TrialJudgement(r::TrialRatio)
    ttol = params(r).time_tolerance
    mtol = params(r).memory_tolerance
    return TrialJudgement(r, judge(time(r), ttol), judge(memory(r), mtol))
end

function Base.:(==)(a::TrialJudgement, b::TrialJudgement)
    return a.ratio == b.ratio &&
           a.time == b.time &&
           a.memory == b.memory
end

Base.copy(t::TrialJudgement) = TrialJudgement(copy(t.params), t.time, t.memory)

Base.time(t::TrialJudgement) = t.time
memory(t::TrialJudgement) = t.memory
ratio(t::TrialJudgement) = t.ratio
params(t::TrialJudgement) = params(ratio(t))

judge(a::TrialEstimate, b::TrialEstimate; kwargs...) = judge(ratio(a, b); kwargs...)

function judge(r::TrialRatio; kwargs...)
    newr = copy(r)
    newr.params = Parameters(params(r); kwargs...)
    return TrialJudgement(newr)
end

function judge(ratio::Real, tolerance::Float64)
    if isnan(ratio) || (ratio - tolerance) > 1.0
        return :regression
    elseif (ratio + tolerance) < 1.0
        return :improvement
    else
        return :invariant
    end
end

isimprovement(f, t::TrialJudgement) = f(t) == :improvement
isimprovement(t::TrialJudgement) = isimprovement(time, t) || isimprovement(memory, t)

isregression(f, t::TrialJudgement) = f(t) == :regression
isregression(t::TrialJudgement) = isregression(time, t) || isregression(memory, t)

isinvariant(f, t::TrialJudgement) = f(t) == :invariant
isinvariant(t::TrialJudgement) = isinvariant(time, t) && isinvariant(memory, t)

const colormap = (
    regression = :red,
    improvement = :green,
    invariant = :normal,
)

printtimejudge(io, t::TrialJudgement) =
    printstyled(io, time(t); color=colormap[time(t)])
printmemoryjudge(io, t::TrialJudgement) =
    printstyled(io, memory(t); color=colormap[memory(t)])

###################
# Pretty Printing #
###################

prettypercent(p) = string(@sprintf("%.2f", p * 100), "%")

function prettydiff(p)
    diff = p - 1.0
    return string(diff >= 0.0 ? "+" : "", @sprintf("%.2f", diff * 100), "%")
end

function prettytime(t)
    if t < 1e3
        value, units = t, "ns"
    elseif t < 1e6
        value, units = t / 1e3, "μs"
    elseif t < 1e9
        value, units = t / 1e6, "ms"
    else
        value, units = t / 1e9, "s"
    end
    return string(@sprintf("%.3f", value), " ", units)
end

function prettymemory(b)
    if b < 1024
        return string(b, " bytes")
    elseif b < 1024^2
        value, units = b / 1024, "KiB"
    elseif b < 1024^3
        value, units = b / 1024^2, "MiB"
    else
        value, units = b / 1024^3, "GiB"
    end
    return string(@sprintf("%.2f", value), " ", units)
end

function withtypename(f, io, t)
    needtype = get(io, :typeinfo, Nothing) !== typeof(t)
    if needtype
        print(io, nameof(typeof(t)), '(')
    end
    f()
    if needtype
        print(io, ')')
    end
end

function bindata(sorteddata, nbins, min, max)
    Δ = (max - min) / nbins
    bins = zeros(nbins)
    lastpos = 0
    for i ∈ 1:nbins
        pos = searchsortedlast(sorteddata, min + i * Δ)
        bins[i] = pos - lastpos
        lastpos = pos
    end
    bins
end
bindata(sorteddata, nbins) = bindata(sorteddata, nbins, first(sorteddata), last(sorteddata))

function asciihist(bins, height=1)
    histbars = ['▁', '▂', '▃', '▄', '▅', '▆', '▇', '█']
    if minimum(bins) == 0
        barheights = 2 .+ round.(Int, (height * length(histbars) - 2) * bins ./ maximum(bins))
        barheights[bins .== 0] .= 1
    else
        barheights = 1 .+ round.(Int, (height * length(histbars) - 1) * bins ./ maximum(bins))
    end
    heightmatrix = [min(length(histbars), barheights[b] - (h-1) * length(histbars))
                    for h in height:-1:1, b in 1:length(bins)]
    map(height -> if height < 1; ' ' else histbars[height] end, heightmatrix)
end

_summary(io, t, args...) = withtypename(() -> print(io, args...), io, t)

Base.summary(io::IO, t::Trial) = _summary(io, t, prettytime(time(t)))
Base.summary(io::IO, t::TrialEstimate) = _summary(io, t, prettytime(time(t)))
Base.summary(io::IO, t::TrialRatio) = _summary(io, t, prettypercent(time(t)))
Base.summary(io::IO, t::TrialJudgement) = withtypename(io, t) do
    print(io, prettydiff(time(ratio(t))), " => ")
    printtimejudge(io, t)
end

_show(io, t) =
    if get(io, :compact, true)
        summary(io, t)
    else
        show(io, MIME"text/plain"(), t)
    end

Base.show(io::IO, t::Trial) = _show(io, t)
Base.show(io::IO, t::TrialEstimate) = _show(io, t)
Base.show(io::IO, t::TrialRatio) = _show(io, t)
Base.show(io::IO, t::TrialJudgement) = _show(io, t)

function Base.show(io::IO, ::MIME"text/plain", t::Trial)

    pad = get(io, :pad, "")
    print(io, "BenchmarkTools.Trial: ", length(t), " sample", if length(t) > 1 "s" else "" end,
          " with ", t.params.evals, " evaluation", if t.params.evals > 1 "s" else "" end ,".\n")

    perm = sortperm(t.times)
    times = t.times[perm]
    gctimes = t.gctimes[perm]

    if length(t) > 1
        med = median(t)
        avg = mean(t)
        std = Statistics.std(t)
        min = minimum(t)
        max = maximum(t)

        medtime, medgc = prettytime(time(med)), prettypercent(gcratio(med))
        avgtime, avggc = prettytime(time(avg)), prettypercent(gcratio(avg))
        stdtime, stdgc = prettytime(time(std)), prettypercent(Statistics.std(gctimes ./ times))
        mintime, mingc = prettytime(time(min)), prettypercent(gcratio(min))
        maxtime, maxgc = prettytime(time(max)), prettypercent(gcratio(max))

        memorystr = string(prettymemory(memory(min)))
        allocsstr = string(allocs(min))
    elseif length(t) == 1
        print(io, pad, " Single result which took ")
        printstyled(io, prettytime(times[1]); color=:blue)
        print(io, " (", prettypercent(gctimes[1]/times[1]), " GC) ")
        print(io, "to evaluate,\n")
        print(io, pad, " with a memory estimate of ")
        printstyled(io, prettymemory(t.memory[1]); color=:yellow)
        print(io, ", over ")
        printstyled(io, t.allocs[1]; color=:yellow)
        print(io, " allocations.")
        return
    else
        print(io, pad, " No results.")
        return
    end

    lmaxtimewidth = maximum(length.((medtime, avgtime, mintime)))
    rmaxtimewidth = maximum(length.((stdtime, maxtime)))
    lmaxgcwidth = maximum(length.((medgc, avggc, mingc)))
    rmaxgcwidth = maximum(length.((stdgc, maxgc)))

    # Main stats

    print(io, pad, " Range ")
    printstyled(io, "("; color=:light_black)
    printstyled(io, "min"; color=:cyan, bold=true)
    print(io, " … ")
    printstyled(io, "max"; color=:magenta)
    printstyled(io, "):  "; color=:light_black)
    printstyled(io, lpad(mintime, lmaxtimewidth); color=:cyan, bold=true)
    print(io, " … ")
    printstyled(io, lpad(maxtime, rmaxtimewidth); color=:magenta)
    print(io, "  ")
    printstyled(io, "┊"; color=:light_black)
    print(io, " GC ")
    printstyled(io, "("; color=:light_black)
    print(io, "min … max")
    printstyled(io, "): "; color=:light_black)
    print(io, lpad(mingc, lmaxgcwidth), " … ", lpad(maxgc, rmaxgcwidth))

    print(io, "\n", pad, " Time  ")
    printstyled(io, "("; color=:light_black)
    printstyled(io, "median"; color=:blue, bold=true)
    printstyled(io, "):     "; color=:light_black)
    printstyled(io, lpad(medtime, lmaxtimewidth), rpad(" ", rmaxtimewidth + 5); color=:blue, bold=true)
    printstyled(io, "┊"; color=:light_black)
    print(io, " GC ")
    printstyled(io, "("; color=:light_black)
    print(io, "median")
    printstyled(io, "):    "; color=:light_black)
    print(io, lpad(medgc, lmaxgcwidth))

    print(io, "\n", pad, " Time  ")
    printstyled(io, "("; color=:light_black)
    printstyled(io, "mean"; color=:green, bold=true)
    print(io, " ± ")
    printstyled(io, "σ"; color=:green)
    printstyled(io, "):   "; color=:light_black)
    printstyled(io, lpad(avgtime, lmaxtimewidth); color=:green, bold=true)
    print(io, " ± ")
    printstyled(io, lpad(stdtime, rmaxtimewidth); color=:green)
    print(io, "  ")
    printstyled(io, "┊"; color=:light_black)
    print(io, " GC ")
    printstyled(io, "("; color=:light_black)
    print(io, "mean ± σ")
    printstyled(io, "):  "; color=:light_black)
    print(io, lpad(avggc, lmaxgcwidth), " ± ", lpad(stdgc, rmaxgcwidth))

    # Histogram

    histquantile = 0.99
    # The height and width of the printed histogram in characters.
    histheight = 2
    histwidth = 42 + lmaxtimewidth + rmaxtimewidth

    histtimes = times[1:round(Int, histquantile*end)]
    histmin = get(io, :histmin, first(histtimes))
    histmax = get(io, :histmax, last(histtimes))
    logbins = get(io, :logbins, nothing)
    bins = bindata(histtimes, histwidth - 1, histmin, histmax)
    append!(bins, [1, floor((1-histquantile) * length(times))])
    # if median size of (bins with >10% average data/bin) is less than 5% of max bin size, log the bin sizes
    if logbins === true || (logbins === nothing && median(filter(b -> b > 0.1 * length(times) / histwidth, bins)) / maximum(bins) < 0.05)
        bins, logbins = log.(1 .+ bins), true
    else
        logbins = false
    end
    hist = asciihist(bins, histheight)
    hist[:,end-1] .= ' '
    maxbin = maximum(bins)

    delta1 = (histmax - histmin) / (histwidth - 1)
    if delta1 > 0
        medpos = 1 + round(Int, (histtimes[length(times) ÷ 2] - histmin) / delta1)
        avgpos = 1 + round(Int, (mean(times) - histmin) / delta1)
    else
        medpos, avgpos = 1, 1
    end

    print(io, "\n")
    for r in axes(hist, 1)
        print(io, "\n", pad, "  ")
        for (i, bar) in enumerate(view(hist, r, :))
            color = :default
            if i == avgpos color = :green end
            if i == medpos color = :blue end
            printstyled(io, bar; color=color)
        end
    end

    remtrailingzeros(timestr) = replace(timestr, r"\.?0+ " => " ")
    minhisttime, maxhisttime = remtrailingzeros.(prettytime.(round.([histmin, histmax], sigdigits=3)))

    print(io, "\n", pad, "  ", minhisttime)
    caption = "Histogram: " * ( logbins ? "log(frequency)" : "frequency" ) * " by time"
    if logbins
        printstyled(io, " " ^ ((histwidth - length(caption)) ÷ 2 - length(minhisttime)); color=:light_black)
        printstyled(io, "Histogram: "; color=:light_black)
        printstyled(io, "log("; bold=true, color=:light_black)
        printstyled(io, "frequency"; color=:light_black)
        printstyled(io, ")"; bold=true, color=:light_black)
        printstyled(io, " by time"; color=:light_black)
    else
        printstyled(io, " " ^ ((histwidth - length(caption)) ÷ 2 - length(minhisttime)), caption; color=:light_black)
    end
    print(io, lpad(maxhisttime, ceil(Int, (histwidth - length(caption)) / 2) - 1), " ")
    printstyled(io, "<"; bold=true)

    # Memory info

    print(io, "\n\n", pad, " Memory estimate")
    printstyled(io, ": "; color=:light_black)
    printstyled(io, memorystr; color=:yellow)
    print(io, ", allocs estimate")
    printstyled(io, ": "; color=:light_black)
    printstyled(io, allocsstr; color=:yellow)
    print(io, ".")
end

function Base.show(io::IO, ::MIME"text/plain", t::TrialEstimate)
    println(io, "BenchmarkTools.TrialEstimate: ")
    pad = get(io, :pad, "")
    println(io, pad, "  time:             ", prettytime(time(t)))
    println(io, pad, "  gctime:           ", prettytime(gctime(t)), " (", prettypercent(gctime(t) / time(t)),")")
    println(io, pad, "  memory:           ", prettymemory(memory(t)))
    print(io,   pad, "  allocs:           ", allocs(t))
end

function Base.show(io::IO, ::MIME"text/plain", t::TrialRatio)
    println(io, "BenchmarkTools.TrialRatio: ")
    pad = get(io, :pad, "")
    println(io, pad, "  time:             ", time(t))
    println(io, pad, "  gctime:           ", gctime(t))
    println(io, pad, "  memory:           ", memory(t))
    print(io,   pad, "  allocs:           ", allocs(t))
end

function Base.show(io::IO, ::MIME"text/plain", t::TrialJudgement)
    println(io, "BenchmarkTools.TrialJudgement: ")
    pad = get(io, :pad, "")
    print(io, pad, "  time:   ", prettydiff(time(ratio(t))), " => ")
    printtimejudge(io, t)
    println(io, " (", prettypercent(params(t).time_tolerance), " tolerance)")
    print(io,   pad, "  memory: ", prettydiff(memory(ratio(t))), " => ")
    printmemoryjudge(io, t)
    println(io, " (", prettypercent(params(t).memory_tolerance), " tolerance)")
end


text:
# Tests for Generalized DF DEA Models
@testset "GeneralizedDFDEAModel" begin

    ## Test Generalized DF DEA Model with Zofío and Prieto (2006) data
    X = [5 3; 2 4; 4 2; 4 8; 7 9]
    Y = [7 4; 10 8; 8 10; 5 4; 3 6]

    # alpha = 0.5 CRS equals Input Oriented CRS
    deaio = dea(X, Y, orient = :Input, rts = :CRS)
    deagdf05crs = deagdf(X, Y, alpha = 0.5, rts = :CRS)

    @test typeof(deagdf05crs) == GeneralizedDFDEAModel

    @test nobs(deagdf05crs) == 5
    @test ninputs(deagdf05crs) == 2
    @test noutputs(deagdf05crs) == 2
    @test efficiency(deaio) ≈ efficiency(deagdf05crs) atol = 1e-7

    @test efficiency(deagdf(targets(deagdf05crs, :X), targets(deagdf05crs, :Y), alpha = 0.5, rts = :CRS, slack = false)) ≈ ones(5)
    @test efficiency(deaadd(targets(deagdf05crs, :X), targets(deagdf05crs, :Y))) ≈ zeros(5) atol=1e-5

    # alphpa = 0.5 VRS
    deagdf05vrs = deagdf(X, Y, alpha = 0.5, rts = :VRS)

    @test nobs(deagdf05vrs) == 5
    @test ninputs(deagdf05vrs) == 2
    @test noutputs(deagdf05vrs) == 2
    @test efficiency(deagdf05vrs) ≈ [0.682;
                               1;
                               1;
                               0.250;
                               0.360] atol = 1e-3
    @test slacks(deagdf05vrs, :X) ≈ [0.605935 0;
                                     0 0;
                                     0 0;
                                     0 0;
                                     0.2 3.4] atol = 1e-5

    @test slacks(deagdf05vrs, :Y) ≈ [0 4.67865;
                                     0 0;
                                     0 0;
                                     0 0;
                                     3.0 0] atol = 1e-5

    @test efficiency(deagdf(targets(deagdf05vrs, :X), targets(deagdf05vrs, :Y), alpha = 0.5, rts = :CRS, slack = false)) ≈ ones(5)
    @test efficiency(deaadd(targets(deagdf05vrs, :X), targets(deagdf05vrs, :Y))) ≈ zeros(5) atol=1e-5

    # Test no slacks
    deagdfnoslack = deagdf(X, Y, alpha = 0.5, rts = :VRS, slack = false)
    @test efficiency(deagdfnoslack) == efficiency(deagdf05vrs)
    @test isempty(slacks(deagdfnoslack, :X)) == 1
    @test isempty(slacks(deagdfnoslack, :Y)) == 1

    @test efficiency(deagdf(targets(deagdfnoslack, :X), targets(deagdfnoslack, :Y), alpha = 0.5, rts = :VRS, slack = false)) ≈ ones(5)
    @test efficiency(deaadd(targets(deagdfnoslack, :X), targets(deagdfnoslack, :Y))) != zeros(5) # Different as there is no slacks in first model

    # Test default alpha is 0.5
    @test efficiency(deagdf(X, Y, rts = :CRS)) == efficiency(deagdf(X, Y, alpha = 0.5, rts = :CRS))

    ## Test if one-by-one DEA using evaluation and reference sets match initial results
    deagdf05crs_ref_eff = zeros(size(X, 1))

    deagdf05vrs_ref_eff = zeros(size(X, 1))

    Xref = X[:,:]
    Yref = Y[:,:]

    for i = 1:size(X, 1)
        Xeval = X[i:i,:]
        Xeval = Xeval[:,:]
        Yeval = Y[i:i,:]
        Yeval = Yeval[:,:]

        deagdf05crs_ref_eff[i] = efficiency(deagdf(Xeval, Yeval, alpha = 0.5, rts = :CRS, Xref = Xref, Yref = Yref, slack = false))[1]

        deagdf05vrs_ref_eff[i] = efficiency(deagdf(Xeval, Yeval, alpha = 0.5, rts = :VRS, Xref = Xref, Yref = Yref, slack = false))[1]
    end

    @test deagdf05crs_ref_eff ≈ efficiency(deagdf05crs)
    @test deagdf05vrs_ref_eff ≈ efficiency(deagdf05vrs)

    # Test with another data
    X = [2; 4; 8; 12; 6; 14; 14; 9.412];
    Y = [1; 5; 8; 9; 3; 7; 9; 2.353];

    gdf2 = deagdf(X, Y, alpha = 0.7, rts = :VRS)

    @test efficiency(gdf2) ≈ [1; 1; 1; 1; 0.423216; 0.69836; 1; 0.2315] atol = 1e-4
    @test slacks(gdf2, :X) ≈ [0; 0; 0; 0; 0; 0.570479; 2.0; 0]  atol = 1e-5
    @test slacks(gdf2, :Y) ≈ [0; 0; 0; 0; 0; 0; 0; 0] atol = 1e-6

    # Print
    show(IOBuffer(), deagdf05crs)
    show(IOBuffer(), deagdfnoslack)

    # Test errors
    @test_throws DimensionMismatch deagdf([1; 2 ; 3], [4 ; 5], alpha = 0.5) #  Different number of observations
    @test_throws DimensionMismatch deagdf([1; 2], [4 ; 5], alpha = 0.5, Xref = [1; 2; 3; 4]) # Different number of observations in reference sets
    @test_throws DimensionMismatch deagdf([1 1; 2 2], [4 4; 5 5], alpha = 0.5, Xref = [1 1 1; 2 2 2]) # Different number of inputs
    @test_throws DimensionMismatch deagdf([1 1; 2 2], [4 4; 5 5], alpha = 0.5, Yref = [4 4 4; 5 5 5]) # Different number of inputs
    @test_throws ArgumentError deagdf([1; 2; 3], [4; 5; 6], alpha = 0.5, rts = :Error) # Invalid returns to scale

    # ------------------
    # Test Vector and Matrix inputs and outputs
    # ------------------

    # Inputs is Matrix, Outputs is Vector
    X = [2 2; 1 4; 4 1; 4 3; 5 5; 6 1; 2 5; 1.6	8]
    Y = [1; 1; 1; 1; 1; 1; 1; 1]

    @test efficiency(deagdf(X, Y, rts = :VRS, slack = false)) ≈ [1; 1; 1; 1; 1; 1; 1; 1] atol = 1e-5

    # Inputs is Vector, Output is Matrix
    X = [1; 1; 1; 1; 1; 1; 1; 1]
    Y = [7 7; 4 8; 8 4; 3 5; 3 3; 8 2; 6 4; 1.5 5]

    @test efficiency(deagdf(X, Y, rts = :VRS, slack = false)) ≈ [1; 1; 1; 1; 1; 1; 1; 1] atol = 1e-5

    # Inputs is Vector, Output is Vector
    X = [2; 4; 8; 12; 6; 14; 14; 9.412]
    Y = [1; 5; 8; 9; 3; 7; 9; 2.353]

    @test efficiency(deagdf(X, Y, rts = :VRS)) ≈ [1; 1; 1; 1; 0.410097; 0.634489; 1; 0.20504] atol = 1e-5

end


text:
using MISAlgorithms, EliminateGraphs
using Test

@testset "adjacent45" begin
    eg = K_eg(4,5)
    v1, v2 = EliminateGraphs.adjacent45(eg)
    @test Set([degree(eg,v2), degree(eg, v1)]) == Set([4,5])
end

@testset "mis2" begin
    @test EliminateGraphs.minx(x->x^2, [2,-1, 3]) == -1

    g1 = EliminateGraph([0 1 1 0 0;
                         1 0 1 1 0;
                         1 1 0 1 0;
                         0 1 1 0 0;
                         0 0 0 0 0])

    @test MISAlgorithms.mis2(g1) == 3
    g2 = disconnected_cliques_eg(6, 6)
    @test MISAlgorithms.mis2(g2) == 2
    g3 = K_eg(7, 2)
    @test MISAlgorithms.mis2(g3) == 7
    g4 = K_eg(7, 7)
    @test MISAlgorithms.mis2(g4) == 7
    g5 = ring_eg(4)
    @test MISAlgorithms.mis2(g5) == 2
    @test MISAlgorithms.mis2(c60_graph) == 24
    g6 = K_eg(4, 5)
    unsafe_disconnect!(g6, 1,5)
    unsafe_connect!(g6, 2,3)
    @test MISAlgorithms.mis2(g6) == 5
    g7 = K_eg(4, 5)
    @test MISAlgorithms.mis2(g7) == 5

    for i = 1:100
        eg = rand_eg(40, 0.2)
        @test MISAlgorithms.mis1(eg) == MISAlgorithms.mis2(eg)
    end
end


text:
function SeisFXDecon(in,param=Dict())

	dt = get(param,"dt",0.002)
	filter_length = get(param,"filter_length",10)
	mu = get(param,"mu",0.01)
	nt = size(in,1)
	nx = size(in,2)
	nf = nt
	dw = 2.*pi/nf/dt
	nw = int(nf/2) + 1
	fmax = get(param,"fmax",int(floor(0.5/dt)))
	if(fmax*dt*nf < nw)
		iw_max = int(floor(fmax*dt*nf))
	else
		iw_max = int(floor(0.5/dt))
	end
	D = fft(in,1)
	Df = zeros(Complex64,size(D))
	Db = zeros(Complex64,size(D))
	for iw=1:iw_max
		y = vec(D[iw,:])
		yf = fxdecon(y,"f",filter_length,mu);
		yb = fxdecon(y,"b",filter_length,mu);
		Df[iw,:] = yf
		Db[iw,:] = yb
	end
	# symmetries
	for iw=nw+1:nf
		Df[iw,:] = conj(Df[nf-iw+2,:])
		Db[iw,:] = conj(Db[nf-iw+2,:])
	end
	d = ( ifft(Df,1) + ifft(Db,1) )/ 2
	out = real(d[1:nt,:])

	return out
end

function fxdecon(d,direction,filter_length,mu)
	nx = length(d);
	if (direction=="b")
		y  = d[1:nx-filter_length];
		C  = d[2:nx-filter_length+1];
		R  = d[nx-filter_length+1:nx];
		M = hankel(C,R);
	else
		y  = d[filter_length+1:nx,1];
		C  = d[filter_length:nx-1];
		R = flipud(d[1:filter_length]);
		M = toeplitz(C,R);
	end
	B = M'*M;
	beta = B[1,1]*mu/100
	ab = (B + beta*eye(filter_length))\M'*y
	temp = M*ab
	if (direction=="b")
		d = [temp ; zeros(filter_length)]
	else
		d = [zeros(filter_length) ; temp]
	end

	return d
end

function hankel(c, r)

	nc = length(r)
	nr = length(c)
	if (nc > 1)
		c = [ c ; r[2:nc] ]
	end
	m = c[ones(nr) * [1:nc]' + [0:nr-1]*ones(nc)']
	m = reshape(m,nr,nc)

	return m
end

function toeplitz(c, r)

	nc = length(r)
	nr = length(c)
	m = zeros(Complex64,nr, nc);
	for i = 1:min(nc, nr)
		m[i:nr,i] = c[1:nr-i+1]
	end
	for i = 1:min(nr, nc-1)
		m[i,i+1:nc] = r[2:nc-i+1];
	end

	return m
end

max_stars_repo_name: fercarozzi/myseismicjulia

max_issues_repo_name: fercarozzi/myseismicjulia

max_forks_repo_name: fercarozzi/myseismicjulia


text:
@testset "Fixed fields" begin
  Qx,x = FlintQQ["x"]

  @testset "FixedField: absolut" begin
    f = x^8+236*x^6+12158*x^4+201788*x^2+779689
    K, a = NumberField(f, "a");
    s = hom(K, K, 43//339072*a^7 + 9265//339072*a^5 + 336481//339072*a^3 + 2429659//339072*a);
    t = hom(K, K, -115//1356288*a^7 - 26257//1356288*a^5 - 1190665//1356288*a^3 - 13355755//1356288*a);

    @test degree(fixed_field(K,t)[1])     == 4
    @test degree(fixed_field(K,s)[1])     == 4
    @test degree(fixed_field(K,[s,t])[1]) == 2
  end

  @testset "FixedField: realtiv" begin
    K,b     = NumberField(x-1,"b")
    Ky,y    = K["y"]
    f = y^8+236*y^6+12158*y^4+201788*y^2+779689
    L,a     = NumberField(f,"a")
    s = hom(L, L, 43//339072*a^7 + 9265//339072*a^5 + 336481//339072*a^3 + 2429659//339072*a);
    t = hom(L, L, -115//1356288*a^7 - 26257//1356288*a^5 - 1190665//1356288*a^3 - 13355755//1356288*a);

    @test degree(fixed_field(L,t)[1])     == 4
    @test degree(fixed_field(L,s)[1])     == 4
    @test degree(fixed_field(L,[s,t])[1]) == 2
  end
end

@testset "Subfields" begin
  Qx, x = FlintQQ["x"]
  f_1 = x^6 + 108;
  K_1,_ = NumberField(f_1);
  #(C_2)^3
  f_2 = x^8 - 12*x^6 + 23*x^4 - 12*x^2 + 1;
  K_2,_ = NumberField(f_2);
  #2xD_4
  f_3 = x^8 - 10*x^4 + 1;
  K_3,_ = NumberField(f_3);
  #[2^4]_4
  f_4 = x^8 + 4*x^6 + 10*x^4 + 12*x^2 + 7;
  K_4,_ = NumberField(f_4);
  #(3^2):2
  f_5 = x^9 - 18*x^8 + 117*x^7 - 348*x^6 + 396*x^5 + 288*x^4 + 3012*x^3 +576*x^2 + 576*x - 512;
  K_5,_ = NumberField(f_5);
  #(2^4):5
  f_6 = x^10 + 38*x^9 - 99*x^8 + 1334*x^7 - 4272*x^6 + 9244*x^5 - 8297*x^4 + 1222*x^3 + 1023*x^2 - 74*x + 1;
  K_6,_ = NumberField(f_6);
  #(2^4):D_5
  f_7 = x^10 -20*x^9 + 80*x^8 + 200*x^7 - 3770*x^6 + 872*x^5 + 29080*x^4 + 36280*x^3 - 456615*x^2 + 541260*x - 517448;
  K_7,_ = NumberField(f_7);
  #((5^2):4):2
  f_8 = x^10 - 10*x^8 + 20*x^7 + 235*x^6 + 606*x^5 +800*x^4 + 600*x^3 + 270*x^2 + 70*x + 16;
  K_8,_ = NumberField(f_8);
  #S_3xS_4
  f_9 = x^12 + 6*x^9 + 4*x^8 + 8*x^6 - 4*x^5 - 12*x^4 + 8*x^3 - 8*x + 8;
  K_9,_ = NumberField(f_9);
  #3xD_4
  f_10 = x^12 + 9*x^11 + 3*x^10 - 73*x^9 - 177*x^8 - 267*x^7 - 315*x^6 - 267*x^5 - 177*x^4 - 73*x^3 + 3*x^2 + 9*x + 1;
  K_10,_ = NumberField(f_10);
  #A_4
  f_11 = x^12 - 34734*x^11 + 401000259*x^10 - 1456627492885*x^9 - 2537142937228035*x^8 + 18762072755679375516*x^7 - 812368636358864062944*x^6 - 70132863629758257512231931*x^5 + 25834472514893102332821062085*x^4 + 76623280610352450247247939584745*x^3 - 45080885015422662132515763499758450*x^2 - 2070499552240812214288316981071818900*x - 550505759097778545485364826246753544;
  K_11,_ = NumberField(f_11);
  #((5^2):4)S_3
  f_12 = x^15 + 20*x^12 + 125*x^11 + 503*x^10 + 1650*x^9 + 3430*x^8 + 4690*x^7 + 4335*x^6 + 2904*x^5 + 1400*x^4 + 485*x^3 +100*x^2 + 15*x + 1;
  K_12,_ = NumberField(f_12);

  ### Ex.Tommy ###

  # Normal mit vielen Teilkoerpern #
  #(C_2^3)
  f_13 = x^8 - x^4 + 1;
  K_13,_ = NumberField(f_13);
  #(C_2^4)
  f_14 = x^16 - 7*x^12 + 48*x^8 - 7*x^4 + 1;
  K_14,_ = NumberField(f_14);
  #(C_2^5)              <--- cmpr runtime fact()
  f_15 = x^32+64*x^30+3072*x^28+78336*x^26+957696*x^24+17534976*x^22+72056832*x^20-10020814848*x^18+102596173824*x^16-271505948672*x^14-53842389499904*x^12-211197100032000*x^10+6793926729007104*x^8+34848900900716544*x^6+79559795278872576*x^4-159748047218147328*x^2+62856539035140096;
  K_15,_ = NumberField(f_15);
  # Normal mit wenig Teilkoerpern #
  #(Q_8)
  f_16 = x^8 + 12*x^6 + 36*x^4 + 36*x^2 + 9;
  K_16,_ = NumberField(f_16);
  #(Q_8)
  f_17 = x^8 + 3*x^7 + 22*x^6 + 60*x^5 + 201*x^4 + 450*x^3 + 1528*x^2 + 3069*x + 4561;
  K_17,_ = NumberField(f_17);
  #(Q_8)
  f_18 = x^8 + 60*x^6 + 900*x^4 + 4500*x^2 + 5625;
  K_18,_ = NumberField(f_18);
  #(Q_8)
  f_19 = x^8 + 84*x^6 + 1764*x^4 + 12348*x^2 + 21609;
  K_19,_ = NumberField(f_19);
  #(Q_8)
  f_20 = x^8 + 44*x^6 + 308*x^4 + 484*x^2 + 121;
  K_20, _ = NumberField(f_20);
  #(C_16)
  f_21 = x^16 - x^15 + x^14 - x^13 + x^12 - x^11 + x^10 - x^9 + x^8 - x^7 + x^6 - x^5 + x^4 - x^3 + x^2 - x + 1;
  K_21, _ = NumberField(f_21);
  #(C_32)        <---- cmpr runtime fact() 28sek
  f_22 = x^32+32*x^30+464*x^28+4032*x^26+23400*x^24+95680*x^22+283360*x^20+615296*x^18+980628*x^16+1136960*x^14+940576*x^12+537472*x^10+201552*x^8+45696*x^6+5440*x^4+256*x^2+2;
  K_22, _ = NumberField(f_22);
  # Nicht normal mit vielen Teilkoerpern #
  #(C_2^2 x S_3)
  f_23 = x^12 - 4*x^6 + 1;
  K_23, _ = NumberField(f_23);
  # Nicht normal mit wenig Teilkoerpern #
  #(S_8)
  f_24 = x^8 - x^7 + 2*x^6 - 3*x^5 + 3*x^4 - 3*x^3 + 3*x^2 - 2*x + 1;
  K_24, _ = NumberField(f_24);
  #(S_8)
  f_25 = x^8 - 2*x^7 + 3*x^6 - 3*x^5 + x^4 + 1;
  K_25, _ = NumberField(f_25);
  #(S_18)
  f_26 = x^18 - 2*x^17 + x^16 + 3*x^15 - 5*x^14 + x^13 + 5*x^12 - 8*x^11 + 2*x^10 + 7*x^9 - 11*x^8 + 2*x^7 + 9*x^6 - 8*x^5 + x^4 + 5*x^3 - 2*x^2 + 1;
  K_26, _ = NumberField(f_26);
  #(18T111)
  f_27 = x^18 - x^17 - x^16 + x^14 - 2*x^13 + x^12 + 5*x^11 + 6*x^10 - 8*x^9 - 6*x^8 - 3*x^7 + 5*x^6 + 6*x^5 + 2*x^4 - 5*x^3 - x^2 + 1;
  K_27, _ = NumberField(f_27);

  ### Ex.Hoeij ###

  #r=6
  #vH, K, N 2011: Total time: 30.690 seconds, Total memory usage: 93.44MB
  #Klueners 1998: Total time:  1.340 seconds, Total memory usage: 12.72MB
  f_28 = x^36-30*x^30+1197*x^24-25218*x^18+256770*x^12+15552*x^6+729;

  #r=7
  # vH, K, N 2011: Total time: 453.240 seconds, Total memory usage: 50.31MB
  # Klueners 1998: Total time:  75.230 seconds, Total memory usage: 23.66MB
  f_29 =  x^75+5*x^70-2535*x^65+17325*x^60+588650*x^55-456974*x^50+9905000*x^45-750608135*x^40-21267217530*x^35+
  3806950040*x^30-26676882954*x^25-3032520830*x^20-685818980*x^15-23381430*x^10-11997545*x^5-1;

  #explicit subfields
  K0 = subfields(K_4);
  L1,_ = NumberField(x^4 + 4*x^3 + 10*x^2 + 12*x + 7, "l1");
  L2,_ = NumberField(x^2 + 6*x + 7, "l2");

  @testset "Fields [2^4]_4" begin
    @test any(L -> isisomorphic(L[1],L1)[1], K0)
    @test any(L -> isisomorphic(L[1],L2)[1], K0)
  end

  ###testing .degree of fields##

  @testset "S_3" begin
    @test length(subfields(K_1)) == 6
    @test length(subfields(K_1, degree =6)) == 1
    @test length(subfields(K_1, degree =3)) == 3
    @test length(subfields(K_1, degree =2)) == 1
    @test length(subfields(K_1, degree =1)) == 1
  end

  @testset "(C_2)^3" begin
    @test length(subfields(K_2)) == 16
    @test length(subfields(K_2, degree =8)) == 1
    @test length(subfields(K_2, degree =4)) == 7
    @test length(subfields(K_2, degree =2)) == 7
    @test length(subfields(K_2, degree =1)) == 1
  end

  @testset "2xD_4" begin
    @test length(subfields(K_3)) == 8
    @test length(subfields(K_3, degree =8)) == 1
    @test length(subfields(K_3, degree =4)) == 3
    @test length(subfields(K_3, degree =2)) == 3
    @test length(subfields(K_3, degree =1)) == 1
  end

  @testset "[2^4]_4" begin
    @test length(subfields(K_4)) == 4
    @test length(subfields(K_4, degree =8)) == 1
    @test length(subfields(K_4, degree =4)) == 1
    @test length(subfields(K_4, degree =2)) == 1
    @test length(subfields(K_4, degree =1)) == 1
  end

  @testset "(3^2):2" begin
    @test length(subfields(K_5)) == 6
    @test length(subfields(K_5, degree =9)) == 1
    @test length(subfields(K_5, degree =3)) == 4
    @test length(subfields(K_5, degree =1)) == 1
  end

  @testset "(2^4):5" begin
    @test length(subfields(K_6)) == 3
    @test length(subfields(K_6, degree =10)) == 1
    @test length(subfields(K_6, degree =5)) == 1
    @test length(subfields(K_6, degree =1)) == 1
  end

  @testset "(2^4):D_5" begin
    @test length(subfields(K_7)) == 3
    @test length(subfields(K_7, degree =10)) == 1
    @test length(subfields(K_7, degree =5)) == 1
    @test length(subfields(K_7, degree =1)) == 1
  end

  @testset "((5^2):4):2" begin
    @test length(subfields(K_8)) == 3
    @test length(subfields(K_8, degree =10)) == 1
    @test length(subfields(K_8, degree =2)) == 1
    @test length(subfields(K_8, degree =1)) == 1
  end

  @testset "S_3xS_4" begin
    @test length(subfields(K_9)) == 4
    @test length(subfields(K_9, degree =12)) == 1
    @test length(subfields(K_9, degree =3)) == 1
    @test length(subfields(K_9, degree =4)) == 1
    @test length(subfields(K_9, degree =1)) == 1
  end

  @testset "3xD_4" begin
    @test length(subfields(K_10)) == 6
    @test length(subfields(K_10, degree =12)) == 1
    @test length(subfields(K_10, degree =6)) == 1
    @test length(subfields(K_10, degree =4)) == 1
    @test length(subfields(K_10, degree =3)) == 1
    @test length(subfields(K_10, degree =2)) == 1
    @test length(subfields(K_10, degree =1)) == 1
  end

  @testset "A_4" begin
    @test length(subfields(K_11)) == 10
    @test length(subfields(K_11, degree =12)) == 1
    @test length(subfields(K_11, degree =6)) == 3
    @test length(subfields(K_11, degree =4)) == 4
    @test length(subfields(K_11, degree =3)) == 1
    @test length(subfields(K_11, degree =1)) == 1
  end

  @testset "((5^2):4)S_3" begin
    @test length(subfields(K_12)) == 3
    @test length(subfields(K_12, degree =15)) == 1
    @test length(subfields(K_12, degree =3)) == 1
    @test length(subfields(K_12, degree =1)) == 1
  end
  #Normal many#
  @testset "(C_2^3)" begin
    @test length(subfields(K_13)) == 16
    @test length(subfields(K_13, degree =8)) == 1
    @test length(subfields(K_13, degree =4)) == 7
    @test length(subfields(K_13, degree =2)) == 7
    @test length(subfields(K_13, degree =1)) == 1
  end

  @testset "(C_2^4)" begin
    @test length(subfields(K_14)) == 67
    @test length(subfields(K_14, degree =16)) == 1
    @test length(subfields(K_14, degree =8)) == 15
    @test length(subfields(K_14, degree =4)) == 35
    @test length(subfields(K_14, degree =2)) == 15
    @test length(subfields(K_14, degree =1)) == 1
  end

  # @testset "(C_2^5)" begin      <--- cmpr runtime fact()
  #     @test length(subfields(f_15,32)) == 1
  #     @test length(subfields(f_15,16)) == 31
  #     @test length(subfields(f_15,8)) == 155
  #     @test length(subfields(f_15,4)) == 155
  #     @test length(subfields(f_15,2)) == 31
  #     @test length(subfields(f_15,1)) == 1
  # end
  #Normal few#
  @testset "(Q_8)" begin
    @test length(subfields(K_16)) == 6
    @test length(subfields(K_16, degree =8)) == 1
    @test length(subfields(K_16, degree =4)) == 1
    @test length(subfields(K_16, degree =2)) == 3
    @test length(subfields(K_16, degree =1)) == 1
  end

  @testset "(Q_8)" begin
    @test length(subfields(K_17)) == 6
    @test length(subfields(K_17, degree =8)) == 1
    @test length(subfields(K_17, degree =4)) == 1
    @test length(subfields(K_17, degree =2)) == 3
    @test length(subfields(K_17, degree =1)) == 1
  end

  @testset "(Q_8)" begin
    @test length(subfields(K_18)) == 6
    @test length(subfields(K_18, degree =8)) == 1
    @test length(subfields(K_18, degree =4)) == 1
    @test length(subfields(K_18, degree =2)) == 3
    @test length(subfields(K_18, degree =1)) == 1
  end

  @testset "(Q_8)" begin
    @test length(subfields(K_19)) == 6
    @test length(subfields(K_19, degree =8)) == 1
    @test length(subfields(K_19, degree =4)) == 1
    @test length(subfields(K_19, degree =2)) == 3
    @test length(subfields(K_19, degree =1)) == 1
  end

  @testset "(Q_8)" begin
    @test length(subfields(K_20)) == 6
    @test length(subfields(K_20, degree =8)) == 1
    @test length(subfields(K_20, degree =4)) == 1
    @test length(subfields(K_20, degree =2)) == 3
    @test length(subfields(K_20, degree =1)) == 1
  end

  @testset "(C_16)" begin
    @test length(subfields(K_21)) == 5
    @test length(subfields(K_21, degree =16)) == 1
    @test length(subfields(K_21, degree =8)) == 1
    @test length(subfields(K_21, degree =4)) == 1
    @test length(subfields(K_21, degree =2)) == 1
    @test length(subfields(K_21, degree =1)) == 1
  end

  @testset "(C_32)" begin
    @test length(subfields(K_22)) == 6
    # @test length(subfields(f_22,32)) == 1
    # @test length(subfields(f_22,16)) == 1
    # @test length(subfields(f_22,8)) == 1
    # @test length(subfields(f_22,4)) == 1
    # @test length(subfields(f_22,2)) == 1
    # @test length(subfields(f_22,1)) == 1
  end
  #Non normal many#
  @testset "(C_2^2 x S_3)" begin
    @test length(subfields(K_23)) == 10
    @test length(subfields(K_23, degree =12)) == 1
    @test length(subfields(K_23, degree =6)) == 3
    @test length(subfields(K_23, degree =4)) == 1
    @test length(subfields(K_23, degree =3)) == 1
    @test length(subfields(K_23, degree =2)) == 3
    @test length(subfields(K_23, degree =1)) == 1
  end
  #Non normal few#
  @testset "(S_8)" begin
    @test length(subfields(K_24)) == 2
    @test length(subfields(K_24, degree =8)) == 1
    @test length(subfields(K_24, degree =4)) == 0
    @test length(subfields(K_24, degree =2)) == 0
    @test length(subfields(K_24, degree =1)) == 1
  end

  @testset "(S_8)" begin
    @test length(subfields(K_25)) == 2
    @test length(subfields(K_25, degree =8)) == 1
    @test length(subfields(K_25, degree =4)) == 0
    @test length(subfields(K_25, degree =2)) == 0
    @test length(subfields(K_25, degree =1)) == 1
  end

  @testset "(S_18)" begin
    @test length(subfields(K_26)) == 2
    @test length(subfields(K_26, degree =18)) == 1
    @test length(subfields(K_26, degree =9)) == 0
    @test length(subfields(K_26, degree =6)) == 0
    @test length(subfields(K_26, degree =4)) == 0
    @test length(subfields(K_26, degree =3)) == 0
    @test length(subfields(K_26, degree =2)) == 0
    @test length(subfields(K_26, degree =1)) == 1
  end

  @testset "(18T111)" begin
    @test length(subfields(K_27)) == 6
    @test length(subfields(K_27, degree =18)) == 1
    @test length(subfields(K_27, degree =9)) == 1
    @test length(subfields(K_27, degree =6)) == 1
    @test length(subfields(K_27, degree =3)) == 2
    @test length(subfields(K_27, degree =2)) == 0
    @test length(subfields(K_27, degree =1)) == 1
  end

  @testset "Trivial fields" begin
    f = x - 1
    K, a = NumberField(f, "a")
    l = subfields(K)
    @test length(l) == 1
    @test degree(l[1][1]) == 1
    @test isisomorphic(l[1][1], K)[1]
  end

  @testset "relative fields" begin
    K, a = NumberField(x^2 +1, "a")
    Kx, x = K["x"]
    f = x^12 - 34734*x^11 + 401000259*x^10 - 1456627492885*x^9 - 2537142937228035*x^8 + 18762072755679375516*x^7 - 812368636358864062944*x^6 - 70132863629758257512231931*x^5 + 25834472514893102332821062085*x^4 + 76623280610352450247247939584745*x^3 - 45080885015422662132515763499758450*x^2 - 2070499552240812214288316981071818900*x - 550505759097778545485364826246753544;
    L, u = NumberField(f)

    @test length(subfields(L)) == 10
    @test length(subfields(L, degree = 6)) == 3

  end

end


text:
# RemoteCall
# ==========

struct RemoteCall
    # The task must return a future as the result value.
    call::Task
end

function Base.wait(rc::RemoteCall)
    return wait(rc.call)::Future
end

function Base.fetch(rc::RemoteCall)
    return fetch(wait(rc))
end

function Base.map(f, rc::RemoteCall)
    return f(rc)
end

function Base.serialize(::AbstractSerializer, rc::RemoteCall)
    throw(ArgumentError("cannot serialize a RemoteCall object"))
end


text:
using Test, Cassette, LinearAlgebra
using Cassette: @context, @pass, @overdub, overdub, recurse, hasmetadata,
                metadata, hasmetameta, metameta, untag, tag, enabletagging,
                untagtype, istagged, istaggedtype, Tagged, fallback, canrecurse,
                similarcontext, disablehooks

println("running unit tests")
@time @testset "unit tests" begin include("unittests.jl") end

println("running misc. tests (w/o tagging)")
@time @testset "misc. tests (w/o tagging)" begin include("misctests.jl") end

println("running misc. tests (w/ tagging)")
@time @testset "misc. tests (w/ tagging)" begin include("misctaggingtests.jl") end


text:
ComponentArray{A}(::UndefInitializer, ax::Axes) where {A<:StaticArrays.StaticArray,Axes<:Tuple} =
    ComponentArray(similar(A), ax...)

max_forks_repo_name: jonniedie/ComponentArrays


text:
module GadflySupport

using Stochy
import Stochy: ERP, support

export bar, density, histogram

gadflyloaded() = isdefined(Main, :Gadfly)
del(names...) = [Main.Gadfly.(n) for n in names]

function samples(erp::Discrete)
    @assert isa(erp.dict, Dict{Any,Int64})
    n = sum(values(erp.dict))
    ret = Array(Any, n)
    i = 1
    for (x,c) in erp.dict
        for _ in 1:c
            ret[i] = x
            i += 1
        end
    end
    ret
end

function bar(erp::ERP)
    @assert gadflyloaded()
    plot, Geom, Scale, Theme, mm, Guide = del(:plot, :Geom, :Scale, :Theme, :mm, :Guide)
    xs = sort(collect(support(erp)))
    ys = map(x->exp(score(erp, x)), xs)
    plot(x=xs, y=ys, Geom.bar(), Scale.x_discrete, Theme(bar_spacing=2mm), Guide.ylabel("p(x)"), Scale.y_continuous(minvalue=0, maxvalue=1))
end

function plotcontinuous(erp::Discrete, geometry, ylabel)
    @assert gadflyloaded()
    plot, Guide, Scale = del(:plot, :Guide, :Scale)
    plot(x=samples(erp), geometry, Scale.y_continuous(scalable=false), Guide.ylabel(ylabel))
end

function density(erp::Discrete)
    Geom, = del(:Geom)
    plotcontinuous(erp, Geom.density, "pdf(x)")
end

function histogram(erp::Discrete; bins=nothing)
    Geom, = del(:Geom)
    kwargs = Dict()
    bins != nothing && (kwargs[:bincount] = bins)
    plotcontinuous(erp, Geom.histogram(;kwargs...), "count")
end

end # module


text:
# BioStructures.jl
# ================
#
# A Julia package to read, write and manipulate macromolecular structures.
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/BioStructures.jl/blob/master/LICENSE.md

module BioStructures

using LinearAlgebra: dot, cross, norm, svd
using Statistics: mean

using CodecZlib
using Format
using RecipesBase
using LightGraphs
using MetaGraphs
using DataFrames
import BioCore
import BioCore.distance
import BioSymbols
import BioSequences.AminoAcidSequence
using BioAlignments
import MMTF: parsemmtf, writemmtf

include("model.jl")
include("pdb.jl")
include("mmcif.jl")
include("mmtf.jl")
include("spatial.jl")

end # BioStructures


text:
import CompEcon
using PyPlot
using PyCall

function main()
    f1(x) = exp(-2x)
    g1(x) = -2exp(-2x)

    function univariate_interp(f, g)
        # fit approximant
        n = 10; a = -1; b = 1
        basis = CompEcon.fundefn(:cheb, n, a, b)
        c = CompEcon.funfitf(basis, f1)

        # graph approximation error for function and derivative
        x = CompEcon.nodeunif(1001, a, b)[1]
        yact = f(x)
        dact = g(x)
        yfit = CompEcon.funeval(c, basis, x)[1]
        dfit = CompEcon.funeval(c, basis, x, 1)[1]

        # Nice plot of function approximation error
        fig, ax = plt.subplots()
        ax[:plot](x, yfit-yact)
        ax[:plot](x, 0*x, "k--", linewidth=2)
        ax[:set_xlabel]("x")
        ax[:set_ylabel]("Error")
        ax[:set_title](L"Error approximating function using Chebyshev")

        # Nice plot of derivative approximation error
        fig, ax = plt.subplots()
        ax[:plot](x, dfit-dact);
        ax[:plot](x, 0*x, "k--", linewidth=2)
        ax[:set_xlabel]("x")
        ax[:set_ylabel]("Error")
        ax[:set_title](L"Error approximating derivative using Chebyshev")
    end

    f2(x) = cos(x[:, 1]) ./ exp(x[:, 2])
    function bivaraite_interp(f)
        n = [7 7]
        a = [0 0]
        b = [1 1]
        basis = CompEcon.fundefn(:cheb, n, a, b)
        c = CompEcon.funfitf(basis, f2)

        fig = plt.figure()
        nplot = [101 101]
        x, xcoord = CompEcon.nodeunif(nplot, a, b)
        yfit = CompEcon.funeval(c, basis, x)[1]
        error = reshape(yfit - f2(x), nplot...)
        poly3d = surf(xcoord[1], xcoord[2], error)
        ax = poly3d[:get_axes]()
        ax[:set_title]("Chebyshev approximation error")

    end
end

main()


text:
include("../src/libjsp.jl")

import .Libjsp
using Plots

function fx(t)
    sin(2*pi*5500*t) + 0.5*sin(2*pi*2300*t + 3*pi/4) + 0.8*cos(2*pi*8150*t + 2*pi/7)
end

fsample = 24000
sampling_num = 256

ts = 1/fsample

x_n = zeros(Float64,sampling_num)

for i_sampling=0:sampling_num-1
    x_n[i_sampling+1] = fx(ts*i_sampling)
end

trangle_window = Libjsp.window("hanning", sampling_num)

x_window = x_n .* trangle_window

X_n = Libjsp.dft(sampling_num, x_window, fsample)

x0 = range(0, stop=sampling_num-1)

y0 = x_n.^2
y0 = 10 * broadcast(log10, y0)

y1 = x_window

y2 = broadcast(abs, X_n)
y2 = 10 * broadcast(log10, y2)

p0 = plot(x0, y0, marker = 2, linecolor= :blue)
p1 = plot(x0, y1, marker = 2, linecolor= :green)
p2 = plot(x0, y2, marker = 3, linecolor= :red)

plot(p0, p1, p2, layout=(3,1))


text:
AbstractTrees.children(invs::InvariantList) = invs.invariants
Base.show(io::IO, inv::AbstractInvariant) = AbstractTrees.print_tree(io, inv)
AbstractTrees.printnode(io::IO, inv::AbstractInvariant) = print(io, nameof(typeof(inv)), "(\"", title(inv), "\")")
AbstractTrees.children(::AbstractInvariant) = ()


text:
export TwoParabolasConstrained, TwoParabolasUnconstrained;

function two_parabolas_doc_str(symb)
    return """
    The canonical two parobola problem with two objectives taking two
    decision variables.

    ```math
        \\min_{𝐱 ∈ X }
        \\begin{bmatrix}
        (x₁ - 1)² + (x₂ - 1)² \\\\
        (x₁ + 1)² + (x₂ + 1)²
        \\end{bmatrix},
    ```

    where ``X = $(symb)``.
    """
end

@doc two_parabolas_doc_str("ℝ²")
struct TwoParabolasUnconstrained <: MOP end;

@doc two_parabolas_doc_str("[\\mathrm{lb}, \\mathrm{ub}]^2")
@with_kw struct TwoParabolasConstrained <: MOP
    lb :: Union{R,Vector{R}} where R<:Real = -4.0;
    ub :: Union{R,Vector{R}} where R<:Real = 4.0;
 end;

num_vars( ::Union{TwoParabolasUnconstrained, TwoParabolasConstrained} ) = 2;
num_objectives( ::Union{TwoParabolasUnconstrained, TwoParabolasConstrained} ) = 2;

function get_objectives( ::Union{TwoParabolasUnconstrained, TwoParabolasConstrained} )
    f1( x :: Vector{R} where R<:Real ) = sum( ( x .- 1.0 ).^2 );
    f2( x :: Vector{R} where R<:Real ) = sum( ( x .+ 1.0 ).^2 );
    return [f1, f2]
end

function constraints( mop::TwoParabolasConstrained )
    Box(
        isa( mop.lb, Real) ? fill( mop.lb, num_vars(mop) ) : mop.lb,
        isa( mop.ub, Real) ? fill( mop.ub, num_vars(mop) ) : mop.ub,
    )
end

function get_pareto_set( mop ::Union{TwoParabolasUnconstrained, TwoParabolasConstrained} )
    return SamplingFunction( mop, :ParetoSet, [:regular, :random] );
end

function get_pareto_front( mop ::Union{TwoParabolasUnconstrained, TwoParabolasConstrained} )
    return SamplingFunction( mop, :ParetoFront, [:regular, :random] );
end

function get_points( sample_func :: SamplingFunction, :: Union{TwoParabolasUnconstrained, TwoParabolasConstrained},
        ::Val{:ParetoSet}, num_points :: Int; method = :regular )

    if method == :regular
        x_range = num_points == 1 ? [0.0,] : range(-1.0, 1.0; length = num_points );
    elseif method == :random
        x_range = [ -1.0 + 2 * rand() for i = 1 : num_points ]
    end
    return [ [x;x] for x in x_range ]
end

function get_points( sample_func :: SamplingFunction, mop :: Union{TwoParabolasUnconstrained, TwoParabolasConstrained},
        ::Val{:ParetoFront}, num_points :: Int; method = :regular )

    pset = get_points( sample_func, mop, Val(:ParetoSet), num_points; method );
    F = get_vector_objective( mop );
    return F.( pset );
end


text:
# run the simulation

import GaussianMixtureTest
using Distributed
using DelimitedFiles

include(joinpath(dirname(pathof(GaussianMixtureTest)), "..",  "test/CompareAcuracy.jl"))
@everywhere Ctrue = 3
@everywhere C_max = max(5, (2*Ctrue - 1))
@everywhere B=200
@everywhere n_vec = [80, 100, 200, 300, 500, 800, 1000]
@everywhere nn = length(n_vec)

for ik in 1:nn
    @everywhere n = $(n_vec[ik])
    teststat = pmap((b) -> compareBIC(Ctrue, n, b), 1:B)
    writedlm("compare_$(Ctrue)_$(n).csv", teststat)
end


text:
ThreadsX.copy!(dest::AbstractVector, src::AbstractVector; kwargs...) =
    ThreadsX.copyto!(resize!(dest, length(src)), src; kwargs...)

function ThreadsX.copyto!(
    dest::AbstractArray,
    src::AbstractArray;
    basesize::Integer = default_copyto_basesize(dest, src),
)
    if length(dest) <= basesize
        copyto!(dest, src)
    elseif IndexStyle(dest) isa IndexLinear && IndexStyle(src) isa IndexLinear
        linear_copyto!(dest, src, basesize)
    else
        cartesian_copyto!(dest, src, basesize)
    end
    return dest
end

function linear_copyto!(dest, src, basesize)
    # TODO: support size-compatible but index-incompatible arrays
    @sync for p in _partition(eachindex(dest, src), basesize)
        @spawn if p isa AbstractUnitRange
            copyto!(dest, first(p), src, first(p), length(p))
        else
            copyto!(view(dest, p), view(src, p))
        end
    end
end

function cartesian_copyto!(dest, src, basesize)
    ThreadsX.foreach(eachindex(dest, src); basesize = basesize) do i
        @inbounds dest[i] = src[i]
    end
end

# TODO: Take into account more properties like: sizeof(T), boxed?, union?,
# need conversion?
function default_copyto_basesize(dest::AbstractArray{T}, ::AbstractArray) where {T}
    # 2^19 for 64 bit T
    basesize = 4194304 ÷ elsizeof(T)
    return max(cld(length(dest), Threads.nthreads()), basesize)
end


text:
using Flux, Statistics
using Flux: onehotbatch, onecold, logitcrossentropy
using Printf, BSON

include("model-commons.jl")
include("myplots.jl")

batch_size = 100
img_n = 48
img_m = 48

function make_batch_indices(i, n, batch_size=batch_size)
  istart = (i - 1) * batch_size + 1
  iend = if (istart+batch_size-1 > n) n else istart+batch_size-1 end
  convert(Int,istart):convert(Int,iend)
end

function train_minibatches(X, Y, batch_size=batch_size, img_n=img_n, img_m=img_m)
  n = size(X,2)
  bcount = ceil(n/batch_size)
  batches = []
  for i in 1:bcount
    inds = make_batch_indices(i, n)
    push!(batches, (reshape(X[:,inds], (img_n, img_m, 1, length(inds))), Y[:,inds]))
  end
  batches
end

function test_minibatch(X, Y, n=img_n, m=img_m)
  bsize = size(X,2)
  inds = make_batch_indices(1, bsize, bsize)
  (reshape(X[:,inds], (n, m, 1, bsize)), Y[:,inds])
end

dataset = create_data_splits()
X_train = reshape(dataset.X_train, (img_n, img_m, 1, size(dataset.X_train,2)))
X_test = reshape(dataset.X_test, (img_n, img_m, 1, size(dataset.X_test,2)))
train_set = train_minibatches(dataset.X_train, dataset.Y_train)
test_set = test_minibatch(dataset.X_test, dataset.Y_test)

@info("Constructing model...")
model = Chain(
    # First convolution, operating upon a 48x48 images
    Conv((8, 8), 1=>16, pad=(1,1), stride=2, relu),

    # Second convolution, operating upon a 22x22 image
    Conv((8, 8), 16=>32, pad=(1,1), stride=2, relu),

    # Reshape 3d tensor into a 2d one, at this point it should be (9, 9, 32, N)
    # which is where we get the 2592 in the `Dense` layer below:
    x -> reshape(x, :, size(x, 4)),
    Dense(2592, 10)
) |> gpu

# Make sure our model is nicely precompiled before starting our training loop
model(train_set[1][1])

loss(x, y) = logitcrossentropy(model(x), y)

function accuracy(x, y)
  logits = model(x)
  probs = [logsoftmax(logits[:,i]) for i in 1:size(logits,2)]
  mean(onecold.(probs) .== onecold(y))
end

opt = ADAM(0.001)

@info("Beginning training loop...")
best_acc = 0.0
last_improvement = 0
for epoch_idx in 1:100
    global best_acc, last_improvement
    # Train for a single epoch
    Flux.train!(loss, params(model), train_set, opt)

    # Calculate metrics:
    acc = accuracy(test_set...)
    err = loss(X_train, dataset.Y_train).data
    @info(@sprintf("[%d]: Test accuracy: %.4f", epoch_idx, acc))
    @info(@sprintf("[%d]: Current loss: %.4f", epoch_idx, err))

    # If our accuracy is good enough, quit out.
    if acc >= 0.999
        @info(" -> Early-exiting: We reached our target accuracy of 99.9%")
        break
    end

    # If this is the best accuracy we've seen so far, save the model out
    if acc >= best_acc
        @info(" -> New best accuracy! Saving model out to sound_cnn2.bson")
        BSON.@save "sound_cnn2.bson" model epoch_idx acc
        best_acc = acc
        last_improvement = epoch_idx
    end

    # If we haven't seen improvement in 5 epochs, drop our learning rate:
    if epoch_idx - last_improvement >= 3 && opt.eta > 1e-5
        opt.eta /= 10.0
        @warn(" -> Haven't improved in a while, dropping learning rate to $(opt.eta)!")

        # After dropping learning rate, give it a few epochs to improve
        last_improvement = epoch_idx
    end

    if epoch_idx - last_improvement >= 10
        @warn(" -> We're calling this converged.")
        break
    end
end

# Best Accuracy: 77.4% - File sound_cnn.bson
#   batch_size       = 100
#   spectrogram size = 48x48
#   eta              = 0.001
#   model:
#       Conv((8, 8), 1=>16, pad=(1,1), stride=2, relu)
#       Conv((8, 8), 16=>32, pad=(1,1), stride=2, relu),
#       Dense(2592, 10)
#       softmax


text:
using LowRankApprox, ApproxFun, SingularIntegralEquations, Test
    import LowRankApprox: _LowRankMatrix

@testset "Hierarchical Solve" begin
    M = 50
    L = Array{LowRankMatrix{Float64}}(undef,14)
    [L[i] = _LowRankMatrix(0.001rand(4M,2),0.001rand(4M,2)) for i=1:2]
    [L[i] = _LowRankMatrix(0.01rand(2M,2),0.01rand(2M,2)) for i=[3:4;9:10]]
    [L[i] = _LowRankMatrix(0.1rand(M,2),0.1rand(M,2)) for i=[5:8;11:14]]
    D = Array{Diagonal{Float64}}(undef,8)
    [D[i] = Diagonal(1 ./collect(1:M)) for i=1:8]

    H = HierarchicalMatrix(D,L)

    @test rank(H) == 400
    #@test blockrank(H) == fill(2,8,8)+diagm(fill(48,8))
    @test isfactored(H) == false

    B = Array{Vector{Float64}}(undef,8)
    [B[i] = rand(M) for i=1:8]
    b = rand(size(H,1))#HierarchicalVector(B)

    Matrix(H)\b
    @time x = Matrix(H)\b
    H\b # includes compile time
    H = HierarchicalMatrix(D,L)
    @time xw = H\b # includes precomputation
    @time xw = H\b # H.A's precomputed and factored :)

    @test isfactored(H) == true

    CH = cond(H)

    @test norm(x-xw) ≤ 10CH^2*eps()
    @test norm(H*x-b) ≤ 10CH*eps()
    @test norm(H*xw-b) ≤ 10CH*eps()

    H+H-(H-H)
    @test norm(Matrix(H-H)) < 10norm(Matrix(H))*eps()

    # Hierarchical operators case

    i = 2

    dom = cantor(Segment(),i)

    ⨍ = DefiniteLineIntegral(dom)
    f = Fun(x->logabs(x-5im),dom)
    sp = Space(dom)

    G = GreensFun((x,y)->0.5,CauchyWeight(sp⊗sp,0);method=:Cholesky)

    @time u1 = ⨍[G]\transpose(f)
    println("Adaptive QR  forward error norm is: ",norm(⨍[G]*u1-f))

    @test norm(⨍[G]*u1-f) < 10eps()

    hdom = clustertree(dom)
    hsp = Space(hdom)

    G1 = GreensFun((x,y)->0.5,CauchyWeight(hsp⊗hsp,0);method=:Cholesky)
    @test SingularIntegralEquations.blockrank(G1) == [1 8 9 9; 8 1 9 9; 9 9 1 8; 9 9 8 1]

    H = ⨍[G1]
    @test SingularIntegralEquations.blocksize(H) == (4,4)

    @time u2 = H\f
    @time u2 = H\f

    @test condest(H) ≤ 20.0

    import ApproxFun: ∞
    @test blockrank(H) == [∞ 8 9 9; 8 ∞ 9 9; 9 9 ∞ 8; 9 9 8 ∞]

    println("The hierarchical forward error norm is: ",norm(⨍[G]*u2-f))

    @test norm(⨍[G]*u2-f) < 10eps()

    # Test three domains

    dom = Segment(-2.0,-1.0+0im) ∪ Segment(1.0,2.0+0im) ∪ Segment(-0.5-2im,0.5-2im)
    hdom = clustertree(dom)

    @test convert(typeof(hdom),hdom)  == hdom

    dom = UnionDomain(SingularIntegralEquations.collectdata(hdom))
    ⨍ = DefiniteLineIntegral(dom)
    f = Fun(x->logabs(x-5im),dom)
    sp = Space(dom)

    G = GreensFun((x,y)->0.5,CauchyWeight(sp⊗sp,0);method=:Cholesky)

    @time u1 = ⨍[G]\transpose(f)

    println("Adaptive QR  forward error norm is: ",norm(⨍[G]*u1-f))

    @test norm(⨍[G]*u1-f) < 300eps()

    hsp = Space(hdom)

    G1 = GreensFun((x,y)->0.5,CauchyWeight(hsp⊗hsp,0);method=:Cholesky)
    H = ⨍[G1]

    @time u2 = H\f
    @time u2 = H\f

    println("The hierarchical forward error norm is: ",norm(⨍[G]*u2-f))

    @test norm(⨍[G]*u2-f) < 10000eps()
end


text:
Pkg.update()
Pkg.add("Primes")
Pkg.add("BenchmarkTools")
Pkg.add("Combinatorics")
Pkg.add("Lazy")
## For whatever reason this is necessary
Pkg.clone("https://github.com/scheinerman/Permutations.jl.git")

max_stars_repo_name: vmchale/julia-problem-solving

max_issues_repo_name: vmchale/julia-problem-solving

max_forks_repo_name: vmchale/julia-problem-solving


text:
using HepApi
using HTTP
using Test

# https://inspirehep.net/api/literature?sort=mostcited&page=3&q=a E.Witten.1

@testset "using all fields of Query" begin
	record_type = "literature"
	query = Dict(
							 "q" => "a E.Witten.1",
							 "sort" => "mostcited",
							 "size" => 10,
							 "page" => 1,
							 "fields" => "authors"
							 )

	# query = HepApi.Query(q, sort, size, page, fields)
	# @test true
	# @info "Made the Query"

	sch = HepApi.Search(record_type, query)
	@test true
	@info "Made the Search"

	uri = HepApi.URI(sch)
	@test true
	@info "URI: "*string(uri)
	@info "Made a URI"

	resp = HepApi.get(sch)
	@info "Successfully ran through the GET request"
	@info "The status code was $(resp.status)"
	@test true

	HepApi.resptojson(resp)
	@info "Converted a response to JSON"
	@test true
end;

@testset "testing search_inspirehep" begin

	record_type = "literature"
	query = Dict{String, Union{String, Integer}}(
							 "q" => "a E.Witten.1",
							 "sort" => "mostcited",
							 "size" => 10,
							 "page" => 1,
							 "fields" => "authors"
							 )

	sch = search_inspirehep(record_type, query)
	@test true
end;


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

using Markdown
using InteractiveUtils

# ╔═╡ 5ad5c202-20f8-11eb-23f1-4f38b687c285
using STMOZOO.Example

# ╔═╡ d2d007b8-20f8-11eb-0ddd-1181d4565a85
using Plots

# ╔═╡ 45189a82-20fa-11eb-0423-05ce1b84639d
using Zygote

# ╔═╡ 171fee18-20f6-11eb-37e5-2d04caea8c35
md"""
# Example: solving quadratic systems

By Michiel Stock

## Introduction

In this notebook, I have chosen to write code so minimize *quadratic functions* of the form:

$$f(\mathbf{x}) =  \frac{1}{2} \mathbf{x}^\intercal P\mathbf{x} + \mathbf{q} \cdot \mathbf{x} + r\,,$$

where $P \succ 0$.

Hence, solving

$$\min_\mathbf{x}\, f(\mathbf{x})\,.$$
"""

# ╔═╡ fb4aeb8c-20f7-11eb-0444-259de7b76883
md"""

## Derivation

The gradient of the quadratic function is

$$\nabla f(\mathbf{x})=P\mathbf{x} +\mathbf{q}\,.$$

Setting this to zero gives

$$\mathbf{x}^\star=-P^{-1}\mathbf{q}\,.$$

We also know that $\nabla^2f(\mathbf{x}) = P \succ 0$, so $\mathbf{x}^\star$ is a minimizer.
"""

# ╔═╡ 52e87238-20f8-11eb-2ea0-27ee1208d3c3
md"""
## Illustration

We will use a simple 2D example to illustrate this code.
"""

# ╔═╡ e9a9d69e-20f8-11eb-0d9d-330ee5e9cf25
P = [7 -1; -1 2]

# ╔═╡ fdd4e550-20f8-11eb-227b-25f36708484d
q = [-14.0, -2.0]

# ╔═╡ 025fd6e8-20f9-11eb-3e7d-3519f3c4b58f
r = 1

# ╔═╡ 096eff98-20f9-11eb-1e61-99d5714895ba
md"We have a function to generate the quadratic function."

# ╔═╡ 165509ca-20f9-11eb-107c-550cbba0f0e9
f_quadr = quadratic_function(P, q, r)

# ╔═╡ 1fffc82a-20f9-11eb-198c-c160d7dac87d
f_quadr([2, 1])

# ╔═╡ 26ab6ce2-20f9-11eb-1836-1756b290e5e3
md"No more need to remember the formulla for the minimizer! Just use `solve_quadratic_system`!"

# ╔═╡ 49832a8e-20f9-11eb-0841-19a40a12db18
x_star = solve_quadratic_system(P, q, r)

# ╔═╡ 55e0e274-20f9-11eb-36c0-753f228f7e9b
begin
	contourf(-20:0.1:20, -20:0.1:20, (x, y) -> f_quadr([x,y]), color=:speed)
	scatter!([x_star[1]], [x_star[2]], label="minimizer")
end

# ╔═╡ b1551758-20f9-11eb-3e8f-ff9a7127d7f8
md"""
## Approximating non-quadratic functions

We can approximate non-quadratic functions by a quadratic function: The second order Taylor approximation $\hat{f}$ of a function $f$ at $\mathbf{x}$ is
$$f(\mathbf{x}+\mathbf{v})\approx\hat{f}(\mathbf{x}+\mathbf{v}) = f(\mathbf{x}) + \nabla f(\mathbf{x})^\top \mathbf{v} + \frac{1}{2} \mathbf{v}^\top \nabla^2 f(\mathbf{x}) \mathbf{v}\,.$$

Let us use this idea for the Rosenbrock function.
"""

# ╔═╡ 41d8f1dc-20fa-11eb-3586-a989427c1fd6
f_nonquadr((x1, x2); a=1, b=5) = (a-x1)^2 + b * (x2 - x1^2)^2

# ╔═╡ 4ed4215e-20fa-11eb-11ee-f7741591163c
x = [0.0, 0.0]

# ╔═╡ 56af99ee-20fa-11eb-0240-69c675efb78c
fx = f_nonquadr(x)

# ╔═╡ 6c5473b4-20fa-11eb-327b-51ac560530eb
∇fx = f_nonquadr'(x)

# ╔═╡ 7518c2c0-20fa-11eb-32c0-a9db2a91cbc5
∇²fx = Zygote.hessian(f_nonquadr, x)

# ╔═╡ 34027942-20fb-11eb-261e-3b991ce4c9f8
v = solve_quadratic_system(∇²fx, ∇fx, fx)

# ╔═╡ 3bbeb85c-20fc-11eb-04d0-fb12d8ace50a
f̂(x′) = quadratic_function(∇²fx, ∇fx, fx)(x′ .- x)

# ╔═╡ 8623ac1a-20fa-11eb-2d45-49cce0fdac86
begin
	plot_nonquadr = contourf(-2:0.01:2, -2:0.01:2, (x, y) -> f_nonquadr([x,y]), color=:speed, title="non-quadratic function")
	scatter!(plot_nonquadr, [x[1]], [x[2]], label="x")
	scatter!(plot_nonquadr, [x[1]+v[1]], [x[2]+v[2]], label="x + v")

	plot_approx = contourf(-2:0.01:2, -2:0.01:2, (x, y) -> f̂([x,y]), color=:speed,
		title="quadratic approximation")
	scatter!(plot_approx, [x[1]], [x[2]], label="x")
	scatter!(plot_approx, [x[1]+v[1]], [x[2]+v[2]], label="x + v")

	plot(plot_nonquadr, plot_approx, layout=(2,1), size=(600, 800))

end

# ╔═╡ Cell order:
# ╟─171fee18-20f6-11eb-37e5-2d04caea8c35
# ╠═5ad5c202-20f8-11eb-23f1-4f38b687c285
# ╠═d2d007b8-20f8-11eb-0ddd-1181d4565a85
# ╟─fb4aeb8c-20f7-11eb-0444-259de7b76883
# ╟─52e87238-20f8-11eb-2ea0-27ee1208d3c3
# ╠═e9a9d69e-20f8-11eb-0d9d-330ee5e9cf25
# ╠═fdd4e550-20f8-11eb-227b-25f36708484d
# ╠═025fd6e8-20f9-11eb-3e7d-3519f3c4b58f
# ╟─096eff98-20f9-11eb-1e61-99d5714895ba
# ╠═165509ca-20f9-11eb-107c-550cbba0f0e9
# ╠═1fffc82a-20f9-11eb-198c-c160d7dac87d
# ╟─26ab6ce2-20f9-11eb-1836-1756b290e5e3
# ╠═49832a8e-20f9-11eb-0841-19a40a12db18
# ╠═55e0e274-20f9-11eb-36c0-753f228f7e9b
# ╟─b1551758-20f9-11eb-3e8f-ff9a7127d7f8
# ╠═41d8f1dc-20fa-11eb-3586-a989427c1fd6
# ╠═45189a82-20fa-11eb-0423-05ce1b84639d
# ╠═4ed4215e-20fa-11eb-11ee-f7741591163c
# ╠═56af99ee-20fa-11eb-0240-69c675efb78c
# ╠═6c5473b4-20fa-11eb-327b-51ac560530eb
# ╠═7518c2c0-20fa-11eb-32c0-a9db2a91cbc5
# ╠═34027942-20fb-11eb-261e-3b991ce4c9f8
# ╠═3bbeb85c-20fc-11eb-04d0-fb12d8ace50a
# ╟─8623ac1a-20fa-11eb-2d45-49cce0fdac86

max_stars_repo_name: manonmontyne/STMOZOO

max_issues_repo_name: manonmontyne/STMOZOO

max_forks_repo_name: manonmontyne/STMOZOO


text:
#!/usr/bin/env julia

#Start Test Script
using OpenStreetMapPlotter
import Test
using Pkg
# Run tests

println("Testing xml parsing...")
way_arr, bbox = parse_ways(open_file("west_philly.osm"))
test_nodes = way_arr[1].nodes[1:10]
Test.@test test_nodes[1].id == "27149008"
Test.@test test_nodes[2].id == "588195318"
Test.@test test_nodes[3].id == "27149007"
Test.@test test_nodes[4].id == "27149006"
Test.@test test_nodes[5].id == "27149005"
Test.@test test_nodes[6].id == "1288458928"
Test.@test test_nodes[7].id == "27149004"
Test.@test test_nodes[8].id == "27149003"
Test.@test test_nodes[9].id == "1288458923"
Test.@test test_nodes[10].id == "27149002"
Test.@test way_arr[1].tags == Dict{Any,Any}("name"=>"Schuylkill River","waterway"=>"river")
Test.@time save_json(way_arr, "test.geojson")
plot_ways(way_arr, bbox)


text:
#+Extracted_Item_Base
# Removes hash
function remove_hashtag(s::String)::String
    input=split(s,"\n",keepempty=false)
    output=Vector{String}()

    for s in input
        if s=="#"
            push!(output,"\n")
        else
            @assert(s[1:2]=="# ")
            push!(output,s[3:end]*"\n")
        end
    end
    return join(output)
end

# +Extracted_Item_Base,OrgString  L:org_string_comment
#
# *Note:* this function does a lot, it uses links.jl and
# *evaluate.jl. That is the reason why it is defined in its own file
#
function org_string_comment(di::Documented_Item,
                            di_array::Array{Documented_Item,1},
                            di_array_universe::Array{Documented_Item,1},
                            link_prefix::String,
                            boxingModule::String)::String

    s=raw_string_doc(di)

    # Process links
    s=doc_link_substitution(s,di_array,di_array_universe, link_prefix)
    # Execute code snippet
    s=with_hash_evaluate(s,boxingModule)
    # remove #
    s=remove_hashtag(s)

    # Quote string
    if !isempty(s)
        s="#+BEGIN_QUOTE\n$(s)#+END_QUOTE\n"
    end

    return s
end


text:
using YaoBase
export RepeatedBlock, repeat

"""
    RepeatedBlock{D,C,GT<:AbstractBlock} <: AbstractContainer{GT,D}

Repeat the same block on given locations.
"""
struct RepeatedBlock{D,C,GT<:AbstractBlock} <: AbstractContainer{GT,D}
    n::Int
    content::GT
    locs::NTuple{C,Int}
end

function RepeatedBlock(n::Int, block::AbstractBlock{D}, locs::NTuple{C,Int}) where {D,C}
    @assert_locs_safe n Tuple(i:i+nqudits(block)-1 for i in locs)
    nqudits(block) > 1 && throw(
        ArgumentError("RepeatedBlock does not support multi-qubit content for the moment."),
    )
    return RepeatedBlock{D,C,typeof(block)}(n, block, locs)
end

function RepeatedBlock(n::Int, block::AbstractBlock{D}, locs::UnitRange{Int}) where {D}
    (0 < locs.start) && (locs.stop <= n) ||
        throw(LocationConflictError("locations conflict."))
    nqudits(block) > 1 && throw(
        ArgumentError("RepeatedBlock does not support multi-qubit content for the moment."),
    )
    return RepeatedBlock{D,length(locs),typeof(block)}(n, block, Tuple(locs))
end

function RepeatedBlock(n::Int, block::GT) where {M,D,GT<:AbstractBlock{D}}
    return RepeatedBlock{D,n,GT}(n::Int, block, Tuple(1:nqudits(block):n-nqudits(block)+1))
end
YaoBase.nqudits(m::RepeatedBlock) = m.n

"""
    repeat(n, x::AbstractBlock[, locs]) -> RepeatedBlock{n}

Create a [`RepeatedBlock`](@ref) with total number of qubits `n` and the block
to repeat on given location or on all the locations.

# Example

This will create a repeat block which puts 4 X gates on each location.

```jldoctest; setup=:(using YaoBlocks)
julia> repeat(4, X)
nqudits: 4
repeat on (1, 2, 3, 4)
└─ X
```

You can also specify the location

```jldoctest; setup=:(using YaoBlocks)
julia> repeat(4, X, (1, 2))
nqudits: 4
repeat on (1, 2)
└─ X
```

But repeat won't copy the gate, thus, if it is a gate with parameter, e.g a `phase(0.1)`, the parameter
will change simultaneously.

```jldoctest; setup=:(using YaoBlocks)
julia> g = repeat(4, phase(0.1))
nqudits: 4
repeat on (1, 2, 3, 4)
└─ phase(0.1)

julia> g.content
phase(0.1)

julia> g.content.theta = 0.2
0.2

julia> g
nqudits: 4
repeat on (1, 2, 3, 4)
└─ phase(0.2)
```
"""
Base.repeat(n::Int, x::AbstractBlock, locs::Int...) = repeat(n, x, locs)
Base.repeat(n::Int, x::AbstractBlock, locs::NTuple{C,Int}) where {C} =
    RepeatedBlock(n, x, locs)
Base.repeat(n::Int, x::AbstractBlock, locs) = repeat(n, x, locs...)
Base.repeat(n::Int, x::AbstractBlock, locs::UnitRange) = RepeatedBlock(n, x, locs)
Base.repeat(n::Int, x::AbstractBlock) = RepeatedBlock(n, x)
Base.repeat(x::AbstractBlock) = @λ(n -> repeat(n, x))

"""
    repeat(x::AbstractBlock, locs)

Lazy curried version of [`repeat`](@ref).
"""
Base.repeat(x::AbstractBlock, locs) = @λ(n -> repeat(n, x, locs...))

occupied_locs(rb::RepeatedBlock) =
    (vcat([(i:i+nqudits(rb.content)-1) for i in rb.locs]...)...,)
chsubblocks(x::RepeatedBlock{D}, blk::AbstractBlock{D}) where {D} =
    RepeatedBlock(x.n, blk, x.locs)
PropertyTrait(x::RepeatedBlock) = PreserveAll()

mat(::Type{T}, rb::RepeatedBlock{D}) where {T,D} =
    hilbertkron(rb.n, fill(mat(T, rb.content), length(rb.locs)), [rb.locs...]; nlevel=D)
mat(::Type{T}, rb::RepeatedBlock{D,0,GT}) where {T,D,GT} = IMatrix{D^nqudits(rb),T}()

function _apply!(r::AbstractRegister, rp::RepeatedBlock)
    m = mat_matchreg(r, rp.content)
    for addr in rp.locs
        instruct!(r, m, Tuple(addr:addr+nqudits(rp.content)-1))
    end
    return r
end

# specialization
for G in [:X, :Y, :Z, :S, :T, :Sdag, :Tdag]
    GT = Expr(:(.), :ConstGate, QuoteNode(Symbol(G, :Gate)))
    @eval function _apply!(r::AbstractRegister, rp::RepeatedBlock{N,C,$GT}) where {N,C}
        instruct!(r, Val($(QuoteNode(G))), rp.locs)
        return r
    end
end

_apply!(reg::AbstractRegister, rp::RepeatedBlock{D,0}) where D = reg

cache_key(rb::RepeatedBlock) = (rb.locs, cache_key(rb.content))

Base.adjoint(blk::RepeatedBlock{D}) where {D} =
    RepeatedBlock(nqudits(blk), adjoint(blk.content), blk.locs)
Base.copy(x::RepeatedBlock) = RepeatedBlock(nqudits(x), x.content, x.locs)
Base.:(==)(A::RepeatedBlock, B::RepeatedBlock) = A.locs == B.locs && A.content == B.content

function YaoBase.iscommute(x::RepeatedBlock{D}, y::RepeatedBlock{D}) where {D}
    if nqudits(x) != nqudits(y)
        throw(QubitMismatchError("got nqudits = `$(nqudits(x))` and `$(nqudits(y))`"))
    end
    if x.locs == y.locs
        return iscommute(x.content, y.content)
    else
        iscommute_fallback(x, y)
    end
end


text:
# ---
# title: 68. Text Justification
# id: problem68
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: String
# link: <https://leetcode.com/problems/text-justification/description/>
# hidden: true
# ---
#
# Given an array of words and a width  _maxWidth_ , format the text such that
# each line has exactly _maxWidth_ characters and is fully (left and right)
# justified.
#
# You should pack your words in a greedy approach; that is, pack as many words
# as you can in each line. Pad extra spaces `' '` when necessary so that each
# line has exactly _maxWidth_ characters.
#
# Extra spaces between words should be distributed as evenly as possible. If the
# number of spaces on a line do not divide evenly between words, the empty slots
# on the left will be assigned more spaces than the slots on the right.
#
# For the last line of text, it should be left justified and no **extra** space
# is inserted between words.
#
# **Note:**
#
#   * A word is defined as a character sequence consisting of non-space characters only.
#   * Each word's length is guaranteed to be greater than 0 and not exceed _maxWidth_.
#   * The input array `words` contains at least one word.
#
#
#
# **Example 1:**
#
#
#
#     Input: words = ["This", "is", "an", "example", "of", "text", "justification."], maxWidth = 16
#     Output:
#     [
#        "This    is    an",
#        "example  of text",
#        "justification.  "
#     ]
#
# **Example 2:**
#
#
#
#     Input: words = ["What","must","be","acknowledgment","shall","be"], maxWidth = 16
#     Output:
#     [
#       "What   must   be",
#       "acknowledgment  ",
#       "shall be        "
#     ]
#     Explanation: Note that the last line is "shall be    " instead of "shall     be", because the last line must be left-justified instead of fully-justified.
#     Note that the second line is also left-justified becase it contains only one word.
#
# **Example 3:**
#
#
#
#     Input: words = ["Science","is","what","we","understand","well","enough","to","explain","to","a","computer.","Art","is","everything","else","we","do"], maxWidth = 20
#     Output:
#     [
#       "Science  is  what we",
#       "understand      well",
#       "enough to explain to",
#       "a  computer.  Art is",
#       "everything  else  we",
#       "do                  "
#     ]
#
#
#
# **Constraints:**
#
#   * `1 <= words.length <= 300`
#   * `1 <= words[i].length <= 20`
#   * `words[i]` consists of only English letters and symbols.
#   * `1 <= maxWidth <= 100`
#   * `words[i].length <= maxWidth`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
function to_real_field_univariate_function(obj::FieldUnivariateFunction, arg0::Field)
    return jcall(obj, "toRealFieldUnivariateFunction", RealFieldUnivariateFunction, (Field,), arg0)
end

function value(obj::FieldUnivariateFunction, arg0::RealFieldElement)
    return jcall(obj, "value", RealFieldElement, (RealFieldElement,), arg0)
end


text:
module SCEDCCorr

using SeisIO, SeisNoise, AWSS3, AWS, Dates, Glob, Statistics, JLD2
using CSV, DataFrames, SCEDC, Serialization
include("tools.jl")
include("autocorr.jl")
include("LHcorrelation.jl")
end


text:
using GLPK
const GLPK_PRIMAL_FACTORY = MOI.OptimizerWithAttributes(GLPK.Optimizer, MOI.Silent() => true)
const GLPK_DUAL_FACTORY = dual_optimizer(GLPK_PRIMAL_FACTORY)
const GLPK_PRIMAL_OPT = MOI.instantiate(GLPK_PRIMAL_FACTORY)
const GLPK_DUAL_OPT = MOI.instantiate(GLPK_DUAL_FACTORY)

push!(primal_linear_factory, GLPK_PRIMAL_FACTORY)
push!(dual_linear_factory, GLPK_DUAL_FACTORY)
push!(dual_linear_optimizer, GLPK_DUAL_OPT)
push!(primal_linear_optimizer, GLPK_PRIMAL_OPT)

@testset "GLPK Linear Problems" begin
    list_of_linear_problems = [
        lp1_test,
        lp2_test,
        lp3_test,
        lp4_test,
        lp5_test,
        lp6_test,
        lp7_test,
        # lp8_test, Int64 problem, does not work
        # lp9_test, Interval is not implemented
        lp10_test,
        # lp11_test, Feasibility not supported
        lp12_test,
        lp13_test
    ]
    test_strong_duality(list_of_linear_problems, GLPK_PRIMAL_FACTORY)
end


text:
" Define a step change of the mechanical power at t = 1.0"

function change_t_one(x,t,integrator)
  t in [1.0]
end

function Y_change!(integrator)
  integrator.p[2].Ybus = integrator.p[1]
end


text:
using ..SymEngine
import YaoBase: rot_mat

rot_mat(::Type{T}, ::Val{:Rx}, theta::Basic) where {T} =
    Basic[cos(theta / 2) -im*sin(theta / 2); -im*sin(theta / 2) cos(theta / 2)]
rot_mat(::Type{T}, ::Val{:Ry}, theta::Basic) where {T} =
    Basic[cos(theta / 2) -sin(theta / 2); sin(theta / 2) cos(theta / 2)]
rot_mat(::Type{T}, ::Val{:Rz}, theta::Basic) where {T} =
    Diagonal(Basic[exp(-im * theta / 2), exp(im * theta / 2)])
rot_mat(::Type{T}, ::Val{:CPHASE}, theta::Basic) where {T} = Diagonal(Basic[1, 1, 1, exp(im * theta)])
rot_mat(::Type{T}, ::Val{:PSWAP}, theta::Basic) where {T} = rot_mat(Basic, Const.SWAP, theta)

for G in [:Rx, :Ry, :Rz, :CPHASE]
    # forward single gates
    @eval function YaoBase.instruct!(
        state::AbstractVecOrMat{T},
        g::Val{$(QuoteNode(G))},
        locs::NTuple{N3,Int},
        control_locs::NTuple{N1,Int},
        control_bits::NTuple{N2,Int},
        theta::Basic,
    ) where {T,N1,N2,N3}
        m = rot_mat(T, g, theta)
        instruct!(state, m, locs, control_locs, control_bits)
    end

    # forward single gates
    @eval function YaoBase.instruct!(
        state::AbstractVecOrMat{T},
        g::Val{$(QuoteNode(G))},
        locs::NTuple{N1,Int},
        theta::Basic,
    ) where {T,N1}
        instruct!(state, g, locs, (), (), theta)
    end

    @eval function YaoBase.instruct!(
        state::AbstractVecOrMat{T},
        g::Val{$(QuoteNode(G))},
        locs::Tuple{Int},
        theta::Basic,
    ) where {T,N1}
        instruct!(state, g, locs, (), (), theta)
    end
end

@eval function YaoBase.instruct!(
    state::AbstractVecOrMat{T},
    g::Val{:PSWAP},
    locs::Tuple{Int,Int},
    theta::Basic,
) where {T,N1}
    instruct!(state, g, locs, (), (), theta)
end


text:
using Documenter, ANSIColoredPrinters

makedocs(
    clean = false,
    checkdocs = :exports,
    modules=[ANSIColoredPrinters],
    format=Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true",
                           assets = ["assets/default.css"]),
    sitename="ANSIColoredPrinters",
    pages=[
        "Introduction" => "index.md",
        "Output Formats" => "output-formats.md",
        "Supported Codes" => "supported-codes.md",
        "Reference" => "reference.md",
    ]
)

deploydocs(
    repo="github.com/JuliaDocs/ANSIColoredPrinters.jl.git",
    devbranch = "main",
    push_preview = true
)


text:
module VizEKP

using Plots
using JLD
using LinearAlgebra
using Statistics
using BoundingSphere

export ekp_sphere_evol
export plot_bounding_radius
export plot_ekp_params
export plot_error_evolution
export plot_covmat
export plot_output_profiles

"""
    function ekp_sphere_evol(ekp_u::Array{Array{Float64,2},1})

Returns the evolution with EKP iteration of the BoundingSphere
parameters (center, radius) for the EK particle ensemble.
"""
function ekp_sphere_evol(ekp_u::Array{Array{Float64,2},1})
    ekp_center = zeros((length(ekp_u), length(ekp_u[1][1,:])))
    ekp_radius = zeros(length(ekp_u))
    ekp_center_jump = zeros(length(ekp_u))
    #Iterate through EKP stages
    for i in 1:length(ekp_u)
        ekp_center[i,:], ekp_radius[i] = boundingsphere(
            mapslices(x->[x], ekp_u[i], dims=2)[:])
        if i > 1
            ekp_center_jump[i] = sqrt( (ekp_center[i,:]-ekp_center[i-1,:])'*(
                ekp_center[i,:]-ekp_center[i-1,:]) )
        end
    end
    return ekp_radius, ekp_center_jump
end

"""
    function plot_bounding_radius(
        ekp_u::Union{Array{Array{Float64,2},1}, Array{Array{Array{Float64,2},1},1}};
        normalized::Bool=false,
        newplot::Bool=true,
        ylims=nothing,)

Plots the evolution of the BoundingSphere radius for one or several EKPs.
"""
function plot_bounding_radius(
    ekp_u::Union{Array{Array{Float64,2},1}, Array{Array{Array{Float64,2},1},1}};
    normalized::Bool=true,
    newplot::Bool=true,
    ylims=nothing,)

    if ekp_u isa Array{Array{Array{Float64,2},1},1}
        # Wrapper for kwarg broadcasting
        wrapper(ekp_u) = plot_bounding_radius(
                            ekp_u, normalized=normalized,
                            newplot=false, ylims=ylims)
        plot()
        wrapper.(ekp_u)
    else

        ekp_radius, _ = ekp_sphere_evol(ekp_u)
        if normalized
            ekp_radius = ekp_radius./ekp_radius[1]
        end
        if newplot
            plot(ekp_radius, lw=2)
        else
            plot!(ekp_radius, ylims=ylims, lw=2)
        end
        xlabel!("EKP iteration")
        ylabel!("Min sphere radius")
        savefig("ekp_radius.png")
    end
    return
end

"""
    function plot_ekp_params(ekp_u::Array{Array{Float64,2},1};
                         exp_transform::Bool=false,
                         param_names::Union{Vector{String}, Nothing}=nothing)

Plots the evolution of EK parameters with EKP iteration. The ensemble is represented by
the mean, min and max values of each parameter at each EKP iteration.

If exp_transform, the parameter values are exponentiated before plotting.
"""
function plot_ekp_params(
    ekp_u::Array{Array{Float64,2},1};
    exp_transform::Bool=false,
    param_names::Union{Vector{String}, Nothing}=nothing)

    if exp_transform
        ekp_u_ = deepcopy( map(x->exp.(x), ekp_u) )
    else
        ekp_u_ = ekp_u
    end
    #One-dimensional measures
    ekp_mean = zeros( (length(ekp_u_), length(ekp_u_[1][1,:])) )
    ekp_min =  zero(ekp_mean)
    ekp_max =  zero(ekp_mean)
    #Iterate through EKP stages
    for i in 1:length(ekp_u_)
        ekp_mean[i,:] = mean(ekp_u_[i], dims=1)
        ekp_min[i,:] = minimum(ekp_u_[i], dims=1)
        ekp_max[i,:] = maximum(ekp_u_[i], dims=1)
    end

    for j in 1:length(ekp_u_[1][1,:])
        plot(ekp_mean[:,j], ribbon=(
            -ekp_min[:,j].+ekp_mean[:,j], ekp_max[:,j].-ekp_mean[:,j]),
            label=string("parameter ", j), lw=2)
        xlabel!("EKP iteration")
        if !isnothing(param_names)
            ylabel!(param_names[j])
        end
        xlabel!("EKP iteration")
        savefig(string("evol_param_",j,".png"))
    end
    return
end

function plot_ekp_params(
    ekp_u::Array{Array{Array{Float64,2},1},1};
    exp_transform::Bool=false,
    param_names::Union{Vector{String}, Nothing}=nothing)

    ekp_u_ = hcat(ekp_u...)
    if exp_transform
        ekp_u_ = deepcopy( map(x->exp.(x), ekp_u_) )
    end

    #One-dimensional measures
    ekp_mean = zeros( (length(ekp_u_[:,1]), length(ekp_u_[1,1][1,:])) )
    ekp_min =  zero(ekp_mean)
    ekp_max =  zero(ekp_mean)
    for j in 1:length(ekp_u_[1,1][1,:]) # For each parameter
        plot()
        for k in 1:length(ekp_u_[1,:]) # For each simulation
            #Iterate through EKP stages
            for i in 1:length(ekp_u_[:,1]) #For each iteration
                ekp_mean[i,:] = mean(ekp_u_[i,k], dims=1)
                ekp_min[i,:] = minimum(ekp_u_[i,k], dims=1)
                ekp_max[i,:] = maximum(ekp_u_[i,k], dims=1)
            end

            plot!(ekp_mean[:,j], ribbon=(
                -ekp_min[:,j].+ekp_mean[:,j], ekp_max[:,j].-ekp_mean[:,j]),
                label=string("parameter ", j), lw=2)
        end
        if !isnothing(param_names)
            ylabel!(param_names[j])
        end
        xlabel!("EKP iteration")
        savefig(string("evol_param_",j,".png"))
    end
    return
end

"""
    function plot_error_evolution(ekp_err::Union{Vector{Float64}, Array{Vector{Float64},1} };
                              ekp_std_scale::Union{Float64, Array{Float64,1}}=1.0 )

Plots the observation covariance-weighted EK error for each iteration. If Arrays are provided,
the function returns a single plot with the errors from all EK processes provided.

ekp_std_scale defines a proportionality constant for the errors from each EK process.

"""
function plot_error_evolution(
    ekp_err::Union{Vector{Float64}, Array{Vector{Float64},1} };
    ekp_std_scale::Union{Float64, Array{Float64,1}}=1.0,
    newplot::Bool=true,
    ylims=nothing, plt_scale=:identity)

    if ekp_err isa Array{Vector{Float64},1}
        # Wrapper for kwarg broadcasting
        wrapper(ekp_err, ekp_std_scale) = plot_error_evolution(
                            ekp_err, ekp_std_scale=ekp_std_scale,
                            newplot=false, plt_scale=plt_scale,
                            ylims=ylims)
        plot()
        wrapper.(ekp_err, ekp_std_scale)
    else
        ekp_err_ = deepcopy(ekp_err.*(ekp_std_scale*ekp_std_scale))

        if newplot
            plot(ekp_err_, lw=2)
        else
            plot!(ekp_err_, ylims=ylims, scale = plt_scale, lw=2)
        end
        xlabel!("EKP iteration")
        ylabel!("Error")
        savefig("ekp_error.png")
    end
    return
end

"""
    function plot_covmat(cov::Array{Float64, 2};
                     figname::Union{String, Nothing}=nothing)

Plots the given covariance matrix as a heat map.

"""
function plot_covmat(
    cov::Array{Float64, 2};
    figname::Union{String, Nothing}=nothing)

    # Scale by diagonal elements
    heatmap(cov, clim=(0, minimum(diag(cov)) ), yflip=true)
    if !isnothing(figname)
        figname_ = string(figname, ".png")
    else
        figname_ = "covmat.png"
    end
    savefig(figname_)
    return
end

"""
    function plot_output_profiles(ekp_g::Array{Array{Float64,2},1},
                      true_g::Array{Float64,1},
                      num_var::Int64,
                      num_heights::Int64)

Plots the mean output from the last EK ensemble, as well as the mean of the true output,
for a model whose outputs are vertical profiles of quantities of interest.=, with the
same number of vertical levels.
"""
function plot_output_profiles(
    ekp_g::Array{Array{Float64,2},1},
    true_g::Array{Float64,1},
    num_var::Array{Int64,1},
    num_heights::Array{Int64,1})

    for sim in 1:length(num_var)
        h_r = range(0, 1, length=Integer(num_heights[sim]))
        for k in 1:num_var[sim]
            inf_lim = Integer(
                num_heights[sim]*(k-1) + (num_heights[1:sim-1]'*num_var[1:sim-1]))+1
            sup_lim = Integer( inf_lim + num_heights[sim] - 1)
            y_model_mean = mean(ekp_g[end][:,inf_lim:sup_lim],dims=1)[1,:]

            plot(true_g[inf_lim:sup_lim] , h_r, label="Truth", lw=2)
            plot!(y_model_mean, h_r, label="Model", lw=2, ls=:dash)
            ylabel!("Normalized height")
            savefig(string("sim_",sim,"_output_variable",k,".png"))
        end
    end
    return
end

end #module