text:
using Test
using OpenQuantumSystems
using LinearAlgebra
using SparseArrays

# https://github.com/JuliaIO/Suppressor.jl/blob/master/src/Suppressor.jl
macro suppress(block)
    quote
        if ccall(:jl_generating_output, Cint, ()) == 0
            original_stdout = stdout
            out_rd, out_wr = redirect_stdout()
            out_reader = @async read(out_rd, String)

            original_stderr = stderr
            err_rd, err_wr = redirect_stderr()
            err_reader = @async read(err_rd, String)

            # approach adapted from https://github.com/JuliaLang/IJulia.jl/pull/667/files
            logstate = Base.CoreLogging._global_logstate
            logger = logstate.logger
            if logger.stream == original_stderr
                new_logstate =
                    Base.CoreLogging.LogState(typeof(logger)(err_wr, logger.min_level))
                Core.eval(Base.CoreLogging, Expr(:(=), :(_global_logstate), new_logstate))
            end
        end

        try
            $(esc(block))
        finally
            if ccall(:jl_generating_output, Cint, ()) == 0
                redirect_stdout(original_stdout)
                close(out_wr)

                redirect_stderr(original_stderr)
                close(err_wr)

                if logger.stream == stderr
                    Core.eval(Base.CoreLogging, Expr(:(=), :(_global_logstate), logstate))
                end
            end
        end
    end
end

@testset "aggregateCore" begin

    mode1 = Mode(0.2, 1.0)
    mode2 = Mode(0.3, 2.0)
    Energy = [0.0, 200.0]
    mol1 = Molecule([mode1, mode2], 3, [2.0, 200.0])
    mol2 = Molecule([mode1], 3, [3.0, 300.0])
    aggCore = AggregateCore([mol1, mol2])
    aggCore.coupling[2, 3] = 50
    aggCore.coupling[3, 2] = 50

    @test aggCore.molecules == [mol1, mol2]
    @test aggCore.coupling == [0.0 0.0 0.0; 0.0 0.0 50.0; 0.0 50.0 0.0]
    @test aggCore.molCount == 2

    @test getNvib(aggCore) == [[3, 3], [3]]
    @test getShifts(aggCore) == [[1.0, 2.0], [1.0]]
    @test getFrequencies(aggCore) == [[0.2, 0.3], [0.2]]

    @test getAggStateEnergy(aggCore, [1, 1], [[1, 1], [1]]) == 5.35
    @test getAggStateEnergy(aggCore, [1, 1], [[1, 2], [1]]) == 5.65
    @test getAggStateEnergy(aggCore, [1, 1], [[2, 1], [1]]) == 5.55
    @test getAggStateEnergy(aggCore, [1, 1], [[1, 1], [2]]) == 5.55
    @test getAggStateEnergy(aggCore, [2, 1], [[1, 1], [1]]) == 203.35

    @test OpenQuantumSystems.elIndOrder([1, 1, 1]) == 1
    @test OpenQuantumSystems.elIndOrder([2, 1, 1]) == 2
    @test OpenQuantumSystems.elIndOrder([1, 2, 1]) == 3

    aggCore_ = AggregateCore([mol1, mol2], [0.0 0.0 0.0; 0.0 0.0 50.0; 0.0 50.0 0.0])
    @test aggCore_ == aggCore
end

max_stars_repo_name: detrin/OpenQuantumSystems

max_issues_repo_name: detrin/OpenQuantumSystems

max_forks_repo_name: detrin/OpenQuantumSystems


text:
"""
return true / false based on hard thresholding of magnitudes above bandwidth
"""
function boxcar_accept(proposalDistance::G, bandwidth::G) where G <: Real
    return abs(proposalDistance) < bandwidth
end

function boxcar_accept(proposalDistance::Array{G, 1}, bandwidth::Array{G, 1}) where G <: Real
    for (dist, bw) in zip(proposalDistance, bandwidth)
        if dist >= bw
            return false
        end
    end
    return true
end

"""
return true / false based on gaussian kernel acceptance probability
"""
function gaussian_accept(proposalDistance::G, bandwidth::G) where G <: Real
    scaledDistane = proposalDistance / bandwidth
    return rand() < exp(-0.5 * scaledDistane ^ 2)
end

function gaussian_accept(proposalDistance::Array{G, 1}, bandwidth::Array{G, 1}) where G <: Real
    for (dist, bw) in zip(proposalDistance, bandwidth)
        scaledDistane = dist / bw
        if rand() >= exp(-0.5 * scaledDistane ^ 2)
            return false
        end
    end
    return true
end

ext: jl


text:
# https://github.com/invenia/PkgTemplates.jl/blob/4302ecb0a8f3304ca519796bae70c196794692c9/src/plugin.jl#L306-L316
"""
    gen_file(file, text::AbstractString)

Create a new file containing some given text.
Trailing whitespace is removed, and the file will end with a newline.
"""
function gen_file(file, text::AbstractString)
    mkpath(dirname(file))
    text = strip(join(map(rstrip, split(text, "\n")), "\n")) * "\n"
    return write(file, text)
end

function make_executable(path)
    try
        # Uses FilePathsBase so we can do "+x"
        chmod(Path(path), "+x")
    catch e
        @error "Could not run `chmod +x $path`. Try running this manually." exception = e
    end
    return nothing
end

function render_dir(template_dir, destination, variables; overwrite)
    for name in readdir(template_dir)
        text = render(read(joinpath(template_dir, name), String), variables)
        path = joinpath(destination, name)
        if !overwrite && isfile(path)
            error("$path already exists; set `overwrite=true` to overwrite existing files.")
        end
        @info("Writing $(path)")
        gen_file(path, text)
    end
    return nothing
end

function check_not_file(destination)
    isfile(destination) && throw(ArgumentError("""
            Destination $destination exists and is a file.
            Must be a directory (or nonexistent in which case a directory will be created).
        """))
    return nothing
end

"""
    setup_tensorboard(destination::AbstractString; app::AbstractString,
                      logdir::AbstractString,
                      ecr::AbstractString = default_ecr(),
                      service_account::AbstractString = default_service_account(),
                      namespace::AbstractString = get_current_namespace(),
                      local_port::Int=6006, overwrite=false)

Sets up a tensorboard launch scripts in destination directory `destination`, which
does not have to exist ahead of time.

* `app`: name of project/application/model that you are using tensorboard with
  (used for labelling the pods and keeping track of which is which)
* `logdir`: log directory that tensorboard should be pointing at. May be an `s3` URI
* `ecr`: name of the ECR repo to use to push docker images to. Can set this as a default
  with [`default_ecr`](@ref).
* `service_account`: name of the Kubernetes service account. Can set this as a default
  with [`default_service_account`](@ref).
* `namespace`: name of the Kubernetes namespace to use. Defaults to the current active
  one introspected by [`get_current_namespace`](@ref).
* `local_port`: local port to use when port forwarding.
* `overwrite`: whether or not to overwrite existing configuration scripts.

Creates:

* `\$destination/tensorboard.sh`: shell script that either connects to an existing tensorboard pod or builds and launches a new tensorboard pod.
* `\$destination/tensorboard.dockerfile`: dockerfile to build tensorboard with access to s3 log directories (used by `tensorboard.sh`)
* `\$destination/tensorboard.yaml`: pod spec used to launch pods ((used by `tensorboard.sh`)

These are pre-populated using the keyword arguments passed here.
"""
function setup_tensorboard(destination::AbstractString; app::AbstractString,
                           logdir::AbstractString, ecr::AbstractString=default_ecr(),
                           service_account::AbstractString=default_service_account(),
                           namespace::AbstractString=get_current_namespace(),
                           local_port::Int=6006, overwrite=false)
    check_not_file(destination)
    variables = Dict("app" => app, "logdir" => logdir, "ecr" => ecr,
                     "service_account" => service_account, "local_port" => local_port,
                     "namespace" => namespace)
    render_dir(joinpath(TEMPLATES, "tensorboard"), destination, variables; overwrite)
    make_executable(joinpath(destination, "tensorboard.sh"))
    @info tensorboard_instructions(destination)
    return nothing
end

function tensorboard_instructions(destination)
    return """
Tensorboard has been configured in $destination.

Make sure to add `BeaconK8sUtilities` to your global Julia environment so that it can be used from the `tensorboard.sh` script.

Then running `$destination/tensorboard.sh` in a shell should launch a tensorboard pod,
or give you the option to connect to an existing one.

Note: You can edit these files freely; running `setup_tensorboard` with `overwrite=true` will replace them with the latest defaults.
"""
end

"""
    setup_follow(destination::AbstractString;
                 namespace::AbstractString=get_current_namespace(), overwrite=false)

Sets up a `follow.sh` script in destination directory `destination`, which
does not have to exist ahead of time.

* `namespace`: name of the Kubernetes namespace to use. Defaults to the current active
  one introspected by [`get_current_namespace`](@ref).
* `overwrite`: whether or not to overwrite existing configuration scripts.

Creates:

* `\$destination/follow.sh`: shell script to wait for a pod to be ready (with progress information) and then follow its logs.
"""
function setup_follow(destination::AbstractString;
                      namespace::AbstractString=get_current_namespace(), overwrite=false)
    check_not_file(destination)
    variables = Dict("namespace" => namespace)
    render_dir(joinpath(TEMPLATES, "follow"), destination, variables; overwrite)
    make_executable(joinpath(destination, "follow.sh"))
    @info follow_instructions(destination)
    return nothing
end

function follow_instructions(destination)
    return """
`follow.sh` has been configured in $destination.

Make sure to add `BeaconK8sUtilities` to your global Julia environment so that it can be used from the `follow.sh` script.

Then running `$destination/follow.sh pod_name` should let you follow the progress of a pod.

Note: You can edit these files freely; running `setup_follow` with `overwrite=true` will replace them with the latest defaults.
"""
end


text:
export taper, taper!

# Create increasing part of a cosine taper
function mktaper!(W::Array{T,1}, L::Int64) where T<:Real
  LL = 2*L
  pp = 2.0*pi
  @inbounds for i = 1:L
    W[i] = 0.5*(1.0 - cos(pp*((i-1)/LL)))
  end
  return nothing
end

# Apply taper W to both ends of X
function taper_seg!(X::AbstractVector, W::Array{T,1}, L::Int64, μ::T; rev::Bool=false) where T<:Real
  if rev == true
    j = L
    @inbounds for i = 1:L
      X[j] = (X[j]-μ)*W[i] + μ
      j -= 1
    end
  else
    @inbounds for i = 1:L
      X[i] = (X[i]-μ)*W[i] + μ
    end
  end
  return nothing
end

# Taper a GphysChannel
@doc """
    taper!(C[; t_max::Real=10.0, α::Real=0.05, N_min::Int64=10])

Cosine taper all time-series data in C. Tapers each segment of each channel
that contains at least `N_min` total samples.

    taper!(S[; chans=CC, t_max::Real=10.0, α::Real=0.05, N_min::Int64=10])

Cosine taper each segment of time-series data in GphysChannel object C that
contains at least `N_min` total samples.

Does not modify irregularly-sampled data channels.

Keywords:
* `chans`: Only taper the specified channels.
* `N_min`: Data segments with N < N_min total samples are not tapered.
* `t_max`: Maximum taper edge in seconds.
* `α`: Taper edge area; as for a Tukey window, the first and last 100*α% of
samples in each window are tapered, up to `t_max` seconds of data.

See also: DSP.Windows.tukey
""" taper!
function taper!(C::GphysChannel; t_max::Real=10.0, α::Real=0.05, N_min::Int64=10)
  if !(C.fs > 0.0)
    return nothing
  end
  N_min = max(div(N_min,2), 0)
  # detrend!(C)

  # Reserve a window of length t_max*C.fs for the taper
  L = round(Int, t_max*C.fs)
  T = eltype(C.x)
  W = Array{T,1}(undef, L)

  # Determine window lengths of all segments
  window_lengths = diff(C.t[:,1])
  window_lengths[end] += 1
  n_seg = length(window_lengths)
  N_max = maximum(window_lengths)

  # Sanity check, with requirement α ≤ 0.5
  L_max = round(Int, t_max*C.fs)
  α = min(α, 0.5)
  L = min(L, round(Int, N_max*α))
  mktaper!(W, L)
  μ = T(mean(C.x))

  # Begin tapering by segment
  if n_seg == 1
    Nx = length(C.x)
    Xl = view(C.x, 1:L)
    Xr = view(C.x, Nx-L+1:Nx)
    taper_seg!(Xl, W, L, μ)
    taper_seg!(Xr, W, L, μ, rev=true)

  else
    # Get taper lengths and segment indices
    si = copy(C.t[1:n_seg, 1])
    ei = copy(C.t[2:n_seg, 1]).-1
    push!(ei, C.t[end,1])

    ii = sortperm(window_lengths, rev=true)
    si = si[ii]
    ei = ei[ii]
    Lw = window_lengths[ii]

    for i = 1:n_seg
      s = si[i]
      t = ei[i]
      L_tap = Lw[i]
      L_tap < N_min && break

      if L_tap < L
        L = min(L_max, round(Int, L_tap*α))
        resize!(W, L)
        mktaper!(W, L)
      end

      X = view(C.x, s:t)
      μ = T(sum(X)/(t-s+1))
      Xl = view(C.x, s:s+L-1)
      Xr = view(C.x, t-L+1:t)
      taper_seg!(Xl, W, L, μ)
      taper_seg!(Xr, W, L, μ, rev=true)
    end
  end
  note!(C,  string( "taper!, ",
                    "t_max = ", t_max, ", ",
                    "α = ", α,  ", ",
                    "N_min = ", N_min ) )
  return nothing
end

# This approach leads to heinous-looking code but uses virtually no memory.
# I could probably clean it up by creating one master taper
# and passing/editing views into the taper.

function taper!(S::GphysData;
  chans::Union{Integer, UnitRange, Array{Int64,1}}=Int64[],
  t_max::Real=10.0,
  α::Real=0.05,
  N_min::Int64=10)

  if !any(getfield(S, :fs) .> 0.0)
    return nothing
  end

  if chans == Int64[]
    chans = 1:S.n
  end

  α = min(α, 0.5)
  N_min = max(div(N_min,2), 0)
  T = unique([eltype(i) for i in S.x])
  nT = length(T)
  N_max = zeros(Int64, nT)

  # Arrays of views for each data type; we'll store these and lengths L_taps
  means = Array{Union{[Array{ty,1} for ty in T]...},1}(undef, nT)
  L_taps = Array{Array{Int64,1},1}(undef, nT)
  Xl = Array{Array{SubArray,1},1}(undef, nT)
  Xr = similar(Xl)
  for j = 1:nT
    Xl[j] = Array{SubArray{T[j],1,Array{T[j],1},Tuple{StepRange{Int64,Int64}},true},1}(undef,0)
    Xr[j] = similar(Xl[j])
    L_taps[j] = Array{Int64,1}(undef, 0)
    means[j] = Array{T[j],1}(undef, 0)
  end

  # Loop over channels, pushing views to the appropriate view array for eltype(S.x[i])
  for i = 1:S.n
    if i in chans
      S.fs[i] <= 0.0 && continue
      j = findfirst(T.==eltype(S.x[i]))

      # Compute length
      window_lengths = diff(S.t[i][:,1])
      window_lengths[end] += 1
      n_seg = length(window_lengths)

      # Get taper lengths and segment indices
      L_max = round(Int, t_max*S.fs[i])
      L_seg = min.(L_max, round.(Int, window_lengths*α))
      si = copy(S.t[i][1:n_seg, 1])
      ei = copy(S.t[i][2:n_seg, 1]).-1
      push!(ei, S.t[i][end,1])
      μ = T[j](mean(S.x[i]))

      for k = 1:length(L_seg)
        push!(Xl[j], view(S.x[i], si[k]:si[k]+L_seg[k]-1))
        push!(Xr[j], view(S.x[i], ei[k]-L_seg[k]+1:ei[k]))
        push!(means[j], μ)
        N_max[j] = max(N_max[j], L_seg[k])
      end
      append!(L_taps[j], L_seg)
    end
  end

  # Loop over data type
  for j = 1:nT
    L = N_max[j]
    W = Array{T[j],1}(undef, L)
    mktaper!(W, L)

    ii = sortperm(L_taps[j], rev=true)
    Nw = L_taps[j][ii]
    Xl[j] = Xl[j][ii]
    Xr[j] = Xr[j][ii]
    means[j] = means[j][ii]

    # Loop over left & right X-views
    for k = 1:length(Nw)
      Nw[k] < N_min && break
      # Note here: "break", not "continue", because we've sorted in reverse
      # order; once our taper regions become shorter than N_min, we're done.
      if Nw[k] < L
        L = Nw[k]
        resize!(W, L)
        mktaper!(W, L)
      end
      taper_seg!(Xl[j][k], W, L, means[j][k])
      taper_seg!(Xr[j][k], W, L, means[j][k], rev=true)
    end
  end

  # Annotate
  for i = 1:S.n
    if i in chans
      S.fs[i] == 0.0 && continue
      note!(S, i, string( "taper!, ",
                          "t_max = ", t_max, ", ",
                          "α = ", α,  ", ",
                          "N_min = ", N_min ) )
    end
  end
  return nothing
end
function taper(C::GphysChannel; t_max::Real=10.0, α::Real=0.05, N_min::Int64=10)
  U = deepcopy(C)
  taper!(U, t_max = t_max, α=α, N_min=N_min)
  return U
end
function taper(S::GphysData;
  chans::Union{Integer, UnitRange, Array{Int64,1}}=Int64[],
  t_max::Real=10.0,
  α::Real=0.05,
  N_min::Int64=10)

  U = deepcopy(S)
  taper!(U, chans=chans, t_max=t_max, α=α, N_min=N_min)
  return U
end


text:
function get_adj_list(tree::EvoTrees.Tree)
    adj = Vector{Vector{Int}}()
    for node in tree.nodes
        node.split ? push!(adj, [node.left, node.right]) : push!(adj, [])
    end
    return adj
end

function get_shapes(tree_layout)
    shapes = Vector(undef, length(tree_layout))
    for i in 1:length(tree_layout)
        x, y = tree_layout[i][1], tree_layout[i][2] # center point
        x_buff = 0.45
        y_buff = 0.45
        shapes[i] = [(x - x_buff, y + y_buff), (x + x_buff, y + y_buff), (x + x_buff, y - y_buff), (x - x_buff, y - y_buff)]
    end
    return shapes
end

function get_annotations(tree_layout, tree, var_names)
    # annotations = Vector{Tuple{Float64, Float64, String, Tuple}}(undef, length(tree_layout))
    annotations = []
    for i in 1:length(tree_layout)
        x, y = tree_layout[i][1], tree_layout[i][2] # center point
        if tree.nodes[i].split
            feat = isnothing(var_names) ? "feat: " *  string(tree.nodes[i].feat) : var_names[tree.nodes[i].feat]
            txt = "$feat\n" * string(round(tree.nodes[i].cond, sigdigits=3))
        else
            txt = "pred:\n" *  string(round(tree.nodes[i].pred[1], sigdigits=3))
        end
        # annotations[i] = (x, y, txt, (9, :white, "helvetica"))
        push!(annotations, (x, y, txt, 10))
    end
    return annotations
end

function get_curves(adj, tree_layout, shapes)
    curves = []
    num_curves = sum(length.(adj))
    for i in 1:length(adj)
        for j in 1:length(adj[i])
            # curves is a length 2 tuple: (vector Xs, vector Ys)
            push!(curves, ([tree_layout[i][1], tree_layout[adj[i][j]][1]], [shapes[i][3][2], shapes[adj[i][j]][1][2]]))
        end
    end
    return curves
end

@recipe function plot(tree::EvoTrees.Tree, var_names=nothing)

    adj = EvoTrees.get_adj_list(tree)
    tree_layout = length(adj) == 1 ? [[0.0,0.0]] : NetworkLayout.Buchheim.layout(adj)
    shapes = EvoTrees.get_shapes(tree_layout) # issue with Shape coming from Plots... to be converted o Shape in Receipe?
    annotations = EvoTrees.get_annotations(tree_layout, tree, var_names) # same with Plots.text
    curves = EvoTrees.get_curves(adj, tree_layout, shapes)

    size_base = floor(log2(length(adj)))
    size = (128 * 2^size_base, 96 * (1 + size_base))

    background_color --> :white
    linecolor --> :black
    legend --> nothing
    axis --> nothing
    framestyle --> :none
    size --> size
    annotations --> annotations

    for i in 1:length(shapes)
        @series begin
            fillcolor = length(adj[i]) == 0 ? "#84DCC6" : "#C8D3D5"
            fillcolor --> fillcolor
            seriestype --> :shape
            return shapes[i]
        end
    end

    for i in 1:length(curves)
        @series begin
            seriestype --> :curves
            return curves[i]
        end
    end
end

@recipe function plot(model::EvoTrees.GBTree, n=1, var_names=nothing)

    isnothing(var_names)

    tree = model.trees[n]
    adj = EvoTrees.get_adj_list(tree)
    tree_layout = length(adj) == 1 ? [[0.0,0.0]] : NetworkLayout.Buchheim.layout(adj)
    shapes = EvoTrees.get_shapes(tree_layout) # issue with Shape coming from Plots... to be converted o Shape in Receipe?
    annotations = EvoTrees.get_annotations(tree_layout, tree, var_names) # same with Plots.text
    curves = EvoTrees.get_curves(adj, tree_layout, shapes)

    size_base = floor(log2(length(adj)))
    size = (128 * 2^size_base, 96 * (1 + size_base))

    background_color --> :white
    linecolor --> :black
    legend --> nothing
    axis --> nothing
    framestyle --> :none
    size --> size
    annotations --> annotations

    for i in 1:length(shapes)
        @series begin
            fillcolor = length(adj[i]) == 0 ? "#84DCC6" : "#C8D3D5"
            fillcolor --> fillcolor
            seriestype --> :shape
            return shapes[i]
        end
    end

    for i in 1:length(curves)
        @series begin
            seriestype --> :curves
            return curves[i]
        end
    end
end


text:
import Pkg
using Pkg
Pkg.activate("GTheoryJulEnvNew")
Pkg.add("PkgTemplates")
Pkg.add(Pkg.PackageSpec(name="Cbc", version="0.7.1"))
Pkg.add(Pkg.PackageSpec(name="Clp", version="0.8.1"))
Pkg.add(Pkg.PackageSpec(name="GLPK", version="0.13.0"))
Pkg.add("Ipopt")
Pkg.add("Literate")
Pkg.add("DataStructures")
Pkg.add(Pkg.PackageSpec(name="JuMP", version="0.21.3"))
Pkg.add(Pkg.PackageSpec(name="Plots", version="1.6.3"))
Pkg.add(Pkg.PackageSpec(name="PrettyTables", version="0.9.1"))
Pkg.add(Pkg.PackageSpec(name="VegaLite", version="2.3.0"))
Pkg.add(Pkg.PackageSpec(name="CPLEX", version="0.6.1"))
Pkg.add(Pkg.PackageSpec(name="CSV", version="0.6.0"))
Pkg.add(Pkg.PackageSpec(name="DataArrays", version="0.7.0"))
Pkg.add(Pkg.PackageSpec(name="DataFrames", version="0.20.2"))
Pkg.add(Pkg.PackageSpec(name="Documenter", version="0.24.7"))
Pkg.add(Pkg.PackageSpec(name="DocumenterTools", version="0.1.9"))
Pkg.add(Pkg.PackageSpec(name="Gurobi", version="0.7.6"))
Pkg.add("MathOptInterface")
Pkg.add(Pkg.PackageSpec(name="MathProgBase", version="0.7.8"))
Pkg.add(Pkg.PackageSpec(name="StatsBase", version="0.33.0"))
Pkg.add(Pkg.PackageSpec(name="YAML", version="0.4.3"))
Pkg.add(Pkg.PackageSpec(name="LinearAlgebra"))
Pkg.status()

max_stars_repo_name: sambuddhac/TestOptJulia

max_issues_repo_name: sambuddhac/TestOptJulia

max_forks_repo_name: sambuddhac/TestOptJulia


text:
module TestSetup

using SpikeSorting

include("detect_test.jl")
include("feature_test.jl")

#include("bench_test.jl")
#include("compatibility_test.jl")

end


text:
function Unitful.ustrip(f::AbstractPICDataStructure)
    !hasunit(f) && return f
    data = unwrapdata(f)
    grid = getdomain(f)
    @debug "Stripping units for data of type $(typeof(data)) and grid $(typeof(grid))"
    udata = hasunit(data) ? ustrip_data(data) : data
    ugrid = ustrip(grid)

    newstruct(f, udata, ugrid)
end

# ustrip(A::Array) uses reinterpret and avoids creating a new array
ustrip_data(data) = ustrip(data)
ustrip_data(data::Base.ReinterpretArray) = reinterpret(recursive_unitless_bottom_eltype(data), data)
ustrip_data(data::SubArray) = reinterpret(recursive_unitless_bottom_eltype(data), data)
ustrip_data(data::StructArray) = StructArray(map(ustrip_data, components(data)))
ustrip_data(u, data) = ustrip.((u,), data)
ustrip_data(u, data::StructArray) = StructArray(map(c->ustrip.((u,),c), components(data)))
uconvert_data(u, data) = uconvert.((u,), data)
uconvert_data(u, data::StructArray) = StructArray(map(c->uconvert.((u,),c), components(data)))

for (f, ff) in zip((:uconvert, :ustrip), (:uconvert_data, :ustrip_data))
    @eval begin
        function (Unitful.$f)(u_data::Units, f::AbstractPICDataStructure)
            data = $ff(u_data, unwrapdata(f))
            grid = getdomain(f)

            newstruct(f, data, grid)
        end

        function (Unitful.$f)(u_data::Units, u_grid::Units, f::AbstractPICDataStructure)
            data = $ff(u_data, unwrapdata(f))
            grid = $f(u_grid, getdomain(f))

            newstruct(f, data, grid)
        end
    end
end

recursive_bottom_unit(f) = unit(recursive_bottom_eltype(f))
hasunit(f) = recursive_bottom_unit(f) ≠ NoUnits
hasunit(f::AbstractPICDataStructure) = hasunit(unwrapdata(f)) || hasunit(getdomain(f))
unitname(f) = hasunit(f) ? string(recursive_bottom_unit(f)) : ""


text:
include("../src/processing/solution_read.jl")

# Test de la méthode is_feasable(sol::Solution)

args = Args.parse_commandline()
Args.set(:level, 1)
Args.set(:timing_algo_solver, :earliest)
# Args.set(:level, 3)

# PARTIE 1
# Test sur l'instance instance_build_mini10 (instance alp01)
# Comme elle reste faisable quelque soit la permutation, on force
# l'infaisablilités en modifiant les caractéristiques d'un avion
# dans l'instance
#

lg1("faisabilité d'une solution PARTIE 1 : test sur alp01... ")
ln2("\nGénération d'une instance instance_build_mini10.")

inst = instance_build_mini10()
@test inst.nb_planes == 10

# Création de la solution vide
sol = Solution(inst)

ln2("La solution initiale de alp10 doit être faisable.")
@test is_feasable(sol) == true

if lg3()
    println(to_s_long(sol))
end

###############################
# on rend la sol infaisable en posant le UB du troisième avion égale
# au LB du premier avion (i.e ON MODIFIE L'INSTANCE) !
ln2("On modifie alp10 par p3.ub := p1.lb pour rendre infaisable.")

sol.planes[3].ub = sol.planes[1].lb
@test is_feasable(sol) == false
if lg3()
    println(to_s_long(sol))
end

ln1("fait.")

######################################################
# PARTIE 2
# Test sur l'instance alp_07 car très peu de solutions faisables.
# PRINCIPE
# - on teste l'instance trié par target croissant
#   => doit être faisable
# - on teste l'instance trié par target décroissant
#   => doit être infaisable
# - on teste NB_SHUFFLE mélanges (e.g 10_000 mélanges)
# - on doit trouver emte 97% et 98% de solutions infaisable !
#

# inst = instance_build_mini10()
# inst = Instance("data/07.alpx")

lg1("faisabilité d'une solution PARTIE 2 : test sur alp07... ")

DATA_DIR = dirname(@__FILE__) * "/data/"
NB_SHUFFLES = 10_000

inst = Instance("$DATA_DIR/alp_07.alp")

@test inst.nb_planes == 44

# Création de la solution vide
sol = Solution(inst)

initial_sort!(sol, presort=:target)
@test is_feasable(sol) == true

initial_sort!(sol, presort=:rtarget)
@test is_feasable(sol) == false

nb_infeasable = 0
if lg3()
    println("\nAvant $(NB_SHUFFLES) mélanges ")
end

for i in 1:NB_SHUFFLES
    shuffle!(sol)
    if !is_feasable(sol)
        global nb_infeasable += 1
        if lg4()
            println("shuffle! NOT FEASABLE : ",  sol )
        end
    end
end
if lg2()
    println("Après $(NB_SHUFFLES) mélanges => nb_infeasable=", nb_infeasable,
            " soit ", nb_infeasable/NB_SHUFFLES , " (test in [0.97, 0.98] ?")
end
@test  0.97 <= nb_infeasable/NB_SHUFFLES <= 0.98
ln1("fait.\n")


text:
using TensorNetworkCodes
using TensorNetworkCodes: _code_to_tensor
using Test
using ITensors

@testset "_code_to_tensor" begin
    # Test on the five-qubit code
    tensor = _code_to_tensor(five_qubit_code())
    # Check total number of code operators is correct
    @test sum(tensor) == 64.0
    # Check dimensions are correct
    @test size(tensor) == (4, 4, 4, 4, 4, 4)

    #Repeat for Steane code
    tensor = _code_to_tensor(steane_code())
    @test sum(tensor) == 256.0
    @test size(tensor) == (4, 4, 4, 4, 4, 4, 4, 4)
end

@testset "code_to_Itensor" begin
    logical_indices = [Index(4,"logical")]
    physical_indices = [Index(4,"physical") for _ in 1:5]
    tensor = code_to_Itensor(five_qubit_surface_code(),logical_indices,physical_indices)
    # Check dimensions
    @test dims(tensor) == (4, 4, 4, 4, 4, 4)
end

@testset "identity_coset" begin
    logical_indices = [Index(4,"logical")]
    physical_indices = [Index(4,"physical") for _ in 1:7]
    tensor = code_to_Itensor(steane_code(),logical_indices,physical_indices)
    id_coset = identity_coset(tensor)
    # Check number of stabilizers is correct
    @test sum(id_coset) == 64.0
end

@testset "all_cosets" begin
    logical_indices = [Index(4,"logical")]
    physical_indices = [Index(4,"physical") for _ in 1:5]
    tensor = code_to_Itensor(five_qubit_surface_code(),logical_indices,physical_indices)
    all_tensor = all_cosets(tensor)
    # Check number of stabilizers + logicals is correct
    @test sum(all_tensor) == 64.0
end

@testset "physical_tensor" begin
    index = Index(4,"logical")
    error_prob = 0.3
    pauli = 3
    tensor = physical_tensor(index,error_prob,pauli)
    arr = array(tensor)
    @test round.(arr,digits=1) == [0.1, 0.1, 0.1, 0.7]
end

@testset "create_virtual_tensor" begin
    code = five_qubit_surface_code()
    tn_code = TensorNetworkCode(code)
    tn_code = TensorNetworkCodes.contract(tn_code,deepcopy(tn_code),[[1,1]])
    tensor = create_virtual_tensor(tn_code,-1)
    # Check to see if tensor has all the right kind of indices
    @test hastags(tensor,"physical")
    @test hastags(tensor,"logical")
    @test hastags(tensor,"bond")
end


text:
"""
Fit Distributions to MC sim data for secondary infections per index case as a function of infection age

# Arguments
- `df::DataFrame` - DataFrame for csv containing MC simulations for daily individual infections.
- `horizon::Int=14` - Number of days in infection age before individual is considered naturally recovered and completely uninfectious.
- `sample_size::Int=50` - Sample size for `infections_path` csv where row entry is average infections for given sample size.
"""
function FitInfectionDistributions(df::DataFrame, horizon::Int=14, sample_size::Int=50)
    Λ = 1_000 # variance scaling parameter
    distributions = Vector{Gamma{Float64}}(undef, horizon)
    for day in 1:horizon
        try
            k, θ = Distributions.params(fit(Gamma, df[!,day]))
            distributions[day] = Gamma(k/(sample_size*Λ), Λ*θ*sample_size)
        catch e
            if e isa DomainError
                distributions[day] = Gamma(1e-100,1e-100)
            else
                throw(e)
            end
        end
    end
    return distributions
end

"""
Proportion of infectious population above some limit of detection from MC simulations

# Arguments
- `df::DataFrame` - DataFrame for csv containing MC simulations for daily individual infections.
- `day::Int` - Infection age (day)
- `LOD::Real` - Limit of detection (Log scale: ``10^x \\rightarrow x``)
"""
function prop_above_LOD(df::DataFrame, day::Int, LOD::Real)::Float64
    sum(df[!,day] .> LOD)/size(df,1)
end

"""
Return distribution resulting from sum of i.i.d RV's characterized by Gamma distribution

# Arguments
- `dist::Gamma` - Distribution characterizing random variable
- `N::Int` - Number of i.i.d RV's summed
"""
function RVsum(dist::Gamma, N::Int)
    k, θ = Distributions.params(dist)
    if k*N > 0
        return Gamma(k*N, θ)
    else
        return Gamma(1e-100, 1e-100)
    end
end


text:
using Documenter, BridgeSPDE

makedocs(;
    modules=[BridgeSPDE],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/mschauer/BridgeSPDE.jl/blob/{commit}{path}#L{line}",
    sitename="BridgeSPDE.jl",
    authors="mschauer <moritzschauer@web.de>",
    assets=String[],
)

deploydocs(;
    repo="github.com/mschauer/BridgeSPDE.jl",
)


text:
cx = (TruncateStd(1), TruncateStd(1))
cy = (TruncateStd(1), TruncateStd(1))
cs = ConstrainedIndexValueResampling(cx, cy)

# If n is not provided, set it to 1
ctest = ConstrainedTest(CrossMappingTest(), cs)
@test ctest isa ConstrainedTest
@test ctest.n == 1

ctest = ConstrainedTest(CrossMappingTest(), cs, 10)
@test ctest isa ConstrainedTest
@test ctest.n == 10


text:
using CSV

# Utility file I/O functions. Not fit for general purpose use.
function load_dipha(filename)
    open(filename, "r") do f
        magic = read(f, Int)
        @assert magic == 8067171840
        type = read(f, Int)
        @assert type == 1
        len = read(f, Int)
        dim = read(f, Int)
        if dim == 2
            m = read(f, Int)
            n = read(f, Int)
            result = Array{Float64}(undef, (m, n))
            read!(f, result)

            return result
        elseif dim == 3
            m = read(f, Int)
            n = read(f, Int)
            o = read(f, Int)
            result = Array{Float64}(undef, (m, n, o))
            read!(f, result)

            return result
        else
            error("not implemented for dim=$dim")
        end
    end
end

function save_dipha(filename, data)
    open(filename, "w") do f
        write(f, 8067171840)
        write(f, 1)
        write(f, length(data))
        write(f, length(size(data)))
        for s in size(data)
            write(f, s)
        end
        for i in data
            write(f, Float64(i))
        end
    end
end

function load_points(filename)
    # ripser uses Float32
    file = CSV.File(filename, header=0, type=Float32, delim=' ')
    nrow = length(file)
    ncol = length(file[1])
    result = Vector{NTuple{ncol, Float32}}(undef, nrow)
    for i in 1:nrow
        result[i] = Tuple(file[i, :]...)
    end
    return result
end

function load_dist(filename)
    file = CSV.File(filename, header=0, type=Float32, delim=' ')
    result = zeros(Float32, (length(file), length(file)))
    for (i, row) in enumerate(file)
        result[i, :] .= collect(row)
    end
    return result
end

function save_sparse(filename, data)
    I, J, V = findnz(sparse(LowerTriangular(data)))
    return CSV.write(joinpath(filename), (I=I.-1, J=J.-1, V=V), writeheader=false, delim=' ')
end

function load_sparse(filename)
    table = CSV.File(filename, header=0, type=Float32, delim=' ')
    I = table[:, 1] .+ 1
    J = table[:, 2] .+ 1
    V = table[:, 3] .+ 1
    sparse([I; J], [J; I], [V; V])
end

function load_data(filename)
    ext = splitext(filename)[2]
    if ext == ".pts" || ext == ".alpha"
        return load_points(filename)
    elseif ext == ".dist"
        return load_dist(filename)
    elseif ext == ".spdist"
        return load_sparse(filename)
    elseif ext == ".dipha"
        return load_dipha(filename)
    else
        error("unsupported filetype \"$(ext)\"")
    end
end

using BenchmarkTools: prettytime

function print_table_eirene(res)
    for k in keys(res)
        r = res[k]
        println(
            k, " & & & & ", prettytime(minimum(r["cohomology"]).time),
            " & ", prettytime(minimum(r["involuted"]).time),
            " & ", prettytime(minimum(r["eirene"]).time),
            " & ", minimum(r["involuted"]).time / minimum(r["eirene"]).time
        )
    end
end

function print_table_matrix(res)
    for k in keys(res)
        r = res[k]
        if haskey(r, "homology")
            println(
                k, " & & & & ", prettytime(minimum(r["cohomology"]).time),
                " & ", prettytime(minimum(r["homology"]).time),
                " & ", prettytime(minimum(r["involuted"]).time),
            )
        else
            println(
                k, " & & & & ", prettytime(minimum(r["cohomology"]).time),
                " & & ", prettytime(minimum(r["involuted"]).time),
            )
        end
    end
end


text:
## stats.jl  Various ways of computing Baum Welch statistics for a GMM
## (c) 2013--2014 David A. van Leeuwen

mem=0.1                          # Working memory for stats extraction, in Gig

## you can set the available working memory for stats calculation---this should be
## quite a bit less than the available memory, since there is (generally unknown) overhead
function setmem(gig::Float64)
    global mem=gig
end

## stats(gmm, x, order) computes zero, first, ... upto order (≤2) statistics of
## a feature file aligned to the gmm.  The statistics are ordered (ng
## * d), as by the general rule for dimension order in types.jl.
## Note: these are _uncentered_ statistics.

## you can dispatch this routine by specifying 3 parameters,
## i.e., an unnamed explicit parameter order

## For reasons of accumulation, this function returns a tuple
## (nx, loglh, N, F [S]).

## This function is admittedly hairy: in Octave this is much more
## efficient than a straightforward calculation.  I don't know if this
## holds for Julia.  We'd have to re-implement using loops and less
## memory.  I've done this now in several ways, it seems that the
## matrix implementation is always much faster.

## The shifting in dimensions (for Gaussian index k) is a nightmare.

## The memory footprint is sizeof(T) * ((2d + 2) ng + (d + ng + 1) nx, and
## results take an additional (2d +1) ng
## This is not very efficient, since this is designed for speed, and
## we don't want to do too much in-memory yet.
## Currently, I don't use a logsumexp implementation because of speed considerations,
## this might turn out numerically less stable for Float32

## diagonal covariance
function stats(gmm::GMM{GT,DCT}, x::Matrix{T}, order::Int) where DCT <: DiagCov{GT} where {GT, T}
    RT = promote_type(GT,T)
    (nₓ, d) = size(x)
    ng = gmm.n
    gmm.d == d || error("dimension mismatch for data")
    1 ≤ order ≤ 2 || error("order out of range")
    prec::Matrix{RT} = 1 ./ gmm.Σ             # ng × d
    mp::Matrix{RT} = gmm.μ .* prec          # mean*precision, ng × d
    ## note that we add exp(-sm2p/2) later to pxx for numerical stability
    a::Matrix{RT} = gmm.w ./ ((2π)^(d/2) * sqrt.(prod(gmm.Σ, dims=2))) # ng × 1
    sm2p::Matrix{RT} = dot(mp, gmm.μ, 2)    # sum over d mean^2*precision, ng × 1
    xx = x .* x                            # nₓ × d
##  γ = broadcast(*, a', exp(x * mp' .- 0.5xx * prec')) # nₓ × ng, Likelihood per frame per Gaussian
    γ = x * mp'                            # nₓ × ng, nₓ * d * ng multiplications
    LinearAlgebra.BLAS.gemm!('N', 'T', -one(RT)/2, xx, prec, one(RT), γ)
    for j = 1:ng
        la = log(a[j]) - 0.5sm2p[j]
        for i = 1:nₓ
            @inbounds γ[i,j] += la
        end
    end
    for i = 1:length(γ) @inbounds γ[i] = exp(γ[i]) end
    lpf=sum(γ,dims=2)                           # nₓ × 1, Likelihood per frame
    broadcast!(/, γ, γ, lpf .+ (lpf .== 0)) # nₓ × ng, posterior per frame per gaussian
    ## zeroth order
    N = vec(sum(γ, dims=1))          # ng, vec()
    ## first order
    F =  γ' * x                           # ng × d, Julia has efficient a' * b
    llh = sum(log.(lpf))                   # total log likeliood
    if order==1
        return (nₓ, llh, N, F)
    else
        ## second order
        S = γ' * xx                       # ng × d
        return (nₓ, llh, N, F, S)
    end
end

## Full covariance
## this is a `slow' implementation, based on posterior()
function stats(gmm::GMM{GT,FCT}, x::Array{T,2}, order::Int) where FCT <: FullCov{GT} where {GT, T}
    RT = promote_type(GT,T)
    (nₓ, d) = size(x)
    ng = gmm.n
    gmm.d == d || error("dimension mismatch for data")
    1 ≤ order ≤ 2 || error("order out of range")
    γ, ll = gmmposterior(gmm, x) # nₓ × ng, both
    llh = sum(logsumexp(ll .+ log.(gmm.w)', 2))
    ## zeroth order
    N = vec(sum(γ, dims=1))
    ## first order
    F = γ' * x
    if order == 1
        return nₓ, llh, N, F
    end
    ## S_k = Σ_i γ _ik x_i' * x
    S = Matrix{RT}[]
    γx = Array{RT}(undef, nₓ, d)
    @inbounds for k=1:ng
        #broadcast!(*, γx, γ[:,k], x) # nₓ × d mults
        for j = 1:d for i=1:nₓ
            γx[i,j] = γ[i,k]*x[i,j]
        end end
        push!(S, x' * γx)            # nₓ * d^2 mults
    end
    return nₓ, llh, N, F, S
end

## ## reduction function for the plain results of stats(::GMM)
## function accumulate(r::Vector{Tuple})
##     res = {r[1]...}           # first stats tuple, as array
##     for i=2:length(r)
##         for j = 1:length(r[i])
##             res[j] += r[i][j]
##         end
##     end
##     tuple(res...)
## end

## split computation up in parts, either because of memory limitations
## or because of parallelization
## You dispatch this by only using 2 parameters
function stats(gmm::GMM, x::Matrix{T}; order::Int=2, parallel=false) where T <: AbstractFloat
    parallel &= nworkers() > 1
    ng = gmm.n
    (nₓ, d) = size(x)
    if kind(gmm) == :diag
        bytes = sizeof(T) * ((2d +2)ng + (d + ng + 1)nₓ)
    elseif kind(gmm) == :full
        bytes = sizeof(T) * ((d + d^2 + 5nₓ + nₓ*d)ng + (2d + 2)nₓ)
    end
    blocks = ceil(Integer, bytes / (mem * (1<<30)))
    if parallel
        blocks= min(nₓ, max(blocks, nworkers()))
    end
    l = nₓ / blocks     # chop array into smaller pieces xx
    xx = Matrix{T}[x[round(Int, i*l+1):round(Int, (i+1)l),:] for i=0:(blocks-1)]
    if parallel
        r = pmap(x->stats(gmm, x, order), xx)
        reduce(+, r)                # get +() from BigData.jl
    else
        r = stats(gmm, popfirst!(xx), order)
        for x in xx
            r += stats(gmm, x, order)
        end
        r
    end
end
## the reduce above needs the following
Base.zero(x::Array{Matrix{T}}) where T = [zero(z) for z in x]

## This function calls stats() for the elements in d::Data, irrespective of the size, or type
function stats(gmm::GMM, d::Data; order::Int=2, parallel=false)
    if parallel
        return dmapreduce(x->stats(gmm, x, order=order, parallel=false), +, d)
    else
        r = stats(gmm, d[1], order=order, parallel=false)
        for i=2:length(d)
            r += stats(gmm, d[i], order=order, parallel=false)
        end
        return r
    end
end

## Same, but UBM centered+scaled stats
## f and s are ng * d
function csstats(gmm::GMM, x::DataOrMatrix{T}, order::Int=2) where T<:AbstractFloat
    kind(gmm) == :diag || error("Can only do centered and scaled stats for diag covariance")
    if order==1
        nₓ, llh, N, F = stats(gmm, x, order)
    else
        nₓ, llh, N, F, S = stats(gmm, x, order)
    end
    Nμ = N .* gmm.μ
    f = (F - Nμ) ./ gmm.Σ
    if order==1
        return(N, f)
    else
        s = (S + (Nμ-2F).*gmm.μ) ./ gmm.Σ
        return(N, f, s)
    end
end

## You can also get centered+scaled stats in a Cstats structure directly by
## using the constructor with a GMM argument
CSstats(gmm::GMM, x::DataOrMatrix) = CSstats(csstats(gmm, x, 1))

## centered stats, but not scaled by UBM covariance
## check full covariance...
function cstats(gmm::GMM, x::DataOrMatrix{T}, parallel=false) where T <: AbstractFloat
    nₓ, llh, N, F, S = stats(gmm, x, order=2, parallel=parallel)
    Nμ =  N .* gmm.μ
    ## center the statistics
    gmmkind = kind(gmm)
    if gmmkind == :diag
        S += (Nμ-2F) .* gmm.μ
    elseif gmmkind == :full
        for i in 1:length(S)
            μi = gmm.μ[i,:]
            Fμi = F[i,:]' * μi
            S[i] += N[i] * μi' * μi - Fμi' - Fμi
        end
    else
        error("Unknown kind")
    end
    F -= Nμ
    return N, F, S
end

Cstats(gmm::GMM, x::DataOrMatrix, parallel=false) = Cstats(cstats(gmm, x, parallel))

## conversion from Cstats to CSstats:
CSstats(gmm::GMM, cstats::Cstats) = CSstats(cstats.N, cstats.F ./ gmm.Σ)

## some convenience functions
Base.eltype(stats::Cstats{T}) where {T} = T
Base.size(stats::Cstats) = size(stats.F)
kind(stats::Cstats) = typeof(stats.S) <: Vector ? :full : :diag
Base.:+(a::Cstats, b::Cstats) = Cstats(a.N + b.N, a.F + b.F, a.S + b.S)
Base.zero(x::Cstats) = Cstats(zero(x.N), zero(x.F), zero(x.S))


text:
include("ieee1788-constructors.jl")
include("libieeep1788_tests_bool.jl")
include("libieeep1788_tests_cancel.jl")
include("libieeep1788_tests_elem.jl")
include("libieeep1788_tests_mul_rev.jl")
include("libieeep1788_tests_num.jl")
include("libieeep1788_tests_overlap.jl")
include("libieeep1788_tests_rec_bool.jl")
include("libieeep1788_tests_rev.jl")
include("libieeep1788_tests_set.jl")
include("mpfi.jl")
include("fi_lib.jl")
include("c-xsc.jl")


text:
import Pkg
Pkg.activate(".")
Pkg.instantiate()

using DelimitedFiles
using LinearAlgebra
using Printf
using SparseArrays

using DiffEqCallbacks
using ForwardDiff
using JSON
using OrdinaryDiffEq
using UnPack
using RecursiveArrayTools
using Roots
using StaticArrays
using SummationByPartsOperators
import SymPy; sp=SymPy

import FFTW; FFTW.set_num_threads(1)

using PyCall, LaTeXStrings; import PyPlot; plt=PyPlot
inset_locator = pyimport("mpl_toolkits.axes_grid.inset_locator")

cycler = pyimport("cycler").cycler
line_cycler   = (cycler(color=["#E69F00", "#56B4E9", "#009E73", "#0072B2", "#D55E00", "#CC79A7", "#F0E442"]) +
                 cycler(linestyle=["-", "--", "-.", ":", "-", "--", "-."]))
marker_cycler = (cycler(color=["#E69F00", "#56B4E9", "#009E73", "#0072B2", "#D55E00", "#CC79A7", "#F0E442"]) +
                 cycler(linestyle=["none", "none", "none", "none", "none", "none", "none"]) +
                 cycler(marker=["4", "2", "3", "1", "+", "x", "."]))

plt.rc("axes", prop_cycle=line_cycler)
plt.rc("text", usetex=true)
plt.rc("text.latex", preamble="\\usepackage{newpxtext}\\usepackage{newpxmath}\\usepackage{commath}\\usepackage{mathtools}")
plt.rc("font", family="serif", size=18.)
plt.rc("savefig", dpi=100)
plt.rc("legend", loc="best", fontsize="medium", fancybox=true, framealpha=0.5)
plt.rc("lines", linewidth=2.5, markersize=10, markeredgewidth=2.5)

function fit_error_growth(t, err)
    @assert axes(t) == axes(err)

    A = hcat(one.(t), t)
    c1, c2 = A \ log.(err)
    exp(c1), c2
end

relative_error(a, b) = (a - b) / max(abs(a), abs(b))

max_stars_repo_name: trixi-framework/paper-EC-KEP-PEP

max_issues_repo_name: trixi-framework/paper-EC-KEP-PEP

max_forks_repo_name: trixi-framework/paper-EC-KEP-PEP


text:
# Once we drop 0.4 support, we want the following to return depwarns
# info("Beginning of tests with deprecation warnings")

v = 1
d = Dict{Int,Int}()
d[1] = 1
@test Compat.@Dict(1 => 1) == d
@test Compat.@Dict(1 => v) == d

@test typeof(@compat(Dict(1 => 1))) == Dict{Int,Int}
@test @compat(Dict(1 => 1)) == d
@test @compat(Dict(1 => v)) == d

ad = Dict{Any,Any}()
ad[1] = 1
@test Compat.@AnyDict(1 => 1) == ad
@test Compat.@AnyDict(1 => v) == ad

@test typeof(@compat(Dict{Any,Any}(1 => 1))) == Dict{Any,Any}
@test @compat(Dict{Any,Any}(1 => 1)) == ad
@test @compat(Dict{Any,Any}(1 => v)) == ad

td = Dict{Int,Float64}()
td[1] = 1.0

@test typeof(@compat(Dict{Int,Float64}(1 => 1))) == Dict{Int,Float64}
@test @compat(Dict{Int,Float64}(1 => 1)) == td
@test @compat(Dict{Int,Float64}(1 => v)) == td

@test @compat(Dict()) == Dict()
@test @compat(Dict{Any,Any}()) == Dict{Any,Any}()
@test @compat(Dict([(1, 1)])) == d
@test @compat(Dict(:Void => :Nothing)) == Dict([(:Void, :Nothing)])

d2 = Dict{Symbol,Dict{Symbol,Int}}()
d2[:a] = Dict{Symbol,Int}()
d2[:a][:b] = 1
@test @compat(Dict(:a => Dict(:b => 1))) == d2

d = Dict(zip([1, 2], [3, 4]))
@test d == @compat Dict(1=>3, 2=>4)

@compat function f()
    a = :a
    b = Dict(:b => 1)
    Dict(a => b)
end
@test f() == d2

module CartesianTest
    using Base.Cartesian, Compat
    @ngenerate N NTuple{N,Int} function f(X::NTuple{N,Int}...)
        @ncall N tuple X
    end
end

@test CartesianTest.f(1) == (1,)
@test CartesianTest.f(1,2) == (1,2)
@test CartesianTest.f(1,2,3) == (1,2,3)
@test CartesianTest.f(1,2,3,4) == (1,2,3,4)
@test CartesianTest.f(1,2,3,4,5) == (1,2,3,4,5)


text:
module ConstantFieldsBenchs

using InplaceArrays.Arrays
using InplaceArrays.Fields

@inline function repeat(n,f,args...)
  for i in 1:n
    f(args...)
  end
  nothing
end

function bench1(n)
  v = 3.0
  f = v
  xi = Point(2,1)
  np = 4
  x = fill(xi,np)
  cf = field_cache(f,x)
  @time repeat(n,evaluate_field!,cf,f,x)
  @time repeat(n,field_gradient,f)
end

function bench2(n)
  v = 3.0
  f = [v,2*v,3*v]
  xi = Point(2,1)
  np = 4
  x = fill(xi,np)
  cf = field_cache(f,x)
  @time repeat(n,evaluate_field!,cf,f,x)
  @time repeat(n,field_gradient,f)
end

for n in (1,1,10,1000,100000)
  @eval begin
    println("+++ runing suite for n = $($n) +++")
    bench1($n)
    bench2($n)
  end
end

end # module


text:
# COUNTS LINE NUMBERS IN LIBRARY AND UPDATES THEM IN THE MAIN FILE

#
#   HEADER PART, DEFINITION OF MAIN FILE
#

# the main file
main_file_filename = "/home/jattig/PackageDevelopement/LatticePhysics.jl/src/LatticePhysics.jl"

# find the folder
main_file_folder = main_file_filename[1:findlast(main_file_filename, '/')]

# get all files within that folder
all_files = readdir(main_file_folder)

# open the main file
main_file = open(main_file_filename, "r")
# read the complete main file
main_file_lines = String[rstrip(line) for line in readlines(main_file)]
# close the main file
close(main_file)

#
#   EXTRACT NUMBER OF LINES FROM SUBFILES
#

# list of all line numbers
line_numbers = Int64[]

# list of all include lines (line numbers where include(...) is written)
include_lines = Int64[]

# process all lines and find include lines
for l in 5:length(main_file_lines)
    # look at line lines preceeding line l to find include lines
    if  main_file_lines[l-3] == "################################################################################" &&
        main_file_lines[l-2] == "" &&
        length(main_file_lines[l-1]) > 10 &&
        main_file_lines[l-1][1:10] == "# included"
        # this line contains the include
        push!(include_lines, l)
    end
end

# extract the included files from these lines
include_files = String[]
for line in [main_file_lines[i] for i in include_lines]
    # extract the filename
    filename = line[findfirst(line, '"')+1:findlast(line, '"')-1]
    # check if it is in the directory
    if !(filename in all_files)
        println("file not found: $(filename)")
        push!(include_files, "NOTHING.jl")
    else
        # push the filename into the list
        push!(include_files, filename)
    end
end

# extract line number information from the files
for include_filename in include_files
    # get the number of lines in the corresponding file
    f = open("$(main_file_folder)$(include_filename)", "r")
    push!(line_numbers, length(readlines(f)))
    close(f)
    # print the success
    println("$(include_filename)  ($(line_numbers[end]) lines)")
end

# overwrite the lines in the main file specifying the include comments with line number counts
for i in 1:length(include_lines)
    # get the line
    l = include_lines[i]
    count = line_numbers[i]
    # compose the new line
    line = "# included in subfile ($(count) lines)"
    # substitute the line
    main_file_lines[l-1] = line
end

# overwrite the total count string at the end of the file
main_file_lines[end-2] = "# total lines: $(length(main_file_lines))"
for l in line_numbers
    main_file_lines[end-2] = "$(main_file_lines[end-2]) + $(l)"
end

# add the total number
main_file_lines[end-1] = "# = $(sum(line_numbers) + length(main_file_lines)) lines"
# print the total number
println("  --> total: $(sum(line_numbers) + length(main_file_lines)) lines")

#
#   SET COMMENT BLOCKS CORRECTLY
#

# find all comment blocks in main file, included file and specify their data
# main_file_line_start, main_file_line_end, sub_file_lines
comment_blocks = Array{Any,1}[]

# find all comment blocks in subfiles
for (i,include_filename) in enumerate(include_files)
    # get the number of lines in the corresponding file
    f = open("$(main_file_folder)$(include_filename)", "r")
    subfile_lines = readlines(f)
    close(f)
    # find the extent of the comment block
    comment_length = 0
    while length(subfile_lines[comment_length+1]) > 0 && subfile_lines[comment_length+1][1] == '#'
        comment_length += 1
    end
    # add a new entry to the comment block list
    push!(comment_blocks, Any[-1, include_lines[i]-3, subfile_lines[1:comment_length]])
end

# find all comment blocks in main file
for (i,cb) in enumerate(comment_blocks)
    # check the extent
    extent = 0
    while length(main_file_lines[cb[2]-extent]) > 0 && main_file_lines[cb[2]-extent][1] == '#'
        extent += 1
    end
    # save into comment block info
    cb[1] = cb[2] - extent + 1
end

# print all comment blocks
for cb in comment_blocks
    println("from line $(cb[1]) to $(cb[2]):")
    for l in cb[3]
        println(l)
    end
end

# build a new list of all main_file lines
main_file_lines_new = String[]
for l in main_file_lines[1:comment_blocks[1][1]-1]
    push!(main_file_lines_new, l)
end
for i in 1:length(comment_blocks)-1
    for l in comment_blocks[i][3]
        push!(main_file_lines_new, l)
    end
    for l in main_file_lines[comment_blocks[i][2]+1:comment_blocks[i+1][1]-1]
        push!(main_file_lines_new, l)
    end
end
for l in comment_blocks[end][3]
    push!(main_file_lines_new, l)
end
for l in main_file_lines[comment_blocks[end][2]+1:end]
    push!(main_file_lines_new, l)
end

# overwrite original main file lines
main_file_lines = main_file_lines_new

# overwrite the total count string at the end of the file
main_file_lines[end-2] = "# total lines: $(length(main_file_lines))"
for l in line_numbers
    main_file_lines[end-2] = "$(main_file_lines[end-2]) + $(l)"
end

# add the total number
main_file_lines[end-1] = "# = $(sum(line_numbers) + length(main_file_lines)) lines"
# print the total number
println("  --> total: $(sum(line_numbers) + length(main_file_lines)) lines (corrected)")

#
#   LAST STEP: WRITE MAIN FILE AGAIN TO FILE
#

# write to file
open(main_file_filename, "w") do main_file
    for line in main_file_lines
        write(main_file, "$(line)\n")
    end
end


text:
module ObservablePmap

export opmap, ologpmap

include("opmap.jl")
include("ologpmap.jl")

end


text:
"""
    makedemos(source::String;
              root = "<current-directory>",
              src = "src",
              build = "build",
              branch = "gh-pages",
              edit_branch = "master",
              credit = true,
              throw_error = false) -> page_path, postprocess_cb

Make a demo page for `source` and return the path to the generated index file.

Processing pipeline:

1. analyze the folder structure `source` and loading all available configs.
2. copy assets
3. preprocess demo files and save it
4. save/copy cover images
5. generate postprocess callback function, which includes url-redirection.

!!! warning
    By default, `makedemos` generates all the necessary files to `docs/src/`, this means that the
    data you pass to `makedemos` should not be placed at `docs/src/`. A recommendation is to
    put folders in `docs/`. For example, `docs/quickstart` is a good choice.

# Inputs

* `source`: dir path to the root page of your demos; relative to `docs`.

# Outputs

* `page_path`: path to demo page's index. You can directly pass it to `makedocs`.
* `postprocess_cb`: callback function for postprocess. You can call `postprocess_cb()` _after_
  `makedocs`.
* `theme_assets`: the stylesheet assets that you may need to pass to `Documenter.HTML`. When the
  demo page has no theme, it will return `nothing`.

# Keywords

* `root::String`: should be equal to `Documenter`'s setting. Typically it is `"docs"` if this
  function is called in `docs/make.jl` file.
* `src::String`: should be equal to `Documenter`'s setting. By default it's `"src"`.
* `build::String`: should be equal to `Documenter`'s setting. By default it's `"build"`.
* `edit_branch::String`: should be equal to `Documenter`'s setting. By default it's `"master"`.
* `branch::String`: should be equal to `Documenter`'s setting. By default it's `"gh-pages"`.
* `credit::Bool`: `true` to show a "This page is generated by ..." info. By default it's `true`.
* `throw_error::Bool`: `true` to throw an error when the julia demo build fails; otherwise it will
  continue the build with warnings.

# Examples

The following is the simplest example for you to start with:

```julia
# 1. preprocess and generate demo files to docs/src
examples, postprocess_cb, demo_assets = makedemos("examples")

assets = []
isnothing(demo_assets) || (push!(assets, demo_assets))

# 2. do the standard Documenter pipeline
makedocs(format = Documenter.HTML(assets = assets),
         pages = [
             "Home" => "index.md",
             examples
         ])

# 3. postprocess after makedocs
postprocess_cb()
```

By default, it won't generate the index page for your demos. To enable it and configure how index
page is generated, you need to create "examples/config.json", whose contents should looks like the
following:

```json
{
    "theme": "grid",
    "order": [
        "basic",
        "advanced"
    ]
}
```

For more details of how `config.json` is configured, please check [`DemoCards.DemoPage`](@ref).
"""
function makedemos(source::String, templates::Union{Dict, Nothing} = nothing;
                   root::String = Base.source_dir(),
                   src::String = "src",
                   build::String = "build",
                   branch::String = "gh-pages",
                   edit_branch::String = "master",
                   credit = true,
                   throw_error = false)

    if !(basename(pwd()) == "docs" || basename(root) == "docs" || root == preview_build_dir())
        # special cases that warnings are not printed:
        # 1. called from `docs/make.jl`
        # 2. called from `preview_demos`
        # 3. pwd() in `docs` -- REPL
        @warn "Suspicious `root` setting: typically `makedemos` should be called from `docs/make.jl`. " pwd=pwd() root
    end

    root = abspath(root)
    page_root = abspath(root, source)

    if page_root == source
        # reach here when source is absolute path
        startswith(source, root) || throw(ArgumentError("invalid demo page source $page_root"))
    end

    page = DemoPage(page_root)

    relative_root = basename(page)

    # Where files are generated by DemoCards for Documenter.makedocs
    # e.g "$PROJECT_ROOT/docs/src/quickstart"
    absolute_root = abspath(root, src, relative_root)

    if page_root == absolute_root
        throw(ArgumentError("Demo page $page_root should not be placed at \"docs/$src\" folder."))
    end

    config_file = joinpath(page.root, config_filename)
    if !isnothing(templates)
        Base.depwarn(
            "It is not recommended to pass `templates`, please configure the index theme in the config file $(config_file). For example `\"theme\": \"grid\"`",
            :makedemos)
    end

    if !isnothing(page.theme)
        page_templates, theme_assets = cardtheme(page.theme, root = root)
        theme_assets = something(page.stylesheet, theme_assets)

        if templates != page_templates && templates != nothing
            @warn "use configured theme from $(config_file)"
        end

        templates = page_templates
    else
        theme_assets = page.stylesheet
    end

    # we can directly pass it to Documenter.makedocs
    if isnothing(templates)
        out_path = walkpage(page; flatten=false) do item
            splitext(joinpath(relative_root, relpath(item.path, page_root)))[1] * ".md"
        end
    else
        out_path = joinpath(relative_root, "index.md")
    end

    @info "SetupDemoCardsDirectory: setting up \"$(source)\" directory."
    if isdir(absolute_root)
        # a typical and probably safe case -- that we're still in docs/ folder
        trigger_prompt = !endswith(dirname(absolute_root), joinpath("docs", "src"))

        if trigger_prompt
            @warn "DemoCards build dir $absolute_root already exists.\nThis should only happen when you are using DemoCards incorrectly."
            # Should never reach this in CI environment
            clean_builddir = input_bool("It might contain important data, do you still want to remove it?")
        else
            clean_builddir = true
        end
        if clean_builddir
            @info "Remove DemoCards build dir: $absolute_root"
            # Ref: https://discourse.julialang.org/t/find-what-has-locked-held-a-file/23278/2
            Base.Sys.iswindows() && GC.gc()
            rm(absolute_root; force=true, recursive=true)
        else
            error("Stopped demos build.")
        end
    end

    mkpath(absolute_root)
    try
        # hard coded "covers" should be consistant to card template
        isnothing(templates) || mkpath(joinpath(absolute_root, "covers"))

        # make a copy before pipeline because `save_democards` modifies card path
        source_files = [x.path for x in walkpage(page)[2]]

        # pipeline
        # for themeless version, we don't need to generate covers and index page
        copy_assets(absolute_root, page)
        # WARNING: julia cards are reconfigured here
        save_democards(absolute_root, page;
                    project_root = root,
                    src = src,
                    credit = credit,
                    nbviewer_root_url = get_nbviewer_root_url(branch),
                    throw_error=throw_error)
        isnothing(templates) || save_cover(joinpath(absolute_root, "covers"), page)
        isnothing(templates) || generate(joinpath(absolute_root, "index.md"), page, templates)

        # pipeline: generate postprocess callback function
        postprocess_cb = ()->begin
            @info "Redirect page URL: redirect docs-edit-link for demos in \"$(source)\" directory."
            isnothing(templates) || push!(source_files, joinpath(page_root, "index.md"))
            foreach(source_files) do source_file
                redirect_link(source_file, source, root, src, build, edit_branch)
            end

            @info "Clean up DemoCards build dir: \"$source\""
            # Ref: https://discourse.julialang.org/t/find-what-has-locked-held-a-file/23278/2
            Base.Sys.iswindows() && GC.gc()
            rm(absolute_root; force=true, recursive=true)

            if !isnothing(theme_assets)
                assets_path = abspath(root, src, theme_assets)
                # Ref: https://discourse.julialang.org/t/find-what-has-locked-held-a-file/23278/2
                Base.Sys.iswindows() && GC.gc()
                rm(assets_path, force=true)
                if(isempty(readdir(dirname(assets_path))))
                    # Ref: https://discourse.julialang.org/t/find-what-has-locked-held-a-file/23278/2
                    Base.Sys.iswindows() && GC.gc()
                    rm(dirname(assets_path))
                end
            end
        end

        return out_path, postprocess_cb, theme_assets

    catch err
        # clean up build dir if anything wrong happens here so that we _hopefully_ never trigger the
        # user prompt logic before the build process.
        # Ref: https://discourse.julialang.org/t/find-what-has-locked-held-a-file/23278/2
        Base.Sys.iswindows() && GC.gc()
        rm(absolute_root; force=true, recursive=true)
        @error "Errors when building demo dir" pwd=pwd() source root src
        rethrow(err)
    end
end

function generate(file::String, page::DemoPage, args...)
    check_ext(file, :markdown)
    open(file, "w") do f
        generate(f, page::DemoPage, args...)
    end
end
generate(io::IO, page::DemoPage, args...) = write(io, generate(page, args...))

function generate(page::DemoPage, templates)
    items = Dict("democards" => generate(page.sections, templates; properties=page.properties))
    Mustache.render(page.template, items)
end

function generate(cards::AbstractVector{<:AbstractDemoCard}, template; properties=Dict{String, Any}())
    # for those hidden cards, only generate the necessary assets and files, but don't add them into
    # the index.md page
    foreach(filter(x->x.hidden, cards)) do x
        generate(x, template; properties=properties)
    end

    mapreduce(*, filter(x->!x.hidden, cards); init="") do x
        generate(x, template; properties=properties)
    end
end

function generate(secs::AbstractVector{DemoSection}, templates; level=1, properties=Dict{String, Any}())
    mapreduce(*, secs; init="") do x
        properties = merge(properties, x.properties) # sec.properties has higher priority
        generate(x, templates; level=level, properties=properties)
    end
end

function generate(sec::DemoSection, templates; level=1, properties=Dict{String, Any}())
    header = repeat("#", level) * " " * sec.title * "\n"
    footer = "\n"
    properties = merge(properties, sec.properties) # sec.properties has higher priority

    # either cards or subsections are empty
    # recursively generate the page contents
    if isempty(sec.cards)
        body = generate(sec.subsections, templates; level=level+1, properties=properties)
    else
        items = Dict(
            "cards" => generate(sec.cards, templates["card"], properties=properties),
            "description" => sec.description
        )
        body = Mustache.render(templates["section"], items)
    end
    header * body * footer
end

function generate(card::AbstractDemoCard, template; properties=Dict{String, Any})
    covername = get_covername(card)

    if isnothing(covername)
        # `generate` are called after `save_cover`, we assume that this default cover file is
        # already generated
        coverpath = "covers/" * basename(get_logopath())
    else
        coverpath = is_remote_url(card.cover) ? covername : "covers/" * covername
    end

    description = card.description
    cut_idx = 500
    if length(card.description) >= cut_idx
        # cut descriptions into ~500 characters
        offset = findfirst(' ', description[cut_idx:end])
        offset === nothing && (offset = 0)
        offset = cut_idx + offset - 2
        description = description[1:cut_idx] * "..."
    end

    items = Dict(
        "coverpath" => coverpath,
        "id" => card.id,
        "title" => card.title,
        "description" => description,
    )
    Mustache.render(template, items)
end

### save demo card covers

save_cover(path::String, page::DemoPage) = save_cover.(path, page.sections)
function save_cover(path::String, sec::DemoSection)
    save_cover.(path, sec.subsections)
    save_cover.(path, sec.cards)
end

"""
    save_cover(path::String, card::AbstractDemoCard)

process the cover image and save it.
"""
function save_cover(path::String, card::AbstractDemoCard)
    covername = get_covername(card)
    if isnothing(covername)
        default_coverpath = get_logopath()
        cover_path = joinpath(path, basename(default_coverpath))
        !isfile(cover_path) && cp(default_coverpath, cover_path)

        # now card uses default cover file as a fallback
        card.cover = basename(default_coverpath)
        return nothing
    end

    is_remote_url(card.cover) && return nothing

    # only save cover image if it is a existing local file
    src_path = joinpath(dirname(card.path), card.cover)
    if !isfile(src_path)
        @warn "cover file doesn't exists" cover_path=src_path

        # reset it back to nothing and fallback to use default cover
        card.cover = nothing
        save_cover(path, card)
    else
        cover_path = joinpath(path, covername)
        if isfile(cover_path)
            @warn "cover file already exists, perhaps you have demos of the same filename" cover_path
        end
        !isfile(cover_path) && cp(src_path, cover_path)
    end
end

function get_covername(card::AbstractDemoCard)
    isnothing(card.cover) && return nothing
    is_remote_url(card.cover) && return card.cover

    default_covername = basename(get_logopath())
    card.cover == default_covername && return default_covername

    return splitext(basename(card))[1] * splitext(card.cover)[2]
end

get_logopath() = joinpath(pkgdir(DemoCards), "assets", "democards_logo.svg")

### save markdown files

"""
    save_democards(root::String, page::DemoPage; credit, nbviewer_root_url)

recursively process and save source demo file
"""
function save_democards(root::String, page::DemoPage; kwargs...)
    save_democards.(root, page.sections; properties=page.properties, kwargs...)
end
function save_democards(root::String, sec::DemoSection; properties, kwargs...)
    properties = merge(properties, sec.properties) # sec.properties has higher priority
    save_democards.(joinpath(root, basename(sec.root)), sec.subsections;
                    properties=properties, kwargs...)
    save_democards.(joinpath(root, basename(sec.root)), sec.cards;
                    properties=properties, kwargs...)
end

### copy assets

function copy_assets(path::String, page::DemoPage)
    _copy_assets(dirname(path), page.root)
    copy_assets.(path, page.sections)
end
function copy_assets(path::String, sec::DemoSection)
    _copy_assets(path, sec.root)
    copy_assets.(joinpath(path, basename(sec.root)), sec.subsections)
end

function _copy_assets(dest_root::String, src_root::String)
    # copy assets of this section
    assets_dirs = filter(x->isdir(x)&&(basename(x) in ignored_dirnames),
                         joinpath.(src_root, readdir(src_root)))
    map(assets_dirs) do src
        dest = joinpath(dest_root, basename(src_root), basename(src))
        mkpath(dest)
        cp(src, dest; force=true)
    end
end

### postprocess

"""
    redirect_link(src_file, source, root, src, build, edit_branch)

Redirect the "Edit On GitHub" link of generated demo files to its original url, without
this a 404 error is expected.
"""
function redirect_link(source_file, source, root, src, build, edit_branch)
    build_file = get_build_file(source_file, source, build)
    if !isfile(build_file)
        @warn "$build_file doesn't exists, skip"
        return nothing
    end
    contents = read(build_file, String)

    if !isfile(source_file)
        # reach here when user doesn't create a page template index.md
        # just remove the whole button so that user don't click and get 404
        new_contents = replace(contents, regex_edit_on_github=>"")
    else
        # otherwise, redirect the url links
        m = match(regex_edit_on_github, contents)
        isnothing(m) && return nothing
        build_url = m.captures[1]

        src_url = get_source_url(build_url, source, basename(source_file), src)
        new_contents = replace(contents, build_url=>src_url)
    end
    write(build_file, new_contents)
end

function get_source_url(build_url, source, cardname, src)
    # given input:
    #   - projct_root:          "$REPO/blob/$edit_branch"
    #   - build_root:           "$projct_root/$docs_root/$src"
    #   - build_dir:            "$build_root/$prefix/$page/$section/$subsection"
    #   - build_url:            "$build_dir/$cardfile"
    # example of build_url:
    #  "https://github.com/johnnychen94/DemoCards.jl/blob/master/docs/src/quickstart/usage_example/julia_demos/2.cover_on_the_fly.md"
    # we need to generate:
    #   - src_root:             "$projct_root/$docs_root/$src"
    #   - src_dir:              "$src_root/$prefix/$page/$section/$subsection"
    #   - src_url:              "$src_dir/$cardfile"
    # example of src_url:
    #   "https://github.com/johnnychen94/DemoCards.jl/blob/master/docs/quickstart/usage_example/julia_demos/2.cover_on_the_fly.jl"
    source = replace(source, Base.Filesystem.path_separator => "/")

    repo, path = strip.(split(build_url, "/blob/"; limit=2), '/')
    root_to_subsection = replace(splitdir(path)[1], "$(src)/" => ""; count=1)
    root_to_subsection = replace(root_to_subsection, "/$(basename(source))/" => "/$(source)/"; count=1)
    src_url = join([repo, "blob", root_to_subsection, cardname], "/")

    return src_url
end

function get_build_file(source_file, source, build)
    # given inputs:
    #   - source_file: "$source_root/$prefix/$page/$section/$subsection/$card.md"
    #   - source:      "$prefix/$page
    #   - build:       "build"
    # we need to generate:
    #   - build_root: "$source_root/$build"
    #   - build_dir: "$build_root/$page/$section/$subsection
    #   - build_file: "$build_dir/$card.html" or "$build_dir/$card/index.html"

    sep = Base.Filesystem.path_separator
    # add trailing / to avoid incorrect substring match
    source_root = first(split(source_file, source * sep; limit=2))
    build_root = joinpath(source_root, build)
    prefix, page = splitdir(source)

    source_dir, name = splitdir(source_file)
    card, ext = splitext(name)
    _, prefix_to_subsection = split(source_dir, source_root; limit=2)
    if !isempty(prefix)
        # add trailing / to remove leading / for prefix_to_subsection
        # otherwise, joinpath of two absolute path would simply drop the first one
        _, prefix_to_subsection = split(prefix_to_subsection, prefix * sep; limit=2)
    end
    build_dir = joinpath(build_root, prefix_to_subsection)

    prettyurls = isdir(joinpath(build_dir, card))
    # Documenter.HTML behaves differently on prettyurls
    if prettyurls
        build_file = joinpath(build_dir, card, "index.html")
    else
        build_file = joinpath(build_dir, card * ".html")
    end
    return build_file
end

# modified from https://github.com/fredrikekre/Literate.jl to replace the use of @__NBVIEWER_ROOT_URL__
function get_nbviewer_root_url(branch)
    if haskey(ENV, "HAS_JOSH_K_SEAL_OF_APPROVAL") # Travis CI
        repo_slug = get(ENV, "TRAVIS_REPO_SLUG", "unknown-repository")
        deploy_folder = if get(ENV, "TRAVIS_PULL_REQUEST", nothing) == "false"
            tag = ENV["TRAVIS_TAG"]
            isempty(tag) ? "dev" : tag
        else
            "previews/PR$(get(ENV, "TRAVIS_PULL_REQUEST", "##"))"
        end
        return "https://nbviewer.jupyter.org/github/$(repo_slug)/blob/$(branch)/$(deploy_folder)"
    elseif haskey(ENV, "GITHUB_ACTIONS")
        repo_slug = get(ENV, "GITHUB_REPOSITORY", "unknown-repository")
        deploy_folder = if get(ENV, "GITHUB_EVENT_NAME", nothing) == "push"
            if (m = match(r"^refs\/tags\/(.*)$", get(ENV, "GITHUB_REF", ""))) !== nothing
                String(m.captures[1])
            else
                "dev"
            end
        elseif (m = match(r"refs\/pull\/(\d+)\/merge", get(ENV, "GITHUB_REF", ""))) !== nothing
            "previews/PR$(m.captures[1])"
        else
            "dev"
        end
        return "https://nbviewer.jupyter.org/github/$(repo_slug)/blob/$(branch)/$(deploy_folder)"
    elseif haskey(ENV, "GITLAB_CI")
        if (url = get(ENV, "CI_PROJECT_URL", nothing)) !== nothing
            cfg["repo_root_url"] = "$(url)/blob/$(devbranch)"
        end
        if (url = get(ENV, "CI_PAGES_URL", nothing)) !== nothing &&
           (m = match(r"https://(.+)", url)) !== nothing
            return "https://nbviewer.jupyter.org/urls/$(m[1])"
        end
    end
    return ""
end


text:
#!/usr/bin/env julia --project

using ArgParse
using LinearAlgebra
using StaticArrays
using Test

using ClimateMachine
using ClimateMachine.Atmos
using ClimateMachine.Orientations
using ClimateMachine.ConfigTypes
using ClimateMachine.Diagnostics
using ClimateMachine.GenericCallbacks
using ClimateMachine.ODESolvers
using ClimateMachine.TurbulenceClosures
using ClimateMachine.SystemSolvers: ManyColumnLU
using ClimateMachine.Mesh.Filters
using ClimateMachine.Mesh.Grids
using ClimateMachine.TemperatureProfiles
using ClimateMachine.Thermodynamics:
    air_density, air_temperature, total_energy, internal_energy, PhasePartition
using ClimateMachine.TurbulenceClosures
using ClimateMachine.VariableTemplates
import ClimateMachine.DGMethods.FVReconstructions: FVLinear

using CLIMAParameters
using CLIMAParameters.Planet: MSLP, R_d, day, grav, Omega, planet_radius
struct EarthParameterSet <: AbstractEarthParameterSet end
const param_set = EarthParameterSet()

function init_baroclinic_wave!(problem, bl, state, aux, localgeo, t)
    FT = eltype(state)

    # parameters
    param_set = parameter_set(bl)
    _grav::FT = grav(param_set)
    _R_d::FT = R_d(param_set)
    _Ω::FT = Omega(param_set)
    _a::FT = planet_radius(param_set)
    _p_0::FT = MSLP(param_set)

    k::FT = 3
    T_E::FT = 310
    T_P::FT = 240
    T_0::FT = 0.5 * (T_E + T_P)
    Γ::FT = 0.005
    A::FT = 1 / Γ
    B::FT = (T_0 - T_P) / T_0 / T_P
    C::FT = 0.5 * (k + 2) * (T_E - T_P) / T_E / T_P
    b::FT = 2
    H::FT = _R_d * T_0 / _grav
    z_t::FT = 15e3
    λ_c::FT = π / 9
    φ_c::FT = 2 * π / 9
    d_0::FT = _a / 6
    V_p::FT = 1
    M_v::FT = 0.608
    p_w::FT = 34e3                 # Pressure width parameter for specific humidity
    η_crit::FT = p_w / _p_0        # Critical pressure coordinate
    q_0::FT = 0.018                # Maximum specific humidity (default: 0.018)
    q_t::FT = 1e-12                # Specific humidity above artificial tropopause
    φ_w::FT = 2π / 9               # Specific humidity latitude wind parameter

    # grid
    φ = latitude(bl.orientation, aux)
    λ = longitude(bl.orientation, aux)
    z = altitude(bl.orientation, param_set, aux)
    r::FT = z + _a
    γ::FT = 1 # set to 0 for shallow-atmosphere case and to 1 for deep atmosphere case

    # convenience functions for temperature and pressure
    τ_z_1::FT = exp(Γ * z / T_0)
    τ_z_2::FT = 1 - 2 * (z / b / H)^2
    τ_z_3::FT = exp(-(z / b / H)^2)
    τ_1::FT = 1 / T_0 * τ_z_1 + B * τ_z_2 * τ_z_3
    τ_2::FT = C * τ_z_2 * τ_z_3
    τ_int_1::FT = A * (τ_z_1 - 1) + B * z * τ_z_3
    τ_int_2::FT = C * z * τ_z_3
    I_T::FT =
        (cos(φ) * (1 + γ * z / _a))^k -
        k / (k + 2) * (cos(φ) * (1 + γ * z / _a))^(k + 2)

    # base state virtual temperature, pressure, specific humidity, density
    T_v::FT = (τ_1 - τ_2 * I_T)^(-1)
    p::FT = _p_0 * exp(-_grav / _R_d * (τ_int_1 - τ_int_2 * I_T))

    # base state velocity
    U::FT =
        _grav * k / _a *
        τ_int_2 *
        T_v *
        (
            (cos(φ) * (1 + γ * z / _a))^(k - 1) -
            (cos(φ) * (1 + γ * z / _a))^(k + 1)
        )
    u_ref::FT =
        -_Ω * (_a + γ * z) * cos(φ) +
        sqrt((_Ω * (_a + γ * z) * cos(φ))^2 + (_a + γ * z) * cos(φ) * U)
    v_ref::FT = 0
    w_ref::FT = 0

    # velocity perturbations
    F_z::FT = 1 - 3 * (z / z_t)^2 + 2 * (z / z_t)^3
    if z > z_t
        F_z = FT(0)
    end
    d::FT = _a * acos(sin(φ) * sin(φ_c) + cos(φ) * cos(φ_c) * cos(λ - λ_c))
    c3::FT = cos(π * d / 2 / d_0)^3
    s1::FT = sin(π * d / 2 / d_0)
    if 0 < d < d_0 && d != FT(_a * π)
        u′::FT =
            -16 * V_p / 3 / sqrt(3) *
            F_z *
            c3 *
            s1 *
            (-sin(φ_c) * cos(φ) + cos(φ_c) * sin(φ) * cos(λ - λ_c)) /
            sin(d / _a)
        v′::FT =
            16 * V_p / 3 / sqrt(3) * F_z * c3 * s1 * cos(φ_c) * sin(λ - λ_c) /
            sin(d / _a)
    else
        u′ = FT(0)
        v′ = FT(0)
    end
    w′::FT = 0
    u_sphere = SVector{3, FT}(u_ref + u′, v_ref + v′, w_ref + w′)
    u_cart = sphr_to_cart_vec(bl.orientation, u_sphere, aux)

    if moisture_model(bl) isa DryModel
        q_tot = FT(0)
    else
        ## Compute moisture profile
        ## Pressure coordinate η
        ## η_crit = p_t / p_w ; p_t = 10000 hPa, p_w = 340 hPa
        η = p / _p_0
        if η > η_crit
            q_tot = q_0 * exp(-(φ / φ_w)^4) * exp(-((η - 1) * _p_0 / p_w)^2)
        else
            q_tot = q_t
        end
    end
    phase_partition = PhasePartition(q_tot)

    ## temperature & density
    T::FT = T_v / (1 + M_v * q_tot)
    ρ::FT = air_density(param_set, T, p, phase_partition)

    ## potential & kinetic energy
    e_pot::FT = gravitational_potential(bl.orientation, aux)
    e_kin::FT = 0.5 * u_cart' * u_cart
    e_tot::FT = total_energy(param_set, e_kin, e_pot, T, phase_partition)

    ## Assign state variables
    state.ρ = ρ
    state.ρu = ρ * u_cart
    state.energy.ρe = ρ * e_tot

    if !(moisture_model(bl) isa DryModel)
        state.moisture.ρq_tot = ρ * q_tot
    end

    nothing
end

function config_baroclinic_wave(
    FT,
    poly_order,
    fv_reconstruction,
    resolution,
    with_moisture,
)
    # Set up a reference state for linearization of equations
    temp_profile_ref =
        DecayingTemperatureProfile{FT}(param_set, FT(290), FT(220), FT(8e3))
    ref_state = HydrostaticState(temp_profile_ref; subtract_off = false)

    # Set up the atmosphere model
    exp_name = "BaroclinicWave"
    domain_height::FT = 30e3 # distance between surface and top of atmosphere (m)
    if with_moisture
        # hyperdiffusion = EquilMoistBiharmonic(FT(8 * 3600))
        moisture = EquilMoist()
        source = (Gravity(), Coriolis())
    else
        # hyperdiffusion = DryBiharmonic(FT(8 * 3600))
        moisture = DryModel()
        source = (Gravity(), Coriolis())
    end
    physics = AtmosPhysics{FT}(
        param_set;
        ref_state = ref_state,
        turbulence = ConstantKinematicViscosity(FT(0)),
        # hyperdiffusion = hyperdiffusion,
        moisture = moisture,
    )
    model = AtmosModel{FT}(
        AtmosGCMConfigType,
        physics;
        init_state_prognostic = init_baroclinic_wave!,
        source = source,
    )

    config = ClimateMachine.AtmosGCMConfiguration(
        exp_name,
        poly_order,
        resolution,
        domain_height,
        param_set,
        init_baroclinic_wave!;
        model = model,
        numerical_flux_first_order = RoeNumericalFlux(),
        fv_reconstruction = HBFVReconstruction(model, fv_reconstruction),
    )

    return config
end

function main()
    # add a command line argument to specify whether to use a moist setup
    # TODO: this will move to the future namelist functionality
    bw_args = ArgParseSettings(autofix_names = true)
    add_arg_group!(bw_args, "BaroclinicWave")
    @add_arg_table! bw_args begin
        "--with-moisture"
        help = "use a moist setup"
        action = :store_const
        constant = true
        default = false
    end

    cl_args = ClimateMachine.init(parse_clargs = true, custom_clargs = bw_args)
    with_moisture = cl_args["with_moisture"]

    # Driver configuration parameters
    FT = Float64                             # floating type precision
    poly_order = (5, 0)                      # discontinuous Galerkin polynomial order
    fv_reconstruction = FVLinear()           # finite volume reconstruction scheme
    n_horz = 8                               # horizontal element number
    n_vert = 20                              # vertical element number
    timestart::FT = 0                        # start time (s)
    timeend::FT = 3600                       # end time (s)

    # Set up driver configuration
    driver_config = config_baroclinic_wave(
        FT,
        poly_order,
        fv_reconstruction,
        (n_horz, n_vert),
        with_moisture,
    )

    # Set up experiment

    ode_solver_type = ClimateMachine.ExplicitSolverType(
        solver_method = LSRK54CarpenterKennedy,
    )

    CFL = FT(0.2) # target acoustic CFL number

    # time step is computed such that the horizontal acoustic Courant number is CFL
    solver_config = ClimateMachine.SolverConfiguration(
        timestart,
        timeend,
        driver_config,
        Courant_number = CFL,
        ode_solver_type = ode_solver_type,
        CFL_direction = EveryDirection(),
        diffdir = HorizontalDirection(),
    )

    # Set up diagnostics
    dgn_config = config_diagnostics(FT, driver_config)

    #TODO enable filter
    # # Set up user-defined callbacks
    filterorder = 20
    filter = ExponentialFilter(solver_config.dg.grid, 0, filterorder)
    cbfilter = GenericCallbacks.EveryXSimulationSteps(1) do
        Filters.apply!(
            solver_config.Q,
            AtmosFilterPerturbations(driver_config.bl),
            solver_config.dg.grid,
            filter,
            state_auxiliary = solver_config.dg.state_auxiliary,
            direction = HorizontalDirection(),
        )
        nothing
    end

    # cbtmarfilter = GenericCallbacks.EveryXSimulationSteps(1) do
    #     Filters.apply!(
    #         solver_config.Q,
    #         ("moisture.ρq_tot",),
    #         solver_config.dg.grid,
    #         TMARFilter(),
    #     )
    #     nothing
    # end

    # Run the model
    result = ClimateMachine.invoke!(
        solver_config;
        diagnostics_config = dgn_config,
        user_callbacks = (cbfilter,),
        #user_callbacks = (cbtmarfilter, cbfilter),
        check_euclidean_distance = true,
    )
end

function config_diagnostics(FT, driver_config)
    interval = "12shours" # chosen to allow diagnostics every 12 simulated hours

    _planet_radius = FT(planet_radius(param_set))

    info = driver_config.config_info
    boundaries = [
        FT(-90.0) FT(-180.0) _planet_radius
        FT(90.0) FT(180.0) FT(_planet_radius + info.domain_height)
    ]
    resolution = (FT(2), FT(2), FT(1000)) # in (deg, deg, m)
    interpol = ClimateMachine.InterpolationConfiguration(
        driver_config,
        boundaries,
        resolution,
    )

    dgngrp = setup_atmos_default_diagnostics(
        AtmosGCMConfigType(),
        interval,
        driver_config.name,
        interpol = interpol,
    )

    return ClimateMachine.DiagnosticsConfiguration([dgngrp])
end

main()


text:
include("product_of_sets.jl")
const GeometricConicForm{T} = MOI.Utilities.GenericModel{
    T,
    MOI.Utilities.ObjectiveContainer{Float64},
    MOI.Utilities.VariablesContainer{Float64},
    MOI.Utilities.MatrixOfConstraints{
        Float64,
        MOI.Utilities.MutableSparseMatrixCSC{
            Float64,
            Int,
            # We use `OneBasedIndexing` as it is the same indexing as used
            # by `SparseMatrixCSC` so we can do an allocation-free conversion to
            # `SparseMatrixCSC`.
            MOI.Utilities.OneBasedIndexing,
        },
        Vector{Float64},
        ProductOfSets{Float64},
    },
}

function build_conic_diff_cache!(model)
    # For theoretical background, refer Section 3 of Differentiating Through a Cone Program, https://arxiv.org/abs/1904.09043

    vis_src = MOI.get(model.optimizer, MOI.ListOfVariableIndices())
    cone_types = unique!([S for (F, S) in MOI.get(model.optimizer, MOI.ListOfConstraintTypesPresent())])
    conic_form = GeometricConicForm{Float64}()
    cones = conic_form.constraints.sets
    set_set_types(cones, cone_types)
    index_map = MOI.copy_to(conic_form, model)

    A = -convert(SparseMatrixCSC{Float64, Int}, conic_form.constraints.coefficients)
    b = conic_form.constraints.constants

    if MOI.get(model, MOI.ObjectiveSense()) == MOI.FEASIBILITY_SENSE
        c = spzeros(length(vis_src))
    else
        obj = MOI.get(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}())
        c = sparse_array_representation(obj, length(vis_src), index_map).terms
        if MOI.get(model, MOI.ObjectiveSense()) == MOI.MAX_SENSE
            c = -c
        end
    end

    # programs in tests were cross-checked against `diffcp`, which follows SCS format
    # hence, some arrays saved during `MOI.optimize!` are not same across all optimizers
    # specifically there's an extra preprocessing step for `PositiveSemidefiniteConeTriangle` constraint for SCS/Mosek

    # get x,y,s
    x = zeros(length(vis_src))
    for vi in vis_src
        i = index_map[vi].value
        x[i] = MOI.get(model, MOI.VariablePrimal(), vi)
    end
    s = map_rows((ci, r) -> MOI.get(model, MOI.ConstraintPrimal(), ci),
        model.optimizer, cones, index_map, Flattened{Float64}())
    y = map_rows((ci, r) -> MOI.get(model, MOI.ConstraintDual(), ci),
        model.optimizer, cones, index_map, Flattened{Float64}())

    # pre-compute quantities for the derivative
    m = A.m
    n = A.n
    N = m + n + 1
    # NOTE: w = 1.0 systematically since we asserted the primal-dual pair is optimal
    (u, v, w) = (x, y - s, 1.0)

    # find gradient of projections on dual of the cones
    Dπv = Dπ(v, model.optimizer, cones, index_map)

    # Q = [
    #      0   A'   c;
    #     -A   0    b;
    #     -c' -b'   0;
    # ]
    # M = (Q- I) * B + I
    # with B =
    # [
    #  I    .   .    # Πx = x because x is a solution and hence satistfies the constraints
    #  .  Dπv   .
    #  .    .   1    # w >= 0, but in the solution x = 1
    # ]
    # see: https://stanford.edu/~boyd/papers/pdf/cone_prog_refine.pdf
    # for the definition of Π and why we get I and 1 for x and w respectectively
    # K is defined in (5), Π in sect 2, and projection sin sect 3

    M = [
        spzeros(n,n)     (A' * Dπv)    c
        -A               -Dπv + I      b
        -c'              -b' * Dπv     0.0
    ]
    # find projections on dual of the cones
    vp = π(v, model.optimizer, cones, index_map)

    model.gradient_cache = ConicCache(
        M = M,
        vp = vp,
        Dπv = Dπv,
        xys = (x, y, s),
        A = A,
        b = b,
        c = c,
        index_map = index_map,
        cones = cones,
    )
    return nothing
end


text:
using RecipesBase

@recipe function f(model::PSFModel, inds...)
    seriestype := :heatmap
    aspect_ratio --> 1
    xlims --> extrema(last(inds))
    ylims --> extrema(first(inds))

    arr = model[inds...]

    return reverse(inds)..., arr
end

@recipe f(model::PSFModel, inds=axes(model)) = model, inds...


text:
# 10001st prime

using Primes

function solve()
    return prime(10001)
end

print(Int(solve()))

max_stars_repo_name: wephy/project-euler

max_issues_repo_name: wephy/project-euler

max_forks_repo_name: wephy/project-euler


text:
using BinaryBuilder

# Collection of sources required to build LLVM
sources = [
    "patches"
]

# Bash recipe for building across all platforms
script = raw"""
cd /workspace/destdir
rm bin/* || true
cp tools/llc* bin
cp tools/lld* bin
cp tools/LLVM.dll bin || true
mkdir newlib
cp -av lib/libLLVM*.so newlib || true
cp -av lib/libLLVM*.dylib newlib || true
rm -r lib
rm -r tools
rm -r include
rm -r share
rm -r libexec
rm -r logs
mv newlib lib
# empty out the manifests
> manifests/LLVM.${target}.list

"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line.
platforms = [
    Windows(:x86_64),
    Linux(:x86_64, :glibc),
    MacOS(),
    Windows(:i686),
    Linux(:i686, :glibc),
    Linux(:aarch64, :glibc),
    Linux(:armv7l, :glibc),
    # ppc64le isn't working right now....
    #Linux(:powerpc64le, :glibc),
]

# The products that we will ensure are always built
products(prefix) = [
    LibraryProduct(prefix, ["libLLVM"], :libLLVM),
    ExecutableProduct(prefix, "llc", :llc),
    ExecutableProduct(prefix, "lld", :lld),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    # "https://github.com/staticfloat/LLVMBuilder/releases/download/v6.0.0-3/build.jl"
    BinaryBuilder.InlineBuildDependency(raw"""
using BinaryProvider # requires BinaryProvider 0.3.0 or later

# Parse some basic command-line arguments
const verbose = "--verbose" in ARGS
const prefix = Prefix(get([a for a in ARGS if a != "--verbose"], 1, joinpath(@__DIR__, "usr")))
products = [
    LibraryProduct(prefix, String["libLLVM"], :libLLVM),
    LibraryProduct(prefix, String["libLTO"], :libLTO),
    LibraryProduct(prefix, String["libclang"], :libclang),
    ExecutableProduct(prefix, "llvm-config", :llvm_config),
]

# Download binaries from hosted location
bin_prefix = "https://github.com/staticfloat/LLVMBuilder/releases/download/v6.0.0-3"

# Listing of files generated by BinaryBuilder:
download_info = Dict(
    Linux(:aarch64, :glibc) => ("$bin_prefix/LLVM.aarch64-linux-gnu.tar.gz", "2f04a47e0b3b8205df963886aeebcfd1d1b7c985f9332a267f66af40339c1f35"),
    Linux(:armv7l, :glibc, :eabihf) => ("$bin_prefix/LLVM.arm-linux-gnueabihf.tar.gz", "a18f0d2c633995fcc68aeecb9a405f9f40293c812556057a1b8125063eea303c"),
    Linux(:i686, :glibc) => ("$bin_prefix/LLVM.i686-linux-gnu.tar.gz", "14f5dff6d707c4a2b71252d00e01a88723f0d14a89c971a7857f005e678b1507"),
    Windows(:i686) => ("$bin_prefix/LLVM.i686-w64-mingw32.tar.gz", "e106a2211236083d6d8aa7be378832f33c9c00dc857db26aad083b97d5272cec"),
    Linux(:powerpc64le, :glibc) => ("$bin_prefix/LLVM.powerpc64le-linux-gnu.tar.gz", "63c8456c20f8cfd9d9ff0b2ea75ee7bb97c87929415a2d7d94d90c710fdc3c2e"),
    MacOS(:x86_64) => ("$bin_prefix/LLVM.x86_64-apple-darwin14.tar.gz", "3fd47e1d52416f930667e98bdddf93e46ec242326c87887fa295a7981d692e25"),
    Linux(:x86_64, :glibc) => ("$bin_prefix/LLVM.x86_64-linux-gnu.tar.gz", "10428267da56b92089e212badf67ccf4c64097904344f3c68024383e65615398"),
    Windows(:x86_64) => ("$bin_prefix/LLVM.x86_64-w64-mingw32.tar.gz", "1241e6862b9c2797263d5337522c9937e0fcc67d63bfe0a1d108e9106e4cc29c"),
)

# Install unsatisfied or updated dependencies:
unsatisfied = any(!satisfied(p; verbose=verbose) for p in products)
if haskey(download_info, platform_key())
    url, tarball_hash = download_info[platform_key()]
    if unsatisfied || !isinstalled(url, tarball_hash; prefix=prefix)
        # Download and install binaries
        install(url, tarball_hash; prefix=prefix, force=true, verbose=verbose)
    end
elseif unsatisfied
    # If we don't have a BinaryProvider-compatible .tar.gz to download, complain.
    # Alternatively, you could attempt to install from a separate provider,
    # build from source or something even more ambitious here.
    error("Your platform $(triplet(platform_key())) is not supported by this package!")
end

# Write out a deps.jl file that will contain mappings for our products
write_deps_file(joinpath(@__DIR__, "deps.jl"), products)
""")
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, "lld", sources, script, platforms, products, dependencies)

max_stars_repo_name: tshort/LLDBuilder

max_issues_repo_name: tshort/LLDBuilder

max_forks_repo_name: tshort/LLDBuilder


text:
using Statistics
using Plots
pyplot()

using SIMD
using BenchmarkTools
using ChipSort
using Unrolled
using JLD

data = load("chip-medium-bench.jld")
bmk = data["bmk"]

T=Dict("Int32"=>UInt32, "Int64"=>UInt64, "UInt32"=>UInt32, "UInt64"=>UInt64)[data["T"]]
V=data["V"]
J=data["J"]

plot()
sizes = sort!(collect(keys(bmk)))
for srt in [:CombSortMix, :JuliaStd, :ChipSortM, :ChipSortL, :InsertionSort]
    xa = hcat(
        filter(
            x->x[2]<0.3 && x[1]<=2^20,
        [[l, median(bmk[l][srt].times * 1e-9)] for l in sizes if srt in keys(bmk[l])]
    )...)

    @show xa
    if length(xa)==0
        continue
    end

    lab = if srt == :ChipSortM "ChipSort-M"
    elseif srt == :ChipSortL "ChipSort-L"
    elseif srt == :CombSortMix "ChipSort"
    else srt end
    # plot!(xa[1,:], xa[2,:], l=2, m=4, label=srt)
    plot!(xa[1,:], xa[2,:] ./ xa[1,:], l=2, m=4, label=lab)
end

srt= :JuliaStd
xa = hcat(
        filter(
            x->x[2]<0.3 && x[1]<=2^20,
            [[l, median(bmk[l][srt].times * 1e-9)] for l in sizes if srt in keys(bmk[l])]
        )...)
plot!(xa[1,:], log.(xa[1,:]) * 5e-9, l=1, label="n⋅log(n)", color=:black, ls=:dash)

cache = 2^10 * [32, 256, 1024+512, 4*(1024+512)] * (1/sizeof(T))
for c in 1:3
    # plot!([cache[c], cache[c]], [10^4, 10^12], color=:black, ls=:dot)
    plot!([cache[c], cache[c]], [1e-8, 1e-7], color=:black, ls=:dot,label="")
end

plot!(xaxis=:log, yaxis=:log,
      yticks=([(1:10)*1e-8;], [(1:10);]),
      xticks=(2 .^(6:20), 2 .^(6:20)),
      title="Performance of sorting methods relative to input size",
      xlabel="Input size [n]",
      ylabel="Time / input size [1E8 × s/n]",
      size=(800,450)
      )

savefig("chipsort-bench-curves.pdf")
savefig("chipsort-bench-curves.svg")
savefig("chipsort-bench-curves.png")

plot(size=(800,600))
sizes = sort!(collect(keys(bmk)))
for srt in [:CombSortMix, :JuliaStd, :ChipSortM, :ChipSortL, :InsertionSort]
    xa = hcat(
        filter(
            x->x[2]<0.3 && x[1]<=2^20,
        [[l, median(bmk[l][srt].times * 1e-9)] for l in sizes if srt in keys(bmk[l])]
    )...)

    @show xa
    if length(xa)==0
        continue
    end

    lab = if srt == :ChipSortM "ChipSort-M"
    elseif srt == :ChipSortL "ChipSort-L"
    elseif srt == :CombSortMix "ChipSort"
    else srt end
    plot!(xa[1,:], xa[2,:], l=2, m=4, label=lab)
end

srt= :JuliaStd
xa = hcat(
        filter(
            x->x[2]<0.3 && x[1]<=2^20,
            [[l, median(bmk[l][srt].times * 1e-9)] for l in sizes if srt in keys(bmk[l])]
        )...)
plot!(xa[1,:], xa[1,:] .* log.(xa[1,:]) * 5e-9, l=1, label="n⋅log(n)", color=:black, ls=:dash)

plot!(xaxis=:log, yaxis=:log,
      # yticks=([(1:10)*1e-8;], [(1:10);]),
      xticks=(2 .^(6:20), 2 .^(6:20)),
      title="Running time of sorting methods",
      xlabel="Input size [n]",
      ylabel="Time [s]"
      )

savefig("chipsort-bench-time.pdf")
savefig("chipsort-bench-time.png")
savefig("chipsort-bench-time.svg")


text:
import Base: *, -

*(f, g::AbstractGP) = CompositeGP((*, f, g), g.gpc)
*(f::AbstractGP, g) = CompositeGP((*, g, f), f.gpc)
*(f::AbstractGP, g::AbstractGP) = throw(ArgumentError("Cannot multiply two GPs together."))

const prod_args{Tf} = Tuple{typeof(*), Tf, <:AbstractGP}

#
# Scale by a function
#

mean_vector((_, σ, g)::prod_args, x::AV) = σ.(x) .* mean_vector(g, x)

function cov((_, σ, g)::prod_args, x::AV)
    σx = σ.(x)
    return σx .* cov(g, x) .* σx'
end
cov_diag((_, σ, g)::prod_args, x::AV) = σ.(x).^2 .* cov_diag(g, x)

function cov((_, σ, g)::prod_args, x::AV, x′::AV)
    return σ.(x) .* cov(g, x, x′) .* σ.(x′)'
end
function cov_diag((_, σ, g)::prod_args, x::AV, x′::AV)
    return σ.(x) .* cov_diag(g, x, x′) .* σ.(x′)
end

cov((_, σ, f)::prod_args, f′::AbstractGP, x::AV, x′::AV) = σ.(x) .* cov(f, f′, x, x′)
cov(f::AbstractGP, (_, σ, f′)::prod_args, x::AV, x′::AV) = cov(f, f′, x, x′) .* (σ.(x′))'

function cov_diag((_, σ, f)::prod_args, f′::AbstractGP, x::AV, x′::AV)
    return σ.(x) .* cov_diag(f, f′, x, x′)
end
function cov_diag(f::AbstractGP, (_, σ, f′)::prod_args, x::AV, x′::AV)
    return cov_diag(f, f′, x, x′) .* σ.(x′)
end

#
# Scale by a constant
#

mean_vector((_, σ, g)::prod_args{<:Real}, x::AV) = σ .* mean_vector(g, x)

cov((_, σ, g)::prod_args{<:Real}, x::AV) = (σ^2) .* cov(g, x)
cov_diag((_, σ, g)::prod_args{<:Real}, x::AV) = (σ^2) .* cov_diag(g, x)

cov((_, σ, g)::prod_args{<:Real}, x::AV, x′::AV) = (σ^2) .* cov(g, x, x′)
cov_diag((_, σ, g)::prod_args{<:Real}, x::AV, x′::AV) = (σ^2) .* cov_diag(g, x, x′)

cov((_, σ, f)::prod_args{<:Real}, f′::AbstractGP, x::AV, x′::AV) = σ .* cov(f, f′, x, x′)
cov(f::AbstractGP, (_, σ, f′)::prod_args{<:Real}, x::AV, x′::AV) = cov(f, f′, x, x′) .* σ

function cov_diag((_, σ, f)::prod_args{<:Real}, f′::AbstractGP, x::AV, x′::AV)
    return σ .* cov_diag(f, f′, x, x′)
end
function cov_diag(f::AbstractGP, (_, σ, f′)::prod_args{<:Real}, x::AV, x′::AV)
    return cov_diag(f, f′, x, x′) .* σ
end

# Use multiplication to define the negation of a GP
-(f::AbstractGP) = (-1) * f


text:
abstract type WeightedOnlineStat{T} <: OnlineStat{T} end
using OnlineStats: fweights

meanweight(o::WeightedOnlineStat) = o.W
weightsum(o::WeightedOnlineStat) = meanweight(o) * nobs(o)
Base.eltype(o::WeightedOnlineStat{T}) where T = T

##############################################################
# Define our own interface so that it accepts two inputs.
##############################################################

function OnlineStatsBase.fit!(o::WeightedOnlineStat, xi::Number, wi::Number)
    _fit!(o, xi, wi)
    return o
end
OnlineStatsBase.fit!(o::WeightedOnlineStat, xi::Missing, wi::Number) = o
OnlineStatsBase.fit!(o::WeightedOnlineStat, xi::Number, wi::Missing) = o

# The missing cases in x are dealt with in the dispatch of _fit!
function OnlineStatsBase.fit!(o::WeightedOnlineStat{VectorOb}, x::VectorOb, w::Number)
    _fit!(o, x, w)
    return o
end
OnlineStatsBase.fit!(o::WeightedOnlineStat{VectorOb}, xi::VectorOb, wi::Missing) = o

function OnlineStatsBase.fit!(o::WeightedOnlineStat, y, w::AbstractVector)
    for (yy, ww) in zip(y, w)
        fit!(o, yy, ww)
    end
    o
end

OnlineStatsBase.fit!(o::WeightedOnlineStat, x::TwoThings) = fit!(o, x[1], x[2])
OnlineStatsBase.fit!(o::WeightedOnlineStat{VectorOb}, x::AbstractMatrix, w::AbstractVector) =
    fit!(o, eachrow(x), w)

function OnlineStatsBase.merge!(o::WeightedOnlineStat, o2::WeightedOnlineStat)
    (weightsum(o) > 0 || weightsum(o2) > 0) && _merge!(o, o2)
    o
end
function Base.show(io::IO, o::WeightedOnlineStat)
    print(io, name(o, false, false), ": ")
    print(io, "∑wᵢ=")
    show(IOContext(io, :compact => true), weightsum(o))
    print(io, " | value=")
    show(IOContext(io, :compact => true), value(o))
end


text:
using ElectronDisplay
using Gadfly
using Base.Test

@testset "ElectronDisplay" begin

p1 = plot(x=rand(10),y=rand(10),Geom.point)

f = display(p1)

@test f isa Electron.Window

end


text:
# Autogenerated wrapper script for LibSSH2_jll for armv6l-linux-musleabihf
export libssh2

using MbedTLS_jll
JLLWrappers.@generate_wrapper_header("LibSSH2")
JLLWrappers.@declare_library_product(libssh2, "libssh2.so.1")
function __init__()
    JLLWrappers.@generate_init_header(MbedTLS_jll)
    JLLWrappers.@init_library_product(
        libssh2,
        "lib/libssh2.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()