text:
#################### Gelman, Rubin, and Brooks Diagnostics ####################

function _gelmandiag(psi::AbstractArray{<:Real,3}; alpha::Real=0.05)
    niters, nparams, nchains = size(psi)
    nchains > 1 || error("Gelman diagnostic requires at least 2 chains")

    rfixed = (niters - 1) / niters
    rrandomscale = (nchains + 1) / (nchains * niters)

    S2 = map(Statistics.cov, (view(psi, :, :, i) for i in axes(psi, 3)))
    W = Statistics.mean(S2)

    psibar = dropdims(Statistics.mean(psi; dims=1); dims=1)'
    B = niters .* Statistics.cov(psibar)

    w = LinearAlgebra.diag(W)
    b = LinearAlgebra.diag(B)
    s2 = mapreduce(LinearAlgebra.diag, hcat, S2)'
    psibar2 = vec(Statistics.mean(psibar; dims=1))

    var_w = vec(Statistics.var(s2; dims=1)) ./ nchains
    var_b = (2 / (nchains - 1)) .* b .^ 2
    var_wb =
        (niters / nchains) .* (
            LinearAlgebra.diag(Statistics.cov(s2, psibar .^ 2)) .-
            2 .* psibar2 .* LinearAlgebra.diag(Statistics.cov(s2, psibar))
        )

    V = @. rfixed * w + rrandomscale * b
    var_V = rfixed^2 * var_w + rrandomscale^2 * var_b + 2 * rfixed * rrandomscale * var_wb

    df = @. 2 * V^2 / var_V

    B_df = nchains - 1
    W_df = @. 2 * w^2 / var_w

    estimates = Array{Float64}(undef, nparams)
    upperlimits = Array{Float64}(undef, nparams)

    q = 1 - alpha / 2
    for i in 1:nparams
        correction = (df[i] + 3) / (df[i] + 1)
        rrandom = rrandomscale * b[i] / w[i]

        estimates[i] = sqrt(correction * (rfixed + rrandom))

        if !isnan(rrandom)
            rrandom *= Distributions.quantile(Distributions.FDist(B_df, W_df[i]), q)
        end
        upperlimits[i] = sqrt(correction * (rfixed + rrandom))
    end

    return estimates, upperlimits, W, B
end

"""
    gelmandiag(chains::AbstractArray{<:Real,3}; alpha::Real=0.95)

Compute the Gelman, Rubin and Brooks diagnostics.
"""
function gelmandiag(chains::AbstractArray{<:Real,3}; kwargs...)
    estimates, upperlimits = _gelmandiag(chains; kwargs...)

    return (psrf=estimates, psrfci=upperlimits)
end

"""
    gelmandiag_multivariate(chains::AbstractArray{<:Real,3}; alpha::Real=0.05)

Compute the multivariate Gelman, Rubin and Brooks diagnostics.
"""
function gelmandiag_multivariate(chains::AbstractArray{<:Real,3}; kwargs...)
    niters, nparams, nchains = size(chains)
    if nparams < 2
        error(
            "computation of the multivariate potential scale reduction factor requires ",
            "at least two variables",
        )
    end

    estimates, upperlimits, W, B = _gelmandiag(chains; kwargs...)

    # compute multivariate potential scale reduction factor (PSRF)
    # the eigenvalues of `X := W⁻¹ B` and `Y := L⁻¹ B L⁻ᵀ = L⁻¹ Bᵀ L⁻ᵀ = L⁻¹ (L⁻¹ B)ᵀ`,
    # where `W = L Lᵀ`, are identical but `Y` is symmetric whereas `X` is not in general
    # (remember, `W` and `B` are symmetric positive semi-definite matrices)
    # for symmetric matrices specialized implementations for computing eigenvalues are used
    rfixed = (niters - 1) / niters
    rrandomscale = (nchains + 1) / (nchains * niters)
    C = LinearAlgebra.cholesky(LinearAlgebra.Symmetric(W))
    L = C.L
    Y = L \ (L \ LinearAlgebra.Symmetric(B))'
    λmax = LinearAlgebra.eigmax(LinearAlgebra.Symmetric(Y))
    multivariate = rfixed + rrandomscale * λmax

    return (psrf=estimates, psrfci=upperlimits, psrfmultivariate=multivariate)
end


text:
using Random
using Statistics
using StatsBase
using Plots
import Base.+
import Base.-

include("Point_Struct.jl")
include("fun.jl")
include("constants.jl")

# This script reproduce the turbulance distribution function (a gaussian)
# Figure 4 in Errard's article

p_0 = [Point{Float64}(20000, 0.9, i, undef, undef, undef)() for i = (1:25000)*( (2*π) / 25000  )  ] # must be constant
p_1 = [Point{Float64}(20000, i, 0.9, undef, undef, undef)() for i = (1:25000)*( (2*π) / 25000 )  ] # must be constant
p_2 = [Point{Float64}(i, 0.9, 0.9, undef, undef, undef)() for i = (10:1000000)/1000000 * 100000  ] # must be constant

w   = [ω(o) for o in p_2]
BB = [B(i) for i in p_2]

z   = [i.z for i in p_2]

plot(  BB[2:length(BB)] / (2*π*(ω₀^2)) , z[2:length(z)], yaxis=:log)#, xlims=(0,1050))
#plot( w[2:length(z)] .* ( (ω₀ * π) / ( λ * 100000 )), z[2:length(z)])#, yaxis=:log)#, xlims=(0,1050))

ext: jl

max_stars_repo_name: gruffoc/AFO

max_issues_repo_name: gruffoc/AFO

max_forks_repo_name: gruffoc/AFO


text:
# convert_a3d

# Method that allows federation by setting the `output_format`  in the Stanmodel().

"""

# convert_a3d

Convert the output file created by cmdstan to the shape of choice.

### Method
```julia
convert_a3d(a3d_array, cnames, ::Val{Symbol}; start=1)
```
### Required arguments
```julia
* `a3d_array::Array{Float64, 3},`      : Read in from output files created by cmdstan
* `cnames::Vector{AbstractString}`     : Monitored variable names
```

### Optional arguments
```julia
* `::Val{Symbol}`                      : Output format, default is :mcmcchains
* `::start=1`                          : First draw for MCMCChains.Chains
```
Method called is based on the output_format defined in the stanmodel, e.g.:

   stanmodel = Stanmodel(`num_samples`=1200, thin=2, name="bernoulli",
     model=bernoullimodel, `output_format`=:mcmcchains);

Current formats supported are:

1. :array (a3d_array format, the default for CmdStan)
2. :mcmcchains (MCMCChains.Chains object)
3. :dataframes (DataFrames.DataFrame object)

Option 2 is provided by the package StanMamba, option 3 is provided by StanMCMCChains.
```

### Return values
```julia
* `res`                       : Draws converted to the specified format.
```
"""
convert_a3d(a3d_array, cnames, ::Val{:array}; start=1) = a3d_array

"""

# convert_a3d

# Convert the output file created by cmdstan to a DataFrame.

$(SIGNATURES)

"""
function convert_a3d(a3d_array, cnames, ::Val{:dataframes}; start=1)
  snames = [Symbol(cnames[i]) for i in 1: length(cnames)]
  [DataFrame(a3d_array[:,:,i], snames) for i in 1:size(a3d_array, 3)]
end


text:
using StatsBase
include("toolbox.jl")

function augmentation_rule1(phi::Array{T,1},xi::Array{T,1},Xflat::Array{T,1},Y::Array{T,1},W::ActiveSet,Wlen::Int64,lamb::LambdaActiveSet,n::Int64,d::Int64,
        viol::Array{T,1},sort_index::Array{Int64,1},all_entries::Array{Int64,1}, entries::Array{Int64,1};
        TOL=1e-3::T,P=1::Int64,block=0::Int64) where {T <:  AbstractFloat}
    fill!(viol,0)
    fill!(entries,0)
    @simd for i in 2:n
        @inbounds all_entries[i-1]=i
    end
    P1::Int64 = 0;
    len::Int64 = 0;
    p::Int64=0;
    for i in 1:n
        len = 0
        if i != 1
            all_entries[i-1] -= 1
        end
        for j in all_entries
            if ~W[j,i]
                len += 1
                entries[len] = j
            end
        end
        A_dot_phi_plus_B_dot_xi!(viol,phi,xi,Xflat,i,entries,len,n,d,block=block);
        P1 = min(sum(viol.<-TOL),P);
        if P1 == 0
            continue
        elseif P1 == 1
            _,p = findmin(viol);
            W[entries[p],i] = true
            lamb[entries[p],i] = -1e-16
            lamb[entries[p],i] = 0.0
            Wlen += 1
        else
            sortperm!(sort_index,viol;alg=PartialQuickSort(P1))
            for k in 1:P1
                W[entries[sort_index[k]],i] = true
                lamb[entries[sort_index[k]],i] = -1e-16
                lamb[entries[sort_index[k]],i] = 0.0
                Wlen += 1
            end
        end
    end
    return Wlen
end

function augmentation_rule2(phi::Array{T,1},xi::Array{T,1},Xflat::Array{T,1},Y::Array{T,1},W::ActiveSet,Wlen::Int64,lamb::LambdaActiveSet,
        n::Int64,d::Int64,samples::Array{Int64,1};TOL=1e-3,K=n,block=0) where T<:AbstractFloat
    i::Int64 = 0
    j::Int64 = 0
    viol::Float64 = 0.
    @assert length(samples)==K
    sample!(1:n*(n-1), samples,replace=false)
    for k in samples
        (i,j) = idx2pair(k, n)
        if block == 0 && W[j,i]
            continue
        elseif block == 1 && W[i,j]
            continue
        end

        @inbounds viol = phi[j] - phi[i] + dot(Xflat[index2d(i,d)]-Xflat[index2d(j,d)], xi[index2d(i,d)])
        if viol < -TOL
            if block == 0
                W[j,i] = true
                lamb[j,i] = -1e-16
                lamb[j,i] = 0.
                Wlen += 1
            else
                W[i,j] = true
                lamb[i,j] = -1e-16
                lamb[i,j] = 0.
                Wlen += 1
            end
        end
    end
    return Wlen
end

function augmentation_rule3(phi::Array{T,1},xi::Array{T,1},Xflat::Array{T,1},Y::Array{T,1},W::ActiveSet,Wlen::Int64,lamb::LambdaActiveSet,
        n::Int64,d::Int64,all_entries::Array{Int64,1}, samples::Array{Int64,1};TOL=1e-3,P=1,block=0) where T<:AbstractFloat
    i::Int64 = 0
    j::Int64 = 0
    viol::Float64 = 0.
    @simd for i in 2:n
        @inbounds all_entries[i-1] = i
    end
    @assert length(samples) == P
    for i in 1:n
        if i != 1
            all_entries[i-1] -= 1
        end
        sample!(all_entries, samples, replace=false)
        for j in samples
            if block == 0
                if W[j,i]
                    continue
                end
                @inbounds viol = phi[j] - phi[i] + dot(Xflat[index2d(i,d)]-Xflat[index2d(j,d)], xi[index2d(i,d)])
                if viol < -TOL
                    W[j,i] = true
                    lamb[j,i] = -1e-16
                    lamb[j,i] = 0.
                    Wlen += 1
                end
            else
                if W[j,i]
                    continue
                end
                @inbounds viol = phi[i] - phi[j] + dot(Xflat[index2d(j,d)]-Xflat[index2d(i,d)], xi[index2d(j,d)])
                if viol < -TOL
                    W[j,i] = true
                    lamb[j,i] = -1e-16
                    lamb[j,i] = 0.
                    Wlen += 1
                end
            end
        end
    end
    return Wlen
end

function augmentation_rule4(phi::Array{T,1},xi::Array{T,1},Xflat::Array{T,1},Y::Array{T,1},W::ActiveSet,Wlen::Int64,lamb::LambdaActiveSet,
        n::Int64,d::Int64,viol::Array{T,1},sort_index::Array{Int64,1},samples::Array{Int64,1}, I_tmp::Array{Int64,1}, J_tmp::Array{Int64,1};TOL=1e-3,M=4n, K=n,block=0) where T<:AbstractFloat
    i::Int64 = 0
    j::Int64 = 0
    fill!(viol,0);
    len = 0;
    total_viol = 0;
    @assert length(samples) == M
    sample!(1:n*(n-1), samples, replace=false)
    for k in samples
        (i,j) = idx2pair(k, n)
        if block == 0 && W[j,i]
            continue
        elseif block == 1 && W[i,j]
            continue
        end
        len += 1
        @inbounds viol[len] = phi[j] - phi[i] + dot(Xflat[index2d(i,d)]-Xflat[index2d(j,d)], xi[index2d(i,d)])
        I_tmp[len] = i;
        J_tmp[len] = j;
        if viol[len] < -TOL
            total_viol += 1;
        end
    end
    K1 = min(K,total_viol);
    if K1 != 0
        sortperm!(sort_index,viol;alg=PartialQuickSort(K1))
        for k in 1:K1
            i = I_tmp[sort_index[k]]
            j = J_tmp[sort_index[k]]
            Wlen += 1
            if block == 0
                W[j,i] = true
                lamb[j,i] = -1e-16
                lamb[j,i] = 0.0
            else
                W[i,j] = true
                lamb[i,j] = -1e-16
                lamb[i,j] = 0.0
            end
        end
    end
    return Wlen
end

function augmentation_rule5(phi::Array{T,1},xi::Array{T,1},Xflat::Array{T,1},Y::Array{T,1},W::ActiveSet,Wlen::Int64,lamb::LambdaActiveSet,n::Int64,d::Int64,
        viol::Array{T,1},sort_index::Array{Int64,1},blocks::Array{Int64,1}, all_entries::Array{Int64,1}, entries::Array{Int64,1};
        TOL=1e-3::T,G=n÷4::Int64,P=4::Int64,block=0::Int64) where {T <:  AbstractFloat}
    fill!(viol,0)
    fill!(entries,0)
    @assert length(blocks) == G
    sample!(1:n, blocks, replace=false, ordered=true)
    @simd for i in 2:n
        @inbounds all_entries[i-1]=i
    end
    P1::Int64 = 0;
    len::Int64 = 0;
    p::Int64=0;
    last_i::Int64 = 1;
    for i in blocks
        len = 0
        all_entries[last_i:i-1] .-= 1
        last_i = i;
        for j in all_entries
            if ~W[j,i]
                len += 1
                entries[len] = j
            end
        end
        A_dot_phi_plus_B_dot_xi!(viol,phi,xi,Xflat,i,entries,len,n,d,block=block);
        P1 = min(sum(viol.<-TOL),P);
        if P1 == 0
            continue
        elseif P1 == 1
            _,p = findmin(viol);
            W[entries[p],i] = true
            lamb[entries[p],i] = -1e-16
            lamb[entries[p],i] = 0.0
            Wlen += 1
        else
            sortperm!(sort_index,viol;alg=PartialQuickSort(P1))
            for k in 1:P1
                W[entries[sort_index[k]],i] = true
                lamb[entries[sort_index[k]],i] = -1e-16
                lamb[entries[sort_index[k]],i] = 0.0
                Wlen += 1
            end
        end
    end
    return Wlen
end

struct ActiveSetAugmentation
    rule::Int;
    n::Int;
    P::Int;
    M::Int;
    K::Int;
    G::Int;
    block::Int;
    TOL::Float64;
    viol::Union{Nothing,Array{Float64,1}};
    sort_index::Union{Nothing, Array{Int64,1}};
    all_entries::Union{Nothing, Array{Int64,1}};
    entries::Union{Nothing, Array{Int64,1}};
    samples::Union{Nothing, Array{Int64,1}};
    I_tmp::Union{Nothing, Array{Int64,1}};
    J_tmp::Union{Nothing, Array{Int64,1}};
    function ActiveSetAugmentation(rule, n, block, TOL; P=0, M=0,K=0,G=0)
        if rule == 1
            P = (P == 0) ? 1 : P;
            viol = zeros(Float64, n-1)
            sort_index = zeros(Int64, n-1)
            all_entries = zeros(Int64, n-1)
            entries = zeros(Int64, n-1)
            samples = nothing;
            I_tmp = nothing;
            J_tmp = nothing;
        elseif rule == 2
            K = (K == 0) ? n : K;
            viol = nothing
            sort_index = nothing;
            all_entries = nothing;
            entries = nothing;
            samples = zeros(Int64,K);
            I_tmp = nothing;
            J_tmp = nothing;
        elseif rule == 3
            P = (P == 0) ? 1 : P;
            viol = nothing
            sort_index = nothing;
            all_entries = zeros(Int64, n-1);
            entries = nothing;
            samples = zeros(Int64, P);
            I_tmp = nothing;
            J_tmp = nothing;
        elseif rule == 4
            K = (K == 0) ? n : K;
            M = (M == 0) ? min(4*K, n*(n-1)) : max(M,K);
            viol = zeros(Float64,M);
            sort_index = zeros(Int64,M);
            all_entries = nothing;
            entries = nothing;
            samples = zeros(Int64,M);
            I_tmp = zeros(Int64, M);
            J_tmp = zeros(Int64, M);
        elseif rule == 5
            P = (P == 0) ? 4 : P;
            G = (G == 0) ? n÷P : min(n÷P, G);
            viol = zeros(Float64, n-1)
            sort_index = zeros(Int64, n-1)
            all_entries = zeros(Int64, n-1)
            entries = zeros(Int64, n-1)
            samples = zeros(Int64, G);
            I_tmp = nothing;
            J_tmp = nothing;
        end
        new(rule, n, P, M, K, G, block,TOL,viol,sort_index,all_entries,entries, samples,I_tmp,J_tmp);
    end
end

function (aug::ActiveSetAugmentation)(phi::Array{T,1},xi::Array{T,1},Xflat::Array{T,1},Y::Array{T,1},W::ActiveSet,Wlen::Int64,lamb::LambdaActiveSet,n::Int64,d::Int64) where{T<:AbstractFloat}
    if aug.rule == 1
        return augmentation_rule1(phi,xi,Xflat,Y,W,Wlen,lamb,n,d,aug.viol,aug.sort_index,aug.all_entries, aug.entries;TOL=aug.TOL,P=aug.P,block=aug.block)
    elseif aug.rule == 2
        return augmentation_rule2(phi,xi,Xflat,Y,W,Wlen,lamb,n,d,aug.samples;TOL=aug.TOL,K=aug.K,block=aug.block)
    elseif aug.rule == 3
        return augmentation_rule3(phi,xi,Xflat,Y,W,Wlen,lamb,n,d,aug.all_entries, aug.samples;TOL=aug.TOL,P=aug.P,block=aug.block)
    elseif aug.rule == 4
        return augmentation_rule4(phi,xi,Xflat,Y,W,Wlen,lamb,n,d,aug.viol,aug.sort_index,aug.samples,aug.I_tmp, aug.J_tmp;TOL=aug.TOL,M=aug.M,K=aug.K,block=aug.block)
    elseif aug.rule == 5
        return augmentation_rule5(phi,xi,Xflat,Y,W,Wlen,lamb,n,d,aug.viol,aug.sort_index,aug.samples,aug.all_entries,aug.entries;TOL=aug.TOL,G=aug.G,P=aug.P,block=aug.block)
    end
end

function initial_active_set(Xmat,Y)
    n = size(Xmat,1);
    W = ifelse.(Y[1:n-1].>Y[2:n],[(i,i+1) for i in 1:n-1],[(i+1,i) for i in 1:n-1])
    return W
end


text:
push!(LOAD_PATH,"../src/")

using Documenter
using DocumenterTools
using TensorNetworks

makedocs(
    format=     Documenter.HTML(
                    prettyurls= !("local" in ARGS),
                    canonical=  "https://quantum-many-body.github.io/TensorNetworks.jl/latest/",
                    assets=     ["assets/favicon.ico"],
                    analytics=  "UA-89508993-1",
                    ),
    clean=      false,
    sitename=   "TensorNetworks.jl",
    pages=      [
                "Home"      =>  "index.md",
                ]
)

deploydocs(
    repo=       "github.com/Quantum-Many-Body/TensorNetworks.jl.git",
    target=     "build",
    deps=       nothing,
    make=       nothing,
)

max_stars_repo_name: Quantum-Many-Body/TensorNetworks

max_issues_repo_name: Quantum-Many-Body/TensorNetworks

max_forks_repo_name: Quantum-Many-Body/TensorNetworks


text:
import Base.Threads: @spawn
using BenchmarkTools
using TimerOutputs

# THESE ARE ONLY VALID INSIDE THE BENCHMARK
query_size() = parse(Int, ENV["QUERY_SIZE"])
num_queries() = parse(Int, ENV["NUM_QUERIES"])

function main(f, n, N, thread_counts, query_latencies_secs)
    vs = fill(0, N)
    thread_counts[:] = fill(0, Threads.nthreads())[:]
    query_latencies_secs[:] = fill(0.0, N)[:]

    @sync begin
        for i in 1:N
            @spawn begin
                # Record value in vs (so work isn't compiled away)
                # Record latency for this query
                vs[i], query_latencies_secs[i] = @timed f(i, 2, n)
                # Count how many queries are handled on each thread
                thread_counts[Threads.threadid()] += 1
            end
        end
    end
    sum(vs)
end

function measure_work(work)
    n = query_size()
    N = num_queries()

    # Output params (since `@benchmark` hides the return values)
    thread_counts = fill(0, Threads.nthreads())
    query_latencies_secs = fill(0.0, N)

    b = @benchmark main($work, $n, $N, $thread_counts, $query_latencies_secs);

    # RETURN RESULTS TO test_harness.jl BY PRINTING THEM
    println(minimum(b).time)
    println(minimum(b).allocs)
    println(minimum(b).memory)
    println(minimum(b).gctime)
    println(repr(thread_counts))
    println(repr(query_latencies_secs))
end


text:
# Julia wrapper for header: darknet-master/include/darknet.h
# Automatically generated using Clang.jl wrap_c

function load_network(cfg, weights, clear)
    ccall((:load_network, libdarknet), Ptr{network}, (Cstring, Cstring, Cint), cfg, weights, clear)
end

function load_network_custom(cfg, weights, clear, batch)
    ccall((:load_network_custom, libdarknet), Ptr{network}, (Cstring, Cstring, Cint, Cint), cfg, weights, clear, batch)
end

# function load_network(cfg, weights, clear)
#     ccall((:load_network, libdarknet), Ptr{network}, (Cstring, Cstring, Cint), cfg, weights, clear)
# end

function get_base_args(net)
    ccall((:get_base_args, libdarknet), load_args, (Ptr{network},), net)
end

function do_nms_sort(dets, total, classes, thresh)
    ccall((:do_nms_sort, libdarknet), Cvoid, (Ptr{detection}, Cint, Cint, Cfloat), dets, total, classes, thresh)
end

function do_nms_obj(dets, total, classes, thresh)
    ccall((:do_nms_obj, libdarknet), Cvoid, (Ptr{detection}, Cint, Cint, Cfloat), dets, total, classes, thresh)
end

function network_predict(net, input)
    ccall((:network_predict, libdarknet), Ptr{Cfloat}, (network, Ptr{Cfloat}), net, input)
end

function network_predict_ptr(net, input)
    ccall((:network_predict_ptr, libdarknet), Ptr{Cfloat}, (Ptr{network}, Ptr{Cfloat}), net, input)
end

function get_network_boxes(net, w, h, thresh, hier, map, relative, num, letter)
    ccall((:get_network_boxes, libdarknet), Ptr{detection}, (Ptr{network}, Cint, Cint, Cfloat, Cfloat, Ptr{Cint}, Cint, Ptr{Cint}, Cint), net, w, h, thresh, hier, map, relative, num, letter)
end

function free_detections(dets, n)
    ccall((:free_detections, libdarknet), Cvoid, (Ptr{detection}, Cint), dets, n)
end

function fuse_conv_batchnorm(net)
    ccall((:fuse_conv_batchnorm, libdarknet), Cvoid, (network,), net)
end

function calculate_binary_weights(net)
    ccall((:calculate_binary_weights, libdarknet), Cvoid, (network,), net)
end

function detection_to_json(dets, nboxes, classes, names, frame_id, filename)
    ccall((:detection_to_json, libdarknet), Cstring, (Ptr{detection}, Cint, Cint, Ptr{Cstring}, Clonglong, Cstring), dets, nboxes, classes, names, frame_id, filename)
end

function get_network_layer(net, i)
    ccall((:get_network_layer, libdarknet), Ptr{layer}, (Ptr{network}, Cint), net, i)
end

function make_network_boxes(net, thresh, num)
    ccall((:make_network_boxes, libdarknet), Ptr{detection}, (Ptr{network}, Cfloat, Ptr{Cint}), net, thresh, num)
end

function reset_rnn(net)
    ccall((:reset_rnn, libdarknet), Cvoid, (Ptr{network},), net)
end

function network_predict_image(net, im)
    ccall((:network_predict_image, libdarknet), Ptr{Cfloat}, (Ptr{network}, image), net, im)
end

function validate_detector_map(datacfg, cfgfile, weightfile, thresh_calc_avg_iou, iou_thresh, existing_net)
    ccall((:validate_detector_map, libdarknet), Cfloat, (Cstring, Cstring, Cstring, Cfloat, Cfloat, Ptr{network}), datacfg, cfgfile, weightfile, thresh_calc_avg_iou, iou_thresh, existing_net)
end

function train_detector(datacfg, cfgfile, weightfile, gpus, ngpus, clear, dont_show, calc_map, mjpeg_port)
    ccall((:train_detector, libdarknet), Cvoid, (Cstring, Cstring, Cstring, Ptr{Cint}, Cint, Cint, Cint, Cint, Cint), datacfg, cfgfile, weightfile, gpus, ngpus, clear, dont_show, calc_map, mjpeg_port)
end

function test_detector(datacfg, cfgfile, weightfile, filename, thresh, hier_thresh, dont_show, ext_output, save_labels, outfile)
    ccall((:test_detector, libdarknet), Cvoid, (Cstring, Cstring, Cstring, Cstring, Cfloat, Cfloat, Cint, Cint, Cint, Cstring), datacfg, cfgfile, weightfile, filename, thresh, hier_thresh, dont_show, ext_output, save_labels, outfile)
end

function network_width(net)
    ccall((:network_width, libdarknet), Cint, (Ptr{network},), net)
end

function network_height(net)
    ccall((:network_height, libdarknet), Cint, (Ptr{network},), net)
end

function optimize_picture(net, orig, max_layer, scale, rate, thresh, norm)
    ccall((:optimize_picture, libdarknet), Cvoid, (Ptr{network}, image, Cint, Cfloat, Cfloat, Cfloat, Cint), net, orig, max_layer, scale, rate, thresh, norm)
end

function resize_image(im, w, h)
    ccall((:resize_image, libdarknet), image, (image, Cint, Cint), im, w, h)
end

function copy_image_from_bytes(im, pdata)
    ccall((:copy_image_from_bytes, libdarknet), Cvoid, (image, Cstring), im, pdata)
end

function letterbox_image(im, w, h)
    ccall((:letterbox_image, libdarknet), image, (image, Cint, Cint), im, w, h)
end

function rgbgr_image(im)
    ccall((:rgbgr_image, libdarknet), Cvoid, (image,), im)
end

function make_image(w, h, c)
    ccall((:make_image, libdarknet), image, (Cint, Cint, Cint), w, h, c)
end

function load_image_color(filename, w, h)
    ccall((:load_image_color, libdarknet), image, (Cstring, Cint, Cint), filename, w, h)
end

function free_image(m)
    ccall((:free_image, libdarknet), Cvoid, (image,), m)
end

function free_layer(arg1)
    ccall((:free_layer, libdarknet), Cvoid, (layer,), arg1)
end

function free_data(d)
    ccall((:free_data, libdarknet), Cvoid, (data,), d)
end

function load_data()
    ccall((:load_data, libdarknet), Cint, ())
end

function load_data_in_thread()
    ccall((:load_data_in_thread, libdarknet), Cint, ())
end

function cuda_pull_array(x_gpu, x, n)
    ccall((:cuda_pull_array, libdarknet), Cvoid, (Ptr{Cfloat}, Ptr{Cfloat}, Cint), x_gpu, x, n)
end

function cuda_pull_array_async(x_gpu, x, n)
    ccall((:cuda_pull_array_async, libdarknet), Cvoid, (Ptr{Cfloat}, Ptr{Cfloat}, Cint), x_gpu, x, n)
end

function cuda_set_device(n)
    ccall((:cuda_set_device, libdarknet), Cvoid, (Cint,), n)
end

function free_ptrs(ptrs, n)
    ccall((:free_ptrs, libdarknet), Cvoid, (Ptr{Ptr{Cvoid}}, Cint), ptrs, n)
end

function top_k(a, n, k, index)
    ccall((:top_k, libdarknet), Cvoid, (Ptr{Cfloat}, Cint, Cint, Ptr{Cint}), a, n, k, index)
end

function read_tree(filename)
    ccall((:read_tree, libdarknet), Ptr{tree}, (Cstring,), filename)
end

function get_metadata(file)
    ccall((:get_metadata, libdarknet), metadata, (Cstring,), file)
end

function get_time_point()
    ccall((:get_time_point, libdarknet), Cdouble, ())
end

function start_timer()
    ccall((:start_timer, libdarknet), Cvoid, ())
end

function stop_timer()
    ccall((:stop_timer, libdarknet), Cvoid, ())
end

function get_time()
    ccall((:get_time, libdarknet), Cdouble, ())
end

function stop_timer_and_show()
    ccall((:stop_timer_and_show, libdarknet), Cvoid, ())
end

function stop_timer_and_show_name(name)
    ccall((:stop_timer_and_show_name, libdarknet), Cvoid, (Cstring,), name)
end

function show_total_time()
    ccall((:show_total_time, libdarknet), Cvoid, ())
end


text:
module DiffEqBase

using RecipesBase, RecursiveArrayTools,
      Requires, Tables, TableTraits, IteratorInterfaceExtensions, TreeViews,
      IterativeSolvers, RecursiveFactorization, Distributed, ArrayInterface

import Logging

import ZygoteRules, ChainRulesCore, LabelledArrays

using Roots # callbacks

using StaticArrays # data arrays

using LinearAlgebra, Statistics, Printf

using DocStringExtensions

using FunctionWrappers: FunctionWrapper

using MuladdMacro, Parameters

# Problems
"""
$(TYPEDEF)

Base type for all DifferentialEquations.jl problems. Concrete subtypes of
`DEProblem` contain the necessary information to fully define a differential
equation of the corresponding type.
"""
abstract type DEProblem end

"""
$(TYPEDEF)
"""
abstract type DEElement end

"""
$(TYPEDEF)
"""
abstract type DESensitivity end

"""
$(TYPEDEF)

Base for types which define linear systems.
"""
abstract type AbstractLinearProblem{bType,isinplace} <: DEProblem end

"""
$(TYPEDEF)

Base for types which define nonlinear systems.
"""
abstract type AbstractNonlinearProblem{uType,isinplace} <: DEProblem end

"""
$(TYPEDEF)

Base for types which define integrals suitable for quadrature.
"""
abstract type AbstractQuadratureProblem{isinplace} <: DEProblem end

"""
$(TYPEDEF)

Base for types which define equations for optimization
"""
abstract type AbstractOptimizationProblem{isinplace} <: DEProblem end

"""
$(TYPEDEF)

Base for types which define steady state problems for ODE systems.
"""
abstract type AbstractSteadyStateProblem{uType,isinplace} <: DEProblem end

"""
$(TYPEDEF)
"""
abstract type AbstractNoiseProblem <: DEProblem end

"""
$(TYPEDEF)

Base for types which define ODE problems.
"""
abstract type AbstractODEProblem{uType,tType,isinplace} <: DEProblem end

"""
$(TYPEDEF)

Base for types which define discrete problems.
"""
abstract type AbstractDiscreteProblem{uType,tType,isinplace} <:
                      AbstractODEProblem{uType,tType,isinplace} end

"""
$(TYPEDEF)
"""
abstract type AbstractAnalyticalProblem{uType,tType,isinplace} <:
                      AbstractODEProblem{uType,tType,isinplace} end

"""
$(TYPEDEF)

Base for types which define RODE problems.
"""
abstract type AbstractRODEProblem{uType,tType,isinplace,ND} <: DEProblem end

"""
$(TYPEDEF)

Base for types which define SDE problems.
"""
abstract type AbstractSDEProblem{uType,tType,isinplace,ND} <:
                      AbstractRODEProblem{uType,tType,isinplace,ND} end

"""
$(TYPEDEF)

Base for types which define DAE problems.
"""
abstract type AbstractDAEProblem{uType,duType,tType,isinplace} <: DEProblem end

"""
$(TYPEDEF)

Base for types which define DDE problems.
"""
abstract type AbstractDDEProblem{uType,tType,lType,isinplace} <: DEProblem end

"""
$(TYPEDEF)
"""
abstract type AbstractConstantLagDDEProblem{uType,tType,lType,isinplace} <:
                      AbstractDDEProblem{uType,tType,lType,isinplace} end

"""
$(TYPEDEF)
"""
abstract type AbstractSecondOrderODEProblem{uType,tType,isinplace} <: AbstractODEProblem{uType,tType,isinplace} end

"""
$(TYPEDEF)

Base for types which define BVP problems.
"""
abstract type AbstractBVProblem{uType,tType,isinplace} <: AbstractODEProblem{uType,tType,isinplace} end

"""
$(TYPEDEF)

Base for types which define jump problems.
"""
abstract type AbstractJumpProblem{P,J} <: DiffEqBase.DEProblem end

"""
$(TYPEDEF)

Base for types which define SDDE problems.
"""
abstract type AbstractSDDEProblem{uType,tType,lType,isinplace,ND} <: DEProblem end

"""
$(TYPEDEF)
"""
abstract type AbstractConstantLagSDDEProblem{uType,tType,lType,isinplace,ND} <:
                      AbstractSDDEProblem{uType,tType,lType,isinplace,ND} end

"""
$(TYPEDEF)

Base for types which define PDE problems.
"""
abstract type AbstractPDEProblem <: DEProblem end

# Algorithms
"""
$(TYPEDEF)
"""
abstract type DEAlgorithm end

"""
$(TYPEDEF)
"""
abstract type AbstractLinearAlgorithm <: DEAlgorithm end

"""
$(TYPEDEF)
"""
abstract type AbstractNonlinearAlgorithm <: DEAlgorithm end

"""
$(TYPEDEF)
"""
abstract type AbstractQuadratureAlgorithm <: DEAlgorithm end

"""
$(TYPEDEF)
"""
abstract type AbstractSteadyStateAlgorithm <: DEAlgorithm end

"""
$(TYPEDEF)
"""
abstract type AbstractODEAlgorithm <: DEAlgorithm end

"""
$(TYPEDEF)
"""
abstract type AbstractSecondOrderODEAlgorithm <: DEAlgorithm end

"""
$(TYPEDEF)
"""
abstract type AbstractRODEAlgorithm <: DEAlgorithm end

"""
$(TYPEDEF)
"""
abstract type AbstractSDEAlgorithm <: DEAlgorithm end

"""
$(TYPEDEF)
"""
abstract type AbstractDAEAlgorithm <: DEAlgorithm end

"""
$(TYPEDEF)
"""
abstract type AbstractDDEAlgorithm <: DEAlgorithm end

"""
$(TYPEDEF)
"""
abstract type AbstractSDDEAlgorithm <: DEAlgorithm end
"""
$(TYPEDEF)
"""
abstract type EnsembleAlgorithm <: DiffEqBase.DEAlgorithm end

"""
$(TYPEDEF)
"""
abstract type AbstractSensitivityAlgorithm{CS,AD,FDT} <: DiffEqBase.DEAlgorithm end

"""
$(TYPEDEF)
"""
abstract type DAEInitializationAlgorithm end

# PDE Discretizations

"""
$(TYPEDEF)
"""
abstract type AbstractDiscretization end

# Monte Carlo Simulations
"""
$(TYPEDEF)
"""
abstract type AbstractEnsembleProblem <: DEProblem end

"""
$(TYPEDEF)
"""
abstract type AbstractEnsembleEstimator <: DEProblem end

export EnsembleProblem
export EnsembleSolution, EnsembleTestSolution, EnsembleSummary

"""
$(TYPEDEF)
"""
abstract type AbstractDiffEqInterpolation <: Function end

"""
$(TYPEDEF)
"""
abstract type AbstractDEOptions end

"""
$(TYPEDEF)
"""
abstract type DECache end

"""
$(TYPEDEF)
"""
abstract type DECallback end

"""
$(TYPEDEF)
"""
abstract type AbstractContinuousCallback <: DECallback end

"""
$(TYPEDEF)
"""
abstract type AbstractDiscreteCallback <: DECallback end

"""
$(TYPEDEF)
"""
abstract type DEDataArray{T,N} <: AbstractArray{T,N} end
const DEDataVector{T} = DEDataArray{T,1}
const DEDataMatrix{T} = DEDataArray{T,2}

# Integrators
"""
$(TYPEDEF)
"""
abstract type DEIntegrator{Alg, IIP, U, T} end

"""
$(TYPEDEF)
"""
abstract type AbstractSteadyStateIntegrator{Alg, IIP, U} <: DEIntegrator{Alg, IIP, U, Nothing} end

"""
$(TYPEDEF)
"""
abstract type AbstractODEIntegrator{Alg, IIP, U, T} <: DEIntegrator{Alg, IIP, U, T} end

"""
$(TYPEDEF)
"""
abstract type AbstractSecondOrderODEIntegrator{Alg, IIP, U, T} <: DEIntegrator{Alg, IIP, U, T} end

"""
$(TYPEDEF)
"""
abstract type AbstractRODEIntegrator{Alg, IIP, U, T} <: DEIntegrator{Alg, IIP, U, T} end

"""
$(TYPEDEF)
"""
abstract type AbstractSDEIntegrator{Alg, IIP, U, T} <: DEIntegrator{Alg, IIP, U, T} end

"""
$(TYPEDEF)
"""
abstract type AbstractDDEIntegrator{Alg, IIP, U, T} <: DEIntegrator{Alg, IIP, U, T} end

"""
$(TYPEDEF)
"""
abstract type AbstractDAEIntegrator{Alg, IIP, U, T} <: DEIntegrator{Alg, IIP, U, T} end

"""
$(TYPEDEF)
"""
abstract type AbstractSDDEIntegrator{Alg, IIP, U, T} <: DEIntegrator{Alg, IIP, U, T} end

# Solutions
"""
$(TYPEDEF)
"""
abstract type AbstractNoTimeSolution{T,N} <: AbstractArray{T,N} end

"""
$(TYPEDEF)
"""
abstract type AbstractTimeseriesSolution{T,N,A} <: AbstractDiffEqArray{T,N,A} end

"""
$(TYPEDEF)
"""
abstract type AbstractEnsembleSolution{T,N,A} <: AbstractVectorOfArray{T,N,A} end

"""
$(TYPEDEF)
"""
abstract type AbstractNoiseProcess{T,N,A,isinplace} <: AbstractDiffEqArray{T,N,A} end

"""
Union of all base solution types.
"""
const DESolution = Union{AbstractTimeseriesSolution,
                         AbstractNoTimeSolution,
                         AbstractEnsembleSolution,
                         AbstractNoiseProcess}
export DESolution

"""
$(TYPEDEF)
"""
abstract type AbstractLinearSolution{T,N} <: AbstractNoTimeSolution{T,N} end

"""
$(TYPEDEF)
"""
abstract type AbstractNonlinearSolution{T,N} <: AbstractNoTimeSolution{T,N} end

"""
$(TYPEDEF)
"""
abstract type AbstractQuadratureSolution{T,N} <: AbstractNoTimeSolution{T,N} end

"""
$(TYPEDEF)
"""
abstract type AbstractSteadyStateSolution{T,N} <: AbstractNoTimeSolution{T,N} end

"""
$(TYPEDEF)
"""
abstract type AbstractAnalyticalSolution{T,N,S} <: AbstractTimeseriesSolution{T,N,S} end

"""
$(TYPEDEF)
"""
abstract type AbstractODESolution{T,N,S} <: AbstractTimeseriesSolution{T,N,S} end

# Needed for plot recipes
"""
$(TYPEDEF)
"""
abstract type AbstractDDESolution{T,N,S} <: AbstractODESolution{T,N,S} end

"""
$(TYPEDEF)
"""
abstract type AbstractRODESolution{T,N,S} <: AbstractODESolution{T,N,S} end

"""
$(TYPEDEF)
"""
abstract type AbstractDAESolution{T,N,S} <: AbstractODESolution{T,N,S} end

"""
$(TYPEDEF)
"""
abstract type AbstractSensitivitySolution{T,N,S} <: AbstractTimeseriesSolution{T,N,S} end

# Misc
"""
$(TYPEDEF)
"""
abstract type Tableau end

"""
$(TYPEDEF)
"""
abstract type ODERKTableau <: Tableau end

"""
$(TYPEDEF)
"""
abstract type DECostFunction end

"""
$(TYPEDEF)
"""
abstract type AbstractDiffEqOperator{T} end

"""
$(TYPEDEF)
"""
abstract type AbstractDiffEqLinearOperator{T} <: AbstractDiffEqOperator{T} end

"""
$(TYPEDEF)

Base for types defining differential equation functions.
"""
abstract type AbstractDiffEqFunction{iip} <: Function end

"""
$(TYPEDEF)

Base for types which define the history of a delay differential equation.
"""
abstract type AbstractHistoryFunction <: Function end

"""
$(TYPEDEF)
"""
abstract type AbstractReactionNetwork <: Function end

include("fastpow.jl")
include("diffeqfastbc.jl")
include("destats.jl")
include("utils.jl")
include("calculate_residuals.jl")
include("solutions/basic_solutions.jl")
include("solutions/steady_state_solutions.jl")
include("solutions/ode_solutions.jl")
include("solutions/rode_solutions.jl")
include("solutions/dae_solutions.jl")
include("solutions/solution_interface.jl")
include("tableaus.jl")
include("function_wrappers.jl")
include("diffeqfunction.jl")
include("problems/problem_utils.jl")
include("problems/discrete_problems.jl")
include("problems/steady_state_problems.jl")
include("problems/analytical_problems.jl")
include("problems/basic_problems.jl")
include("problems/ode_problems.jl")
include("problems/rode_problems.jl")
include("problems/sde_problems.jl")
include("problems/noise_problems.jl")
include("problems/bvp_problems.jl")
include("problems/dae_problems.jl")
include("problems/dde_problems.jl")
include("problems/sdde_problems.jl")
include("problems/pde_problems.jl")
include("problems/problem_traits.jl")
include("ensemble/ensemble_solutions.jl")
include("ensemble/ensemble_problems.jl")
include("ensemble/basic_ensemble_solve.jl")
include("ensemble/ensemble_analysis.jl")
include("nlsolve/type.jl")
include("nlsolve/newton.jl")
include("nlsolve/functional.jl")
include("nlsolve/utils.jl")
include("operators/diffeq_operator.jl")
include("operators/common_defaults.jl")
include("operators/basic_operators.jl")
include("interpolation.jl")
include("callbacks.jl")
include("integrator_interface.jl")
include("linear_nonlinear.jl")
include("common_defaults.jl")
include("data_array.jl")
include("solve.jl")
include("internal_euler.jl")
include("tabletraits.jl")
include("alg_traits.jl")
include("remake.jl")
include("init.jl")
include("zygote.jl")

"""
$(TYPEDEF)
"""
abstract type AbstractParameterizedFunction{iip} <: AbstractODEFunction{iip} end

"""
$(TYPEDEF)
"""
struct ConvergenceSetup{P,C}
    probs::P
    convergence_axis::C
end

const AbstractMonteCarloProblem = AbstractEnsembleProblem
const AbstractMonteCarloSolution = AbstractEnsembleSolution
const MonteCarloAlgorithm = EnsembleAlgorithm
const MonteCarloProblem = EnsembleProblem
const MonteCarloSolution = EnsembleSolution
const MonteCarloSummary = EnsembleSummary
@deprecate MonteCarloProblem(args...) EnsembleProblem(args...)
@deprecate MonteCarloSolution(args...) EnsembleSolution(args...)
@deprecate MonteCarloSummary(args...) EnsembleSummary(args...)
@deprecate calculate_monte_errors(args...;kwargs...) calculate_ensemble_errors(args...;kwargs...)

export isinplace

export solve, solve!, init, step!, discretize, concrete_solve

export tuples, intervals, TimeChoiceIterator

export resize!,deleteat!,addat!,get_tmp_cache,
       full_cache,user_cache,u_cache,du_cache,
       rand_cache,ratenoise_cache,
       resize_non_user_cache!,deleteat_non_user_cache!,addat_non_user_cache!,
       terminate!,
       add_tstop!,add_saveat!,set_abstol!,
       set_reltol!,get_du, get_du!, get_dt,get_proposed_dt,set_proposed_dt!,
       u_modified!, savevalues!,reinit!, auto_dt_reset!, set_t!,
       set_u!, check_error, change_t_via_interpolation!, addsteps!,
       isdiscrete, reeval_internals_due_to_modification!

export LinearProblem, NonlinearProblem, QuadratureProblem, OptimizationProblem

export DiscreteProblem
export SteadyStateProblem, SteadyStateSolution
export NoiseProblem
export ODEProblem, ODESolution
export DynamicalODEFunction, DynamicalODEProblem,
       SecondOrderODEProblem, SplitFunction, SplitODEProblem
export SplitSDEProblem
export RODEProblem, RODESolution, SDEProblem
export DAEProblem, DAESolution
export DDEProblem
export SDDEProblem
export PDEProblem
export IncrementingODEProblem

export BVProblem, TwoPointBVProblem

export remake

export DEDataArray, DEDataVector, DEDataMatrix

export ODEFunction, DiscreteFunction, SplitFunction, DAEFunction, DDEFunction,
       SDEFunction, SplitSDEFunction, RODEFunction, SDDEFunction, IncrementingODEFunction

export ContinuousCallback, DiscreteCallback, CallbackSet, VectorContinuousCallback

export initialize!

export LinSolveFactorize, LinSolveGPUFactorize, DefaultLinSolve, DEFAULT_LINSOLVE,
       LinSolveGMRES, LinSolveCG, LinSolveBiCGStabl, LinSolveChebyshev,
       LinSolveMINRES, LinSolveIterativeSolvers

export AffineDiffEqOperator, update_coefficients!, update_coefficients,
       has_expmv!, has_expmv, has_exp, has_mul, has_mul!, has_ldiv, has_ldiv!

export DiffEqScalar, DiffEqArrayOperator, DiffEqIdentity

export NLNewton, NLFunctional, NLAnderson

export EnsembleThreads, EnsembleDistributed, EnsembleSplitThreads, EnsembleSerial

export EnsembleAnalysis

export SensitivityADPassThrough

end # module


text:
## Aircraft analysis benchmarking
using BenchmarkTools
using TimerOutputs

# All subsequent analyses use no symmetry tricks for performance as AeroMDAO hasn't implemented them and apples must be compared to apples.

## AeroMDAO tests: https://github.com/GodotMisogi/AeroMDAO.jl
#=======================================================#
using AeroMDAO

println("AeroMDAO Aircraft -")
function aeromdao_steady_vlm()
    # Wing
    wing = Wing(foils     = fill(naca4((0,0,1,2)), 2),
                chords    = [1.0, 0.6],
                twists    = [2.0, 2.0],
                spans     = [10.0],
                dihedrals = [11.31],
                sweeps    = [2.29]);

    # Horizontal tail
    htail = Wing(foils     = fill(naca4((0,0,1,2)), 2),
                 chords    = [0.7, 0.42],
                 twists    = [0.0, 0.0],
                 spans     = [2.5],
                 dihedrals = [0.],
                 sweeps    = [6.39],
                 position  = [4., 0, 0],
                 angle     = 0.,
                 axis      = [0., 1., 0.])

    # Vertical tail
    vtail = HalfWing(foils     = fill(naca4((0,0,0,9)), 2),
                     chords    = [0.7, 0.42],
                     twists    = [0.0, 0.0],
                     spans     = [1.0],
                     dihedrals = [0.],
                     sweeps    = [7.97],
                     position  = [4., 0, 0],
                     angle     = 90.,
                     axis      = [1., 0., 0.])

    wing_panels, wing_normals   = panel_wing(wing, 20, 10,
                                             spacing = Cosine()
                                            )

    htail_panels, htail_normals = panel_wing(htail, 12, 6;
                                             spacing  = Cosine()
                                            )

    vtail_panels, vtail_normals = panel_wing(vtail, 10, 5;
                                             spacing  = Cosine()
                                            )

    # Aircraft assembly
    aircraft = [
                Horseshoe.(wing_panels,  wing_normals)[:];
                Horseshoe.(htail_panels, htail_normals)[:];
                Horseshoe.(vtail_panels, vtail_normals)[:];
               ]

    α, β = 5.0, 0.0
    Ω    = [0.0, 0.0, 0.0]
    fs   = AeroMDAO.Freestream(α, β, Ω)
    refs = References(
                         speed     = 10.0,
                         density   = 1.225,
                         viscosity = 1.5e-5,
                         area      = 9.0,
                         span      = 10.0,
                         chord     = 0.9,
                         location  = [0.5, 0., 0.]
                     )

    # reset_timer!()
    # @timeit "Solving Case"
    data = solve_case(aircraft, fs, refs;)

    # @timeit "Computing Forces" Fs = surface_forces(data)
    CFs, CMs = surface_coefficients(data)
    # FFs      = farfield_coefficients(data)

    sum(CFs), sum(CMs)
    # nearfield(data), farfield(data)
    # print_timer()
end

t1 = @benchmark aeromdao_steady_vlm()

##
println("AeroMDAO Aircraft ComponentArrays.jl -")
function aeromdao_steady_vlm_components()
    # Wing
    wing = Wing(foils     = fill(naca4((0,0,1,2)), 2),
                chords    = [1.0, 0.6],
                twists    = [2.0, 2.0],
                spans     = [10.0],
                dihedrals = [11.31],
                sweeps    = [2.29]);

    # Horizontal tail
    htail = Wing(foils     = fill(naca4((0,0,1,2)), 2),
                 chords    = [0.7, 0.42],
                 twists    = [0.0, 0.0],
                 spans     = [2.5],
                 dihedrals = [0.],
                 sweeps    = [6.39],
                 position  = [4., 0, 0],
                 angle     = 0.,
                 axis      = [0., 1., 0.])

    # Vertical tail
    vtail = HalfWing(foils     = fill(naca4((0,0,0,9)), 2),
                     chords    = [0.7, 0.42],
                     twists    = [0.0, 0.0],
                     spans     = [1.0],
                     dihedrals = [0.],
                     sweeps    = [7.97],
                     position  = [4., 0, 0],
                     angle     = 90.,
                     axis      = [1., 0., 0.])

    wing_mesh = WingMesh(wing, [20], 10,
                         span_spacing = Cosine()
                        )

    htail_mesh = WingMesh(htail, [12], 6;
                          span_spacing  = Cosine()
                         )

    vtail_mesh = WingMesh(vtail, [10], 5;
                          span_spacing  = Cosine()
                         )

    # Aircraft assembly
    aircraft = ComponentArray(
                    wing  = make_horseshoes(wing_mesh),
                    htail = make_horseshoes(htail_mesh),
                    vtail = make_horseshoes(vtail_mesh),
                )

    α, β = 5.0, 0.0
    Ω    = [0.0, 0.0, 0.0]
    fs   = AeroMDAO.Freestream(α, β, Ω)
    refs = References(
                         speed     = 10.0,
                         density   = 1.225,
                         viscosity = 1.5e-5,
                         area      = 9.0,
                         span      = 10.0,
                         chord     = 0.9,
                         location  = [0.5, 0., 0.]
                     )

    # reset_timer!()
    # @timeit "Solving Case"
    data = solve_case(aircraft, fs, refs;)

    # @timeit "Computing Forces" Fs = surface_forces(data)
    CFs, CMs = surface_coefficients(data)
    # FFs      = farfield_coefficients(data)

    sum(CFs), sum(CMs)
    # nearfield(data), farfield(data)
    # print_timer()
end

t2 = @benchmark aeromdao_steady_vlm_components()

## BYU FLOW Lab tests: https://github.com/byuflowlab/VortexLattice.jl
#=======================================================#

# Example taken from: https://flow.byu.edu/VortexLattice.jl/dev/examples/

# Pulled these out to avoid compilation for proper benchmarking
fc   = fill((xc) -> 0, 2) # camberline function for each section
fc_h = fill((xc) -> 0, 2) # camberline function for each section
fc_v = fill((xc) -> 0, 2) # camberline function for each section

##
using VortexLattice

println("BYU FLOW Lab VortexLattice.jl - ")
function byu_steady_vlm()
    # wing
    xle = [0.0, 0.2]
    yle = [0.0, 5.0]
    zle = [0.0, 1.0]
    chord = [1.0, 0.6]
    theta = [2.0*pi/180, 2.0*pi/180]
    phi = [0.0, 0.0]
    ns = 20
    nc = 10
    spacing_s = VortexLattice.Cosine()
    spacing_c = VortexLattice.Cosine()
    mirror = true

    # horizontal stabilizer
    xle_h = [0.0, 0.14]
    yle_h = [0.0, 1.25]
    zle_h = [0.0, 0.0]
    chord_h = [0.7, 0.42]
    theta_h = [0.0, 0.0]
    phi_h = [0.0, 0.0]
    ns_h = 12
    nc_h = 6
    spacing_s_h = VortexLattice.Cosine()
    spacing_c_h = VortexLattice.Cosine()
    mirror_h = true

    # vertical stabilizer
    xle_v = [0.0, 0.14]
    yle_v = [0.0, 0.0]
    zle_v = [0.0, 1.0]
    chord_v = [0.7, 0.42]
    theta_v = [0.0, 0.0]
    phi_v = [0.0, 0.0]
    ns_v = 10
    nc_v = 5
    spacing_s_v = VortexLattice.Cosine()
    spacing_c_v = VortexLattice.Cosine()
    mirror_v = false

    Sref = 9.0
    cref = 0.9
    bref = 10.0
    rref = [0.5, 0.0, 0.0]
    Vinf = 1.0
    ref = Reference(Sref, cref, bref, rref, Vinf)

    alph = 5.0*pi/180
    beta = 0.0
    Omega = [0.0; 0.0; 0.0]
    fs = VortexLattice.Freestream(Vinf, alph, beta, Omega)

    symmetric = [false, false, false]

    # generate surface panels for wing
    wgrid, wing = wing_to_surface_panels(xle, yle, zle, chord, theta, phi, ns, nc;
        mirror=mirror, fc = fc, spacing_s=spacing_s, spacing_c=spacing_c)

    # generate surface panels for horizontal tail
    hgrid, htail = wing_to_surface_panels(xle_h, yle_h, zle_h, chord_h, theta_h, phi_h, ns_h, nc_h;
        mirror=mirror_h, fc=fc_h, spacing_s=spacing_s_h, spacing_c=spacing_c_h)
    VortexLattice.translate!(hgrid, [4.0, 0.0, 0.0]) # Plots also exports translate! -_-
    VortexLattice.translate!(htail, [4.0, 0.0, 0.0]) # Plots also exports translate! -_-

    # generate surface panels for vertical tail
    vgrid, vtail = wing_to_surface_panels(xle_v, yle_v, zle_v, chord_v, theta_v, phi_v, ns_v, nc_v;
        mirror=mirror_v, fc=fc_v, spacing_s=spacing_s_v, spacing_c=spacing_c_v)
    VortexLattice.translate!(vgrid, [4.0, 0.0, 0.0]) # Plots also exports translate! -_-
    VortexLattice.translate!(vtail, [4.0, 0.0, 0.0]) # Plots also exports translate! -_-

    # display(size.([ wing, htail, vtail ])) # Checking sizes

    grids = [wgrid, hgrid, vgrid]
    surfs = [wing, htail, vtail]
    surface_id = [1, 2, 3]

    system = steady_analysis(surfs, ref, fs; symmetric=symmetric, surface_id=surface_id)

    CF, CM = body_forces(system; frame=VortexLattice.Wind())

    # CDiff = far_field_drag(system)

    # dCF_b, dCM_b = body_derivatives(system)
    # dCF_s, dCM_s = stability_derivatives(system)

    # CF, CM, CDiff #, dCF_b, dCM_b, dCF_s, dCM_s
end

t2 = @benchmark byu_steady_vlm()


text:
using JSON3

include("../src/parser.jl")

folder_for_generated_types = "../examples/acme_parsed_schemas/"

filename = "../examples/acme_schemas/acme_flatschema.json"
json_schema = JSON3.read(read(filename, String))
generate_types(json_schema;
    gen_folder=folder_for_generated_types,
    module_name = "AcmeFlat")

filename = "../examples/acme_schemas/acme_nestedschema.json"
json_schema = JSON3.read(read(filename, String))
generate_types(json_schema;
    gen_folder=folder_for_generated_types,
    module_name = "AcmeNested")

filename = "../examples/acme_schemas/acme_nestedschema_mod.json"
json_schema = JSON3.read(read(filename, String))
generate_types(json_schema;
    gen_folder=folder_for_generated_types,
    module_name = "AcmeNestedMod")


text:
#
# Time analysis
#

struct TimeSeriesAnalysis{T, D, A}
        time :: T
        data :: D
    analysis :: A
end

TimeSeriesAnalysis(time, analysis) = TimeSeriesAnalysis(time, zeros(length(time)), analysis)

struct TimeAverage end

# Use trapz integral to compute time average of data in case times are not evenly spaced
@inline (::TimeAverage)(data, time) = trapz(data, time) / (time[end] - time[1])

#
# Ensemble time analysis
#

struct EnsembleTimeSeriesAnalysis{T, D, A}
        time :: T
        data :: D
    analysis :: A
end

# Ensemble approach assumes evenly spaced time steps
EnsembleTimeSeriesAnalysis(time, ensemble_size) = EnsembleTimeSeriesAnalysis(time, zeros(ensemble_size, length(time)), data -> mean(data, dims=2))


text:
using Histograms
using Base.Test

NENTRIES = 10000
x = randn(NENTRIES)
h1 = H1D(100, -10, 10)
for i in x
    hfill!(h1, i)
end
@test sum(h1.entries) == NENTRIES

h2 = H2D(100, -10, 10, 100, -10, 10)
y = randn(NENTRIES)
for i = 1:NENTRIES
    hfill!(h2, x[i], y[i])
end
@test sum(h2.entries) == NENTRIES


text:
module AnalogMTK

#==Base simulation engine
===============================================================================#
using ModelingToolkit
import ModelingToolkit.IfElse: ifelse
const MTK_PATH = abspath(joinpath(dirname(pathof(ModelingToolkit)), ".."))
@info("Using MTK_PATH: $MTK_PATH")

#==Plotting infrastructure
===============================================================================#
using CMDimData
using CMDimData.MDDatasets
using CMDimData.EasyPlot
using CMDimData.EasyPlot.Colors #Should be ok with whatever version it needs
import CMDimData.EasyPlot: BoundingBox #Use this one to avoid version conflicts with Graphics.
using InspectDR
CMDimData.@includepkg EasyPlotInspect

#Start using electrical components defined in MTK_PATH/examples/:
include(joinpath(MTK_PATH, "examples", "electrical_components.jl"))
include("base.jl")
include("signals.jl")
include("sources.jl")
include("diodes.jl")
include("display.jl")
include("plot_base.jl")
include("plot_transient.jl")

#Signal & device types
export Ideal, Shockley
export Sine

#Circuit element definitions:
export Resistor, Capacitor, Inductor
export Diode
export VSource
export Ground

#Plotting:
export plot_transient, plot_transient!

#Displaying plots:
export inlinedisp, saveimage, displaygui

#Inform user of the existance of `t`:
@info("NOTE: AnalogMTK also defines parameter `t`")

end # module


text:
function FieldTransform(arg0::Field, arg1::Transform)
    return FieldTransform((Field, Transform), arg0, arg1)
end

function FieldTransform(arg0::FieldAbsoluteDate, arg1::FieldAngularCoordinates)
    return FieldTransform((FieldAbsoluteDate, FieldAngularCoordinates), arg0, arg1)
end

function FieldTransform(arg0::FieldAbsoluteDate, arg1::FieldPVCoordinates)
    return FieldTransform((FieldAbsoluteDate, FieldPVCoordinates), arg0, arg1)
end

function FieldTransform(arg0::FieldAbsoluteDate, arg1::FieldRotation)
    return FieldTransform((FieldAbsoluteDate, FieldRotation), arg0, arg1)
end

function FieldTransform(arg0::FieldAbsoluteDate, arg1::FieldRotation, arg2::FieldVector3D)
    return FieldTransform((FieldAbsoluteDate, FieldRotation, FieldVector3D), arg0, arg1, arg2)
end

function FieldTransform(arg0::FieldAbsoluteDate, arg1::FieldRotation, arg2::FieldVector3D, arg3::FieldVector3D)
    return FieldTransform((FieldAbsoluteDate, FieldRotation, FieldVector3D, FieldVector3D), arg0, arg1, arg2, arg3)
end

function FieldTransform(arg0::FieldAbsoluteDate, arg1::FieldTransform, arg2::FieldTransform)
    return FieldTransform((FieldAbsoluteDate, FieldTransform, FieldTransform), arg0, arg1, arg2)
end

function FieldTransform(arg0::FieldAbsoluteDate, arg1::FieldVector3D)
    return FieldTransform((FieldAbsoluteDate, FieldVector3D), arg0, arg1)
end

function FieldTransform(arg0::FieldAbsoluteDate, arg1::FieldVector3D, arg2::FieldVector3D)
    return FieldTransform((FieldAbsoluteDate, FieldVector3D, FieldVector3D), arg0, arg1, arg2)
end

function FieldTransform(arg0::FieldAbsoluteDate, arg1::FieldVector3D, arg2::FieldVector3D, arg3::FieldVector3D)
    return FieldTransform((FieldAbsoluteDate, FieldVector3D, FieldVector3D, FieldVector3D), arg0, arg1, arg2, arg3)
end

function freeze(obj::FieldTransform)
    return jcall(obj, "freeze", FieldTransform, ())
end

function get_acceleration(obj::FieldTransform)
    return jcall(obj, "getAcceleration", FieldVector3D, ())
end

function get_angular(obj::FieldTransform)
    return jcall(obj, "getAngular", FieldAngularCoordinates, ())
end

function get_cartesian(obj::FieldTransform)
    return jcall(obj, "getCartesian", FieldPVCoordinates, ())
end

function get_date(obj::FieldTransform)
    return jcall(obj, "getDate", AbsoluteDate, ())
end

function get_field_date(obj::FieldTransform)
    return jcall(obj, "getFieldDate", FieldAbsoluteDate, ())
end

function get_identity(::Type{FieldTransform}, arg0::Field)
    return jcall(FieldTransform, "getIdentity", FieldTransform, (Field,), arg0)
end

function get_inverse(obj::FieldTransform)
    return jcall(obj, "getInverse", FieldTransform, ())
end

function get_jacobian(obj::FieldTransform, arg0::CartesianDerivativesFilter, arg1::Vector{Vector{RealFieldElement}})
    return jcall(obj, "getJacobian", void, (CartesianDerivativesFilter, Vector{Vector{RealFieldElement}}), arg0, arg1)
end

function get_rotation(obj::FieldTransform)
    return jcall(obj, "getRotation", FieldRotation, ())
end

function get_rotation_acceleration(obj::FieldTransform)
    return jcall(obj, "getRotationAcceleration", FieldVector3D, ())
end

function get_rotation_rate(obj::FieldTransform)
    return jcall(obj, "getRotationRate", FieldVector3D, ())
end

function get_translation(obj::FieldTransform)
    return jcall(obj, "getTranslation", FieldVector3D, ())
end

function get_velocity(obj::FieldTransform)
    return jcall(obj, "getVelocity", FieldVector3D, ())
end

function interpolate(::Type{FieldTransform}, arg0::FieldAbsoluteDate, arg1::CartesianDerivativesFilter, arg2::AngularDerivativesFilter, arg3::Collection)
    return jcall(FieldTransform, "interpolate", FieldTransform, (FieldAbsoluteDate, CartesianDerivativesFilter, AngularDerivativesFilter, Collection), arg0, arg1, arg2, arg3)
end

function interpolate(::Type{FieldTransform}, arg0::FieldAbsoluteDate, arg1::CartesianDerivativesFilter, arg2::AngularDerivativesFilter, arg3::Stream)
    return jcall(FieldTransform, "interpolate", FieldTransform, (FieldAbsoluteDate, CartesianDerivativesFilter, AngularDerivativesFilter, Stream), arg0, arg1, arg2, arg3)
end

function interpolate(::Type{FieldTransform}, arg0::FieldAbsoluteDate, arg1::Collection)
    return jcall(FieldTransform, "interpolate", FieldTransform, (FieldAbsoluteDate, Collection), arg0, arg1)
end

function shifted_by(obj::FieldTransform, arg0::RealFieldElement)
    return jcall(obj, "shiftedBy", FieldTransform, (RealFieldElement,), arg0)
end

function shifted_by(obj::FieldTransform, arg0::jdouble)
    return jcall(obj, "shiftedBy", FieldTransform, (jdouble,), arg0)
end

function transform_line(obj::FieldTransform, arg0::FieldLine)
    return jcall(obj, "transformLine", FieldLine, (FieldLine,), arg0)
end

function transform_line(obj::FieldTransform, arg0::Line)
    return jcall(obj, "transformLine", FieldLine, (Line,), arg0)
end

function transform_position(obj::FieldTransform, arg0::FieldVector3D)
    return jcall(obj, "transformPosition", FieldVector3D, (FieldVector3D,), arg0)
end

function transform_position(obj::FieldTransform, arg0::Vector3D)
    return jcall(obj, "transformPosition", FieldVector3D, (Vector3D,), arg0)
end

function transform_pv_coordinates(obj::FieldTransform, arg0::FieldPVCoordinates)
    return jcall(obj, "transformPVCoordinates", FieldPVCoordinates, (FieldPVCoordinates,), arg0)
end

function transform_pv_coordinates(obj::FieldTransform, arg0::PVCoordinates)
    return jcall(obj, "transformPVCoordinates", FieldPVCoordinates, (PVCoordinates,), arg0)
end

function transform_pv_coordinates(obj::FieldTransform, arg0::TimeStampedFieldPVCoordinates)
    return jcall(obj, "transformPVCoordinates", TimeStampedFieldPVCoordinates, (TimeStampedFieldPVCoordinates,), arg0)
end

function transform_pv_coordinates(obj::FieldTransform, arg0::TimeStampedPVCoordinates)
    return jcall(obj, "transformPVCoordinates", TimeStampedFieldPVCoordinates, (TimeStampedPVCoordinates,), arg0)
end

function transform_vector(obj::FieldTransform, arg0::FieldVector3D)
    return jcall(obj, "transformVector", FieldVector3D, (FieldVector3D,), arg0)
end

function transform_vector(obj::FieldTransform, arg0::Vector3D)
    return jcall(obj, "transformVector", FieldVector3D, (Vector3D,), arg0)
end


text:
using Gcvspl, Test

# Dummy data
x = collect(0.:1:100)
y = 2*x

spline_coeff = 1.0
SplOrder = 3

c, WK, IER = Gcvspl.gcv_spl(x,y,ese,spline_coeff;SplineOrder = Int32(SplOrder-1)) # with cubic spline as the default
y_calc = Gcvspl.spl_der(x,x,c,SplineOrder= Int32(SplOrder-1))

@test isapprox(y,y_calc,atol=1e-1)

max_stars_repo_name: charlesll/gcvspl

max_issues_repo_name: charlesll/gcvspl

max_forks_repo_name: charlesll/gcvspl


text:
abstract type SKCE <: CalibrationErrorEstimator end

# kernel used by the estimator
kernel(skce::SKCE) = skce.kernel

# result type
function skce_result_type(skce::SKCE, predictions::AbstractMatrix{<:Real})
    s = zero(eltype(predictions))
    eltype(s^2 * zero(kernel_result_type(kernel(skce), s, s)))
end

"""
    skce_kernel(k, p, y, p̃, ỹ)

Evaluate kernel function
```math
(e_y - p)^T k(p, p̃) (e_ỹ - p̃)
```
of the squared kernel calibration error for the matrix-valued kernel `k` and predictions
`p` and `p̃` with corresponding labels `y` and `ỹ`.

This method assumes that `p`, `p̃`, `y`, and `ỹ` are valid and specified correctly, and
does not perform any checks.
"""
function skce_kernel(kernel::UniformScalingKernel, p::AbstractVector{<:Real},
                     y::Integer, p̃::AbstractVector{<:Real}, ỹ::Integer)
    @unpack λ = kernel

    # compute scalar products
    @inbounds s = dot(p, p̃) - p[ỹ] - p̃[y]
    y == ỹ && (s += 1)

    # multiply with scalar-valued kernel
    λ * s * kernel.kernel(p, p̃)
end

# evaluation of kernel function for diagonal kernels
function skce_kernel(kernel::DiagonalKernel, p::AbstractVector{<:Real},
                     y::Integer, p̃::AbstractVector{<:Real}, ỹ::Integer)
    @unpack diag = kernel

    # compute scalar products
    @inbounds s = dot(p, diag .* p̃) - p[ỹ] * diag[ỹ] - diag[y] * p̃[y]
    y == ỹ && (@inbounds s += diag[y])

    # multiply with scalar-valued kernel
    s * kernel.kernel(p, p̃)
end


text:
"""
# module HyperGraphTools

- Julia version:
- Author: bramb
- Date: 2019-02-13

# Examples

```jldoctest
julia>
```
"""
module HyperGraphTools
    include("simple_tagml_tokenizer.jl")
    include("visualization.jl")
    include("parser.jl")
    include("partiallyorderedparser.jl")
end

max_stars_repo_name: rhdekker/HyperGraphTools

max_issues_repo_name: rhdekker/HyperGraphTools

max_forks_repo_name: rhdekker/HyperGraphTools


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

using ..TestRequest

@testset "HTTP Stack Inserting" begin
    @testset "Insert - Beginning" begin
        expected = TestLayer{TopLayer{RedirectLayer{BasicAuthLayer{MessageLayer{RetryLayer{ExceptionLayer{ConnectionPoolLayer{StreamLayer{Union{}}}}}}}}}}
        result = insert(stack(), TopLayer, TestLayer)

        @test expected == result
    end

    @testset "Insert - Middle" begin
        expected = TopLayer{RedirectLayer{BasicAuthLayer{MessageLayer{RetryLayer{TestLayer{ExceptionLayer{ConnectionPoolLayer{StreamLayer{Union{}}}}}}}}}}
        result = insert(stack(), ExceptionLayer, TestLayer)

        @test expected == result
    end

    @testset "Insert - End" begin
        expected = TopLayer{RedirectLayer{BasicAuthLayer{MessageLayer{RetryLayer{ExceptionLayer{ConnectionPoolLayer{StreamLayer{TestLayer{Union{}}}}}}}}}}
        result = insert(stack(), Union{}, TestLayer)

        @test expected == result
    end

    @testset "Insert - Non-existant layer" begin
        @test_throws HTTP.Layers.LayerNotFoundException insert(stack(), AWS4AuthLayer, TestLayer)
    end

    @testset "Insert - Multiple Same layer" begin
        test_stack = insert(stack(), RetryLayer, ExceptionLayer)

        expected = TopLayer{RedirectLayer{BasicAuthLayer{MessageLayer{TestLayer{ExceptionLayer{RetryLayer{ExceptionLayer{ConnectionPoolLayer{StreamLayer{Union{}}}}}}}}}}}
        result = insert(test_stack, ExceptionLayer, TestLayer)

        @test expected == result
    end

    @testset "Inserted final layer runs handler" begin
        TestRequest.FLAG[] = false
        request(insert(stack(), Union{}, LastLayer), "GET", "https://httpbin.org/anything")
        @test TestRequest.FLAG[]
    end

    @testset "Insert/remove default layers" begin
        top = HTTP.top_layer(stack())
        insert_default!(top, TestLayer)
        @test HTTP.top_layer(stack()) <: TestLayer
        remove_default!(top, TestLayer)
        @test HTTP.top_layer(stack()) <: top
        insert_default!(Union{}, TestLayer)
        remove_default!(Union{}, TestLayer)
    end
end


text:
# A Simple Computational Algebra System

# scale
Base.:(-)(x::AbstractBlock{N}) where {N} = Scale(Val(-1), x)
Base.:(-)(x::Scale{Val{-1}}) = content(x)
Base.:(-)(x::Scale{Val{S}}) where {S} = Scale(Val(-S), content(x))
Base.:(-)(x::Scale) = Scale(-x.alpha, content(x))

Base.:(*)(x::AbstractBlock, α::Number) = α * x

# NOTE: ±,±im should be identical
Base.:(*)(α::Val{S}, x::AbstractBlock) where {S} = Scale(α, x)

function Base.:(*)(α::T, x::AbstractBlock) where {T<:Number}
    return α == one(T) ? x :
           α == -one(T) ? Scale(Val(-1), x) :
           α == im ? Scale(Val(im), x) : α == -im ? Scale(Val(-im), x) : Scale(α, x)
end

Base.:(*)(α::T, x::Scale) where {T<:Number} = α == one(T) ? x : Scale(x.alpha * α, content(x))
Base.:(*)(α::T, x::Scale{Val{S}}) where {T<:Number,S} = α * S * content(x)

Base.:(*)(α::Val{S}, x::Scale) where {S} = (S * x.alpha) * content(x)
Base.:(*)(α::Val{S1}, x::Scale{Val{S2}}) where {S1,S2} = (S1 * S2) * content(x)

Base.:(*)(x::Scale, y::Scale) = (x.alpha * y.alpha) * (content(x) * content(y))
Base.:(*)(x::Scale{Val{S1}}, y::Scale{Val{S2}}) where {S1,S2} = (S1 * S2) * (content(x) * content(y))
Base.:(*)(x::Scale, y::Scale{Val{S}}) where {S} = (x.alpha * S) * (content(x) * content(y))
Base.:(*)(x::Scale{Val{S}}, y::Scale) where {S} = (S * y.alpha) * (content(x) * content(y))
Base.:(*)(x::Scale, y::AbstractBlock) = x.alpha * chain(y, content(x))
Base.:(*)(y::AbstractBlock, x::Scale) = x.alpha * chain(content(x), y)

Base.:(+)(xs::AbstractBlock...) = Add(xs...)
Base.:(*)(xs::AbstractBlock...) = chain(Iterators.reverse(xs)...)
Base.:(/)(A::AbstractBlock, x::Number) = (1 / x) * A
# reduce
Base.prod(blocks::AbstractVector{<:AbstractBlock{N}}) where {N} = chain(Iterators.reverse(blocks)...)
Base.sum(blocks::AbstractVector{<:AbstractBlock{N}}) where {N} = +(blocks...)

Base.:(-)(lhs::AbstractBlock, rhs::AbstractBlock) = Add(lhs, -rhs)
Base.:(^)(x::AbstractBlock, n::Int) = chain((copy(x) for k in 1:n)...)


text:
# General approach to computing tracks from entity edge data.
using JLD

# Load the data file
file_dir = joinpath(Base.source_dir(),"../1EntityAnalysis/Entity.jld")
E = load(file_dir)["E"]
#E = loadassoc(file_dir)
Es = E
E = logical(E)

# Show general purpose method for building tracks.
p = StartsWith("PERSON/,")      # Set entity range.
t = StartsWith("TIME/,")        # Set time range.
x = StartsWith("LOCATION/,")    # Set spatial range.
a = StartsWith("PERSON/,TIME/,LOCATION/,")

# Limit to edges with all three.
E3 = E[getrow( sum(E[:,p],2) & sum(E[:,t],2) & sum(E[:,x],2) ),a]
printFull(E3[1:5,:])

# Collapse to get unique time and space for each edge and get triples.
edge,timeCols  = find(  val2col(col2type(E3[:,t],"/"),"/")  )
edge,space = find(  val2col(col2type(E3[:,x],"/"),"/")  )

Etx = Assoc(edge,timeCols,space)     # Combine edge, time and space.
Ext = Assoc(edge,space,timeCols)     # Combine edge, space and time.

# Construct time tracks with matrix multiply.
At = CatValMul(transpose(Etx),E3[:,p])
figure()
spy(At')
axis("auto")

# Construct space tracks with matrix multiply.
Ax = CatValMul(transpose(Ext),E3[:,p])
figure()
spy(Ax')
axis("auto")


text:
using Distances, LinearAlgebra
using Test
using KernelFunctions
using PDMats

dims = [10,5]

A = rand(dims...)
B = rand(dims...)
C = rand(8,9)
K = [zeros(dims[1],dims[1]),zeros(dims[2],dims[2])]
Kdiag = [zeros(dims[1]),zeros(dims[2])]
k = SqExponentialKernel()
@testset "Kernel Matrix Operations" begin
    @testset "Inplace Kernel Matrix" begin
        for obsdim in [1,2]
            @test kernelmatrix!(K[obsdim],k,A,B,obsdim=obsdim) == kernelmatrix(k,A,B,obsdim=obsdim)
            @test kernelmatrix!(K[obsdim],k,A,obsdim=obsdim) == kernelmatrix(k,A,obsdim=obsdim)
            @test kerneldiagmatrix!(Kdiag[obsdim],k,A,obsdim=obsdim) == kerneldiagmatrix(k,A,obsdim=obsdim)
            @test_throws DimensionMismatch kernelmatrix!(K[obsdim],k,A,C,obsdim=obsdim)
            @test_throws DimensionMismatch kernelmatrix!(K[obsdim],k,C,obsdim=obsdim)
            @test_throws DimensionMismatch kerneldiagmatrix!(Kdiag[obsdim],k,C,obsdim=obsdim)
        end
    end
    @testset "Kernel matrix" begin
        for obsdim in [1,2]
            @test kernelmatrix(k,A,B,obsdim=obsdim) == kappa.(k,pairwise(KernelFunctions.metric(k),A,B,dims=obsdim))
            @test kernelmatrix(k,A,obsdim=obsdim) == kappa.(k,pairwise(KernelFunctions.metric(k),A,dims=obsdim))
            @test kerneldiagmatrix(k,A,obsdim=obsdim) == diag(kernelmatrix(k,A,obsdim=obsdim))
            @test k(A,B,obsdim=obsdim) == kernelmatrix(k,A,B,obsdim=obsdim)
            @test k(A,obsdim=obsdim) == kernelmatrix(k,A,obsdim=obsdim)
            @test KernelFunctions._kernel(k,1.0,2.0) == KernelFunctions._kernel(k,[1.0],[2.0])
            @test_throws DimensionMismatch kernelmatrix(k,A,C,obsdim=obsdim)
        end
    end
    @testset "KernelSum" begin
        k1 = SqExponentialKernel()
        k2 = LinearKernel()
        ks = k1 + k2
        w1 = 0.4; w2 = 1.2;
        ks2 = KernelSum([k1,k2],weights=[w1,w2])
        @test all(kernelmatrix(ks,A) .== kernelmatrix(k1,A) + kernelmatrix(k2,A))
        @test all(kernelmatrix(ks+k1,A) .≈ 2*kernelmatrix(k1,A) + kernelmatrix(k2,A))
        @test all(kernelmatrix(k1+ks,A) .≈ 2*kernelmatrix(k1,A) + kernelmatrix(k2,A))
        @test all(kernelmatrix(ks,A,B) .== kernelmatrix(k1,A,B) + kernelmatrix(k2,A,B))
        @test all(kerneldiagmatrix(ks,A) .== kerneldiagmatrix(k1,A) + kerneldiagmatrix(k2,A))
        @test all(kernelmatrix(ks2,A) .== w1*kernelmatrix(k1,A) + w2*kernelmatrix(k2,A))
    end
    @testset "KernelProduct" begin
        k1 = SqExponentialKernel()
        k2 = LinearKernel()
        k3 = RationalQuadraticKernel()
        kp = k1 * k2
        kp2 = k1 * k3
        @test all(KernelFunctions.metric(kp).==[KernelFunctions.metric(k1),KernelFunctions.metric(k2)])
        @test all(kernelmatrix(kp,A) .≈ kernelmatrix(k1,A) .* kernelmatrix(k2,A))
        @test all(kernelmatrix(kp*k1,A) .≈ kernelmatrix(k1,A).^2 .* kernelmatrix(k2,A))
        @test all(kernelmatrix(k1*kp,A) .≈ kernelmatrix(k1,A).^2 .* kernelmatrix(k2,A))
        @test all(kernelmatrix(kp,A) .≈ kernelmatrix(k1,A) .* kernelmatrix(k2,A))
        @test all(kernelmatrix(kp,A,B) .≈ kernelmatrix(k1,A,B) .* kernelmatrix(k2,A,B))
        @test all(kernelmatrix(kp,A) .≈ kernelmatrix(k1,A) .* kernelmatrix(k2,A))
        @test all(kerneldiagmatrix(kp,A) .== kerneldiagmatrix(k1,A) .* kerneldiagmatrix(k2,A))
    end
    @testset "PDMat" begin
        for obsdim in [1,2]
            @test all(Matrix(kernelpdmat(k,A,obsdim=obsdim)) .≈ Matrix(PDMat(kernelmatrix(k,A,obsdim=obsdim))))
            # @test_throws ErrorException kernelpdmat(k,ones(100,100),obsdim=obsdim)
        end
    end
end


text:
# GaussDCA
# ========

"""
Wrapper function to `GaussDCA.gDCA`.
You need to install GaussDCA:
```julia
Pkg.clone("https://github.com/carlobaldassi/GaussDCA.jl")
```
And you need to load it before using this wrapper:
```julia
using GaussDCA
```
Look into [GaussDCA.jl README](https://github.com/carlobaldassi/GaussDCA.jl) for further information.
If you use this wrapper, **please cite the GaussDCA publication and the package's doi**.

**GaussDCA Publication:**
Baldassi, Carlo, Marco Zamparo, Christoph Feinauer, Andrea Procaccini, Riccardo Zecchina, Martin Weigt, and Andrea Pagnani.
"Fast and accurate multivariate Gaussian modeling of protein families: predicting residue contacts and protein-interaction partners."
PloS one 9, no. 3 (2014): e92721.
"""
function gaussdca(msa, args...; kargs...)
    filename = tempname()
    write(filename, msa, FASTA)
    try
        pairedvalues = Main.GaussDCA.gDCA(filename, args...; kargs...)
        plm = fill!(columnpairsmatrix(msa), NaN)
        for (i,j,value) in pairedvalues
           plm[i,j] = value
        end
        plm
    finally
        rm(filename)
    end
end


text:
function solution()
  n = 600851475143

  factor = 1
  while true
    factor += 2

    while n % factor == 0
      n = div(n, factor)
    end

    if factor * factor > n
      return n
    elseif n == 1
      return factor
    end
  end
end

if length(ARGS) == 1 && ARGS[1] == "-a"
  println(solution())
  exit(0)
end

for line in eachline(STDIN)
  iters = int(line)

  start = time_ns()
  for i in 1:iters
    solution()
  end
  end_ = time_ns()

  println(end_ - start)
end


text:
# MadNLP.jl
# Created by Sungho Shin (sungho.shin@wisc.edu)

module MadNLPRichardson

import ..MadNLP:
    @kwdef, Logger, @debug, @warn, @error,
    norm,
    AbstractOptions, AbstractIterator, set_options!, @sprintf, StrideOneVector,
    solve_refine!

@kwdef mutable struct Options <: AbstractOptions
    richardson_max_iter::Int = 10
    richardson_tol::Float64 = 1e-10
    richardson_acceptable_tol::Float64 = 1e-5
end

mutable struct Solver <: AbstractIterator
    res::Vector{Float64}
    mul!::Function
    div!::Function
    opt::Options
    logger::Logger
end
function Solver(res::Vector{Float64},mul!,div!;
                opt=Options(),logger=Logger(),
                option_dict::Dict{Symbol,Any}=Dict{Symbol,Any}(),kwargs...)
    !isempty(kwargs) && (for (key,val) in kwargs; option_dict[key]=val; end)
    set_options!(opt,option_dict)
    return Solver(res,mul!,div!,opt,logger)
end

function solve_refine!(x::StrideOneVector{Float64},
                       IS::Solver,
                       b::AbstractVector{Float64})
    @debug(IS.logger,"Iterative solver initiated")
    norm_b = norm(b,Inf)
    x.=0
    IS.res.=.-b
    norm_res = norm(IS.res,Inf)
    residual_ratio = norm_res/(1+norm_b)

    iter   = 0
    residual_ratio_old = Inf
    noprogress = 0

    while true
        mod(iter,10)==0 &&
            @debug(IS.logger,"iter ||res||")
        @debug(IS.logger,@sprintf("%4i %6.2e",iter,residual_ratio))
        iter += 1
        iter > IS.opt.richardson_max_iter && break
        residual_ratio < IS.opt.richardson_tol && break
        # residual_ratio>=residual_ratio_old && residual_ratio>IS.opt.richardson_acceptable_tol && (noprogress+=1)
        # if noprogress >= 3
        #     @debug(logger,@sprintf(
        #         "Iterative solver terminated with %4i refinement steps and residual = %6.2e",
        #         iter,residual_ratio))
        #     return :Singular
        # end

        IS.div!(IS.res)
        x.-=IS.res
        IS.mul!(IS.res,x)
        IS.res.-=b
        norm_res = norm(IS.res,Inf)

        residual_ratio_old = residual_ratio
        residual_ratio = norm_res/(1+norm_b)
    end

    @debug(IS.logger,@sprintf(
        "Iterative solver terminated with %4i refinement steps and residual = %6.2e",
        iter,residual_ratio))

    return (residual_ratio < IS.opt.richardson_acceptable_tol ? :Solved : :Failed)
end
end # module


text:
module Infiltrator

using REPL
using REPL.LineEdit

export @infiltrate

"""
    @infiltrate ex = true

Sets a "breakpoint" in a local context. Accepts an optional argument that must evaluate to a boolean,
and will drop you into an interactive REPL session that lets you inspect local variables and the call
stack as well as execute aribtrary statements in the context of the current function's module.

Usage:
```
julia> function f(x)
         x *= 2
         y = rand(3)
         @infiltrate
         x += 2
       end
f (generic function with 1 method)

julia> f(3)
Hit `@infiltrate` in f(::Int64) at none:4:

debug> ?
  Code entered is evaluated in the current function's module. Note that you cannot change local
  variables.
  The following commands are special cased:
    - `@trace`: Print the current stack trace.
    - `@locals`: Print local variables.
    - `@stop`: Stop infiltrating at this `@infiltrate` spot.

  Exit this REPL mode with `Ctrl-D`, and clear the effect of `@stop` with `Infiltrator.clear_stop()`.

debug> @trace
[1] f(::Int64) at none:4
[2] top-level scope at none:0

debug> @locals
- y::Array{Float64,1} = [0.187253, 0.145958, 0.183677]
- x::Int64 = 6

debug> x.+y
3-element Array{Float64,1}:
 6.187252565686353
 6.145958004935359
 6.1836766675450034

debug>

8

julia>
```
"""
macro infiltrate(ex = true)
  quote
    if $(esc(ex))
      start_prompt(@__MODULE__, Base.@locals, $(QuoteNode(ex)),
                   $(String(__source__.file)), $(__source__.line))
    end
  end
end

const prompt_color = if Sys.iswindows()
  "\e[33m"
else
  "\e[38;5;166m"
end

const TEST_TERMINAL_REF = Ref{Any}(nothing)
const TEST_REPL_REF = Ref{Any}(nothing)
const TEST_NOSTACK = Ref{Any}(false)

const STOP_SPOTS = Set()
clear_stop() = (empty!(STOP_SPOTS); nothing)

function start_prompt(mod, locals, ex, file, fileline;
                        terminal = TEST_TERMINAL_REF[],
                        repl = TEST_REPL_REF[],
                        nostack = TEST_NOSTACK[]
                      )
  (file, fileline) in STOP_SPOTS && return
  if terminal === nothing || repl === nothing
    if isdefined(Base, :active_repl) && isdefined(Base.active_repl, :t)
      repl = Base.active_repl
      terminal = Base.active_repl.t
    else
      println("Infiltrator.jl needs a proper Julia REPL.")
      return
    end
  end
  io = Base.pipe_writer(terminal)

  trace = stacktrace()
  start = something(findlast(x -> x.func === Symbol("start_prompt"), trace), 0) + 2
  last = something(findlast(x -> x.func === Symbol("eval"), trace),
                   length(trace)) - 2
  trace = trace[start:last]
  current = trace[1]
  println(io, "Hit `@infiltrate", ex == true ? "" : " $(ex)", "` in ", current, ":")
  println(io)
  debugprompt(mod, locals, trace, terminal, repl, nostack, file=file, fileline=fileline)
  println(io)
end

function show_help(io)
  println(io, """
    Code entered is evaluated in the current function's module. Note that you cannot change local
    variables.
    The following commands are special cased:
      - `@trace`: Print the current stack trace.
      - `@locals`: Print local variables.
      - `@stop`: Stop infiltrating at this `@infiltrate` spot.

    Exit this REPL mode with `Ctrl-D`, and clear the effect of `@stop` with `Infiltrator.clear_stop()`.
  """)
end

function show_trace(io, trace)
  for (i, frame) in enumerate(trace)
    println(io, "[", i, "] ", frame)
  end
  println(io)
end

function show_locals(io, locals)
  for (var, val) in locals
    println(io, "- ", var, "::", typeof(val), " = ", repr(val))
  end
  println(io)
end

function debugprompt(mod, locals, trace, terminal, repl, nostack = false; file, fileline)
  io = Base.pipe_writer(terminal)

  try
    panel = REPL.LineEdit.Prompt("debug> ";
              prompt_prefix = prompt_color,
              prompt_suffix = Base.text_colors[:normal],
              complete = InfiltratorCompletionProvider(mod, locals),
              on_enter = s -> true)

    panel.hist = REPL.REPLHistoryProvider(Dict{Symbol,Any}(:Infiltrator => panel))
    REPL.history_reset_state(panel.hist)

    search_prompt, skeymap = LineEdit.setup_search_keymap(panel.hist)
    search_prompt.complete = REPL.LatexCompletions()

    panel.on_done = (s, buf, ok) -> begin
      if !ok
        LineEdit.transition(s, :abort)
        REPL.LineEdit.reset_state(s)
        return false
      end

      line = String(take!(buf))

      if isempty(line)
        println(io)
        LineEdit.reset_state(s)
        return true
      end

      sline = strip(line)
      if sline == "?"
        show_help(io)
        LineEdit.reset_state(s)
        return true
      elseif sline == "@trace"
        show_trace(io, trace)
        LineEdit.reset_state(s)
        return true
      elseif sline == "@locals"
        show_locals(io, locals)
        LineEdit.reset_state(s)
        return true
      elseif sline == "@stop"
        push!(STOP_SPOTS, (file, fileline))
        println(io)
        LineEdit.reset_state(s)
        return true
      end
      ok = true
      result = nothing

      @static if VERSION >= v"1.2.0-DEV.253"
        try
          result = interpret(line, mod, locals)
        catch err
          ok = false
          result = Base.catch_stack()
          nostack && (result = map(r -> Any[first(r), []], result))
        end
        REPL.print_response(repl, (result, !ok), true, true)
      else
        try
          result = interpret(line, mod, locals)
        catch err
          ok = false
          result = (err, nostack ? Any[] : catch_backtrace())
        end
        REPL.print_response(repl, ok ? result : result[1], ok ? nothing : result[2], true, true)
      end
      println(io)
      LineEdit.reset_state(s)

      return true
    end

    panel.keymap_dict = LineEdit.keymap(Dict{Any,Any}[skeymap, LineEdit.history_keymap, LineEdit.default_keymap, LineEdit.escape_defaults])

    REPL.run_interface(terminal, REPL.LineEdit.ModalInterface([panel, search_prompt]))
  catch e
    e isa InterruptException || rethrow(e)
  end
end

maybe_quote(x) = (isa(x, Expr) || isa(x, Symbol)) ? QuoteNode(x) : x

function interpret(command::AbstractString, mod, locals)
    expr = Base.parse_input_line(command)
    Base.Meta.isexpr(expr, :toplevel) && (expr = expr.args[end])
    res = gensym()
    eval_expr = Expr(:let,
        Expr(:block, map(x->Expr(:(=), x...), [(k, maybe_quote(v)) for (k, v) in locals])...),
        Expr(:block,
            Expr(:(=), res, expr),
            Expr(:tuple, res, Expr(:tuple, [k for (k, v) in locals]...))
        ))
    eval_res, res = Core.eval(mod, eval_expr)
    eval_res
end

# completions

mutable struct InfiltratorCompletionProvider <: REPL.CompletionProvider
  mod::Module
  locals::Dict{Symbol, Any}
end

function LineEdit.complete_line(c::InfiltratorCompletionProvider, s)
  partial = REPL.beforecursor(s.input_buffer)
  full = LineEdit.input_string(s)

  ret, range, should_complete = completions(c, full, partial)
  return ret, partial[range], should_complete
end

function completions(c::InfiltratorCompletionProvider, full, partial)
  # repl backend completions
  comps, range, should_complete = REPL.REPLCompletions.completions(full, lastindex(partial), c.mod)
  ret = map(REPL.REPLCompletions.completion_text, comps)

  # local completions
  prepend!(ret, filter!(v -> startswith(v, partial), string.(keys(c.locals))))

  # Infiltrator commands completions
  commands = ["?", "@trace", "@locals", "@stop"]
  prepend!(ret, filter!(c -> startswith(c, partial), commands))

  unique!(ret), range, should_complete
end

end # module


text:
module Estimation

using LinearAlgebra: eigen, svd

export esprit

"""
    esprit(x::AbstractArray, M::Integer, p::Integer, Fs::Real=1.0)

ESPRIT [^Roy1986] algorithm for frequency estimation.
Estimation of Signal Parameters via Rotational Invariance Techniques

Given length N signal "x" that is the sum of p sinusoids of unknown frequencies,
estimate and return an array of the p frequencies.

# Arguments
- `x::AbstractArray`: complex length N signal array
- `M::Integer`: size of correlation matrix, must be <= N.
      The signal subspace is computed from the SVD of an M x (N-M+1) signal matrix
      formed from N-M+1 length-M shifts of the signal x in its columns.
      For best performance for 1 sinusoid, use M = (N+1)/3 (according to van der Veen and Leus).
      For faster execution (due to smaller SVD), use small M or small N-M
- `p::Integer`: number of sinusoids to estimate.
- `Fs::Float64`: sampling frequency, in Hz.

# Returns
length p real array of frequencies in units of Hz.

[^Roy1986]: R Roy, A Paulraj and T Kailath, ESPRIT - A subspace approach to estimation of parameters of cisoids in noise, IEEE Trans. Acoustics, Speech, Signal Process., 34, 1340-1342 (1986). [url](http://ieeexplore.ieee.org/abstract/document/1164935/).
"""
function esprit(x::AbstractArray, M::Integer, p::Integer, Fs::Real=1.0)
    count(v->v != 1, size(x)) <= 1 || error("`x` must be a 1D signal")
    N = length(x)
    X = x[ (1:M) .+ (0:N-M)' ]
    U,s,V = svd(X)
    D,_ = eigen( U[1:end-1,1:p] \ U[2:end,1:p] )

    angle.(D)*Fs/2π
end

end # end module definition


text:
#-----------------------------------------------------------------------
# MetadataTools
# https://github.com/IainNZ/MetadataTools.jl
#-----------------------------------------------------------------------
# Copyright (c) 2015: Iain Dunning
# Licensed under the MIT License
#-----------------------------------------------------------------------
# src/pkg_graph.jl
# Converts METADATA to a graph and provides some operations on that
# graph such as extracting [reverse] dependency graphs for a package.
#----------------------------------------------------------------------

export make_dep_graph, get_pkg_dep_graph

"""
    PkgGraph

A lightweight wrapper around a LightGraphs.jl representation of a
package dependency graph.
"""
type PkgGraph
    adjlist::Vector{Vector{Int}}
    pkgnames::Vector{Compat.UTF8String}
    pkgname_idx::Dict{Compat.UTF8String,Int}
end

"""
    size(pg::PkgGraph)

Returns the number of packages in a package dependency graph.
"""
Base.size(pg::PkgGraph) = length(pg.adjlist)

"""
    packages(pg::PkgGraph)

Returns the names of the packages in a package dependency graph
as a vector.
"""
packages(pg::PkgGraph) = copy(pkgnames)

"""
    adjlist(pkg::PkgGraph)

Returns a representation of a package dependency graph as a simple
adjacency list, using integer indices
"""
adjlist(pg::PkgGraph) = pg.adjlist

"""
    make_dep_graph(pkgs::Dict{UTF8String,PkgMeta}; reverse=false)

Given a Dict{UTF8String,PkgMeta} (e.g., from `get_all_pkg`), build a
directed graph with an edge from PkgA to PkgB iff PkgA directly requires
PkgB. Alternatively reverse the direction of the edges if `reverse` is true.
"""
function make_dep_graph(pkgs::Dict{Compat.UTF8String,PkgMeta}; reverse=false)
    pkgnames = collect(keys(pkgs))
    pkgname_idx = Dict{Compat.UTF8String,Int}(
        [(pkgname, i) for (i,pkgname) in enumerate(pkgnames)])  # TODO: Dict(pkgname=>i for ...)
    numpkg   = length(pkgnames)
    adjlist  = Vector{Int}[Vector{Int}() for pkgname in pkgnames]

    # Build up the edges of the graph
    for (pkg_name, pkg_meta) in pkgs
        if length(pkg_meta.versions) == 0
            # This package doesn't have any tagged versions, so
            # we can't figure out anything about what it depends
            # on from METADATA alone. Give up.
            continue
        end
        # We use the most recent version
        # TODO: add an option to modify this behaviour.
        pkg_ver = pkg_meta.versions[end]
        for req in pkg_ver.requires  # Requirements are plain strings
            if contains(req, "julia")
                # Julia version dependency, skip it
                continue
            end
            # Strip OS-specific conditions
            other_pkg = (req[1] == '@') ? split(req," ")[2] : split(req," ")[1]
            # Special case the the cycle:
            # LibCURL -> WinRPM (on Windows only)
            # WinRPM -> HTTPClient (on Unix only)
            # HTTPClient -> LibCurl
            # by breaking WinRPM -> HTTPClient
            if pkg_name == "WinRPM" && other_pkg == "HTTPClient"
                continue
            end
            # Map names to indices
            src, dst = pkgname_idx[pkg_name], pkgname_idx[other_pkg]
            # Add the directed edge
            if reverse
                dst, src = src, dst
            end
            if dst ∉ adjlist[src]
                push!(adjlist[src], dst)
            end
        end
    end

    return PkgGraph(adjlist, pkgnames, pkgname_idx)
end

"""
    get_pkg_dep_graph(pkg::PkgMeta, pg::PkgGraph; depth=Inf)
    get_pkg_dep_graph(pkgname::AbstractString, pg::PkgGraph; depth=Inf)

Returns a subgraph of the package dependency graph starting at a given
package. Optional keyword argument controls depth - depth=2 is just the
immediate dependencies.
"""
function get_pkg_dep_graph(pkg::PkgMeta, pg::PkgGraph; depth=Inf)
    get_pkg_dep_graph(pkg.name, pg, depth=depth)
end
function get_pkg_dep_graph(pkgname::AbstractString, pg::PkgGraph; depth=Inf)
    # Run a DFS to find the connected component
    # While doing so, building a mapping from old indices to
    # new indices in the subgraph
    n = size(pg)
    visited = zeros(Bool, n)
    stack = Int[pg.pkgname_idx[pkgname]]
    depth_stack = Int[1]
    old_idx_to_new_idx = Dict{Int,Int}()
    new_n = 0
    while length(stack) > 0
        cur_idx = pop!(stack)
        cur_depth = pop!(depth_stack)
        visited[cur_idx] && continue
        cur_depth > depth && continue
        visited[cur_idx] = true
        new_n += 1
        old_idx_to_new_idx[cur_idx] = new_n
        for dep_idx in pg.adjlist[cur_idx]
            if !visited[dep_idx]
                push!(stack, dep_idx)
                push!(depth_stack, cur_depth + 1)
            end
        end
    end
    # Build subgraph
    new_pkgnames = Compat.UTF8String["" for i in 1:new_n]
    new_pkgname_idx = Dict{Compat.UTF8String,Int}()
    new_adjlist  = Vector{Int}[Vector{Int}() for i in 1:new_n]
    for old_idx in 1:n
        !visited[old_idx] && continue
        new_idx = old_idx_to_new_idx[old_idx]
        new_pkgnames[new_idx] = pg.pkgnames[old_idx]
        new_pkgname_idx[pg.pkgnames[old_idx]] = new_idx
        for old_dep_idx in pg.adjlist[old_idx]
            if visited[old_dep_idx]  # might not have due to depth
                new_dep_idx = old_idx_to_new_idx[old_dep_idx]
                push!(new_adjlist[new_idx], new_dep_idx)
            end
        end
    end
    return PkgGraph(new_adjlist, new_pkgnames, new_pkgname_idx)
end


text:
@testset "Internal format strings" begin
    shortp(x) = repr(x, context = :color=>true)
    longp(x) = repr(:"text/plain", x, context = :color=>true)
    shortp_bw(x) = repr(x)
    longp_bw(x) = repr(:"text/plain", x)

    spuri = SpUri()
    spid = SpId()
    spcategoryid = SpCategoryId()
    spuserid = SpUserId()
    spurl = SpUrl()
    spidvec = [SpId(), SpId()]

    @test shortp(spuri)[1] == '\e'
    @test longp(spuri)[1] == '\e'
    @test shortp_bw(spuri)[1] == 's'
    @test longp_bw(spuri)[1] == '"'

    @test shortp(spid)[1] == '\e'
    @test longp(spid)[1] == '\e'
    @test shortp_bw(spid)[1] == '6'
    @test longp_bw(spid)[1] == '"'

    @test shortp(spcategoryid)[1] == '\e'
    @test longp(spcategoryid)[1] == '\e'
    @test shortp_bw(spcategoryid)[1] == 'p'
    @test longp_bw(spcategoryid)[1] == '"'

    @test shortp(spuserid)[1] == '\e'
    @test longp(spuserid)[1] == '\e'

    # This test case depends on the user ID info from .ini file
    #@test shortp_bw(spuserid)[1] == 'S'
    @test longp_bw(spuserid)[1] == '"'

    @test shortp(spurl)[1] == '\e'
    @test longp(spurl)[1] == '\e'
    @test shortp_bw(spurl)[1] == 'h'
    @test longp_bw(spurl)[1] == '"'

    @test "AA|$(spidvec)|BB" ==  "AA|6rqhFgbbKwnb9MLmUQDhG6, 6rqhFgbbKwnb9MLmUQDhG6|BB"
end


text:
mutable struct SceneBoot <: Ranger.AbstractScene
    base::Nodes.NodeData

    replacement::Ranger.AbstractScene

    function SceneBoot(world::Ranger.World, name::String, replacement::Ranger.AbstractScene)
        obj = new()

        obj.base = Nodes.NodeData(name, Nodes.NodeNil(), world)
        obj.replacement = replacement

        obj
    end
end

# --------------------------------------------------------
# Transitioning
# --------------------------------------------------------
function Nodes.transition(node::SceneBoot)
    Scenes.REPLACE_TAKE
end

function Nodes.get_replacement(node::SceneBoot)
    node.replacement
end

max_stars_repo_name: wdevore/Ranger-Julia-SDL

max_issues_repo_name: wdevore/Ranger-Julia-SDL

max_forks_repo_name: wdevore/Ranger-Julia-SDL


text:
using Revise
using Gadfly
using ANOVADDPTest
using DataFrames
using Distributions
using Random

# Here, I simulate data
rng = MersenneTwister(1); # seed
N = 1000; # number of observations
D = 2; # number of covariates
X = rand(rng, 1:3, N, D);
y = randn(rng, N); # clearly, all groups are equal
for i in 1:N
    ((X[i, 1] != 1) && (X[i, 2] == 1)) && (y[i] += 1.0)
end

# Organize the data
data = NormalData(X, y);

# Create some data for prediction (useful for density plots)
pred = expandgrid([1:3 for _ in 1:D]...,  range(-2.0, stop = 2.0, length = 50));
pred = NormalData(hcat(pred[1:D]...), pred[D+1]);

# Create the model
m = NormalDDP(rng, N, length(data.Xunique));

# Fit the model
ch = train(rng, m, data, pred);

# Compute the effects
simple_effect_probabilities(ch, data)
interaction_effect_probabilities(ch, data)


text:
## Extensions to Polynomials.jl

" Return the variable `x` of a polynomial as a polynomial."
variable(p::Poly) = poly(zeros(eltype(p),1), p.var)

"Create a monic polynomial from `p`"
monic(p::Poly) = Poly(p.a/p[degree(p)], p.var)

## Reverse coefficients
rcoeffs(p::Poly) = reverse(coeffs(p))

function Base.convert{T<:Integer}(::Type{Poly{T}}, p::Poly{Rational{T}})
    l = reduce(lcm, [x.den for x in p.a])
    q = l * p
    Poly(T[x for x in q.a], p.var)
end

## convert Poly <--> function
Base.convert(::Type{Function}, p::Poly) = x -> Polynomials.polyval(p,x)
## convert a function to a polynomial with error if conversion is not possible
## This needs a hack, as Polynomials.jl defines `/` to return a `div` and not an error for p(x)/q(x)
immutable PolyTest
   x
end
import Base: +,-,*,/,^
+(a::PolyTest, b::PolyTest) = PolyTest(a.x + b.x)
+{T<:Number}(a::T, b::PolyTest) = PolyTest(a + b.x)
+{T<:Number}(a::PolyTest,b::T) = PolyTest(a.x + b)
-(a::PolyTest,b::PolyTest) = PolyTest(a.x - b.x)
-{T<:Number}(a::T, b::PolyTest) = PolyTest(a - b.x)
-{T<:Number}(a::PolyTest,b::T) = PolyTest(a.x - b)
-(a::PolyTest) = PolyTest(-a.x)
*(a::PolyTest, b::PolyTest) = PolyTest(a.x * b.x)
*(a::Bool, b::PolyTest) = PolyTest(a * b.x)
*{T<:Number}(a::T, b::PolyTest) = PolyTest(a * b.x)
*{T<:Number}(b::PolyTest, a::T) = PolyTest(a * b.x)
/{T<:Number}(b::PolyTest, a::T) = PolyTest(b.x / a)
^{T<:Integer}(b::PolyTest, a::T) = PolyTest(b.x ^ a)

typealias QQR @compat Union{Int, BigInt, Rational{Int}, Rational{BigInt}, Float64}
function Base.convert{T<:QQR}(::Type{Poly{T}}, f::Function)
    try
        f(PolyTest(0))                  # error if not a poly
        x = poly(zeros(T,1))
        out = f(x)
        if !isa(out, Poly)
            out = Poly([out])   # maybe a constant
        end
        out
    catch e
        rethrow(e)
    end
end
function Base.convert(::Type{Poly}, f::Function)
    ## try integers first, then float
    for T in [BigInt, Int, Float64]
        try
            fn = convert(Poly{T}, f)
            return(fn)
        catch e
        end
    end
    DomainError()
end

*{T, S}(A::Array{T,2}, p::Poly{S}) = Poly(A * rcoeffs(p))


text:
using CairoMakie
using Random, Distributions
using SpecialFunctions

α = 1.0
β = 3.0

function beta_distribution(x, α, β)
    return gamma(α+β) / (gamma(α) * gamma(β)) * x^(α-1) * (1 - x)^(β - 1)
end

x = range(0, 1, length=100)
y = beta_distribution.(x, α, β)
l = lines(x, y, color = :grey, axis = (xminorticks = IntervalsBetween(4),xminorticksvisible = true,
yminorticksvisible = true,
xminorgridvisible = true,
yminorgridvisible = true,))
xlims!(0, nothing)
ylims!(0, nothing)

save("beta-alpha$(α)-beta$β.pdf", l)


text:
using CartesianGP
using Base.Test

# TODO: Add correctness tests

# Half adder

function parameters_half_adder()
    numperlevel = 1
    numlevels = 2
    numlevelsback = numlevels
    numinputs = 2
    numoutputs = 2

    mu = 1
    lambda = 4
    mutrate = 0.05
    targetfitness = 1.0
    funcs = [AND, OR, XOR]
    fitfunc = hamming_max

    return Parameters(mu, lambda, mutrate, targetfitness, numinputs, numoutputs, numperlevel, numlevels, numlevelsback, funcs, fitfunc)
end

function goal_half_adder()
    return Goal(2, (0b0110, 0b1000))
end

p = parameters_half_adder()
g = goal_half_adder()
n = 1000
r = mu_lambda(p, g, n)

println("running mu_lambda on half adder")
println("fitness:",fitness(r[1], g),"  gens:",r[2])
@assert fitness(r[1], g) == p.targetfitness
@assert r[2] <= n
print_chromosome(r[1])

#Half (single-output) full adder
function parameters_half_full_adder()
    numperlevel = 1
    numlevels = 10
    numlevelsback = numlevels
    numinputs = 3
    numoutputs = 1

    mu = 1
    lambda = 4
    mutrate = 0.12
    targetfitness = 1.0
    funcs = [AND, OR, XOR]
    fitfunc = hamming_max

    return Parameters(mu, lambda, mutrate, targetfitness, numinputs, numoutputs, numperlevel, numlevels, numlevelsback, funcs, fitfunc)
end

function goal_half_full_adder()
    return Goal(3, (0xe8, ))
end

p = parameters_half_full_adder()
g = goal_half_full_adder()
n = 100000
println("runing mu_lambda on half full adder")
r = mu_lambda(p, g, n)
println("generations: ",r[2])
#println("trying evolve on half full adder")
#r = evolve(p, g, n)

println("fitness:",fitness(r[1], g),"  gens:",r[2])
@assert fitness(r[1], g) == p.targetfitness
@assert r[2] <= n
print_chromosome(r[1])

# Full adder

function parameters_full_adder()
    numperlevel = 1
    numlevels = 4
    numlevelsback = numlevels
    numinputs = 3
    numoutputs = 2

    mu = 1
    lambda = 4
    mutrate = 0.1
    targetfitness = 1.0
    funcs = [AND, OR, XOR]
    fitfunc = hamming_max

    return Parameters(mu, lambda, mutrate, targetfitness, numinputs, numoutputs, numperlevel, numlevels, numlevelsback, funcs, fitfunc)
end

function goal_full_adder()
    return Goal(3, (0b10010110, 0b11101000))
end

p = parameters_full_adder()
g = goal_full_adder()
n = 1000000
println("runing mu_lambda on full adder with n: ",n)
r = mu_lambda(p, g, n)
println("generations: ",r[2])
@assert fitness(r[1], g) == p.targetfitness
@assert r[2] <= n
print_chromosome(r[1])


text:
abstract type InfoDist end

"""
    observe!(d, args...)
"""
observe!(::InfoDist, args...)

"""
    clear!(d)
"""
clear!(::InfoDist)

"""
    estimate(d)
"""
estimate(::InfoDist)


text:
"""
`Basics.merge()`  ... to merge data from different instances into a single one.

  + `("atomic bases", bases::Array{Basis,1})`
    ... to merge two (or more) atomic bases into a single basis::Basis. This method assumes the same number of electrons in all basis and
        that the subshell lists are the same or can be made 'consistent' to each other. Two bases have a consistent subshell list if all
        subshells, that appear in any of the two lists appear always in the same sequence (if they are not missing at all). In the merged
        basis, the radial orbitals are taken from the basis (in the bases-array}, from where they are found first.
"""
function Basics.merge(sa::String, bases::Array{Basis,1})
    sa != "atomic bases"   &&   error("Unsupported keystring = $sa")

    if  length(bases) > 2
        println("Number of bases = $(length(bases)) ...")
        bs       = Basics.merge("atomic bases", [ bases[1], bases[2] ])
        basesNew = [bs]
        for i = 3:length(bases)    push!( basesNew, bases[i])    end
        Basics.merge("atomic bases", basesNew )
    elseif  length(bases) == 1    return( basis[1] )
    elseif  length(bases) == 2
        # This is the essential step here; first create a consistent subshell list
        basisA = bases[1];    basisB = bases[2]
        wa = deepcopy( basisA.subshells );   wb = [ 100i  for i = 1:length(wa) ]
        # Add the subshells of basisB
        for  i = 1:length( basisB.subshells )
            if    basisB.subshells[i] in wa
            else  push!( wa, basisB.subshells[i] );    push!( wb, wb[i-1]+1 )
            end
        end
        # Now sort wb and create a new subshell list from wa due to this sorting
        wc = sortperm( wb )
        newSubshells = Subshell[]
        for  i in wc   push!( newSubshells, wa[i] )    end
        println("Sorted subshell list: $(newSubshells) ")
        error("Not yet fully implemented, see source code ")
        # Transform all CSF from basisA and basisB to newSubshells ... and apply CsfRExcludeDouble()
        # return new basis
        #
    else    error("stop a")
    end
end

"""
  + `("multiplets", multiplets::Array{Multiplet,1})`
    ... to merge two (or more) atomic multiplets into a single multiplet::Multiplet. This method assumes (and checks) that all
        levels have level.hasStateRep = true and that all levels refer to the same basis.
"""
function Basics.merge(sa::String, multiplets::Array{Multiplet,1})
    sa != "multiplets"   &&   error("Unsupported keystring = $sa")
    if      length(multiplets) == 1    return( multiplets[1] )
    elseif  length(multiplets) == 2
        levels = copy(multiplets[1].levels);   basis = multiplets[1].levels[1].basis
        for  lev  in  multiplets[2].levels
            if  !(lev.hasStateRep)   ||   lev.basis != basis    error("Levels of multiplets do not refer to the same basis.")   end
            push!(levels, lev)
        end
        mp = Multiplet( multiplets[1].name * "+" * multiplets[2].name, levels )
        return( mp )
    elseif  length(multiplets) > 2
        mp = Basics.merge("multiplets", [ multiplets[1], multiplets[2] ])
        for  k = 3:length(multiplets)
            mp = Basics.merge("multiplets", [ mp, multiplets[k] ])
        end
        return( mp )
    else    error("stop a")
    end
end