text:
using Test, Parameters

include(joinpath(dirname(@__DIR__), "src", "io.jl"))

# --- Test Utility ---

@test flatten([1, 2]) == [1, 2]
@test flatten([[1, 2], [3, 4, 5]]) == [1, 2, 3, 4, 5]

@test shape([1, 2]) == [2]
@test shape([[1, 2, 3], [4, 5, 6]]) == [3, 2]

@test mdim([1, 2]) == [1, 2]
@test mdim([[1, 2, 3], [4, 5, 6]]) == reshape([1, 2, 3, 4, 5, 6], 3, 2)

# --- Test IO ---

@with_kw struct TestStruct
    b::Bool
    i::Int
    f::Float64
    ia1::Array{Int, 1}
    fa1::Array{Float64, 1}
    ia2::Array{Int, 2}
    fa2::Array{Float64, 2}
    iaa::Array{Array{Int, 1}, 1}
    faa::Array{Array{Float64, 1}, 1}
end

test_struct = TestStruct(
    b=true,
    i=1,
    f=1.1,
    ia1=[1, 2, 3],
    fa1=[1.1, 2.2, 3.3],
    ia2=reshape([1, 2, 3, 4, 5, 6], 2, 3),
    fa2=reshape([1.1, 2.2, 3.3, 4.4, 5.5, 6.6], 2, 3),
    iaa=[[0], [1, 2], [3, 4, 5]],
    faa=[[0.0], [1.1, 2.2], [3.3, 4.4, 5.5]]
)

output_dir = tempdir()
filepath = joinpath(output_dir, "test_struct.json")
save_json(test_struct, filepath)
test_struct2 = load_json(TestStruct, filepath)

for field in fieldnames(TestStruct)
    @test getfield(test_struct, field) == getfield(test_struct2, field)
end

max_stars_repo_name: nnhjy/JuMP-project-template

max_issues_repo_name: nnhjy/JuMP-project-template

max_forks_repo_name: nnhjy/JuMP-project-template


text:
function check_folder(folder::String)
    !isdir(folder) && throw(IS.ConflictingInputsError("Specified folder is not valid"))
    try
        mkdir(joinpath(folder, "fake"))
        rm(joinpath(folder, "fake"))
    catch e
        throw(IS.ConflictingInputsError("Specified folder does not have write access [$e]"))
    end
end

ext: jl


text:
"""
    Raibert hopper policy for Hopper2D
"""

ContactImplicitMPC.@with_kw mutable struct RaibertOptions
    live_plotting::Bool=false # Use the live plotting tool to debug
end

mutable struct Raibert3D <: ContactImplicitMPC.Policy
	kr_c::Real
	kr_p::Real
	kr_v_stance::Real
	kr_v_flight::Real
	kθ_c::Real
	kθ_p::Real
	kθ_v::Real
	u::AbstractVector
	contact::Bool
	v0::Vector{Real}
	Tstance::Real
	Tflight::Real
	q0::AbstractVector
	q1::AbstractVector
	opts::RaibertOptions
end

function raibert_policy(model::ContactImplicitMPC.Hopper3D;
		v0=[0.0; 0.0], Tstance=0.13, Tflight=0.62, h = 0.1,
		kr_c =  8e1,
		kr_p = -1e3,
		kr_v_stance = -1e-2,
		kr_v_flight = -1e1,
		kθ_c =  0.0,
		kθ_p = -6e1,
		kθ_v = -1e1,
		mpc_opts = RaibertOptions(),
		)

	@assert v0[1] == 0.0 || v0[2] == 0.0

	u  = zeros(model.nu)
	q0 = zeros(model.nq)
	q1 = zeros(model.nq)
	contact = false

	Raibert3D(kr_c, kr_p, kr_v_stance, kr_v_flight, kθ_c, kθ_p, kθ_v,
		u, contact, v0, Tstance, Tflight, copy(q0), copy(q1), mpc_opts)
end

function ContactImplicitMPC.policy(p::Raibert3D, x, traj, t)
	# Initialization
	h = traj.h
	p.q0 = copy(traj.q[t])
	p.q1 = copy(x)

	# get Euler angles
	rot0 = RotXYZ(MRP(p.q0[4:6]...))
	eul0 = [rot0.theta1; rot0.theta2; rot0.theta3]

	rot1 = RotXYZ(MRP(p.q1[4:6]...))
	eul1 = [rot1.theta1; rot1.theta2; rot1.theta3]

	# Detect contact
	if any(traj.γ[max(1,t-1)] .> 1.5e-2)
		p.contact = true
	else
		p.contact = false
	end

	# Velocities
	qv = (p.q1 - p.q0)/h
	θv = (eul1 - eul0)/h
	rv = qv[7]

	# References
	θ1 = eul1
	r1 = p.q1[7]
	rref = 0.5

	# td = touchdown
	θtd  = [asin(p.v0[1]*Tstance/(2*rref))/2; asin(p.v0[2]*Tstance/(2*rref))/2]

	# Gains
	kr_c = p.kr_c
	kr_p = p.kr_p
	kr_v_stance = p.kr_v_stance
	kr_v_flight = p.kr_v_flight

	kθ_c = p.kθ_c
	kθ_p = p.kθ_p
	kθ_v = p.kθ_v

	dir = qv[1:2] ./ norm(qv[1:2]) # TODO: change to binary switch

	if p.contact
		# regulate around θtd using same gain kθ_p but taking into account
		# that stance time is much lower than flight time
		p.u[1] = (dir[2] * kθ_c + kθ_p*(θ1[1] + dir[2] * θtd[1]) * p.Tflight/p.Tstance)
		p.u[2] = (dir[1] * kθ_c + kθ_p*(θ1[2] + dir[1] * θtd[2]) * p.Tflight/p.Tstance)
		p.u[3] = kr_c + kr_p*(r1 - rref) + kr_v_stance*rv
	else
		# regulate around θtd
		p.u[1] = (kθ_p*(θ1[1] - dir[2] * θtd[1]))  + kθ_v*θv[1]
		p.u[2] = (kθ_p*(θ1[2] - dir[1] * θtd[2]))  + kθ_v*θv[2]
		p.u[3] = kr_p*(r1 - rref) + kr_v_flight*rv
	end

	p.u .*= h

    return p.u
end


text:
function prob(q)
    p = zeros(Float64, 21, 50)

    p[1, 50] = 50/q
    p[2, 50] = 1 - 50/q

    for k in 49:-1:1
        for n in 0:min(51-k, 20)
            if n == 0
                p[n+1, k] = (k/q) * p[n+1, k+1]
            elseif n + k == 51
                p[n+1, k] = (1 - k/q) * p[n, k+1]
            else
                p[n+1, k] = (k/q) * p[n+1, k+1] + (1 - k/q) * p[n, k+1]
            end
        end
    end

    return p[21, 1]
end

function main()
    l = 50
    h = 55
    e = 10^-12

    while h - l > e
        p = prob((h + l)/2)
        if p > 0.02
            l = (h + l)/2
        else
            h = (h + l)/2
        end
    end

    println(round((h + l)/2, digits=10))
end

main()

max_stars_repo_name: HarrisonGreen/Project-Euler-Solutions

max_issues_repo_name: HarrisonGreen/Project-Euler-Solutions

max_forks_repo_name: HarrisonGreen/Project-Euler-Solutions


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

using Documenter, DocumenterMarkdown, Turing, AdvancedHMC, Bijectors, AdvancedMH
using LibGit2: clone, tag_list, GitRepo
using Documenter: GitHubActions

# Include the utility functions.
include("make-utils.jl")

# Get the version number
trim_version(x) = x[1:findlast('.', x) - 1]

version, is_dev = if haskey(ENV, "TURING_VERSION")
    if ENV["TURING_VERSION"] == "dev"
        ENV["TURING_VERSION"], true
    else
        ENV["TURING_VERSION"], false
    end
else
    if length(ARGS) > 0
        trim_version(ARGS[1]), false
    else
        "dev", true
    end
end

# set default baseurl for the master branch
baseurl = "/" * version
@info "" baseurl

# Make a temporary folder to build from
tmp_path = mktempdir(cleanup=false)

# Paths
## The location of the package to build (Turing())
package_directory = dirname(dirname(pathof(Turing)))
## The location of the docs folder inside that package
docs_path = joinpath(package_directory, "docs")
## The src files for the docs -- markdown documents, typically.
source_path = joinpath(package_directory, "docs", "src")
## The path of turing.ml that we are running this code from.
local_path = @__DIR__
## The place to put the files from source_path after they go through Documenter.jl
build_path = joinpath(tmp_path, "_docs")

# Get any files from Turing's directory
for (root, dirs, files) in walkdir(docs_path)
    new_root = replace(root, docs_path => tmp_path)

    for file in files
        old_file = joinpath(root, file)
        new_file = joinpath(new_root, file)
        @debug "" old_file new_file
        if !isdir(dirname(new_file))
            mkpath(dirname(new_file))
        end
        cp(old_file, new_file)
    end
end

# Copy all the local files to the temporary path
paths = readdir(local_path, join=true)
filter!(x -> !(basename(x) in ["make.jl", "make-utils.jl"]), paths)
for path in paths
    new_path = replace(path, local_path => tmp_path)
    @info "" path new_path

    cp(path, new_path, force=true)
    # if path != new_path
    #     cp(path, new_path, force=true)
    # end
end

# Build docs
with_clean_docs(source_path, build_path) do source, build
    makedocs(
        sitename = "Turing.jl",
        source = source,
        build = build,
        format = Markdown(),
        checkdocs = :all,
    )
end

# You can skip this part if you are on a metered
# connection by calling `julia make.jl no-tutorials`
tutorial_path = joinpath(tmp_path, "_tutorials")
in("no-tutorials", ARGS) || copy_tutorial(tutorial_path)

# deploy
old_jekyll_build = joinpath(local_path, "jekyll-build")
new_jekyll_build = joinpath(tmp_path, "jekyll-build")

# Move jekyll-build to the temporary path
cp(old_jekyll_build, new_jekyll_build, force=true)
with_baseurl(() -> run(`$new_jekyll_build`), baseurl, joinpath(local_path, "_config.yml"))

# Copy assets to folder
cp(joinpath(tmp_path, "assets"), joinpath(tmp_path, "_site", "assets"), force=true)

repo = "github.com:TuringLang/turing.ml.git"

deploy_config = GitHubActions(
    "TuringLang/turing.ml", #github_repository::String
    "push", #github_event_name::String
    is_dev ? "refs/heads/master" : "refs/tags/$(ARGS[1])" #github_ref::String
)

deploydocs(
    target = joinpath(tmp_path, "_site"),
    repo = repo,
    branch = "gh-pages",
    devbranch = "master",
    devurl = "dev",
    versions = ["stable" => "v^", "v#.#", "dev" => "dev"],
    deploy_config = deploy_config
)

@info "" get(ENV, "GITHUB_REF", missing)


text:
const DATA_DIR = joinpath(MODULE_DIR, "..", "data")

const COERCE_ADULT = (
    :age => Continuous,
    :workclass => Multiclass,
    :fnlwgt => Continuous,
    :education => Multiclass,
    :education_num => Continuous,
    :marital_status => Multiclass,
    :occupation => Multiclass,
    :relationship => Multiclass,
    :race => Multiclass,
    :sex => Multiclass,
    :capital_gain => Continuous,
    :capital_loss => Continuous,
    :hours_per_week => Continuous,
    :native_country => Multiclass,
    :income_per_year => Multiclass,
)

"""
Checks whether the dataset is already present in data directory. Downloads it if not present.
"""
function ensure_download(url::String, file::String)
    cd(DATA_DIR) # This is to ensue that the dataset is not downloaded to /tmp instead of ./data
    fpath = joinpath(DATA_DIR, file)
    if !isfile(fpath)
        download(url, file)
    end
end

"
Macro to read csv file of job data and convert columns to categorical.
Returns the tuple (X, y, ŷ)
"
macro load_toydata()
    quote
        fpath = joinpath(DATA_DIR, "jobs.csv")
        data = DataFrame!(CSV.File(fpath))
        categorical!(data, names(data)[1:6])
        (data[!, names(data)[1:4]], data[!, :Class], data[!, :Pred])
    end
end

"
Macro to create fairness Tensor for jobs.csv
The fairness tensor will be created on the basis of the column Job Type.
This column has 3 different values for job types.
"
macro load_toyfairtensor()
    quote
        X, y, ŷ = @load_toydata
        fair_tensor(ŷ, y, X[!, names(X)[4]])
    end
end

"
Macro to load COMPAS dataset.
https://github.com/propublica/compas-analysis/blob/master/Compas%20Analysis.ipynb
"
macro load_compas()
    quote
        url = "https://raw.githubusercontent.com/propublica/compas-analysis/master/compas-scores-two-years.csv"
        fname = "compas-scores-two-years.csv"
        ensure_download(url, fname)

        fpath = joinpath(DATA_DIR, fname)
        data = DataFrame!(CSV.File(fpath))
        data = data[!, ["sex", "age", "age_cat", "race", "c_charge_degree", "priors_count", "days_b_screening_arrest", "decile_score", "is_recid"]]
        dropmissing!(data, disallowmissing=true)
        coerce!(data, Textual => Multiclass)
        coerce!(data, :is_recid => Multiclass)
        y, X = unpack(data, ==(:is_recid), col -> true)

        X = data[!, ["sex", "age", "age_cat", "race", "c_charge_degree", "priors_count", "days_b_screening_arrest", "decile_score"]]
        y = data[!, "is_recid"]
        (X, y)
    end
end

"Macro to load Adult dataset."
macro load_adult()
    quote
        url = "https://archive.ics.uci.edu/ml/machine-learning-databases/adult/adult.data"
        fname = "adult.data"
        cols = ["age", "workclass", "fnlwgt", "education",
            "education_num", "marital_status", "occupation",
            "relationship", "race", "sex", "capital_gain",
            "capital_loss", "hours_per_week", "native_country",
            "income_per_year"
        ]
        ensure_download(url, fname)
        fpath = joinpath(DATA_DIR, fname)
        data = DataFrame!(CSV.File(fpath, header=cols))

        data = dropmissing(data, names(data))
        data.income_per_year = map(data.income_per_year) do η
            η == " <=50K" ? 0 : 1
        end

        coerce!(data, COERCE_ADULT...)
        coerce!(data, :income_per_year => Multiclass)
        y, X = unpack(data, ==(:income_per_year), col -> true)
        (X, y)
    end
end

"""
Load the full version of [German credit dataset](https://archive.ics.uci.edu/ml/datasets/Statlog+%28German+Credit+Data%29).
This dataset has 20 features and 1000 rows. The protected attributes are gender_status and age (>25 is priviledged)
Using the 20 features, it classifies the credit decision to a person as good or bad credit risks.

Returns (X, y)
"""
macro load_german()
    quote
        url = "https://archive.ics.uci.edu/ml/machine-learning-databases/statlog/german/german.data"
        fname = "german.data"
        cols = [
            "status", "duration", "credit_history",
            "purpose", "credit_amount", "savings", "employment",
            "installment_rate", "gender_status",
            "other_debtors", "residence_since", "property", "age",
            "installment_plans", "housing", "existing_credits",
            "skill_level", "people_liable", "telephone",
            "foreign_worker", "label"
        ]
        ensure_download(url, fname)
        fpath = joinpath(DATA_DIR, fname)
        df = DataFrame!(CSV.File(fpath, header=cols))

        # As of now the Package supports only labels which are subset of (0, 1)
        # But the german dataset has labels (1, 2). So, we need to map the labels {1=>0, 2=>1}
        # This restriction will soon be removed and levels(label) shall be used to decide favourable label
        df.label = map(df.label) do η
            η == 2 ? 1 : 0
        end

        gender_status = Dict(
            "A91" => "male_divorced_separated",
            "A92" => "female_divorced_separated_married",
            "A93" => "male_single",
            "A94" => "male_married_widowed",
            "A95" => "female_single", # There seems to be no row with this value in data
        )

        df.gender_status = map(df.gender_status) do η
            gender_status[η]
        end

        coerce!(df, Textual => Multiclass)
        coerce!(df, :label => Multiclass)
        y, X = unpack(df, ==(:label), col -> true);
        (X, y)
    end
end


text:
#=
Afft.jl
2019-07-07, Jeff Fessler, University of Michigan
2020-06-30 in-place fft!
=#

export Afft

# using MIRT: embed!, getindex!
using LinearMapsAA: LinearMapAA, LinearMapAM, LinearMapAO
using FFTW: fft!, bfft!

"""
    A = Afft(samp::AbstractArray{Bool} ; T::DataType = ComplexF32)
Make a LinearMapAO object for (under-sampled) FFT, of type `T`,
using given sampling pattern `samp`.
Especially for compressed sensing MRI with Cartesian sampling.

Option:
- `operator::Bool=true` set to `false` to return a `LinearMapAM`
- `work::AbstractArray` work space for in-place fft operations
"""
function Afft(samp::AbstractArray{Bool,D} ;
    T::DataType = ComplexF32,
    operator::Bool = true, # !
    work::AbstractArray{S,D} = Array{T}(undef, size(samp)),
) where {D, S <: Number}

    forw! = (y,x) -> getindex!(y, fft!(copyto!(work,x)), samp)
    back! = (x,y) -> bfft!(embed!(x,y,samp))

    dim = size(samp)

    if operator
        return LinearMapAA(forw!, back!,
            (sum(samp), prod(dim)), (name="fft",) ; T=T,
            idim = dim, # odim is 1D
            operator = operator,
        )
    end

    return LinearMapAA(
        (y,x) -> forw!(y, reshape(x,dim)), # fft(reshape(x,dim))[samp]
        (x,y) -> vec(back!(reshape(x,dim),y)), # vec(bfft(embed(y,samp)))
        (sum(samp), prod(dim)), (name="fft",) ; T=T,
    )
end


text:
using ZMQ
using JSON
using Zlib
using Nettle

ctx = Context(1)
s = Socket(ctx, ZMQ.REQ)
ZMQ.connect(s, "tcp://localhost:7001")
h = HashState(SHA1)
update!(h, "test")
jobhash = hexdigest!(h)
update!(h, "test2")
paramhash = hexdigest!(h)

fullhash = jobhash * "|" * paramhash
wiredata = bytestring(compress(json((Int=>Int)[1=>3]), 7))

#ZMQ.send(s, fullhash, SNDMORE)
#ZMQ.send(s, wiredata)
#ZMQ.recv(s)

#ZMQ.send(s, fullhash, SNDMORE)
#ZMQ.send(s, wiredata)
#ZMQ.recv(s)

#t1 = time()
for i=1:20
    tic()
    ZMQ.send(s, fullhash, SNDMORE)
    ZMQ.send(s, wiredata)
    msg = ZMQ.recv(s)
    #println(bytestring(msg))
    toc()
end
#println(time()-t1)
#println("^^ for 10000 messages")

max_stars_repo_name: binarybana/jobmon

max_issues_repo_name: binarybana/jobmon

max_forks_repo_name: binarybana/jobmon


text:
using FountainCodes

function test_encode_decode()
    for n in 1:100
        if graydecode(grayencode(n)) != n
            error("gray-code failure: graydecode(grayencode($n)) != $n")
        end
    end
    return true
end
@test test_encode_decode()

function test_nextgray()
    for hw in 1:10
        g = 0
        for n in 1:10
            g = nextgray(g, hw)
            if hamming_weight(g) != hw
                error("nextgray failure: $n-th number has weight $(hamming_weight(g)), but should have weight $hw")
            end
        end
    end
    return true
end
@test test_nextgray()

max_stars_repo_name: severinson/RaptorCodes

max_issues_repo_name: severinson/RaptorCodes

max_forks_repo_name: severinson/RaptorCodes


text:
# This file was generated, do not modify it. # hide
controls = [
    Step(1),              # to increment iteration parameter (`pipe.nrounds`)
    NumberSinceBest(4),   # main stopping criterion
    TimeLimit(2/3600),    # never train more than 2 sec
    InvalidValue()        # stop if NaN or ±Inf encountered
]

max_stars_repo_name: alan-turing-institute/MLJTutorials

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials


text:
using HorizonSideRobots
include("../basics.jl")

md"""
        # Найти площадь наибольшей прямоугольной перегородки
"""

"""Решение задачи 24"""
function solve_24!(r)

end


text:
include("/home/jefter66/Workspace/lablib/src/calcs/Uncertainty.jl")

function I1()
    Δr = 0.01
    r1 = 0.0335
    r2 = 0.06

    m = 2.2429
    Δm = 0.0001

    I1 = (1/2)*m*(r1^2 + r2^2)
    ΔI1 = (1/2)*((Δm*(r1^2 + r2^2) + 2*m*(r1*Δr + r2*Δr)))
    println(" Ig ± Δ Ig = ", I1, " ± ", ΔI1)
#    return [I1,ΔI1]
end

function I2(ω1, ω)
   Δω=1e-3
   I2 = I1()[1]*( ω1/(ω1 - ω))
   ΔI2   =  (I1()[2]  + (I1()[2]*ω +  Δω*I1()[1])/(ω^2))
    println("I2 ± Δ I2 = ", I2, " ± ", ΔI2)
    return [I2, ΔI2]
end

function equivalencia(v1,v2,σ1, σ2)
    abs(v1-v2) < 2*(σ1 + σ2)
end

function ΔE(ω1, ω,I1, I2)
    Δω = 0.001
    Ef = (1/2)*(I1[1] + I2[1])*ω^2
    Ei = (1/2)*(I1[1])*ω1^2

    δE = ((1/2)*(I1[2] + I2[2])*ω^2 + (I1[1] + I2[1])*ω*Δω) + (1/2)*(I1[2]*ω1^2 + 2*ω1*Δω*I1[1])

    println("E ± ΔE = ", Ef - Ei, " ± ", δE)
end

function results()
ω1 = [14.444, 22.765,24.353]
ω   = [7.453,12.777,13.629]
print("(1) : ")
I2_1 = I2(ω1[1], ω[1])
print("(2) : ")
I2_2 = I2(ω1[2],ω[2])
print("(3) : ")
I2_3 = I2(ω1[3],ω[3])

println("==============================================")

println("(1) & (2) : ", equivalencia(I2_1[1], I2_1[2], I2_2[1] , I2_2[2]))
println("(1) & (3) : ", equivalencia(I2_1[1],I2_1[2], I2_3[1],I2_3[2]))
println("(2) & (3) : ", equivalencia(I2_2[1],I2_2[2], I2_3[1], I2_3[2]))

println("==========================")

E1 = ΔE(ω1[1], ω[1], I1(), I2_1)
E2 = ΔE(ω1[2], ω[2], I1(), I2_2)
E3 = ΔE(ω1[3], ω[3], I1(), I2_3)
#println("ΔE1 ± δ E1 = ", E1[1], " ± ", E1[2])
#println("ΔE2 ± δ E2 =",  E2[1], " ± ", E2[2])
#println("ΔE3 ± δ E3 = ", E3[1] , " ± ", E3[2])

end

max_stars_repo_path: Fisica experimental/02/01/src/exp2.jl

max_issues_repo_path: Fisica experimental/02/01/src/exp2.jl

max_forks_repo_path: Fisica experimental/02/01/src/exp2.jl


text:
# prints the classification result produce by running mil_data.jl
# as an input, it takes the folder with saved outputs and possibly a number
# of an output file that only needs to be read
# run from the command line, e.g.
# julia path/to/file/plot_output.jl path/with/outputs path/to/save - this plots all the files in output path
# and saves them to given folder
# or
# julia path/to/file/plot_output.jl path/with/outputs path/to/save 1 3 7- this plots first, third and seventh file in the folder
# and saves them to given folder
using JLD
include("mil_util.jl")

output_path = ARGS[1]
if size(ARGS)[1] > 1
    save_path = ARGS[2]
    mkpath(save_path)
else
    save_path =""
end
files = readdir(output_path)
nfiles = size(files)[1]

if size(ARGS)[1] > 2
    inds = map(x->parse(Float64,x),ARGS[3:end])
else
    inds = 1:nfiles
end

for n in inds
    n = convert(Int, n)
    file = files[n]
    println("File n = $n")
    try
        data = load("$output_path/$file")
        plot_statistics(data, verb = false, save_path = save_path)
    catch e
        println("Something bad happened, probably wrong file type.")
    end
    n+=1
end

#merge all the output files?
#if save_path != ""
    #try
        #cd(save_path)
        #run(pipeline(`gs -sDEVICE=pdfwrite -dNOPAUSE -dBATCH -dSAFER -sOutputFile=out.pdf *.eps`))
        #println("All output files merged to out.pdf")
    #end
#end


text: const LIBC_RLIB = "/Users/ericdavies/.julia/v0.5/Rust/deps/libcshell/target/release/deps/liblibc-d9dba2869ce64308.rlib"


text:
function test_spec(f, value, args...)
    spec = f(value)
    @test spec == typeof(spec)(value, args...)
end

test_handle(args...) = test_spec(handle_by_name, args...)
test_bitmask(args...) = test_spec(bitmask_by_name, args...)
test_enum(args...) = test_spec(enum_by_name, args...)
test_constant(args...) = test_spec(constant_by_name, args...)
test_struct(name, type, is_returnedonly, extends, members) = test_spec(
    struct_by_name,
    name,
    type,
    is_returnedonly,
    extends,
    StructVector(map(x -> SpecStructMember(name, x...), members)),
)
test_union(name, types, selectors, is_returnedonly) = test_spec(union_by_name, name, types, selectors, is_returnedonly)
test_func(name, type, return_type, render_pass_compatibility, queue_compatibility, params, success_codes, error_codes) =
    test_spec(
        func_by_name,
        name,
        type,
        return_type,
        render_pass_compatibility,
        queue_compatibility,
        StructVector(map(x -> SpecFuncParam(name, x...), params)),
        success_codes,
        error_codes,
    )

function test_create_func(name, handle, _struct, index, batch)
    cf = create_func(name)
    f = func_by_name(name)
    @test cf == typeof(cf)(f, handle_by_name(handle), struct_by_name(_struct), children(f)[index], batch)
end

function test_destroy_func(name, handle, index, batch)
    df = destroy_func(name)
    f = func_by_name(name)
    @test df == typeof(df)(f, handle_by_name(handle), children(f)[index], batch)
end

@testset "Vulkan Specification" begin
    @testset "Handles" begin
        @testset "Definition" begin
            test_handle(:VkInstance, nothing, true)
            test_handle(:VkDevice, :VkPhysicalDevice, true)
            test_handle(:VkDeviceMemory, :VkDevice, false)
            @test isnothing(handle_by_name(:VkDescriptorUpdateTemplateKHR))
            test_handle(:VkDescriptorUpdateTemplate, :VkDevice, false)
        end

        @testset "Creation" begin
            test_create_func(
                :vkCreateInstance,
                :VkInstance,
                :VkInstanceCreateInfo,
                1,
                false,
            )
            test_create_func(
                :vkCreateGraphicsPipelines,
                :VkPipeline,
                :VkGraphicsPipelineCreateInfo,
                4,
                true,
            )
            test_create_func(
                :vkAllocateCommandBuffers,
                :VkCommandBuffer,
                :VkCommandBufferAllocateInfo,
                2,
                true,
            )
            test_create_func(
                :vkAllocateMemory,
                :VkDeviceMemory,
                :VkMemoryAllocateInfo,
                2,
                false,
            )
        end

        @testset "Destruction" begin
            test_destroy_func(:vkDestroyInstance, :VkInstance, 1, false)
            test_destroy_func(:vkDestroyDevice, :VkDevice, 1, false)
            test_destroy_func(:vkFreeDescriptorSets, :VkDescriptorSet, 4, true)
            test_destroy_func(:vkDestroyPipeline, :VkPipeline, 2, false)
            test_destroy_func(:vkFreeMemory, :VkDeviceMemory, 2, false)
        end
    end

    @testset "Defaults" begin
        @test default(handle_by_name(:VkInstance)) == :C_NULL
        @test all(
            default.(struct_by_name(:VkInstanceCreateInfo).members) .==
            [0, :C_NULL, 0, :C_NULL, 0, :C_NULL, 0, :C_NULL],
        )
    end

    @testset "Constants" begin
        test_constant(:VkBool32, :UInt32)
        test_constant(:VkFlags, :UInt32)
        test_constant(:VkPipelineCreateFlags, :VkFlags)

        @test follow_constant(constant_by_name(:VkPipelineCreateFlags)) ==
              follow_constant(:VkPipelineCreateFlags) ==
              :UInt32
        @test follow_constant(:UInt32) == :UInt32
        @test follow_constant(:(Ptr{UInt32})) == :(Ptr{UInt32})
    end

    @testset "Structs" begin
        test_struct(
            :VkApplicationInfo,
            GENERIC_INFO,
            false,
            [],
            [
                (:sType, :VkStructureType, false, false, REQUIRED, nothing, [], true),
                (:pNext, :(Ptr{Cvoid}), true, false, OPTIONAL, nothing, [], true),
                (:pApplicationName, :Cstring, true, false, OPTIONAL, nothing, [], true),
                (:applicationVersion, :UInt32, false, false, REQUIRED, nothing, [], true),
                (:pEngineName, :Cstring, true, false, OPTIONAL, nothing, [], true),
                (:engineVersion, :UInt32, false, false, REQUIRED, nothing, [], true),
                (:apiVersion, :UInt32, false, false, REQUIRED, nothing, [], true),
            ],
        )
        test_struct(
            :VkPipelineViewportCoarseSampleOrderStateCreateInfoNV,
            CREATE_INFO,
            false,
            [:VkPipelineViewportStateCreateInfo],
            [
                (:sType, :VkStructureType, false, false, REQUIRED, nothing, [], true),
                (:pNext, :(Ptr{Cvoid}), true, false, OPTIONAL, nothing, [], true),
                (:sampleOrderType, :VkCoarseSampleOrderTypeNV, false, false, REQUIRED, nothing, [], true),
                (:customSampleOrderCount, :UInt32, false, false, OPTIONAL, nothing, [:pCustomSampleOrders], true),
                (
                    :pCustomSampleOrders,
                    :(Ptr{VkCoarseSampleOrderCustomNV}),
                    true,
                    false,
                    REQUIRED,
                    :customSampleOrderCount,
                    [],
                    true,
                ),
            ],
        )
        test_struct(
            :VkPhysicalDeviceShadingRateImageFeaturesNV,
            DATA,
            false,
            [:VkPhysicalDeviceFeatures2, :VkDeviceCreateInfo],
            [
                (:sType, :VkStructureType, false, false, REQUIRED, nothing, [], true),
                (:pNext, :(Ptr{Cvoid}), false, false, OPTIONAL, nothing, [], true),
                (:shadingRateImage, :VkBool32, false, false, REQUIRED, nothing, [], true),
                (:shadingRateCoarseSampleOrder, :VkBool32, false, false, REQUIRED, nothing, [], true),
            ],
        )
        test_struct(
            :VkDescriptorSetLayoutBindingFlagsCreateInfo,
            CREATE_INFO,
            false,
            [:VkDescriptorSetLayoutCreateInfo],
            [
                (:sType, :VkStructureType, false, false, REQUIRED, nothing, [], true),
                (:pNext, :(Ptr{Cvoid}), true, false, OPTIONAL, nothing, [], true),
                (:bindingCount, :UInt32, false, false, OPTIONAL, nothing, [:pBindingFlags], true),
                (:pBindingFlags, :(Ptr{VkDescriptorBindingFlags}), true, false, POINTER_REQUIRED, :bindingCount, [], true),
            ],
        ),
        test_struct(
            :VkInstanceCreateInfo,
            CREATE_INFO,
            false,
            [],
            [
                (:sType, :VkStructureType, false, false, REQUIRED, nothing, [], true),
                (:pNext, :(Ptr{Cvoid}), true, false, OPTIONAL, nothing, [], true),
                (:flags, :VkInstanceCreateFlags, false, false, OPTIONAL, nothing, [], true),
                (:pApplicationInfo, :(Ptr{VkApplicationInfo}), true, false, OPTIONAL, nothing, [], true),
                (:enabledLayerCount, :UInt32, false, false, OPTIONAL, nothing, [:ppEnabledLayerNames], true),
                (:ppEnabledLayerNames, :(Ptr{Cstring}), true, false, REQUIRED, :enabledLayerCount, [], true),
                (:enabledExtensionCount, :UInt32, false, false, OPTIONAL, nothing, [:ppEnabledExtensionNames], true),
                (:ppEnabledExtensionNames, :(Ptr{Cstring}), true, false, REQUIRED, :enabledExtensionCount, [], true),
            ],
        )
        test_struct(
            :VkDisplayPlaneInfo2KHR,
            GENERIC_INFO,
            false,
            [],
            [
                (:sType, :VkStructureType, false, false, REQUIRED, nothing, [], true),
                (:pNext, :(Ptr{Cvoid}), true, false, OPTIONAL, nothing, [], true),
                (:mode, :VkDisplayModeKHR, false, true, REQUIRED, nothing, [], true),
                (:planeIndex, :UInt32, false, false, REQUIRED, nothing, [], true),
            ],
        )
        test_struct(
            :VkTransformMatrixKHR,
            DATA,
            false,
            [],
            [(:matrix, :(NTuple{3,NTuple{4,Float32}}), false, false, REQUIRED, nothing, [], true)],
        )
    end

    @testset "Unions" begin
        test_union(:VkClearColorValue, [:(NTuple{4,Float32}), :(NTuple{4,Int32}), :(NTuple{4,UInt32})], [], false)

        test_union(
            :VkPerformanceValueDataINTEL,
            [:UInt32, :UInt64, :Float32, :VkBool32, :Cstring],
            [
                :VK_PERFORMANCE_VALUE_TYPE_UINT32_INTEL,
                :VK_PERFORMANCE_VALUE_TYPE_UINT64_INTEL,
                :VK_PERFORMANCE_VALUE_TYPE_FLOAT_INTEL,
                :VK_PERFORMANCE_VALUE_TYPE_BOOL_INTEL,
                :VK_PERFORMANCE_VALUE_TYPE_STRING_INTEL,
            ],
            false,
        )
    end

    @testset "Functions" begin
        test_func(
            :vkCreateInstance,
            CREATE,
            :VkResult,
            [],
            [],
            [
                (:pCreateInfo, :(Ptr{VkInstanceCreateInfo}), true, false, REQUIRED, nothing, [], true),
                (:pAllocator, :(Ptr{VkAllocationCallbacks}), true, false, OPTIONAL, nothing, [], true),
                (:pInstance, :(Ptr{VkInstance}), false, false, REQUIRED, nothing, [], true),
            ],
            [:VK_SUCCESS],
            [
                :VK_ERROR_OUT_OF_HOST_MEMORY,
                :VK_ERROR_OUT_OF_DEVICE_MEMORY,
                :VK_ERROR_INITIALIZATION_FAILED,
                :VK_ERROR_LAYER_NOT_PRESENT,
                :VK_ERROR_EXTENSION_NOT_PRESENT,
                :VK_ERROR_INCOMPATIBLE_DRIVER,
            ],
        )
        test_func(
            :vkCmdBindPipeline,
            COMMAND,
            :Cvoid,
            [RenderPassInside(), RenderPassOutside()],
            [QueueGraphics(), QueueCompute()],
            [
                (:commandBuffer, :VkCommandBuffer, false, true, REQUIRED, nothing, [], true),
                (:pipelineBindPoint, :VkPipelineBindPoint, false, false, REQUIRED, nothing, [], true),
                (:pipeline, :VkPipeline, false, false, REQUIRED, nothing, [], true),
            ],
            [],
            [],
        )
    end

    @testset "Aliases" begin
        @test isalias(:VkDescriptorUpdateTemplateKHR)
        @test hasalias(:VkDescriptorUpdateTemplate)
        @test isalias(:VkImageStencilUsageCreateInfoEXT)
        @test !isalias(:VkImageStencilUsageCreateInfo)

        @test follow_alias(:VkDescriptorUpdateTemplateKHR) == :VkDescriptorUpdateTemplate
        @test follow_alias(:VkPhysicalDeviceMemoryProperties2KHR) == :VkPhysicalDeviceMemoryProperties2
        @test follow_alias(:VkPhysicalDeviceMemoryProperties2) == :VkPhysicalDeviceMemoryProperties2
    end

    @testset "Miscellaneous" begin
        @test length_chain(children(func_by_name(:vkAllocateDescriptorSets))[3], "pAllocateInfo->descriptorSetCount") ==
              [
            children(func_by_name(:vkAllocateDescriptorSets))[2],
            children(struct_by_name(:VkDescriptorSetAllocateInfo))[4],
        ]
    end
end


text:
abstract AbstractMeta
abstract AbstractΨMeta <: AbstractMeta

"""A `ψMeta` object contains all the meta-information to identify
an oriented wavelet within a filter bank. Fields:
* `γ` log-scale. `2^(-γ)` is proportional to center frequency
* `θ` orientation, i.e. angle (in 2d) or sign of center frequency (in 1d)
* `χ` chroma. `χ` is `mod(γ, n_filters_per_octave)`
* `j` octave. `j` is `div(γ, n_filters_per_octave)`
* `bandwidth` ∈]0,1] the width at -3dB, expressed in fraction of signal length
* `centerfrequency` ∈]0,1] is expressed in fraction of signal length
* `qualityfactor` ∈[1,max_qualityfactor] is equal to `centerfrequency/bandwidth`
* `scale` is the FWTM (full width at tenth maximum) in spatial domain"""
immutable ΨMeta1D <: AbstractΨMeta
    γ::Int16
    θ::Int8
    χ::Int8
    bandwidth::Float64
    centerfrequency::Float64
    j::Int8
    qualityfactor::Float64
    scale::Float64
end

immutable ΨMeta2D <: AbstractΨMeta
    γ::Int16
    θ::Int8
    χ::Int8
    aspectratio::Float64
    bandwidth::Float64
    centerfrequency::Float64
    j::Int8
    qualityfactor::Float64
    scale::Float64
end

"""A `ΦMeta` object contains all the meta-information to identify a low-pass
filter within a wavelet filter bank. Fields:
* `bandwidth` ∈]0,1] the width at -3dB, expressed in fraction of signal length
* `scale` ∈]0,1] is the FWTM (full width at tenth maximum) in spatial domain."""
immutable ΦMeta <: AbstractMeta
    bandwidth::Float64
    scale::Float64
end

get_γ(meta::AbstractΨMeta) = meta.γ
get_θ(meta::AbstractΨMeta) = meta.θ
get_χ(meta::AbstractΨMeta) = meta.χ
get_aspectratio(meta::ΨMeta2D) = meta.aspectratio
get_bandwidth(meta::AbstractMeta) = meta.bandwidth
get_centerfrequency(meta::AbstractΨMeta) = meta.centerfrequency
get_j(meta::AbstractΨMeta) = meta.j
get_qualityfactor(meta::AbstractΨMeta) = meta.qualityfactor
get_scale(meta::AbstractMeta) = meta.scale


text:
# Define plotting recipe
@recipe function f(s::SynthControlModel; kind = "overall")
    legend --> :topleft

    if kind == "weights"

        wp = sort(DataFrame(comp = s.treatment_panel.comp_labels,
                  weight = s.w), :weight, rev = true)

        @series begin
            seriestype := :bar
            label --> ""
            xguide --> "Donor"
            yguide --> "Weight"
            linecolor := "white"
            wp[wp.weight .> 0.05, :comp], wp[wp.weight .> 0.05, :weight]
        end

    elseif kind == "overall"

        @series begin
            label --> s.treatment_panel.treatment[1]
            sort!(unique(s.treatment_panel.data[!, s.treatment_panel.t_var])), [s.treatment_panel.y₁₀; s.treatment_panel.y₁₁]
        end

        @series begin
            label --> "Control"
            sort!(unique(s.treatment_panel.data[!, s.treatment_panel.t_var])), s.ŷ₁
        end

        @series begin
            label --> "Impact"
            seriestype --> :bar
            seriesalpha --> 0.5
            linecolor --> "white"
            sort!(unique(s.treatment_panel.data[!, s.treatment_panel.t_var]))[findfirst(x -> x >= s.treatment_panel.treatment[2],
                                    sort!(unique(s.treatment_panel.data[!, s.treatment_panel.t_var]))):end], s.τ̂
        end

        @series begin
            label --> "Intervention"
            seriestype := :vline
            linestyle := :dash
            seriescolor := "black"
            xguide := "Time"
            yguide := "Outcome"
            [s.treatment_panel.treatment[2]]
        end

    elseif kind == "diffplot"
        @series begin
            label --> s.treatment_panel.treatment[1]
            sort!(unique(s.treatment_panel.data[!, s.treatment_panel.t_var])), [s.treatment_panel.y₁₀; s.treatment_panel.y₁₁] .- s.ŷ₁
        end

        @series begin
            label --> ""
            seriestype := :scatter
            seriescolor := 1
            markerstrokecolor := "white"
            sort!(unique(s.treatment_panel.data[!, s.treatment_panel.t_var])), [s.treatment_panel.y₁₀; s.treatment_panel.y₁₁] .- s.ŷ₁
        end

        @series begin
            label --> "Intervention"
            seriestype := :vline
            linestyle := :dash
            seriescolor := "black"
            xguide := "Time"
            yguide := "Outcome"
            [s.treatment_panel.treatment[2]]
        end
    end

     nothing
end


text:
using Distributions
using LinearAlgebra
using Random
using Test

using EnsembleKalmanProcesses.EnsembleKalmanProcessModule
using EnsembleKalmanProcesses.ParameterDistributionStorage
@testset "EnsembleKalmanProcessModule" begin

    # Seed for pseudo-random number generator
    rng_seed = 41
    Random.seed!(rng_seed)

    ### Generate data from a linear model: a regression problem with n_par parameters
    ### and n_obs. G(u) = A \times u, where A : R^n_par -> R
    n_obs = 10                  # Number of synthetic observations from G(u)
    n_par =  2                  # Number of parameteres
    u_star = [-1., 2.]          # True parameters
    noise_level = .1            # Defining the observation noise level
    Γy = noise_level * Matrix(I, n_obs, n_obs) # Independent noise for synthetic observations
    noise = MvNormal(zeros(n_obs), Γy)

    C = [1 -.9; -.9 1]          # Correlation structure for linear operator
    A = rand(MvNormal(zeros(2,), C), n_obs)'    # Linear operator in R^{n_obs \times n_par}

    @test size(A) == (n_obs, n_par)

    y_star = A * u_star
    y_obs  = y_star + rand(noise)

    @test size(y_star) == (n_obs,)

    #### Define linear model
    function G(u)
        A * u
    end

    @test norm(y_star - G(u_star)) < n_obs * noise_level^2

    #### Define prior information on parameters
    prior_distns = [Parameterized(Normal(1., sqrt(2))),
              Parameterized(Normal(1., sqrt(2)))]
    constraints = [[no_constraint()], [no_constraint()]]
    prior_names = ["u1","u2"]
    prior = ParameterDistribution(prior_distns, constraints, prior_names)

    prior_mean = get_mean(prior)

    # Assuming independence of u1 and u2
    prior_cov = get_cov(prior)#convert(Array, Diagonal([sqrt(2.), sqrt(2.)]))

    ###
    ###  Calibrate (1): Ensemble Kalman Inversion
    ###

    N_ens = 50 # number of ensemble members
    N_iter = 20 # number of EKI iterations
    initial_ensemble = EnsembleKalmanProcessModule.construct_initial_ensemble(prior,N_ens;
                                                    rng_seed=rng_seed)
    @test size(initial_ensemble) == (n_par, N_ens)

    ekiobj = EnsembleKalmanProcessModule.EnsembleKalmanProcess(initial_ensemble,
                       y_obs, Γy, Inversion())

    # Find EKI timestep
    g_ens = G(get_u_final(ekiobj))
    @test size(g_ens) == (n_obs, N_ens)
    # as the columns of g are the data, this should throw an error
    g_ens_t = permutedims(g_ens, (2,1))
    @test_throws DimensionMismatch find_ekp_stepsize(ekiobj, g_ens_t)
    Δ = find_ekp_stepsize(ekiobj, g_ens)
    @test Δ ≈ 0.0625

    # EKI iterations
    params_i_vec = []
    g_ens_vec = []
    for i in 1:N_iter
        params_i = get_u_final(ekiobj)
        push!(params_i_vec,params_i)
        g_ens = G(params_i)
        push!(g_ens_vec, g_ens)
        if i == 1
            g_ens_t = permutedims(g_ens, (2,1))
            @test_throws DimensionMismatch EnsembleKalmanProcessModule.update_ensemble!(ekiobj, g_ens_t)
        end
        EnsembleKalmanProcessModule.update_ensemble!(ekiobj, g_ens)
    end
    push!(params_i_vec, get_u_final(ekiobj))

    @test get_u_prior(ekiobj) == params_i_vec[1]
    @test get_u(ekiobj) == params_i_vec
    @test get_g(ekiobj) == g_ens_vec
    @test get_g_final(ekiobj) == g_ens_vec[end]
    @test get_error(ekiobj) == ekiobj.err

    # EKI results: Test if ensemble has collapsed toward the true parameter
    # values
    eki_final_result = vec(mean(get_u_final(ekiobj), dims=2))
    # @test norm(u_star - eki_final_result) < 0.5

    # Plot evolution of the EKI particles
    eki_final_result = vec(mean(get_u_final(ekiobj), dims=2))

    if TEST_PLOT_OUTPUT
        gr()
        p = plot(get_u_prior(ekiobj)[1,:], get_u_prior(ekiobj)[2,:], seriestype=:scatter)
        plot!(get_u_final(ekiobj)[1, :],  get_u_final(ekiobj)[2,:], seriestype=:scatter)
        plot!([u_star[1]], xaxis="u1", yaxis="u2", seriestype="vline",
            linestyle=:dash, linecolor=:red)
        plot!([u_star[2]], seriestype="hline", linestyle=:dash, linecolor=:red)
        savefig(p, "EKI_test.png")
    end

    ###
    ###  Calibrate (2): Ensemble Kalman Sampler
    ###
    eksobj = EnsembleKalmanProcessModule.EnsembleKalmanProcess(initial_ensemble,
                      y_obs, Γy,
                      Sampler(prior_mean, prior_cov))

    # EKS iterations
    for i in 1:N_iter
        params_i = get_u_final(eksobj)
        g_ens = G(params_i)
        if i == 1
            g_ens_t = permutedims(g_ens, (2,1))
            @test_throws DimensionMismatch EnsembleKalmanProcessModule.update_ensemble!(eksobj, g_ens_t)
        end

        EnsembleKalmanProcessModule.update_ensemble!(eksobj, g_ens)
    end

    # Plot evolution of the EKS particles
    eks_final_result = vec(mean(get_u_final(eksobj), dims=2))

    if TEST_PLOT_OUTPUT
        gr()
        p = plot(get_u_prior(eksobj)[1,:], get_u_prior(eksobj)[2,:], seriestype=:scatter)
        plot!(get_u_final(eksobj)[1, :],  get_u_final(eksobj)[2,:], seriestype=:scatter)
        plot!([u_star[1]], xaxis="u1", yaxis="u2", seriestype="vline",
            linestyle=:dash, linecolor=:red)
        plot!([u_star[2]], seriestype="hline", linestyle=:dash, linecolor=:red)
        savefig(p, "EKS_test.png")
    end

    posterior_cov_inv = (A' * (Γy\A) + 1 * Matrix(I, n_par, n_par)/prior_cov)
    ols_mean          = (A'*(Γy\A)) \ (A'*(Γy\y_obs))
    posterior_mean    = posterior_cov_inv \ ( (A'*(Γy\A)) * ols_mean  + (prior_cov\prior_mean))

    #### This tests correspond to:
    # EKI provides a solution closer to the ordinary Least Squares estimate
    @test norm(ols_mean - eki_final_result) < norm(ols_mean - eks_final_result)
    # EKS provides a solution closer to the posterior mean
    @test norm(posterior_mean - eks_final_result) < norm(posterior_mean - eki_final_result)

    ##### I expect this test to make sense:
    # In words: the ensemble covariance is still a bit ill-dispersed since the
    # algorithm employed still does not include the correction term for finite-sized
    # ensembles.
    @test abs(sum(diag(posterior_cov_inv\cov(get_u_final(eksobj),dims=2))) - n_par) > 1e-5

    ###
    ###  Calibrate (3): Unscented Kalman Inversion
    ###

    N_iter = 20 # number of UKI iterations
    α_reg =  1.0
    update_freq = 0
    process = Unscented(prior_mean, prior_cov, length(y_star),  α_reg, update_freq)
    ukiobj = EnsembleKalmanProcessModule.EnsembleKalmanProcess(y_star, Γy, process)

    # UKI iterations
    params_i_vec = []
    g_ens_vec = []
    for i in 1:N_iter
        params_i = get_u_final(ukiobj)
        push!(params_i_vec,params_i)
        g_ens = G(params_i)
        push!(g_ens_vec, g_ens)
        if i == 1
            g_ens_t = permutedims(g_ens, (2,1))
            @test_throws DimensionMismatch EnsembleKalmanProcessModule.update_ensemble!(ukiobj, g_ens_t)
        end
        EnsembleKalmanProcessModule.update_ensemble!(ukiobj, g_ens)
    end
    push!(params_i_vec,get_u_final(ukiobj))

    @test get_u_prior(ukiobj) == params_i_vec[1]
    @test get_u(ukiobj) == params_i_vec
    @test get_g(ukiobj) == g_ens_vec
    @test get_g_final(ukiobj) == g_ens_vec[end]
    @test get_error(ukiobj) == ukiobj.err

    # UKI results: Test if ensemble has collapsed toward the true parameter
    # values
    uki_final_result = get_u_mean_final(ukiobj)
    @test norm(u_star - uki_final_result) < 0.5

    if TEST_PLOT_OUTPUT
        gr()
        θ_mean_arr = hcat(ukiobj.process.u_mean...)
        N_θ = length(ukiobj.process.u_mean[1])
        θθ_std_arr = zeros(Float64, (N_θ, N_iter+1))
        for i = 1:N_iter+1
            for j = 1:N_θ
                θθ_std_arr[j, i] = sqrt(ukiobj.process.uu_cov[i][j,j])
            end
        end

        ites = Array(LinRange(1, N_iter+1, N_iter+1))
        p = plot(ites,grid=false, θ_mean_arr[1,:], yerror=3.0*θθ_std_arr[1,:],  label="u1")
        plot!(ites, fill(u_star[1], N_iter+1), linestyle=:dash, linecolor=:grey,label=nothing)
        plot!(ites,grid=false, θ_mean_arr[2,:], yerror=3.0*θθ_std_arr[2,:], label="u2", xaxis="Iterations")
        plot!(ites, fill(u_star[2], N_iter+1), linestyle=:dash, linecolor=:grey,label=nothing)
        savefig(p, "UKI_test.png")
    end

end


text:
# This file was generated, do not modify it. # hide
using RDatasets, MLJ
iris = dataset("datasets", "iris")

first(iris, 3) |> pretty


text:
"""
    nt(s::Symbol, v = 0.0)

Creates a named tuple `(s=v,)` where `s` is a provided symbol, and `v` is the value.
"""
nt(s::Symbol, v = 0.0) = NamedTuple{(s,)}([v])

inrange(x,r) = r[1] ≤ x ≤ r[2]


text:
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:light
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.4'
#       jupytext_version: 1.2.1
#   kernelspec:
#     display_name: Julia 1.1.0
#     language: julia
#     name: julia-1.1
# ---

# # Ocean transports and maps
#
# Ocean transport is represented as a vector field defined over the Global Ocean. These case be integrated along a grid path to compute transports between any two points.
#
# For more about how these methods, please refer to **Forget et al, 2015.** ECCO version 4: An integrated framework for non-linear inverse modeling and global ocean state estimation. https://doi.org/10.5194/gmd-8-3071-2015
#
# Key functions:
# - `LatitudeCircles` computes integration paths that follow latitude circles
# - `ThroughFlow` computes transports through these integration paths

# +
using MeshArrays, Plots
include(joinpath(dirname(pathof(MeshArrays)),"Plots.jl"))
include("prepare_transports.jl")

if !isdir("GRID_LLC90")
    run(`git clone https://github.com/gaelforget/GRID_LLC90`)
end

mygrid=GridSpec("LLC90");
GridVariables=GridLoad(mygrid);
# -

# ### Time average and vertically integrate transports
#
# _Note: these intermediate results have pre-computed for you._

# +
#using Statistics
#using FortranFiles
#!isdir("nctiles_climatology") ? error("missing files") : nothing
#include(joinpath(dirname(pathof(MeshArrays)),"gcmfaces_nctiles.jl"))
#(TrspX, TrspY, TauX, TauY, SSH)=trsp_prep(mygrid,GridVariables,"GRID_LLC90/");
# -

# Let's read the intermediate results from file using `trsp_read`.

(TrspX, TrspY, TauX, TauY, SSH)=trsp_read(mygrid,"GRID_LLC90/");

# Alternatively these can be recomputed from three-dimensional time-varying velocity fields via `trsp_prep`.

# ### Transports between latitude bands

# +
UVmean=Dict("U"=>TrspX,"V"=>TrspY,"dimensions"=>["x","y"]);
LC=LatitudeCircles(-89.0:89.0,GridVariables);

T=Array{Float64,1}(undef,length(LC));
for i=1:length(LC)
   T[i]=ThroughFlow(UVmean,LC[i],GridVariables)
end
# -

# ### Plot result

lat=-89.0:89.0
plot(lat,T/1e6,xlabel="latitude",ylabel="Sverdrup (10^6 m^3 s^-1)",
    label="ECCOv4r2",title="Northward transport of seawater (Global Ocean)")

# ### Plot transport maps
#
# _Note that vector field orientations differ amongst the arrays._

heatmap(1e-6*TrspX,clims=(-20.0,20.0))
#heatmap(1e-6*TrspY,clims=(-20.0,20.0))

# Convert to `Sv` and mask out land

u=1e-6 .*UVmean["U"]; v=1e-6 .*UVmean["V"];
u[findall(GridVariables["hFacW"][:,1].==0)].=NaN;
v[findall(GridVariables["hFacS"][:,1].==0)].=NaN;

# `x/y` transport at cell center

# +
using Statistics
nanmean(x) = mean(filter(!isnan,x))
nanmean(x,y) = mapslices(nanmean,x,dims=y)

(u,v)=exch_UV(u,v);
uC=similar(u); vC=similar(v);
for iF=1:mygrid.nFaces;
    tmp1=u[iF][1:end-1,:]; tmp2=u[iF][2:end,:];
    uC[iF]=reshape(nanmean([tmp1[:] tmp2[:]],2),size(tmp1));
    tmp1=v[iF][:,1:end-1]; tmp2=v[iF][:,2:end];
    vC[iF]=reshape(nanmean([tmp1[:] tmp2[:]],2),size(tmp1));
end
# -

# `Eastward/Northward` transport
#
# _Note: compare vector field orientations with previous plot._

cs=GridVariables["AngleCS"];
sn=GridVariables["AngleSN"];
u=uC.*cs-vC.*sn;
v=uC.*sn+vC.*cs;
heatmap(u,clims=(-20.0,20.0))

# Transport as a function of longitude and latitude

# +
#cbp="https://github.com/gaelforget/CbiomesProcessing.jl"
#using Pkg; Pkg.add(PackageSpec(url=cbp))
using CbiomesProcessing

SPM,lon,lat=CbiomesProcessing.read_SPM("GRID_LLC90/")
uI=CbiomesProcessing.MatrixInterp(write(u),SPM,size(lon))
heatmap(vec(lon[:,1]),vec(lat[1,:]),transpose(uI))
# -

max_stars_repo_name: lmilechin/MeshArrayNotebooks

max_issues_repo_name: lmilechin/MeshArrayNotebooks

max_forks_repo_name: lmilechin/MeshArrayNotebooks


text:
######### Stan program example  ###########

using AxisKeys, DataFrames
using StanSample
using Test

binorm_model = "
  transformed data {
      matrix[2,2] Sigma;
      vector[2] mu;

      mu[1] <- 0.0;
      mu[2] <- 0.0;
      Sigma[1,1] <- 1.0;
      Sigma[2,2] <- 1.0;
      Sigma[1,2] <- 0.10;
      Sigma[2,1] <- 0.10;
  }
  parameters {
      vector[2] y;
  }
  model {
        y ~ multi_normal(mu,Sigma);
  }
  "

sm = SampleModel("binormal", binorm_model);

rc = stan_sample(sm)

if success(rc)
  samples = read_samples(sm)

  # Fetch cmdstan summary data frame
  df = read_summary(sm)
  @test df[df.parameters .== Symbol("y[1]"), :mean][1] ≈ 0.0 atol=2.0
  @test df[df.parameters .== Symbol("y[2]"), :mean][1] ≈ 0.0 atol=2.0

end


text:
module SIGPROC

include("read.jl")
include("fake.jl")
include("plotting.jl")

end


text:
using Test
using BinaryBuilderBase: extract_kwargs, extract_fields, strip_backslash, ArchiveSource

@testset "Compat functions" begin
    foo(; kwargs...) = collect(extract_kwargs(kwargs, (:bar, :qux)))
    @test foo(; a = 1) == Pair[]
    @test foo(; bar = 2, b = 3) == [:bar => 2]
    @test foo(; c = 4, qux = 5) == [:qux => 5]
    @test Set(foo(; qux = 6, d = 7, bar = 8)) == Set([:bar => 8, :qux => 6])

    @test extract_fields(ArchiveSource("http://example.org", "this is the hash"; unpack_target = "target")) == Dict("url" => "http://example.org", "hash" => "this is the hash", "unpack_target" => "target")

    @test strip_backslash("/home/wizard") == "/home/wizard"
    @test strip_backslash("/home/wizard/") == "/home/wizard"
    @test strip_backslash("/home//wizard///") == "/home//wizard"
    @test strip_backslash("wizard.jl") == "wizard.jl"
    # Test integration with `basename`
    @test basename(strip_backslash("/home/wizard")) == "wizard"
    @test basename(strip_backslash("/home/wizard/")) == "wizard"
    @test basename(strip_backslash("/home//wizard///")) == "wizard"
    @test basename(strip_backslash("wizard.jl")) == "wizard.jl"
end


text:
#!/usr/bin/env julia
using Base: LineEdit, REPL
using Base.LineEdit: Prompt

function stripblock(block::Vector{SubString{String}})
    length(block) > 0 && block[1] == "begin" || return block
    shift!(block)

    length(block) > 0 && block[end] == "end" || return :incomplete
    pop!(block)
    block
end

splitstring(x) = split(String(x), "\n"; keep=false)

function return_callback(leps::LineEdit.PromptState)
    # this could probably be more efficient
    tmp = splitstring(LineEdit.buffer(leps))
    stripblock(tmp) !== :incomplete
end

function RunShell(lang::AbstractString = "mdl")
    t = Base.Terminals.TTYTerminal("xterm", STDIN, STDOUT, STDERR)
    repl = REPL.LineEditREPL(t)

    replc = REPL.REPLCompletionProvider()
    render_prompt = Prompt(rpad(lang, 5) * "> ";
                           prompt_prefix = repl.hascolor ? Base.text_colors[:magenta] : "",
                           prompt_suffix = repl.hascolor ? Base.input_color : "",
                           keymap_func_data = repl,
                           on_enter = return_callback)

    # unified history
    hp = REPL.REPLHistoryProvider(Dict{Symbol,Any}(:render => render_prompt))

    REPL.history_reset_state(hp)
    render_prompt.hist = hp

    render_prompt.on_done = (s, buf, ok) -> begin
        ok || return transition(s, :abort)
        try
            # this is awfully complicated
            toparse = String.(take!(buf) |> splitstring |> stripblock)
            ps = Render.renderfile(toparse, lang)
            length(toparse) == 1 && Render.display(ps, `display`)
        catch err
            REPL.print_response(repl, err, catch_backtrace(), true, Base.have_color)
        end
        REPL.reset(repl)
        REPL.prepare_next(repl)
        REPL.reset_state(s)
        s.current_mode.sticky || transition(s, render_prompt)
    end

    search_prompt, skeymap = LineEdit.setup_search_keymap(hp)

    repl_transition(symb::Char, mode, s, o...) = begin
        if isempty(s) || position(LineEdit.buffer(s)) == 0
            buf = copy(LineEdit.buffer(s))
            LineEdit.transition(s, mode) do
                LineEdit.state(s, mode).input_buffer = buf
            end
        else
            edit_insert(s, symb)
        end
    end

    prefix_prompt, prefix_keymap = LineEdit.setup_prefix_keymap(hp, render_prompt)

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

    b = Dict{Any,Any}[skeymap, REPL.mode_keymap(render_prompt), prefix_keymap, LineEdit.history_keymap, LineEdit.default_keymap, LineEdit.escape_defaults]

    modal = LineEdit.ModalInterface([render_prompt, search_prompt, prefix_prompt])

    println("""
            Welcome to interactive mode!
            Every block is parsed as a separate file.
            The special keywords `begin` and `end` are available for declaring multi-line blocks.
            A line outside of a block will be run, and then displayed.
            """)

    REPL.run_interface(t, modal)
end
RunShell()

max_stars_repo_name: AnAverageHuman/jlgraphics

max_issues_repo_name: AnAverageHuman/jlgraphics

max_forks_repo_name: AnAverageHuman/jlgraphics


text:
using DrWatson
using DataFrames
using BSON
using CSV

raw = collect_results(datadir())

cols = eachcol(raw)

vals = map(values(cols)) do col
  map(col) do v
    if isa(v,NamedTuple)
      v.max
    elseif isa(v,Tuple)
      prod(v)
    else
      v
    end
  end
end

dfmax = DataFrame(vals,keys(cols))
sort!(dfmax,[:np,:ir])

rows = eachrow(copy(dfmax))
dict = Dict{Tuple{Int,Int}}{Any}()
for row in rows
  key = (row.ngdofs,row.np)
  if haskey(dict,key)
    prevrow = dict[key]
    for (k,v) in pairs(row)
      if v < prevrow[k]
        prevrow[k] = v
      end
    end
  else
    dict[key] = row
  end
end

df = DataFrame(collect(values(dict)))
sort!(df,[order(:ngdofs,rev=false),order(:np,rev=false)])

mkpath(plotsdir())
fn = plotsdir("summary.csv")
CSV.write(fn,df)
println(df)
df

max_stars_repo_name: santiagobadia/GridapGeosciences

max_issues_repo_name: santiagobadia/GridapGeosciences

max_forks_repo_name: santiagobadia/GridapGeosciences


text:
horizontal, vertical, aim = 0, 0, 0
for line in eachline()
    direction, size_string = split(line)
    size_int = parse(Int, size_string)
    if direction == "forward"
        global horizontal += size_int
        global vertical += aim * size_int
    elseif direction == "down"
        global aim += size_int
    elseif direction == "up"
        global aim -= size_int
    end
end
println(horizontal * aim)
println(horizontal * vertical)


text:
"""

```
forces!(x::Vector{T},f::Vector{T},opt::Options) where T
```

Subroutine that computes the force. It modifies the input `f` vector.

"""
function forces!(f::Vector{T},x::Vector{T},sys::System{T},opt::Options) where T
  @unpack n, sides = sys
  @unpack eps, sig = opt
  @. f = zero(T)
  for i in 1:n-1
    for j in i+1:n
      dx = image(x[j] - x[i], sides)
      r = norm(dx)
      drdx = -dx/r
      @fastpow dudx = 4*eps*(12*sig^12/r^13 - 6*sig^6/r^7)*drdx
      f[i] = f[i] + dudx
      f[j] = f[j] - dudx
    end
  end
  return f
end


text:
# Tests conversion among basis-structure representations

@testset "Test Basis Structure Representations" begin
    mb = Basis(SplineParams(linspace(-3, 3, 9), 0, 1),
               ChebParams(7, 1e-4, 10.0),
               LinParams(linspace(-4, 2, 11), 0))

    # construct evaluation points
    X, x123 = nodes(mb)

    # construct expanded, direct, and tensor basis-structure representation
    Φ_expanded = BasisMatrices.BasisMatrix(mb, BasisMatrices.Expanded(), X)
    Φ_direct = BasisMatrices.BasisMatrix(mb, BasisMatrices.Direct(), X)
    Φ_tensor = BasisMatrices.BasisMatrix(mb, BasisMatrices.Tensor(), x123)

    # construct expanded, direct, and tensor basis-structure representation with 1D basis
    Φ_expanded_1d = BasisMatrices.BasisMatrix(mb[1], BasisMatrices.Expanded(), X[:,1])
    Φ_direct_1d = BasisMatrices.BasisMatrix(mb[1], BasisMatrices.Direct(), X[:,1])

    @testset "test standard Base methods" begin
        # test == and ndims, multiD
        for Φ in (Φ_expanded, Φ_direct, Φ_tensor)
            @test  Φ  ==  Φ
            @test ndims(Φ) == 3
        end

        for Φ in (Φ_direct, Φ_tensor)
            @test ( Φ_expanded  ==  Φ ) == false
        end

        # test == and ndims, 1D
        for Φ in (Φ_expanded_1d, Φ_direct_1d)
            @test  Φ  ==  Φ
            @test ndims(Φ) == 1
        end

        @test  Φ_expanded_1d  ==  Φ_direct_1d

    end

    @testset "test convert methods" begin
        @test  Φ_direct  ==  convert(Direct, Φ_tensor)
        @test  Φ_expanded  ==  convert(Expanded, Φ_direct)
        @test  Φ_expanded  ==  convert(Expanded, Φ_tensor)
        @test ==(Φ_expanded,
                 convert(Expanded, convert(Direct, Φ_tensor)))

        # test the no-op covert method
        @test convert(Expanded, Φ_expanded) == Φ_expanded
        @test convert(Direct, Φ_direct) == Φ_direct
        @test convert(Tensor, Φ_tensor) == Φ_tensor

        # test that we can't do expanded -> (direct|tensor) or direct -> tensor
        @test_throws MethodError convert(Direct, Φ_expanded)
        @test_throws MethodError convert(Tensor, Φ_expanded)
        @test_throws MethodError convert(Tensor, Φ_direct)
    end

    @testset "test internal tools" begin
        ## test _vals_type
        for (TF, TM) in [(Spline, SparseMatrixCSC{Float64,Int}),
                         (Lin, SparseMatrixCSC{Float64,Int}),
                         (Cheb, Matrix{Float64})]
            @test BasisMatrices._vals_type(TF) == TM
            @test BasisMatrices._vals_type(TF()) == TM
        end

        @test BasisMatrices._vals_type(BasisMatrices.BasisFamily) == AbstractMatrix{Float64}

        ## test _checkx
        # create test data
        xm = rand(1, 2)
        xv = rand(2)
        xvv = [rand(2) for i=1:2]

        @test BasisMatrices._checkx(2, xm) == xm
        @test BasisMatrices._checkx(2, xv) == reshape(xv, 1, 2)
        @test BasisMatrices._checkx(2, xvv) == xvv
        @test BasisMatrices._checkx(1, xv) == xv

        @test_throws ErrorException BasisMatrices._checkx(2, rand(1, 3))
        @test_throws ErrorException BasisMatrices._checkx(2, rand(3))
        @test_throws ErrorException BasisMatrices._checkx(2, [rand(2) for i=1:3])

        ## test check_convert
        for Φ in (Φ_expanded, Φ_direct, Φ_tensor)
            @test BasisMatrices.check_convert(Φ, zeros(1, 3)) == (3, 1, 3)
            @test_throws ErrorException BasisMatrices.check_convert(Φ, zeros(1, 2))
            @test_throws ErrorException BasisMatrices.check_convert(Φ, fill(-1, 1, 3))
        end

        ## test check_basis_structure
        @testset "test check_basis_structure" begin
            order1 = zeros(Int, 1, 2)
            out1 = BasisMatrices.check_basis_structure(2, xm, order1)
            @test out1 == (1, order1, order1, [1 1], xm)

            # test when order::Int
            @test BasisMatrices.check_basis_structure(2, xm, 0) == out1

            # test the reshape(order, 1, N) branch (isa(order, Vector))
            @test BasisMatrices.check_basis_structure(2, xm, [0, 0]) == out1

            # check N=1 --> order = fill(order, 1, 1) branch
            out_1d = BasisMatrices.check_basis_structure(1, xv, 0)
            order_1d_out = fill(0, 1, 1)
            @test out_1d == (1, order_1d_out, order_1d_out, fill(1, 1, 1), xv)

            # check m > 1 branch
            order_m2 = [0 0; 1 1; 0 -1]
            out_m2 = BasisMatrices.check_basis_structure(2, xm, order_m2)
            @test out_m2 == (3, order_m2, [0 -1], [2 3], xm)

            order2 = zeros(Int, 1, 5)  # should throw error
            @test_throws ErrorException BasisMatrices.check_basis_structure(2, X,
                                                                       order2)
        end

    end

    @testset "constructors" begin
        nothing
    end

    @testset "Test from PR #25" begin
        basisμ = BasisMatrices.Basis(Cheb, 20, 0.0, 1.0)
        basisσ = BasisMatrices.Basis(Cheb, 20, 0.0, 1.0)
        basis = BasisMatrices.Basis(basisμ, basisσ)
        S, (μs, σs) = BasisMatrices.nodes(basis)
        bs = BasisMatrices.BasisMatrix(basis, BasisMatrices.Expanded(), S, [0 2])
        @test isa(bs, BasisMatrices.BasisMatrix{Expanded}) == true
    end

    # call show (which calls writemime) so we can get 100% coverage
    @testset "Printing" begin
        iob = IOBuffer()
        show(iob, Φ_tensor)
    end

end  # testset


text:
"""
    SelectTransform(dims)

Select the dimensions `dims` that the kernel is applied to.
```
    dims = [1,3,5,6,7]
    tr = SelectTransform(dims)
    X = rand(100,10)
    transform(tr,X,obsdim=2) == X[dims,:]
```
"""
struct SelectTransform{T} <: Transform
    select::T
end

set!(t::SelectTransform, dims) = t.select .= dims

duplicate(t::SelectTransform,θ) = t

(t::SelectTransform)(x::AbstractVector) = _maybe_unwrap(view(x, t.select))

_maybe_unwrap(x) = x
_maybe_unwrap(x::AbstractArray{<:Any, 0}) = x[]

_map(t::SelectTransform, x::ColVecs) = _wrap(view(x.X, t.select, :), ColVecs)
_map(t::SelectTransform, x::RowVecs) = _wrap(view(x.X, :, t.select), RowVecs)

_wrap(x::AbstractVector{<:Real}, ::Any) = x
_wrap(X::AbstractMatrix{<:Real}, ::Type{T}) where {T} = T(X)

Base.show(io::IO, t::SelectTransform) = print(io, "Select Transform (dims: ", t.select, ")")


text:
@testset "Length" begin
    for len in [1, 16, 32, 64]
        if len <= 32
            @test length(DNAMer(random_dna_kmer(len))) == len
            @test length(RNAMer(random_rna_kmer(len))) == len
        end
        @test length(BigDNAMer(random_dna_kmer(len))) == len
        @test length(BigRNAMer(random_rna_kmer(len))) == len
    end
end


text:
# debug info collection module
module DebugInfo

mutable struct StackInfo
  funcName::AbstractString
  filepath::AbstractString
  lineno::Int64
end

mutable struct VarInfo
  name::AbstractString
  vtype::AbstractString
  value::AbstractString
end

# Tree structure to log variable
Domain = Union{VarInfo, AbstractString}
mutable struct TreeNode
  # var value could be: VarInfo or Module name (at top level)
  content::Domain
  parent::Int64
  children::Array{Int64}
end
mutable struct Tree
  nodes::Array{TreeNode}
end
Tree() = Tree([])
function addRoot!(tree::Tree, node::TreeNode)
  push!(tree.nodes, node)
end
function addChild!(tree::Tree, id::Int64, var::Domain)
  1 <= id <= length(tree.nodes) || throw(BoundsError(tree, id))
  push!(tree.nodes, TreeNode(var, id, []))
  child = length(tree.nodes)
  push!(tree.nodes[id].children, child)
  return length(tree.nodes)
end
children(tree, id) = tree.nodes[id].children
parent(tree,id) = tree.nodes[id].parent

# Vars: Vars[ModuleName] = VarTree
Vars = Tree()
# stack frame
stacks = []

# should be called when hit a breakpoint
# collect stack info, debugger function will be
# deleted at getStackInfo()
function collectStackInfo()
  global stacks
  stacks = []
  for frame in stacktrace()
    funcName = String(frame.func)
    filepath = String(frame.file)
    lineno = frame.line
    stackInfo = StackInfo(funcName, filepath, lineno)
    push!(stacks, stackInfo)
  end
end

# return stack info list (ele: dict)
# Should be called when vs-code request stack
# debugger stack will be removed
function getStackInfo()
  global stacks
  results = []
  cnt = 0
  for frame in stacks
    debugger_file = r"(EventHandler\.jl)|(judy\.jl)|(MsgHandler\.jl)|(DebugInfo\.jl)|(RunTime\.jl)"
    m = match(debugger_file, frame.filepath)
    if !isa(m, RegexMatch)
      cnt += 1
      push!(results, Dict("frameId" => cnt,
                          "name" => frame.funcName,
                          "path" => frame.filepath,
                          "line" => frame.lineno))
    end
  end
  return results
end

# collect variable info
# currently we only collect Main module
function collectVarInfo()
  global Vars
  Vars = Tree()
  addRoot!(Vars, TreeNode("Main", 1, []))
  for var in names(Main)[5:end]
    # collect var name
    var_name = string(var)
    # collect var type
    type_ast = Meta.parse("typeof($(var))")
    vtype = Core.eval(Main, type_ast)
    var_type = string(vtype)
    if var_type == "DataType" || var_type == "Module"
      continue
    end
    # collect var value
    var_value = string(Core.eval(Main, var))
    # collect var ref
    var_info = VarInfo(var_name, var_type, var_value)
    ref = addChild!(Vars, 1, var_info)
    # decomposite var
    if length(fieldnames(vtype)) != 0
      decomposeVar(vtype, var, ref)
    end
  end
end

# decompose var and log them into Vars recursivly
function decomposeVar(vtype::DataType, var::Symbol, parent::Int64)
  global Vars
  for field in fieldnames(vtype)
    var_name = string(field)
    ast = Meta.parse("getfield($(var), :$(field))")
    var_value = string(Core.eval(Main, ast))
    ast = Meta.parse("typeof($(var).$(field))")
    vtype = Core.eval(Main, ast)
    var_type = string(vtype)
    var_info = VarInfo(var_name, var_type, var_value)
    ref = addChild!(Vars, parent, var_info)
    if length(fieldnames(vtype)) != 0
      devar = Core.eval(Main, Meta.parse("$(var).$(field)"))
      decomposeVar(vtype, Symbol(devar), ref)
    end
  end
end

function getVarInfo(ref)
  global Vars
  result = []
  var_refs = children(Vars, ref)
  for ref in var_refs
    var_ref = 0
    if length(children(Vars, ref)) != 0
      var_ref = ref
    end
    var = Vars.nodes[ref].content
    push!(result, Dict("name" => var.name,
                       "value" => var.value,
                       "type" => var.vtype,
                       "variablesReference" => var_ref))
  end
  return result
end

end # module DebugInfo

max_stars_repo_name: judy-vscode/judy

max_issues_repo_name: judy-vscode/judy

max_forks_repo_name: judy-vscode/judy


text:
@with_kw struct DeepPoly <: Solver
    max_iter::Int64 = 100
end

function init_deep_poly_symbolic_interval(domain)
    # VF = Vector{HalfSpace{Float64, Vector{Float64}}}
    # domain = HPolytope(VF(constraints_list(domain)))

    n = dim(domain)
    I = Matrix{Float64}(LinearAlgebra.I(n))
    Z = zeros(n) # for the constant terms
    return SymbolicInterval([I Z], [I Z], domain)
end

struct LayerNegPos{F<:ActivationFunction, N<:Number}
    weights::Matrix{N}
    bias::Vector{N}
    W_neg::Matrix{N}
    W_pos::Matrix{N}
    activation::F
end

LayerNegPos(L::Layer) = LayerNegPos(L.weights, L.bias, min.(L.weights, 0),
                                    max.(L.weights, 0), L.activation)
Layer(L::LayerNegPos) = Layer(L.weights, L.bias, L.activation)

n_nodes(L::LayerNegPos) = length(L.bias)
affine_map(L::LayerNegPos, x) = L.weights*x + L.bias

struct NetworkNegPos <: AbstractNetwork
    layers::Vector{LayerNegPos}
end

NetworkNegPos(net::Network) = NetworkNegPos([LayerNegPos(l) for l in net.layers])

function compute_output(nnet::NetworkNegPos, input)
    curr_value = input
    for layer in nnet.layers
        curr_value = layer.activation(affine_map(layer, curr_value))
    end
    return curr_value
end

function interval_map(W_neg::AbstractMatrix{N}, W_pos::AbstractMatrix{N},
                      l::AbstractVecOrMat, u::AbstractVecOrMat) where N
    l_new = W_pos * l .+ W_neg * u
    u_new = W_pos * u .+ W_neg * l

    return (l_new, u_new)
end

function solve(solver::DeepPoly, problem::Problem)
    isbounded(problem.input) || throw(UnboundedInputError("DeepPoly can only handle bounded input sets."))

    nnet, output = problem.network, problem.output
    domain = init_deep_poly_symbolic_interval(problem.input)

    reach = forward_network(solver, nnet, domain, collect=true)
    return reach
end

function forward_linear(solver::DeepPoly, L::LayerNegPos, input::AbstractHyperrectangle)
    # as Hyperrectangle is no subtype of SymbolicInterval, we need separate
    # method for the input layer
    domain = init_deep_poly_symbolic_interval(input)
    return forward_linear(solver, L, domain)
end

function forward_linear(solver::DeepPoly, L::LayerNegPos, input::SymbolicInterval)
    #output_Low, output_Up = interval_map(L.weights, input.Low, input.Up)
    output_Low, output_Up = interval_map(L.W_neg, L.W_pos, input.Low, input.Up)
    output_Up[:, end] += L.bias
    output_Low[:, end] += L.bias
    return SymbolicInterval(output_Low, output_Up, domain(input))
end

function forward_act(solver::DeepPoly, L::LayerNegPos{ReLU}, input::SymbolicInterval)
    n_node = n_nodes(L)
    output_Low, output_Up = copy(input.Low), copy(input.Up)
    los, his = bounds(input)

    @inbounds for j in 1:n_node
        # up_low, up_up = bounds(upper(input), j)
        # low_low, low_up = bounds(lower(input), j)

        #up_up = upper_bound(upper(input), j)
        #low_low = lower_bound(lower(input), j)
        low_low, up_up = los[j], his[j]

        slope = relaxed_relu_gradient(low_low, up_up)

        output_Up[j, :] .*= slope
        # only -low_low, if ReLU is not fixed
        output_Up[j, end] += slope * max(-low_low, 0)

        if up_up <= 0 || low_low > 0
            output_Low[j, :] .*= slope
        elseif up_up <= abs(low_low)
            output_Low[j, :] .*= 0
        end
        # if up_up > abs(low_low) do nothing, since output_Low is already a
        # copy of input.sym.Low
    end

    return SymbolicInterval(output_Low, output_Up, domain(input))
end

function forward_act(solver::DeepPoly, L::LayerNegPos{Id}, input::SymbolicInterval)
    return input
end

function bounds(sym::SymbolicInterval{<:Hyperrectangle})
    input_set = domain(sym)::Hyperrectangle
    input_los = input_set.center - input_set.radius
    input_his = input_set.center + input_set.radius

    # TODO: maybe use views here instead of copying the subarray into W, b
    W = sym.Low[:, 1:dim(input_set)]
    b = sym.Low[:, end]

    W_pos = max.(W, 0)
    W_neg = min.(W, 0)

    lb = W_pos * input_los .+ W_neg * input_his .+ b

    W = sym.Up[:, 1:dim(input_set)]
    b = sym.Up[:, end]

    # .= because W_pos, W_neg in appropriate size were already allocated before
    W_pos .= max.(W, 0)
    W_neg .= min.(W, 0)

    ub = W_pos * input_his .+ W_neg * input_los .+ b

    return lb, ub
end

# ρ(d, M) is max(bounds(d^T M))
function LazySets.ρ(d::AbstractArray{T,1} where T, sym::NeuralVerification.SymbolicInterval)
    # applying direction to symbolic interval is the same as passing the interval through a linear layer
    # (without bias), the result should be a 1-dimensional
    d_Low, d_Up = NeuralVerification.interval_map(reshape(d, 1, length(d)), sym.Low, sym.Up)
    sym_prime = NeuralVerification.SymbolicInterval(d_Low, d_Up, NeuralVerification.domain(sym))
    up_lo, up_up = bounds(upper(sym_prime), 1)
    return up_up
end


text:
using AbInitioSoftwareBase
using Test

@testset "AbInitioSoftwareBase.jl" begin
    # Write your tests here.
    include("fileops.jl")
    include("format.jl")
    include("Commands.jl")
end


text:
module CUDNN

issimilar(x,y) = (typeof(x) === typeof(y) && (x === nothing || size(x) === size(y)))

include("autograd.jl")
include("activation.jl")
include("softmax.jl")
include("pooling.jl")
include("dropout.jl")
include("rnn.jl")
include("convolution.jl")
include("multiheadattn.jl")
include("normalization.jl")

end

## grep cudnnHandle_t
# ?add8: cudnnBackendExecute
# ?add8: cudnnBuildRNNDynamic
# add8: cudnnCTCLoss_v8
# add8: cudnnGetFoldedConvBackwardDataDescriptors
# add8: cudnnGetNormalizationBackwardWorkspaceSize
# add8: cudnnGetNormalizationForwardTrainingWorkspaceSize
# add8: cudnnGetNormalizationTrainingReserveSpaceSize
# +add8: cudnnGetRNNTempSpaceSizes
# add8: cudnnGetRNNWeightParams
# +add8: cudnnGetRNNWeightSpaceSize
# add8: cudnnNormalizationBackward
# add8: cudnnNormalizationForwardInference
# add8: cudnnNormalizationForwardTraining
# +add8: cudnnRNNBackwardData_v8
# +add8: cudnnRNNBackwardWeights_v8
# +add8: cudnnRNNForward
# +cudnnActivationBackward
# +cudnnActivationForward
# +cudnnAddTensor
# cudnnBatchNormalizationBackward
# cudnnBatchNormalizationBackwardEx
# cudnnBatchNormalizationForwardInference
# cudnnBatchNormalizationForwardTraining
# cudnnBatchNormalizationForwardTrainingEx
# +cudnnConvolutionBackwardBias
# +cudnnConvolutionBackwardData
# +cudnnConvolutionBackwardFilter
# +cudnnConvolutionBiasActivationForward
# +cudnnConvolutionForward
# -cudnnCreate
# cudnnCTCLoss
# -cudnnDestroy
# cudnnDivisiveNormalizationBackward
# cudnnDivisiveNormalizationForward
# +cudnnDropoutBackward
# +cudnnDropoutForward
# +cudnnDropoutGetStatesSize
# +cudnnFindConvolutionBackwardDataAlgorithm
# +cudnnFindConvolutionBackwardDataAlgorithmEx
# +cudnnFindConvolutionBackwardFilterAlgorithm
# +cudnnFindConvolutionBackwardFilterAlgorithmEx
# +cudnnFindConvolutionForwardAlgorithm
# +cudnnFindConvolutionForwardAlgorithmEx
# cudnnFusedOpsExecute
# cudnnGetBatchNormalizationBackwardExWorkspaceSize
# cudnnGetBatchNormalizationForwardTrainingExWorkspaceSize
# cudnnGetBatchNormalizationTrainingExReserveSpaceSize
# +cudnnGetConvolutionBackwardDataAlgorithm_v7
# +cudnnGetConvolutionBackwardDataAlgorithmMaxCount
# +cudnnGetConvolutionBackwardDataWorkspaceSize
# +cudnnGetConvolutionBackwardFilterAlgorithm_v7
# +cudnnGetConvolutionBackwardFilterAlgorithmMaxCount
# +cudnnGetConvolutionBackwardFilterWorkspaceSize
# +cudnnGetConvolutionForwardAlgorithm_v7
# +cudnnGetConvolutionForwardAlgorithmMaxCount
# +cudnnGetConvolutionForwardWorkspaceSize
# cudnnGetCTCLossWorkspaceSize
# cudnnGetCTCLossWorkspaceSize_v8
# +cudnnGetDropoutDescriptor
# cudnnGetMultiHeadAttnBuffers
# cudnnGetMultiHeadAttnWeights
# +cudnnGetReductionIndicesSize
# +cudnnGetReductionWorkspaceSize
# -cudnnGetStream
# -cudnnIm2Col
# cudnnLRNCrossChannelBackward
# cudnnLRNCrossChannelForward
# cudnnMakeFusedOpsPlan
# cudnnMultiHeadAttnBackwardData
# cudnnMultiHeadAttnBackwardWeights
# cudnnMultiHeadAttnForward
# +cudnnOpTensor
# +cudnnPoolingBackward
# +cudnnPoolingForward
# -cudnnQueryRuntimeError
# +cudnnReduceTensor
# ?cudnnReorderFilterAndBias
# +cudnnRestoreDropoutDescriptor
# +cudnnScaleTensor
# +cudnnSetDropoutDescriptor
# -cudnnSetStream
# +cudnnSetTensor
# +cudnnSoftmaxBackward
# +cudnnSoftmaxForward
# -cudnnSpatialTfGridGeneratorBackward
# -cudnnSpatialTfGridGeneratorForward
# -cudnnSpatialTfSamplerBackward
# -cudnnSpatialTfSamplerForward
# -cudnnTransformFilter
# -cudnnTransformTensor
# -cudnnTransformTensorEx
# dep8: cudnnFindRNNBackwardDataAlgorithmEx
# dep8: cudnnFindRNNBackwardWeightsAlgorithmEx
# dep8: cudnnFindRNNForwardInferenceAlgorithmEx
# dep8: cudnnFindRNNForwardTrainingAlgorithmEx
# dep8: cudnnGetAlgorithmSpaceSize
# dep8: cudnnGetRNNBackwardDataAlgorithmMaxCount
# dep8: cudnnGetRNNBackwardWeightsAlgorithmMaxCount
# dep8: cudnnGetRNNDescriptor_v6
# dep8: cudnnGetRNNForwardInferenceAlgorithmMaxCount
# dep8: cudnnGetRNNForwardTrainingAlgorithmMaxCount
# dep8: cudnnGetRNNLinLayerBiasParams
# dep8: cudnnGetRNNLinLayerMatrixParams
# dep8: cudnnGetRNNParamsSize
# dep8: cudnnGetRNNProjectionLayers
# dep8: cudnnGetRNNTrainingReserveSize
# dep8: cudnnGetRNNWorkspaceSize
# dep8: cudnnRestoreAlgorithm
# dep8: cudnnRNNBackwardData
# dep8: cudnnRNNBackwardDataEx
# dep8: cudnnRNNBackwardWeights
# dep8: cudnnRNNBackwardWeightsEx
# dep8: cudnnRNNForwardInference
# dep8: cudnnRNNForwardInferenceEx
# dep8: cudnnRNNForwardTraining
# dep8: cudnnRNNForwardTrainingEx
# dep8: cudnnRNNGetClip
# dep8: cudnnRNNSetClip
# dep8: cudnnSaveAlgorithm
# dep8: cudnnSetRNNAlgorithmDescriptor
# dep8: cudnnSetRNNDescriptor_v6
# dep8: cudnnSetRNNProjectionLayers


text:
module MPSensitivity

using UnROOT, Plots, StatsBase

include("MPHist.jl")
export get_bin_content
export MP_heatmap

include("MPEfficiency.jl")
export get_efficiency

include("MPStoB.jl")
export get_SB_ratio

include("MPROI.jl")
export get_best_ROI_index
export get_best_ROI

end # module MPSensitivity


text:
function bubblesort!(x::AbstractArray)
    n = length(x)
    for i ∈ 1:n
        for j ∈ (i+1):n
            @inbounds if x[j] ≤ x[i]
                @inbounds tmp  = x[i]
                @inbounds x[i] = x[j]
                @inbounds x[j] = tmp
            end
        end
    end
    return nothing
end


text:
#=
A specialization of the
the resynthesizer algorithm,
inpainting from context that is a nearby frisket.
The results are better (than resynthesizing from the entire corpus)
since pragmatically, nearby corpus is usually a "background".

Takes one image and a mask.

Returns new image with the masked region synthesized from its *nearby* surroundings.

The one image is first divided into a target and corpus image.
The target image is the input.
The algorithm mutates the target image in the masked region, called the "synth" (a noun) region.

The corpus is inverse of mask applied to input.
The corpus is immutable.
=#

include("../resynthesizer.jl")

include("../image/boundingBox.jl")
include("../image/mask.jl")

function healSelection(
    image::AbstractArray{T,N},      # target
    mask::AbstractArray{Bool,N},    # synth region
    frisketDepth::Int64
    ) where {N,T}

    # Assert image and mask are similar

    # Target for resynthesize() is as passed
    targetImage = MaskedImage(image, mask)         # image is mutable

    #=
    The corpus image is a copy of the in image.
    The corpus is not mutated.

    The corpus mask is a frisket.
    Slightly larger than,
    and surrounding the masked (selected) area of the target,
    with a hole in it where is the selected area of the target.
    =#

    copyImage = copy(image)

    corpusMask = frisketMaskAroundMask(image, mask, frisketDepth)

    # An optimization would be to subarray the image and the mask

    corpusImage = copyImage

    corpusImage = MaskedImage(corpusImage, corpusMask)

    resynthesize(targetImage, corpusImage)
end

max_stars_repo_name: bootchk/resynthesizerJulia

max_issues_repo_name: bootchk/resynthesizerJulia

max_forks_repo_name: bootchk/resynthesizerJulia


text:
using Echidna
using Plots

gr()

f1(θ) = [cos(θ), sin(θ)]

pts = [f1(π/2 * x) for x in rand(10)]
x = [p[1] for p in pts]
y = [p[2] for p in pts]
scatter(x, y)

zmin = [1.0, 1.0]
extremepts = find_extreme_points(pts, zmin)

xx = [x[i] for i in extremepts]
yy = [y[i] for i in extremepts]
scatter!(xx, yy)

f2(x) = [x, exp(-x)]
pts = [f2(5*x) for x in rand(10)]
x = [p[1] for p in pts]
y = [p[2] for p in pts]
scatter(x, y)

extremepts = find_extreme_points(pts)

xx = [x[i] for i in extremepts]
yy = [y[i] for i in extremepts]
scatter!(xx, yy)

zpts = [pts[i] for i in extremepts]
find_axis_intercepts(zpts)

#### TESTING POINT ASSOCIATION
canpoints = [[0.5, 1.0], [2.0, 0.0]]
refpoints = [[1.0, 0.0], [0.5, 0.5], [0.0, 1.0]]

index, dists = associate_points(canpoints, refpoints)

max_stars_repo_name: phawthorne/echidna

max_issues_repo_name: phawthorne/echidna

max_forks_repo_name: phawthorne/echidna


text:
using PooksoftAlphaVantageDataStore

# ok, so lets build by user object -
my_current_dir = pwd()   # where am I?
path_to_config_file = my_current_dir*"/test/configuration/Configuration.json"

# build the api user model -
user_model = build_api_user_model(path_to_config_file) |> checkresult

# setup call -
requestDictionary = Dict{String,Any}()
requestDictionary["apikey"] = user_model.alphavantage_api_key
requestDictionary["symbol"] = "ALLY"
requestDictionary["function"] = "TIME_SERIES_INTRADAY"
requestDictionary["interval"] = "30min"
requestDictionary["adjusted"] = "true"
requestDictionary["outputsize"] = "compact"
requestDictionary["datatype"] = "json"
api_call_result = execute_sts_intraday_api_call(requestDictionary)
df = checkresult(api_call_result)


text:
using IntervalArithmetic
using Test

setformat(:full)

# NB: Due to "world age" problems, do *not* wrap everything in another @testset
# The `setrounding` commands must be at top level

setrounding(Interval, :slow)
x = Interval(0.5)
@testset "Tight, slow rounding by changing rounding mode" begin
    @test rounding(Interval) == :slow
    @test sin(x) == Interval(0.47942553860420295, 0.479425538604203)
end

setrounding(Interval, :accurate)
@testset "Accurate, fast rounding using prevfloat and nextfloat" begin
    @test rounding(Interval) == :accurate
    @test sin(x) == Interval(0.47942553860420295, 0.47942553860420306)
end

setrounding(Interval, :none)
@testset "No rounding" begin
    @test rounding(Interval) == :none
    @test sin(x) == Interval(0.479425538604203, 0.479425538604203)
    w = 0 ± big(1)
    @test isa(w+w, Interval)
end

setrounding(Interval, :tight)
@testset "Tight, fast rounding using FastRounding.jl" begin
    @test rounding(Interval) == :tight
    @test sin(x) == Interval(0.47942553860420295, 0.479425538604203)
end

setformat(:standard)

# end
# end