text:
@testset "Matmul" begin

    # Test multiplocations with Double64 and Float64
    A1 = [1.0 2 3; 4 5 6; 7 8 9]
    B1 = [7.0 8; 9 10; 11 12]
    A1B1 = A1*B1

    DA1DB1 = Double64.(A1)*Double64.(B1)
    DA1FB1 = Double64.(A1)*B1
    FA1DB1 = A1*Double64.(B1)

    @test DA1DB1 == Double64.(A1B1)
    @test DA1FB1 == Double64.(A1B1)
    @test FA1DB1 == Double64.(A1B1)
    @test DA1DB1 isa Array{Double64,2}
    @test DA1FB1 isa Array{Double64,2}
    @test FA1DB1 isa Array{Double64,2}

    # Also test transpose
    @test Double64.(B1)'Double64.(A1)' == Double64.(A1B1')
    # Incompatible Dimension
    @test_throws DimensionMismatch Double64.(B1)*Double64.(A1)

    #Test with zero inner dimension (zero result)
    A2 = Array{Float64,2}(undef, 2, 0)
    B2 = Array{Float64,2}(undef, 0, 3)
    A2B2 = A2*B2

    DA2DB2 = Double64.(A2)*Double64.(B2)
    DA2FB2 = Double64.(A2)*B2
    FA2DB2 = A2*Double64.(B2)

    @test DA2DB2 == Double64.(A2B2)
    @test DA2FB2 == Double64.(A2B2)
    @test FA2DB2 == Double64.(A2B2)
    @test DA2DB2 isa Array{Double64,2}
    @test DA2FB2 isa Array{Double64,2}
    @test FA2DB2 isa Array{Double64,2}

    #Test with zero outer dimensions (empty result)
    A3 = Array{Float64,2}(undef, 0, 3)
    B3 = Array{Float64,2}(undef, 3, 0)
    A3B3 = A3*B3

    DA3DB3 = Double64.(A3)*Double64.(B3)
    DA3FB3 = Double64.(A3)*B3
    FA3DB3 = A3*Double64.(B3)

    @test DA3DB3 == Double64.(A3B3)
    @test DA3FB3 == Double64.(A3B3)
    @test FA3DB3 == Double64.(A3B3)
    @test DA3DB3 isa Array{Double64,2}
    @test DA3FB3 isa Array{Double64,2}
    @test FA3DB3 isa Array{Double64,2}

end


text:
"""
Given a dataset, generate a closure that computes the negative log likelihood
for a model and takes only a vector of parameters as an argument.
"""
function simulate(
    m::Model,
    inputs::Vector,
    n_trials::Integer,
)
    sim_inputs = similar(inputs, 0)
    sim_weights = Int[]
    n = length(inputs)

    for i in 1:n
        @inbounds x2, t2, x1, t1, y = inputs[i]
        p = predict(m, x2, t2, x1, t1)
        w = rand(Distributions.Binomial(n_trials, p))
        push!(sim_inputs, (x2, t2, x1, t1, 1.0))
        push!(sim_weights, w)
        push!(sim_inputs, (x2, t2, x1, t1, 0.0))
        push!(sim_weights, n_trials - w)
    end

    return sim_inputs, sim_weights
end

ext: jl


text:
# AUTO GENERATED FILE - DO NOT EDIT

export ''_leafletmaplistener

"""
    ''_leafletmaplistener(;kwargs...)

A LeafletMapListener component.

Keyword arguments:
- `id` (String; optional)
- `_bounds` (optional): . _bounds has the following type: lists containing elements 'minx', 'miny', 'maxx', 'maxy'.
Those elements have the following types:
  - `minx` (Real; optional)
  - `miny` (Real; optional)
  - `maxx` (Real; optional)
  - `maxy` (Real; optional)
- `_center` (optional): . _center has the following type: lists containing elements 'lng', 'lat'.
Those elements have the following types:
  - `lng` (Real; optional)
  - `lat` (Real; optional)
- `_clickedLatLng` (optional): . _clickedLatLng has the following type: lists containing elements 'lng', 'lat'.
Those elements have the following types:
  - `lng` (Real; optional)
  - `lat` (Real; optional)
- `_zoom` (Real; optional)
- `loading_state` (optional): . loading_state has the following type: lists containing elements 'is_loading', 'prop_name', 'component_name'.
Those elements have the following types:
  - `is_loading` (Bool; optional): Determines if the component is loading or not
  - `prop_name` (String; optional): Holds which property is loading
  - `component_name` (String; optional): Holds the name of the component that is loading
"""
function ''_leafletmaplistener(; kwargs...)
        available_props = Symbol[:id, :_bounds, :_center, :_clickedLatLng, :_zoom, :loading_state]
        wild_props = Symbol[]
        return Component("''_leafletmaplistener", "LeafletMapListener", "feffery_leaflet_components", available_props, wild_props; kwargs...)
end

max_stars_repo_name: CNFeffery/feffery-leaflet-components

max_issues_repo_name: CNFeffery/feffery-leaflet-components

max_forks_repo_name: CNFeffery/feffery-leaflet-components


text:
export CompileError

# Generate a temporary file with specific suffix
# NOTE: mkstemps is glibc 2.19+, so emulate its behavior
function mkstemps(suffix::AbstractString)
    base = tempname()
    filename = base * suffix

    # make sure the filename is unique
    i = 0
    while isfile(filename)
        i += 1
        filename = base * ".$i" * suffix
    end

    return (filename, Base.open(filename, "w"))
end

macro compile(dev, kernel, code)
    kernel_name = string(kernel)
    containing_file = @__FILE__

    return Expr(:toplevel,
        Expr(:export,esc(kernel)),
        :($(esc(kernel)) = _compile($(esc(dev)), $kernel_name, $code, $containing_file)))
end

immutable CompileError <: Exception
    message::String
end

const builddir = joinpath(@__DIR__, ".cache")

function _compile(dev, kernel, code, containing_file)
    arch = CUDArt.architecture

    if !isdir(builddir)
        println("Writing build artifacts to $builddir")
        mkpath(builddir)
    end

    # check if we need to compile
    codehash = hex(hash(code))
    output = "$builddir/$(kernel)_$(codehash)-$(arch).ptx"
    if isfile(output)
        need_compile = (stat(containing_file).mtime > stat(output).mtime)
    else
        need_compile = true
    end

    # compile the source, if necessary
    if need_compile
        # write the source to a compilable file
        (source, io) = mkstemps(".cu")
        write(io, """
extern "C"
{
$code
}
""")
        Base.close(io)

        compiler = CUDArt.cuda_compiler
        flags = ["--compiler-bindir", CUDArt.host_compiler, "--gpu-architecture", arch]

        err = Pipe()
        cmd = `$compiler $flags -ptx -o $output $source`
        result = success(pipeline(cmd; stdout=DevNull, stderr=err))
        Base.close(err.in)
        rm(source)

        errors = readstring(err)
        if !result
            throw(CompileError("compilation of kernel $kernel failed\n$errors"))
        elseif !isempty(errors)
            warn("during compilation of kernel $kernel:\n$errors")
        end

        if !isfile(output)
            error("compilation of kernel $kernel failed (no output generated)")
        end
    end

    mod = CUDAdrv.CuModuleFile(output)
    return CUDAdrv.CuFunction(mod, kernel)
end

function clean_cache()
    if ispath(builddir)
        @assert isdir(builddir)
        rm(builddir; recursive=true)
    end
end


text:
function isvalid(sim::Chains)
    if !all(gelmandiag(sim).value[:,1,1] .< 1.01)
        return false
    end
    for i in 1:size(sim.value)[2]
        data = mean(sim.value[:,i,:], dims = [1])[:]
        # mu = mean(data)
        sig = std(data)
        if sig > 0.05
            return false
        end
    end
    return true
end

function getstats(value::AbstractArray;standardize::Bool = false)
    # get statistics or natural parameter of given MCMC chains
    # in this function, return the first four moments
    m1 = mean(value, dims = [1,3])[:]
    m2 = var(value, dims = [1,3])[:]
    if standardize
        m3 = map(1:size(value)[2]) do i
            return skewness(value[:,i,:][:])
        end
        m3 = m3[:]
        m4 = map(1:size(value)[2]) do i
            return kurtosis(value[:,i,:][:])
        end
        m4 = m4[:]
    else
        m3 = mean(value.^3, dims = [1,3])[:]
        m4 = mean(value.^4, dims = [1,3])[:]
    end
    res = [m1 m2 m3 m4]
    @assert size(res) == (size(value)[2], 4) "size incorrect"
    return res
end


text:
"""
Estimate if the number is prime
"""
function check_composite(a::BigInt, d::BigInt, n::BigInt, s::BigInt)
    x = exp!(a, d, n)

    if (x == 1 || x == (n - 1))
        return (false)
    end

    for _ in 1:s
        x = (x * x) % n
        if (x == (n - 1))
            return (false)
        end
    end
    true
end

check_composite(a::T, d::T, n::T, s::T) where {T <: Integer} =
check_composite(big(a), big(d), big(n), big(s))

"""Miller Rabin primality test
Return False if n is composite, True(probably prime) otherwise.

n: Number to be tested for primality
k: number of iterations to run the test
"""
function miller_rabin(n::BigInt, k::BigInt = big(10))
    if (n == 2)
        return (true)
    end
    if (n == 1 || n % 2 == 0)
        return (false)
    end

    s::BigInt = 0
    d = n - 1

    while d % 2 == 0
        s += 1
        d /= 2
    end

    @assert floor(BigInt, 2 ^ s * d) == n - 1

    for _ in 1:k
        a = rand(2:n - 1)
        d = floor(BigInt, d)
        if (check_composite(a, d, n, s))
            return (false)
        end
    end
    true
end

miller_rabin(n::T, k::T = 10) where {T <: Integer} =
miller_rabin(big(n), big(k))

"""
Solovay Strassen Primality Test
Returns False is n is composite, True(probably prime) otherwise

n: Number to be tested for primality
k: number of iterations to run the test
"""
function solovay_strassen(n::BigInt, k::BigInt = big(10))
    if (n == 2)
        return (true)
    end
    if (n == 1 || n & 1 == 0)
        return (false)
    end

    for _ in 1:k
        a = rand(2:n-1)
        x = (jacobi!(a, n) + n) % n
        if (x == 0 || exp!(a, floor(BigInt, (n - 1) / 2), n) != x)
            return (false)
        end
    end
    true
end

solovay_strassen(n::T, k::T = 10) where {T <: Integer} =
solovay_strassen(big(n), big(k))

"""
Fermat's Primality test
Returns True(probably prime) if n is prime, Flase if n is composite

n: Number to be tested for primality
k: number of iterations of the Fermat's Primality Test
"""
function fermat_test(n::BigInt, k::BigInt = big(10))
    if (n == 2)
        return (true)
    end
    if (n & 1 == 0)
        return (false)
    end

    for _ in 1:k
        a = rand(2:n-1)
        if (exp!(a, n - 1, n) != 1)
            return (false)
        end
    end
    true
end

fermat_test(n::T, k::T = 10) where {T <: Integer} =
fermat_test(big(n), big(k))


text:
"""
    warmup_from_variability(
        model::MetabolicModel,
        optimizer,
        n_points::Int;
        kwargs...
    )

Generates FVA-like warmup points for samplers, by selecting random points by
minimizing and maximizing reactions. Can not return more than 2 times the
number of reactions in the model.
"""
function warmup_from_variability(model::MetabolicModel, optimizer, n_points::Int; kwargs...)
    nr = n_reactions(model)

    n_points > 2 * nr && throw(
        DomainError(
            n_points,
            "Variability method can not generate more than $(2*nr) points from this model",
        ),
    )

    sample = shuffle(vcat(1:nr, -(1:nr)))[begin:n_points]
    warmup_from_variability(
        model,
        optimizer,
        -filter(x -> x < 0, sample),
        filter(x -> x > 0, sample);
        kwargs...,
    )
end

"""
    function warmup_from_variability(
        model::MetabolicModel,
        optimizer,
        min_reactions::Vector{Int}=1:n_reactions(model),
        max_reactions::Vector{Int}=1:n_reactions(model);
        modifications = [],
        workers::Vector{Int} = [myid()],
    )::Tuple{Matrix{Float64}, Vector{Float64}, Vector{Float64}}

Generate FVA-like warmup points for samplers, by minimizing and maximizing the
specified reactions. The result is returned as a matrix, each point occupies as
single column in the result.
"""
function warmup_from_variability(
    model::MetabolicModel,
    optimizer,
    min_reactions::AbstractVector{Int} = 1:n_reactions(model),
    max_reactions::AbstractVector{Int} = 1:n_reactions(model);
    modifications = [],
    workers::Vector{Int} = [myid()],
)::Tuple{Matrix{Float64},Vector{Float64},Vector{Float64}}

    # create optimization problem at workers, apply modifications
    save_model = :(
        begin
            model = $model
            optmodel = $COBREXA.make_optimization_model(model, $optimizer)
            for mod in $modifications
                mod(model, optmodel)
            end
            optmodel
        end
    )

    map(fetch, save_at.(workers, :cobrexa_sampling_warmup_optmodel, Ref(save_model)))

    fluxes = hcat(
        dpmap(
            rid -> :($COBREXA._FVA_optimize_reaction(
                cobrexa_sampling_warmup_optmodel,
                $rid,
                optmodel -> $COBREXA.JuMP.value.(optmodel[:x]),
            )),
            CachingPool(workers),
            vcat(-min_reactions, max_reactions),
        )...,
    )

    # snatch the bounds from whatever worker is around
    lbs, ubs = get_val_from(
        workers[1],
        :($COBREXA.get_bound_vectors(cobrexa_sampling_warmup_optmodel)),
    )

    # free the data on workers
    map(fetch, remove_from.(workers, :cobrexa_sampling_warmup_optmodel))

    return fluxes, lbs, ubs
end


text:
@testset "Sanity/String" begin
    @test lcs("", "", ==) == ( [ ], [ ] )

    @test lcs("a", "", ==) == ( [ ( true, 1:1:1 ) ], [ ] )
    @test lcs("", "a", ==) == ( [ ], [ ( true, 1:1:1 ) ] )

    @test lcs("ab", "", ==) == ( [ ( true, 1:1:2 ) ], [ ] )
    @test lcs("abc", "", ==) == ( [ ( true, 1:1:3 ) ], [ ] )

    @test lcs("a", "a", ==) == ( [ ( false, 1:1:1 ) ], [ ( false, 1:1:1 ) ] )
    @test lcs("a", "b", ==) == ( [ ( true, 1:1:1 ) ], [ ( true, 1:1:1 ) ] )

    @test lcs("a", "ab", ==) == ( [ ( false, 1:1:1 ) ], [ ( false, 1:1:1 ), ( true, 2:1:2 ) ] )
    @test lcs("a", "bb", ==) == ( [ ( true, 1:1:1 ) ], [ ( true, 1:1:2 ) ] )

    @test lcs("ab", "bb", ==) == ( [ ( true, 1:1:1 ), ( false, 2:1:2 ) ], [ ( true, 1:1:1 ), ( false, 2:1:2 ) ] )

    @test lcs("abc", "bbc", ==) == ( [ ( true, 1:1:1 ), ( false, 2:1:3 ) ], [ ( true, 1:1:1 ), ( false, 2:1:3 ) ] )
    @test lcs("abc", "abb", ==) == ( [ ( false, 1:1:2 ), ( true, 3:1:3 ) ], [ ( false, 1:1:2 ), ( true, 3:1:3 ) ] )
end

max_stars_repo_name: jlegare/Longest-Common-Subsequence

max_issues_repo_name: jlegare/Longest-Common-Subsequence

max_forks_repo_name: jlegare/Longest-Common-Subsequence


text:
# © 2016 Massachusetts Institute of Technology.  See LICENSE file for details.
VALID_ID = Union{Integer, AbstractString, Symbol}
node_name(n::AbstractString, k) = "$(n)_$(k)"

abstract Flow
get_params{T<:Flow}(creek::T) = error("not implemented ($T)")
get_params{T<:Flow}(creeks::AbstractArray{T}) = vcat([get_params(c) for c in creeks]...)
reg_params{T<:Flow}(creek::T) = error("not implemented ($T)")
reg_params{T<:Flow}(creeks::AbstractArray{T}) = vcat([reg_params(c) for c in creeks]...)
#@compat (creek::T){T<:Flow}(args...) = error("not implemented")
Base.size(creek::Flow, args...) = error("not implemented")
Base.eltype(creek::Flow) = error("not implemented")

in_size{T<:Flow}(creek::T) = error("not implemented ($T)")
out_size{T<:Flow}(creek::T) = error("not implemented ($T)")

function scanl{T}(f, ::Type{T}, a0::T, xa)
  y = Array{T}(size(xa))
  a = a0
  for i in 1:length(xa)
    y[i] = a = f(a, xa[i])
  end
  return y
end

function scanl!(f, a0, xa)
  a = a0
  for i in 1:length(xa)
    xa[i] = a = f(a, xa[i])
  end
  return xa
end

# An array of flows behaves as composition
@compat (layers::Array{T}){T<:Flow}(x; kw...) = foldl((x, l)->l(x), x, layers)

# Any unary function (e.g. tanh, sigm) may be used as a flow
immutable Activation{F} <: Flow f::F end
Base.convert{F<:Function}(::Type{Flow}, f::F) = Seep.Activation(f)
Base.convert{F<:Function}(::Type{Activation}, f::F) = Seep.Activation{F}(f)
@compat (a::Activation)(x) = a.f(x)
@compat (a::Activation)(xs::AbstractArray) = map(a.f, xs)

get_params(::Activation) = ANode[]
reg_params(::Activation) = ANode[]

# Uniform Xavier Initialization
function xavier{T}(::Type{T}, shape...; scale=shape[end])
    @assert scale>0
    val = map(T, rand(shape...) / sqrt(scale))
    !haskey(ENV, "SEEP_NO_GPU") ?  CudaArray(val) : val
end
xavier{T}(::Type{T}, nout) = xavier(T, nout, 1)

# Gaussian (Normal) Xavier Initialization
function xaviern{T}(::Type{T}, shape...; scale=shape[end])
    @assert scale>0
    val = map(T, randn(shape...) / sqrt(scale))
    !haskey(ENV, "SEEP_NO_GPU") ?  CudaArray(val) : val
end
xaviern{T}(::Type{T}, nout) = xaviern(T, nout, 1)

function zero_init{T}(::Type{T}, s...)
    val = map(T, zeros(s...))
    !haskey(ENV, "SEEP_NO_GPU") ?  CudaArray(val) : val
end

type Linear <: Flow
    name::VALID_ID
    W::ANode
    b::ANode
end
function Linear{T}(::Type{T}, k, input_size::Integer, output_size::Integer;
    W=xaviern(T, output_size, input_size), b=zero_init(T, output_size))
    Linear(k, ANode(node_name("W", k), W), ANode(node_name("b", k), b))
end
Base.size(l::Linear, args...) = size(l.W, args...)
Base.eltype(l::Linear) = eltype(l.W)
get_params(l::Linear) = ANode[l.W, l.b]
reg_params(l::Linear) = ANode[l.W]
in_size(l::Linear) = size(l.W, 2)
out_size(l::Linear) = size(l.W, 1)
@compat (l::Linear)(x) = l.W * x .+ l.b
@compat (l::Linear)(xs::AbstractArray) = map(l, xs)

type Dropout <: Flow
    a::Real
end
get_params(::Dropout) = ANode[]
reg_params(::Dropout) = ANode[]
@compat (d::Dropout)(x) = (rand!(ANode(size(x))) .> d.a).*x./(1-d.a)
@compat (d::Dropout)(x::AbstractArray) = map(d, x)

type RandomNoise <: Flow
    a::Real
end
get_params(::RandomNoise) = ANode[]
reg_params(::RandomNoise) = ANode[]
@compat (d::RandomNoise)(x) = x.+(randn!(ANode(size(x))))
@compat (d::RandomNoise)(x::AbstractArray) = map(d, x)

type Conv2D <: Flow
    name::VALID_ID
    W::ANode
    b::ANode
    filter_shape
    stride
end
function Conv2D{T}(::Type{T}, k, filter_shape::Tuple{Vararg{Int}}; stride=(),
    W=xaviern(T, filter_shape...; scale=prod(filter_shape[2:end])),
    b=zero_init(T, filter_shape[1]))

    Conv2D(k, ANode(node_name("W", k), W), ANode(node_name("b", k), b), filter_shape, stride)
end
Base.size(l::Conv2D, args...) = size(l.W, args...)
Base.eltype(l::Conv2D) = eltype(l.W)
get_params(l::Conv2D) = ANode[l.W, l.b]
reg_params(l::Conv2D) = ANode[]
@compat (l::Conv2D)(x) = conv2(l.W, x; stride=l.stride) .+ l.b
@compat (l::Conv2D)(xs::AbstractArray) = map(x->l(x), xs)

type Pool <: Flow
    filter_shape
    stride
end
function Pool(filter_shape::Tuple{Vararg{Int}}; stride=filter_shape)
    Pool(filter_shape, stride)
end
get_params(l::Pool) = ANode[]
reg_params(l::Pool) = ANode[]
@compat (l::Pool)(x) = pool(x, l.filter_shape; stride=l.stride)
@compat (l::Pool)(xs::AbstractArray) = map(x->l(x), xs)

type RNN <: Flow
    name::VALID_ID
    Wx::ANode
    Wh::ANode
    b::ANode
end
function RNN{T}(::Type{T}, k, input_size::Integer, output_size::Integer;
    Wx=xaviern(T, output_size, input_size),
    Wh=xaviern(T, output_size, input_size),
    b=zero_init(T, output_size))

    RNN(k, ANode(node_name("Wx", k), Wx), ANode(node_name("Wh", k), Wh), ANode(node_name("b", k), b))
end
Base.size(l::RNN, args...) = size(l.Wx, args...)
Base.eltype(l::RNN) = eltype(l.Wx)
get_params(l::RNN) = ANode[l.Wx, l.Wh, l.b]
reg_params(l::RNN) = ANode[l.Wx, l.Wh]
in_size(l::RNN) = size(l.Wx, 2)
out_size(l::RNN) = size(l.Wx, 1)
@compat (l::RNN)(x, h) = l.Wx*x .+ l.Wh*h .+ l.b
@compat (l::RNN)(xs::AbstractArray, h0=h0_default(l, xs)) = scanl(l, ANode, h0, xs)
h0_default(l::Flow, xs::AbstractArray) = ANode(:zeros, (), (out_size(l), size(xs[1], 2)))

type LSTM <: Flow
    name::VALID_ID
    # Weights
    Wxa::ANode
    Wha::ANode
    ba::ANode
    Wxi::ANode
    Whi::ANode
    bi::ANode
    Wxf::ANode
    Whf::ANode
    bf::ANode
    Wxo::ANode
    Who::ANode
    bo::ANode
end
function LSTM{T}(::Type{T}, k, input_size, output_size;
    Wx=xaviern(T, 4*output_size, input_size), Wh=xaviern(T, 4*output_size, output_size), b=zero_init(T, 4*output_size))

    ind = 1:output_size
    Wxa = ANode(node_name("Wxa", k), Wx[ind, :])
    Wha = ANode(node_name("Wha", k), Wh[ind, :])
    ba = ANode(node_name("ba", k), b[ind])

    ind = ind[end] + (1:output_size)
    Wxi = ANode(node_name("Wxi", k), Wx[ind, :])
    Whi = ANode(node_name("Whi", k), Wh[ind, :])
    bi = ANode(node_name("bi", k), b[ind])

    ind = ind[end] + (1:output_size)
    Wxf = ANode(node_name("Wxf", k), Wx[ind, :])
    Whf = ANode(node_name("Whf", k), Wh[ind, :])
    bf = ANode(node_name("Wbf", k), b[ind])

    ind = ind[end] + (1:output_size)
    Wxo = ANode(node_name("Wxo", k), Wx[ind, :])
    Who = ANode(node_name("Who", k), Wh[ind, :])
    bo = ANode(node_name("bo", k), b[ind])

    LSTM(k, Wxa, Wha, ba, Wxi, Whi, bi, Wxf, Whf, bf, Wxo, Who, bo)
end
Base.size(l::LSTM, args...) = size(l.Wxa, args...)
Base.eltype(l::LSTM) = eltype(l.Wxa)
get_params(l::LSTM) = begin
    params = ANode[]
    push!(params, l.Wxa, l.Wha, l.ba)
    push!(params, l.Wxi, l.Whi, l.bi)
    push!(params, l.Wxf, l.Whf, l.bf)
    push!(params, l.Wxo, l.Who, l.bo)
    params
end
reg_params(l::LSTM) = begin
    params = ANode[]
    push!(params, l.Wxa, l.Wha)
    push!(params, l.Wxi, l.Whi)
    push!(params, l.Wxf, l.Whf)
    push!(params, l.Wxo, l.Who)
    params
end
in_size(l::LSTM) = size(l.Wxa, 2)
out_size(l::LSTM) = size(l.Wxa, 1)

function apply_lstm(l::LSTM, x, c, h)
    a = tanh(l.Wxa * x .+ l.Wha * h .+ l.ba)
    i = sigm(l.Wxi * x .+ l.Whi * h .+ l.bi)
    f = sigm(l.Wxf * x .+ l.Whf * h .+ l.bf)
    o = sigm(l.Wxo * x .+ l.Who * h .+ l.bo)
    c = a .* i + f.* c
    h = o .* tanh(c)
    c, h
end

apply_lstm(l::LSTM, xs::AbstractArray, ch=(h0_default(l, xs), h0_default(l, xs))) =
  scanl((ch,x)->apply_lstm(l, x, ch[1], ch[2]), Tuple{ANode,ANode}, ch, xs)

@compat (l::LSTM)(x) = error("LSTM should be applied to an array of input nodes")
@compat (l::LSTM)(xs::AbstractArray) = map(x->x[2], apply_lstm(l, xs))
@compat (l::LSTM)(xs::AbstractArray, ch) = map(x->x[2], apply_lstm(l, xs, ch))

type GRU <: Flow
    name::VALID_ID
    Wr::ANode
    Ur::ANode
    br::ANode
    Wz::ANode
    Uz::ANode
    bz::ANode
    W::ANode
    U::ANode
    b::ANode
end
function GRU{T}(::Type{T}, k, n, args...)
    Wr = ANode(node_name("Wr", k), xaviern(T, n, n, 2*n))
    Ur = ANode(node_name("Ur", k), xaviern(T, n, n, 2*n))
    br = ANode(node_name("br", k), zero_init(T, n))
    Wz = ANode(node_name("Wz", k), xaviern(T, n, n, 2*n))
    Uz = ANode(node_name("Uz", k), xaviern(T, n, n, 2*n))
    bz = ANode(node_name("bz", k), zero_init(T, n))
    W = ANode(node_name("W", k), xaviern(T, n, n, 2*n))
    U = ANode(node_name("U", k), xaviern(T, n, n, 2*n))
    b = ANode(node_name("b", k), zero_init(T, n))

    GRU(k, Wr, Ur, br, Wz, Uz, bz, W, U, b)
end
Base.size(l::GRU, args...) = size(l.Wxr, args...)
Base.eltype(l::GRU) = eltype(l.Wr)
get_params(gru::GRU) = begin
    params = ANode[]
    push!(params, gru.Wr, gru.Ur, gru.br)
    push!(params, gru.Wz, gru.Uz, gru.bz)
    push!(params, gru.W, gru.U, gru.b)
    params
end
reg_params(gru::GRU) = begin
    params = ANode[]
    push!(params, gru.Wr, gru.Ur)
    push!(params, gru.Wz, gru.Uz)
    push!(params, gru.W, gru.U)
    params
end
in_size(l::GRU) = out_size(l)
out_size(l::GRU) = size(l.Wr, 1)
@compat (g::GRU)(x, h) = begin
    r = sigm(g.Wr * x .+ g.Ur * h .+ g.br)
    z = sigm(g.Wz * x .+ g.Uz * h .+ g.bz)
    hbar = tanh(g.W * x .+ g.U * (r .* h) .+ g.b)
    (1-z) .* g + z .* hbar
end
@compat (g::GRU)(xs::AbstractArray, h=h0_default(g, xs)) = scanl(g, ANode, h, xs)


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule ephemeralpg_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("ephemeralpg")
JLLWrappers.@generate_main_file("ephemeralpg", UUID("1b42a6b1-a0b2-50eb-a0be-acc79e4b08a0"))
end  # module ephemeralpg_jll


text:
const DCHyperParamConfig = NamedTuple{(:λ, :γ),Tuple{Float64,Float64}}

struct DifferenceOfConvex{T <: Algorithm} <: Algorithm
    optimizer_factory
    initialization_algorithm::T
    hyperparameter_set::Array{DCHyperParamConfig}
end

function DifferenceOfConvex(optimizer_factory, initialization_algorithm::T=RandomGuessing()) where {T <: Algorithm}
    return DifferenceOfConvex{T}(
        optimizer_factory,
        initialization_algorithm,
        [(λ = 2.0^λ_exp_val, γ = 10.0^γ_exp_val) for λ_exp_val in -5:5, γ_exp_val in -3:-1]
    )
end

function _get_initializer_learner(existing_learner::Learner{DifferenceOfConvex{T}}) where {T <: Algorithm}
    return Learner(
        existing_learner.train_data,
        existing_learner.test_data,
        existing_learner.validation_data,
        existing_learner.algorithm.initialization_algorithm,
        existing_learner.param_lb,
        existing_learner.param_ub,
        existing_learner.offset_lb,
        existing_learner.offset_ub,
        existing_learner.objective_penalties
    )
end

function δv(r, bid_price, pay_price, γ)
    if r < pay_price
        return -1.0
    elseif r > (1 + γ) * bid_price
        return 1.0
    else
        return 0.0
    end
end

function δV(dataset::DataSet, weights, γ)
    num_features = feature_count(dataset)
    num_samples = sample_count(dataset)
    grads = zeros(num_features)
    for sample in dataset.samples
        linear_value = dot(weights, sample.features)
        for i in 1:num_features
            grads[i] += sample.features[i] * δv(linear_value, sample.bid_price, sample.pay_price, γ)
        end
    end
    return grads
end

function DCA(learner::Learner{DifferenceOfConvex{T}}, w_prev::Vector{Float64}, λ::Float64, γ::Float64) where {T <: Algorithm}
    dataset = learner.train_data
    num_features = feature_count(learner.train_data)
    num_samples = sample_count(learner.train_data)
    jump_model = JuMP.Model(JuMP.with_optimizer(learner.algorithm.optimizer_factory))
    JuMP.@variable(jump_model, w[1:num_features])
    JuMP.@variable(jump_model, s[1:num_samples])
    JuMP.@objective(jump_model, Min, λ * sum(w[i]^2 for i in 1:num_features) + sum(s) - dot(δV(dataset, w_prev, γ), w))
    for i in 1:num_samples
        sample = dataset.samples[i]
        features = sample.features
        bid_price = sample.bid_price
        JuMP.@constraint(jump_model, s[i] ≥ - dot(w, features))
        JuMP.@constraint(jump_model, s[i] ≥ 1 / γ * (dot(w, features) - (1 + γ) * bid_price))
    end
    JuMP.optimize!(jump_model)
    @assert JuMP.primal_status(jump_model) == MOI.FEASIBLE_POINT
    return JuMP.value.(w)
end

function line_search(dataset::DataSet, u::Vector{Float64}, γ::Float64)
    best_r = 0.0
    best_r_reward = Inf
    for _sample in dataset.samples
        _ux_val = dot(u, _sample.features)
        _ux_val <= 0 && continue
        r = _sample.bid_price / _ux_val
        r_reward = 0.0
        for sample in dataset.samples
            ux_val = dot(u, sample.features)
            ux_val <= 0 && continue
            pay_price = sample.pay_price / ux_val
            bid_price = sample.bid_price / ux_val
            if r <= pay_price
                r_reward += -pay_price * ux_val
            elseif pay_price < r <= bid_price
                r_reward += -r * ux_val
            elseif bid_price < r <= (1 + γ) * bid_price
                r_reward += 1 / γ * (r - (1 + γ) * bid_price) * ux_val
            else
                @assert r > (1 + γ) * bid_price
            end
        end
        if r_reward < best_r_reward
            best_r = r
            best_r_reward = r_reward
        end
    end
    return best_r
end

function train(learner::Learner{DifferenceOfConvex{T}}) where {T <: Algorithm}
    num_features = feature_count(learner.train_data)
    num_samples = sample_count(learner.train_data)
    # Subtlety of algorithm: Presume that a fixed λ attains the maximum for
    # any choice of w_{t-1} (i.e. it can stay the same after restart)
    hyperparam_results = TrainingResults[]
    for (λ, γ) in learner.algorithm.hyperparameter_set
        _initialization_learner = _get_initializer_learner(learner)
        _initialization_results = train(_initialization_learner)
        _linear_model = _initialization_results.model
        _audit_weights(_initialization_learner, _linear_model.weights)
        linear_model = _linear_model
        max_iter = 1000
        for iter in 1:max_iter
            iis_init_reward = _evaluate(learner.train_data, linear_model).reward
            v = DCA(learner, linear_model.weights, λ, γ)
            dca_linear_model = LinearModel(v)
            iis_dca_reward = _evaluate(learner.train_data, dca_linear_model).reward
            linear_model = dca_linear_model
            u = v / norm(v)
            η_star = line_search(learner.train_data, u, γ)
            # DIFFERENT THAN IN MM PAPER: Update should be η * u, not η * v
            w = η_star * u
            ls_linear_model = LinearModel(w)
            iis_ls_reward = _evaluate(learner.train_data, ls_linear_model).reward
            # New model is not strictly better than old one
            if !(iis_ls_reward > iis_init_reward)
                break
            end
            linear_model = ls_linear_model
        end
        push!(hyperparam_results, TrainingResults(learner, linear_model))
    end
    best_validation_reward = -Inf
    best_validation_index = 0
    for (i, result) in enumerate(hyperparam_results)
        reward = result.final_validation_eval.reward
        if reward > best_validation_reward
            best_validation_reward = reward
            best_validation_index = i
        end
    end
    return hyperparam_results[best_validation_index]
end

max_stars_repo_name: joehuchette/reserve-price-optimization

max_issues_repo_name: joehuchette/reserve-price-optimization

max_forks_repo_name: joehuchette/reserve-price-optimization


text:
using CharacterBuffer
using Test

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


text:
@defcomp ClimateTemperature begin
    region = Index()

    # Basic parameters
    area = Parameter(index=[region], unit="km^2")
    y_year_0 = Parameter(unit="year")
    y_year = Parameter(index=[time], unit="year")

    # Climate sensitivity calculations
    tcr_transientresponse = Parameter(unit="degreeC", default=1.70)
    frt_warminghalflife = Parameter(unit="year", default=35.00)

    sens_climatesensitivity = Variable(unit="degreeC")

    # Unadjusted temperature calculations
    fslope_CO2forcingslope = Parameter(unit="W/m2", default=5.5)

    ft_totalforcing = Parameter(index=[time], unit="W/m2")
    fs_sulfateforcing = Parameter(index=[time, region], unit="W/m2")

    et_equilibriumtemperature = Variable(index=[time, region], unit="degreeC")
    rt_realizedtemperature = Variable(index=[time, region], unit="degreeC") # unadjusted temperature

    # Adjusted temperature calculations
    pole_polardifference = Parameter(unit="degreeC", default=1.50) # near 1 degC, the temperature increase difference between equator and pole
    lat_latitude = Parameter(index=[region], unit="degreeLatitude")
    lat_g_meanlatitude = Parameter(unit="degreeLatitude", default=30.21989459076828) # Area-weighted average latitude
    rlo_ratiolandocean = Parameter(unit="unitless", default=1.40) # near 1.4, the ratio between mean land and ocean temperature increases

    rtl_0_realizedtemperature = Parameter(index=[region], unit="degreeC")
    rtl_realizedtemperature = Variable(index=[time, region], unit="degreeC")

    # Global outputs
    rtl_g_landtemperature = Variable(index=[time], unit="degreeC")
    rto_g_oceantemperature = Variable(index=[time], unit="degreeC")
    rt_g_globaltemperature = Variable(index=[time], unit="degreeC")
    rt_g0_baseglobaltemp=Variable(unit="degreeC") #needed for feedback in CO2 cycle component
    rtl_g0_baselandtemp=Variable(unit="degreeC") #needed for feedback in CH4 and N2O cycles

    function init(p, v, d)
        #calculate global baseline temperature from initial regional temperatures

        ocean_prop_ortion = 1. - sum(p.area) / 510000000.

        # Equation 21 from Hope (2006): initial global land temperature
        v.rtl_g0_baselandtemp = sum(p.rtl_0_realizedtemperature' .* p.area') / sum(p.area)

        # initial ocean and global temperatures
        rto_g0_baseoceantemp = v.rtl_g0_baselandtemp/ p.rlo_ratiolandocean
        v.rt_g0_baseglobaltemp = ocean_prop_ortion * rto_g0_baseoceantemp + (1. - ocean_prop_ortion) * v.rtl_g0_baselandtemp
    end

    function run_timestep(p, v, d, tt)

        # Inclusion of transient climate response from Hope (2009)
        if is_first(tt) # only calculate once
            v.sens_climatesensitivity = p.tcr_transientresponse / (1. - (p.frt_warminghalflife / 70.) * (1. - exp(-70. / p.frt_warminghalflife)))
        end

        ## Adjustment for latitude and land
        ocean_prop_ortion = 1. - (sum(p.area) / 510000000.)
        rt_adj_temperatureadjustment = (p.pole_polardifference / 90.) * (abs.(p.lat_latitude) .- p.lat_g_meanlatitude)

        ## Unadjusted realized temperature

        # Equation 19 from Hope (2006): equilibrium temperature estimate
        for rr in d.region
            v.et_equilibriumtemperature[tt, rr] = (v.sens_climatesensitivity / log(2.0)) * (p.ft_totalforcing[tt] + p.fs_sulfateforcing[tt, rr]) / p.fslope_CO2forcingslope
        end

        # Equation 20 from Hope (2006): realized temperature estimate
        # Hope (2009) replaced OCEAN with FRT
        if is_first(tt)
            # Calculate baseline realized temperature by subtracting off adjustment
            rt_0_realizedtemperature = (p.rtl_0_realizedtemperature - rt_adj_temperatureadjustment) * (1. + (ocean_prop_ortion / p.rlo_ratiolandocean) - ocean_prop_ortion)
            for rr in d.region
                v.rt_realizedtemperature[tt, rr] = rt_0_realizedtemperature[rr] + (1 - exp(-(p.y_year[tt] - p.y_year_0) / p.frt_warminghalflife)) * (v.et_equilibriumtemperature[tt, rr] - rt_0_realizedtemperature[rr])
            end
        else
            for rr in d.region
                v.rt_realizedtemperature[tt, rr] = v.rt_realizedtemperature[tt-1, rr] + (1 - exp(-(p.y_year[tt] - p.y_year[tt-1]) / p.frt_warminghalflife)) * (v.et_equilibriumtemperature[tt, rr] - v.rt_realizedtemperature[tt-1, rr])
            end
        end

        ## Adjusted realized temperature

        # Adding adjustment, from Hope (2009)
        for rr in d.region
            v.rtl_realizedtemperature[tt, rr] = v.rt_realizedtemperature[tt, rr] / (1. + (ocean_prop_ortion / p.rlo_ratiolandocean) - ocean_prop_ortion) + rt_adj_temperatureadjustment[rr]
        end

        # Equation 21 from Hope (2006): global realized temperature estimate
        v.rtl_g_landtemperature[tt] = sum(v.rtl_realizedtemperature[tt, :]' .* p.area') / sum(p.area)

        # Ocean and global average temperature from Hope (2009)
        v.rto_g_oceantemperature[tt] = v.rtl_g_landtemperature[tt] / p.rlo_ratiolandocean
        v.rt_g_globaltemperature[tt] = ocean_prop_ortion * v.rto_g_oceantemperature[tt] + (1. - ocean_prop_ortion) * v.rtl_g_landtemperature[tt]
    end
end


text:
mutable struct SquareNeighbors
    E::Bool
    N::Bool
    W::Bool
    S::Bool
end

struct SquareGrid
    vertices::Array{Bool,2}
    neighbors::Array{SquareNeighbors,2}
    location::Tuple{<:Real,<:Real}
    spacing::Real
end

import Base.size
Base.size(G::SquareGrid) = size(G.vertices)

import Base.show
function Base.show(io::IO,G::SquareGrid)
    n = sum(G.vertices)
    E = sum([nb.N+nb.E for nb in G.neighbors])
    vert = n == 1 ? "vertex" : "vertices"
    edge = E == 1 ? "edge" : "edges"
    print(io,"SquareGrid($n $vert, $E $edge)")
end

function vertices(G::SquareGrid)
    m,n = size(G)
    [(i,j) for i=1:m for j=1:n if G.vertices[i,j]]
end

function edges(G::SquareGrid)
    m,n = size(G)
    ([((i,j),(i,j+1)) for i=1:m,j=1:n-1 if G.neighbors[i,j].N] ∪
     [((i,j),(i+1,j)) for i=1:m-1,j=1:n if G.neighbors[i,j].E])
end

import Base.position
function position(G::SquareGrid,i,j)
    return G.location .+ G.spacing.*(i-1,j-1)
end

function rectangulargrid(m::Integer,n::Integer;location=(0,0),spacing=1)
    V = [true for i=1:m,j=1:n]
    N = [SquareNeighbors(true,true,true,true) for i=1:m,j=1:n]
    for i=1:m
        N[i,1].S = false
        N[i,end].N = false
    end
    for j=1:n
        N[1,j].W = false
        N[end,j].E = false
    end
    return SquareGrid(V,N,location,spacing)
end

const STEP = Dict(
    :E => (1,0),
    :N => (0,1),
    :W => (-1,0),
    :S => (0,-1)
)

const OPP = Dict(
    :E => :W,
    :W => :E,
    :N => :S,
    :S => :N
)

function neighbors(G::SquareGrid,i,j)
    [(i,j) .+ STEP[s] for s in (:E,:N,:W,:S)
                            if getfield(G.neighbors[i,j],s)]
end

function between(a::Real,b::Real)
    if a == b
        return 1:0 # empty range
    elseif a < b
        return ceil(Int,a):floor(Int,b)
    else
        return ceil(Int,b):floor(Int,a)
    end
end

function domainapprox(loop::Array{<:AsyPlots.Vec2,1},N::Integer)
    mx = minimum(v.x for v in loop)
    my = minimum(v.y for v in loop)
    Mx = maximum(v.x for v in loop)
    My = maximum(v.y for v in loop)
    location = (mx,my)
    δ = max(Mx-mx,My-my)/(N-1)
    xvals = mx:δ:Mx
    yvals = my:δ:My
    m,n = map(length,(xvals,yvals))
    V = [AsyPlots.iswellinside(AsyPlots.Vec2(xvals[i],yvals[j]),
                loop;epsilon=1e-3) for i=1:m,j=1:n]
    northEdgesHit = NTuple{2,Int64}[]
    eastEdgesHit = NTuple{2,Int64}[]
    for i=1:length(loop)
        P,Q = i < length(loop) ? loop[i:i+1] : [loop[end],loop[1]]
        a,b,c,d = [t/δ for t in (P.x-mx,P.y-my,Q.x-mx,Q.y-my)]
        for i = between(a,c)
            j = floor(Int,(d-b)/(c-a)*(i-a)+b)
            push!(northEdgesHit,(i+1,j+1))
        end
        for j = between(b,d)
            i = floor(Int,(c-a)/(d-b)*(j-b)+a)
            push!(eastEdgesHit,(i+1,j+1))
        end
    end
    northHitSet, eastHitSet = map(Set,(northEdgesHit,eastEdgesHit))
    N = [SquareNeighbors(false,false,false,false) for i=1:m,j=1:n]
    for i=1:m
        for j=1:n
            for dir in (:E,:N)
                nb = (i,j) .+ STEP[dir]
                if (nb[1] ≤ m && nb[2] ≤ n &&
                       V[i,j] && V[nb...]  &&
                        ~(dir == :N && (i,j) in northHitSet) &&
                        ~(dir == :E && (i,j) in eastHitSet))
                    setfield!(N[i,j],dir,true)
                    setfield!(N[nb...],OPP[dir],true)
                end
            end
        end
    end
    return SquareGrid(V,N,location,δ)
end

function isboundary(G::SquareGrid,i,j)
    m,n = size(G)
    if ~G.vertices[i,j]
        return false
    elseif i == 1 || i == m || j == 1 || j == n
        return true
    elseif ~all(G.vertices[i+k,j+l] for (k,l) in values(STEP))
        return true
    else
        return false
    end
end

function boundary(G::SquareGrid)
    m,n = size(G)
    [(i,j) for i=1:m,j=1:n if isboundary(G,i,j)]
end

function solvelaplace(G::SquareGrid,
                      Δvals::Function,
                      bvals::Function,
                      bvertices::Set=Set(boundary(G)))
    m,n = size(G)
    function newentry!(I,J,V,i,j,v)
        push!(I,i); push!(J,j); push!(V,v);
    end
    # number the vertices in the graph, for purposes of
    # setting up a linear system
    vertexmap = Dict(map(reverse,enumerate([(i,j) for j=1:n
                                for i=1:m if G.vertices[i,j]])))
    # I, J are index lists and V a value list for a sparse matrix
    # used to construct the linear system describing the Laplace
    # equation
    I = Int64[]
    J = Int64[]
    V = Float64[]
    # b is a vector representing the right-hand side of the
    # linear system describing the Laplace equation
    b = zeros(length(vertexmap))
    for (v,k) in vertexmap
        if v in bvertices
            # ensure boundary values are satisfied
            newentry!(I,J,V,k,k,1)
            b[k] = bvals(v...)
        else
            nbs = neighbors(G,v...)
            l = length(nbs)
            for nb in nbs
                # row for Laplacian of vertex v
                j = vertexmap[nb]
                newentry!(I,J,V,k,j,1/l)
                newentry!(I,J,V,k,k,-1/l)
            end
            b[k] = Δvals(v...)
        end
    end
    A = SparseArrays.sparse(I,J,V)
    soln = LinearAlgebra.lu(A) \ b
    h = zeros(m,n)
    for (k,v) in vertexmap
        h[k...] = soln[v]
    end
    return h
end

function draw(G::SquareGrid;kwargs...)
    grlist = AsyPlots.GraphicElement2D[]
    for v in vertices(G)
        push!(grlist,AsyPlots.Point(position(G,v...);kwargs...))
    end
    for e in edges(G)
        push!(grlist,AsyPlots.Path([position(G,v...) for v in e];kwargs...))
    end
    AsyPlots.Plot(grlist)
end

function laplacesurface(G::SquareGrid,
                        Δvals::Function,
                        bvals::Function,
                        bvertices::Set=Set(boundary(G)))
    h = solvelaplace(G,Δvals,bvals,bvertices)
    AsyPlots.Plot(AsyPlots.Surface(h,
                                  clip=G.vertices,
                                  spline=false),
                  ignoreaspect=true,zmax=2)
end


text:
#!/usr/bin/env -S julia --compile=min

using VegaLite, VegaDatasets

# p = dataset("cars") |> @vlplot(:bar, x="count()", y=:Origin)

p = dataset("cars") |>
       @vlplot(
           :point,
           x=:Horsepower,
           y=:Miles_per_Gallon,
           color=:Origin,
           width=400,
           height=400
       )

save("myplot.svg", p)

#  Time:
# ./hello.jl  12.42s

# @time using VegaLite, VegaDatasets
# 3.722158 seconds. using
# 2.087367 seconds	vplot
# 3.784914 seconds	save
# 9.4 seconds

# 2nd run:
# 0.004162 seconds vplot (includes near instantaneous open of webpage with plot)
# 2.814209 seconds save
#

# And who can even explain this?:
# back to back:
# julia> @time save("myplot2.svg", p)
#   2.814209 seconds (234.42 k allocations: 7.033 MiB)
#
# julia> @time save("myplot3.svg", p)
#   3.642300 seconds (234.41 k allocations: 7.033 MiB)
#
# julia> @time save("myplot4.svg", p)
#   2.916682 seconds (234.41 k allocations: 7.033 MiB)

max_stars_repo_name: jdurbin/sandbox

max_issues_repo_name: jdurbin/sandbox

max_forks_repo_name: jdurbin/sandbox


text:
# type conversions

## SparseChar conversions

function Base.convert(::Type{cusparseOperation_t}, trans::SparseChar)
    if trans == 'N'
        CUSPARSE_OPERATION_NON_TRANSPOSE
    elseif trans == 'T'
        CUSPARSE_OPERATION_TRANSPOSE
    elseif trans == 'C'
        CUSPARSE_OPERATION_CONJUGATE_TRANSPOSE
    else
        throw(ArgumentError("Unknown operation $trans"))
    end
end

function Base.convert(::Type{cusparseMatrixType_t}, mattype::SparseChar)
    if mattype == 'G'
        CUSPARSE_MATRIX_TYPE_GENERAL
    elseif mattype == 'T'
        CUSPARSE_MATRIX_TYPE_TRIANGULAR
    elseif mattype == 'S'
        CUSPARSE_MATRIX_TYPE_SYMMETRIC
    elseif mattype == 'H'
        CUSPARSE_MATRIX_TYPE_HERMITIAN
    else
        throw(ArgumentError("Unknown matrix type $mattype"))
    end
end

function Base.convert(::Type{cusparseFillMode_t}, uplo::SparseChar)
    if uplo == 'U'
        CUSPARSE_FILL_MODE_UPPER
    elseif uplo == 'L'
        CUSPARSE_FILL_MODE_LOWER
    else
        throw(ArgumentError("Unknown fill mode $uplo"))
    end
end

function Base.convert(::Type{cusparseDiagType_t}, diag::SparseChar)
    if diag == 'U'
        CUSPARSE_DIAG_TYPE_UNIT
    elseif diag == 'N'
        CUSPARSE_DIAG_TYPE_NON_UNIT
    else
        throw(ArgumentError("Unknown diag mode $diag"))
    end
end

function Base.convert(::Type{cusparseIndexBase_t}, index::SparseChar)
    if index == 'Z'
        CUSPARSE_INDEX_BASE_ZERO
    elseif index == 'O'
        CUSPARSE_INDEX_BASE_ONE
    else
        throw(ArgumentError("Unknown index base"))
    end
end

function Base.convert(::Type{cusparseDirection_t}, dir::SparseChar)
    if dir == 'R'
        CUSPARSE_DIRECTION_ROW
    elseif dir == 'C'
        CUSPARSE_DIRECTION_COLUMN
    else
        throw(ArgumentError("Unknown direction $dir"))
    end
end


text:
# This file is a part of Julia. License is MIT: http://julialang.org/license

@testset "Core" begin
    m = [1+im 2; 2 4-im]
    cm = ConjMatrix(m)
    @test cm[1,1] == 1-im
    @test cm[1] == 1-im
    @test trace(cm*m) == 27
    @test cm' == m

    cm[:,2] = [3; 3-im] #setindex! with a vector
    @test conj(cm) == [1+im 3; 2 3+im]

    v = [[1+im], [1-im]]
    cv = ConjVector(v)
    @test cv[1] == [1-im]
end

@testset "RowVector conjugates" begin
    v = [1+im, 1-im]
    rv = v'
    @test (parent(rv) isa ConjArray)
    @test rv' === v

    # Currently, view behavior defaults to only RowVectors.
    @test isa((v').', Vector)
    @test isa((v.')', Vector)
end


text:
using BenchmarkTools
using LocalSearchSolvers

const suite = BenchmarkGroup()

suite["adaptive"] = BenchmarkGroup(["integer", "discrete", "all different"])

# bench for the error functions and predicate of usual constraints
suite["constraints"] = BenchmarkGroup(["error", "predicate", "automatic", "handmade"])
for c in [all_different]
    for i in 0:10
        n = 2^i
        values = rand(1:2n, n)
        suite["constraints"][string(c), length(values)] = @benchmarkable $(c)($values...)
    end
end

# bench for the different problems modeling
suite["problems"] = BenchmarkGroup(["generation"])
for p in [sudoku]
    for size in 2:10
        suite["problems"][string(p), size] = @benchmarkable $(p)($size)
    end
end

## commands to store the tuning parameters
# tune!(suite)
# BenchmarkTools.save("benchmark/params.json", params(suite));

## syntax is loadparams!(group, paramsgroup, fields...)
# loadparams!(suite, BenchmarkTools.load("benchmark/params.json")[1], :evals, :samples);

results = run(suite, verbose = true, seconds = 1)


text:
using AlgamesDriving
using Algames
using StaticArrays
using LinearAlgebra
using MeshCat

# Visualize Scenario
vis = Visualizer()
open(vis)

mutable struct DataCollectionParams
    num_agents::Int
    # One element per agent in each of the arrays below
    lanes::Array
    init_vel::Array
    target_vel::Array
    collision_params::Array
end
# Generate Autonomous Driving Scenario

# Create Roadway
roadway_opts = FourIntersectionRoadwayOptions()
roadway_opts.lane_length = 4.8
roadway_opts.lane_width = 1.2
roadway_opts.turn_radius = 0.175
roadway = build_roadway(roadway_opts)

# Create Dynamics Model
p = 2
model = UnicycleGame(p=p)
params = DataCollectionParams(p, [1,10], [0.3, 0.4], [0.6,0.8], [0.2, 0.25])

players = Vector{Player}(undef, p)
# For each player specify the start and goal state
# west_straight
players[1] = Player(model, roadway.lane[params.lanes[1]],
    Q=Diagonal(SVector{model.ni[1]}([1e-2,1,0.1,5])),
    x0=VehicleState(-2.4, -0.3, 0.0, params.init_vel[1]),
    xf=VehicleState( 2.2, -0.3, 0.0, params.target_vel[1]),
    u_min=SVector{2, Float64}([-0.36, -1.0]),
    u_max=SVector{2, Float64}([0.36, 1.0]),
    v_min=-Inf, # no constraints on the velocity
    v_max=Inf, # no constraints on the velocity
    r_col=params.collision_params[1])

# south_straight
players[2] = Player(model, roadway.lane[params.lanes[2]],
    Q=Diagonal(SVector{model.ni[2]}([1,1e-2,0.1,5])),
    x0=VehicleState(0.3, -3.7, π/2, params.init_vel[2]),
    xf=VehicleState(0.3, 2.2, π/2, params.target_vel[2]),
    u_min=SVector{2, Float64}([-0.36, -1.0]),
    u_max=SVector{2, Float64}([0.36, 1.0]),
    v_min=0.0, # constrained velocity
    v_max=1.05*params.target_vel[2], # constrained velocity
    r_col=params.collision_params[2])

# Create Scenario
sce = Scenario(model, roadway, players)

# # Create Dynamic Game Problem
N = 71 # number of time steps
dt = 0.1 # time step
solver_opts = Options() # solver options
prob = GameProblem(N, dt, sce, solver_opts) # autonomous driving problem
@time newton_solve!(prob) # solve the problem

# Build environment
player_vis_opts = PlayerVisualizationOptions(car_scale=2.5)
set_scenario!(vis, sce, player_vis_opts=player_vis_opts)
# Set camera
set_env!(vis, VehicleState(0.,0.,0.,0.))
# Animate the vehicles
set_traj!(vis, model, sce, prob.pdtraj.pr)

# Plotting vehicles' trajectories and velocities and the solver's progress
plot!(prob.stats)
plot!(model, prob.pdtraj.pr)
plot_velocity_!(model, prob.pdtraj.pr)


text:
module Common
export sqerr,
       sqerr_deriv,
       nlogsig,
       nlogsig_deriv,
       heaviside,
       sigmoid

"""Squared error"""
sqerr(p::Number, y::Number)       = (p - y)^2
sqerr_deriv(p::Number, y::Number) = 2(p - y)

"""Negative logistic sigmoid"""
nlogsig(p::Number, y::Number)       = -log(sigmoid(p * y))
nlogsig_deriv(p::Number, y::Number) = y * sigmoid(p * y) - y

"""Heaviside step function"""
heaviside(p::Number, y::Number) = p == y ? zero(p) : one(p)

"""Sigmoid"""
sigmoid(x::Number) = one(x) / (one(x) + exp(-x))

end


text:
using Test
using Documenter

@test isdefined(Documenter.Writers, :enable_backend)
@test isdefined(Documenter.Writers, :backends_enabled)

@test Documenter.Writers.backends_enabled[:markdown] === false

using DocumenterMarkdown

@test Documenter.Writers.backends_enabled[:markdown] === true


text:
const julia_exts = [".jl",]
const nbviewer_badge = "https://img.shields.io/badge/show-nbviewer-579ACA.svg"
const download_badge = "https://img.shields.io/badge/download-julia-brightgreen.svg"
const julia_footer = raw"""

---

*This page was generated using [DemoCards.jl](https://github.com/johnnychen94/DemoCards.jl) and [Literate.jl](https://github.com/fredrikekre/Literate.jl).*

"""

"""
    struct JuliaDemoCard <: AbstractDemoCard
    JuliaDemoCard(path::String)

Constructs a julia-format demo card from existing julia file `path`.

The julia file is written in [Literate syntax](https://fredrikekre.github.io/Literate.jl/stable/fileformat).

# Fields

Besides `path`, this struct has some other fields:

* `path`: path to the source julia file
* `cover`: path to the cover image
* `id`: cross-reference id
* `title`: one-line description of the demo card
* `author`: author(s) of this demo.
* `date`: the update date of this demo.
* `description`: multi-line description of the demo card
* `julia`: Julia version compatibility
* `hidden`: whether this card is shown in the generated index page
* `notebook`: enable or disable the jupyter notebook generation. Valid values are `true` or `false`.

# Configuration

You can pass additional information by adding a YAML front matter to the julia file.
Supported items are:

* `cover`: an URL or a relative path to the cover image. If not specified, it will use the first available image link, or all-white image if there's no image links.
* `description`: a multi-line description to this file, will be displayed when the demo card is hovered. By default it uses `title`.
* `id`: specify the `id` tag for cross-references. By default it's infered from the filename, e.g., `simple_demo` from `simple demo.md`.
* `title`: one-line description to this file, will be displayed under the cover image. By default, it's the name of the file (without extension).
* `author`: author name. If there are multiple authors, split them with semicolon `;`.
* `date`: any string contents that can be passed to `Dates.DateTime`. For example, `2020-09-13`.
* `julia`: Julia version compatibility. Any string that can be converted to `VersionNumber`
* `hidden`: whether this card is shown in the layout of index page. The default value is `false`.

An example of the front matter (note the leading `#`):

```julia
# ---
# title: passing extra information
# cover: cover.png
# id: non_ambiguious_id
# author: Jane Doe; John Roe
# date: 2020-01-31
# description: this demo shows how you can pass extra demo information to DemoCards package. All these are optional.
# julia: 1.0
# hidden: false
# ---
```

See also: [`MarkdownDemoCard`](@ref DemoCards.MarkdownDemoCard), [`DemoSection`](@ref DemoCards.DemoSection), [`DemoPage`](@ref DemoCards.DemoPage)
"""
mutable struct JuliaDemoCard <: AbstractDemoCard
    path::String
    cover::Union{String, Nothing}
    id::String
    title::String
    description::String
    author::String
    date::DateTime
    julia::VersionNumber
    hidden::Bool
    notebook::Union{Nothing, Bool}
end

function JuliaDemoCard(path::String)::JuliaDemoCard
    # first consturct an incomplete democard, and then load the config
    card = JuliaDemoCard(path, "", "", "", "", "", DateTime(0), JULIA_COMPAT, false, nothing)

    config = parse(card)
    card.cover = load_config(card, "cover"; config=config)
    card.title = load_config(card, "title"; config=config)
    card.date = load_config(card, "date"; config=config)
    card.author = load_config(card, "author"; config=config)
    card.julia = load_config(card, "julia"; config=config)
    # default id requires a title
    card.id    = load_config(card, "id"; config=config)
    # default description requires a title
    card.description = load_config(card, "description"; config=config)
    card.hidden = load_config(card, "hidden"; config=config)

    # Because we want bottom-level cards configuration to have higher priority, figuring
    # out the default `notebook` option needs a reverse walk from top-level page to the
    # bottom-level card.
    if haskey(config, "notebook")
        notebook = config["notebook"]
        if notebook isa Bool
            card.notebook = notebook
        else
            card.notebook = try
                parse(Bool, lowercase(notebook))
            catch err
                @warn "`notebook` option should be either `\"true\"` or `\"false\"`, instead it is: $notebook. Fallback to unconfigured."
                nothing
            end
        end
    end

    return card
end

"""
    save_democards(card_dir::String, card::JuliaDemoCard;
                   project_dir,
                   src,
                   credit,
                   nbviewer_root_url)

process the original julia file and save it.

The processing pipeline is:

1. preprocess and copy source file
2. generate ipynb file
3. generate markdown file
4. insert header and footer to generated markdown file
"""
function save_democards(card_dir::String,
                        card::JuliaDemoCard;
                        credit,
                        nbviewer_root_url,
                        project_dir=Base.source_dir(),
                        src="src",
                        throw_error = false,
                        properties = Dict{String, Any}(),
                        kwargs...)
    isdir(card_dir) || mkpath(card_dir)
    cardname = splitext(basename(card.path))[1]
    md_path = joinpath(card_dir, "$(cardname).md")
    nb_path = joinpath(card_dir, "$(cardname).ipynb")
    src_path = joinpath(card_dir, "$(cardname).jl")

    if VERSION < card.julia
        # It may work, it may not work; I hope it would work.
        @warn "The running Julia version `$(VERSION)` is older than the declared compatible version `$(card.julia)`. You might need to upgrade your Julia."
    end

    card.notebook = if isnothing(card.notebook)
        # Backward compatibility: we used to generate notebooks for all jl files
        op = get(properties, "notebook", "true")
        op = isa(op, Bool) ? op : try
            Base.parse(Bool, lowercase(op))
        catch err
            @warn "`notebook` option should be either `\"true\"` or `\"false\"`, instead it is: $op. Fallback to \"true\"."
            true
        end
    else
        card.notebook
    end
    @assert !isnothing(card.notebook)

    # 1. generating assets
    cp(card.path, src_path; force=true)
    cd(card_dir) do
        try
            foreach(ignored_dirnames) do x
                isdir(x) || mkdir(x)
            end

            # run scripts in a sandbox
            m = Module(gensym())
            # modules created with Module() does not have include defined
            # abspath is needed since this will call `include_relative`
            Core.eval(m, :(include(x) = Base.include($m, abspath(x))))
            gen_assets() = Core.eval(m, :(include($cardname * ".jl")))
            verbose_mode() ? gen_assets() : @suppress gen_assets()

            # WARNING: card.path is modified here
            card.path = cardname*".jl"
            card.cover = load_config(card, "cover")
            card.path = src_path

            foreach(ignored_dirnames) do x
                isempty(readdir(x)) && rm(x; recursive=true, force=true)
            end
        catch err
            if err isa LoadError
                # TODO: provide more informative debug message
                @warn "something wrong during the assets generation for demo $cardname.jl"
                if throw_error
                    rethrow(err)
                else
                    @warn err
                end
            else
                rethrow(err)
            end
        end
    end

    # remove YAML frontmatter
    src_header, _, body = split_frontmatter(read(src_path, String))
    isempty(strip(src_header)) || (src_header *= "\n\n")

    # insert header badge
    badges = make_badges(card;
                         src=src,
                         card_dir=card_dir,
                         nbviewer_root_url=nbviewer_root_url,
                         project_dir=project_dir,
                         build_notebook=card.notebook)
    write(src_path, badges * "\n\n", body)

    # 2. notebook
    if card.notebook
        try
            @suppress Literate.notebook(src_path, card_dir; credit=credit)
        catch err
            nothing # There are already early warning in the assets generation
        end
    end

    # 3. markdown
    @suppress Literate.markdown(src_path, card_dir; credit=false) # manually add credit later
    # remove meta info generated by Literate.jl
    contents = readlines(md_path)
    offsets = findall(x->startswith(x, "```"), contents)
    body = join(contents[offsets[2]+1:end], "\n")

    # 4. insert header and footer to generated markdown file
    config = parse(Val(:Julia), body)
    need_header = !haskey(config, "title")
    # @ref syntax: https://juliadocs.github.io/Documenter.jl/stable/man/syntax/#@ref-link-1
    header = need_header ? "# [$(card.title)](@id $(card.id))\n" : "\n"
    footer = credit ? julia_footer : "\n"
    write(md_path, header, body, footer)

    # 5. filter out source file
    mktempdir(card_dir) do tmpdir
        @suppress Literate.script(src_path, tmpdir; credit=credit)
        write(src_path, src_header, read(joinpath(tmpdir, basename(src_path)), String))
    end
end

function make_badges(card::JuliaDemoCard; src, card_dir, nbviewer_root_url, project_dir, build_notebook)
    cardname = splitext(basename(card.path))[1]
    badges = ["#md #"]
    push!(badges, "[![Source code]($download_badge)]($(cardname).jl)")

    if build_notebook
        if !isempty(nbviewer_root_url)
            # Note: this is only reachable in CI environment
            nbviewer_folder = normpath(relpath(card_dir, "$project_dir/$src"))
            nbviewer_url = replace("$(nbviewer_root_url)/$(nbviewer_folder)/$(cardname).ipynb", Base.Filesystem.path_separator=>'/')
        else
            nbviewer_url = "$(cardname).ipynb"
        end
        push!(badges, "[![notebook]($nbviewer_badge)]($nbviewer_url)")
    end
    if card.julia != JULIA_COMPAT
        # It might be over verbose to insert a compat badge for every julia card, only add one for
        # cards that users should care about
        # Question: Is this too conservative?
        push!(badges, "![compat](https://img.shields.io/badge/julia-$(card.julia)-blue.svg)")
    end

    push!(badges, invoke(make_badges, Tuple{AbstractDemoCard}, card))

    join(badges, " ")
end


text:
using DataFrames
using GraphDataFrameBridge
using MetaGraphs
using CSV
using LightGraphs
using Random
using DataFramesMeta

"""    load_graph_data(input_data::Array{Any,1})

Helper function that allows a synthetic set of vertices or edges to be represented as an array for the purpose of serialization.
"""
function load_graph_data(input_data)

    return input_data
end

"""    gen_rand_vertex_name(vertex_type::String, tag::String)

This function outputs a vertex name that reflects the provided type and tag, and includes a random component to ensure uniqueness.

see also: [`Graphs.generate_synthetic_vertices`](@ref), [`gen_vertex_hash`](@ref)
"""
function gen_rand_vertex_name(vertex_type::String, tag="synth")
    return string(tag, "_", vertex_type, "_", randstring(MersenneTwister(), 'a':'z', 6))
end

"""    gen_vertex_hash(vertex_name::String, vertex_type::String)

This function computes the hash of a vertex's name and type; this combination is assumed to be unique within the graph.

see also: [`generate_synthetic_vertices`](@ref), [`gen_rand_vertex_name`](@ref)
"""
function gen_vertex_hash(vertex_name::String, vertex_type::String)
    return hash(string(vertex_name, vertex_type))
end

# +
"""    generate_synthetic_vertices(vertex_type_defs::String, output_path::String)

Generate synthetic test data. The synthetic vertices are returned as a dataframe
that can be used for testing/debugging/developing the knowledge graph.

see also: [`generate_synthetic_edges`](@ref)
"""
function generate_synthetic_vertices(vertex_type_defs::String, output_path::String)
    vertex_types = CSV.read(vertex_type_defs)
    synth_vertices_df = []
    for vertex in eachrow(vertex_types)
        # Generate a synthetic name
        v_name = gen_rand_vertex_name(vertex.v_type, "synth")
        # Get an id for that name
        v_hash = gen_vertex_hash(v_name, vertex.v_type)
        synth_v_attrs = (v_hash="$v_hash", v_name=v_name, v_type=v_type)
        push!(synth_vertices_df, synth_v_attrs)

    end
    open(output_path, "w") do io
        print(io, repr(load_graph_data(synth_vertices_df)))
    end
    @info("Synthetic vertex dataframe generated and saved as a Julia file.")
    return DataFrame(synth_vertices_df)

end
# -

# """    generate_synthetic_edges(edge_type_defs::String, synth_vertex_df::DataFrame, output_path::String)
#
# Generate synthetic test data. The synthetic edges are returned as a dataframe
# that can be used for testing/debugging/developing the knowledge graph.
#
# see also: [`generate_synthetic_vertices`](@ref)
# """
# function generate_synthetic_edges(edge_type_defs::String, synth_vertex_df::DataFrame, output_path::String)
#
#     edge_types = CSV.read(edge_type_defs)
#
#     synth_edges_df = []
#
#     for edge_row in eachrow(edge_types)
#
#         src_row = @linq synth_vertex_df |>
#             where(:v_type .== edge_row.src_type) |>
#             select(:v_hash, :v_name, :v_type)
#
#         dst_row = @linq synth_vertex_df |>
#             where(:v_type .== edge_row.dst_type) |>
#             select(:v_hash, :v_name, :v_type)
#
#         if size(src_row)[1] >=1 && size(dst_row)[1] >= 1
#             src_vhash = src_row.v_hash[1]
#             src_name = src_row.v_name[1]
#             src_vtype = src_row.v_type[1]
#
#             dst_vhash = dst_row.v_hash[1]
#             dst_name = dst_row.v_name[1]
#             dst_vtype = dst_row.v_type[1]
#
#             if edge_row.value_field != "nothing"
#                 edge_val = ("placeholder")
#             else
#                 edge_val = nothing
#             end
#
#             synth_edge_attrs = (src_vhash="$src_vhash",
#                                 src_name=src_name,
#                                 src_vtype=src_vtype,
#                                 dst_vhash="$dst_vhash",
#                                 dst_name=dst_name,
#                                 dst_vtype=dst_vtype,
#                                 edge_relation=edge_row.edge_relation,
#                                 edge_description=edge_row.description,
#                                 value=edge_val),
#             push!(synth_edges_df, synth_edge_attrs)
#         else
#             continue
#         end
#     end
#
#     open(output_path, "w") do io
#         print(io, repr(load_graph_data(synth_edges_df)))
#     end
#     @info("Synthetic edge dataframe generated", file=output_path)
#     return DataFrame(synth_edges_df)
#
# end

# +
"""    insert_vertices_from_jl(vertices_file::String, input_graph::Nothing)

Ingests and evaulates a Julia file containing vertex information; instantiates an empty knowledge graph and inserts each unique vertex into this graph.

see also: [`insert_edges_from_jl`](@ref)
"""
function insert_vertices_from_jl(vertices_file::String, input_graph::Nothing)
    # The graph is currently empty;
    # objective is to create it from the input vertex and edge files.
    v_df = DataFrame(include(vertices_file))

    G = MetaDiGraph()
    counter = 1::Int64
    set_indexing_prop!(G, :v_hash)
    for v in eachrow(v_df)
        v_hash = gen_vertex_hash(v.v_name, v.v_type)
        v_attrs = Dict(
            :v_name=>v.v_name,
            :v_type=>v.v_type)
        add_vertex!(G)
        set_indexing_prop!(G, counter, :v_hash, "$v_hash")
        set_props!(G, counter, v_attrs)
        @info("Inserting new vertex into graph:", vertex=v)
        counter += 1
    end
    n = nv(G)
    @info("Returning graph G", nv=nv(G))
    return G

end
# -

"""    insert_vertices_from_jl(vertices_file::String, input_graph::MetaDiGraph)

Ingests and evaulates a Julia file containing vertex information; instantiates
an empty knowledge graph and inserts each unique vertex into this graph.

see also: [`insert_edges_from_jl`](@ref)
"""
function insert_vertices_from_jl(vertices_file::String, input_graph::MetaDiGraph)
    v_df = DataFrame(include(vertices_file))
    n_vertices_input = nv(input_graph)
    n_vertices = size(v_df,1) + n_vertices_input

    counter = nv(input_graph) + 1

    G = MetaDiGraph()

    # insert vertices (with properties) from the input graph
    for i in 1:nv(input_graph)
        v_props = props(input_graph, i)
        add_vertex!(G)
        set_indexing_prop!(G, i, :v_hash, v_props[:v_hash])
        set_props!(G, i,
            Dict(:v_name=>v_props[:v_name],
                 :v_type=>v_props[:v_type]))
    end
    @info("input graph have been inserted.", nv=n_vertices_input)
    for v in eachrow(v_df)
        v_hash = gen_vertex_hash(v.v_name, v.v_type)
        try
            G["$v_hash", :v_hash]
            @info("A vertex with this hash value already exists in the graph.", vhash=v_hash)
        catch error
            if isa(error, KeyError) || isa(error, ExceptionError)
            v_attrs = Dict(
                    :v_name=>v.v_name,
                    :v_type=>v.v_type)
                add_vertex!(G)
                set_indexing_prop!(G, counter, :v_hash, "$v_hash")
                set_props!(G, counter, v_attrs)
                @info("Inserting new vertex into graph", v=v)
                counter += 1
            end
        end
    end
    @info("Returning graph G.", nv=nv(G))
    return G
end

"""    copy_input_graph_to_new_graph(input_graph::MetaDiGraph)

Helper function that instantiates a new MetaDiGraph and inserts vertices/edges
from an (existing) input graph.

see also: [`insert_edges_from_jl`](@ref)
"""
function copy_input_graph_to_new_graph(input_graph::MetaDiGraph)
    G = MetaDiGraph()
    # Insert vertices (with properties) from the input graph
    for i in 1:nv(input_graph)
        v_props = props(input_graph, i)
        add_vertex!(G)
        set_indexing_prop!(G, i, :v_hash, v_props[:v_hash])
        set_props!(G, i,
            Dict(:v_name=>v_props[:v_name],
                 :v_type=>v_props[:v_type]))
    end

    # Insert edges (with properties) from the input graph
    for edge in edges(input_graph)
        e_props = props(input_graph, edge)
        add_edge!(G, edge.src, edge.dst)
        set_props!(G, Edge(edge.src, edge.dst), e_props)
    end
    return G
end

# +
"""    insert_edges_from_jl(edges_file::String, input_graph::MetaDiGraph)

Takes as input an existing graph and an edge file. Each edge in the file is
either inserted (if new) or (if already in G), an associated integer weight is
incremented.

see also: [`insert_vertices_from_jl`](@ref), [`copy_input_graph_to_new_graph`](@ref)
"""
function insert_edges_from_jl(edges_file::String, input_graph::MetaDiGraph)
    # Edges may contain vertices that are \notin G (e.g. v \in E_1 \bigcup E_0 may be \emptyset)
    e_df = DataFrame(include(edges_file))

    if nv(input_graph) > 0
        vertices_already_in_G = [get_prop(input_graph, v, :v_hash) for v in vertices(input_graph)]
        G = copy_input_graph_to_new_graph(input_graph)
        set_indexing_prop!(G, :v_hash)
    else
        vertices_already_in_G = []
        G = copy(input_graph)
    end

    vertices_in_edgelist = cat(dims=1, unique(e_df.src_vhash), unique(e_df.dst_vhash))
    refs_to_existing_vertices  = intersect(Set(vertices_already_in_G), Set(vertices_in_edgelist))

    num_vertices_g = nv(input_graph)
    num_vertices_el = length(vertices_in_edgelist)
    inter_cardinality = length(refs_to_existing_vertices)
    union_cardinality = length(union(Set(vertices_already_in_G), Set(vertices_in_edgelist)))

    if num_vertices_g == 0
        counter = 1
    elseif inter_cardinality != num_vertices_g
        counter = num_vertices_g + 1
    end

    @info("The input graph contains $num_vertices_g unique vertices")
    @info("The input edge list refers to $num_vertices_el unique vertices.", nv=num_vertices_el)
    @info("The size of the intersection of these two sets is: $inter_cardinality.", nv=inter_cardinality)

    set_indexing_prop!(G, :v_hash)

    for e in eachrow(e_df)
        edge_attrs = Dict(
                        :e_rel=>e.edge_relation,
                        :e_desc=>e.edge_description,
                        :e_value=>e.value,
                        :weight=>1)

        # The src vertex is \notin G; we need to insert it before we can insert the edge
        if !(e.src_vhash in vertices_already_in_G)
            src_v_hash = e.src_vhash
            src_v_attrs = Dict(
                        :v_name=>e.src_name,
                        :v_type=>e.src_vtype
            )

            add_vertex!(G)
            set_indexing_prop!(G, counter, :v_hash, "$src_v_hash")
            set_props!(G, counter, src_v_attrs)

            counter += 1
            vertices_already_in_G = cat(dims=1, "$src_v_hash", vertices_already_in_G)
            vname = e.src_name
            @info("src vertex $vname was not in G, and has been inserted.", vname=vname)
        end

        # The dst vertex is \notin G; we need to insert it before we can insert the edge
        if !(e.dst_vhash in vertices_already_in_G)
            dst_v_hash = e.dst_vhash
            dst_v_attrs = Dict(
                :v_name=>e.dst_name,
                :v_type=>e.dst_vtype
            )

            add_vertex!(G)
            set_indexing_prop!(G, counter, :v_hash, "$dst_v_hash")
            set_props!(G, counter, dst_v_attrs)

            counter += 1
            vertices_already_in_G = cat(dims=1, "$dst_v_hash", vertices_already_in_G)
            vname = e.dst_name
            @info("dst vertex $vname was not in G, and has been inserted.", vname=vname)
        end

        # All vertices referenced are now \in G; we can insert the edge or increment an existing edge's counter
        src_int_id = G[e.src_vhash, :v_hash]
        dst_int_id = G[e.dst_vhash, :v_hash]

        edge_type = e.edge_relation
        src_name = e.src_name
        dst_name = e.dst_name

        # Add a new edge to the graph
        if !(has_edge(G, src_int_id, dst_int_id))
            add_edge!(G, src_int_id, dst_int_id)
            set_props!(G, Edge(src_int_id, dst_int_id), edge_attrs)
            @info("Inserting directed edge of type $edge_type from $src_name to $dst_name.")
        #  This edge already exists \in G; increment the (frequency) counter
        #  (currently used for weight; we will probably want to think about how
        #  to best normalize these counts)
        else
            new_weight = get_prop(G, Edge(src_int_id, dst_int_id), :weight) + 1
            set_prop!(G, Edge(src_int_id, dst_int_id), :weight, new_weight)
            @info("Incrementing weight of existing directed edge",
                  edge_type=edge_type, weight=new_weight,
                  type=edge_type, src=src_name, dst=dst_name)
        end

    end

    weightfield!(G, :weight)
    num_edges_final = ne(G)
    @info("Returning graph G", nedges=num_edges_final, cardinality=union_cardinality)
    return G

end

# +
"""    insert_edges_from_jl(edges_file::String, input_graph::MetaDiGraph)

Takes as input an existing graph and an edge file. Each edge in the file is
either inserted (if new) or (if already in G), an associated integer weight is
incremented.

see also: [`insert_vertices_from_jl`](@ref), [`copy_input_graph_to_new_graph`](@ref)
"""
function insert_edges_from_jl(edges_file::DataFrame, input_graph::MetaDiGraph)
    # Edges may contain vertices that are \notin G (e.g. v \in E_1 \bigcup E_0 may be \emptyset)
    e_df = DataFrame(edges_file)

        if nv(input_graph) > 0
        vertices_already_in_G = [get_prop(input_graph, v, :v_hash) for v in vertices(input_graph)]
        G = copy_input_graph_to_new_graph(input_graph)
        set_indexing_prop!(G, :v_hash)
    else
        vertices_already_in_G = []
        G = copy(input_graph)
    end

    vertices_in_edgelist = cat(dims=1, unique(e_df.src_vhash), unique(e_df.dst_vhash))
    refs_to_existing_vertices  = intersect(Set(vertices_already_in_G), Set(vertices_in_edgelist))

    num_vertices_g = nv(input_graph)
    num_vertices_el = length(vertices_in_edgelist)
    inter_cardinality = length(refs_to_existing_vertices)
    union_cardinality = length(union(Set(vertices_already_in_G), Set(vertices_in_edgelist)))

    if num_vertices_g == 0
        counter = 1
    elseif inter_cardinality != num_vertices_g
        counter = num_vertices_g + 1
    end

    @info("The input graph contains $num_vertices_g unique vertices")
    @info("The input edge list refers to $num_vertices_el unique vertices.", nv=num_vertices_el)
    @info("The size of the intersection of these two sets is: $inter_cardinality.", nv=inter_cardinality)

    set_indexing_prop!(G, :v_hash)

    for e in eachrow(e_df)
        edge_attrs = Dict(
                        :e_rel=>e.edge_relation,
                        :e_desc=>e.edge_description,
                        :e_value=>e.value,
                        :weight=>1)

        # The src vertex is \notin G; we need to insert it before we can insert the edge
        if !(e.src_vhash in vertices_already_in_G)
            src_v_hash = e.src_vhash
            src_v_attrs = Dict(
                        :v_name=>e.src_name,
                        :v_type=>e.src_vtype
            )

            add_vertex!(G)
            set_indexing_prop!(G, counter, :v_hash, "$src_v_hash")
            set_props!(G, counter, src_v_attrs)

            counter += 1
            vertices_already_in_G = cat(dims=1, "$src_v_hash", vertices_already_in_G)
            vname = e.src_name
            @info("src vertex $vname was not in G, and has been inserted.", vname=vname)
        end

        # The dst vertex is \notin G; we need to insert it before we can insert the edge
        if !(e.dst_vhash in vertices_already_in_G)
            dst_v_hash = e.dst_vhash
            dst_v_attrs = Dict(
                :v_name=>e.dst_name,
                :v_type=>e.dst_vtype
            )

            add_vertex!(G)
            set_indexing_prop!(G, counter, :v_hash, "$dst_v_hash")
            set_props!(G, counter, dst_v_attrs)

            counter += 1
            vertices_already_in_G = cat(dims=1, "$dst_v_hash", vertices_already_in_G)
            vname = e.dst_name
            @info("dst vertex $vname was not in G, and has been inserted.", vname=vname)
        end

        # All vertices referenced are now \in G; we can insert the edge or increment an existing edge's counter
        src_int_id = G[e.src_vhash, :v_hash]
        dst_int_id = G[e.dst_vhash, :v_hash]

        edge_type = e.edge_relation
        src_name = e.src_name
        dst_name = e.dst_name

        # Add a new edge to the graph
        if !(has_edge(G, src_int_id, dst_int_id))
            add_edge!(G, src_int_id, dst_int_id)
            set_props!(G, Edge(src_int_id, dst_int_id), edge_attrs)
            @info("Inserting directed edge of type $edge_type from $src_name to $dst_name.")
        #  This edge already exists \in G; increment the (frequency) counter
        #  (currently used for weight; we will probably want to think about how
        #  to best normalize these counts)
        else
            new_weight = get_prop(G, Edge(src_int_id, dst_int_id), :weight) + 1
            set_prop!(G, Edge(src_int_id, dst_int_id), :weight, new_weight)
            @info("Incrementing weight of existing directed edge",
                  edge_type=edge_type, weight=new_weight,
                  type=edge_type, src=src_name, dst=dst_name)
        end

    end

    weightfield!(G, :weight)
    num_edges_final = ne(G)
    @info("Returning graph G", nedges=num_edges_final, cardinality=union_cardinality)
    return G

end

# +
function main(path_to_v_types="../examples/knowledge_graph/data/kg_vertex_types.csv", path_to_edge_types="../examples/knowledge_graph/data/kg_edge_types.csv")

    @info "Generating dataframe of synthetic vertices (vertex types reflect our schema)"
    synth_vdf = generate_synthetic_vertices(path_to_v_types,"../examples/knowledge_graph/data/synth_kg_vertices.jl")

    @info "Generating dataframe of synthetic edges (edge types reflect our schema)"
    synth_edf = generate_synthetic_edges(path_to_edge_types, synth_vdf, "../examples/knowledge_graph/data/synth_kg_edges.jl")

    @info "Example of inserting vertices to instantiate a KG from scratch"
    g1 = insert_vertices_from_jl("../examples/knowledge_graph/data/synth_kg_vertices.jl", nothing)

    @info "Example of inserting vertices to an existing KG"
    g2 = insert_vertices_from_jl("../examples/knowledge_graph/data/synth_kg_vertices.jl", g1)

    @info "Example of inserting edges to an existing KG"
    g3 = insert_edges_from_jl("../examples/knowledge_graph/data/synth_kg_edges.jl", g2)
end

main()


text:
using GLPlot, GLAbstraction, ModernGL, GeometryTypes
n = 100
h = 1./n
r = h:h:1.
t = (-1:h:1+h)*π
x = map(Float32, r*cos(t)')
y = map(Float32, r*sin(t)')

f(x,y) = exp(-10x.^2-20y.^2)  # arbitrary function of f
z      = Float32[(f(x[k,j],y[k,j])) for k=1:size(x,1),j=1:size(x,2)]
obj    = glplot((x, y, z), :surface)


text:
module TestFinRelations
using Test

using Catlab.Theories, Catlab.CategoricalAlgebra.FinRelations
using Catlab.CategoricalAlgebra.Matrices: MatrixDom

# Category of finite ordinals and relations
###########################################

R = FinRelation(Matrix{BoolRig}([true false true; false true false]))
@test dom(R) == FinRel(3)
@test codom(R) == FinRel(2)

# Evaluation.
@test map(R, [1,3], [1,2]) == [true, false]
@test map(id(FinRel(3)), [1,1], [1,2]) == [true, false]
@test map(FinRelation((x,y) -> 2x == y, 3, 6), [1,1], [2,1]) == [true, false]

# Compatibility of function and matrix representations.
R = FinRelation((x,y) -> x == y || x == 2y || x == 3y, 10, 5)
S = FinRelation((x,y) -> (x+y) % 2 == 0, 5, 10)
A, B = dom(R), dom(S)
m, n = length(A), length(B)
R_mat, S_mat = force(R), force(S)
A_mat, B_mat = MatrixDom{Matrix{BoolRig}}(m), MatrixDom{Matrix{BoolRig}}(n)

@test force(R ⋅ S) == R_mat ⋅ S_mat
@test force(R ⋅ id(codom(R))) == R_mat
@test force(id(dom(R)) ⋅ R) == R_mat
@test force(dagger(R)) == dagger(R_mat)

@test force(R ⊗ S) == R_mat ⊗ S_mat
@test force(braid(A, B)) == FinRelation(braid(A_mat, B_mat))

@test force(R ⊕ S) == R_mat ⊕ S_mat
@test force(swap(A, B)) == FinRelation(swap(A_mat, B_mat))

S = FinRelation((x,y) -> (x+y) % 2 == 0, 10, 5)
S_mat = force(S)
@test force(meet(R, S)) == meet(R_mat, S_mat)
@test force(join(R, S)) == join(R_mat, S_mat)
@test force(mcopy(dom(R))⋅(R⊗S)⋅mmerge(codom(R))) == meet(R_mat, S_mat)
@test force(coplus(dom(R))⋅(R⊕S)⋅plus(codom(R))) == join(R_mat, S_mat)

@test force(top(A, B)) == FinRelation(ones(BoolRig, n, m))
@test force(bottom(A, B)) == FinRelation(zeros(BoolRig, n, m))
@test force(delete(A)⋅create(B)) == FinRelation(ones(BoolRig, n, m))
@test force(cozero(A)⋅zero(B)) == FinRelation(zeros(BoolRig, n, m))

@test force(pair(R, S)) == pair(R_mat, S_mat)
@test force(copair(R, S)) == copair(R_mat, S_mat)

end


text:
using ImageMethodReverb, Random, WAV, LinearAlgebra
Fs = 16e3                     # sampling frequency
T60 = 0.9                     # reverberation time

xs = [1, 2, 2]          #src pos (in meters)
xr = [2, 1.5, 1] #mic pos
L = Lx, Ly, Lz =  4.,4.,4.    # room dimensions in meters
Nt = round(Int,Fs)            #time samples 1 sec

h_se, = rim(xs,xr,L,0.93.*ones(6),Nt,Fs; Rd=0)  # with sweeping echo
h, = rim(xs,xr,L,0.93.*ones(6),Nt,Fs)  # without sweeping echo

wavplay(0.8 .* h_se ./ norm(h,Inf),Fs)
wavplay(0.8 .* h ./ norm(h,Inf),Fs)

using Plots, FFTW
f = rfftfreq(Nt,Fs)
t = range(0; length = Nt, step = 1/Fs )
p1 = plot(t, h_se, xlabel="Time (t)", ylabel="Sound Pressure (Pa)")
plot!(p1, t, h)
p2 = plot(f, 10 .* log10.(abs.(rfft(h_se))),
          xlabel="Frequency (Hz)",
          ylabel="Sound Pressure (dB)",
          xlim=[50;500],
         )
plot!(p2, f, 10 .* log10.(abs.(rfft(h))))
plot(p1,p2)

max_stars_repo_name: nantonel/JuRIM

max_issues_repo_name: nantonel/JuRIM

max_forks_repo_name: nantonel/JuRIM


text:
using LinearAlgebra: cholesky

@testset "Factorization" begin
    @testset "Cholesky" begin
        A = [1 0 0.5;
             0 2 0;
             0.5 0 3]

        B = rand(4)
        B *= B'
        B += I

        K = A ⊗ B

        @test isposdef(K)

        KC = cholesky(K)
        @test isposdef(K)
        @test size(K) == (12, 12)

        C = cholesky(collect(K))

        @test det(KC) ≈ det(C)
        @test logdet(KC) ≈ logdet(C)
        @test collect(inv(KC)) ≈ inv(C)
    end
end


text:
# Align structures and ensembles

export
    kabschalignment,
    displacements,
    align!,
    alignsimple!,
    alignatoms!,
    alignensemble!,
    alignensemblesimple!,
    centroid,
    rmsd,
    ensemblermsds,
    selfalignensemble!

"""
Find the transformations to map the first set of coordinates onto the second.
See Kabsch algorithm for details.
Returns the translation to center the first set, the translation to center the second set
and the rotation to map the first set onto the second.
"""
function kabschalignment(coords_one::Array{Float64}, coords_two::Array{Float64})
    n_coords = size(coords_one, 2)
    @assert n_coords == size(coords_two, 2) "Number of coordinates is different"
    @assert n_coords > 0 "The coordinate arrays are empty"
    # Find the translation to move the coordinate centroids to the origin
    trans_one = sum(coords_one, dims=2) / n_coords
    trans_two = sum(coords_two, dims=2) / n_coords
    p = coords_one - repeat(trans_one, inner=[1,n_coords])
    q = coords_two - repeat(trans_two, inner=[1,n_coords])
    # Find the rotation that maps the coordinates
    cov = p * transpose(q)
    svd_res = svd(cov)
    rotation = svd_res.V * transpose(svd_res.U)
    return trans_one, trans_two, rotation
end

"Return the distances between each point for two sets of aligned coordinates."
function displacements(coords_one::Array{Float64}, coords_two::Array{Float64})
    n_coords = size(coords_one, 2)
    @assert n_coords == size(coords_two, 2) "Number of coordinates is different"
    @assert n_coords > 0 "The coordinate arrays are empty"
    diff = coords_one - coords_two
    sq_diff = diff .* diff
    dists = sqrt.(sum(sq_diff, dims=1))
    return dists[:]
end

"""
Align the first set of atomic coordinates to the second using an iterative procedure.
At each step atoms with a high structural deviation are removed.
If the standard deviation of the structural deviations is greater than `cutoff`, then structural
deviations further than 2 standard deviations from the median are discarded.
This is repeated `n_cycles` times.
If `calpha` is true, only considers the C-alpha atoms.
"""
function align!(coords_one::Array{Float64},
                    coords_two::Array{Float64},
                    atoms::Array{Atom,1};
                    calpha::Bool=true,
                    n_cycles::Integer=defaults["align_cycles"],
                    cutoff::Real=defaults["align_cutoff"])
    n_coords = size(coords_one, 2)
    @assert n_coords == size(coords_two, 2) "Number of coordinates is different"
    @assert n_coords > 0 "The coordinate arrays are empty"
    n_atoms = length(atoms)
    @assert n_coords == n_atoms "Number of coordinates and number of atoms are different"
    if calpha
        inds_to_use = calphaindices(atoms)
    else
        inds_to_use = collect(1:n_atoms)
    end
    n_starting = length(inds_to_use)
    @assert n_starting > 0 "Not aligning on any atoms"
    new_coords = zeros(3, n_atoms)
    for i in 1:n_cycles
        # Find the alignment
        trans_one, trans_two, rot = kabschalignment(getindex(coords_one, collect(1:3), inds_to_use), getindex(coords_two, collect(1:3), inds_to_use))
        # Move the first set to the origin, rotate it then move it to the centroid of the second set
        new_coords = coords_one - repeat(trans_one, inner=[1,n_atoms])
        new_coords = rot * new_coords + repeat(trans_two, inner=[1,n_atoms])
        #println("Iteration ", i, " RMSD: ", round(rmsd(getindex(new_coords, collect(1:3), inds_to_use), getindex(coords_two, collect(1:3), inds_to_use)), 1), " A")
        struc_devs = displacements(getindex(new_coords, collect(1:3), inds_to_use), getindex(coords_two, collect(1:3), inds_to_use))
        sd = sqrt(var(struc_devs))
        dev_cutoff = median(struc_devs) + 2*sd
        inds_to_remove = Int[]
        for j in 1:length(struc_devs)
            if struc_devs[j] > dev_cutoff
                push!(inds_to_remove, j)
            end
        end
        if length(inds_to_remove) > 0 && sd > cutoff && i != n_cycles
            #println("Removing ", length(inds_to_remove), " atoms")
            deleteat!(inds_to_use, inds_to_remove)
            @assert length(inds_to_use) > 0 "Not aligning on any atoms"
        elseif i == n_cycles
            #println("Finishing after ", n_cycles, " iterations")
        else
            #println("Stopping early")
            break
        end
    end
    #println("Alignment has RMSD ", round(rmsd(getindex(new_coords, collect(1:3), inds_to_use), getindex(coords_two, collect(1:3), inds_to_use)), 1), " A based on ", length(inds_to_use), " out of ", n_starting, " atoms")
    coords_one[:] = new_coords[:]
end

"""
Align the first set of atomic coordinates to the second using a single Kabsch rotation.
If `calpha` is true, only considers the C-alpha atoms.
"""
function alignsimple!(coords_one::Array{Float64},
                    coords_two::Array{Float64},
                    atoms::Array{Atom,1};
                    calpha::Bool=true)
    n_coords = size(coords_one, 2)
    @assert n_coords == size(coords_two, 2) "Number of coordinates is different"
    @assert n_coords > 0 "The coordinate arrays are empty"
    n_atoms = length(atoms)
    @assert n_coords == n_atoms "Number of coordinates and number of atoms are different"
    if calpha
        inds_to_use = calphaindices(atoms)
    else
        inds_to_use = collect(1:n_atoms)
    end
    @assert length(inds_to_use) > 0 "Not aligning on any atoms"
    # Find the alignment
    trans_one, trans_two, rot = kabschalignment(getindex(coords_one, collect(1:3), inds_to_use), getindex(coords_two, collect(1:3), inds_to_use))
    # Move the first set to the origin, rotate it then move it to the centroid of the second set
    new_coords = coords_one - repeat(trans_one, inner=[1,n_atoms])
    new_coords = rot * new_coords + repeat(trans_two, inner=[1,n_atoms])
    coords_one[:] = new_coords[:]
end

"""
Align an array of atoms to a set of coordinates or another set of atoms.
Input is list of atoms and coords to align to.
Changes the atom coords.
"""
function alignatoms!(atoms::Array{Atom,1}, coords_ref::Array{Float64})
    n_atoms = length(atoms)
    @assert n_atoms > 0 "There are no atoms in the atom list"
    @assert n_atoms == size(coords_ref, 2) "The number of atoms and reference coordinates differ"
    coords = atomcoords(atoms)
    align!(coords, coords_ref, atoms)
    atomcoords!(atoms, coords)
end

alignatoms!(atoms::Array{Atom,1}, atoms_ref::Array{Atom,1}) = alignatoms!(atoms, atomcoords(atoms_ref))

"Align an ensemble of coordinates to a coordinate set using an iterative procedure for each structure."
function alignensemble!(ensemble::ModelledEnsemble, coords_ref::Array{Float64})
    strucs = ensemble.strucs
    @assert length(strucs) > 0 "No structures in the ensemble"
    n_coords_ref = size(coords_ref, 2)
    @assert n_coords_ref > 0 "No coords in reference coordinate array"
    @assert size(strucs[1].coords, 2) == n_coords_ref "Number of atoms in ensemble and reference structure are different"
    atoms = ensemble.atoms
    for struc in strucs
        align!(struc.coords, coords_ref, atoms)
    end
end

"Align an ensemble of coordinates to a coordinate set using a single Kabsch rotation for each structure."
function alignensemblesimple!(ensemble::ModelledEnsemble, coords_ref::Array{Float64})
    strucs = ensemble.strucs
    @assert length(strucs) > 0 "No structures in the ensemble"
    n_coords_ref = size(coords_ref, 2)
    @assert n_coords_ref > 0 "No coords in reference coordinate array"
    @assert size(strucs[1].coords, 2) == n_coords_ref "Number of atoms different in ensemble and reference structure"
    atoms = ensemble.atoms
    for struc in strucs
        alignsimple!(struc.coords, coords_ref, atoms)
    end
end

"""
Find the average structure in an ensemble.
The ensemble should be aligned.
"""
function centroid(ensemble::ModelledEnsemble)
    strucs = ensemble.strucs
    @assert length(strucs) > 0 "No structures in the ensemble"
    n_coords = size(strucs[1].coords, 2)
    @assert n_coords > 0 "No coordinates for one or more of the structures"
    sum_coords = zeros(3, n_coords)
    for struc in strucs
        @assert n_coords == size(struc.coords, 2) "There are varying numbers of coordinates between structures"
        sum_coords += struc.coords
    end
    sum_coords /= length(strucs)
    return sum_coords
end

"""
Get the RMSD between two identical structures or two sets of coordinates.
Assumes the structures are aligned.
If `calpha` is true, only considers the C-alpha atoms.
"""
function rmsd(coords_one::Array{Float64},
                    coords_two::Array{Float64},
                    atoms::Array{Atom,1};
                    calpha::Bool=true)
    n_coords = size(coords_one, 2)
    @assert n_coords == size(coords_two, 2) "Number of coordinates is different"
    @assert n_coords > 0 "The coordinate arrays are empty"
    n_atoms = length(atoms)
    @assert n_coords == n_atoms "Number of coordinates and number of atoms are different"
    if calpha
        inds_to_use = calphaindices(atoms)
    else
        inds_to_use = collect(1:n_atoms)
    end
    return rmsd(getindex(coords_one, collect(1:3), inds_to_use), getindex(coords_two, collect(1:3), inds_to_use))
end

function rmsd(coords_one::Array{Float64}, coords_two::Array{Float64})
    n_coords = size(coords_one, 2)
    @assert n_coords == size(coords_two, 2) "Number of coordinates is different"
    @assert n_coords > 0 "The coordinate arrays are empty"
    diff = coords_one - coords_two
    sq_diff = diff .* diff
    return sqrt(sum(sq_diff) / n_coords)
end

"""
Get the RMSD between each structure in an ensemble and a reference structure.
"""
function ensemblermsds(ensemble::ModelledEnsemble, atoms_ref::Array{Atom,1})
    @assert length(atoms_ref) == length(ensemble.atoms) "Number of atoms in ensemble and reference are different"
    rmsds = Float64[]
    coords_ref = atomcoords(atoms_ref)
    for struc in ensemble.strucs
        push!(rmsds, rmsd(struc.coords, coords_ref, atoms_ref))
    end
    return rmsds
end

"""
Align an ensemble to itself.
See Bakan and Bahar, PNAS, 2009 for algorithm.
The single Kabsch alignment is used during iteration, then the iterative alignment is done at the end.
Returns the structure the final alignment was on, not the ensemble average (but they will be close).
The cutoff parameter determines the RMSD between two average structures below which iteration is terminated.
"""
function selfalignensemble!(ensemble::ModelledEnsemble; cutoff::Real=defaults["ens_align_cutoff"])
    strucs = ensemble.strucs
    n_strucs = length(strucs)
    @assert n_strucs > 0 "The ensemble does not contain any structures"
    @assert n_strucs > 1 "The ensemble contains 1 structure - need multiple for alignment"
    atoms = ensemble.atoms
    @assert length(atoms) > 0 "The atom list for the ensemble is empty"
    chosen_struc = strucs[rand(1:n_strucs)]
    coords_old = chosen_struc.coords
    alignensemblesimple!(ensemble, coords_old)
    current_rmsd = cutoff+1.0 # Initialised with arbitrary value greater than the cutoff
    counter = 0
    coords_new = Float64[]
    while current_rmsd >= cutoff
        if counter > 1000
            error("Could not align ensemble to itself")
        end
        coords_new = centroid(ensemble)
        alignensemblesimple!(ensemble, coords_new)
        current_rmsd = rmsd(coords_new, coords_old, atoms)
        coords_old = coords_new
        counter += 1
    end
    alignensemble!(ensemble, coords_new)
    println("Aligned ensemble to itself after ", counter, " iteration(s)")
    return coords_new
end


text:
# This code is part of QuantumCircuits.
#
# (C) Copyright Rafał Pracht 2022.
#
# This code is licensed under the Apache License, Version 2.0. You may
# obtain a copy of this license in the LICENSE.txt file in the root directory
# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
#
# Any modifications or derivative works of this code must retain this
# copyright notice, and modified files need to carry a notice indicating
# that they have been altered from the originals.

module Math

using LinearAlgebra

export dagger, unitary_error, matrix_norm, eye, observe_unitary_error,
       min_observe_unitary_error

"The dagger => conj transpose"
dagger(x) = transpose(conj(x))

"Calculate the unitary matrix error"
function unitary_error(exp_unit, mat_unit)
    dif = exp_unit - mat_unit
    return matrix_norm(dif)
end

"Calculate the observe unitary matrix error, the error which we can obserwe
from quantum state"
function observe_unitary_error(exp_unit, mat_unit, index=1)
    # Remove phase
    exp_unit = exp_unit[:] .* exp(-im * angle(exp_unit[index]))
    mat_unit = mat_unit[:] .* exp(-im * angle(mat_unit[index]))

    return unitary_error(exp_unit, mat_unit)
end

"Calculate the minimum observe unitary matrix error."
function min_observe_unitary_error(exp_unit, mat_unit)
    l = length(exp_unit)
    ret = [observe_unitary_error(exp_unit, mat_unit, i) for i in 1:l]
    return minimum(ret)
end

"Calculate the matrix norm"
function matrix_norm(mat)
    return abs(tr(dagger(mat) * mat))
end

"Return the identity matrix."
eye(n) = Matrix(I, n, n)

end  # module Math