text:
# Classifies MNIST digits with a convolutional network.
# Writes out saved model to the file "mnist_conv.bson".
# Demonstrates basic model construction, training, saving,
# conditional early-exit, and learning rate scheduling.
#
# This model, while simple, should hit around 99% test
# accuracy after training for approximately 20 epochs.

using Flux, Flux.Data.MNIST, Statistics
using Flux: onehotbatch, onecold, crossentropy, throttle
using Base.Iterators: repeated, partition
using Printf, BSON
using Base.Iterators
using LinearAlgebra
using ToyFHE
using Primes
using GaloisFields
using Test

using OffsetArrays

# Prepare test set in batches of 64
function make_minibatch(X, Y, idxs)
    X_batch = Array{Float32}(undef, size(X[1])..., 1, length(idxs))
    for i in 1:length(idxs)
        X_batch[:, :, :, i] = Float32.(X[idxs[i]])
    end
Y_batch = onehotbatch(Y[idxs], 0:9)
    return (X_batch, Y_batch)
end
batch_size = 64
test_imgs = MNIST.images(:test)
test_labels = MNIST.labels(:test)
mb_idxs = partition(1:length(test_labels), batch_size)
test_set =  [make_minibatch(test_imgs, test_labels, i) for i in mb_idxs]

# Load up the model we trained in train.jl
BSON.@load joinpath(dirname(@__FILE__), "mnist_conv.bson") model

function reshape_and_vcat(x)
    let y=reshape(x, 64, 4, size(x, 4))
        vcat((y[:,i,:] for i=axes(y,2))...)
    end
end

model = Chain(model.layers[1], reshape_and_vcat, model.layers[3:end]...)

function encrypted_matmul(weights, x)
    sum(diag(circshift(weights, (0,(k-1)))).*circshift(x, (k-1,0)) for k = 1:size(x, 1))
end

function naive_rectangular_matmul(weights, x)
    @assert size(weights, 1) < size(weights, 2)
    weights = vcat(weights, zeros(eltype(weights), size(weights, 2)-size(weights, 1), size(weights, 2)))
    encrypted_matmul(weights, x)
end

function public_preprocess(batch)
    ka = OffsetArray(0:7, 0:7)
    # Create feature extracted matrix
    I = [[batch[i′*3 .+ (1:7), j′*3 .+ (1:7), 1, k] for i′=ka, j′=ka] for k = 1:64]

    # Reshape into the ciphertext
    Iᵢⱼ = [[I[k][l...][i,j] for k=1:64, l=product(ka, ka)] for i=1:7, j=1:7]
end

function do_encrypted_inference(model, batch)
    Iᵢⱼ = public_preprocess(batch)

    # Evaluate the convolution
    weights = model.layers[1].weight
    conv_weights = reverse(reverse(weights, dims=1), dims=2)
    conved = [sum(Iᵢⱼ[i,j]*conv_weights[i,j,1,channel] for i=1:7, j=1:7) for channel = 1:4]
    conved = map(((x,b),)->x .+ b, zip(conved, model.layers[1].bias))
    sqed1 = map(x->x.^2, conved)
    sqed1 = map(x->reshape(x, 64, 64)', sqed1)
    fq1_weights = model.layers[3].W
    fq1 = sum(enumerate(partition(1:256, 64))) do (i,range)
        encrypted_matmul(fq1_weights[:, range], sqed1[i])
    end
    fq1 .+= model.layers[3].b
    sqed2 = fq1.^2
    fq2_weights = model.layers[4].W
    result = naive_rectangular_matmul(fq2_weights, sqed2)[1:10, :]
    result .+= model.layers[4].b
    result
end

batch = test_set[1][1]

# First test that the plaintext implementation above is correct
@test model(batch) ≈ do_encrypted_inference(model, batch)

# Set up the crypto parameters
# For now, we match the parameters of [JKLS19] and see how far we get

N = 2^13
q₀ = nextprime(Int128(2)^60 + 1, 1; interval=2N)
ps = nextprime(q₀ + 2N, 1; interval=2N)

q₁ = nextprime(2^40 + 1, 1; interval=2N)
q₂ = nextprime(q₁ + 2N, 1; interval=2N)
q₃ = nextprime(q₂ + 2N, 1; interval=2N)
q₄ = nextprime(q₃ + 2N, 1; interval=2N)
q₅ = nextprime(q₄ + 2N, 1; interval=2N)

ℛ = let CT = ToyFHE.CRTEncoded{7, Tuple{GaloisField.((q₀,q₁,q₂,q₃,q₄,q₅,ps))...}}
    ζ₂n = GaloisFields.minimal_primitive_root(CT, 2N)
    ToyFHE.NegacyclicRing{CT, N}(ζ₂n)
end

ckks_params = ModulusRaised(CKKSParams(ℛ, 0, 3.2))
kp = keygen(ckks_params)

Iᵢⱼ = public_preprocess(batch)

scale = Int128(2^40)
Tscale = FixedRational{scale}

C_Iij = map(Iᵢⱼ) do Iij
    plain = CKKSEncoding{Tscale}(zero(plaintext_space(ckks_params)))
    plain .= OffsetArray(vec(Iij), 0:(N÷2-1))
    encrypt(kp, plain)
end

weights = model.layers[1].weight
conv_weights = reverse(reverse(weights, dims=1), dims=2)
conved3 = [sum(C_Iij[i,j]*conv_weights[i,j,1,channel] for i=1:7, j=1:7) for channel = 1:4]
conved2 = map(((x,b),)->x .+ b, zip(conved3, model.layers[1].bias))
conved1 = map(ToyFHE.modswitch, conved2)

ek = keygen(EvalMultKey, kp.priv)
gk = keygen(GaloisKey, kp.priv; steps=64)

Csqed1 = map(x->x*x, conved1)
Csqed1 = map(x->keyswitch(ek, x), Csqed1)
Csqed1 = map(ToyFHE.modswitch, Csqed1)

decrypt(kp, Csqed1[1])

function encrypted_matmul(gk, weights, x::ToyFHE.CipherText)
    result = repeat(diag(weights), inner=64).*x
    rotated = x
    for k = 2:64
        @show k
        rotated = ToyFHE.rotate(gk, rotated)
        result += repeat(diag(circshift(weights, (0,(k-1)))), inner=64) .* rotated
    end
    result
end

decrypt_matrix(kp, xx) = reshape(collect(decrypt(kp, xx)), (64, 64))'

fq1_weights = model.layers[3].W
reshape(collect(decrypt(kp, encrypted_matmul(gk, fq1_weights[:, 1:64], Csqed1[1]))), (64, 64))'
fq1_weights[:, 1:64] * reshape(collect(decrypt(kp, Csqed1[1])), (64, 64))'

Cfq1 = sum(enumerate(partition(1:256, 64))) do (i,range)
    encrypted_matmul(gk, fq1_weights[:, range], Csqed1[i])
end

Cfq1 = Cfq1 .+ OffsetArray(repeat(model.layers[3].b, inner=64), 0:4095)
Cfq1 = modswitch(Cfq1)

Csqed2 = Cfq1*Cfq1
Csqed2 = keyswitch(ek, Csqed2)
Csqed2 = modswitch(Csqed2)

function naive_rectangular_matmul(gk, weights, x)
    @assert size(weights, 1) < size(weights, 2)
    weights = vcat(weights, zeros(eltype(weights), size(weights, 2)-size(weights, 1), size(weights, 2)))
    encrypted_matmul(gk, weights, x)
end

fq2_weights = model.layers[4].W
Cresult = naive_rectangular_matmul(gk, fq2_weights, Csqed2)
Cresult = Cresult .+ OffsetArray(repeat(vcat(model.layers[4].b, zeros(54)), inner=64), 0:4095)

plain_result = model(batch)
enc_result = real.(decrypt_matrix(kp, Cresult))

# N.B.: To see the images in your terminal, make sure to have `Images` and `TerminalExtensions` loaded,
# and make sure your terminal supports inline images (currently iTerm2 on OS X).
function compare_models(test_imgs, test_labels, plain_result, enc_result, plain_label, enc_label)
    print(" "^14)
    display(reduce(hcat, test_imgs[1:64]))
    println()
    print("Ground truth:  ")
    print(join(test_labels[1:64], " "))
    println()
    print(rpad(plain_label, 15))
    for idx = 1:64
        label = argmax(plain_result[:,idx])-1
        printstyled(label, color = label == test_labels[idx] ? :green : red)
        print(" ")
    end
    println()
    print(rpad(enc_label, 15))
    for idx = 1:64
        label = argmax(enc_result[:,idx])-1
        printstyled(label, color = label == test_labels[idx] ? :green : red)
        print(" ")
    end
    println()
end

compare_models(test_imgs, test_labels, model(batch), real.(decrypt_matrix(kp, Cresult)), "Model (Plain):", "Model (Enc):")


text:
function noop_sampler(rawdata)
    @warn "Noop data sampler, returning nothing."
    return nothing
end

ext: jl


text:
###############################################################################
#
# Initialize CFCache
#
###############################################################################
"""
    create_cf_cache(FT, rt_dim::RTDimensions)

Create a [`CFCache`](@ref) type struct, given
- `FT` Floating number type
- `rt_dim` [`RTDimensions`](@ref) type struct
"""
function create_cf_cache(FT, rt_dim::RTDimensions)
    @unpack nAzi, nLayer, nPAR, nWL = rt_dim;

    abs_wave    = zeros(FT, nWL   );
    absfs_lidf  = zeros(FT, nAzi  );
    E_all       = zeros(FT, nWL   );
    E_iPAR      = zeros(FT, nPAR  );
    lPs         = zeros(FT, nLayer);
    kChlrel     = zeros(FT, nPAR  );
    PAR_diff    = zeros(FT, nPAR  );
    PAR_diffCab = zeros(FT, nPAR  );
    PAR_dir     = zeros(FT, nPAR  );
    PAR_dirCab  = zeros(FT, nPAR  );

    return CFCache{FT}(abs_wave    = abs_wave   ,
                       absfs_lidf  = absfs_lidf ,
                       E_all       = E_all      ,
                       E_iPAR      = E_iPAR     ,
                       lPs         = lPs        ,
                       kChlrel     = kChlrel    ,
                       PAR_diff    = PAR_diff   ,
                       PAR_diffCab = PAR_diffCab,
                       PAR_dir     = PAR_dir    ,
                       PAR_dirCab  = PAR_dirCab )
end

###############################################################################
#
# Initialize CGCache
#
###############################################################################
"""
    create_cg_cache(FT, rt_dim::RTDimensions)

Create a [`CGCache`](@ref) type struct, given
- `FT` Floating number type
- `rt_dim` [`RTDimensions`](@ref) type struct
"""
function create_cg_cache(FT, rt_dim::RTDimensions)
    @unpack nAzi, nIncl = rt_dim;

    _Co  = zeros(FT, nIncl);
    _Cs  = zeros(FT, nIncl);
    _So  = zeros(FT, nIncl);
    _Ss  = zeros(FT, nIncl);

    _1s  =  ones(FT, (1    , nAzi));
    _2d  = zeros(FT, (nIncl, nAzi));
    _cdo = zeros(FT, (nIncl, nAzi));
    _cds = zeros(FT, (nIncl, nAzi));

    return CGCache{FT}(_Co  = _Co ,
                       _Cs  = _Cs ,
                       _So  = _So ,
                       _Ss  = _Ss ,
                       _1s  = _1s ,
                       _2d  = _2d ,
                       _cdo = _cdo,
                       _cds = _cds)
end

###############################################################################
#
# Initialize SWCache
#
###############################################################################
"""
    create_sf_cache(FT, rt_dim::RTDimensions)

Create a [`SFCache`](@ref) type struct, given
- `FT` Floating number type
- `rt_dim` [`RTDimensions`](@ref) type struct
"""
function create_sf_cache(FT, rt_dim::RTDimensions)
    @unpack nAzi, nIncl, nLayer, nLevel, nWLE, nWLF = rt_dim;

    τ_dd            = zeros(FT, (nWLF,nLayer));
    Xdd             = zeros(FT, (nWLF,nLayer));
    Rdd             = zeros(FT, (nWLF,nLevel));
    Y               = zeros(FT, (nWLF,nLayer));
    U               = zeros(FT, (nWLF,nLevel));
    dnorm           = zeros(FT, nWLF);
    ρ_dd            = zeros(FT, (nWLF,nLayer));
    S⁻              = zeros(FT, (nWLF,nLayer));
    S⁺              = zeros(FT, (nWLF,nLayer));
    piLs            = zeros(FT, (nWLF,nLayer));
    piLd            = zeros(FT, (nWLF,nLayer));
    Fsmin           = zeros(FT, (nWLF,nLayer));
    Fsplu           = zeros(FT, (nWLF,nLayer));
    Fdmin           = zeros(FT, (nWLF,nLayer));
    Fdplu           = zeros(FT, (nWLF,nLayer));
    Femo            = zeros(FT, (nWLF,nLayer));
    M⁺              = zeros(FT, (nWLF,nWLE));
    M⁻              = zeros(FT, (nWLF,nWLE));
    M⁻_sun          = zeros(FT, nWLF);
    M⁺_sun          = zeros(FT, nWLF);
    wfEs            = zeros(FT, nWLF);
    sfEs            = zeros(FT, nWLF);
    sbEs            = zeros(FT, nWLF);
    M⁺⁻             = zeros(FT, nWLF);
    M⁺⁺             = zeros(FT, nWLF);
    M⁻⁺             = zeros(FT, nWLF);
    M⁻⁻             = zeros(FT, nWLF);
    sun_dwl_iWlE    = zeros(FT, nWLE);
    tmp_dwl_iWlE    = zeros(FT, nWLE);
    ϕ_cosΘ          = zeros(FT, (nIncl,nAzi));
    ϕ_cosΘ_lidf     = zeros(FT, nAzi);
    vfEplu_shade    = zeros(FT, nWLF);
    vbEmin_shade    = zeros(FT, nWLF);
    vfEplu_sun      = zeros(FT, nWLF);
    vbEmin_sun      = zeros(FT, nWLF);
    sigfEmin_shade  = zeros(FT, nWLF);
    sigbEmin_shade  = zeros(FT, nWLF);
    sigfEmin_sun    = zeros(FT, nWLF);
    sigbEmin_sun    = zeros(FT, nWLF);
    sigfEplu_shade  = zeros(FT, nWLF);
    sigbEplu_shade  = zeros(FT, nWLF);
    sigfEplu_sun    = zeros(FT, nWLF);
    sigbEplu_sun    = zeros(FT, nWLF);
    zeroB           = zeros(FT, nWLF);
    F⁻              = zeros(FT, (nWLF,nLevel));
    F⁺              = zeros(FT, (nWLF,nLevel));
    net_diffuse     = zeros(FT, (nWLF,nLayer));
    tmp_1d_nWlF     = zeros(FT, nWLF);
    tmp_1d_nLayer   = zeros(FT, nLayer);
    tmp_2d_nWlF_nLayer   = zeros(FT, (nWLF,nLayer));
    tmp_2d_nWlF_nLayer_2 = zeros(FT, (nWLF,nLayer));

    return SFCache{FT}(τ_dd                 = τ_dd                ,
                       Xdd                  = Xdd                 ,
                       Rdd                  = Rdd                 ,
                       Y                    = Y                   ,
                       U                    = U                   ,
                       dnorm                = dnorm               ,
                       ρ_dd                 = ρ_dd                ,
                       S⁻                   = S⁻                  ,
                       S⁺                   = S⁺                  ,
                       piLs                 = piLs                ,
                       piLd                 = piLd                ,
                       Fsmin                = Fsmin               ,
                       Fsplu                = Fsplu               ,
                       Fdmin                = Fdmin               ,
                       Fdplu                = Fdplu               ,
                       Femo                 = Femo                ,
                       M⁺                   = M⁺                  ,
                       M⁻                   = M⁻                  ,
                       M⁻_sun               = M⁻_sun              ,
                       M⁺_sun               = M⁺_sun              ,
                       wfEs                 = wfEs                ,
                       sfEs                 = sfEs                ,
                       sbEs                 = sbEs                ,
                       M⁺⁻                  = M⁺⁻                 ,
                       M⁺⁺                  = M⁺⁺                 ,
                       M⁻⁺                  = M⁻⁺                 ,
                       M⁻⁻                  = M⁻⁻                 ,
                       sun_dwl_iWlE         = sun_dwl_iWlE        ,
                       tmp_dwl_iWlE         = tmp_dwl_iWlE        ,
                       ϕ_cosΘ               = ϕ_cosΘ              ,
                       ϕ_cosΘ_lidf          = ϕ_cosΘ_lidf         ,
                       vfEplu_shade         = vfEplu_shade        ,
                       vbEmin_shade         = vbEmin_shade        ,
                       vfEplu_sun           = vfEplu_sun          ,
                       vbEmin_sun           = vbEmin_sun          ,
                       sigfEmin_shade       = sigfEmin_shade      ,
                       sigbEmin_shade       = sigbEmin_shade      ,
                       sigfEmin_sun         = sigfEmin_sun        ,
                       sigbEmin_sun         = sigbEmin_sun        ,
                       sigfEplu_shade       = sigfEplu_shade      ,
                       sigbEplu_shade       = sigbEplu_shade      ,
                       sigfEplu_sun         = sigfEplu_sun        ,
                       sigbEplu_sun         = sigbEplu_sun        ,
                       zeroB                = zeroB               ,
                       F⁻                   = F⁻                  ,
                       F⁺                   = F⁺                  ,
                       net_diffuse          = net_diffuse         ,
                       tmp_1d_nWlF          = tmp_1d_nWlF         ,
                       tmp_1d_nLayer        = tmp_1d_nLayer       ,
                       tmp_2d_nWlF_nLayer   = tmp_2d_nWlF_nLayer  ,
                       tmp_2d_nWlF_nLayer_2 = tmp_2d_nWlF_nLayer_2)
end

###############################################################################
#
# Initialize SWCache
#
###############################################################################
"""
    create_sw_cache(FT, rt_dim::RTDimensions)

Create a [`CGCache`](@ref) type struct, given
- `FT` Floating number type
- `rt_dim` [`RTDimensions`](@ref) type struct
"""
function create_sw_cache(FT, rt_dim::RTDimensions)
    @unpack nLayer, nWL = rt_dim;

    dnorm  = zeros(FT, nWL);
    piLo   = zeros(FT, nWL);
    piLoc  = zeros(FT, nWL);
    piLos  = zeros(FT, nWL);

    piLoc2 = zeros(FT, (nWL,nLayer));
    ρ_dd   = zeros(FT, (nWL,nLayer));
    ρ_sd   = zeros(FT, (nWL,nLayer));
    τ_dd   = zeros(FT, (nWL,nLayer));
    τ_sd   = zeros(FT, (nWL,nLayer));

    return SWCache{FT}(dnorm  = dnorm ,
                       piLo   = piLo  ,
                       piLoc  = piLoc ,
                       piLos  = piLos ,
                       piLoc2 = piLoc2,
                       ρ_dd   = ρ_dd  ,
                       ρ_sd   = ρ_sd  ,
                       τ_dd   = τ_dd  ,
                       τ_sd   = τ_sd  )
end

###############################################################################
#
# Initialize RTCache
#
###############################################################################
"""
    create_rt_cache(FT, rt_dim::RTDimensions)

Create an [`RTCache`](@ref), given
- `FT` Floating number type
- `rt_dim` [`RTDimensions`](@ref) type struct
"""
function create_rt_cache(FT, rt_dim::RTDimensions)
    cf_con = create_cf_cache(FT, rt_dim);
    cg_con = create_cg_cache(FT, rt_dim);
    sf_con = create_sf_cache(FT, rt_dim);
    sw_con = create_sw_cache(FT, rt_dim);

    return RTCache{FT}(cf_con = cf_con,
                       cg_con = cg_con,
                       sf_con = sf_con,
                       sw_con = sw_con)
end

max_stars_repo_name: bbuman/Land

max_issues_repo_name: bbuman/Land

max_forks_repo_name: bbuman/Land


text:
#!/usr/bin/env julia

using Pkg, SHA, BinaryBuilder, Pkg.BinaryPlatforms

include("common.jl")

if Sys.which("mksquashfs") == nothing
    error("Must install squashfs-tools!")
end

# Parse CLI args
verbose = "--verbose" in ARGS
debug = "--debug" in ARGS

# First, figure out the path to BB and see if we can modify Artifacts.toml
BB_DIR = dirname(dirname(pathof(BinaryBuilder)))
if uperm(joinpath(BB_DIR,"Artifacts.toml")) & 0x02 != 0x02
    error("Unable to write to $(BB_DIR)/Artifacts.toml: dev BinaryBuilder first!")
end

# Let's start by building the Rootfs. We don't want to deploy this one quite yet,
# since it doesn't get a jll package
build_tarballs("Rootfs")

# Next, build the PlatformSupport packages.
for target in supported_platforms()
    build_tarballs("PlatformSupport", triplet(target))
end

# Next, bulid the GCC packages
for gcc_version in ("4.8.5", "5.2.0", "6.1.0", "7.1.0", "8.1.0")
    for target in supported_platforms()
        build_tarballs("GCCBootstrap", "--gcc-version", gcc_version, triplet(target))
    end
end

# LLVM, Rust, etc...
build_tarballs("LLVMBootstrap")
build_tarballs("Rust")
build_tarballs("Go")


text:
# Maybe this should know about t0 as well.
struct SystemTrajectory{h, ΔT, nx, nu, TX<:SizedVector{h,<:SVector{nx}},
                        TU<:SizedVector{h,<:SVector{nu}}}
    "The sequence of states."
    x::TX
    "The sequence of controls."
    u::TU
    "The start time of this trajectory."
    t0::Float64
end

function SystemTrajectory{ΔT}(x::TX, u::TU, t0::Float64) where {h, ΔT, nx, nu,
                                                                TX<:SizedVector{h,<:SVector{nx}},
                                                                TU<:SizedVector{h,<:SVector{nu}}}
    return SystemTrajectory{h, ΔT, nx, nu, TX, TU}(x, u, t0)
end

samplingtime(::SystemTrajectory{h, ΔT}) where {h, ΔT} = ΔT
horizon(::SystemTrajectory{h}) where {h} = h
initialtime(traj::SystemTrajectory) = traj.t0

function time_disc2cont(traj::SystemTrajectory, k::Int)
    @assert 1 <= k <= horizon(traj)
    return initialtime(traj) + (k-1)*samplingtime(traj)
end

function time_cont2disc(traj::SystemTrajectory, t::Float64)
    @assert t >= initialtime(traj)
    k = 1 + round(Int, (t-initialtime(traj))/samplingtime(traj))
    @assert k <= horizon(traj)
    return k
end

timepoints(traj::SystemTrajectory) = (time_disc2cont(traj, k) for k in
                                      1:horizon(traj))

# thin interface with Base for convenience
@inline function Base.zero(::Type{<:SystemTrajectory{h,ΔT,nx,nu}},
                           t0::Float64=0.) where{h,ΔT,nx,nu}
    return SystemTrajectory{ΔT}(zero(SizedVector{h,SVector{nx,Float64}}),
                                zero(SizedVector{h,SVector{nu,Float64}}), t0)
end

Base.zero(traj::SystemTrajectory) = zero(typeof(traj), initialtime(traj))

function Base.zero(::Type{SystemTrajectory}, g)
    h = horizon(g)
    ΔT = samplingtime(g)
    nx = n_states(g)
    nu = n_controls(g)
    return zero(SystemTrajectory{h,ΔT,nx,nu})
end

function zero!(traj::SystemTrajectory)
    for k in 1:horizon(traj)
        traj.x[k] = zero(traj.x[k])
        traj.u[k] = zero(traj.u[k])
    end
    return traj
end

@inline function Base.copy(traj::SystemTrajectory)
    return SystemTrajectory{samplingtime(traj)}(copy(traj.x), copy(traj.u),
                                                initialtime(traj))
end

@inline @inbounds function Base.copyto!(dest::SystemTrajectory, src::SystemTrajectory)
    @assert initialtime(dest) == initialtime(src)
    @assert horizon(dest) == horizon(src)
    copyto!(dest.x, src.x)
    copyto!(dest.u, src.u)
    return dest
end


text:
using BenchmarkTools, CartesianGrids

BenchmarkTools.DEFAULT_PARAMETERS.gcsample = true

SUITE = BenchmarkGroup()
SUITE["field ops"] = BenchmarkGroup()
SUITE["ddf"] = BenchmarkGroup()
SUITE["regularization"] = BenchmarkGroup()

include("setupproblem.jl")

w .= randn(size(w));
f .= randn(size(f));

SUITE["field ops"]["inverse Laplacian"] = @benchmarkable L\w;
SUITE["field ops"]["curl of dual nodes"] = @benchmarkable curl!(q,w);
SUITE["field ops"]["divergence of primal edges"] = @benchmarkable divergence!(p,q);

for kernel in (:Roma,:Yang3,:M4prime)
  name = string(kernel)
  @eval SUITE["ddf"][$name] = @benchmarkable evaluate_ddf(DDF(ddftype=CartesianGrids.$kernel))
end

SUITE["regularization"]["create node regularization matrix"] = @benchmarkable RegularizationMatrix(reg,f,w);
SUITE["regularization"]["create node interpolation matrix"] = @benchmarkable InterpolationMatrix(reg,w,f);

Rn = RegularizationMatrix(reg,f,w);
En = InterpolationMatrix(reg,w,f);
SUITE["regularization"]["evaluate node regularization"] = @benchmarkable Rn*f;
SUITE["regularization"]["evaluate node interpolation"] = @benchmarkable En*w;

SUITE["regularization"]["create edge regularization matrix"] = @benchmarkable RegularizationMatrix(reg,τ,q);
SUITE["regularization"]["create edge interpolation matrix"] = @benchmarkable InterpolationMatrix(reg,q,τ);

Re = RegularizationMatrix(reg,τ,q);
Ee = InterpolationMatrix(reg,q,τ);
SUITE["regularization"]["evaluate edge regularization"] = @benchmarkable Re*τ;
SUITE["regularization"]["evaluate edge interpolation"] = @benchmarkable Ee*q;

run(SUITE)

max_stars_repo_name: JuliaIBPM/CartesianGrids

max_issues_repo_name: JuliaIBPM/CartesianGrids

max_forks_repo_name: JuliaIBPM/CartesianGrids


text:
using RandomMatrices
using Base.Test

#Test far outside support
#Tracy-Widom has support on all x>0, but the integration won't pick it up
@test pdf(TracyWidom, -10) == pdf(TracyWidom, 10) == 0
@test cdf(TracyWidom, -10) == 0
@test cdf(TracyWidom, 10) == 1

if isdefined(:ODE) && isa(ODE, Module)
    t = rand()
    @test pdf(TracyWidom, t) > 0
    @test 0 < cdf(TracyWidom, t) < 1
end

@test isfinite(rand(TracyWidom, 10))
@test isfinite(rand(TracyWidom, 100))


text:
# Autogenerated wrapper script for gb_jll for x86_64-apple-darwin
export libgb

using GMP_jll
using CompilerSupportLibraries_jll
JLLWrappers.@generate_wrapper_header("gb")
JLLWrappers.@declare_library_product(libgb, "@rpath/libgb.0.dylib")
function __init__()
    JLLWrappers.@generate_init_header(GMP_jll, CompilerSupportLibraries_jll)
    JLLWrappers.@init_library_product(
        libgb,
        "lib/libgb.0.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
struct Scaler
    channel::Int8
    scaler::Int32
    Scaler() = new(0, 0)
    Scaler(x, y) = new(x, y)
end

mutable struct ADCTDC
    channel::Int8
    adc::Int32
    tdc::Int32
    ADCTDC() = new(0, 0, 0)
    ADCTDC(x, y, z) = new(x, y, z)
end

mutable struct ADC
    channel::Int8
    data::Int32
    associated_channel::Int8
    ADC() = new(0, 0, 0)
end

mutable struct Event
    # Channel scaler - Not used
    # scaler::Scaler
    #
    # The SiRi front detector ADC values
    enum::Int8
    e::Vector{ADC}

    # The SiRi back detector ADC values
    Δenum::Int8
    Δe::Vector{ADC}

    # OSCAR TDC and ADC values
    # Uses a poor man's dict to avoid Julia's costly
    # hashing. labrnum keeps the length of the
    # labrindices, and labrindices keeps the indices
    # /channels used in labr. Labr has the actual channels
    # and data
    labrnum::Int8 # Keeps track of how many labr are used
    labrindices::Vector{Int8}
    labr::Vector{ADCTDC}
    # labr::Dict{UInt32, ADCTDC}

    # Wall clock time
    # time::Union{Int32, Nothing}
    # The TPU pattern
    # pattern::Union{Int32, Nothing}
    function Event()
        e = [ADC() for _ in 1:16]
        Δe = [ADC() for _ in 1:16]
        labr = [ADCTDC() for _ in 1:32]
        labrind = zeros(Int8, 32)
        new(0, e, 0, Δe, 0, labrind, labr)
    end
end

function Base.show(io::IO, e::Event)
    # println("LaBr: $(e.labr)")
    println("# of back: ", e.enum)
    println("# of front: ", e.Δenum)
    println("# of LaBr: ", e.labrnum)
    println("SiRi back(e):")
    [println(adc) for adc in e.e]
    println("\nSiRi front(Δe):")
    [println(adc) for adc in e.Δe]
    println("\nLaBr: ")
    for i in 1:e.labrnum
        chn = e.labrindices[i]
        println("chn: ", chn, " = ", e.labr[chn])
    end
    # println("\nWall clock time: $(repr(e.time))")
    # println("TPU Pattern: $(repr(e.pattern))")
end

function Base.show(io::IO, adc::ADC)
    print("\tchl: $(adc.channel), adc: $(adc.data)")
end

function Base.show(io::IO, adc::ADCTDC)
    print("\tadc: $(adc.adc), tdc: $(adc.tdc)")
end

function setadc!(event::Event, channel, adc)::Nothing
    event.labr[channel+1].adc = adc
    found = false
    for i ∈ 1:event.labrnum
        if event.labrindices[i] == channel+1
            found = true
            break
        end
    end
    if !found
        event.labrnum += 1
        event.labrindices[event.labrnum] = channel+1
        event.labr[channel+1].channel = channel
    end
    nothing
end

function settdc!(event::Event, channel, tdc)::Nothing
    event.labr[channel+1].tdc = tdc
    found = false
    for i ∈ 1:event.labrnum
        if event.labrindices[i] == channel+1
            found = true
            break
        end
    end
    if !found
        event.labrnum += 1
        event.labrindices[event.labrnum] = channel+1
        event.labr[channel+1].channel = channel
    end
    nothing
end

function adde!(event, channel, data)
    event.enum += 1 # First to account for 1-indexing
    event.e[event.enum].channel = channel >> 1
    event.e[event.enum].data = data
end

function addΔe!(event, channel, data)
    # HACK This is a bug. Sometimes it gets greater than 17 and crashes
    if event.Δenum == 16
        return
    end
    event.Δenum += 1
    event.Δe[event.Δenum].channel = channel
    event.Δe[event.Δenum].associated_channel = div(channel, 8)
    event.Δe[event.Δenum].data = data
end

function reset!(e::Event)
    #e.labr = Dict{UInt32, ADCTDC}()
    # Can skip this step, as the
    # counters decide what is being read
    # for i in 1:16
        # e.e[i].data = e.e[i].channel = 0
        # e.Δe[i].data = e.Δe[i].channel = 0
    # end
    e.enum = 0
    e.Δenum = 0
    for channel in e.labrindices[1:e.labrnum]
        e.labr[channel].adc = 0
        e.labr[channel].tdc = 0
    end
    e.labrnum = 0
end

function eachlabr(e::Event)
    # Returns channel and labr in that channel
    return [(e.labrindices[i], e.labr[e.labrindices[i]]) for i in 1:e.labrnum]
end

function islabrgood(adctdc::ADCTDC)
    return adctdc.adc > 0 &&  adctdc.tdc > 0
end

struct MiniEvent
    front::Int8
    back::Int8
    frontchannel::Int8
    backchannel::Int8
    Δe::Float32
    e::Float32
end

function Base.show(io::IO, event::MiniEvent)
    println("b$(event.back)f$(event.front): Δe = $(event.Δe), e = $(event.e)")
end

struct LaBrEvent
    mevent::MiniEvent
    labr::Vector{ADCTDC}
    function LaBrEvent(event::Event, mevent::MiniEvent)
        new(mevent, filter(x -> !islabrgood(x), event.labr[event.labrindices[1:event.labrnum]]))
    end
    function LaBrEvent(mevent::MiniEvent)
        new(mevent, ADCTDC[])
    end
end

function serialize(event::LaBrEvent)
    return [serialize(event.mevent)..., convert(UInt32, length(event.labr))], [serialize(labr) for labr in event.labr]
end

function serialize(event::MiniEvent)
    return (event.front, event.back, event.Δe, event.e) .|> x -> convert(UInt32, x)
end

function serialize(adc::ADCTDC)
    return (adc.channel, adc.adc, adc.tdc) .|> x -> convert(UInt32, x)
end

struct TrueEvent
    f::Int8
    b::Int8
    Δe::Float32
    e::Float32
    labr::Vector{ADCTDC}
    function TrueEvent(f, b, de, e, labr)
        new(f, b, de, e, labr)
    end
end

function save(labrs::AbstractArray{LaBrEvent}, path::AbstractString)
    # Terrible performance. Concoct some smarter way to save data
    !isdir(path) && mkdir(path)
    fhandle = open(joinpath(path, "labr.bin"), "w")
    for labr in labrs
        header, body = serialize(labr)
        write(fhandle, header)
        for b in body
            write(fhandle, b...)
        end
    end
    close(fhandle)
end

max_stars_repo_name: Caronthir/JSort

max_issues_repo_name: Caronthir/JSort

max_forks_repo_name: Caronthir/JSort


text:
let
    gw = LegacyGridWorld(sx=2, sy=2, rs=[GridWorldState(1,1)], rv=[10.0])

    pvec = fill(GridWorldAction(:left), 5)

    solver = VectorSolver(pvec)

    p = solve(solver, gw)

    io = IOBuffer()
    d = TextDisplay(io)
    display(d, p)
    @test length(split(String(take!(io)), '\n')) == 6

    for s1 in states(gw)
        @test action(p, s1) == GridWorldAction(:left)
    end

    p2 = VectorPolicy(gw, pvec)
    for s2 in states(gw)
        @test action(p2, s2) == GridWorldAction(:left)
    end

    p3 = ValuePolicy(gw)
    for s2 in states(gw)
        @inferred(action(p3, s2)) isa GridWorldAction
    end

    io = IOBuffer()
    d = TextDisplay(io)
    display(d, p3)
    @test length(split(String(take!(io)), '\n')) == 6
end

max_stars_repo_name: JuliaPOMDP/POMDPPolicies

max_issues_repo_name: JuliaPOMDP/POMDPPolicies

max_forks_repo_name: JuliaPOMDP/POMDPPolicies


text:
macro timeit(f)
      start = time()
      out = eval(f)
      elapsed = time() - start
      print("Time elapsed: ", elapsed, "\n")
      (output = out, time = elapsed)
end

function std_aux_laws(AuxLaw, param, θ, obs, obs_time, idx=1:length(obs[1]))
    if param === nothing # a temporary solution before redefining FHN
        laws = [AuxLaw(θ..., t₀, u[idx], T, v[idx]) for (t₀,T,u,v)
                in zip(obs_time[1:end-1], obs_time[2:end], obs[1:end-1], obs[2:end])]
    else
        laws = [AuxLaw(param, θ..., t₀, u[idx], T, v[idx]) for (t₀,T,u,v)
                in zip(obs_time[1:end-1], obs_time[2:end], obs[1:end-1], obs[2:end])]
    end
    display(laws[1])
    laws
end


text:
mutable struct BaseData{T <: Unsigned}
    id::Int64
    input::T
    output::T

    function BaseData{T}() where {T <: Unsigned}
        new(0, 0, 0)
    end
end


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md

using JuliaFEM, Test

# Test stiffness matrix of geometrically nonlinear problem

X = Dict(1 => [0.0, 0.0], 2 => [1.0, 0.0], 3 => [1.0, 1.0], 4 => [0.0, 1.0])
u = Dict(1 => [0.1, 0.2], 2 => [0.3, 0.4], 3 => [0.5, 0.6], 4 => [0.7, 0.8])
T = Dict(3 => [0.0, 288.0], 4 => [0.0, 288.0])

element = Element(Quad4, (1, 2, 3, 4))
update!(element, "geometry", X)
update!(element, "displacement", u)
update!(element, "youngs modulus", 288.0)
update!(element, "poissons ratio", 1/3)

traction = Element(Seg2, (3, 4))
update!(traction, "geometry", X)
update!(traction, "displacement", u)
update!(traction, "displacement traction force", T)

# field problem
block = Problem(Elasticity, "block", 2)
block.properties.formulation = :plane_stress
block.properties.finite_strain = true
block.properties.geometric_stiffness = true
add_element!(block, element)
#add_elements!(block, traction)

assemble!(block, 0.0)
Km = Matrix(block.assembly.K)
Kg = Matrix(block.assembly.Kg)
K = Km + Kg
f = Vector(block.assembly.f)

K_expected = [
  401.76   200.88  -123.12    -5.76  -191.52  -117.36  -87.12   -77.76
  200.88   473.76    -5.76    28.08  -117.36  -205.92  -77.76  -295.92
 -123.12    -5.76   197.28    -2.16   -12.24   -25.92  -61.92    33.84
   -5.76    28.08    -2.16   298.08   -25.92  -163.44   33.84  -162.72
 -191.52  -117.36   -12.24   -25.92   240.48   120.24  -36.72    23.04
 -117.36  -205.92   -25.92  -163.44   120.24   312.48   23.04    56.88
  -87.12   -77.76   -61.92    33.84   -36.72    23.04  185.76    20.88
  -77.76  -295.92    33.84  -162.72    23.04    56.88   20.88   401.76]
#   f_expected = [142.272, 214.272, -13.824, 58.176, -75.456, 19.584, -52.992, -4.032]
f_expected = [142.272, 214.272, -13.824, 58.176, -75.456, -124.416, -52.992, -148.032]
@test isapprox(K, K_expected)
@test isapprox(f, f_expected)


text:
############################################################################################
"""
$(TYPEDEF)
Markov Kernel for particle propagation.

# Fields
$(TYPEDFIELDS)
"""
struct Markov{A,B,C} <: ParticleKernel
    "Initial distribution, function of iter only."
    initial::A
    "Transition distribution, function of full particle trajectory and current iteration count."
    transition::B
    "Data distribution to weight particles. Function of full data, particle trajectory and current iteration count."
    evidence::C
    function Markov(initial::A, transition::B, evidence::C) where {A,B,C}
        return new{A,B,C}(initial, transition, evidence)
    end
end

############################################################################################
function initial(_rng::Random.AbstractRNG, kernel::Markov)
    return rand(_rng, kernel.initial)
end

############################################################################################
function transition(
    _rng::Random.AbstractRNG, kernel::Markov, val::AbstractArray{P}, iter::Integer
) where {P}
    return rand(_rng, kernel.transition(val, iter))
end

############################################################################################
function ℓtransition(
    valₜ::Union{P,AbstractArray{P}}, kernel::Markov, val::AbstractArray{P}, iter::Integer
) where {P}
    return logpdf(kernel.transition(val, iter), valₜ)
end

############################################################################################
export Markov


text:
## ------------------------------------------------------------------
function _new_section!(callinfo::LineNumberNode, label)

    # book
    doc = _check_currdoc()

    # TODO: Update book but just if necessary

    # checks
    _check_call_inbook(callinfo, doc)
    _check_bang_call(callinfo, doc)

    # new section
    sec = RBSection(label)
    setdoc!(sec, doc)
    srcfile!(sec, srcfile(callinfo))
    srcline!(sec, srcline(callinfo))

    # add to doc
    add_sec!(doc, sec)

    currsec!(sec)
end

## ------------------------------------------------------------------
macro new_section!(ex...)
    label = _insert_label!("new_section!", ex, __source__)
    _new_section!(__source__, label)
end


text:
using Knet
using ReinforcementLearning
using ArgParse
using JLD

const GAMMA = 0.95
const WINIT = 0.1
const EPISODES = 1000

function main(args)
    s = ArgParseSettings()
    s.description = ""

    @add_arg_table s begin
        ("--actor"; nargs='*'; arg_type=Int64)
        ("--critic"; nargs='*'; arg_type=Int64)
        ("--atype"; default=(gpu()>=0 ? "KnetArray{Float32}":"Array{Float32}"))
        ("--optim"; default="Sgd(;lr=0.001)")
        ("--discount"; default=GAMMA; arg_type=Float64)
        ("--winit"; default=WINIT; arg_type=Float64)
        ("--loadfile"; default=nothing)
        ("--episodes"; default=EPISODES; arg_type=Int64)
        ("--seed"; default=-1; arg_type=Int64)
        ("--monitor"; default=nothing)
        ("--generate"; default=0; arg_type=Int64)
        ("--period"; default=100; arg_type=Int64)
        ("--pole"; default=1; arg_type=Int64)
    end

    isa(args, AbstractString) && (args=split(args))
    o = parse_args(args, s; as_symbols=true); println(o); flush(STDOUT)
    o[:atype] = eval(parse(o[:atype]))
    sr = o[:seed] > 0 ? srand(o[:seed]) : srand()

    pole = o[:pole]
    threshold = pole == 0 ? 200 : 500
    env = GymEnv("CartPole-v$pole")
    w = θ = opts = nothing
    if o[:loadfile] == nothing
        θ = initweights(o[:atype], o[:actor])
        w = initweights(o[:atype], o[:critic])
    else
        w = load(o[:loadfile], "w")
        w = map(wi->convert(o[:atype]), w)
        θ = load(o[:loadfile], "w")
        θ = map(θi->convert(o[:atype]), θ)
    end
    aopt = map(θi->eval(parse(o[:optim])), θ)
    copt = map(wi->eval(parse(o[:optim])), w)

    if o[:monitor] != nothing
        monitor_start(env, o[:monitor])
    end

    best_episode = -Inf
    avgreward = 0
    for k = 1:o[:episodes]
        s0 = getInitialState(env); state = s0;
        episode_rewards = 0
        episode_length  = 0
        history = []

        while !state.done
            s = convert(o[:atype], state.data)
            probs = softmax(predict(θ,s))

            # this is a, action
            a = action = sample_action(probs)[1]

            # this is transition sampled
            next_state, reward = transfer(env, state, env.actions[action])
            push!(history, Transition(state, action, next_state, reward))

            # compute TD target
            sp = history[end].state.data
            ap = sample_action(softmax(predict(θ,s)))[1]
            td = reward + o[:discount] * predict(w,sp)[ap]

            if next_state.done && episode_length < threshold
                td -= threshold
            end

            values = []
            s = convert(o[:atype], s)
            ∇w = criticgrad(w, s, a, td; values=values)
            qw = values[1]

            # train actor
            ∇θ = actorgrad(θ, s, a, qw)
            update!(θ, ∇θ, aopt)

            # train critic
            update!(w, ∇w, copt)

            state = next_state
            episode_length += 1
            episode_rewards += reward
        end

        # average_reward = episode_rewards / episode_length
        if k < 100
            avgreward = avgreward * (k-1) + episode_rewards
            avgreward = avgreward / k
        else
            avgreward = 0.01 * episode_rewards + avgreward * 0.99
        end

        if k % o[:period] == 0
            println("(episode:$k,avgreward:$avgreward)")
        end

        if avgreward > best_episode
            best_episode = avgreward
        end

        empty!(history)
    end

    for k = 1:o[:generate]
        done = false
        state = getInitialState(env)
        rewards = 0
        while !done
            s = convert(o[:atype], state.data)
            pred = predict(θ,s)
            probs = softmax(pred)
            # action = sample_action(probs)[1]
            action = select_action(probs)
            render_env(env)
            next_state, reward = transfer(env,state,env.actions[action])
            state = next_state
            done = next_state.done
            rewards += reward
        end
        println("(generation:$k,reward:$rewards")
    end
end

type Transition
    state
    action
    next_state
    reward
end

function initweights(atype, hidden, xsize=4, ysize=2, winit=0.1)
    w = []
    x = xsize
    for y in [hidden..., 2]
        push!(w, xavier(y,x))
        push!(w, zeros(y,1))
        x = y
    end
    return map(wi->convert(atype,wi), w)
end

function predict(w,x)
    for i = 1:2:length(w)
        x = w[i] * x .+ w[i+1]
        if i < length(w)-1
            x = relu(x)
        end
    end
    return x
end

function softmax(x0)
    x1 = maximum(x0,1)
    x2 = x0 .- x1
    x3 = exp(x2)
    x4 = x3 ./ sum(x3,1)
end

function actor(w, x, actions, values)
    ypred = predict(w,x)
    return -sum(logprob(actions,ypred) .* values) / size(x,2)
end

actorgrad = grad(actor)

function critic(w,s,a,ret; values=[])
    qsa = predict(w,s)[a]
    val = sumabs2(ret-qsa)
    push!(values, val)
    return val
end

criticgrad = grad(critic)

function sample_action(probs)
    probs = Array(probs)

    cprobs = cumsum(probs)
    sampled = cprobs .> rand()
    actions = mapslices(indmax, sampled, 1)

    # an alternative?
    # sampled = probs .* rand(size(probs))
    # actions = mapslices(indmax, sampled, 1)

    return actions
end

function select_action(probs)
    return indmax(probs)
end

function logprob(output, ypred)
    nrows,ncols = size(ypred)
    index = output + nrows*(0:(length(output)-1))
    o1 = logp(ypred,1)
    o2 = o1[index]
    return o2
end

!isinteractive() && !isdefined(Core.Main, :load_only) && main(ARGS)

max_stars_repo_name: ilkerkesen/rl

max_issues_repo_name: ilkerkesen/rl

max_forks_repo_name: ilkerkesen/rl


text:
"""
    PointNode(coordinate, properties)

## Examples
```jldoctest
julia> using NeuroCore

julia> PointNode((1,2))
```
"""
struct PointNode{Dim,T,C<:AbstractPoint{Dim,T},P} <: AbstractPoint{Dim,T}
    coordinates::C
    properties::P

    PointNode{Dim,T,C,P}(coord::C, props::P) where {Dim,T,C,P} = new{Dim,T,C,P}(coord, props)
end

function PointNode{Dim,T,C}(coord::C, props::P=nothing) where {Dim,T,C<:AbstractPoint{Dim,T},P}
    return PointNode{Dim,T,C,P}(coord, props)
end

function PointNode{Dim,T}(coord::AbstractPoint{Dim,T}, props=nothing) where {Dim,T}
    return PointNode{Dim,T,typeof(coord)}(coord, props)
end

function PointNode{Dim}(coord::AbstractPoint{Dim,T}, props=nothing) where {Dim,T}
    return PointNode{Dim,T}(coord, props)
end

function PointNode(coord::AbstractPoint{Dim}, prop=nothing) where {Dim}
    return PointNode{Dim}(coord, prop)
end

PointNode(coord::Tuple) = PointNode(Point(coord))

GeometryBasics.coordinates(n::PointNode) = getfield(n, :coordinates)

FieldProperties.properties(n::PointNode) = getfield(n, :properties)

#FieldProperties.@hasfieldprop PointNode

Base.@propagate_inbounds Base.getindex(n::PointNode, i::Int) = coordinates(n)[i]

Base.show(io::IO, n::PointNode) = print_node(io, n)

Base.show(io::IO, ::MIME"text/plain", n::PointNode) = print_node(io, n)

function print_node(io::IO, n::PointNode{Dim,T}) where {Dim,T}
    return print(io, "PointNode($(coordinates(n).data), $(properties(n)))")
end

"""
    PointNodeList{Dim,T,V,P}

"""
const PointNodeList{Dim,T<:PointNode{Dim},V<:AbstractVector{T},P} = MetadataArray{T,1,P,V}

#FieldProperties.@hasfieldprop PointNodeList

function PointNodeList{Dim,T}(v::AbstractVector{T}=T[], props=nothing) where {Dim,T<:PointNode{Dim}}
    return MetadataArray(v, props)
end

PointNodeList(v::AbstractVector{T}, props=nothing) where {T<:PointNode} = PointNodeList{length(T),T}(v, props)

"""
    MultiPointNodeList{T,V,P}
"""
const MultiPointNodeList{T<:PointNodeList,V<:AbstractVector{T},P} = MetadataArray{T,1,P,V}

#FieldProperties.@hasfieldprop PointNodeListVec


text:
using Statistics
include("CKlib.jl")
function Analyze(chain)
    lb,ub = PriorSupport()
    θtrue = TrueParameters()
    nParams = size(lb,1)
    inci01 = zeros(nParams)
    inci05 = zeros(nParams)
    inci10 = zeros(nParams)
    lower = zeros(nParams)
    upper = zeros(nParams)
    for i = 1:nParams
        lower[i] = quantile(chain[:,i],0.005)
        upper[i] = quantile(chain[:,i],0.995)
        inci01[i] = θtrue[i] >= lower[i] && θtrue[i] <= upper[i]
        lower[i] = quantile(chain[:,i],0.025)
        upper[i] = quantile(chain[:,i],0.975)
        inci05[i] = θtrue[i] >= lower[i] && θtrue[i] <= upper[i]
        lower[i] = quantile(chain[:,i],0.05)
        upper[i] = quantile(chain[:,i],0.95)
        inci10[i] = θtrue[i] >= lower[i] && θtrue[i] <= upper[i]
    end
    return vcat(mean(chain,dims=1)[:], inci01[:], inci05[:], inci10[:])
end

max_stars_repo_name: mcreel/ILSNM

max_issues_repo_name: mcreel/ILSNM

max_forks_repo_name: mcreel/ILSNM


text:
peek{T}(v::Vector{T}) = v[end]

function shove!{T}(a::Vector{T}, i::Integer)
  ln = length(a)
  i = min(max(1, ln - i), ln + 1)
  pop!(insert!(a, i, a[end]))
  return a
end

function yank!{T}(a::Vector{T}, i::Integer)
  ln = length(a)
  i = min(max(1, ln - i), ln)
  push!(a, splice!(a, i))
end

function yankdup!{T}(a::Vector{T}, i::Integer)
  ln = length(a)
  i = min(max(1, ln - i), ln)
  push!(a, a[i])
end


text:
function collect_eachmatch(re, str, overlap=false)
    [m.match for m in collect(eachmatch(re, str, overlap))]
end

for f in [matchall, collect_eachmatch]
    @test f(r"a?b?", "asbd") == ["a","","b","",""] == f(r"""a?b?""", "asbd")
    @test f(r"a?b?", "asbd", true) == ["a","","b","",""]
    @test f(r"\w+", "hello", true) == ["hello","ello","llo","lo","o"]
    @test f(r".\s", "x \u2200 x \u2203 y") == ["x ", "∀ ", "x ", "∃ "]
    @test f(r"(\w+)(\s*)", "The dark side of the moon") ==
          ["The ", "dark ", "side ", "of ", "the ", "moon"]
    @test f(r"", "") == [""]
    @test f(r"", "", true) == [""]
    @test f(r"aa", "aaaa") == ["aa", "aa"]
    @test f(r"aa", "aaaa", true) == ["aa", "aa", "aa"]
    @test f(r"", "aaa") == ["", "", "", ""]
    @test f(r"", "aaa", true) == ["", "", "", ""]
end

max_stars_repo_name: jakebolewski/julia

max_issues_repo_name: jakebolewski/julia

max_forks_repo_name: jakebolewski/julia


text:
using MyPkg
using Test

@testset "MyPkg.jl" begin
    # Write your own tests here.
end

max_stars_repo_name: jannittayao/MyPkg

max_issues_repo_name: jannittayao/MyPkg

max_forks_repo_name: jannittayao/MyPkg


text:
@testset "Test quadratures" begin

    using jFEMTools
    jF = jFEMTools

    # Triangle Quadratures
    ###############################################

    # Test Strang Quadrature
    for i in 1:5
        quad_rule = jF.QuadratureRule{jF.Triangle}(jF.Strang(),i)
        @test sum(quad_rule.weights) ≈ 0.5
    end

end

max_stars_repo_name: Paulms/jFEMToools

max_issues_repo_name: Paulms/jFEMToools

max_forks_repo_name: Paulms/jFEMToools


text:
# # Using IterationControl.jl to train a tree-booster on the iris data set

# In this demonstration we show how to the controls in
# [IterationControl.jl](https://github.com/ablaom/IterationControl.jl)
# with an iterative
# [MLJ](https://github.com/alan-turing-institute/MLJ.jl) model, using
# our bare hands. (MLJ will ultimately provide its own canned
# `IteratedModel` wrapper to make this more convenient and
# compositional.)

using Pkg
Pkg.activate(@__DIR__)
Pkg.instantiate()

import MLJ
using IterationControl

using Statistics
using Random
Random.seed!(123)

MLJ.color_off()

# Loading some data and splitting observation indices into test/train:

X, y = MLJ.@load_iris;
train, test = MLJ.partition(eachindex(y), 0.7, shuffle=true)

# Import an model type:

Booster = MLJ.@load EvoTreeClassifier verbosity=0

# Note that in MLJ a "model" is just a container for
# hyper-parameters. The objects we will iterate here are MLJ
# [*machines*](https://alan-turing-institute.github.io/MLJ.jl/dev/machines/);
# these bind the model to train/test data and, in the case of
# iterative models, can be trained using a warm-restart.

# Creating a machine:

mach = MLJ.machine(Booster(nrounds=1), X, y);

# Lifting MLJ's `fit!(::Machine)` method to `IterativeControl.train!`:

function IterationControl.train!(mach::MLJ.Machine{<:Booster}, n::Int)
    mlj_model = mach.model
    mlj_model.nrounds = mlj_model.nrounds + n
    MLJ.fit!(mach, rows=train, verbosity=0)
end

# Lifting the out-of-sample loss:

function IterationControl.loss(mach::MLJ.Machine{<:Booster})
    mlj_model = mach.model
    yhat = MLJ.predict(mach, rows=test)
    return MLJ.log_loss(yhat, y[test]) |> mean
end

# Iterating with controls:

logging(mach) = "loss: $(IterationControl.loss(mach))"

IterationControl.train!(mach,
                        Step(5),
                        GL(),
                        Info(logging))

# Continuing iteration with a different stopping criterion:

IterationControl.train!(mach,
                        Step(5),
                        NumberLimit(10),
                        Info(logging))

using Literate #src
Literate.markdown(@__FILE__, @__DIR__, execute=false) #src
Literate.notebook(@__FILE__, @__DIR__, execute=true) #src


text:
include("engine.jl")

##############
#AI
##############

#DEBUG
global abN = 0
global gLast=nothing
#END DEBUG

moves = Dict{Card, Int64}()

#egymas melletti lapok (pl. tF,tU,t9) eleg az egyiket valasztanunk
#mondjuk a legnagyobbat.
#Kivetel1: zT,zK nem mindegy. Kivetel2: adu 7-es vagy 8-as nem mindegy (ulti-repulonel)

#a legnagyobbat valasztjuk - lehetne a legkisebb vagy random is
pickFromStreak(streak::CardSet32) = length(streak) > 1 ? first(streak) : streak
adu8(g::GameState) = Card(g.contract.suit, _8)
adu7(g::GameState) = Card(g.contract.suit, _7)

extendStreak(result, streak, card) = (result, union(streak, card))
breakStreak(result, streak, card, inV) = (union(result, pickFromStreak(streak)), inV ? card : CardSet32())

#BUG TODO: check 494c2399f96cc0af5508336cab0f7415b5708f44 vs. d51d5d75f51ee3a1ada61c989eccfaa7df8bc0e3 2x time and abN
function removeNeighbours(g::GameState, vm::CardSet32)
    if length(vm) <= 1 return vm end

    (result, streak) = (CardSet32(), CardSet32())

    remaining = union(ps(g, p1).hand, ps(g, p2).hand, ps(g, p3).hand)
    @assert (vm in remaining)
    # remaining = setdiff(s1, s2)(remaining, vm)

    for card in remaining
        #DEBUG show(result); print(" | "); show(streak); print(" | "); show(card);println()
        if card in vm
            if !(card in AT) && hanyAT(streak) > 1 || #tizes es nemtizes
                    !sameSuit(card, streak, suitof(card)) || #mas szin
                    card == adu7(g)  #adu7 tori, mert 1. adu 8-7 - mindketto fontos 2. csendesultinal szamit | Note: 8-as nem, mert ugysem fordul elo, hogy 8-9 teheto

                (result, streak) = breakStreak(result, streak, card, true)
            else
                (result, streak) = extendStreak(result, streak, card)
            end
        else #in remaining
            (result, streak) = breakStreak(result, streak, card, false)
        end
    end
    #DEBUG show(result); print(" | "); show(streak); print(" | ");println()

    (result, streak) = breakStreak(result, streak, CardSet32(), false)  #register last streak
    #DEBUG show(result); print(" | "); show(streak); print(" | ");println()

    return result
end

#     #terrible spaghetti - but a bit faster :(
# function removeNeighbours(g::GameState, vm::CardSet32)
#     #terrible spaghetti - maybe rewrite one day?
#     if length(vm) <= 1 return vm end

#     remaining = union(ps(g, p1).hand, ps(g, p2).hand, ps(g, p3).hand)
#     @assert (length(remaining) > length(vm))

#     streak = CardSet32()
#     result = CardSet32()

#     stateV = start(vm)
#     stateR = start(remaining)
#     (cV, stateV) = next(vm, stateV)
#     (cR, stateR) = next(remaining, stateR)
#     oneMoreRound = 2 #to make sure the last is recorded
#     while !done(vm, stateV) || !done(vm, stateR) || oneMoreRound > 0
#         #DEBUG: show(streak);print(" | "); show(result); print(" | "); show(cV); show(cR); println()
#         if cR == cV
#             #TK boundary, suit boundary, or adu8,7
#             if (cV in K && !isempty(streak)) || !sameSuit(cV, streak, suitof(cV)) || (cV == adu7(g) && adu8(g) in streak)
#                 result = union(result, pickFromStreak(streak)) #a legnagyobbat valasztjuk - lehetne a legkisebb vagy random is
#                 streak = cV
#             end
#             streak = union(streak, cV) #continue streak

#             #move both
#             if !done(vm, stateV)
#                 (cV, stateV) = next(vm, stateV)
#             end
#             if !done(vm, stateR)
#                 (cR, stateR) = next(remaining, stateR)
#             end
#         else
#             if !isempty(streak)
#                 result = union(result, pickFromStreak(streak)) #a legnagyobbat valasztjuk - lehetne a legkisebb vagy random is
#                 streak = CardSet32()
#             end
#             if !done(vm, stateR)
#                 (cR, stateR) = next(remaining, stateR)
#             end
#         end

#         if done(vm, stateV)
#             oneMoreRound -= 1
#         end
#     end

#     if !isempty(streak)
#         result = union(result, pickFromStreak(streak))
#     end

#     return result
# end

#heuristics of ordering the right moves
#TODO smallestCards
function sortMoves(g::GameState, vm::CardSet32)
    if length(vm) <= 1 return vm end
    valid = vm

    result = Vector{Card}() #ordered

    #sort suit by length
    shorterThan(x,y) = length(valid * x)>length(valid * y)
    suitsInvalid = filter(x->!isempty(valid * x), [t,z,m,p]) #TODO: randomize maybe later?
    suitsInLengthOrder = sort(suitsInvalid, lt=shorterThan)

    function largestCardsInsuitLengthOrder(result, valid, what=:largest, rev::Bool=false)
        suitList = rev ? suitsInLengthOrder : reverse(suitsInLengthOrder)
        for suit in suitList
            if isempty(valid * suit)
                continue
            end

            if what == :largest
                cards = largestCard(valid * suit, suit)
                push!(result, cards)
            elseif what == :smallest
                # cards = last(valid * suit) #TODO: properly for betli/duri
                cards = smallestCard(valid * suit, suit)
                push!(result, cards)
            else #:all
                cards = valid * suit
                append!(result, toArray(cards))
            end
            valid -= cards
        end
        #show(result); print(" | "); show(valid); println() #DEBUG
        return (result, valid)
    end

    smallestCardsInsuitLengthOrder(result, valid) = largestCardsInsuitLengthOrder(result, valid, :smallest)
    smallestCardsInReverseSuitLengthOrder(result, valid) = largestCardsInsuitLengthOrder(result, valid, :smallest, true)
    allCardsInsuitLengthOrder(result, valid) = largestCardsInsuitLengthOrder(result, valid, :all)
    allCardsInReverseSuitLengthOrder(result, valid) = largestCardsInsuitLengthOrder(result, valid, :all)
    largestCardsInReverseSuitLengthOrder(result, valid) = largestCardsInsuitLengthOrder(result, valid, :largest, true)

    if g.contract.suit <= p #adu jatek
        if length(g.talon) < 2 # talonozas
            (result, valid) = smallestCardsInReverseSuitLengthOrder(result, valid)
            (result, valid) = allCardsInReverseSuitLengthOrder(result, valid)
        elseif isempty(g.asztal) #elso helyrol
            #largest cards in suit length order
            (result, valid) = largestCardsInsuitLengthOrder(result, valid)
            #smallest cards in suit length order
            (result, valid) = smallestCardsInsuitLengthOrder(result, valid)
            #the rest in suit length order
            (result, valid) = allCardsInsuitLengthOrder(result, valid)
        else # masodik vagy harmadik helyrol - legkisebb dobas vagy felulutes
            #smallest card in suit length order
            (result, valid) = smallestCardsInsuitLengthOrder(result, valid)
            #then largest card in suit length order
            (result, valid) = largestCardsInsuitLengthOrder(result, valid)
            #the rest in suit length order
            (result, valid) = allCardsInsuitLengthOrder(result, valid)
        end
    else #TODO szintelen jatek, foleg betli

    end

    if !isempty(valid) push!(result, toArray(valid)) end

    # show(result); print(" | "); show(valid); println() #DEBUG
    @assert (length(vm) == length(result)) #make sure we have all cards sorted
    return result
end

#heuristics of picking the right moves
function pickMoves(g::GameState, vm::CardSet32, depth::Int, α::Int, β::Int)
    vm = removeNeighbours(g, vm)

    vm = sortMoves(g, vm)
    return vm
end

#TODO
#this is the main (recursive loop for the minimax evaluation)
#note: may be implemented partially on the client side in javascript
# type MiniMaxTree
#tree for evaluation
#TODO memoise: make sure states are identical if situation is identical - maybe define custom hash, use IntSets everywhere, and maybe-maybe (doublecheck) leave alpha and beta out of the memoise hash
#BUG players are not in turn. Make sure max-min players are called right (if maximisingPlayer ...). Eventually there will be three types of players with different score goals (csendes ulti miatt)
#TODO: memoize better?
@memoize Dict function alfabeta(g::GameState, depth::Int, α::Int, β::Int)
    #DEBUG
    global abN += 1
    # global gLast = g
    # info(depth, " ", α, " ", β)
    #END DEBUG
    if  depth == 0 || isempty(ps(g, g.currentPlayer).hand)
        return score(g), CardSet32()
    end
    bestMoves = CardSet32()
    if  g.currentPlayer == p1
        for card in pickMoves(g, validMoves(g), depth, α, β)
            #create new state
            gNew = newState(g, card)

            newAlpha, bm = alfabeta(gNew, depth-1, α, β)
            if newAlpha > α
                bestMoves = card
                α = newAlpha
            elseif newAlpha == α
                bestMoves += card
            end

            # # #DEBUG
            # if depth > -2
            #     moves[card] = α
            # end
            # # #END DEBUG

            if β ≤ α
                break                            # (* Beta cut-off *)
            end
        end
        return α, bestMoves
    else
        for card in pickMoves(g, validMoves(g), depth, α, β)
            #create new state
            gNew = newState(g, card)

            newBeta, bm = alfabeta(gNew, depth-1, α, β)
            if newBeta < β
                bestMoves = card
                β = newBeta
            elseif newBeta == β
                bestMoves += card
            end

            # # #DEBUG
            # if depth > -2
            #     moves[card] = α
            # end
            # # #END DEBUG

            if β ≤ α
                break                             # (* Alpha cut-off *)
            end
        end
        return β, bestMoves
    end
end
# Notes:
#
# -'for each child of node' - Rather than editing the state of the current board, create an entirely new board that is the result of applying the move. By using immutable objects, your code will be less prone to bugs and quicker to reason about in general.
#
# -To use this method, call it for every possible move you can make from the current state, giving it depth -1, -Infinity for alpha and +Infinity for beta, and it should start by being the non-moving player's turn in each of these calls - the one that returns the highest value is the best one to take.
#
# It's very very conceptually simple. If you code it right then you never instantiate more than (depth) boards at once, you never consider pointless branches and so on.

#generate (sub)hands from the current state, and run
#Monte Carlo simulation on it
#use existing probabilities for hand distribution
#and hypergeometric probabilities otherwise
function MonteCarlo()

    #@parallel (vcat) for ... end | or use SharedArray to accumulate results in parallel for
end

##############
#Tests
##############

##############
#Erdekes partik
##############
erdekesPartik = [
"pAFZ8 mAFXX zTKYY | pTX7 mTX zAXX tYY | tA xXXXXXXXXX", #vege: asszal ut, F-vel kihajtja a 10-est, 8-asra adu hetes marad bennt
"pA p[FK] p[78] b[FK]YY | p0 bAT cA dAY | xXXXXX", #az indokolatlan kontra
"pAU8 mAFUX tAT zU| zATXX7 pT7 mAX tX | xXXXXX | mY tY",
]

max_stars_repo_name: albapa/ulti

max_issues_repo_name: albapa/ulti

max_forks_repo_name: albapa/ulti


text:
# regular prior for Pose2

export PriorPose3, PackedPriorPose3

"""
$(TYPEDEF)

Direct observation information of `Pose3` variable type.
"""
mutable struct PriorPose3{T <: IIF.SamplableBelief} <: IncrementalInference.AbstractPrior
  Zi::T
  # empty constructor
  PriorPose3{T}() where T = new{T}()
  # regular constructor
  PriorPose3{T}(z::T) where {T <: IIF.SamplableBelief} = new{T}(z)
end
# convenience and default object helper
PriorPose3(z::T=MvNormal(zeros(6), LinearAlgebra.diagm([0.01*ones(3);0.0001*ones(3)]) )) where {T <: SamplableBelief} = PriorPose3{T}(z)

# Standardized sampling function
function getSample(cfo::CalcFactor{<:PriorPose3}, N::Int=1)
  return (rand(cfo.factor.Zi, N),)
end

"""
$(TYPEDEF)

Serialization type for PriorPose3.
"""
mutable struct PackedPriorPose3  <: IncrementalInference.PackedInferenceType
    Zi::String
    PackedPriorPose3() = new()
    PackedPriorPose3(x::AbstractString) = new(x)
end
function convert(::Type{PriorPose3}, packed::PackedPriorPose3)
  return PriorPose3( convert(SamplableBelief, packed.Zi) )
end
function convert(::Type{PackedPriorPose3}, obj::PriorPose3)
  return PackedPriorPose3(convert(PackedSamplableBelief, obj.Zi))
end


text:
using PowerDynamics
using NetworkDynamics
using LaTeXStrings
include("plotting.jl")
include("ThirdOrderEq.jl")
include("ThirdOrderEqGovernorIEEEG1.jl")
include("ThirdOrderEqExciterIEEEDC1A.jl")

line_list = []
append!(line_list, [StaticLine(from = 1, to = 2, Y = -1im / 0.02)])
append!(line_list, [StaticLine(from = 2, to = 3, Y = -1im / 0.02)])
append!(line_list, [StaticLine(from = 3, to = 4, Y = -1im / 0.02)])
append!(line_list, [StaticLine(from = 1, to = 4, Y = -1im / 0.02)])

node_list = []
append!(node_list, [SlackAlgebraic(U = 1, Y_n = 1)])
append!(node_list, [ThirdOrderEq(H = 5.0, P = 0.6, D = 0.01, Ω = 50, E_f = 0.2, X_d_dash = 0.111, T_d_dash = 0.1, X_q_dash = 0.103, X_d = 0.1)])
append!(node_list, [ThirdOrderEq(H = 5.0, P = 0.6, D = 0.01, Ω = 50, E_f = 0.2, X_d_dash = 0.111, T_d_dash = 0.1, X_q_dash = 0.103, X_d = 0.1)])
#append!(node_list, [ThirdOrderEq(H = 5.0, P = 0.6, D = 0.01, Ω = 50, E_f = 0.2, X_d_dash = 0.111, T_d_dash = 0.1, X_q_dash = 0.103, X_d = 0.1)])
append!(node_list, [ThirdOrderEqGovernorIEEEGone(H = 5.0, D = 0.1, Ω = 50, E_f = 0.2, T_d_dash = 0.1, X_d_dash = 0.5,  X_q_dash = 0.103, X_d = 0.5, P0 = 0.5, Pmax = 1.03, Pmin = 0, Pup = 1.38, Pdown = -1.25, T1 = 0.82, T2 = 0, T3 = 0.97, K = 20.0)])

powergrid = PowerGrid(node_list, line_list)

operationpoint = find_steady_state(powergrid)

result = simulate(
    Perturbation(4, :ω, Inc(10)),
    powergrid,
    operationpoint,
    timespan = (0.0, 28.0),
)

plot_res(result, powergrid, 4)

max_stars_repo_name: luap-pik/PowerDynIC

max_issues_repo_name: luap-pik/PowerDynIC

max_forks_repo_name: luap-pik/PowerDynIC


text:
# using MCMC to fit a Bayesian model

# inverse
using KissMCMC
using StatsBase
using OnlineStats
import JLD2
const JLD=JLD2

export thick_err, iv2d_err, iv1d_err, make_loglikelihood, make_logposterior,
    MCMCPara, mcmc, Theta0

# misc tools
include("tools.jl")

include("probabilistic-model.jl")
include("priors.jl")

# numerical fitting
include("numerics.jl")

include("plotting-inverse.jl")
include("postproc.jl")
include("evaluation.jl")


text:
using CSV
using DataFrames

# Random ECG of length 40s with a known number of 52 Beats
ecg = CSV.read(joinpath(dirname(@__FILE__), "data", "ecg.csv"), DataFrame)[!, :ECG]

sf = 250
peaks = 53
@testset "Test R peak detection" begin
    @test length(detectRPeaks(ecg, sf)) == peaks
    @test length(detectRPeaks(ecg, 249)) == peaks
    @test length(detectRPeaks(ecg, 249.141)) == peaks
    @test length(detectRPeaks(ecg, 249.141, minPeakDist = 0.2)) == peaks
    @test_throws DomainError detectRPeaks(ecg, 40)
end

@testset "Test baseline detection" begin
    @test_throws  DomainError getECGBaseline(collect(1:10), 5)
    @test length(getECGBaseline(ecg, 250)) == length(ecg)
end

sf = 1000
peaks = 19

bp = CSV.read(joinpath(dirname(@__FILE__), "data", "BP.csv"), DataFrame)[!, :BP]
@testset "Test PW peak detection" begin
    @test length(detectPWPeaks(bp, sf)) == peaks
    @test length(detectPWPeaks(bp, sf-1)) == peaks
    @test length(detectPWPeaks(bp, sf - 0.891)) == peaks
    @test_throws DomainError detectRPeaks(bp, 10)
end


text:
# Basic Linear algebra
# https://www.juliabox.com/notebook/notebooks/tutorials/intro-to-julia/11.%20Basic%20linear%20algebra.ipynb

A = [1 2 3;
     4 5 6;
     7 8 9]

@show A
@show A' # conjuate transpose
@show transpose(A)

x = fill(1.0, (3,))

@show typeof(A), typeof(x)

@show A * x
@show A' * A
@show A'A #allows us to write this wiout `*`

# Solving linear systems
# The problem `Ax=b` for square A is solved by `\` function

A = [1 3 2;
     3 2 4;
     1 9 3]

b = A*x
@show A\b # expected the value of `x`

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
if GPU
    import CUDNN: cudnnConvolutionForward, cudnnConvolutionBackwardFilter, cudnnConvolutionBackwardData, cudnnPoolingForward, cudnnPoolingBackward
    using CUDNN: CUDNN_CONVOLUTION, CUDNN_CROSS_CORRELATION, CUDNN_CONVOLUTION_FWD_ALGO_IMPLICIT_GEMM, CUDNN_CONVOLUTION_FWD_ALGO_IMPLICIT_PRECOMP_GEMM, CUDNN_CONVOLUTION_FWD_ALGO_GEMM, CUDNN_CONVOLUTION_FWD_ALGO_DIRECT, CUDNN_CONVOLUTION_FWD_ALGO_FFT, CUDNN_POOLING_MAX, CUDNN_POOLING_AVERAGE_COUNT_INCLUDE_PADDING, CUDNN_POOLING_AVERAGE_COUNT_EXCLUDE_PADDING
else
    const CUDNN_CONVOLUTION = (UInt32)(0)
    const CUDNN_CROSS_CORRELATION = (UInt32)(1)
    const CUDNN_CONVOLUTION_FWD_ALGO_IMPLICIT_GEMM = (UInt32)(0)
    const CUDNN_CONVOLUTION_FWD_ALGO_IMPLICIT_PRECOMP_GEMM = (UInt32)(1)
    const CUDNN_CONVOLUTION_FWD_ALGO_GEMM = (UInt32)(2)
    const CUDNN_CONVOLUTION_FWD_ALGO_DIRECT = (UInt32)(3)
    const CUDNN_CONVOLUTION_FWD_ALGO_FFT = (UInt32)(4)
    const CUDNN_POOLING_MAX = (UInt32)(0)
    const CUDNN_POOLING_AVERAGE_COUNT_INCLUDE_PADDING = (UInt32)(1)
    const CUDNN_POOLING_AVERAGE_COUNT_EXCLUDE_PADDING = (UInt32)(2)
end

function _conv2_gemm{T}(x0::Array{T,2}, w::Array{T,2}; pad=0, stride=1, xcorr=false)
    if pad > 0
        x=zeros(eltype(x0),map(m->2pad+m,size(x0)))
        x[pad+1:end-pad,pad+1:end-pad] = x0
    else
        x=x0
    end
    rwindow, cwindow = size(w)
    row_extend = size(x,1)-rwindow+1
    col_extend = size(x,2)-cwindow+1

    widx = [(j-1)*size(x,1)+i for i in 1:row_extend, j in 1:col_extend]
    oidx = [(j-1)*size(x,1)+i for i in 1:rwindow, j in 1:cwindow]
    destidx = [i+(j-1) for i in widx, j in oidx]
    return reshape(x[destidx]*(xcorr ? w[:] : reverse(w[:])),row_extend,col_extend)
end

function _conv2{T}(x::Array{T,2}, w::Array{T,2}; pad=0, stride=1, xcorr=false)
    max_pad = map(x->x-1-pad,size(w))
    y = conv2(x, xcorr ? rot180(w) : w)
    return y[1+max_pad[1]:stride:end-max_pad[1], 1+max_pad[2]:stride:end-max_pad[2]]
end

function cudnnConvolutionForward{T}(x::Array{T,4}, w::Array{T,4}, y::Array{T,4}; padding=0, stride=1,
                                    upscale=1, mode=CUDNN_CONVOLUTION, cd=nothing,
                                    algorithm=CUDNN_CONVOLUTION_FWD_ALGO_IMPLICIT_GEMM,
                                    workSpace=0, workSpaceSizeInBytes=0, alpha=1, beta=1,im2col=1)
    # x: (W,H,C,N)
    # w: (W,H,C,K)
    # y: (W,H,K,N)
    fill!(y,0)
    @assert (padding==0 && stride==1 && upscale==1 && mode==CUDNN_CONVOLUTION && algorithm == CUDNN_CONVOLUTION_FWD_ALGO_IMPLICIT_GEMM) "$((padding,stride,upscale,mode,algorithm))"
    Wx,Hx,Cx,N = size(x)
    Ww,Hw,Cw,K = size(w)
    @assert (Cx==Cw && Hx>=Hw && Wx>=Ww) "$((Wx,Hw,Ww,Hw))"

    @inbounds for n in 1:N, k in 1:K, c in 1:Cx
        y[:,:,k,n] += _conv2_gemm(x[:,:,c,n], w[:,:,c,k]; pad=padding, stride=stride, xcorr=mode!=0)
    end
    return y
end

# dw = rot180(xcorr(x,dy))
function cudnnConvolutionBackwardFilter{T}(x::Array{T,4}, dy::Array{T,4}, dw::Array{T,4}; padding=0, stride=1, upscale=1, mode=CUDNN_CONVOLUTION)
    # x:    (Wx,Hx,Cx,N)
    # dy:   (Wy,Hy,K,N)
    # dw:    (Ww,Hw,Cw,K)
    fill!(dw,0)
    @assert (padding==0&& stride==1&& upscale==1&& mode==CUDNN_CONVOLUTION)
    Wx,Hx,C,Nx = size(x)
    Wy,Hy,K,Ny = size(dy)
    @inbounds for c in 1:C, k in 1:K, n in 1:Ny
        dw[:,:,c,k] += rot180(_conv2_gemm(x[:,:,c,n], dy[:,:,k,n]; pad=padding, stride=stride, xcorr=true))
    end
    return dw
end

# dx = xcorr(dy, w, 'full')
function cudnnConvolutionBackwardData{T}(w::Array{T,4}, dy::Array{T,4}, dx::Array{T,4}; padding=0, stride=1, upscale=1, mode=CUDNN_CONVOLUTION)
    fill!(dx,0)
    @assert (padding==0&& stride==1&& upscale==1&& mode==CUDNN_CONVOLUTION)
    Wy,Hy,Ky,N = size(dy)
    Ww,Hw,C,Kw = size(w)
    @assert Ky==Kw
    @inbounds for n in 1:N, c in 1:C, k in 1:Kw
        t = _conv2_gemm(dy[:,:,k,n], w[:,:,c,k]; xcorr=true, pad=Ww-1)
        # t = conv2(dy[:,:,k,n], rot180(w[:,:,c,k]))
        # t = _conv2(dy[:,:,k,n], w[:,:,c,k]; pad=Ww-1, stride=stride, xcorr=true)
        dx[:,:,c,n] += t
    end
    return dx
end

function cudnnPoolingForward{T}(x::Array{T,4}, y; window=2, padding=0, stride=window, mode=CUDNN_POOLING_MAX)
    fill!(y,0)
    @assert (padding==0 &&  mode==CUDNN_POOLING_MAX &&  stride==window)
    # x: (W,H,C,N)
    Wx,Hx,C,Nx = size(x);
    Wy,Hy,K,Ny = size(y);
    @assert (Nx == Ny && C==K)
    @inbounds for n in 1:Nx, c in 1:C, j in 1:stride:Hx, i in 1:stride:Wx
        iy, jy = div(i,stride)+1, div(j,stride)+1
        hx_end = j+window-1 > Hx ? Hx : j+window-1
        wx_end = i+window-1 > Hx ? Hx : i+window-1
        y[iy,jy,c,n] = maximum(x[i:wx_end,j:hx_end,c,n])
    end
    return y
end

function cudnnPoolingBackward{T}(y::Array{T,4}, dy::Array{T,4}, x::Array{T,4}, dx::Array{T,4}; window=2, padding=0, stride=1, mode=CUDNN_POOLING_MAX)
    fill!(dx,0)
    @assert (padding==0 && stride==window && mode==CUDNN_POOLING_MAX)
    # x: (W,H,C,N)
    Wx,Hx,C,Nx = size(x);
    Wy,Hy,K,Ny = size(y);
    @assert (Nx == Ny && C==K)
    @inbounds for n in 1:Nx, c in 1:C, j in 1:stride:Hx, i in 1:stride:Wx
        iy, jy = div(i,stride)+1, div(j,stride)+1
        hx_end = j+window-1 > Hx ? Hx : j+window-1
        wx_end = i+window-1 > Hx ? Hx : i+window-1
        a = x[i:wx_end,j:hx_end,c,n]
        di,dj = ind2sub(a,indmax(a))
        dx[i+di-1,j+dj-1,c,n] += dy[iy,jy,c,n]
    end
    return dx
end

function cudnnGetConvolutionNdForwardOutputDim{T}(x::Array{T,4}, w::Array{T,4}; padding=padding,stride=stride)
    Wx,Hx,Cx,N = size(x)
    Ww,Hw,Cw,K = size(w)
    @assert Cx==Cw
    Wy,Hy = floor(Int, 1 + (Int[Wx,Hx] + 2*padding - Int[Ww,Hw]) / stride)
    return (Wy,Hy,K,N)
end

function cudnnGetPoolingNdForwardOutputDim{T}(x::Array{T,4}; window=2, padding=0, stride=1, mode=CUDNN_POOLING_MAX)
    # @assert padding==0 && stride==1 && mode==CUDNN_POOLING_MAX
    dims = [size(x)...]
    # (mode, pdims, window, padding, stride) = cudnnGetPoolingNdDescriptor(pd)
    for i=1:length(dims)-2
        # dims[i] = 1 + floor((dims[i] + 2*padding - window) / stride)
        dims[i] = 1 + ceil((dims[i] + 2*padding - window) / stride)
    end
    tuple(dims...)
end


text:
module TestExamples2DEuler

using Test
using Trixi

include("test_trixi.jl")

# pathof(Trixi) returns /path/to/Trixi/src/Trixi.jl, dirname gives the parent directory
EXAMPLES_DIR = joinpath(pathof(Trixi) |> dirname |> dirname, "examples", "tree_2d_dgsem")

@testset "Compressible Euler" begin
  @trixi_testset "elixir_euler_source_terms.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_source_terms.jl"),
      l2   = [8.517808508019351e-7, 1.2350203856098537e-6, 1.2350203856728076e-6, 4.277886946638239e-6],
      linf = [8.357848139128876e-6, 1.0326302096741458e-5, 1.0326302101404394e-5, 4.496194024383726e-5])
  end

  @trixi_testset "elixir_euler_convergence_pure_fv.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_convergence_pure_fv.jl"),
      l2   = [0.026440292358506527, 0.013245905852168414, 0.013245905852168479, 0.03912520302609374],
      linf = [0.042130817806361964, 0.022685499230187034, 0.022685499230187922, 0.06999771202145322])
  end

  @trixi_testset "elixir_euler_density_wave.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_density_wave.jl"),
      l2   = [0.0010600778457964775, 0.00010600778457634275, 0.00021201556915872665, 2.650194614399671e-5],
      linf = [0.006614198043413566, 0.0006614198043973507, 0.001322839608837334, 0.000165354951256802],
      tspan = (0.0, 0.5))
  end

  @trixi_testset "elixir_euler_source_terms_nonperiodic.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_source_terms_nonperiodic.jl"),
      l2   = [2.3653424742011065e-6, 2.1388875092652046e-6, 2.1388875093792834e-6, 6.010896863407163e-6],
      linf = [1.4080465934984687e-5, 1.7579850582816192e-5, 1.757985059525069e-5, 5.95689353266593e-5])
  end

  @trixi_testset "elixir_euler_ec.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_ec.jl"),
      l2   = [0.061751715597716854, 0.05018223615408711, 0.05018989446443463, 0.225871559730513],
      linf = [0.29347582879608825, 0.31081249232844693, 0.3107380389947736, 1.0540358049885143])
  end

  @trixi_testset "elixir_euler_shockcapturing.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_shockcapturing.jl"),
      l2   = [0.053797946432602085, 0.04696120828935379, 0.04696384063506395, 0.19685320969570913],
      linf = [0.18540158860112732, 0.24029373364236004, 0.23267525584314722, 0.6874555954921888])
  end

  @trixi_testset "elixir_euler_blast_wave.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_blast_wave.jl"),
      l2   = [0.13938437665766198, 0.11552277347912725, 0.11552281282617058, 0.3391365873985103],
      linf = [1.4554615795909327, 1.3256380647458283, 1.3256380647457329, 1.8161954000082436],
      maxiters = 30)
  end

  @trixi_testset "elixir_euler_blast_wave_pure_fv.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_blast_wave_pure_fv.jl"),
      l2   = [0.39957047631960346, 0.21006912294983154, 0.21006903549932, 0.6280328163981136],
      linf = [2.20417889887697, 1.5487238480003327, 1.5486788679247812, 2.4656795949035857],
      tspan = (0.0, 0.5))
  end

  @trixi_testset "elixir_euler_blast_wave_amr.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_blast_wave_amr.jl"),
      l2   = [0.6777595378750056, 0.2813201698283611, 0.2813199969265205, 0.720079968985097],
      linf = [2.889020548335417, 1.8035429384705923, 1.8034830438685023, 3.033907232631001],
      tspan = (0.0, 1.0))
  end

  @trixi_testset "elixir_euler_sedov_blast_wave.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_sedov_blast_wave.jl"),
      l2   = [0.4819231474652106, 0.1655453843075888, 0.16554538430758944, 0.6182451347074112],
      linf = [2.4849541392059202, 1.2816207838889486, 1.2816207838888944, 6.4743208586529875],
      tspan = (0.0, 1.0))
  end

  @trixi_testset "elixir_euler_positivity.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_positivity.jl"),
      l2   = [0.48408821703857413, 0.16610607713651496, 0.16610607713651482, 0.6181627089877778],
      linf = [2.4926053782999267, 1.2903378547503723, 1.2903378547503739, 6.472338570555728],
      tspan = (0.0, 1.0))
  end

  @trixi_testset "elixir_euler_blob_mortar.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_blob_mortar.jl"),
      l2   = [0.22112279110928998, 0.6274680581327996, 0.24324346433572555, 2.925884837102588],
      linf = [10.209709958634505, 28.427728124038527, 9.45255862469211, 97.24756873689934],
      tspan = (0.0, 0.5))
  end

  @trixi_testset "elixir_euler_blob_amr.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_blob_amr.jl"),
      l2   = [0.20227473253519843, 1.1852146938118358, 0.10163211952508579, 5.24084854294007],
      linf = [14.241296885123061, 71.64475965226391, 7.3265332346563605, 293.17007294795917],
      tspan = (0.0, 0.12))
  end

  @trixi_testset "elixir_euler_kelvin_helmholtz_instability.jl" begin
    if Threads.nthreads() == 1
      # This example uses random numbers to generate the initial condition.
      # Hence, we can only check "errors" if everything is made reproducible.
      # However, that's not enough to ensure reproducibility since the stream
      # of random numbers is not guaranteed to be the same across different
      # minor versions of Julia.
      # See https://github.com/trixi-framework/Trixi.jl/issues/232#issuecomment-709738400
      @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_kelvin_helmholtz_instability.jl"),
        l2   = [0.002046113073936985, 0.002862623943300569, 0.001971116879236713, 0.004816623657677065],
        linf = [0.024375050653856478, 0.01803061241763637, 0.009938942915093363, 0.02097211774984231],
        tspan = (0.0, 0.2))
    else
      @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_kelvin_helmholtz_instability.jl"),
        tspan = (0.0, 0.2))
    end
  end

  @trixi_testset "elixir_euler_kelvin_helmholtz_instability_amr.jl" begin
    if Threads.nthreads() == 1
      # This example uses random numbers to generate the initial condition.
      # Hence, we can only check "errors" if everything is made reproducible.
      # However, that's not enough to ensure reproducibility since the stream
      # of random numbers is not guaranteed to be the same across different
      # minor versions of Julia.
      # See https://github.com/trixi-framework/Trixi.jl/issues/232#issuecomment-709738400
      @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_kelvin_helmholtz_instability_amr.jl"),
        l2   = [0.0015773943210539198, 0.0023012411276933705, 0.0013793453564426015, 0.003239622733845669],
        linf = [0.025791632274568155, 0.015170620391733003, 0.00958159135384009, 0.013265593557157729],
        tspan = (0.0, 0.2))
    else
      @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_kelvin_helmholtz_instability_amr.jl"),
        tspan = (0.0, 0.2))
    end
  end

  @trixi_testset "elixir_euler_vortex.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_vortex.jl"),
      l2   = [3.6343141788303172e-6, 0.0032111379945554378, 0.0032111482803927763, 0.004545715899533244],
      linf = [7.901851921288117e-5, 0.030561510550305093, 0.03050266851613559, 0.042876690674344076])
  end

  @trixi_testset "elixir_euler_vortex_mortar.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_vortex_mortar.jl"),
      l2   = [2.1203067597905396e-6, 2.792922844344066e-5, 3.7593422075434476e-5, 8.813644239472934e-5],
      linf = [5.932046017209647e-5, 0.0007491265252463908, 0.0008165690091972433, 0.0022122638482677814])
  end

  @trixi_testset "elixir_euler_vortex_mortar_split.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_vortex_mortar_split.jl"),
      l2   = [2.1203686464224497e-6, 2.8053511778944724e-5, 3.7617944213324846e-5, 8.840785095671702e-5],
      linf = [5.900566824046383e-5, 0.0007554116580206216, 0.0008166047918783947, 0.0022090163480150693])
  end

  @trixi_testset "elixir_euler_vortex_mortar_split.jl with flux_central" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_vortex_mortar_split.jl"),
      l2   = [2.1203067597988456e-6, 2.7929228443420677e-5, 3.7593422075440256e-5, 8.813644239501579e-5],
      linf = [5.932046017187442e-5, 0.0007491265252450585, 0.0008165690091986866, 0.0022122638482677814],
      volume_flux = flux_central)
  end

  @trixi_testset "elixir_euler_vortex_mortar_split.jl with flux_shima_etal" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_vortex_mortar_split.jl"),
      l2   = [2.120102590290442e-6, 2.8056524989269332e-5, 3.7595003938956103e-5, 8.84137216636431e-5],
      linf = [5.9341032821902395e-5, 0.000755231667019185, 0.0008152449093754566, 0.0022069874183223703],
      volume_flux = flux_shima_etal)
  end

  @trixi_testset "elixir_euler_vortex_mortar_split.jl with flux_ranocha" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_vortex_mortar_split.jl"),
      l2   = [2.1201025794676757e-6, 2.805652743675726e-5, 3.759500922486255e-5, 8.841377002453846e-5],
      linf = [5.93402785843411e-5, 0.0007552314167595942, 0.0008152450162266511, 0.002206997654994325],
      volume_flux = flux_ranocha)
  end

  @trixi_testset "elixir_euler_vortex_shockcapturing.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_vortex_shockcapturing.jl"),
      l2   = [3.803468104603101e-6, 5.561125881446469e-5, 5.564051742708959e-5, 0.00015706487208149976],
      linf = [8.49343523420254e-5, 0.0009602338002179245, 0.0009669941278601657, 0.0030758222779923017])
  end

  @trixi_testset "elixir_euler_vortex_mortar_shockcapturing.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_vortex_mortar_shockcapturing.jl"),
      l2   = [2.120587333728832e-6, 2.8053708180265217e-5, 3.7618048362912055e-5, 8.841529146530808e-5],
      linf = [5.900528113544912e-5, 0.0007554176274675584, 0.000816603085540657, 0.002208999078852969])
  end

  @trixi_testset "elixir_euler_vortex_amr.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_vortex_amr.jl"),
      l2   = [2.1325938318497577e-6, 0.003281545066125044, 0.0032806321457482615, 0.00464605975938406],
      linf = [4.667437348149228e-5, 0.03175420871507906, 0.0318039789241531, 0.04561735256198318])
  end
end

end # module


text:
using Base: @propagate_inbounds, OneTo

import LinearAlgebra: lmul!, rmul!
import Base: Matrix

@propagate_inbounds is_offdiagonal_small(H::AbstractMatrix{T}, i::Int, tol = eps(real(T))) where {T} =
    abs(H[i+1,i]) ≤ tol*(abs(H[i,i]) + abs(H[i+1,i+1]))

function double_shift_schur!(H::AbstractMatrix{Tv}, from::Int, to::Int,
                             μ::Complex, Q = NotWanted()) where {Tv<:Real}
    m, n = size(H)

    # Compute the nonzero entries of p = (H - μI)(H - μ'I)e₁.
    # Because of the Hessenberg structure we only need H[min:min+2,min:min+1] to form p.
    @inbounds H₁₁ = H[from+0,from+0]
    @inbounds H₂₁ = H[from+1,from+0]

    @inbounds H₁₂ = H[from+0,from+1]
    @inbounds H₂₂ = H[from+1,from+1]
    @inbounds H₃₂ = H[from+2,from+1]

    p₁ = abs2(μ) - 2real(μ) * H₁₁ + H₁₁ * H₁₁ + H₁₂ * H₂₁
    p₂ = -2real(μ) * H₂₁ + H₂₁ * H₁₁ + H₂₂ * H₂₁
    p₃ = H₃₂ * H₂₁

    # Map that column to a mulitiple of e₁ via two Given's rotations
    G₁, nrm = get_rotation(p₁, p₂, p₃, from)

    # Apply the Given's rotations
    lmul!(G₁, view(H, :, from:n))
    rmul!(view(H, OneTo(min(from+3, m)), :), G₁)
    rmul!(Q, G₁)

    # Bulge chasing. First step of the for-loop below looks like:
    #  from           to
    #     ↓           ↓
    #     x x x x x x x     x x x x x x x     x + + + x x x
    # i → x x x x x x x     + + + + + + +     x + + + x x x
    #     x x x x x x x     o + + + + + +       + + + x x x
    #     x x x x x x x  ⇒  o + + + + + +  ⇒    + + + x x x
    #       |   x x x x           x x x x       + + + x x x
    #       |     x x x             x x x             x x x
    #       |       x x               x x               x x
    #       ↑
    #       i
    #
    # Last iterations looks like:
    #  from           to
    #     ↓           ↓
    #     x x x x x x x     x x x x x x x     x x x x + + +
    #     x x x x x x x     x x x x x x x     x x x x + + +
    #       x x x x x x       x x x x x x       x x x + + +
    #         x x x x x  ⇒    x x x x x x  ⇒      x x + + +
    # i → ----- x x x x           + + + +           x + + +
    #           x x x x           o + + +             + + +
    #           x x x x           o + + +             + + +
    #             ↑
    #             i

    @inbounds for i = from + 1 : to - 2
        p₁ = H[i+0,i-1]
        p₂ = H[i+1,i-1]
        p₃ = H[i+2,i-1]

        G, nrm = get_rotation(p₁, p₂, p₃, i)

        # First column is done by hand
        H[i+0,i-1] = nrm
        H[i+1,i-1] = zero(Tv)
        H[i+2,i-1] = zero(Tv)

        # Rotate remaining columns
        lmul!(G, view(H, :, i:n))

        # Create a new bulge
        rmul!(view(H, OneTo(min(i + 3, m)), :), G)
        rmul!(Q, G)
    end

    # Last bulge is just one Given's rotation
    #     from          to
    #       ↓           ↓
    #from → x x x x x x x    x x x x x x x    x x x x x + +
    #       x x x x x x x    x x x x x x x    x x x x x + +
    #         x x x x x x      x x x x x x      x x x x + +
    #           x x x x x  ⇒     x x x x x  ⇒     x x x + +
    #             x x x x          x x x x          x x + +
    #               x x x            + + +            x + +
    # to  → ------- x x x            o + +              + +

    @inbounds Gₙ, nrm = get_rotation(H[to-1,to-2], H[to,to-2], to-1)
    @inbounds H[to-1,to-2] = nrm
    @inbounds H[to,to-2] = zero(Tv)

    lmul!(Gₙ, view(H, :, to-1:n))
    rmul!(view(H, OneTo(to), :), Gₙ)
    rmul!(Q, Gₙ)

    H
end

function single_shift_schur!(H::AbstractMatrix{Tv}, from::Int, to::Int,
                             μ::Number, Q = NotWanted()) where {Tv<:Number}
    m, n = size(H)

    # Compute the nonzero entries of p = (H - μI)e₁.
    @inbounds H₁₁ = H[from+0,from+0]
    @inbounds H₂₁ = H[from+1,from+0]

    p₁ = H₁₁ - μ
    p₂ = H₂₁

    # Map that column to a mulitiple of e₁ via two Given's rotations
    G₁, nrm = get_rotation(p₁, p₂, from)

    # Apply the Given's rotations
    lmul!(G₁, view(H, :, from:n))
    rmul!(view(H, OneTo(min(from + 2, m)), :), G₁)
    rmul!(Q, G₁)

    # Bulge chasing. First step of the for-loop below looks like:
    #  from           to
    #     ↓           ↓
    #     x x x x x x x     x x x x x x x     x + + x x x x
    # i → x x x x x x x     + + + + + + +     x + + x x x x
    #     x x x x x x x     o + + + + + +       + + x x x x
    #         x x x x x  ⇒      x x x x x  ⇒    + + x x x x
    #       |   x x x x           x x x x           x x x x
    #       |     x x x             x x x             x x x
    #       |       x x               x x               x x
    #       ↑
    #       i
    #
    # Last iterations looks like:
    #  from           to
    #     ↓           ↓
    #     x x x x x x x     x x x x x x x     x x x x x + +
    #     x x x x x x x     x x x x x x x     x x x x x + +
    #       x x x x x x       x x x x x x       x x x x + +
    #         x x x x x  ⇒      x x x x x  ⇒      x x x + +
    #           x x x x           x x x x           x x + +
    # i → ------- x x x             + + +             x + +
    #             x x x             0 + +               + +
    #               ↑
    #               i

    @inbounds for i = from + 1 : to - 1
        p₁ = H[i+0,i-1]
        p₂ = H[i+1,i-1]

        G, nrm = get_rotation(p₁, p₂, i)

        # First column is done by hand
        H[i+0,i-1] = nrm
        H[i+1,i-1] = zero(Tv)

        # Rotate remaining columns
        lmul!(G, view(H, :, i:n))

        # Create a new bulge
        rmul!(view(H, OneTo(min(i + 2, m)), :), G)
        rmul!(Q, G)
    end

    H
end

###
### Real arithmetic
###
function local_schurfact!(H::AbstractMatrix{T}, start::Int, to::Int,
                          Q = NotWanted(), tol = eps(T),
                          maxiter = 100*size(H, 1)) where {T<:Real}
    # iteration count
    iter = 0

    @inbounds while true
        iter += 1

        if iter > maxiter
            throw("QR algorithm did not converge")
        end

        # Indexing
        # `to` points to the column where the off-diagonal value was last zero.
        # while `from` points to the smallest index such that there is no small off-diagonal
        # value in columns from:end-1. Sometimes `from` is just 1. Cartoon of a split
        # with from != 1:
        #
        #  + + | | | | + +
        #  + + | | | | + +
        #    o X X X X = =
        #      X X X X = =
        #      . X X X = =
        #      .   X X = =
        #      .     o + +
        #      .     . + +
        #      ^     ^
        #   from   to
        # The X's form the unreduced Hessenberg matrix we are applying QR iterations to,
        # the | and = values get updated by Given's rotations, while the + values remain
        # untouched! The o's are zeros -- or numerically considered zeros.

        # We keep `from` one column past the zero off-diagonal value, so we check whether
        # the `from - 1` column has a small off-diagonal value.
        from = to
        while from > start && !is_offdiagonal_small(H, from - 1, tol)
            from -= 1
        end

        if from == to
            # This just means H[to, to-1] == 0, so one eigenvalue converged at the end
            H[from,from-1] = zero(T)
            to -= 1
        else
            # Now we are sure we can work with a 2×2 block H[to-1:to,to-1:to]
            # We check if this block has a conjugate eigenpair, which might mean we have
            # converged w.r.t. this block if from + 1 == to.
            # Otherwise, if from + 1 < to, we do either a single or double shift, based on
            # whether the H[to-1:to,to-1:to] part has real eigenvalues or a conjugate pair.

            H₁₁, H₁₂ = H[to-1,to-1], H[to-1,to]
            H₂₁, H₂₂ = H[to  ,to-1], H[to  ,to]

            # Scaling to avoid losing precision in the case where we have nearly
            # repeated eigenvalues.
            scale = abs(H₁₁) + abs(H₁₂) + abs(H₂₁) + abs(H₂₂)
            H₁₁ /= scale; H₁₂ /= scale; H₂₁ /= scale; H₂₂ /= scale

            # Trace and discriminant of small eigenvalue problem.
            t = (H₁₁ + H₂₂) / 2
            d = (H₁₁ - t) * (H₂₂ - t) - H₁₂ * H₂₁
            sqrt_discr = sqrt(abs(d))

            # Very important to have a strict comparison here!
            if d < zero(T)
                # Real eigenvalues.
                # Note that if from + 1 == to in this case, then just one additional
                # iteration is necessary, since the Wilkinson shift will do an exact shift.

                # Determine the Wilkinson shift -- the closest eigenvalue of the 2x2 block
                # near H[to,to]

                λ₁ = t + sqrt_discr
                λ₂ = t - sqrt_discr
                λ = abs(H₂₂ - λ₁) < abs(H₂₂ - λ₂) ? λ₁ : λ₂
                λ *= scale

                # Run a bulge chase
                single_shift_schur!(H, from, to, λ, Q)
            else
                # Conjugate pair
                if from + 1 == to
                    # A conjugate pair has converged apparently!
                    if from != 1
                        H[from,from-1] = zero(T)
                    end
                    to -= 2
                else
                    # Otherwise we do a double shift!
                    complex_shift = scale * (t + sqrt_discr * im)
                    double_shift_schur!(H, from, to, complex_shift, Q)
                end
            end
        end

        # Converged!
        to ≤ start && break
    end

    return true
end

###
### Generic implementation
###
function local_schurfact!(H::AbstractMatrix{T}, start::Int, to::Int,
                          Q = NotWanted(), tol = eps(real(T)),
                          maxiter = 100*size(H, 1)) where {T}
    # iteration count
    iter = 0

    @inbounds while true
        iter += 1
        iter > maxiter && return false

        # Indexing, see the real arithmetic version!

        from = to
        while from > start && !is_offdiagonal_small(H, from - 1, tol)
            from -= 1
        end

        if from == to
            # This just means H[to, to-1] == 0, so one eigenvalue converged at the end
            H[from,from-1] = zero(T)
            to -= 1
        else
            # Compute Wilkinson shift
            H₁₁, H₁₂ = H[to-1,to-1], H[to-1,to]
            H₂₁, H₂₂ = H[to  ,to-1], H[to  ,to]
            d = H₁₁ * H₂₂ - H₂₁ * H₁₂
            t = H₁₁ + H₂₂
            sqr = sqrt(t * t - 4d)
            λ₁ = (t + sqr) / 2
            λ₂ = (t - sqr) / 2
            λ = abs(H₂₂ - λ₁) < abs(H₂₂ - λ₂) ? λ₁ : λ₂

            # Run a bulge chase
            single_shift_schur!(H, from, to, λ, Q)
        end

        # Converged!
        to ≤ start && break
    end

    return true
end

local_schurfact!(H::AbstractMatrix{T}, Q = NotWanted(), tol = eps(real(T)), maxiter = 100*size(H, 1)) where {T} =
    local_schurfact!(H, 1, size(H, 2), Q, tol, maxiter)


text:
# UnitCommitment.jl: Optimization Package for Security-Constrained Unit Commitment
# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
# Released under the modified BSD license. See COPYING.md for more details.

using PackageCompiler

using DataStructures
using JSON
using JuMP
using MathOptInterface
using SparseArrays

pkg = [:DataStructures,
       :JSON,
       :JuMP,
       :MathOptInterface,
       :SparseArrays,
       ]

@info "Building system image..."
create_sysimage(pkg,
                precompile_statements_file="build/precompile.jl",
                sysimage_path="build/sysimage.so")


text:
"""
    Readme(;
        file="$(contractuser(default_file("README.md")))",
        destination="README.md",
        inline_badges=false,
    )

Creates a `README` file that contains badges for other included plugins.

## Keyword Arguments
- `file::AbstractString`: Template file for the `README`.
- `destination::AbstractString`: File destination, relative to the repository root.
  For example, values of `"README"` or `"README.rst"` might be desired.
- `inline_badges::Bool`: Whether or not to put the badges on the same line as the package name.
"""
@plugin struct Readme <: FilePlugin
    file::String = default_file("README.md")
    destination::String = "README.md"
    inline_badges::Bool = false
    badge_order::Vector{DataType} = default_badge_order()
end

source(p::Readme) = p.file
destination(p::Readme) = p.destination

function view(p::Readme, t::Template, pkg::AbstractString)
    # Explicitly ordered badges go first.
    strings = String[]
    done = DataType[]
    foreach(p.badge_order) do T
        if hasplugin(t, T)
            append!(strings, badges(getplugin(t, T), t, pkg))
            push!(done, T)
        end
    end
    # And the rest go after, in no particular order.
    foreach(setdiff(map(typeof, t.plugins), done)) do T
        append!(strings, badges(getplugin(t, T), t, pkg))
    end

    return Dict(
        "BADGES" => strings,
        "HAS_CITATION" => hasplugin(t, Citation) && getplugin(t, Citation).readme,
        "HAS_INLINE_BADGES" => !isempty(strings) && p.inline_badges,
        "PKG" => pkg,
    )
end

default_badge_order() = [
    Documenter{GitHubActions},
    Documenter{GitLabCI},
    Documenter{TravisCI},
    GitHubActions,
    GitLabCI,
    TravisCI,
    AppVeyor,
    DroneCI,
    CirrusCI,
    Codecov,
    Coveralls,
    subtypes(BadgePlugin)...
]


text:
#computeOffsets
rowPtrs = Array{Int, 1}(undef, 30)
JuliaPetra.computeOffsets(rowPtrs, 9)
@test collect(1:9:9*30) == rowPtrs

rowPtrs = Array{Int, 1}(undef, 30)
numEnts = collect(1:29)
@test sum(1:29) == JuliaPetra.computeOffsets(rowPtrs, numEnts)
@test [sum(numEnts[1:i-1])+1 for i = 1:30] == rowPtrs

@test_throws InvalidArgumentError JuliaPetra.computeOffsets(rowPtrs, Array{Int, 1}(undef, 30))
@test_throws InvalidArgumentError JuliaPetra.computeOffsets(rowPtrs, Array{Int, 1}(undef, 31))
@test_throws InvalidArgumentError JuliaPetra.computeOffsets(rowPtrs, Array{Int, 1}(undef, 42))


text:
using Documenter
using Bootstrap

makedocs(modules = [Bootstrap],
    sitename = "Bootstrap.jl",
    format = Documenter.HTML(disable_git = true,  # disable source and edit links to github
        canonical = "https://juliangehring.github.io/Bootstrap.jl/stable/"),
    doctest = true,
    checkdocs = :exports,
    linkcheck = true,
    pages = [
        "Home" => "index.md",
        "Library" => "library.md",
        "Changelog" => "NEWS.md",
        "License" => "LICENSE.md"
    ])

deploydocs(
    repo = "github.com/juliangehring/Bootstrap.jl.git",
    target = "build",
    devbranch = "develop",
    push_preview = true
)


text:
@testset "Angles" begin
  @test ∠(P2(0,1),P2(0,0),P2(1,0)) ≈ T(-π/2)
  @test ∠(P2(1,0),P2(0,0),P2(0,1)) ≈ T(π/2)
  @test ∠(P2(-1,0),P2(0,0),P2(0,1)) ≈ T(-π/2)
  @test ∠(P2(0,1),P2(0,0),P2(-1,0)) ≈ T(π/2)
  @test ∠(P2(0,-1),P2(0,0),P2(1,0)) ≈ T(π/2)
  @test ∠(P2(1,0),P2(0,0),P2(0,-1)) ≈ T(-π/2)
  @test ∠(P2(0,-1),P2(0,0),P2(-1,0)) ≈ T(-π/2)
  @test ∠(P2(-1,0),P2(0,0),P2(0,-1)) ≈ T(π/2)
end


text:
# AUTO GENERATED FILE - DO NOT EDIT

export 'feffery'_antdcheckboxgroup

"""
    'feffery'_antdcheckboxgroup(;kwargs...)

An AntdCheckboxGroup component.

Keyword arguments:
- `id` (String; optional)
- `className` (String; optional)
- `disabled` (Bool; optional)
- `loading_state` (optional): . loading_state has the following type: lists containing elements 'is_loading', 'prop_name', 'component_name'.
Those elements have the following types:
  - `is_loading` (Bool; optional): Determines if the component is loading or not
  - `prop_name` (String; optional): Holds which property is loading
  - `component_name` (String; optional): Holds the name of the component that is loading
- `options` (optional): . options has the following type: Array of lists containing elements 'label', 'value', 'disabled'.
Those elements have the following types:
  - `label` (String; optional)
  - `value` (String; optional)
  - `disabled` (Bool; optional)s
- `style` (Dict; optional)
- `value` (Array of Strings; optional)
"""
function 'feffery'_antdcheckboxgroup(; kwargs...)
        available_props = Symbol[:id, :className, :disabled, :loading_state, :options, :style, :value]
        wild_props = Symbol[]
        return Component("'feffery'_antdcheckboxgroup", "AntdCheckboxGroup", "feffery_antd_components", available_props, wild_props; kwargs...)
end

max_stars_repo_name: limmyforget/feffery-antd-components

max_issues_repo_name: limmyforget/feffery-antd-components

max_forks_repo_name: limmyforget/feffery-antd-components


text:
@rule typeof(+)(:out, Marginalisation) (m_in1::PointMass, m_in2::PointMass) = PointMass(mean(m_in1) + mean(m_in2))

@rule typeof(+)(:out, Marginalisation) (m_in1::PointMass, m_in2::UnivariateNormalDistributionsFamily) = NormalMeanVariance(mean(m_in1) + mean(m_in2), var(m_in2))
@rule typeof(+)(:out, Marginalisation) (m_in1::UnivariateNormalDistributionsFamily, m_in2::PointMass) = NormalMeanVariance(mean(m_in1) + mean(m_in2), var(m_in1))

@rule typeof(+)(:out, Marginalisation) (m_in1::PointMass, m_in2::MultivariateNormalDistributionsFamily) = MvNormalMeanCovariance(mean(m_in1) + mean(m_in2), cov(m_in2))
@rule typeof(+)(:out, Marginalisation) (m_in1::MultivariateNormalDistributionsFamily, m_in2::PointMass) = MvNormalMeanCovariance(mean(m_in1) + mean(m_in2), cov(m_in1))

# TODO: Do we need it?
# Specific case for precision parametrisation
@rule typeof(+)(:out, Marginalisation) (m_in1::NormalMeanPrecision, m_in2::PointMass) = NormalMeanPrecision(mean(m_in1) + mean(m_in2), precision(m_in1))
@rule typeof(+)(:out, Marginalisation) (m_in1::PointMass, m_in2::NormalMeanPrecision) = NormalMeanPrecision(mean(m_in1) + mean(m_in2), precision(m_in2))

@rule typeof(+)(:out, Marginalisation) (m_in1::NormalMeanPrecision, m_in2::NormalMeanPrecision) = begin
    min1, pin1 = mean_precision(m_in1)
    min2, pin2 = mean_precision(m_in2)
    return NormalMeanPrecision(min1 + min2, pin1 * pin2 / (pin1 + pin2))
end

@rule typeof(+)(:out, Marginalisation) (m_in1::UnivariateNormalDistributionsFamily, m_in2::UnivariateNormalDistributionsFamily) = begin
    m1, v1 = mean_var(m_in1)
    m2, v2 = mean_var(m_in2)
    return NormalMeanVariance(m1 + m2, v1 + v2)
end

@rule typeof(+)(:out, Marginalisation) (m_in1::MultivariateNormalDistributionsFamily, m_in2::MultivariateNormalDistributionsFamily) = begin
    m1, v1 = mean_cov(m_in1)
    m2, v2 = mean_cov(m_in2)
    return MvNormalMeanCovariance(m1 + m2, v1 + v2)
end