text:
# ## [`DatasetRecipe`] implementations

"""
    ImageFolders(; labelfn = parentname, split = false)

Recipe for loading a single-label image classification dataset
stored in a hierarchical folder format. If `split == true`, split
the data container on the name of the grandparent folder. The label
defaults to the name of the parent folder but a custom function can
be passed as `labelfn`.

```julia
julia> recipeblocks(ImageFolders)
Tuple{Image{2}, Label}
```
"""
Base.@kwdef struct ImageFolders <: Datasets.DatasetRecipe
    labelfn = parentname
    split::Bool = false
    filefilterfn = _ -> true
end

function Datasets.loadrecipe(recipe::ImageFolders, path)
    isdir(path) || error("$path is not a directory")
    data = loadfolderdata(
        path,
        filterfn=f -> isimagefile(f) && recipe.filefilterfn(f),
        loadfn=(loadfile, recipe.labelfn),
        splitfn=recipe.split ? grandparentname : nothing)

    (recipe.split ? length(data) > 0 : nobs(data) > 0) || error("No image files found in $path")

    labels = recipe.split ? first(values(data))[2] : data[2]
    blocks = Image{2}(), Label(unique(eachobs(labels)))
    length(blocks[2].classes) > 1 || error("Expected multiple different labels, got: $(blocks[2].classes))")
    return data, blocks
end

Datasets.recipeblocks(::Type{ImageFolders}) = Tuple{Image{2},Label}

"""
    ImageSegmentationFolders(; imagefolder="images", maskfolder="labels", labelfile="codes.txt")

Dataset recipe for loading 2D image segmentation datasets from a common format
where images and masks are stored as images in two different subfolders
"<root>/<imagefolder>" and "<root>/<maskfolder>"
The class labels should be in a newline-delimited file "<root>/<labelfile>".
"""
Base.@kwdef struct ImageSegmentationFolders <: Datasets.DatasetRecipe
    imagefolder::String = "images"
    maskfolder::String = "labels"
    labelfile::String = "codes.txt"
end

function Datasets.loadrecipe(recipe::ImageSegmentationFolders, path)
    isdir(path) || error("$path is not a directory")
    imagepath = joinpath(path, recipe.imagefolder)
    maskpath = joinpath(path, recipe.maskfolder)
    classespath = joinpath(path, recipe.labelfile)

    isdir(imagepath) || error("Image folder $imagepath is not a directory")
    isdir(maskpath) || error("Mask folder $maskpath is not a directory")

    isfile(classespath) || error("Classes file $classespath does not exist")
    classes = readlines(open(joinpath(path, recipe.labelfile)))
    length(classes) > 1 || error("Expected multiple different labels, got: $(blocks[2].classes))")

    images = loadfolderdata(imagepath, filterfn=isimagefile, loadfn=loadfile)
    masks = loadfolderdata(maskpath, filterfn=isimagefile, loadfn=f -> loadmask(f, classes))
    nobs(images) == nobs(masks) || error("Expected the same number of images and masks, but found $(nobs(images)) images and $(nobs(masks)) masks")
    nobs(images) > 0 || error("No images or masks found in folders $imagepath and $maskpath")

    blocks = Image{2}(), Mask{2}(classes)
    return (images, masks), blocks
end

Datasets.recipeblocks(::Type{ImageSegmentationFolders}) = Tuple{Image{2},Mask{2}}

# ImageTableMultiLabel

Base.@kwdef struct ImageTableMultiLabel <: Datasets.DatasetRecipe
    csvfile::String = "train.csv"
    imagefolder::String = "train"
    filecol::Symbol = :fname
    fileext::String = ""
    labelcol::Symbol = :labels
    split::Bool = false
    splitcol::Symbol = :is_valid
    labelsep::String = " "
end

function Datasets.loadrecipe(recipe::ImageTableMultiLabel, path)
    csvpath = joinpath(path, recipe.csvfile)
    isfile(csvpath) || error("File $csvpath does not exist")
    df = loadfile(csvpath)
    images = mapobs(
        f -> loadfile(joinpath(path, recipe.imagefolder, f * recipe.fileext)),
        df[:, recipe.filecol])
    labels = map(
        str -> split(str, recipe.labelsep),
        df[:,recipe.labelcol])
    data = (images, labels)
    blocks = Image{2}(), LabelMulti(unique(Iterators.flatten(labels)))
    if recipe.split
        idxs = 1:nobs(data)
        splits = df[:, recipe.splitcol]
        data = Dict(
            "train" => datasubset(data, idxs[splits]),
            "valid" => datasubset(data, idxs[(!).(splits)])
        )
    end
    return data, blocks
end

Datasets.recipeblocks(::Type{ImageTableMultiLabel}) = Tuple{Image{2},LabelMulti}

# ## Registering recipes for fastai datasets

const RECIPES = Dict{String,Vector{Datasets.DatasetRecipe}}(
    [name => [ImageFolders()] for name in (
        "imagenette", "imagenette-160", "imagenette-320",
        "imagenette2", "imagenette2-160", "imagenette2-320",
        "imagewoof", "imagewoof-160", "imagewoof-320",
        "imagewoof2", "imagewoof2-160", "imagewoof2-320",
        "cifar10", "cifar100", "caltech_101", "mnist_png",
        "mnist_sample", "CUB_200_2011"
    )]...,
    [name => [ImageFolders(filefilterfn=f -> !(occursin("unsup", f)))]
        for name in ("imagewang-160", "imagewang-320", "imagewang")]...,
    "camvid" => [ImageSegmentationFolders()],
    "camvid_tiny" => [ImageSegmentationFolders()],
    "pascal_2007" => [ImageTableMultiLabel()],
    "mnist_tiny" => [ImageFolders(filefilterfn=f -> !occursin("test", f))],
    "mnist_var_size_tiny" => [ImageFolders(filefilterfn=f -> !occursin("test", f))],
)

function _registerrecipes()
    for (name, recipes) in RECIPES, recipe in recipes
        Datasets.registerrecipe!(Datasets.FASTAI_DATA_REGISTRY, name, recipe)
    end
end

# ## Tests

@testset "ImageFolders [recipe]" begin
    path = joinpath(datasetpath("mnist_var_size_tiny"), "train")

    @testset "Basic configuration" begin
        recipe = Vision.ImageFolders()
        data, blocks = loadrecipe(recipe, path)
        Datasets.testrecipe(recipe, data, blocks)
        @test blocks[1] isa Image
        @test blocks[2].classes == ["3", "7"]
    end

    @testset "Split configuration" begin
        recipe = Vision.ImageFolders(split=true)
        data, blocks = loadrecipe(recipe, path)
        Datasets.testrecipe(recipe, data["train"], blocks)
    end

    @testset "Error cases" begin
        @testset "Empty directory" begin
            recipe = Vision.ImageFolders(split=true)
            @test_throws ErrorException loadrecipe(recipe, mktempdir())
        end

        @testset "Only one label" begin
            recipe = Vision.ImageFolders(labelfn=x -> "1")
            @test_throws ErrorException loadrecipe(recipe, path)
        end
    end

end

@testset "ImageSegmentationFolders [recipe]" begin
    path = datasetpath("camvid_tiny")

    @testset "Basic configuration" begin
        recipe = ImageSegmentationFolders()
        data, blocks = loadrecipe(recipe, path)
        Datasets.testrecipe(recipe, data, blocks)
        @test blocks[1] isa Image
        @test blocks[2] isa Mask
    end

    @testset "Error cases" begin
        @testset "Empty directory" begin
            recipe = ImageSegmentationFolders()
            @test_throws ErrorException loadrecipe(recipe, mktempdir())
        end

        @testset "Only one label" begin
            recipe = ImageSegmentationFolders(labelfile="idontexist")
            @test_throws ErrorException loadrecipe(recipe, path)
        end
    end
end


text:
#The m-by-n Hilbert matrix has matrix elements
# H_{ij} = 1/(i+j-1)
export Hilbert, InverseHilbert

struct Hilbert{T} <: AbstractMatrix{T}
    m :: Int
    n :: Int
end
Hilbert(T::Type, m::Integer, n::Integer) = Hilbert{T}(m, n)
Hilbert(T::Type{<:Integer}, m::Integer, n::Integer) = Hilbert(Rational{T}, m, n)
Hilbert(m::Integer, n::Integer) = Hilbert(Int, m, n)
Hilbert(n::Integer) = Hilbert(n, n)

# Define its size
size(H::Hilbert, dim::Integer) = dim==1 ? H.m : dim==2 ? H.n : 1
size(H::Hilbert)= size(H,1), size(H,2)

# Index into a Hilbert matrix
function getindex(H::Hilbert{T}, i::Integer, j::Integer) where {T}
    return one(T)/(i+j-1)
end

# Dense version
Matrix(H::Hilbert) = [H[i,j] for i=1:size(H,1), j=1:size(H,2)]

# Some properties
ishermitian(H::Hilbert) = H.m==H.n
isposdef(H::Hilbert) = H.m==H.n
det(H::Hilbert) = inv(det(inv(H)))

# Inverse of a Hilbert matrix
struct InverseHilbert{T} <: AbstractMatrix{T}
    n :: Int
end
InverseHilbert(T::Type, n::Integer) = InverseHilbert{T}(n)
InverseHilbert(n::Integer) = InverseHilbert(Int, n)

# Define its size
size(A::InverseHilbert, dim::Integer) = dim==1 || dim==2 ? A.n : 1
size(A::InverseHilbert) = size(A,1), size(A,2)

# Index into a inverse Hilbert matrix
function getindex(A::InverseHilbert{T}, i::Integer, j::Integer) where {T}
    return T((-1)^(i+j)*(i+j-1)*binomial(A.n+i-1,A.n-j)*
              binomial(A.n+j-1,A.n-i)*binomial(i+j-2,i-1)^2)
end
# Use bigger numbers in case of bigger types
function getindex(A::InverseHilbert{T}, i::Integer, j::Integer) where {T<:Union{BigInt,BigFloat}}
    N = big(A.n)
    return T((-1)^(i+j)*(i+j-1)*binomial(N+i-1,N-j)*
              binomial(N+j-1,N-i)*binomial(big(i+j-2),i-1)^2)
end

# Explicit formula for the determinant
det(A::InverseHilbert{T}) where {T} = prod(T,(2k+1)*binomial(2k,k)^2 for k=1:A.n-1)

# Dense version
Matrix(A::InverseHilbert) = [A[i,j] for i=1:size(A,1), j=1:size(A,2)]

# Define Base.inv
function inv(H::Hilbert{T}) where {T}
    H.m == H.n || throw(ArgumentError("Works only for square Hilbert matrices."))
    return InverseHilbert(T,H.n)
end
function inv(A::InverseHilbert{T}) where {T}
    HT = promote_type(T,typeof(Rational(one(T))))
    return Hilbert(HT,A.n)
end

ext: jl


text:
"""
    ltbtmatrix(a,c)

Build a lower triangular banded Toeplitz (LTBT) matrix corresponding to a univariate polynomial `a` and with the number of columns given by the integer `c`. The matrix is stored as dense.

# Examples

```julia
julia> a = Polynomial([1.0,2.0,3.0,4.0],:s)
julia> ltbtmatrix(a,3)
6×3 Array{Real,2}:
 1.0  0    0
 2.0  1.0  0
 3.0  2.0  1.0
 4.0  3.0  2.0
 0    4.0  3.0
 0    0    4.0
```
"""
function ltbtmatrix(a::Polynomial{T},c::Signed) where {T}
    da = degree(a)          # degree of the polynomial a
    r = da+c                # number of rows of the resulting matrix
    M = zeros(T,r,c)
    for i = 1:r
        for j = 1:c
            if i >= j && i <= (j+da)
                M[i,j] = a[i-j]
            end
        end
    end
    return M
end

"""
    ltbtmatrix(a,c)

Build a lower triangular banded Toeplitz (LTBT) matrix corresponding to a univariate Laurent polynomial `a` and with the number of columns given by the integer `c`. The matrix is stored as dense.

# Examples

```julia
julia> a = LaurentPolynomial([1,2,3,4,5],-2:2,:z)
LaurentPolynomial(z⁻² + 2*z⁻¹ + 3 + 4*z + 5*z²)

julia> A = ltbtmatrix(a,3)
7×3 Array{Float64,2}:
 1.0  0.0  0.0
 2.0  1.0  0.0
 3.0  2.0  1.0
 4.0  3.0  2.0
 5.0  4.0  3.0
 0.0  5.0  4.0
 0.0  0.0  5.0
```
"""
function ltbtmatrix(a::LaurentPolynomial{T},c::Signed) where {T}
    m = a.m[]
    n = a.n[]
    d = n-m+1               # length of the coefficient vector
    r = d+c-1               # number of rows of the resulting matrix
    M = zeros(T,r,c)
    for i = 1:r
        for j = 1:c
            if i >= j && i <= (j+d)
                M[i,j] = a[i-j+m]   # but assumes that m ≦ 0
            end
        end
    end
    return M
end

"""
    sylvestermatrix(a,b)

Build a Sylvester resultant matrix for a pair of univariate polynomials `a` and `b`.

The matrix is obtained by stacking horizontally two lower triangular banded Toeplitz matrices corresponding to the two polynomials. The matrix is square. Its nonsingularity is used to detect coprimeness of the two polynomials.

The matrix is stored as dense.

# Examples

```julia
julia> a = Polynomial([1,2,3],:s)
Polynomial(1 + 2*s + 3*s^2)

julia> b = Polynomial([4,5],:s)
Polynomial(4 + 5*s)

julia> S = sylvestermatrix(a,b)
3×3 Array{Int64,2}:
 1  4  0
 2  5  4
 3  0  5
```
"""
function sylvestermatrix(a::Polynomial{T},b::Polynomial{T}) where {T}
    da = degree(a)
    db = degree(b)
    dx = degree(b)-1
    dy = da+dx-db
    m = 1+da+dx                    # number of rows: da+db
    n = 1+dx+1+dy                  # number of columns: da+db
    S = zeros(T,m,n)
    for i = 1:m
        for j = 1:n
            if j <= dx+1           # if in the "left part" of the Sylvester matrix
                S[i,j] = i >= j && i <= (j+da) ? a[i-j] : 0.0
            else
                k = j-(dx+1)       # index of the column in the "right part" of the Sylvester matrix
                S[i,j] = i >= k && i <= (k+db) ? b[i-k] : 0.0
            end
        end
    end
    return S
end

"""
    sylvestermatrix(a,b,w)

Build a Sylvester matrix for a pair of univariate polynomials `a` and `b`. The matrix is obtained by stacking horizontally two lower triangular banded Toeplitz matrices corresponding to the two polynomials `a` and `b`. The third input argument `w` is a tuple whose two components specify the widths of the two Toeplitz blocks. The matrix is stored as dense.

# Examples

```julia
julia> a = Polynomial([1,2,3],:s)
Polynomial(1 + 2*s + 3*s^2)

julia> b = Polynomial([4,5],:s)
Polynomial(4 + 5*s)

julia> S = sylvestermatrix(a,b,(1,3))
3×3 Array{Int64,2}:
 1  4  0  0
 2  5  4  0
 3  0  5  4
 0  0  0  5
```
"""
function sylvestermatrix(a::Polynomial{T},b::Polynomial{T},w::Tuple{Int,Int}) where {T}
    da = degree(a)
    db = degree(b)
    m = max(da+w[1],db+w[2])       # number of rows of the resulting matrix
    n = w[1]+w[2]                  # number of columns
    S = zeros(T,m,n)
    for i = 1:m
        for j = 1:n
            if j <= w[1]           # if in the "left part" of the Sylvester matrix
                S[i,j] = i >= j && i <= (j+da) ? a[i-j] : 0.0
            else
                k = j-w[1]         # index of the column in the "right part" of the Sylvester matrix
                S[i,j] = i >= k && i <= (k+db) ? b[i-k] : 0.0
            end
        end
    end
    return S
end


text:
using Crayons

function Base.show(io::IO, gw::AbstractGridWorld)
    p, d = get_agent_pos(gw), get_agent_dir(gw)
    w = convert(GridWorldBase, gw)

    println(io, "World:")
    for i in 1:size(w, 2)
        for j in 1:size(w, 3)
            if CartesianIndex(i, j) ∈ get_agent_view_inds(gw)
                bg = :dark_gray
            else
                bg = :black
            end
            if i == p[1] && j == p[2]
                agent = get_agent(gw)
                print(io, Crayon(background=bg, foreground=get_color(agent),reset=true), convert(Char, agent))
            else
                o = get_object(gw)[findfirst(w.world[:, i, j])]
                print(io, Crayon(background=bg, foreground=get_color(o),reset=true),  convert(Char, o))
            end
        end
        println(io, Crayon(reset=true))
    end
    println(io)

    println(io, "Agent's view:")
    v = get_agent_view(gw)
    for i in 1:size(v, 2)
        for j in 1:size(v, 3)
            if i == 1 && j == size(v, 3) ÷ 2 + 1
                print(io, Agent(dir=DOWN))
            else
                x = findfirst(v[:, i, j])
                if isnothing(x)
                    print(io, '_')
                else
                    print(io, get_object(gw)[x])
                end
            end
        end
        println(io)
    end
    println(io)
end


text:
using Revise
using Test, QuHamiltonian, Lattices

# set up something useful
using LinearAlgebra, SparseArrays
⊗ = kron
II = Matrix{Float64}(I, 2, 2)
##########################

@testset "vertex sum" begin

    h = @sum Z[@vertex(i)]
    ltc = Chain(4)
    H = h(ltc)

    A = Z ⊗ II ⊗ II ⊗ II + II ⊗ Z ⊗ II ⊗ II + II ⊗ II ⊗ Z ⊗ II + II ⊗ II ⊗ II ⊗ Z

    lhs = Bit[0, 0, 0, 0]
    rhs = Bit[0, 0, 0, 0]

    for i = 1:2^4
        for j = 1:2^4
            @test H[lhs, rhs] == A[convert(Int, lhs) + 1, convert(Int, rhs) + 1]
            lhs << 1
        end
        rhs << 1
    end

    @test Matrix(H) == A
    @test sparse(H) == sparse(A)
end

@testset "edge sum" begin

    @nearest i, j
    h = @sum Z[i]Z[j]
    ltc = Chain(4)
    H = h(ltc)

    lhs = Bit[0, 0, 0, 0]
    rhs = Bit[0, 0, 0, 0]

    A = Z ⊗ Z ⊗ II ⊗ II + II ⊗ Z ⊗ Z ⊗ II + II ⊗ II ⊗ Z ⊗ Z + Z ⊗ II ⊗ II ⊗ Z
    Matrix(H)
    for i = 1:1<<4
        for j = 1:1<<4
            @test H[lhs, rhs] == A[convert(Int, lhs) + 1, convert(Int, rhs) + 1]
            lhs << 1
        end
        rhs << 1
    end

    @test Matrix(H) == A
    @test sparse(H) == sparse(A)
end

@testset "add/sub" begin

    @nearest i, j
    @vertex l
    h = @sum(X[i]X[j]) - 0.1 * @sum(Z[l])
    ltc = Chain(4)
    H = h(ltc)

    lhs = Bit[0, 0, 0, 0]
    rhs = Bit[0, 0, 0, 0]

    H[lhs, rhs]
end

# @testset "multiple vertex sum" begin
#
#     @vertex i, j
#     h = @sum Z[i]X[j] + @sum (Z * X)[i]
#     ltc = Chain(4)
#     H = h(ltc)
#
#     A =
#     # (2, 1)
#     X ⊗ Z ⊗ II ⊗ II +
#     # (3, 1)
#     X ⊗ II ⊗ Z ⊗ II +
#     # (4, 1)
#     X ⊗ II ⊗ II ⊗ Z +
#
#     # (1, 2)
#     Z ⊗ X ⊗ II ⊗ II +
#     # (3, 2)
#     II ⊗ X ⊗ Z ⊗ II +
#     # (4, 2)
#     II ⊗ X ⊗ II ⊗ Z +
#
#     # (1, 3)
#     Z ⊗ II ⊗ X ⊗ II +
#     # (2, 3)
#     II ⊗ Z ⊗ X ⊗ II +
#     # (4, 3)
#     II ⊗ II ⊗ X ⊗ Z +
#
#     # (1, 4)
#     Z ⊗ II ⊗ II ⊗ X +
#     # (2, 4)
#     II ⊗ Z ⊗ II ⊗ X +
#     # (3, 4)
#     II ⊗ II ⊗ Z ⊗ X +
#
#     # (1, 1)
#     (Z * X) ⊗ II ⊗ II ⊗ II +
#     # (2, 2)
#     II ⊗ (Z * X) ⊗ II ⊗ II +
#     # (3, 3)
#     II ⊗ II ⊗ (Z * X) ⊗ II +
#     # (4, 4)
#     II ⊗ II ⊗ II ⊗ (Z * X)
#
#     Matrix(H) == A
#
# end

max_stars_repo_name: Roger-luo/QuHamiltonian

max_issues_repo_name: Roger-luo/QuHamiltonian

max_forks_repo_name: Roger-luo/QuHamiltonian


text:
using ILog2
using Test

# FIXME: Tests for other floating point types

@testset "ILog2" begin
    bitstypes = (Int8, Int16, Int32, Int64,
                 UInt8, UInt16, UInt32, UInt64, Int128, UInt128)
    largest_Float64_corresponding_to_unique_integer = 2^53 - 1
    for T in  bitstypes
        N = min(typemax(T) - 1, largest_Float64_corresponding_to_unique_integer)
        num_trials_per_type = min(N, 10^4)
        nums = rand(T(1):T(N), num_trials_per_type)
        for n in nums
            @test ilog2(T(n)) == ilog2(float(n))
        end
    end
end

@testset "BigInt" begin
    for expt in rand(10:10^4, 10^3)
        n = big(2)^expt
        np = n + rand(1:100)
        nm = n - rand(1:100)
        @test ilog2(n) == expt
        @test ilog2(np) == expt
        @test ilog2(nm) == expt - 1
        @test ilog2(n) == ilog2(float(n))
    end
end

@testset "exceptions" begin
    @test_throws ArgumentError ILog2.msbindex(BigInt)
    @test_throws InexactError  ilog2(float(typemax(Int)-2^8))
    @test_throws DomainError ilog2(0)
    @test_throws DomainError ilog2(-1)
end


text:
using Compat.LinearAlgebra: copytri!

function state_likelihood_model(Eₘ_list)
    # Generate the A matrix used to calculate likelihoods
    # The A matrix depends on the measurement operators
    sum(abs,data(sum(Eₘ_list))-I)<1E-14 ||
        throw(ArgumentError("Eₘ operators do not form a valid POVM!"))
    return copy(transpose(mapreduce(vec∘full,hcat,Eₘ_list)::Matrix{ComplexF64}))
end

function mle_state_tomo(M,A)
    (nm,d²) = size(A); d = isqrt(d²)
    length(M) == nm || throw(ArgumentError("number of measurements does not match A-matrix size!"))
    d^2 == d² || throw(ArgumentError("invalid number of columns in A-matrix!"))
    f = let ρ = Matrix{ComplexF64}(undef,d,d), X = Vector{ComplexF64}(undef,nm)
        t -> loglikehood_binomial(t,M,A,ρ,X,zero(ρ))
    end
    return optimize(f,[fill(1/√d,d); zeros(d²-d)])
end

function loglikehood_gaussian(t,M,A,ρ,X,T)
    # Gaussian statistics for the measurement count probability
    t_expectation_values!(X,t,A,ρ,T)
    return sum((real(p)-m)^2/2p for (p,m) in zip(X,M))
end

function loglikehood_binomial(t,M,A,ρ,X,T)
    # Binomial statistics for the measurement count probability, up to some irrelevant constant
    t_expectation_values!(X,t,A,ρ,T)
    return sum(-m*log(real(p)) for (p,m) in zip(X,M))
end

function t_expectation_values!(X,t::AbstractVector,A::Matrix,ρ::Matrix,T::Matrix)
    # calculate ρ from the t-vector
    if length(t) == 4
        @inbounds t₁,t₂,t₃,t₄ = t
        build_density_matrix!(ρ,t₁,t₂,t₃,t₄)
    else
        build_density_matrix!(ρ,t,T)
    end
    # calculate the probabilities by multiplying the vectorized density
    # operator with the A matrix
    mul!(X,A,vec(ρ))
end

function build_density_matrix(t::AbstractVector{S}) where S<:Real
    d = isqrt(length(t))
    d^2 == length(t) || throw(ArgumentError("invalid t-vector length!"))
    ρ = Matrix{complex(float(S))}(undef,d,d)
    return build_density_matrix!(ρ,t,zero(ρ))
end

function build_density_matrix!(ρ::Matrix,t₁::Real,t₂::Real,t₃::Real,t₄::Real)
    # reconstruct the density matrix from the t-parameters, special 2-d case
    t₁²,t₂²,t₃²,t₄² = t₁^2,t₂^2,t₃^2,t₄^2
    N = t₁²+t₂²+t₃²+t₄²
    ρ[1,1] = (t₁²+t₃²+t₄²)/N;  ρ[1,2] = t₂*(t₃-1im*t₄)/N
    ρ[2,1] = t₂*(t₃+1im*t₄)/N; ρ[2,2] = t₂²/N
    return ρ
end

function build_density_matrix!(ρ::Matrix,t::AbstractVector,T::Matrix)
    # reconstruct the density matrix from the t-parameters, general case
    # T must have its upper triangle zeroed out
    size(ρ) == size(T) || throw(DimensionMismatch("ρ & T do not have the same dimensions!"))
    d = isqrt(length(t))
    (d,d) == size(T) || throw(DimensionMismatch("ρ & T incompatible with the t-vector!"))
    for i = 1:d
        @inbounds T[i,i] = t[i]
    end
    r,c = 2,1
    for i = d+1:2:d^2
        @inbounds T[r,c] = complex(t[i],t[i+1])
        r += 1
        if r > d
            c = c + 1
            r = c + 1
        end
    end
    α = 1/real(dot(T,T))
    # Ac_mul_B!(ρ,T,T)
    # ρ .= α .* ρ
    copytri!(BLAS.herk!('U', 'C', α, T, zero(α), ρ), 'U', true)
    return ρ
end

M = [103, 633, 621, 115, 377, 354]

# only valid for a certain set of measurement, and for a single qubit
# here, the inputs are count values for |+⟩,|-⟩,|+i⟩,|-i⟩,|0⟩,|1⟩
# measurements, e.g. M[1] = Nm*⟨+|ρ|+⟩
# a physical density matrix can be represented by 4 parameters; assuming
# some measurement statistics, we minimize the loglikehood function
# t_expectation_values(t::Array) = (@inbounds t₁,t₂,t₃,t₄=t; t_expectation_values(t₁,t₂,t₃,t₄))
#
# function t_expectation_values(t₁,t₂,t₃,t₄)
#     # Calculate the expectation values for |+⟩,|-⟩,|+i⟩,|-i⟩,|0⟩,|1⟩ with the
#     # t-parameterized density matrix ρ(t) = T†*T/trace(T†*T), where
#     # T = [  t₁    0
#     #      t₃+it₄  t₂]
#     t₁²,t₂²,t₃²,t₄² = t₁^2,t₂^2,t₃^2,t₄^2
#     N = t₁²+t₂²+t₃²+t₄²
#     # ⟨+|ρ|+⟩, ⟨-|ρ|-⟩
#     Xp, Xm = 1/2 + t₂*t₃/N, 1/2 - t₂*t₃/N
#     # ⟨+i|ρ|+i⟩, ⟨-i|ρ|-i⟩
#     Yp, Ym = 1/2 + t₂*t₄/N, 1/2 - t₂*t₄/N
#     # ⟨0|ρ|0⟩, ⟨1|ρ|1⟩
#     Zp, Zm = (t₁²+t₃²+t₄²)/N, t₂²/N
#     return Xp,Xm,Yp,Ym,Zp,Zm
# end


text:
include("TestFeature.jl")
function testInstance(outInstance::Ref, outGlfw::Ref)
    # Create a Vulkan instance, tell it we need glfw and the validation as
    # extension features
    fs = Vector{features.IFeature}()
    fsValid = features.create(features.Validation)
    glfw = features.create(features.GlfwOutput)
    fsTest = create(TestFeature)
    fs = [fsValid, glfw, fsTest]
    instance = lava.create(lava.Instance, fs)

    @test fsTest.mOnInstanceCreatedCalled
    @test lava.handle(instance) != vk.VK_NULL_HANDLE
    @test glfw.mTempWindow != vk.VK_NULL_HANDLE
    @test glfw.mTempSurface != vk.VK_NULL_HANDLE

    outInstance[] = instance
    outGlfw[] = glfw
    return true
end


text:
#                        .-'''-.                               _..._
#                       '   _    \         _______          .-'_..._''.
#  /|                 /   /` '.   \        \  ___ `'.     .' .'      '.\
#  ||                .   |     \  '         ' |--.\  \   / .'
#  ||        .-,.--. |   '      |  '        | |    \  ' . '                                     .|
#  ||  __    |  .-. |\    \     / /  __     | |     |  '| |                 __                .' |_
#  ||/'__ '. | |  | | `.   ` ..' /.:--.'.   | |     |  || |              .:--.'.         _  .'     |
#  |:/`  '. '| |  | |    '-...-'`/ |   \ |  | |     ' .'. '             / |   \ |      .' |'--.  .-'
#  ||     | || |  '-             `" __ | |  | |___.' /'  \ '.          .`" __ | |     .   | / |  |
#  ||\    / '| |                  .'.''| | /_______.'/    '. `._____.-'/ .'.''| |   .'.'| |// |  |
#  |/\'..' / | |                 / /   | |_\_______|/       `-.______ / / /   | |_.'.'.-'  /  |  '.'
#  '  `'-'`  |_|                 \ \._,\ '/                          `  \ \._,\ '/.'   \_.'   |   /
#                                 `--'  `"                               `--'  `"             `'-'

using Base.Broadcast
using Base.Broadcast: Broadcasted, AbstractArrayStyle, DefaultArrayStyle, broadcasted,
  instantiate, materialize, flatten, combine_eltypes, _broadcast_getindex
using ForwardDiff: Dual

trim(x, Δ) = reshape(Δ, ntuple(i -> size(Δ, i), Val(ndims(x))))

unbroadcast(x::AbstractArray, Δ) =
  size(x) == size(Δ) ? Δ :
  length(x) == length(Δ) ? trim(x, Δ) :
    trim(x, sum(Δ, dims = ntuple(i -> size(x, i) == 1 ? i : ndims(Δ)+1, Val(ndims(Δ)))))

unbroadcast(x::Number, Δ) = sum(Δ)

dual(x, p) = x
dual(x::Real, p) = Dual(x, p)

dualtype(::Type{Dual{G,T,P}}) where {G,T,P} = T

function dual_function(f::F) where F
  function (args::Vararg{Any,N}) where N
    ds = map(args, ntuple(identity,Val(N))) do x, i
      dual(x, ntuple(j -> i==j, Val(N)))
    end
    return f(ds...)
  end
end

dualify(bc::Broadcasted{S}) where S = Broadcasted{S}(dual_function(bc.f), bc.args, bc.axes)

function broadcast_gradient!(bc::Broadcasted, dest::AbstractArray, grads::Vararg{Any})
  @simd for I in eachindex(bc)
    @inbounds begin
      out = bc[I]
      dest[I] = out.value
      map((g, p) -> g[I] = p, grads, out.partials.values)
    end
  end
end

function broadcast_gradient(bc::Broadcasted, ::Type{T}) where T
  dest = similar(bc, T)
  grads = map(_ -> similar(bc, T), bc.args)
  broadcast_gradient!(bc, dest, grads...)
  return dest, grads
end

@inline function ∇broadcast(bc′::Broadcasted) where {F,N}
  bc = dualify(instantiate(flatten(bc′)))
  T = combine_eltypes(bc.f, bc.args)
  y, gs = broadcast_gradient(bc, dualtype(T))
  back(Δ) = map((x, d) -> unbroadcast(x, Δ.*d), bc.args, gs)
  return y, back
end

function ∇broadcast(bc::Broadcasted{<:AbstractArrayStyle{0}})
  out = dualify(instantiate(flatten(bc)))[]
  return out.value, Δ -> map(x -> x*Δ, out.partials.values)
end

using Base: tail

_unflatten(x, xs) = first(xs), tail(xs)

_unflatten(x::Tuple{}, xs) = (), xs

function _unflatten(x::Tuple, xs)
  t1, xs1 = _unflatten(first(x), xs)
  t2, xs2 = _unflatten(tail(x), xs1)
  (t1, t2...), xs2
end

function _unflatten(bc::Broadcasted, xs)
  t, xs′ = _unflatten(bc.args, xs)
  (args=t,f=nothing,axes=nothing), xs′
end

unflatten(x, xs) = _unflatten(x, xs)[1]

@grad function broadcasted(f, args...)
  broadcasted(f, args...), Δ -> (nothing, Δ.args...)
end

@grad function materialize(bc::Broadcasted{<:DefaultArrayStyle})
  let (y, back) = ∇broadcast(bc)
    y, Δ -> (unflatten(bc, back(Δ)),)
  end
end


text:
include("activation.jl")

abstract Layer

"""
a perceptron layer with any number
"""
type PerceptronLayer <: Layer
    weights::Matrix{Float64}
    outputs::Vector{Float64}
    inputs::Vector{Float64}
    bias::Vector{Float64}
    activation::Function

    weights_delta_acc::Matrix{Float64}
    bias_delta_acc::Vector{Float64}

    PerceptronLayer(input_num::Integer, output_num::Integer,
    activation::Function=sigmoid) = begin
        # weights * input = output, therefore the row should be output nodes
        weights = rand(output_num, input_num)

        # normalize for the weights
        for i = 1:output_num
            weights[i, :] ./= sum(weights[i, :])
        end

        new(weights, zeros(output_num), zeros(input_num), rand(output_num),
        activation, zeros(output_num, input_num), zeros(output_num))
    end
end

"""
update parameters of a layer using errors from mini-batch
"""
function batch_update(layer::PerceptronLayer; eta=0.03)
    layer.weights -= eta * layer.weights_delta_acc
    layer.bias -= eta * layer.bias_delta_acc
    layer.weights_delta_acc = zeros(layer.weights_delta_acc)
    layer.bias_delta_acc = zeros(layer.bias_delta_acc)
end

"""
Forward computing in PerceptronLayer
"""
function forward(layer::PerceptronLayer, inputs)
    layer.inputs = inputs
    net = layer.weights * layer.inputs + layer.bias
    layer.outputs = layer.activation(net)
end

"""
Backward error update for PerceptronLayer
"""
function backward(layer::PerceptronLayer, error; eta=0.03)
    # partial error w.r.t. activation function
    activation_derivatives = activation_partial(layer)
    output_num, input_num = size(layer.weights)

    # partial error w.r.t. weighted net sum
    net_error = error .* activation_derivatives

    # build new error to propagate backward
    new_error = [sum([net_error[j] * layer.weights[j, i] for j = 1:output_num])
    for i = 1:input_num ]

    # update weights
    for i = 1:input_num, j = 1:output_num
        gradient = net_error[j] * layer.inputs[i]
        layer.weights_delta_acc[j, i] += gradient
    end

    # update bias
    layer.bias_delta_acc += net_error

    new_error
end

"""
Partial derivative value of activation functions w.r.t. layer_weight
"""
function activation_partial(layer::PerceptronLayer)
    map(layer.outputs) do z
        if layer.activation == sigmoid
            partial_sigmoid(z)
        elseif layer.activation == tanh
            partial_tanh(z)
        else
            warn("partial for activation not implemented")
            z
        end
    end
end

max_stars_repo_name: zxteloiv/julia-ex

max_issues_repo_name: zxteloiv/julia-ex

max_forks_repo_name: zxteloiv/julia-ex


text:
@doc raw"""
`SPDAE`: Split Partitioned Differential Algebraic Equation *EXPERIMENTAL*

Defines a split differential algebraic initial value problem, that is
a canonical Hamiltonian system of equations subject to Dirac constraints,
```math
\begin{aligned}
\dot{q} (t) &= v_1(t, q(t), p(t)) + v_2(t, q(t), p(t), \lambda(t)) + v_3(t, q(t), p(t), \lambda(t), \gamma(t)) , & q(t_{0}) &= q_{0} , \\
\dot{p} (t) &= f_1(t, q(t), p(t)) + f_2(t, q(t), p(t), \lambda(t)) + f_3(t, q(t), p(t), \lambda(t), \gamma(t)) , & p(t_{0}) &= p_{0} , \\
0 &= \phi (t, q(t), p(t)) , \\
0 &= \psi (t, q(t), p(t), \dot{q}(t), \dot{p}(t)) ,
\end{aligned}
```
with vector fields ``v_i`` and ``f_i`` for ``i = 1 ... 3``,
primary constraint ``\phi(q,p)=0`` and secondary constraint ``\psi(q,p,\lambda)=0``,
initial conditions ``(q_{0}, p_{0})``, the dynamical variables ``(q,p)``
taking values in ``\mathbb{R}^{d} \times \mathbb{R}^{d}`` and
the algebraic variables ``(\lambda, \gamma)`` taking values in
``\mathbb{R}^{n} \times \mathbb{R}^{d}``.

### Parameters

* `DT <: Number`: data type
* `TT <: Real`: time step type
* `AT <: AbstractArray{DT}`: array type
* `vType <: Function`: type of `v`
* `fType <: Function`: type of `f`
* `ϕType <: Function`: type of `ϕ`
* `ψType <: Function`: type of `ψ`
* `invType <: OptionalNamedTuple`: invariants type
* `parType <: OptionalNamedTuple`: parameters type
* `perType <: OptionalArray{AT}`: periodicity type

### Fields

* `d`: dimension of dynamical variables ``q`` and ``p`` as well as the vector fields ``v`` and ``f``
* `m`: dimension of algebraic variables ``\lambda`` and ``\gamma`` and the constraints ``\phi`` and ``\psi``
* `v`: tuple of functions computing the vector fields ``v_i``, ``i = 1 ... 3``
* `f`: tuple of functions computing the vector fields ``f_i``, ``i = 1 ... 3``
* `ϕ`: primary constraints
* `ψ`: secondary constraints
* `t₀`: initial time
* `q₀`: initial condition for dynamical variable ``q``
* `p₀`: initial condition for dynamical variable ``p``
* `λ₀`: initial condition for algebraic variable ``λ``
* `μ₀`: initial condition for algebraic variable ``μ`` (optional)
* `invariants`: either a `NamedTuple` containing the equation's invariants or `nothing`
* `parameters`: either a `NamedTuple` containing the equation's parameters or `nothing`
* `periodicity`: determines the periodicity of the state vector `q` for cutting periodic solutions

### Constructors

```julia
SPDAE(v, f, ϕ, ψ, t₀, q₀::StateVector, p₀::StateVector, λ₀::StateVector; kwargs...)
SPDAE(v, f, ϕ, ψ, q₀::StateVector, p₀::StateVector, λ₀::StateVector; kwargs...)
SPDAE(v, f, ϕ, ψ, t₀, q₀::State, p₀::State, λ₀::State=zero(q₀); kwargs...)
SPDAE(v, f, ϕ, ψ, q₀::State, p₀::State, λ₀::State=zero(q₀); kwargs...)
```

### Keyword arguments:

* `invariants = nothing`
* `parameters = nothing`
* `periodicity = nothing`

"""
struct SPDAE{dType <: Number, tType <: Real, arrayType <: AbstractArray{dType},
             vType <: Tuple, fType <: Tuple, ϕType <: Function, ψType <: OptionalFunction,
             invType <: OptionalInvariants,
             parType <: OptionalParameters,
             perType <: OptionalPeriodicity} <: AbstractEquationPDAE{invType,parType,perType,ψType}

    d::Int
    m::Int

    v::vType
    f::fType
    ϕ::ϕType
    ψ::ψType

    t₀::tType
    q₀::Vector{arrayType}
    p₀::Vector{arrayType}
    λ₀::Vector{arrayType}
    μ₀::Vector{arrayType}

    invariants::invType
    parameters::parType
    periodicity::perType

    function SPDAE(v::vType, f::fType, ϕ::ϕType, ψ::ψType,
                   t₀::tType, q₀::Vector{arrayType}, p₀::Vector{arrayType}, λ₀::Vector{arrayType}, μ₀::Vector{arrayType},
                   invariants::invType, parameters::parType, periodicity::perType) where {
                        dType <: Number, tType <: Real, arrayType <: AbstractArray{dType},
                        vType <: Tuple, fType <: Tuple, ϕType <: Function, ψType <: Function,
                        invType <: OptionalInvariants,
                        parType <: OptionalParameters,
                        perType <: OptionalPeriodicity}

        d = length(q₀[begin])
        m = length(λ₀[begin])

        @assert 2d ≥ m

        @assert length(q₀) == length(p₀)
        @assert length(λ₀) == length(μ₀)

        @assert all(length(q) == d for q in q₀)
        @assert all(length(p) == d for p in p₀)
        @assert all(length(λ) == m for λ in λ₀)
        @assert all(length(μ) == m for μ in μ₀)

        @assert all([ndims(q) == ndims(p) == ndims(λ) == ndims(μ) for (q,p,λ,μ) in zip(q₀,p₀,λ₀,μ₀)])

        new{dType, tType, arrayType, vType, fType, ϕType, ψType, invType, parType, perType}(d, m, v, f, ϕ, ψ, t₀, q₀, p₀, λ₀, μ₀, invariants, parameters, periodicity)
    end
end

_SPDAE(v, f, ϕ, ψ, t₀, q₀, p₀, λ₀, μ₀; invariants=NullInvariants(), parameters=NullParameters(), periodicity=NullPeriodicity()) = SPDAE(v, f, ϕ, ψ, t₀, q₀, p₀, λ₀, μ₀, invariants, parameters, periodicity)

SPDAE(v, f, ϕ, ψ, t₀, q₀::StateVector, p₀::StateVector, λ₀::StateVector, μ₀::StateVector=zero(λ₀); kwargs...) = _SPDAE(v, f, ϕ, ψ, t₀, q₀, p₀, λ₀, μ₀; kwargs...)
SPDAE(v, f, ϕ, ψ, q₀::StateVector, p₀::StateVector, λ₀::StateVector, μ₀::StateVector=zero(λ₀); kwargs...) = SPDAE(v, f, ϕ, ψ, 0.0, q₀, p₀, λ₀, μ₀; kwargs...)
SPDAE(v, f, ϕ, ψ, t₀, q₀::State, p₀::State, λ₀::State, μ₀::State=zero(λ₀); kwargs...) = SPDAE(v, f, ϕ, ψ, t₀, [q₀], [p₀], [λ₀], [μ₀]; kwargs...)
SPDAE(v, f, ϕ, ψ, q₀::State, p₀::State, λ₀::State, μ₀::State=zero(λ₀); kwargs...) = SPDAE(v, f, ϕ, ψ, 0.0, q₀, p₀, λ₀, μ₀; kwargs...)

Base.hash(dae::SPDAE, h::UInt) = hash(dae.d, hash(dae.m,
                    hash(dae.v, hash(dae.f, hash(dae.ϕ, hash(dae.ψ,
                    hash(dae.t₀, hash(dae.q₀, hash(dae.p₀,
                    hash(dae.invariants, hash(dae.parameters, hash(dae.periodicity, h))))))))))))

Base.:(==)(dae1::SPDAE, dae2::SPDAE) = (
                                dae1.d == dae2.d
                             && dae1.m == dae2.m
                             && dae1.v == dae2.v
                             && dae1.f == dae2.f
                             && dae1.ϕ == dae2.ϕ
                             && dae1.ψ == dae2.ψ
                             && dae1.t₀ == dae2.t₀
                             && dae1.q₀ == dae2.q₀
                             && dae1.p₀ == dae2.p₀
                             && dae1.λ₀ == dae2.λ₀
                             && dae1.μ₀ == dae2.μ₀
                             && dae1.invariants  == dae2.invariants
                             && dae1.parameters  == dae2.parameters
                             && dae1.periodicity == dae2.periodicity)

function Base.similar(equ::SPDAE, t₀::Real, q₀::StateVector, p₀::StateVector, λ₀::StateVector, μ₀::StateVector=initial_multiplier(q₀, equ.μ₀); parameters=equ.parameters)
    @assert all([length(q) == equ.d for q in q₀])
    @assert all([length(p) == equ.d for p in p₀])
    @assert all([length(λ) == equ.m for λ in λ₀])
    @assert all([length(μ) == equ.m for μ in μ₀])
    _SPDAE(equ.v, equ.f, equ.ϕ, equ.ψ, t₀, q₀, p₀, λ₀, μ₀;
          invariants=equ.invariants, parameters=parameters, periodicity=equ.periodicity)
end

Base.similar(equ::SPDAE, q₀, p₀, λ₀=initial_multiplier(q₀, equ.λ₀), μ₀=initial_multiplier(q₀, equ.μ₀); kwargs...) = similar(equ, equ.t₀, q₀, p₀, λ₀, μ₀; kwargs...)
Base.similar(equ::SPDAE, t₀::Real, q₀::State, p₀::State, λ₀::State=initial_multiplier(q₀, equ.λ₀), μ₀::State=initial_multiplier(λ₀, equ.μ₀); kwargs...) = similar(equ, equ.t₀, [q₀], [p₀], [λ₀], [μ₀]; kwargs...)

@inline Base.axes(equation::SPDAE) = axes(equation.q₀[begin])
@inline Base.ndims(equation::SPDAE) = equation.d
@inline GeometricBase.nsamples(equation::SPDAE) = length(eachindex(equation.q₀))
@inline GeometricBase.nconstraints(equation::SPDAE) = equation.m

@inline GeometricBase.periodicity(equation::SPDAE) = hasperiodicity(equation) ? equation.periodicity : zero(equation.q₀[begin])
@inline initial_conditions(equation::SPDAE) = (equation.t₀, equation.q₀, equation.p₀, equation.λ₀, equation.μ₀)

function functions(equation::SPDAE{DT,TT,AT,VT,FT,ϕT,ψT,IT,PT}) where {DT, TT, AT, VT, FT, ϕT, ψT, IT, PT <: NullParameters}
    NamedTuple{(:v, :f, :ϕ, :ψ)}((equation.v, equation.f, equation.ϕ, equation.ψ))
end

function functions(equation::SPDAE{DT,TT,AT,VT,FT,ϕT,ψT,IT,PT}) where {DT, TT, AT, VT, FT, ϕT, ψT, IT, PT <: NamedTuple}
    vₚ = (t,q,p,v)   -> equation.v(t, q, p, v, equation.parameters)
    fₚ = (t,q,p,f)   -> equation.f(t, q, p, f, equation.parameters)
    ϕₚ = (t,q,p,ϕ)   -> equation.ϕ(t, q, p, ϕ, equation.parameters)
    ψₚ = (t,q,p,v,f,ψ) -> equation.ψ(t, q, p, v, f, ψ, equation.parameters)

    names = (:v, :f, :ϕ, :ψ)
    equs  = (vₚ, fₚ, ϕₚ, ψₚ)

    NamedTuple{names}(equs)
end


text:
"""
**Adam Optimizer**
```julia
    Adam(;α=0.001, β₁=0.9, β₂=0.999, ϵ=10e-8)
```

Algorithm:
```math
\\begin{align*}
    m_t =& \\beta_1 m_{t-1} + (1-\\beta_1)g_t\\\\
    v_t =& \\beta_2 v_{t-1} + (1-\\beta_2)g_t^2\\\\
    \\hat{m}_t =& \\frac{m_t}{1-\\beta_1^t}\\\\
    \\hat{v}_t =& \\frac{v_t}{1-\\beta_2^t}\\\\
    \\Delta x_t =& \\frac{\\alpha}{\\sqrt{\\hat{v}_t}+\\epsilon}\\hat{m}_t\\\\
\\end{align*}
```
[Algorithm Reference](https://arxiv.org/abs/1412.6980)
"""
mutable struct Adam <: Optimizer
    opt_type::String
    t::Int64
    ϵ::Float64
    α::Float64
    β₁::Float64
    β₂::Float64
    m_t::AbstractArray
    v_t::AbstractArray
end

function Adam(;α::Real=0.001, β₁::Real=0.9, β₂::Real=0.999, ϵ::Real=10e-8)
    @assert α > 0.0 "α must be greater than 0"
    @assert β₁ > 0.0 "β₁ must be greater than 0"
    @assert β₂ > 0.0 "β₂ must be greater than 0"
    @assert ϵ > 0.0 "ϵ must be greater than 0"

    Adam("Adam", 0, ϵ, α, β₁, β₂, [], [])
end

params(opt::Adam) = "ϵ=$(opt.ϵ), α=$(opt.α), β₁=$(opt.β₁), β₂=$(opt.β₂)"

function update(opt::Adam, g_t::AbstractArray{T}) where {T<:Real}
    # resize biased moment estimates if first iteration
    if opt.t == 0
        opt.m_t = zero(g_t)
        opt.v_t = zero(g_t)
    end

    # update timestep
    opt.t += 1

    # update biased first moment estimate
    opt.m_t = opt.β₁ * opt.m_t + (one(T) - opt.β₁) * g_t

    # update biased second raw moment estimate
    opt.v_t = opt.β₂ * opt.v_t + (one(T) - opt.β₂) * ((g_t) .^2)

    # compute bias corrected first moment estimate
    m̂_t = opt.m_t / (one(T) - opt.β₁^opt.t)

    # compute bias corrected second raw moment estimate
    v̂_t = opt.v_t / (one(T) - opt.β₂^opt.t)

    # apply update
    ρ = opt.α * m̂_t ./ (sqrt.(v̂_t .+ opt.ϵ))

    return ρ
end


text:
function gwr(data,MAXNUMBEROFNODES)

#the initial parameters for the algorithm:
global maxnodes, at, en, eb, h0, ab, an, tb, tn, amax, A, C

maxnodes = MAXNUMBEROFNODES; #maximum number of nodes/neurons in the gas
at = 0.95; #activity threshold
en = 0.006; #epsilon subscript n
eb = 0.2; #epsilon subscript b
h0 = 1;
ab = 0.95;
an = 0.95;
tb = 3.33;
tn = 3.33;
amax = 50; #greatest allowed age

t0 = time();
STATIC = true;
RANDOMSTART = true;
DOOVER = 1;

#test some algorithm conditions:
if !(0 < en || en < eb || eb < 1)
    error("en and/or eb definitions are wrong. They are as: 0<en<eb<1.")
end
# (1)
# pick n1 and n2 from data
ni1 = 1;
ni2 = 2;
if RANDOMSTART
    n = randperm(size(data,2));
    ni1 = n[1];
    ni2 = n[2];
end
n1 = data[:,ni1]; n2 = data[:,ni2];

#println(strcat("Initial parameters ","n1 =  ",num2str(ni1)," n2 =  ",num2str(ni2)))

A = zeros(size(n1,1),maxnodes);
    A[:,[1,2]] = hcat(n1, n2)

# (2)
# initialize empty set C

C = spzeros(maxnodes,maxnodes);
C_age = C;

r = 3;
h = zeros(1,maxnodes);
datasetsize = size(data,2);

#some variables to display the graphs
activations = [];
nodecount = [];
epoch = 1;

### SPEEDUP CHANGE
if STATIC
    hizero = hi(0)*ones(1,maxnodes);
    hszero = hs(0);
else
    timetime = 0;
end
OldA = A;
# start of the loop
for epochscounter = 1:DOOVER
for k = 1:datasetsize #step 1
	#global A
        eta = data[:,k]; # this the k-th data sample
        ws, (), s, t, () = findnearest(eta, A); #step 2 and 3
    if C[s,t]==0 #step 4
        C = spdi_bind(C,s,t);
    else
        C_age = spdi_del(C_age,s,t);
    end
    a = exp(-norm(eta-ws)); #step 5

    neighbours = findneighbours(s, C);

    num_of_neighbours = size(neighbours,1);

    if a < at && r <= maxnodes #step 6
        wr = 0.5*(ws+eta); #too low activity, needs to create new node r
        A[:,r] = wr;
        C = spdi_bind(C,t,r);
        C = spdi_bind(C,s,r);
        C = spdi_del(C,s,t);
        r = r+1;
    else #step 7
        for j = 1:num_of_neighbours # check this for possible indexing errors
            i = neighbours[j];
            wi = A[:,i];
            A[:,i] = wi + en*h[i]*(eta-wi);
        end
        A[:,s] = ws + eb*h[s]*(eta-ws);
    end
    #step 8 : age edges with end at s
    #first we need to find if the edges connect to s

    for j = 1:num_of_neighbours # check this for possible indexing errors
           i = neighbours[j];
           C_age = spdi_add(C_age,s,i);
    end

    #step 9: again we do it inverted, for loop first
    #### this strange check is a speedup for the case when the algorithm is static
    if STATIC # skips this if algorithm is static
        h = hizero;
        h[s] = hszero;
    else
        for i = 1:r
            h[i] = hi(time);
        end
        h[s] = hs(timetime);
        timetime = (time - t0)*1;
    end

    #step 10: check if a node has no edges and delete them
    if r>2 # don't remove everything
        C, C_age  = removeedge(C, C_age);
        C, A, C_age, h, r  = removenode(C, A, C_age, h, r);  #inverted order as it says on the algorithm to remove points faster
    end
    #activations = vcat(activations, a);
    epoch += 1;

end
end
return A,C,ni1,ni2
end

function spdi_add(sparsemat, a, b) #increases the number
sparsemat[a,b] = sparsemat[a,b] + 1;
sparsemat[b,a] = sparsemat[a,b] + 1;
return sparsemat
end

function spdi_bind(sparsemat, a, b) # adds a 2 way connection
sparsemat[a,b] = 1;
sparsemat[b,a] = 1;
return sparsemat
end

function spdi_del(sparsemat, a, b) # removes a 2 way connection
sparsemat[a,b] = 0;
sparsemat[b,a] = 0;
return sparsemat
end

function S(t)
    X = 1;
return X
end

function hs(t)
global h0,ab,tb
X = h0 - S(t)/ab*(1-exp(-ab*t/tb));
return X
end

function hi(t)
global h0, an, tn
X = h0 - S(t)/an*(1-exp(-an*t/tn));
return X
end


text:
# Window2 ==========================================================
"""Windows2(y,t,nw,noverlap,window_func)
noverlap = -1 sets the noverlap to dpw/2
"""
struct Windows2
    y::AbstractVector
    t::AbstractVector
    nw::Int
    dpw::Int
    noverlap::Int
    W
    function Windows2(y,t,nw,noverlap,window_func)
        N   = length(y)
        @assert N == length(t) "y and t has to be the same length"
        dpw = floor(Int64,N/nw)
        if noverlap == -1
            noverlap = round(Int64,dpw/2)
        end
        W = window_func(dpw)
        new(y,t,nw,dpw,noverlap,W)
    end
end

function Base.iterate(w::Windows2, state=0)
    state >= w.nw && return nothing
    if state < w.nw
        inds =  (1:w.dpw) .+ state*(w.dpw-w.noverlap)
    else
        inds =  state*(w.dpw-w.noverlap):length(w.y)
    end
    ((w.y[inds],w.t[inds]), state+1)
end

Base.length(w::Windows2) = w.nw;

# Window3 ==========================================================
"""Windows3(y,t,v,nw::Integer,noverlap::Integer,window_func::Function)
noverlap = -1 sets the noverlap to dpw/2
"""
struct Windows3
    y::AbstractVector
    t::AbstractVector
    v::AbstractVector
    nw::Int
    dpw::Int
    noverlap::Int
    W
end

function Windows3(y,t,v,nw::Integer,noverlap::Integer,window_func::Function)
    N       = length(y)
    @assert N == length(t) == length(v) "y, t and v has to be the same length"
    dpw     = floor(Int64,N/nw)
    if noverlap == -1
        noverlap = round(Int64,dpw/2)
    end
    W       = window_func(dpw)
    Windows3(y,t,v,nw,dpw,noverlap,W)
end

function Base.iterate(w::Windows3, state=0)
    state >= w.nw && return nothing
    if state < w.nw
        inds =  (1:w.dpw) .+ state*(w.dpw-w.noverlap)
    else
        inds =  state*(w.dpw-w.noverlap):length(w.y)
    end
    ((w.y[inds],w.t[inds],w.v[inds]), state+1)
end

Base.length(w::Windows3) = w.nw;


text:
"""
    struct IdentifierTable <: Any
Hold a pointer to a `clang::IdentifierTable` object.
"""
struct IdentifierTable
    ptr::CXIdentifierTable
end

Base.unsafe_convert(::Type{CXIdentifierTable}, x::IdentifierTable) = x.ptr
Base.cconvert(::Type{CXIdentifierTable}, x::IdentifierTable) = x

"""
    struct IdentifierInfo <: Any
Hold a pointer to a `clang::IdentifierInfo` object.
"""
struct IdentifierInfo
    ptr::CXIdentifierInfo
end

Base.unsafe_convert(::Type{CXIdentifierInfo}, x::IdentifierInfo) = x.ptr
Base.cconvert(::Type{CXIdentifierInfo}, x::IdentifierInfo) = x


text:
""" Backend-agnostic layout of wiring diagrams via morphism expressions.

This module lays out wiring diagrams for visualization, independent of any
specific graphics system. It uses the structure of a morphism expression to
determine the layout. Thus, the first step of the algorithm is to convert the
wiring diagram to a symbolic expression, using the submodule
`WiringDiagrams.Expressions`. Morphism expressions may also be given directly.
"""
module WiringDiagramLayouts
export LayoutOrientation, LeftToRight, RightToLeft, TopToBottom, BottomToTop,
  layout_diagram, layout_box

using Compat
import Base: sign
using LinearAlgebra: dot, normalize
using Logging
using Parameters: @with_kw_noshow
using StaticArrays: StaticVector, SVector
using Statistics: mean

using ...Syntax
using ...Theories: ObExpr, HomExpr, dom, codom, compose, id, otimes, braid
using ...WiringDiagrams

# Data types
############

const AbstractVector2D = StaticVector{2,<:Real}
const Vector2D = SVector{2,Float64}

""" Orientation of wiring diagram.
"""
@enum LayoutOrientation LeftToRight RightToLeft TopToBottom BottomToTop

is_horizontal(orient::LayoutOrientation) = orient in (LeftToRight, RightToLeft)
is_vertical(orient::LayoutOrientation) = orient in (TopToBottom, BottomToTop)
is_positive(orient::LayoutOrientation) = orient in (LeftToRight, TopToBottom)
is_negative(orient::LayoutOrientation) = orient in (RightToLeft, BottomToTop)

sign(orient::LayoutOrientation) = is_positive(orient) ? +1 : -1
port_sign(kind::PortKind, orient::LayoutOrientation) =
  (kind == InputPort ? -1 : +1) * sign(orient)
port_sign(d::WiringDiagram, port::Port, orient::LayoutOrientation) =
  (port.box in outer_ids(d) ? -1 : +1) * port_sign(port.kind, orient)

svector(orient::LayoutOrientation, first, second) =
  is_horizontal(orient) ? SVector(first, second) : SVector(second, first)

""" Internal data type for configurable options of wiring diagram layout.
"""
@with_kw_noshow struct LayoutOptions
  orientation::LayoutOrientation = LeftToRight
  default_box_shape::Symbol = :rectangle
  box_shapes::AbstractDict = Dict()
  box_styles::AbstractDict = Dict()
  outer_ports_layout::Symbol = :isotonic
  anchor_wires::Union{Bool,AbstractSet,AbstractVector} = [:id,:braid]
  base_box_size::Real = 2
  box_stretch::Real = 1
  sequence_pad::Real = 2
  parallel_pad::Real = 1
  junction_size::Real = 0.25
end

svector(opts::LayoutOptions, args...) = svector(opts.orientation, args...)

""" Layout for box in a wiring diagram.

Suggests a shape, size, and position for the box. Also includes the box's
value, typically used for labels, and a style name, which does not affect the
layout but is interpreted by graphics backends to adjust the visual style.

Specific features of the shape are determined by the graphics backend. For
example, a rectangle could be rendered with or without rounded corners or even
as another, similar shape, such as a parallelogram.
"""
@with_kw_noshow mutable struct BoxLayout
  value::Any = nothing
  shape::Symbol = :rectangle
  position::Vector2D = zeros(Vector2D)
  size::Vector2D = zeros(Vector2D)
  style::Union{Symbol,Nothing} = nothing
end

position(layout::BoxLayout) = layout.position
size(layout::BoxLayout) = layout.size
lower_corner(layout::BoxLayout) = position(layout) - size(layout)/2
upper_corner(layout::BoxLayout) = position(layout) + size(layout)/2

position(box::AbstractBox) = position(box.value)
size(box::AbstractBox) = size(box.value)
lower_corner(box::AbstractBox) = lower_corner(box.value)
upper_corner(box::AbstractBox) = upper_corner(box.value)

contents_lower_corner(diagram::WiringDiagram) =
  mapreduce(lower_corner, (c,d) -> min.(c,d), boxes(diagram))
contents_upper_corner(diagram::WiringDiagram) =
  mapreduce(upper_corner, (c,d) -> max.(c,d), boxes(diagram))

""" Layout for port in a wiring diagram.
"""
@with_kw_noshow struct PortLayout
  value::Any = nothing
  position::Vector2D = zeros(Vector2D) # Position relative to box center.
  normal::Vector2D = zeros(Vector2D)   # Outward unit normal vector.
  label_wires::Bool = true             # Show labels for wires into/out of port?
  reverse_wires::Bool = false          # Reverse wires into/out of port?
end

position(layout::PortLayout) = layout.position
normal(layout::PortLayout) = layout.normal

""" Interior point of wire in a wiring diagram.
"""
mutable struct WirePoint
  position::Vector2D     # Position of wire at point.
  tangent::Vector2D      # Unit tangent vector along wire at point.
end

position(point::WirePoint) = point.position
tangent(point::WirePoint) = point.tangent

const WireLayout = Union{Vector{WirePoint},Nothing}

wire_points(points::Vector{WirePoint}) = points
wire_points(::Nothing) = WirePoint[]
wire_points(wire::Wire) = wire_points(wire.value)

# Main entry point
##################

""" Lay out a wiring diagram or morphism expression for visualization.

If a wiring diagram is given, it is first to converted to a morphism expression.

The layout is calculated with respect to a cartesian coordinate system with
origin at the center of the diagram and the positive y-axis pointing
*downwards*. Box positions are relative to their centers. All positions and
sizes are dimensionless (unitless).
"""
function layout_diagram(Syntax::Module, diagram::WiringDiagram; kw...)
  layout_diagram(to_hom_expr(Syntax, diagram); kw...)
end
function layout_diagram(Ob::Type, Hom::Type, diagram::WiringDiagram; kw...)
  layout_diagram(to_hom_expr(Ob, Hom, diagram); kw...)
end

function layout_diagram(expr::HomExpr; kw...)::WiringDiagram
  opts = LayoutOptions(; kw...)
  diagram = layout_hom_expr(expr, opts)
  layout_outer_ports!(diagram, opts, method=opts.outer_ports_layout)
end

# Morphism expression layout
############################

""" Lay out a morphism expression as a wiring diagram.
"""
layout_hom_expr(f::HomExpr, opts) = layout_box(f, opts)

layout_hom_expr(f::HomExpr{:compose}, opts) =
  compose_with_layout!(map(arg -> layout_hom_expr(arg, opts), args(f)), opts)
layout_hom_expr(f::Union{HomExpr{:otimes},HomExpr{:oplus}}, opts) =
  otimes_with_layout!(map(arg -> layout_hom_expr(arg, opts), args(f)), opts)

layout_hom_expr(f::HomExpr{:id}, opts) = layout_pure_wiring(f, opts)
layout_hom_expr(f::Union{HomExpr{:braid},HomExpr{:swap}}, opts) =
  layout_pure_wiring(f, opts)

layout_hom_expr(f::HomExpr{:mcopy}, opts) = layout_supply(f, opts, shape=:junction)
layout_hom_expr(f::HomExpr{:delete}, opts) = layout_supply(f, opts, shape=:junction)
layout_hom_expr(f::HomExpr{:mmerge}, opts) = layout_supply(f, opts, shape=:junction)
layout_hom_expr(f::HomExpr{:create}, opts) = layout_supply(f, opts, shape=:junction)

layout_hom_expr(f::HomExpr{:plus}, opts) =
  layout_supply(f, opts, shape=:junction, style=:variant_junction)
layout_hom_expr(f::HomExpr{:zero}, opts) =
  layout_supply(f, opts, shape=:junction, style=:variant_junction)
layout_hom_expr(f::HomExpr{:coplus}, opts) =
  layout_supply(f, opts, shape=:junction, style=:variant_junction)
layout_hom_expr(f::HomExpr{:cozero}, opts) =
  layout_supply(f, opts, shape=:junction, style=:variant_junction)

layout_hom_expr(f::HomExpr{:dunit}, opts) =
  layout_supply(f, opts, shape=:junction, visible=false, pad=false)
layout_hom_expr(f::HomExpr{:dcounit}, opts) =
  layout_supply(f, opts, shape=:junction, visible=false, pad=false)

layout_port(A::ObExpr{:dual}; kw...) =
  PortLayout(; value=A, reverse_wires=true, kw...)
wire_label(mime::MIME, A::ObExpr{:dual}) = wire_label(mime, first(A))

function layout_box(f::HomExpr, opts::LayoutOptions; kw...)
  layout_box(collect(dom(f)), collect(codom(f)), opts; value=f, kw...)
end

function layout_pure_wiring(f::HomExpr, opts::LayoutOptions; kw...)
  layout_pure_wiring(to_wiring_diagram(f, identity, identity), opts;
    anchor_wires=head(f) in opts.anchor_wires, kw...)
end

function layout_supply(f::HomExpr, opts::LayoutOptions;
                       ob::Union{ObExpr,Nothing}=nothing, kw...)
  obs = collect(isnothing(ob) ? first(f)::ObExpr : ob)
  n = length(obs)
  ndom, ncodom = ndims(dom(f)) ÷ n, ndims(codom(f)) ÷ n
  diagram = otimes_with_layout!(map(obs) do ob
    layout_box(repeat([ob], ndom), repeat([ob], ncodom), opts; kw...)
  end, opts)
  anchor_wires = head(f) in opts.anchor_wires
  if n > 1 && ndom > 1
    σ = reshape(permutedims(reshape(1:ndom*n, n, :)), :)
    σ_diagram = permute(dom(diagram), σ, inverse=true)
    σ_diagram = layout_pure_wiring(σ_diagram, opts, anchor_wires=anchor_wires)
    diagram = compose_with_layout!(σ_diagram, diagram, opts)
  end
  if n > 1 && ncodom > 1
    σ = reshape(permutedims(reshape(1:ncodom*n, n, :)), :)
    σ_diagram = permute(codom(diagram), σ)
    σ_diagram = layout_pure_wiring(σ_diagram, opts, anchor_wires=anchor_wires)
    diagram = compose_with_layout!(diagram, σ_diagram, opts)
  end
  diagram
end

# Diagram layout
################

""" Compose wiring diagrams and their layouts.

Compare with: `WiringDiagram.compose`.
"""
function compose_with_layout!(d1::WiringDiagram, d2::WiringDiagram, opts::LayoutOptions)
  diagram = compose(d1, d2; unsubstituted=true)
  dir = svector(opts, sign(opts.orientation), 0)
  place_adjacent!(d1, d2; dir=dir, pad=-opts.sequence_pad)
  substitute_with_layout!(size_to_fit!(diagram, opts, pad=false))
end
function compose_with_layout!(diagrams::Vector{<:WiringDiagram}, opts::LayoutOptions)
  foldl((d1,d2) -> compose_with_layout!(d1, d2, opts), diagrams)
end

""" Tensor wiring diagrams and their layouts.

Compare with: `WiringDiagram.otimes`.
"""
function otimes_with_layout!(d1::WiringDiagram, d2::WiringDiagram, opts::LayoutOptions)
  # Compare with `WiringDiagrams.otimes`.
  diagram = otimes(d1, d2; unsubstituted=true)
  dir = svector(opts, 0, 1)
  place_adjacent!(d1, d2; dir=dir, pad=-opts.parallel_pad)
  substitute_with_layout!(size_to_fit!(diagram, opts, pad=false))
end
function otimes_with_layout!(diagrams::Vector{<:WiringDiagram}, opts::LayoutOptions)
  foldl((d1,d2) -> otimes_with_layout!(d1, d2, opts), diagrams)
end

""" Size a wiring diagram to fit its contents.

The inner boxes are also shifted to be centered within the new bounds.
"""
function size_to_fit!(diagram::WiringDiagram, opts::LayoutOptions; pad::Bool=true)
  nin, nout = length(input_ports(diagram)), length(output_ports(diagram))
  minimum_size = minimum_diagram_size(nin, nout, opts)

  lower, upper = contents_lower_corner(diagram), contents_upper_corner(diagram)
  content_size = upper - lower
  if pad; content_size += 2 * diagram_padding(opts) end
  size = max.(minimum_size, content_size)

  content_center = (lower + upper)/2
  shift_contents!(diagram, -content_center)
  diagram.value = BoxLayout(size=size)
  diagram
end

""" Substitute sub-wiring diagrams, preserving their layouts.
"""
function substitute_with_layout!(d::WiringDiagram)
  substitute_with_layout!(d, filter(v -> box(d,v) isa WiringDiagram, box_ids(d)))
end
function substitute_with_layout!(d::WiringDiagram, vs::Vector{Int})
  for v in vs
    sub = box(d, v)::WiringDiagram
    shift_contents!(sub, position(sub))
  end
  substitute(d, vs; merge_wire_values=merge_wire_layouts)
end

""" Place one box adjacent to another.

The absolute positions are undefined; only relative positions are guaranteed.
"""
function place_adjacent!(box1::AbstractBox, box2::AbstractBox;
                         dir::AbstractVector2D=SVector(1,0), pad::Real=0)
  layout1, layout2 = box1.value::BoxLayout, box2.value::BoxLayout
  layout1.position = -(pad .+ size(layout1))/2 .* dir
  layout2.position = (pad .+ size(layout2))/2 .* dir
end

""" Shift all boxes and wires within wiring diagram by a fixed offset.
"""
function shift_contents!(diagram::WiringDiagram, offset::AbstractVector2D)
  for box in boxes(diagram)
    layout = box.value::BoxLayout
    layout.position += offset
  end
  for wire in wires(diagram)
    for point in wire_points(wire)
      point.position += offset
    end
  end
  diagram
end

""" Compute the minimum size of a wiring diagram from the number of its ports.
"""
function minimum_diagram_size(nin::Int, nout::Int, opts::LayoutOptions)
  default_box_size(nin, nout, opts, stretch=0) + 2*diagram_padding(opts)
end
function diagram_padding(opts::LayoutOptions)
  svector(opts, opts.sequence_pad, opts.parallel_pad)
end

# Box layout
############

""" Lay out a box and its ports.

By default the box is rectangular, but other shapes are also supported.
"""
function layout_box(inputs::Vector, outputs::Vector, opts::LayoutOptions;
                    shape::Union{Symbol,Nothing}=nothing,
                    style::Union{Symbol,Nothing}=nothing, value=nothing, kw...)
  shape = get(opts.box_shapes, value) do
    isnothing(shape) ? opts.default_box_shape : shape
  end
  style = get(opts.box_styles, value, style)
  layout_box(Val{shape}, inputs, outputs, opts; style=style, value=value, kw...)
end

function layout_box(::Type{Val{:rectangle}}, inputs::Vector, outputs::Vector,
                    opts::LayoutOptions; shape::Symbol=:rectangle, kw...)
  size = default_box_size(length(inputs), length(outputs), opts)
  box = Box(BoxLayout(; shape=shape, size=size, kw...),
            layout_linear_ports(InputPort, inputs, size, opts),
            layout_linear_ports(OutputPort, outputs, size, opts))
  size_to_fit!(singleton_diagram(box), opts)
end

function layout_box(::Type{Val{:circle}}, inputs::Vector, outputs::Vector,
                    opts::LayoutOptions; pad::Bool=true, kw...)
  size = default_box_size(1, 1, opts)
  radius = first(size) / 2
  box = Box(BoxLayout(; shape=:circle, size=size, kw...),
            layout_circular_ports(InputPort, inputs, radius, opts; pad=pad),
            layout_circular_ports(OutputPort, outputs, radius, opts; pad=pad))
  size_to_fit!(singleton_diagram(box), opts)
end

function layout_box(::Type{Val{:ellipse}}, inputs::Vector, outputs::Vector,
                    opts::LayoutOptions; pad::Bool=true, kw...)
  size = default_box_size(length(inputs), length(outputs), opts)
  radii = size / 2
  box = Box(BoxLayout(; shape=:ellipse, size=size, kw...),
            layout_elliptical_ports(InputPort, inputs, radii, opts; pad=pad),
            layout_elliptical_ports(OutputPort, outputs, radii, opts; pad=pad))
  size_to_fit!(singleton_diagram(box), opts)
end

function layout_box(::Type{Val{:triangle}}, inputs::Vector, outputs::Vector,
                    opts::LayoutOptions; kw...)
  @assert length(inputs) <= 1 "Cannot use triangle shape with multiple inputs"
  layout_box(Val{:rectangle}, inputs, outputs, opts; shape=:triangle, kw...)
end
function layout_box(::Type{Val{:invtriangle}}, inputs::Vector, outputs::Vector,
                    opts::LayoutOptions; kw...)
  @assert length(outputs) <= 1 "Cannot use invtriangle shape with multiple outputs"
  layout_box(Val{:rectangle}, inputs, outputs, opts; shape=:invtriangle, kw...)
end

# Although `trapezoid` is the standard term in North American English,
# we use the term `trapezium` because both Graphviz and TikZ do.
function layout_box(::Type{Val{:trapezium}}, inputs::Vector, outputs::Vector,
                    opts::LayoutOptions; kw...)
  layout_box(Val{:rectangle}, inputs, outputs, opts; shape=:trapezium, kw...)
end
function layout_box(::Type{Val{:invtrapezium}}, inputs::Vector, outputs::Vector,
                    opts::LayoutOptions; kw...)
  layout_box(Val{:rectangle}, inputs, outputs, opts; shape=:invtrapezium, kw...)
end

function layout_box(::Type{Val{:junction}}, inputs::Vector, outputs::Vector,
                    opts::LayoutOptions; visible::Bool=true, pad::Bool=true, kw...)
  nin, nout = length(inputs), length(outputs)
  shape, radius = visible ? (:junction, opts.junction_size) : (:invisible, 0)
  size = 2*SVector(radius, radius)
  box = Box(BoxLayout(; shape=shape, size=size, kw...),
            layout_circular_ports(InputPort, inputs, radius, opts;
                                  pad=pad, label_wires = nin == 1 || nout == 0),
            layout_circular_ports(OutputPort, outputs, radius, opts;
                                  pad=pad, label_wires = nin == 0 || nout == 1))
  size_to_fit!(singleton_diagram(box), opts)
end

""" Compute the default size of a rectangular box from the number of its ports.

We use the unique formula consistent with the padding for monoidal products,
ensuring that the size of a product of boxes depends only on the total number of
ports, not on the number of boxes.
"""
function default_box_size(nin::Int, nout::Int, opts::LayoutOptions;
                          stretch::Union{Real,Nothing}=nothing)
  base_size = opts.base_box_size
  stretch = isnothing(stretch) ? opts.box_stretch : stretch
  n = max(1, nin, nout)
  svector(opts, stretch*base_size, n*base_size + (n-1)*opts.parallel_pad)
end

# Port layout
#############

""" Lay out ports linearly, for use with a rectangular box.

Assumes the box is at least as large as `default_box_size()`.
"""
function layout_linear_ports(port_kind::PortKind, port_values::Vector,
                             box_size::AbstractVector2D, opts::LayoutOptions; kw...)
  n = length(port_values)
  coeffs = range(-(n-1), n-1, step=2) / 2 # Always length n
  coords = (opts.base_box_size + opts.parallel_pad) .* collect(coeffs)
  layout_linear_ports(port_kind, port_values, coords, box_size, opts, kw...)
end

function layout_linear_ports(port_kind::PortKind, port_values::Vector,
                             port_coords::Vector, box_size::AbstractVector2D,
                             opts::LayoutOptions; kw...)
  normal_dir = svector(opts, port_sign(port_kind, opts.orientation), 0.0)
  port_dir = svector(opts, 0.0, 1.0)
  start = box_size/2 .* normal_dir
  positions = [ start + coord .* port_dir for coord in port_coords ]
  PortLayout[ layout_port(value, position=pos, normal=normal_dir; kw...)
              for (value, pos) in zip(port_values, positions) ]
end

""" Lay out ports along a circular arc.
"""
function layout_circular_ports(port_kind::PortKind, port_values::Vector,
                               radius::Real, opts::LayoutOptions; kw...)
  radii = SVector(radius, radius)
  layout_elliptical_ports(port_kind, port_values, radii, opts; kw...)
end

""" Lay out ports along an elliptical arc.
"""
function layout_elliptical_ports(port_kind::PortKind, port_values::Vector,
                                 radii::AbstractVector2D, opts::LayoutOptions;
                                 pad::Bool=true, kw...)
  n = length(port_values)
  port_dir = svector(opts, port_sign(port_kind, opts.orientation), 0)
  θ1, θ2 = if port_dir == SVector(1,0); (π/2, -π/2)
    elseif port_dir == SVector(0,1); (-π, 0)
    elseif port_dir == SVector(-1,0); (π/2, 3π/2)
    elseif port_dir == SVector(0,-1); (π, 0) end
  θs = collect(pad ? range(θ1,θ2,length=n+2)[2:n+1] : range(θ1,θ2,length=n))
  dirs = [ SVector(cos(θ),-sin(θ)) for θ in θs ] # positive y-axis downwards
  PortLayout[ layout_port(value, position=radii .* dir, normal=dir; kw...)
              for (value, dir) in zip(port_values, dirs) ]
end

""" Lay out outer ports of wiring diagram.
"""
function layout_outer_ports!(diagram::WiringDiagram, opts::LayoutOptions; kw...)
  original_value = x -> x isa PortLayout ? x.value : x # XXX: This is ugly.
  diagram.input_ports = map(original_value, input_ports(diagram))
  diagram.output_ports = map(original_value, output_ports(diagram))
  diagram.input_ports, diagram.output_ports =
    layout_outer_ports(diagram, opts; kw...)
  diagram
end

""" Lay out outer ports of wiring diagram.
"""
function layout_outer_ports(diagram::WiringDiagram, opts::LayoutOptions;
                            method::Symbol=:fixed, kw...)
  if !has_port_layout_method(method)
    if method == :isotonic
      @warn """
        Both Convex.jl and SCS.jl must be loaded to use isotonic
        regression for port layout. Falling back to fixed layout.
        """ maxlog=1
    else
      error("Unknown port layout method: $method")
    end
    method = :fixed
  end
  layout_outer_ports(diagram, opts, Val{method}; kw...)
end

function layout_outer_ports(diagram::WiringDiagram, opts::LayoutOptions,
                            ::Type{Val{:fixed}})
  (layout_linear_ports(InputPort, input_ports(diagram), size(diagram), opts),
   layout_linear_ports(OutputPort, output_ports(diagram), size(diagram), opts))
end

function layout_outer_ports(diagram::WiringDiagram, opts::LayoutOptions,
                            ::Type{Val{:isotonic}})
  inputs, outputs = input_ports(diagram), output_ports(diagram)
  diagram_size = diagram.value.size
  port_dir = svector(opts, 0, 1)
  port_coord = v -> dot(v, port_dir)
  max_range = port_coord(diagram_size)
  lower, upper = -max_range/2, max_range/2
  pad = opts.base_box_size + opts.parallel_pad

  # Solve optimization problem for input ports.
  # Note that we are minimizing distance only along the port axis, not in the
  # full Euclidean plan.
  ys = [ Float64[] for i in eachindex(inputs) ]
  for wire in out_wires(diagram, input_id(diagram))
    i, layout = wire.source.port, wire.value::WireLayout
    if !(isnothing(layout) || isempty(layout))
      push!(ys[i], port_coord(position(first(layout))))
    elseif wire.target.box != output_id(diagram)
      push!(ys[i], port_coord(position(diagram, wire.target)))
    end
  end
  y = Float64[ isempty(y) ? 0.0 : mean(y) for y in ys ]
  in_coords = solve_isotonic(y; lower=lower, upper=upper, pad=pad)

  # Solve optimization problem for output ports.
  ys = [ Float64[] for i in eachindex(outputs) ]
  for wire in in_wires(diagram, output_id(diagram))
    i, layout = wire.target.port, wire.value::WireLayout
    if !(isnothing(layout) || isempty(layout))
      push!(ys[i], port_coord(position(last(layout))))
    elseif wire.source.box != input_id(diagram)
      push!(ys[i], port_coord(position(diagram, wire.source)))
    end
  end
  y = Float64[ isempty(y) ? 0.0 : mean(y) for y in ys ]
  out_coords = solve_isotonic(y; lower=lower, upper=upper, pad=pad)

  (layout_linear_ports(InputPort, inputs, in_coords, diagram_size, opts),
   layout_linear_ports(OutputPort, outputs, out_coords, diagram_size, opts))
end
solve_isotonic(args...; kw...) = error("Isotonic regression backend not available")

has_port_layout_method(method::Symbol) = has_port_layout_method(Val{method})
has_port_layout_method(::Type{<:Val}) = false
has_port_layout_method(::Type{Val{:fixed}}) = true

layout_port(value; kw...) = PortLayout(; value=value, kw...)

function position(diagram::WiringDiagram, port::Port)
  pos = position(port_value(diagram, port))
  if !(port.box in outer_ids(diagram))
    pos += position(box(diagram, port.box))
  end
  pos
end

function normal(diagram::WiringDiagram, port::Port)
  (port.box in outer_ids(diagram) ? -1 : +1) * normal(port_value(diagram, port))
end

# Wire layout
#############

""" Lay out the wires in a wiring diagram with no boxes.
"""
function layout_pure_wiring(diagram::WiringDiagram, opts::LayoutOptions;
                            anchor_wires::Bool=true)
  @assert nboxes(diagram) == 0
  inputs, outputs = input_ports(diagram), output_ports(diagram)
  size = minimum_diagram_size(length(inputs), length(outputs), opts)
  result = WiringDiagram(BoxLayout(size=size), inputs, outputs)
  if anchor_wires
    result = layout_outer_ports!(result, opts, method=:fixed)
    e1, e2 = svector(opts, 1, 0), svector(opts, 0, 1)
    for wire in wires(diagram)
      src, tgt = wire.source, wire.target
      src_pos, tgt_pos = (position(port_value(result, p)) for p in (src, tgt))
      mid_pos = (src_pos + tgt_pos) / 2
      v = sign.((tgt_pos - src_pos) .* e1) + (tgt.port - src.port) .* e2
      point = WirePoint(mid_pos .* e2, normalize(v))
      add_wire!(result, Wire([point], src, tgt))
    end
  else
    add_wires!(result, wires(diagram))
  end
  result
end

merge_wire_layouts(left::WireLayout, middle::WireLayout, right::WireLayout) =
  vcat(wire_points(left), wire_points(middle), wire_points(right))

# Labels
########

""" Label for box in wiring diagram.
"""
box_label(mime::MIME, value) = diagram_element_label(mime, value)

""" Label for wire in wiring diagram.

Note: This function takes a port value, not a wire value.
"""
wire_label(mime::MIME, value) = diagram_element_label(mime, value)

""" Label for element in wiring diagram.

By default, both `box_label` and `wire_label` fall back to this function.
"""
diagram_element_label(::MIME, value) = string(value)
diagram_element_label(::MIME, ::Nothing) = ""

function diagram_element_label(::MIME"text/latex", expr::GATExpr)
  string('$', sprint(show_latex, expr), '$')
end
function diagram_element_label(::MIME"text/latex", s::Union{String,Symbol})
  string('$', s, '$')
end

end


text:
module Test

  using Base.Test;
  using Promisables;

  @testset "await macro" begin
    @testset "must wait for a promise to resolve" begin
      p0 = Promise();
      p1 = x -> x + 2;
      p2 = x -> Resolve(x * 2);
      p3 = x -> begin
        x ^ 2;
        return x^2;
      end

      final = foldl((prev, next) -> Then(next, prev), p0, [p1, p2, p3]);

      Timer((t) -> begin
        Fulfill(p0, 0)
        close(t)
      end, 2);

      b = @pawait final;
      @test b == 1;
      @test final.value == 16;
    end
  end

  @testset "A promise" begin
    @testset "when already resolved and given to a then" begin
      @testset "must call with the given function" begin
        p1 = Promise();
        Fulfill(p1, 2);
        handler = (value) -> value + 2;
        p2 = Then(handler, p1);
        @pawait p2;
        @test p2.value === 4;
      end
    end
    @testset "when already rejected and given to a then" begin
      @testset "must call with the given error handler function" begin
        p1 = Promise();
        Reject(p1, ErrorException("foo"));
        handler = (value) -> error("failed"); ## This should not be called
        errHandler = (err) -> Resolve(2);
        p2 = Then(handler, errHandler, p1);
        @pawait p2;
        @test p2.value === 2;
      end
    end
    @testset "must resolve when fullfilled" begin
      chan = Channel{Any}(32)
      p = Promise(chan);
      Fulfill(p, "foo");
      result = take!(chan);
      @test result == "foo";
      @test typeof(p.status) == Fulfilled;
      @test p.value == "foo";
    end

    @testset "must reject when rejected" begin
      chan = Channel{Any}(32)
      p = Promise(chan);
      Reject(p, ErrorException("Basic Error"));
      result = take!(chan);
      @test typeof(result) == ErrorException;
      @test typeof(p.status) == Rejected;
      @test p.value == nothing;
    end

    @testset "when resolved with another promise" begin
      @testset "should chain promises" begin
        p1 = Promise();
        p2 = Promise();
        Fulfill(p1, p2);
        @async Fulfill(p2, "a successful value");
        result = take!(p1.channel);
        @test p1.value ==  p2.value;
      end
    end

    @testset "after then is called" begin
      @testset "when fullfilled" begin
        @testset "must call the then's next function" begin
          chan = Channel{Any}(32);
          p = Promise();
          p1 = Then((k) -> put!(chan, k), p);
          Fulfill(p, "foo");
          result = take!(chan);
          @test result == "foo";
          @test typeof(p.status) == Fulfilled;
          @test typeof(p1.status) == Fulfilled;
          @test p.value == "foo";
        end

      end

      @testset "when rejected" begin
        @testset "must call the error handler if one exists with the given error" begin
          chan = Channel{Any}(32);
          tE = () -> error("failed");
          errorH = (err) -> put!(chan, err);
          p = Promise();
          p1 = Then(tE, errorH, p);
          givenError = ErrorException("Blank");
          Reject(p, givenError);
          result = take!(chan);
          @test typeof(p.status) == Rejected;
          @test typeof(p1.status) == Rejected;
          @test result == givenError;
        end
        @testset "if no error handler is present" begin
          @testset "it should not continue" begin
            chan = Channel{Any}(32);
            tE = () -> error("failed");
            p = Promise();
            p1 = Then(tE, p);
            givenError = ErrorException("Blank");
            Reject(p, givenError);
            result = take!(p1.channel);
            @test typeof(p.status) == Rejected;
            @test typeof(p1.status) == Rejected;
            @test result == givenError;
          end
        end
        @testset "error handler" begin
          @testset "when returning a promise" begin
            @testset "should continue with the promise chain" begin
              chan = Channel{Any}(32);
              tE = () -> error("failed"); # Should throw if called
              p = Promise();
              handler = (err) -> Resolve("solved");
              p1 = Then(tE, handler, p);
              givenError = ErrorException("Blank");
              Reject(p, givenError);
              result = take!(p1.channel);
              @test typeof(p.status) == Rejected;
              @test typeof(p1.status) == Fulfilled;
              @test result == "solved";
            end
          end
          @testset "when returning anything else" begin
            @testset "should not continue with the promise chain" begin
              chan = Channel{Any}(32);
              tE = () -> error("failed"); # Should throw if called
              p = Promise();
              handler = (err) -> Reject(ErrorException("Blank"));
              p1 = Then(tE, handler, p);
              givenError = ErrorException("Blank");
              Reject(p, givenError);
              result = take!(p1.channel);
              @test typeof(p.status) == Rejected;
              @test typeof(p1.status) == Rejected;
            end
          end
        end
      end
    end
  end
end


text:
function build_data_dictionary()

    # Initialize -
    data_dictionary = Dict{String,Any}()

    # setup system dimension -

    # setup initial conditions -
    initial_condition_array = [

    ]

    # --- DO NOT EDIT BELOW THIS LINE ----------------------------------------------- #
    data_dictionary["initial_condition_array"] = initial_condition_array

    return data_dictionary
    # --- DO NOT EDIT ABOVE THIS LINE ----------------------------------------------- #
end

max_stars_repo_name: varnerlab/N-Compartment-Chip-Model

max_issues_repo_name: varnerlab/N-Compartment-Chip-Model

max_forks_repo_name: varnerlab/N-Compartment-Chip-Model


text:
"""
    Zeros.@pirate Base

Extend some functions in Base to return `Zero()` or `One()`, even when called without
any arguments of types "owned" by Zeros.jl. (This is sometimes referred to as "type piracy".)

The function calls `+()` and `*()` (without arguments) will return `Zero()` and `One()`
respectively. The same will be true for `zero()` and `one()` without arguments.

The functions `zero(Any)` and `one(Any)` will also return `Zero()` and `One()`, respectively,
which in turn will make `sum([])` and `prod([])` return those values.

Finally, `sin(π)` and `tan(π)` will return `Zero()`, fixing a bug in Base.
"""
macro pirate(m::Symbol)
    esc(pirate_code(Val(m)))
end

pirate_code(::Val{:Base}) = quote
   Base.:+() = Zero()
   Base.:*() = One()
   Base.sum(::Tuple{}) = Zero()
   Base.zero() = Zero()
   Base.one() = One()
   Base.zero(::Type{Any}) = Zero()
   Base.one(::Type{Any}) = One()
   Base.sin(::Irrational{:π}) = Zero()
   Base.tan(::Irrational{:π}) = Zero()
end


text: """WIP Work in progress""

max_stars_repo_name: adamwillisXanax/accountingJulia

max_issues_repo_name: adamwillisXanax/accountingJulia

max_forks_repo_name: adamwillisXanax/accountingJulia


text:
fname = joinpath(dirname(@__FILE__), "..", "data", "test-4d.dat")
sname = joinpath(dirname(@__FILE__), "..", "data", "tmp", "same.dat")

#x, y, z, s, t = read_dat(fname)

#@test size(x) == (30,)
#@test size(y) == (36,)
#@test size(z) == (28,)
#@test size(s) == (30,36,28,2)
#@test size(t) == (2,)

#@test maximum(x) == 72.5
#@test maximum(y) == 71.220001
#@test maximum(z) == 76.809998
#@test maximum(s) == 0.067409396
#@test maximum(t) == 0.24
#
#@test minimum(x) == -72.5
#@test minimum(y) == -103.779999
#@test minimum(z) == -58.189999
#@test minimum(s) == 0.0
#@test minimum(t) == 0.12

# write_dat(sname, x, y, z, s[:,:,:,:], t)

# x2, y2, z2, s2, t2 = read_dat(sname)

# @test x == x2
# @test y == y2
# @test z == z2
# @test s == s2
# @test t == t2

# fname = joinpath(dirname(@__FILE__), "..", "data", "test-3d.dat")
# x, y, z, s, t = read_dat(fname)

#@test size(x) == (30,)
#@test size(y) == (36,)
#@test size(z) == (28,)
#@test size(s) == (30,36,28,1)
#@test size(t) == (1,)

#@test maximum(x) == 72.5
#@test maximum(y) == 71.220001
#@test maximum(z) == 76.809998
#@test maximum(s) == 33.2692985535
#@test maximum(t) == 0

#@test minimum(x) == -72.5
#@test minimum(y) == -103.779999
#@test minimum(z) == -58.189999
#@test minimum(s) == -7.5189352036
#@test minimum(t) == 0


text:
@test NewModule(:Main, nothing, Dict()) != nothing
@test NewTypeName(:TN, MainModule) != nothing
@test NewTypeVar(:X, NewBottom, NewAny) != nothing


text:
# middle

@test middle(3) === 3.0
@test middle(2, 3) === 2.5
@test middle(1:8) === 4.5
@test middle([1:8]) === 4.5

# ensure type-correctness
for T in [Bool,Int8,Int16,Int32,Int64,Int128,Uint8,Uint16,Uint32,Uint64,Uint128,Float16,Float32,Float64]
    @test middle(one(T)) === middle(one(T), one(T))
end

# median
@test median([1.]) === 1.
@test median([1.,3]) === 2.
@test median([1.,3,2]) === 2.

@test median([1,3,2]) === 2.0
@test median([1,3,2,4]) === 2.5

@test median([0.0,Inf]) == Inf
@test median([0.0,-Inf]) == -Inf
@test median([0.,Inf,-Inf]) == 0.0
@test median([1.,-1.,Inf,-Inf]) == 0.0
@test isnan(median([-Inf,Inf]))

@test all(median([2 3 1 -1; 7 4 5 -4], 2) .== [1.5, 4.5])
@test all(median([2 3 1 -1; 7 4 5 -4], 1) .== [4.5 3.5 3.0 -2.5])

@test_throws ErrorException median([])
@test isnan(median([NaN]))
@test isnan(median([0.0,NaN]))
@test isnan(median([NaN,0.0]))
@test isequal(median([NaN 0.0; 1.2 4.5], 2), reshape([NaN; 2.85], 2, 1))

@test mean([1,2,3]) == 2.
@test mean([0 1 2; 4 5 6], 1) == [2.  3.  4.]
@test mean([1 2 3; 4 5 6], 1) == [2.5 3.5 4.5]

# test var & std

# edge case: empty vector
# iterable; this has to throw for type stability
@test_throws ErrorException var(())
@test_throws ErrorException var((); corrected=false)
@test_throws ErrorException var((); mean=2)
@test_throws ErrorException var((); mean=2, corrected=false)
# reduction
@test isnan(var(Int[]))
@test isnan(var(Int[]; corrected=false))
@test isnan(var(Int[]; mean=2))
@test isnan(var(Int[]; mean=2, corrected=false))
# reduction across dimensions
@test_approx_eq var(Int[], 1) [NaN]
@test_approx_eq var(Int[], 1; corrected=false) [NaN]
@test_approx_eq var(Int[], 1; mean=[2]) [NaN]
@test_approx_eq var(Int[], 1; mean=[2], corrected=false) [NaN]

# edge case: one-element vector
# iterable
@test isnan(@inferred(var((1,))))
@test var((1,); corrected=false) === 0.0
@test var((1,); mean=2) === Inf
@test var((1,); mean=2, corrected=false) === 1.0
# reduction
@test isnan(@inferred(var([1])))
@test var([1]; corrected=false) === 0.0
@test var([1]; mean=2) === Inf
@test var([1]; mean=2, corrected=false) === 1.0
# reduction across dimensions
@test_approx_eq @inferred(var([1], 1)) [NaN]
@test_approx_eq var([1], 1; corrected=false) [0.0]
@test_approx_eq var([1], 1; mean=[2]) [Inf]
@test_approx_eq var([1], 1; mean=[2], corrected=false) [1.0]

@test var(1:8) == 6.

@test_approx_eq varm([1,2,3], 2) 1.
@test_approx_eq var([1,2,3]) 1.
@test_approx_eq var([1,2,3]; corrected=false) 2.0/3
@test_approx_eq var([1,2,3]; mean=0) 7.
@test_approx_eq var([1,2,3]; mean=0, corrected=false) 14.0/3

@test_approx_eq varm((1,2,3), 2) 1.
@test_approx_eq var((1,2,3)) 1.
@test_approx_eq var((1,2,3); corrected=false) 2.0/3
@test_approx_eq var((1,2,3); mean=0) 7.
@test_approx_eq var((1,2,3); mean=0, corrected=false) 14.0/3

@test_approx_eq var([1 2 3 4 5; 6 7 8 9 10], 2) [2.5 2.5]'
@test_approx_eq var([1 2 3 4 5; 6 7 8 9 10], 2; corrected=false) [2.0 2.0]'

@test_approx_eq stdm([1,2,3], 2) 1.
@test_approx_eq std([1,2,3]) 1.
@test_approx_eq std([1,2,3]; corrected=false) sqrt(2.0/3)
@test_approx_eq std([1,2,3]; mean=0) sqrt(7.0)
@test_approx_eq std([1,2,3]; mean=0, corrected=false) sqrt(14.0/3)

@test_approx_eq stdm((1,2,3), 2) 1.
@test_approx_eq std((1,2,3)) 1.
@test_approx_eq std((1,2,3); corrected=false) sqrt(2.0/3)
@test_approx_eq std((1,2,3); mean=0) sqrt(7.0)
@test_approx_eq std((1,2,3); mean=0, corrected=false) sqrt(14.0/3)

@test_approx_eq std([1 2 3 4 5; 6 7 8 9 10], 2) sqrt([2.5 2.5]')
@test_approx_eq std([1 2 3 4 5; 6 7 8 9 10], 2; corrected=false) sqrt([2.0 2.0]')

A = Complex128[exp(i*im) for i in 1:10^4]
@test_approx_eq varm(A,0.) sum(map(abs2,A))/(length(A)-1)
@test_approx_eq varm(A,mean(A)) var(A)

# test covariance

function safe_cov(x, y, zm::Bool, cr::Bool)
    n = length(x)
    if !zm
        x = x .- mean(x)
        y = y .- mean(y)
    end
    dot(vec(x), vec(y)) / (n - int(cr))
end

X = [1. 2. 3. 4. 5.; 5. 4. 6. 2. 1.]'
Y = [6. 1. 5. 3. 2.; 2. 7. 8. 4. 3.]'

for vd in [1, 2], zm in [true, false], cr in [true, false]
    # println("vd = $vd: zm = $zm, cr = $cr")
    if vd == 1
        k = size(X, 2)
        Cxx = zeros(k, k)
        Cxy = zeros(k, k)
        for i = 1:k, j = 1:k
            Cxx[i,j] = safe_cov(X[:,i], X[:,j], zm, cr)
            Cxy[i,j] = safe_cov(X[:,i], Y[:,j], zm, cr)
        end
        x1 = vec(X[:,1])
        y1 = vec(Y[:,1])
    else
        k = size(X, 1)
        Cxx = zeros(k, k)
        Cxy = zeros(k, k)
        for i = 1:k, j = 1:k
            Cxx[i,j] = safe_cov(X[i,:], X[j,:], zm, cr)
            Cxy[i,j] = safe_cov(X[i,:], Y[j,:], zm, cr)
        end
        x1 = vec(X[1,:])
        y1 = vec(Y[1,:])
    end

    c = zm ? cov(x1; mean=0, corrected=cr) :
             cov(x1; corrected=cr)
    @test isa(c, Float64)
    @test_approx_eq c Cxx[1,1]

    C = zm ? cov(X; vardim=vd, mean=0, corrected=cr) :
             cov(X; vardim=vd, corrected=cr)
    @test size(C) == (k, k)
    @test_approx_eq C Cxx

    c = zm ? cov(x1, y1; mean=0, corrected=cr) :
             cov(x1, y1; corrected=cr)
    @test isa(c, Float64)
    @test_approx_eq c Cxy[1,1]

    C = zm ? cov(x1, Y; vardim=vd, mean=0, corrected=cr) :
             cov(x1, Y; vardim=vd, corrected=cr)
    @test size(C) == (1, k)
    @test_approx_eq C Cxy[1,:]

    C = zm ? cov(X, y1; vardim=vd, mean=0, corrected=cr) :
             cov(X, y1; vardim=vd, corrected=cr)
    @test size(C) == (k, 1)
    @test_approx_eq C Cxy[:,1]

    C = zm ? cov(X, Y; vardim=vd, mean=0, corrected=cr) :
             cov(X, Y; vardim=vd, corrected=cr)
    @test size(C) == (k, k)
    @test_approx_eq C Cxy
end

# test correlation

function safe_cor(x, y, zm::Bool)
    if !zm
        x = x .- mean(x)
        y = y .- mean(y)
    end
    x = vec(x)
    y = vec(y)
    dot(x, y) / (sqrt(dot(x, x)) * sqrt(dot(y, y)))
end

for vd in [1, 2], zm in [true, false]
    # println("vd = $vd: zm = $zm")
    if vd == 1
        k = size(X, 2)
        Cxx = zeros(k, k)
        Cxy = zeros(k, k)
        for i = 1:k, j = 1:k
            Cxx[i,j] = safe_cor(X[:,i], X[:,j], zm)
            Cxy[i,j] = safe_cor(X[:,i], Y[:,j], zm)
        end
        x1 = vec(X[:,1])
        y1 = vec(Y[:,1])
    else
        k = size(X, 1)
        Cxx = zeros(k, k)
        Cxy = zeros(k, k)
        for i = 1:k, j = 1:k
            Cxx[i,j] = safe_cor(X[i,:], X[j,:], zm)
            Cxy[i,j] = safe_cor(X[i,:], Y[j,:], zm)
        end
        x1 = vec(X[1,:])
        y1 = vec(Y[1,:])
    end

    c = zm ? cor(x1; mean=0) : cor(x1)
    @test isa(c, Float64)
    @test_approx_eq c Cxx[1,1]

    C = zm ? cor(X; vardim=vd, mean=0) : cor(X; vardim=vd)
    @test size(C) == (k, k)
    @test_approx_eq C Cxx

    c = zm ? cor(x1, y1; mean=0) : cor(x1, y1)
    @test isa(c, Float64)
    @test_approx_eq c Cxy[1,1]

    C = zm ? cor(x1, Y; vardim=vd, mean=0) : cor(x1, Y; vardim=vd)
    @test size(C) == (1, k)
    @test_approx_eq C Cxy[1,:]

    C = zm ? cor(X, y1; vardim=vd, mean=0) : cor(X, y1; vardim=vd)
    @test size(C) == (k, 1)
    @test_approx_eq C Cxy[:,1]

    C = zm ? cor(X, Y; vardim=vd, mean=0) : cor(X, Y; vardim=vd)
    @test size(C) == (k, k)
    @test_approx_eq C Cxy
end

# test hist

@test sum(hist([1,2,3])[2]) == 3
@test hist([])[2] == []
@test hist([1])[2] == [1]
@test hist([1,2,3],[0,2,4]) == ([0,2,4],[2,1])
@test hist([1,2,3],0:2:4) == (0:2:4,[2,1])
@test all(hist([1:100]/100,0.0:0.01:1.0)[2] .==1)
@test hist([1,1,1,1,1])[2][1] == 5
@test sum(hist2d(rand(100, 2))[3]) == 100
@test hist([1 2 3 4;1 2 3 4]) == (0.0:2.0:4.0, [2 2 0 0; 0 0 2 2])

@test midpoints(1.0:1.0:10.0) == 1.5:1.0:9.5
@test midpoints(1:10) == 1.5:9.5
@test midpoints(Float64[1.0:1.0:10.0]) == Float64[1.5:1.0:9.5]

@test quantile([1,2,3,4],0.5) == 2.5
@test quantile([1., 3],[.25,.5,.75])[2] == median([1., 3])
@test quantile([0.:100.],[.1,.2,.3,.4,.5,.6,.7,.8,.9])[1] == 10.0

max_stars_repo_name: Aerlinger/julia

max_issues_repo_name: Aerlinger/julia

max_forks_repo_name: Aerlinger/julia


text:
# This file contains functions for the Revenue Efficiency DDF DEA model
"""
    RevenueDDFDEAModel
An data structure representing a revenue DDF DEA model.
"""
struct RevenueDDFDEAModel <: AbstractRevenueDEAModel
    n::Int64
    m::Int64
    s::Int64
    Gy::Symbol
    monetary::Bool
    dmunames::Union{Vector{String},Nothing}
    rts::Symbol
    eff::Vector
    lambda::SparseMatrixCSC{Float64, Int64}
    techeff::Vector
    alloceff::Vector
    normalization::Vector
    Xtarget::Matrix
    Ytarget::Matrix
end

"""
    dearevenueddf(X, Y, P; Gy)
Compute revenue efficiency using directional distance function data envelopment analysis for
inputs `X`, outputs `Y` and price of outputs `P`.

# Direction specification:

The direction `Gy` can be one of the following symbols.
- `:Ones`: use ones.
- `:Observed`: use observed values.
- `:Mean`: use column means.
- `:Monetary`: use direction so that profit inefficiency is expressed in monetary values.

Alternatively, a vector or matrix with the desired directions can be supplied.

# Optional Arguments
- `rts=:VRS`: chooses variable returns to scale. For constant returns to scale choose `:CRS`.

# Examples
```jldoctest
julia> X = [1; 1; 1; 1; 1; 1; 1; 1];

julia> Y = [7 7; 4 8; 8 4; 3 5; 3 3; 8 2; 6 4; 1.5 5];

julia> P = [1 1; 1 1; 1 1; 1 1; 1 1; 1 1; 1 1; 1 1];

julia> dearevenueddf(X, Y, P, Gy = :Monetary)
Revenue DDF DEA Model
DMUs = 8; Inputs = 1; Outputs = 2
Orientation = Output; Returns to Scale = VRS
Gy = Monetary
─────────────────────────────────
   Revenue  Technical  Allocative
─────────────────────────────────
1      0.0       0.0         0.0
2      2.0       0.0         2.0
3      2.0       0.0         2.0
4      6.0       5.0         1.0
5      8.0       8.0         0.0
6      4.0       0.0         4.0
7      4.0       3.0         1.0
8      7.5       5.75        1.75
─────────────────────────────────
```
"""
function dearevenueddf(X::Union{Matrix,Vector}, Y::Union{Matrix,Vector},
    P::Union{Matrix,Vector};
    Gy::Union{Symbol,Matrix,Vector},
    rts::Symbol = :VRS, monetary::Bool = false,
    names::Union{Vector{String},Nothing} = nothing,
    optimizer::Union{DEAOptimizer,Nothing} = nothing)::RevenueDDFDEAModel

    # Check parameters
    nx, m = size(X, 1), size(X, 2)
    ny, s = size(Y, 1), size(Y, 2)

    np, sp = size(P, 1), size(P, 2)

    if nx != ny
        throw(DimensionMismatch("number of rows in X and Y ($nx, $ny) are not equal"));
    end
    if np != ny
        throw(DimensionMismatch("number of rows in P and Y ($np, $ny) are not equal"));
    end
    if sp != s
        throw(DimensionMismatch("number of columns in P and Y ($sp, $s) are not equal"));
    end

    # Build or get user directions
    if typeof(Gy) == Symbol
        Gysym = Gy

        if Gy == :Ones
            Gy = ones(size(Y))
        elseif Gy == :Observed
            Gy = Y
        elseif Gy == :Mean
            Gy = repeat(mean(Y, dims = 1), size(Y, 1))
        elseif Gy == :Monetary
            GxGydollar = 1 ./ (sum(P, dims = 2));
            Gy = repeat(GxGydollar, 1, s);
        else
            throw(ArgumentError("Invalid `Gy`"));
        end

    else
        Gysym = :Custom
    end

    if (size(Gy, 1) != size(Y, 1)) | (size(Gy, 2) != size(Y, 2))
        throw(DimensionMismatch("size of Gy and Y ($(size(Gy)), $(size(Y))) are not equal"));
    end

    # Set input direction to zeros
    Gx = zeros(size(X))

    # Default optimizer
    if optimizer === nothing
        optimizer = DEAOptimizer(:LP)
    end

    # Get maximum revenue targets and lambdas
    n = nx

    Xtarget = X[:,:]
    Ytarget, rlambdaeff = deamaxrevenue(X, Y, P, rts = rts, dispos = :Strong, optimizer = optimizer)

    # Revenue, technical and allocative efficiency
    refficiency  = vec(sum(P .* Ytarget, dims = 2)  .- sum(P .* Y, dims = 2))
    normalization = vec(sum(P .* Gy, dims = 2))
    techefficiency = efficiency(deaddf(X, Y, Gx = Gx, Gy = Gy, rts = rts, slack = false, optimizer = optimizer))

    if monetary
        techefficiency = techefficiency .* normalization
    else
        refficiency = refficiency ./ normalization
    end

    allocefficiency = refficiency .- techefficiency

    return RevenueDDFDEAModel(n, m, s, Gysym, monetary, names, rts, refficiency, rlambdaeff, techefficiency, allocefficiency, normalization, Xtarget, Ytarget)
end

function Base.show(io::IO, x::RevenueDDFDEAModel)
    compact = get(io, :compact, false)

    n = nobs(x)
    m = ninputs(x)
    s = noutputs(x)
    dmunames = names(x)

    eff = efficiency(x)
    techeff = efficiency(x, :Technical)
    alloceff = efficiency(x, :Allocative)

    if !compact
        print(io, "Revenue DDF DEA Model \n")
        print(io, "DMUs = ", n)
        print(io, "; Inputs = ", m)
        print(io, "; Outputs = ", s)
        print(io, "\n")
        print(io, "Orientation = Output")
        print(io, "; Returns to Scale = ", string(x.rts))
        print(io, "\n")
        print(io, "Gy = ", string(x.Gy))
        print(io, "\n")
        show(io, CoefTable(hcat(eff, techeff, alloceff), ["Revenue", "Technical", "Allocative"], dmunames))
    end

end


text:
export q, @q

const q = _Q()

function (f::_Q)(cmd::String, args...)
    x = k(KDB_HANDLE[], cmd, map(K_new, args)...)
    systemerror("k", x == C_NULL)
    K(x)
end

macro q(ex::Expr)
    if ex.head === :let
        f = :(Dict())
        for s in ex.args[2:end]
            push!(f.args, :($(Meta.quot(s)) => $(esc(s))))
        end
        :(q("eval", resolve!($(esc(K(ex.args[1]))),$f)))
    else
        :(q("eval", $(esc(K(ex)))))
    end
end

macro q(x::Symbol)
    cmd = string(x)
    if cmd[1] == '_'
        cmd = "." * cmd[2:end]
    end
    :(q($cmd))
end

raw"""
 An idea for implementing the following syntax:

       @q let x, y
          x + til(y)
       end

The invocation above should treat x and y as Julia variables that
need to be sent to kdb+.

julia> macro frame(ex)
          f = :(Dict())
          for s in ex.args[2:end]
              push!(f.args, :($(Meta.quot(s)) => $(esc(s))))
          end
          f
       end
@frame (macro with 1 method)

julia> @macroexpand @frame(let x, y;end)
:((Main.Dict)(:x => x, :y => y))

The @q(let x,y;ex end) should expand into q("eval", resolve!(K(ex.args[1]), f))
expression, where f is the frame similar to the one generated by the @frame
macro above and resolve() is a function that substitutes the values from f to
a K parse tree.
"""
function resolve!(e, f)
    # TODO: Check e[1] and treat : differently.
    for (i, x) in enumerate(e)
        t = ktypecode(x)
        if t == KK
            resolve!(x, f)
        elseif t == -KS
            sym = Symbol(x)
            if haskey(f, sym)
                e[i] = K(f[sym])
            end
        end
    end
    e
end