text:
#=
This file defines additional functions to return objects of type Distribution. This is
necessary because we specify prior distributions wrt mean and SD
(for beta and gamma-distributed parameters) and ν and σ (for inverse gamma-distributed
parameters). Note these functions are NOT new methods for the Distributions.Beta, etc.
functions, but rather new functions with the same names.
=#
"""
```
BetaAlt(μ::T, σ::T) where {T<:Real}
```

Given μ and σ, calculate α and β and return a Distributions.Beta Distribution object.

### Arguments
`μ`: The mean of the desired distribution
`σ`: The standard deviation of the desired distribution
"""
function BetaAlt(μ::T, σ::T) where {T<:Real}
    α = (1-μ) * μ^2 / σ^2 - μ
    β = α * (1/μ - 1)
    return Distributions.Beta(α, β)
end

"""
```
GammaAlt(μ::T, σ::T) where {T<:Real}
```

Given μ and σ, calculate α and β and return a Distributions.Gamma object.

### Arguments
`μ`: The mean of the desired distribution
`σ`: The standard deviation of the desired distribution
"""
function GammaAlt(μ::T, σ::T) where {T<:Real}
    β = σ^2 / μ
    α = μ / β
    return Distributions.Gamma(α, β)
end

"""
```
mutable struct RootInverseGamma <: Distribution{Univariate, Continuous}
```

If x  ~ RootInverseGamma(ν, τ), then
   x² ~ ScaledInverseChiSquared(ν, τ²)
   x² ~ InverseGamma(ν/2, ντ²/2)

x has mode τ and ν degrees of freedom.
"""
mutable struct RootInverseGamma <: Distribution{Univariate, Continuous}
    ν::Float64
    τ::Float64
end

Distributions.params(d::RootInverseGamma) = (d.ν, d.τ)

"""
```
Distributions.pdf(d::RootInverseGamma, x::T) where {T<:Real}
```

Compute the pdf of a RootInverseGamma distribution at x.
"""
function Distributions.pdf(d::RootInverseGamma, x::T) where {T<:Real}
    (ν, τ) = params(d)
    return 2 * (ν*τ^2/2)^(ν/2) * exp((-ν*τ^2)/(2x^2)) / gamma(ν/2) / x^(ν+1)
end

"""
```
Distributions.logpdf(d::RootInverseGamma, x::T) where {T<:Real}
```

Compute the log pdf of a RootInverseGamma distribution at x.
"""
function Distributions.logpdf(d::RootInverseGamma, x::T) where {T<:Real}
    (ν, τ) = params(d)
    return log(2) - log(gamma(ν/2)) + (ν/2)*log(ν*τ^2/2) - ((ν+1)/2)*log(x^2) - ν*τ^2/(2x^2)
end

"""
```
Distributions.cdf(d::RootInverseGamma, x::T) where {T<:Real}
```

Compute the cdf of a RootInverseGamma distribution at x.
"""
function Distributions.cdf(d::RootInverseGamma, x::T) where {T<:Real}
    (ν, τ) = params(d)
    return 2^(ν/2)# * SpecialFunctions.gamma(ν/2, τ^2*ν/(4x^2)) / gamma(ν/2)
end

"""
```
Distributions.rand(d::RootInverseGamma; cc::T = 1.0) where T <: AbstractFloat
```

Generate a draw from the RootInverseGamma distribution `d`.
"""
function Distributions.rand(d::RootInverseGamma; cc::T = 1.0) where T<:AbstractFloat
    return sqrt(d.ν * d.τ^2 / sum(randn(round(Int,d.ν)).^2))
end

"""
```
DegenerateMvNormal <: Distribution{Multivariate, Continuous}
```

The `DegenerateMvNormal` mutable struct implements a degenerate multivariate normal
distribution. The covariance matrix may not be full rank (hence degenerate).

See [Multivariate normal distribution - Degenerate case](en.wikipedia.org/wiki/Multivariate_normal_distribution#Degenerate_case).
"""
mutable struct DegenerateMvNormal <: Distribution{Multivariate, Continuous}
    μ::Vector        # mean
    σ::Matrix        # standard deviation
    σ_inv::Matrix    # inverse of standard deviation matrix
    λ_vals::Vector   # eigenvalues of σ
    cov::Bool        # indicates if struct contains covariance matrix
    Σ::Matrix        # covariance matrix
end

"""
```
DegenerateMvNormal(μ::Vector, σ::Matrix)
```
Constructor for DegenerateMvNormal type.
"""
function DegenerateMvNormal(μ::Vector, σ::Matrix; stdev::Bool = true)
    if stdev
        return DegenerateMvNormal(μ, σ, Matrix{eltype(μ)}(undef,0,0), Vector{eltype(μ)}(undef,0),
                                  false, Matrix{eltype(μ)}(undef,0,0))
    else
        U, λ_vals, Vt = svd(σ)
        stdevs = Diagonal(sqrt.(λ_vals)) * Vt'
        #Q,R = qr(stdevs)
        return DegenerateMvNormal(μ, stdevs,
                                  Matrix{eltype(μ)}(undef,0,0), Vector{eltype(μ)}(undef,0),
                                  true, σ)
    end
end

"""
```
DegenerateMvNormal(μ::Vector, σ::Matrix,  σ_inv::Matrix, λ_vals::Vector; stdev::Bool = true))
```
Constructor for DegenerateMvNormal type.
"""
function DegenerateMvNormal(μ::Vector, σ::Matrix, σ_inv::Matrix, λ_vals::Vector; stdev::Bool = true)
    return DegenerateMvNormal(μ, σ, σ_inv, λ_vals, false, Matrix{eltype(μ)}(undef,0,0))
end
"""
```
function init_deg_mvnormal(μ::Vector, σ::Matrix)
```
Initializes fields for DegenerateMvNormal type.
"""
function init_deg_mvnormal(μ::Vector, σ::Matrix; stdev::Bool = true)
    if stdev
        U, λ_vals, Vt = svd(σ)
        λ_inv = [λ > 1e-12 ? 1/λ : 0.0 for λ in λ_vals]
        σ_inv = Vt' * Diagonal(λ_inv) * U'
        return DegenerateMvNormal(μ, σ, σ_inv, λ_vals, false, Matrix{eltype(μ)}(undef,0,0))
    else
        return DegenerateMvNormal(μ, σ, Matrix{eltype(μ)}(undef,0,0), Vector{eltype(μ)}(undef,0),
                                  false, Matrix{eltype(μ)}(undef,0,0))
    end
end
"""
```
Distributions.logpdf(d::DegenerateMvNormal)
```

Method bypasses Distributions package implementation of logpdf so as to minimize numerical error.
"""
function Distributions.logpdf(d::DegenerateMvNormal, x::Vector{T}) where T<:Real
    # We need Σ to compute the logpdf, not σ. Here, Σ = σ σ'
    if d.cov
        U, Λ, Vt = svd(d.Σ)
        ind_zero = findall(x -> x > 0.0, Λ)
        return -(length(d.μ[ind_zero]) * log2π + sum(log.(Λ[ind_zero])) + ((x[ind_zero] .- d.μ[ind_zero])'*inv(d.Σ[ind_zero, ind_zero])*(x[ind_zero] .- d.μ[ind_zero]))) / 2.0
    else
        Σ = d.σ * d.σ'

        # b/c Σ may be singular, we need to use a generalized inverse
        # We opt for the Moore-Penrose pseudoinverse
        U, Λ, Vt = svd((Σ + Σ')/2)
        ind_zero = findall(x -> x > 0.0, Λ)
        Σ = Σ[ind_zero, ind_zero]
        return -(length(d.μ[ind_zero]) * log2π + sum(log.(Λ[ind_zero])) + ((x[ind_zero] .- d.μ[ind_zero])'*inv(Σ)*(x[ind_zero] .- d.μ[ind_zero]))) / 2.0
    end
end

"""
```
rank(d::DegenerateMvNormal)
```

Returns the rank of `d.σ`.
"""
function LinearAlgebra.rank(d::DegenerateMvNormal)
    return rank(d.σ)
end

"""
```
length(d::DegenerateMvNormal)
```

Returns the dimension of `d`.
"""
Base.length(d::DegenerateMvNormal) = length(d.μ)

"""
```
Distributions.rand(d::DegenerateMvNormal; cc::T = 1.0) where T<:AbstractFloat
```

Generate a draw from `d` with variance optionally scaled by `cc^2`.
"""
function Distributions.rand(d::DegenerateMvNormal;  cc::T = 1.0) where T<:AbstractFloat
    return d.μ + cc * d.σ * randn(length(d))
end

"""
```
Distributions.rand(d::DegenerateMvNormal, n::Int)
```

Generate `n` draws from `d`. This returns a matrix of size `(length(d), n)`,
where each column is a sample.
"""
function Distributions.rand(d::DegenerateMvNormal, n::Int)
    return d.μ .+ d.σ*randn(length(d), n)
end

"""
```
DegenerateDiagMvTDist <: Distribution{Multivariate, Continuous}
```

The `DegenerateDiagMvTDist` mutable struct implements a degenerate multivariate Student's t
distribution, where the covariance matrix is diagonal. The covariance matrix may
not be full rank (hence degenerate).
"""
mutable struct DegenerateDiagMvTDist <: Distribution{Multivariate, Continuous}
    μ::Vector          # mean
    σ::Matrix          # standard deviation
    ν::Int             # degrees of freedom

    function DegenerateDiagMvTDist(μ::Vector, σ::Matrix, ν::Int)
        ν > 0       ? nothing : error("Degrees of freedom (ν) must be positive")
        isdiag(σ^2) ? nothing : error("Covariance matrix (σ^2) must be diagonal")
        return new(μ, σ, ν)
    end
end

"""
```
LinearAlgebra.rank(d::DegenerateDiagMvTDist)
```

Returns the rank of `d.σ`.
"""
function rank(d::DegenerateDiagMvTDist)
    return rank(d.σ)
end

"""
```
length(d::DegenerateDiagMvTDist)
```

Returns the dimension of `d`.
"""
Base.length(d::DegenerateDiagMvTDist) = length(d.μ)

"""
```
Distributions.rand(d::DegenerateDiagMvTDist)
```

Generate a draw from `d`.
"""
function Distributions.rand(d::DegenerateDiagMvTDist)
    return d.μ + d.σ*rand(TDist(d.ν), length(d))
end

"""
```
Distributions.rand(d::DegenerateDiagMvTDist, n::Int)
```

Generate `n` draws from `d`. This returns a matrix of size `(length(d), n)`,
where each column is a sample.
"""
function Distributions.rand(d::DegenerateDiagMvTDist, n::Int)
    return d.μ .+ d.σ*rand(TDist(d.ν), length(d), n)
end

# Compute the mean μ and standard deviation σ of a RootInverseGamma object.

# A Root Inverse Gamma / Nagasaki Scaled Chi2 Distribution's mean and standard deviation
# can be computed as follows:

#     μ = √(β) * Γ(α - 0.5) / Γ(α)
#     σ = √(β / (α - 1) - μ²)

# where α = ν/2 and β = τ²*ν/2.
function mean(dist::RootInverseGamma)
    α = dist.ν/2
    β = dist.τ^2 * dist.ν/2

    μ = β^(.5) * gamma(α - 0.5) / gamma(α)
    return μ
end

function std(dist::RootInverseGamma)
    μ = mean(dist)
    α = dist.ν/2
    β = dist.τ^2 * dist.ν/2

    σ = (β / (α - 1) - μ^2)^(0.5)
    return σ
end

"""
```
MatrixNormal <: Distribution{Matrixvariate, Continuous}
```

The `MatrixNormal` mutable struct implements a matrixvariate normal
distribution. Note that the matrix must be square.

See [Matrix normal distribution - Degenerate case](en.wikipedia.org/wiki/Matrix_normal_distribution).
"""
mutable struct MatrixNormal <: Distribution{Matrixvariate, Continuous}
    μ::Matrix # mean
    U::Matrix # row variance
    V::Matrix # col variance
    U_sqrt::Matrix # cholesky of U
    V_sqrt::Matrix # choleksy of V
    U_inv::Matrix # inverse of U
    V_inv::Matrix # inverse of V

    function MatrixNormal(μ::Matrix, U::Matrix, V::Matrix)
        size(μ,1) == size(U,1) || error("μ and U must have the same number of rows")
        size(μ,2) == size(V,1) || error("μ and V must have the same number of columns")
        isposdef(U) || error("U is not a positive definite matrix")
        isposdef(V) || error("V is not a positive definite matrix")

        U_sqrt = Matrix(chol(U))
        V_sqrt = Matrix(chol(V))
        U_inv = inv(U)
        V_inv = inv(V)

        return new(μ, U, V, U_sqrt, V_sqrt, U_inv, V_inv)
    end

    function MatrixNormal(μ::Matrix, Σ::Matrix)
        size(μ) == size(Σ) || error("μ and Σ must be the same size")
        isposdef(Σ) || error("Σ is not a positive definite matrix")
        Σ_sqrt = Matrix(chol(U))
        Σ_inv = inv(U)

        return new(μ, Σ, Σ, Σ_sqrt, Σ_sqrt, Σ_inv, Σ_inv)
    end

end

"""
```
size(d::MatrixNormal)
```

Returns the dimension of `d`.
"""
Base.size(d::MatrixNormal) = size(d.μ)

"""
```
Distributions.rand(d::MatrixNormal)
```

Generate a draw from `d`.
"""
function Distributions.rand(d::MatrixNormal)
    return d.μ + d.U_sqrt*randn(size(d.μ))*d.V_sqrt'
end

function mean(d::MatrixNormal)
    return d.μ
end

Distributions.params(d::MatrixNormal) = d.μ, d.U, d.V

"""
```
Distributions.pdf(d::MatrixNormal, x::Matrix)
```

Compute the pdf of a MatrixNormal distribution at x.
"""
function Distributions.pdf(d::MatrixNormal, x::Matrix)
    μ, U, V = params(d)
    n, p = size(μ)
    U_inv = d.U_inv
    V_inv = d.V_inv
    return exp(-.5*trace(V_inv * (x-μ)' * U_inv * (x-μ))) / ((2π)^(n*p/2) * det(U)^(p/2) * det(V)^(n/2))
end

"""
```
Distributions.logpdf(d::MatrixNormal, x::Matrix)
```

Compute the logpdf of a MatrixNormal distribution at x.
"""
function Distributions.logpdf(d::MatrixNormal, x::Matrix)
    μ, U, V = params(d)
    n, p = size(μ)
    U_inv = d.U_inv
    V_inv = d.V_inv
    return -.5*trace(V_inv * (x-μ)' * U_inv * (x-μ)) - log((2π)^(n*p/2) * det(U)^(p/2) * det(V)^(n/2))
end

"""
```
truncmean(Dist::Truncated)
```
Computes the mean of truncated continuous and discrete distributions.
Written by person in this [comment](https://github.com/JuliaStats/Distributions.jl/issues/709).
Thanks, dude!
"""
function truncmean(Dist::Truncated)
    F(x) = x*pdf(Dist,x)
    y = 0.0;

    if typeof(Dist) <: ContinuousDistribution # Continuous distribution
        y = quadgk(F, Dist.lower, Dist.upper)[1]

    else # Discrete distriubtion
        x = ceil(Dist.lower)
        q_max = 1 - 1E-9;
        x_max = min(Dist.upper, quantile(Dist.untruncated, q_max))

        while x < x_max
            y += F(x)
            x += 1
        end
    end
    return y
end

max_stars_repo_name: FRBNY-DSGE/ModelConstructors

max_issues_repo_name: FRBNY-DSGE/ModelConstructors

max_forks_repo_name: FRBNY-DSGE/ModelConstructors


text:
"""
    Develop()

Adds generated packages to the current environment by `dev`ing them.
See the Pkg documentation
[here](https://julialang.github.io/Pkg.jl/v1/managing-packages/#Developing-packages-1)
for more details.
"""
struct Develop <: Plugin end

function posthook(::Develop, ::Template, pkg_dir::AbstractString)
    Pkg.develop(PackageSpec(; path=pkg_dir))
end

ext: jl


text:
module SpatialIndexing

export SpatialIndex, SpatialIndexException, SpatialElem,
    RTree, SimpleSpatialIndex,
    contained_in, intersects_with

include("regions.jl")
include("pool.jl")

include("abstract.jl")

## Spatial Indices

#   SimpleSpatialIndex
include("simple/simpleindex.jl")

#   R-tree
include("rtree/node.jl")
include("rtree/rtree.jl")
include("rtree/check.jl")
include("rtree/query.jl")
include("rtree/split.jl")
include("rtree/insert.jl")
include("rtree/adjust.jl")
include("rtree/delete.jl")
include("rtree/bulk.jl")

end # module


text:
function getModule(::Image{MLImage}) end
name(dataset::Image{MLImage}) = String(nameof(typeof(dataset)))
function url(dataset::Image{MLImage})
    text = "URL address for $dataset can be found in MLDatasets.jl."
    @info text
    return text
end
function checksum(dataset::Image{MLImage})
    text = "Checksum address for $dataset can be found in MLDatasets.jl."
    @info text
    return text
end

target(dataset::Image{MLImage}) = @info "Target for $dataset can be obtained by loading the whole data."

# ---------------------------- Util functions for splitting ---------------------------- */

"""
    getdata(dataset::Image{MLImage}, type::Type{<:Split})

Returns data from dataset. `type` is either `Train` or `Test.`
"""
function getdata(dataset::Image{MLImage}, type::Type{Train})
    datadep = getModule(dataset)
    return datadep.traindata()
end

function getdata(dataset::Image{MLImage}, type::Type{Test})
    datadep = getModule(dataset)
    return datadep.testdata()
end

getdata(dataset::Image{MLImage}, type::Type{Valid}) = error("No valid data in image datasets from MLDatasets.jl.")

"""
    load(dataset::Image, type::Symbol; kwargs...)

Loads dataset of given type.

- `type::Type{<:Split}=Train`: other possible type is `Test`.
"""
load(dataset::Image{<:ImageType}, type::Type{<:Split}=Train) = getdata(dataset, type)

"""
    splits(dataset::Image{MLImage}, data, indeces1, indeces2)

Devides `data` from `dataset` to two parts.
"""
function splits(dataset::Image{MLImage}, data,  indeces1, indeces2)
    datadep = getModule(dataset)
    return datadep.traindata(indeces1), datadep.traindata(indeces2)
end

"""
    postprocess(dataset::Image{MLImage}, data...)

Returns image data.
"""
function postprocess(dataset::Image{MLImage}, data...; kwargs...)
    return data
end

"""
    labels(dataset::MLImage)

Returns train label array for given dataset.
"""
function labels(dataset::MLImage)
    datadep = getModule(dataset)
    return datadep.trainlabels()
end

max_stars_repo_name: GruAna/PreprocessData

max_issues_repo_name: GruAna/PreprocessData

max_forks_repo_name: GruAna/PreprocessData


text:
using AngularMomentumAlgebra
using WignerSymbols
using LinearAlgebra
using AtomicLevels
using HalfIntegers
using Test

include("linear_combinations.jl")
include("wigner_eckart.jl")
include("multipole_expansions.jl")
include("angular_momenta.jl")
include("spherical_tensors.jl")


text:
module ArrayFire

using Cxx
using Base.Meta

import Base: rand, show, randn, ones, diag, eltype, size, elsize,
    sizeof, length, showarray, convert, ndims
import Cxx: CppEnum
export AFArray

# If you have a crash, enable this
const AF_DEBUG = true

function init_library()
	#Libdl.dlopen(Pkg.dir("ArrayFire","deps","build","arrayfire","src","backend","opencl","libafopencl"),Libdl.RTLD_GLOBAL)
	#addHeaderDir(Pkg.dir("ArrayFire","deps","src","arrayfire","include");kind = C_System)
  Libdl.dlopen("/usr/local/lib/libaf.so")
  addHeaderDir("/usr/local/include", kind=C_System)
end
init_library()

__init__() = init_library()

cxx"""
#include <arrayfire.h>
"""

setDevice(i) = icxx"af::setDevice($i);"

const f32 = CppEnum{:af_dtype}(0) # 32-bit floating point values
const c32 = CppEnum{:af_dtype}(1) # 32-bit complex floating point values
const f64 = CppEnum{:af_dtype}(2) # 64-bit complex floating point values
const c64 = CppEnum{:af_dtype}(3) # 64-bit complex floating point values
const b8  = CppEnum{:af_dtype}(4) #  8-bit boolean values
const s32 = CppEnum{:af_dtype}(5) # 32-bit signed integral values
const u32 = CppEnum{:af_dtype}(6) # 32-bit unsigned integral values
const u8  = CppEnum{:af_dtype}(7) #  8-bit unsigned integral values
const s64 = CppEnum{:af_dtype}(8) # 64-bit signed integral values
const u64 = CppEnum{:af_dtype}(9) # 64-bit unsigned integral values

aftype(::Type{Float32})          = f32
aftype(::Type{Complex{Float32}}) = c32
aftype(::Type{Float64})          = f64
aftype(::Type{Complex{Float64}}) = c64
aftype(::Type{Bool})             = b8
aftype(::Type{Int32})            = s32
aftype(::Type{UInt32})           = u32
aftype(::Type{UInt8})            = u8
aftype(::Type{Int64})            = s64
aftype(::Type{UInt64})           = u64

function jltype(af_dtype)
    if af_dtype == f32
        return Float32
    elseif af_dtype == c32
        return Complex{Float32}
    elseif af_dtype == f64
        return Float64
    elseif af_dtype == c64
        return Complex{Float64}
    elseif af_dtype == b8
        return Bool
    elseif af_dtype == s32
        return Int32
    elseif af_dtype == u32
        return UInt32
    elseif af_dtype == u8
        return UInt8
    elseif af_dtype == s64
        return Int64
    elseif af_dtype == u64
        return UInt64
    end
end

function af_promote{T,S}(::Type{T},::Type{S})
    if T == S
        return T
    elseif T == Complex{Float64} || S == Complex{Float64}
        return Complex{Float64}
    elseif T == Complex{Float32} || S == Complex{Float32}
        (T == Float64 || S == Float64) && return Complex{Float64}
        return Complex{Float32}
    elseif T == Float64 || S == Float64
        return Float64
    elseif T == Float32 || S == Float32
        return Float32
    elseif T == UInt64 || S == UInt64
        return UInt64
    elseif T == Int64 || S == Int64
        return Int64
    elseif T == UInt32 || S == UInt32
        return UInt32
    elseif T == Int32 || S == Int32
        return Int32
    elseif T == UInt8 || S == UInt8
        return UInt8
    elseif T == Bool || S == Bool
        return Bool
    else
        return Float32
    end
end

abstract AFAbstractArray{T} <: AbstractArray{T,4}

immutable AFArray{T} <: AFAbstractArray{T}
    array::vcpp"af::array"
    function AFArray(array::vcpp"af::array")
        ret = new(array)
        if AF_DEBUG && backend_eltype(ret) != T
            error("Tried to create AFArray{$T} with array of eltype $(backend_eltype(ret))")
        end
        ret
    end
end
immutable AFSubArray{T} <: AFAbstractArray{T}
    array::vcpp"af::array::array_proxy"
    function AFSubArray(array::vcpp"af::array::array_proxy")
        ret = new(array)
        if AF_DEBUG && backend_eltype(ret) != T
            error("Tried to create AFArray{$T} with array of eltype $(backend_eltype(ret))")
        end
        ret
    end
end

call{T}(::Type{AFAbstractArray{T}}, proxy::vcpp"af::array::array_proxy") =
    AFSubArray{T}(proxy)
call{T}(::Type{AFAbstractArray{T}}, arr::vcpp"af::array") =
    AFArray{T}(arr)

Cxx.cppconvert{T}(x::AFAbstractArray{T}) = x.array

convert{T}(::Type{AFArray{T}}, arr::AFSubArray{T}) =
    AFArray{T}(icxx"(af::array)$arr")

eltype{T}(x::AFAbstractArray{T}) = T
backend_eltype(x::AFAbstractArray) = jltype(icxx"$x.type();")
sizeof{T}(a::AFAbstractArray{T}) = elsize(a) * length(a)

# GPU to Host
function convert{T}(::Type{Array{T}},x::AFAbstractArray{T})
    ret = Array(Uint8,sizeof(x))
    icxx"$x.host($(pointer(ret)));"
    ret = reinterpret(T, ret)
    ret = reshape(ret, size(x)...)
    ret
end
convert{T}(::Type{Array},x::AFAbstractArray{T}) = convert(Array{T},x)

# Host to GPU
function convert{T}(::Type{AFArray{T}}, x::Array{T})
    AFArray{T}(icxx"return af::array{$(dims_to_dim4(size(x))),$(pointer(x)),afHost};")
end
convert{T}(::Type{AFArray}, x::Array{T}) = convert(AFArray{T}, x)

# Show this array, by getting a data pointer to it and using julia's printing
# mechanism. ArrayFire copies internally anyway, so there is nothing to be gained
# by using it's printing
function showarray{T}(io::IO, X::AFAbstractArray{T};
                   header::Bool=true, kwargs...)
    header && print(io, summary(X))
    !isempty(X) && println(io,":")
    showarray(io, convert(Array{T},X); header = false, kwargs...)
end

af_print(X::AFAbstractArray) = icxx"""af::print("",$X);"""

# Translate dimensions. Note that we need to translate between 0 and 1 based
# indexing
ndims(dim4) = icxx"$dim4.ndims();"
ndims(A::AFAbstractArray) = ndims(icxx"$A.dims();")
function dims_to_dim4(dims)
    if length(dims) == 1
        icxx"af::dim4{$(dims[1])};"
    elseif length(dims) == 2
        icxx"af::dim4{$(dims[1]),$(dims[2])};"
    elseif length(dims) == 3
        icxx"af::dim4{$(dims[1]),$(dims[2]),$(dims[3])};"
    elseif length(dims) == 4
        icxx"af::dim4{$(dims[1]),$(dims[2]),$(dims[3]),$(dims[4])};"
    else
        error("Too many dimensions")
    end
end
function dim4_to_dims(dim4)
    d = ndims(dim4)
    t = (Int(icxx"(int)$dim4[0];"),)
    d == 1 && return t
    t = tuple(t...,Int(icxx"(int)$dim4[1];"))
    d == 2 && return t
    t = tuple(t...,Int(icxx"(int)$dim4[2];"))
    d == 3 && return t
    t = tuple(t...,Int(icxx"(int)$dim4[3];"))
    d == 4 && return t
    error("Bad dim4 object")
end

size(x::AFAbstractArray) = dim4_to_dims(icxx"($x).dims();")

rand{T}(::Type{AFArray{T}}, dims...) = AFArray{T}(icxx"af::randu($(dims_to_dim4(dims)),$(aftype(T)));")
randn{T}(::Type{AFArray{T}}, dims...) = AFArray{T}(icxx"af::randn($(dims_to_dim4(dims)),$(aftype(T)));")
eye{T}(::Type{AFArray{T}}, dims...) = AFArray{T}(icxx"af::identity($(dims_to_dim4(dims)),$(aftype(T)));")
diag{T}(x::AFArray{T}, k = 0) = AFArray{T}(icxx"af::diag($(dims_to_dim4(dims)),$k);");

for op in (:sin, :cos, :tan, :asin, :acos, :log, :log1p, :log10, :sqrt, :transpose,
    :exp, :expm1, :erf, :erfc, :cbrt, :lgamma, :transpose)
    @eval Base.($(quot(op))){T}(x::AFAbstractArray{T}) = AFArray{T}(@cxx af::($op)(x.array))
end

Base.gamma{T}(x::AFAbstractArray{T}) = AFArray{T}(@cxx af::tgamma(x))

#
import Base: +, -

# Resolve conflicts
+(x::AFAbstractArray{Bool},y::Bool) = AFArray{Bool}(@cxx +(x.array,y))
+(y::Bool,x::AFAbstractArray{Bool}) = AFArray{Bool}(@cxx +(y,x.array))
-(x::AFAbstractArray{Bool},y::Bool) = AFArray{Bool}(@cxx -(x.array,y))
-(y::Bool,x::AFAbstractArray{Bool}) = AFArray{Bool}(@cxx -(y,x.array))

for (op,cppop) in ((:+,:+),(:(.+),:+),(:-,:-),(:(.-),:-),(:.*,:*),(:./,:/),(:.>>,:>>),(:.<<,:<<))
    @eval function Base.($(quot(op))){T,S}(x::AFAbstractArray{T}, y::AFAbstractArray{S})
        a1 = x.array
        a2 = y.array
        AFArray{af_promote(T,S)}(@cxx ($(cppop))(a1,a2))
    end
    @eval function Base.($(quot(op))){T,S<:Number}(x::AFAbstractArray{T}, y::S)
        a = x.array
        # This is special behavior hardcoded in arrayfire for real floats
        ST = S <: FloatingPoint ? T : S
        AFArray{af_promote(T,ST)}(@cxx ($(cppop))(a, y))
    end
    @eval function Base.($(quot(op))){T,S<:Number}(y::S, x::AFAbstractArray{T})
        a = x.array
        # This is special behavior hardcoded in arrayfire for real floats
        ST = S <: FloatingPoint ? T : S
        AFArray{af_promote(T,ST)}(@cxx ($(cppop))(y, a))
    end
end
# TODO: add! using +=, etc.

import Base: getindex

# Note that we need to translate between 0 and 1 based indexing
to_af_idx(x::Real) = convert(Cint, x-1)
to_af_idx(x::Range) = icxx"af::seq($(first(x))-1,$(last(x))-1,$(step(x)));"
to_af_idx(x::Colon) = icxx"af::span;"
const tis = to_af_idx

IS = Union(Real,Range,Colon)

_getindex{T}(x::AFAbstractArray{T},y::AFAbstractArray)   = icxx"$x($y-1);"
# Avoid an ambiguity
_getindex{T}(x::AFAbstractArray{T},s0::Real)             = icxx"$x($(tis(s0)));"

_getindex{T}(x::AFAbstractArray{T},s0::IS)               = icxx"$x($(tis(s0)));"

_getindex{T}(x::AFAbstractArray{T},s0::IS,s1::IS)        = icxx"$x($(tis(s0)),
                                                                   $(tis(s1)));"

_getindex{T}(x::AFAbstractArray{T},s0::IS,s1::IS,s2::IS) = icxx"$x($(tis(s0)),
                                                                   $(tis(s1)),
                                                                   $(tis(s2)));"
_getindex{T}(x::AFAbstractArray{T},
                    s0::IS,s1::IS,s2::IS, s3::IS)        = icxx"$x($(tis(s0)),
                                                                   $(tis(s1)),
                                                                   $(tis(s2)),
                                                                   $(tis(s3)));"

function getindex{T}(x::AFAbstractArray{T}, idxs...)
    proxy = _getindex(x,idxs...)
    AFSubArray{T}(proxy)
end

function setindex!{T}(x::AFAbstractArray{T}, val, idxs...)
    proxy = _getindex(x,idxs...)
    icxx"$proxy = $val;"
end

# Sorting
import Base: sort

function sort{T}(A::AFArray{T}; rev = false)
    ndims(A) == 1 ||
        error("Must explicitly specify dimension when sorting mutlidimensional array")
    icxx"af::sort($A,0,$(!rev));"
end
sort{T}(A::AFArray{T}, dim; rev = false) = icxx"af::sort($A,$dim,$(!rev));"

# BLAS operations

import Base: dot, A_mul_Bt, At_mul_B, At_mul_Bt, A_mul_Bc,
    Ac_mul_B, Ac_mul_Bc, transpose, ctranspose, transpose!, ctranspose!

dot{T,S}(lhs::AFAbstractArray{T}, rhs::AFAbstractArray{S}) =
    icxx"dot($lhs, $rhs);"

# Matmul
*{T,S}(a::AFAbstractArray{T}, b::AFAbstractArray{S}) =
    AFArray{af_promote(T,S)}(icxx"matmul($a, $b);")
*{T,S,V}(a::AFAbstractArray{T}, b::AFAbstractArray{S}, c::AFAbstractArray{V}) =
    AFArray{af_promote(af_promote(T,S), V)}(icxx"matmul($a, $b, $c);")
*{T,S,V,W}(a::AFAbstractArray{T}, b::AFAbstractArray{S}, c::AFAbstractArray{V}, d::AFAbstractArray{W}) =
    AFArray{af_promte(af_promote(af_promote(T,S), V), W)}(icxx"matmul($a, $b, $c, $d);")

#Image
function loadImage(str::String)
	a = pointer(str)
	return AFArray{Float32}(icxx"loadImage($a);")
end

function _matmul(a::AFAbstractArray, b::AFAbstractArray;
    lhsProp = AF_MAT_NONE, rhsProp = AF_MAT_NONE)
    icxx"matmul($a,$b,$lhsProp,$rhsProp);"
end

# with transpose
A_mul_Bt(a::AFAbstractArray, b::AFAbstractArray) =
    icxx"matmulNT($a,$b);"

At_mul_B(a::AFAbstractArray, b::AFAbstractArray) =
    icxx"matmulTN($a,$b);"

At_mul_Bt(a::AFAbstractArray, b::AFAbstractArray) =
    icxx"matmulTT($a,$b);"

# with complex conjugate
A_mul_Bc(a::AFAbstractArray, b::AFAbstractArray) =
    _matmul(a,b,rhsProp=AF_MAT_CTRANS)
Ac_mul_B(a::AFAbstractArray, b::AFAbstractArray) =
    _matmul(a,b,lhsProp=AF_MAT_CTRANS)
Ac_mul_Bc(a::AFAbstractArray, b::AFAbstractArray) =
    _matmul(a,b,lhsProp=AF_MAT_CTRANS,rhsProp=AF_MAT_CTRANS)

# transpose

transpose(x::AFAbstractArray) = icxx"af::transpose($x);"
ctranspose(x::AFAbstractArray) = icxx"af::transpose($x,true);"

transpose!(x::AFAbstractArray) = icxx"af::transposeInPlace($x);"
ctranspose!(x::AFAbstractArray) = icxx"af::transposeInPlace($x,true);"

# solve

# TODO: The documentation says only AF_MAT_LOWER/AF_MAT_UPPER are supported
# once AF_MAT_(C)TRANS is supported this could be useful for A_rdiv, etc
\(a::AFAbstractArray, b::AFAbstractArray) = icxx"af::solve($a,$b);"

# Factorizations
import Base.LinAlg: chol, chol!, PosDefException, UpperTriangular,
    LowerTriangular, lufact!, lufact

function _chol{T}(a::AFAbstractArray{T}, is_upper::Bool)
    out = similar(a)
    info = icxx"af::cholesky($out,$a,$is_upper);"
    info > 0 && throw(PosDefException(info))
    AFArray{T}(out)
end

function _chol!{T}(a::AFAbstractArray{T}, is_upper::Bool)
    info = icxx"af::choleskyInPlace($a,$is_upper);"
    info > 0 && throw(PosDefException(info))
    a
end

chol(a::AFAbstractArray, ::Type{Val{:U}}) = UpperTriangular(_chol(a), true)
chol(a::AFAbstractArray, ::Type{Val{:L}}) = LowerTriangular(_chol(a), false)
chol(a::AFAbstractArray) = chol(A,Val{:U})

chol!(a::AFAbstractArray, ::Type{Val{:U}}) = UpperTriangular(_chol!(a), true)
chol!(a::AFAbstractArray, ::Type{Val{:L}}) = LowerTriangular(_chol!(a), false)
chol!(a::AFAbstractArray) = chol!(A,Val{:U})

# Fourier Transforms
# TODO: Multidimensional
import Base: fft
function fft{T}(a::AFAbstractArray{T})
	if ndims(a) == 1
		AFArray{Complex{T}}(icxx"fft($a);")
	elseif ndims(a) == 2
		AFArray{Complex{T}}(icxx"fft2($a);")
	end
end

function fft{T<:Complex}(a::AFAbstractArray{T})
	if ndims(a) == 1
		AFArray{T}(icxx"fft($a);")
	elseif ndims(a) == 2
		AFArray{T}(icxx"fft2($a);")
	end
end

# Exception handling
import Base: showerror
@exception function showerror(io::IO, e::rcpp"af::exception")
    println(io, bytestring(icxx"$e.what();"))
end

end # module


text:
module Stats
using Statistics

export mean, max, min, sum, count

mean(values) = Statistics.mean(values)
max(values) = Base.maximum(values)
min(values) = Base.minimum(values)
sum(values) = Base.sum(values)
count(values) = Base.length(values)
end


text:
abstract Shape

immutable Cube <: Shape
    name::String
    size::Tuple{Real,Real,Real}
    panels::Tuple{Int,Int,Int}
    edgewidth::Real
    centered::Bool
    discretized::Bool
    exclude::Array{String} #top: t, bottom: b, sides: p
    function Cube(name;size=(1,1,1),panels=(3,3,3),edgewidth=.1,centered=false,discretized=true,exclude=[])
        return new(name,size,panels,edgewidth,centered,discretized,exclude)
    end
end

immutable Sphere <: Shape
    name::String
    radius::Real
    depth::Int
    function Sphere(name;radius=1,depth=1)
        return new(name,radius,depth)
    end
end

immutable GmshShape <: Shape
    name::String
    file::String
end

function getfile(sh::Cube)
    height  = ["-xh$(sh.size[1])";"-yh$(sh.size[2])";"-zh$(sh.size[3])"]
    panels  = ["-nx$(sh.panels[1])";"-ny$(sh.panels[2])";"-nz$(sh.panels[3])"]
    misc    = ["-e$(sh.edgewidth)";sh.centered?"-o":[];sh.discretized?[]:"-d"]
    exclude = ["-$x" for x in sh.exclude]
    file = strip(readstring(`cubegen_ $height $panels $misc $exclude`))[2:end]
    return "0 cube," * file
end

function getfile(sh::Sphere)
    radius = ["-r$(sh.radius)"]
    misc   = ["-d$(sh.depth)"]
    file = strip(readstring(`spheregen $radius $misc`))[121:end]
    return "0 sphere, $(sh.radius)m radius sphere (d=$(sh.depth))" * file
end

function getfile(sh::GmshShape)
    return sh.file
end

max_stars_repo_name: jebej/FastCap

max_issues_repo_name: jebej/FastCap

max_forks_repo_name: jebej/FastCap


text:
using SciMLBase: issymbollike, sym_to_index, getobserved

Zygote.@adjoint function Base.getindex(VA::ODESolution, sym, j::Int)
    function ODESolution_getindex_pullback(Δ)
        # i = issymbollike(sym) ? sym_to_index(sym, VA) : sym
		i = nothing
        if i === nothing
            getter = getobserved(VA)
            grz = Zygote.pullback(getter, sym, VA.u[j], VA.prob.p, VA.t[j])[2](Δ)
            du = [k == j ? grz[2] : zero(VA.u[1]) for k in 1:length(VA.u)]
            dp = grz[3] # pullback for p
            dprob = remake(VA.prob, p = dp)
            T = eltype(eltype(VA.u))
            N = length(VA.prob.p)
            Δ′ = ODESolution{T, N, typeof(du), Nothing, Nothing, Nothing, Nothing,
                             typeof(dprob), Nothing, Nothing, Nothing}(du, nothing,
                              nothing, nothing, nothing, dprob, nothing, nothing,
                              VA.dense, 0, nothing, VA.retcode)
            (Δ′, nothing, nothing)
        else
            Δ′ = [m == j ? [i == k ? Δ : zero(VA.u[1][1]) for k in 1:length(VA.u[1])] : zero(VA.u[1]) for m in 1:length(VA.u)]
            (Δ′, nothing, nothing)
        end
    end
    VA[sym, j], ODESolution_getindex_pullback
end

@adjoint function Base.getindex(VA::ODESolution, sym)
    function ODESolution_getindex_pullback(Δ)
        i = issymbollike(sym) ? sym_to_index(sym, VA) : sym
        if i === nothing
            throw("Zygote AD of purely-symbolic slicing for observed quantities is not yet supported. Work around this by using `A[sym,i]` to access each element sequentially in the function being differentiated.")
        else
            Δ′ = [ [i == k ? Δ[j] : zero(x[1]) for k in 1:length(x)] for (x, j) in zip(VA.u, 1:length(VA))]
            (Δ′, nothing)
        end
    end
    VA[sym], ODESolution_getindex_pullback
end

# Zygote.@adjoint function Base.getindex(VA::ODESolution, sym, j::Int)
#     function ODESolution_getindex_pullback(Δ)
#         i = SciMLBase.issymbollike(sym) ? SciMLBase.sym_to_index(sym, VA) : sym
#         if i === nothing
#
# 					  zerou = zero(VA.prob.u0)
# 					  _Δ = @. ifelse(Δ == nothing,(zerou,),Δ)
#
# 					  #return (DiffEqBase.build_solution(VA.prob,VA.alg,VA.t,_Δ), nothing, nothing)
# 						return (DiffEqBase.build_solution(VA.prob,VA.alg,VA.t,_Δ), nothing, nothing)
#
# 						getter = SciMLBase.getobserved(VA)
# 						# @show getter
# 						# getter = VA.prob.f.observed
#             grz = Zygote.pullback(getter, sym, VA.u[j], VA.prob.p, VA.t[j])[2](Δ)
# 						# @show grz
#             du = [k == j ? grz[2] : zero(VA.u[1]) for k in 1:length(VA.u)]
#             dp = grz[3] # pullback for p
# 						# @show size(dp)
# 						# dp = dp == nothing ? zeros(eltype(eltype(VA.u)), length(VA.prob.p)) : dp
# 						dprob = remake(VA.prob, p = dp)
#             T = eltype(eltype(VA.u))
#             N = length(VA.prob.p)
#             Δ′ = ODESolution{T, N, typeof(du), Nothing, Nothing, Nothing, Nothing,
#                              typeof(dprob), Nothing, Nothing, Nothing}(du, nothing,
#                               nothing, nothing, nothing, dprob, nothing, nothing,
#                               VA.dense, 0, nothing, VA.retcode)
# 			(Δ′, nothing, nothing)
#         else
#             Δ′ = [m == j ? [i == k ? Δ : zero(VA.u[1][1]) for k in 1:length(VA.u[1])] : zero(VA.u[1]) for m in 1:length(VA.u)]
#             (Δ′, nothing, nothing)
#         end
#     end
#     VA[sym, j], ODESolution_getindex_pullback
# end
#
# Zygote.@adjoint function Base.getindex(VA::DiffEqBase.ODESolution, sym)
#     function ODESolution_getindex_pullback(Δ)
#         i = SciMLBase.issymbollike(sym) ? SciMLBase.sym_to_index(sym, VA) : sym
#         if i === nothing
#             throw("Zygote AD of purely-symbolic slicing for observed quantities is not yet supported. Work around this by using `A[sym,i]` to access each element sequentially in the function being differentiated.")
#         else
#             Δ′ = [ [i == k ? Δ[j] : zero(x[1]) for k in 1:length(x)] for (x, j) in zip(VA.u, 1:length(VA))]
#             (Δ′, nothing)
#         end
#     end
#     VA[sym], ODESolution_getindex_pullback
# end
#

Zygote.@adjoint function Base.getindex(sim::DiffEqBase.EnsembleSolution, i::Int) #where {T,N,S}
	function EnsembleSolution_getindex_pullback(Δ::ODESolution)
		# prob = sim[1].prob
		# du = zeros(eltype(Δ.u[1]),size(Δ.u))
		# dp = zeros(eltype(prob.p),length(sim[1].prob.p))
		# eprob = remake(prob, u=du, p=dp)
		# empty_sol = ODESolution{T, N, typeof(du), Nothing, Nothing, Nothing, Nothing,
		# 										 typeof(eprob), Nothing, Nothing, Nothing}(du, nothing,
		# 											nothing, nothing, nothing, eprob, nothing, nothing,
		# 											Δ.dense, 0, nothing, Δ.retcode)
    # arr = [t == i ? Δ : empty_sol for t in 1:length(sim)]
		# arr = [t == i ? Δ : sim[t] for t in 1:length(sim)]
		arr = [t == i ? Δ : Δ for t in 1:length(sim)]
	  (EnsembleSolution(arr, 0.0, true), nothing)
	end
	sim[i], EnsembleSolution_getindex_pullback
end


text:
module CMBuffers

using Requires
using MPI
using CUDAapi

export CMBuffer
export SingleCMBuffer, DoubleCMBuffer
export get_stage, get_transfer, prepare_transfer!, prepare_stage!

@enum CMBufferKind begin
    SingleCMBuffer
    DoubleCMBuffer
end

###
# Note: We use pinned and device-mapped hostbuffers in double staging
# Potential improvements
# - Investigate if we need to device-map
# - Implement single buffered with pinned + device-mapped hostbuffers

"""
    CMBuffer{T}(::Type{Arr}, kind, dims...; pinned = true)

CUDA/MPI buffer abstracts storage for MPI communication. The buffer is
used for staging data and for MPI transfers. When running on:

  - CPU -- a single buffer is used for staging and MPI transfers can be
    initiated directly to/from it.
  - CUDA -- either:
    - MPI is CUDA-aware: a single buffer on the device for staging and MPI
      transfers can be initiated directly to/from it, or
    - MPI is not CUDA-aware: a double buffering scheme with the staging
      buffer on the device and a transfer buffer on the host

# Arguments
- `T`: element type
- `Arr::Type`: what kind of array to allocate for `stage`
- `kind::CMBufferKind`: either `SingleCMBuffer` or `DoubleCMBuffer`
- `dims...`: dimensions of the array

# Keyword Arguments
- `pin`:  register the `transfer` buffer with CUDA
"""
struct CMBuffer{T, Arr, Buff}
    stage :: Arr     # Same type as Q.data, used for staging
    transfer :: Buff # Union{Nothing,Buff}

    function CMBuffer{T}(::Type{Arr}, kind, dims...; pin=true) where {T, Arr}
        if kind == SingleCMBuffer
            transfer = nothing
        elseif kind == DoubleCMBuffer
            transfer = zeros(T, dims...)
            if pin
                # XXX: make this
                # transfer = register(transfer) to obtain a HostCMBuffer
                # or try:
                # ```
                #   transfer = register(transfer)::HostCMBuffer
                #   stage    = # CuArray(transfer)
                # ```
                register(transfer)
                # XXX: Should this finalizer be attached to the CMBuffer struct?
                finalizer(unregister, transfer)
            end
        else
            error("CMBufferkind $kind is not implemented yet")
        end

        stage = similar(Arr, T, dims...)
        buffer = new{T, typeof(stage), typeof(transfer)}(stage, transfer)

        return buffer
    end
end

function get_stage(buf::CMBuffer)
    return buf.stage
end

function get_transfer(buf::CMBuffer)
    if buf.transfer === nothing
        return buf.stage
    else
        return buf.transfer
    end
end

function prepare_transfer!(buf::CMBuffer)
    if buf.transfer === nothing
        # nothing to do here
    else
        copybuffer!(buf.transfer, buf.stage, async=false)
    end
end

function prepare_stage!(buf::CMBuffer)
    if buf.transfer === nothing
        # nothing to do here
    else
        copybuffer!(buf.stage, buf.transfer, async=false)
    end
end

######
# Internal methods
######
register(x) = nothing
unregister(x) = nothing

"""
    copybuffer!(A, B; async=true)

Copy a buffer from device to host or vice-versa. Internally this uses
`cudaMemcpyAsync` on the `CuDefaultStream`. The keyword argument
`async` determines whether it is asynchronous with regard to the host.
"""
function copybuffer! end

function copybuffer!(A::AbstractArray, B::AbstractArray; async=true)
    copy!(A, B)
end

@init @require CUDAdrv = "c5f51814-7f29-56b8-a69c-e4d8f6be1fde" @require CuArrays = "3a865a2d-5b23-5a0f-bc46-62713ec82fae" begin
    using .CUDAdrv
    using .CuArrays

    import CUDAdrv.Mem

    function register(arr)
        if sizeof(arr) == 0
            #@warn "Array is of size 0. Can't pin register array with CUDA" size(arr) typeof(arr)
            return arr
        end
        GC.@preserve arr begin
            # XXX: is HOSTREGISTER_DEVICEMAP useful?
            Mem.register(Mem.HostBuffer, pointer(arr), sizeof(arr), Mem.HOSTREGISTER_DEVICEMAP)
        end
    end

    function unregister(arr)
        if sizeof(arr) == 0
            return
        end
        GC.@preserve arr begin
            Mem.unregister(Mem.HostBuffer(pointer(arr), sizeof(arr), CuCurrentContext(), true))
        end
    end

    # CUDAdrv.jl throws on CUDA_ERROR_NOT_READY
    function queryStream(stream)
        err = CUDAapi.@runtime_ccall((:cuStreamQuery, CUDAdrv.libcuda), CUDAdrv.CUresult,
                                     (CUDAdrv.CUstream,), stream)

        if err === CUDAdrv.CUDA_ERROR_NOT_READY
            return false
        elseif err === CUDAdrv.CUDA_SUCCESS
            return true
        else
            CUDAdrv.throw_api_error(err)
        end
    end

    """
        friendlysynchronize(stream)

    MPI defines a notion of progress which means that MPI operations
    need the program to call MPI functions (potentially multiple times)
    to make progress and eventually complete. In some implementations,
    progress on one rank may need MPI to be called on another rank.

    As a result blocking by for example calling cudaStreamSynchronize,
    may create a deadlock in some cases because not calling MPI will
    not make other ranks progress.
    """
    function friendlysynchronize(stream)
        status = false
        while !status
            status = queryStream(stream)
            MPI.Iprobe(MPI.MPI_ANY_SOURCE, MPI.MPI_ANY_TAG, MPI.COMM_WORLD)
        end
        return
    end

    function async_copy!(A, B, N, stream)
        GC.@preserve A B begin
            #copyto!(A, B)
            ptrA = pointer(A)
            ptrB = pointer(B)
            unsafe_copyto!(ptrA, ptrB, N, async=true, stream=stream)
        end
    end
    function copybuffer!(A::Array, B::CuArray; async=true)
        @assert sizeof(A) == sizeof(B)
        stream = CuDefaultStream()
        async_copy!(A, B, length(A), stream)
        if !async
            friendlysynchronize(stream)
        end
    end
    function copybuffer!(A::CuArray, B::Array; async=true)
        @assert sizeof(A) == sizeof(B)
        stream = CuDefaultStream()
        async_copy!(A, B, length(A), stream)
        if !async
            friendlysynchronize(stream)
        end
    end
end

end # module

max_stars_repo_name: favba/CLIMA

max_issues_repo_name: favba/CLIMA

max_forks_repo_name: favba/CLIMA


text:
# CodeGen
function start_llvm_module(x::CodeGenerator, ctx::Context, mod_name::String)
    return StartModule(x, ctx, mod_name)
end

get_llvm_module(x::CodeGenerator) = GetModule(x)

release_llvm_module(x::CodeGenerator) = ReleaseModule(x)

get_decl(x::CodeGenerator, s::String) = GetDeclForMangledName(x, s)


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:3] <= 1)

                     @variable(m, -19.736879506320445 <= q <= 22.217442927624916)

                     add_NL_constraint(m, :((-0.3136061325511519 + 0.2821341474705279*(0.7833824073251248 + 0.1390648465658577*$(x[1]) + -0.6617731040629309*$(x[2]) + -0.5974574624055107*$(x[3]))/(1 + exp(-(0.7833824073251248 + 0.1390648465658577*$(x[1]) + -0.6617731040629309*$(x[2]) + -0.5974574624055107*$(x[3])))) + 0.4396178028948987*(-0.7125949557702831 + 0.9109805514170546*$(x[1]) + 0.42115237771853087*$(x[2]) + 0.7134573642319126*$(x[3]))/(1 + exp(-(-0.7125949557702831 + 0.9109805514170546*$(x[1]) + 0.42115237771853087*$(x[2]) + 0.7134573642319126*$(x[3])))) + -0.24849525966922226*(-0.36270963442709325 + 0.07407412170431016*$(x[1]) + -0.6239570881148437*$(x[2]) + -0.6542866380333012*$(x[3]))/(1 + exp(-(-0.36270963442709325 + 0.07407412170431016*$(x[1]) + -0.6239570881148437*$(x[2]) + -0.6542866380333012*$(x[3])))) + 0.5972646324288218*(-0.5451890615285553 + -0.19626106399785304*$(x[1]) + 0.5284606976702149*$(x[2]) + 0.2694706847237054*$(x[3]))/(1 + exp(-(-0.5451890615285553 + -0.19626106399785304*$(x[1]) + 0.5284606976702149*$(x[2]) + 0.2694706847237054*$(x[3])))) + -0.8620507270647968*(0.35158020081654806 + 0.8300503955389886*$(x[1]) + -0.8477304064925764*$(x[2]) + 0.1840810006485709*$(x[3]))/(1 + exp(-(0.35158020081654806 + 0.8300503955389886*$(x[1]) + -0.8477304064925764*$(x[2]) + 0.1840810006485709*$(x[3])))))/(1 + exp(-(-0.3136061325511519 + 0.2821341474705279*(0.7833824073251248 + 0.1390648465658577*$(x[1]) + -0.6617731040629309*$(x[2]) + -0.5974574624055107*$(x[3]))/(1 + exp(-(0.7833824073251248 + 0.1390648465658577*$(x[1]) + -0.6617731040629309*$(x[2]) + -0.5974574624055107*$(x[3])))) + 0.4396178028948987*(-0.7125949557702831 + 0.9109805514170546*$(x[1]) + 0.42115237771853087*$(x[2]) + 0.7134573642319126*$(x[3]))/(1 + exp(-(-0.7125949557702831 + 0.9109805514170546*$(x[1]) + 0.42115237771853087*$(x[2]) + 0.7134573642319126*$(x[3])))) + -0.24849525966922226*(-0.36270963442709325 + 0.07407412170431016*$(x[1]) + -0.6239570881148437*$(x[2]) + -0.6542866380333012*$(x[3]))/(1 + exp(-(-0.36270963442709325 + 0.07407412170431016*$(x[1]) + -0.6239570881148437*$(x[2]) + -0.6542866380333012*$(x[3])))) + 0.5972646324288218*(-0.5451890615285553 + -0.19626106399785304*$(x[1]) + 0.5284606976702149*$(x[2]) + 0.2694706847237054*$(x[3]))/(1 + exp(-(-0.5451890615285553 + -0.19626106399785304*$(x[1]) + 0.5284606976702149*$(x[2]) + 0.2694706847237054*$(x[3])))) + -0.8620507270647968*(0.35158020081654806 + 0.8300503955389886*$(x[1]) + -0.8477304064925764*$(x[2]) + 0.1840810006485709*$(x[3]))/(1 + exp(-(0.35158020081654806 + 0.8300503955389886*$(x[1]) + -0.8477304064925764*$(x[2]) + 0.1840810006485709*$(x[3]))))))) + (-0.9705953045797062 + -0.6259685253929317*(0.7833824073251248 + 0.1390648465658577*$(x[1]) + -0.6617731040629309*$(x[2]) + -0.5974574624055107*$(x[3]))/(1 + exp(-(0.7833824073251248 + 0.1390648465658577*$(x[1]) + -0.6617731040629309*$(x[2]) + -0.5974574624055107*$(x[3])))) + 0.41918630672941015*(-0.7125949557702831 + 0.9109805514170546*$(x[1]) + 0.42115237771853087*$(x[2]) + 0.7134573642319126*$(x[3]))/(1 + exp(-(-0.7125949557702831 + 0.9109805514170546*$(x[1]) + 0.42115237771853087*$(x[2]) + 0.7134573642319126*$(x[3])))) + -0.9388772793667424*(-0.36270963442709325 + 0.07407412170431016*$(x[1]) + -0.6239570881148437*$(x[2]) + -0.6542866380333012*$(x[3]))/(1 + exp(-(-0.36270963442709325 + 0.07407412170431016*$(x[1]) + -0.6239570881148437*$(x[2]) + -0.6542866380333012*$(x[3])))) + -0.18849282727638839*(-0.5451890615285553 + -0.19626106399785304*$(x[1]) + 0.5284606976702149*$(x[2]) + 0.2694706847237054*$(x[3]))/(1 + exp(-(-0.5451890615285553 + -0.19626106399785304*$(x[1]) + 0.5284606976702149*$(x[2]) + 0.2694706847237054*$(x[3])))) + -0.7819718693068158*(0.35158020081654806 + 0.8300503955389886*$(x[1]) + -0.8477304064925764*$(x[2]) + 0.1840810006485709*$(x[3]))/(1 + exp(-(0.35158020081654806 + 0.8300503955389886*$(x[1]) + -0.8477304064925764*$(x[2]) + 0.1840810006485709*$(x[3])))))/(1 + exp(-(-0.9705953045797062 + -0.6259685253929317*(0.7833824073251248 + 0.1390648465658577*$(x[1]) + -0.6617731040629309*$(x[2]) + -0.5974574624055107*$(x[3]))/(1 + exp(-(0.7833824073251248 + 0.1390648465658577*$(x[1]) + -0.6617731040629309*$(x[2]) + -0.5974574624055107*$(x[3])))) + 0.41918630672941015*(-0.7125949557702831 + 0.9109805514170546*$(x[1]) + 0.42115237771853087*$(x[2]) + 0.7134573642319126*$(x[3]))/(1 + exp(-(-0.7125949557702831 + 0.9109805514170546*$(x[1]) + 0.42115237771853087*$(x[2]) + 0.7134573642319126*$(x[3])))) + -0.9388772793667424*(-0.36270963442709325 + 0.07407412170431016*$(x[1]) + -0.6239570881148437*$(x[2]) + -0.6542866380333012*$(x[3]))/(1 + exp(-(-0.36270963442709325 + 0.07407412170431016*$(x[1]) + -0.6239570881148437*$(x[2]) + -0.6542866380333012*$(x[3])))) + -0.18849282727638839*(-0.5451890615285553 + -0.19626106399785304*$(x[1]) + 0.5284606976702149*$(x[2]) + 0.2694706847237054*$(x[3]))/(1 + exp(-(-0.5451890615285553 + -0.19626106399785304*$(x[1]) + 0.5284606976702149*$(x[2]) + 0.2694706847237054*$(x[3])))) + -0.7819718693068158*(0.35158020081654806 + 0.8300503955389886*$(x[1]) + -0.8477304064925764*$(x[2]) + 0.1840810006485709*$(x[3]))/(1 + exp(-(0.35158020081654806 + 0.8300503955389886*$(x[1]) + -0.8477304064925764*$(x[2]) + 0.1840810006485709*$(x[3]))))))) + (0.8285261526161638 + 0.6863822914423507*(0.7833824073251248 + 0.1390648465658577*$(x[1]) + -0.6617731040629309*$(x[2]) + -0.5974574624055107*$(x[3]))/(1 + exp(-(0.7833824073251248 + 0.1390648465658577*$(x[1]) + -0.6617731040629309*$(x[2]) + -0.5974574624055107*$(x[3])))) + -0.7060930690780838*(-0.7125949557702831 + 0.9109805514170546*$(x[1]) + 0.42115237771853087*$(x[2]) + 0.7134573642319126*$(x[3]))/(1 + exp(-(-0.7125949557702831 + 0.9109805514170546*$(x[1]) + 0.42115237771853087*$(x[2]) + 0.7134573642319126*$(x[3])))) + -0.2956343589188948*(-0.36270963442709325 + 0.07407412170431016*$(x[1]) + -0.6239570881148437*$(x[2]) + -0.6542866380333012*$(x[3]))/(1 + exp(-(-0.36270963442709325 + 0.07407412170431016*$(x[1]) + -0.6239570881148437*$(x[2]) + -0.6542866380333012*$(x[3])))) + 0.32353098534638214*(-0.5451890615285553 + -0.19626106399785304*$(x[1]) + 0.5284606976702149*$(x[2]) + 0.2694706847237054*$(x[3]))/(1 + exp(-(-0.5451890615285553 + -0.19626106399785304*$(x[1]) + 0.5284606976702149*$(x[2]) + 0.2694706847237054*$(x[3])))) + 0.6963022628505415*(0.35158020081654806 + 0.8300503955389886*$(x[1]) + -0.8477304064925764*$(x[2]) + 0.1840810006485709*$(x[3]))/(1 + exp(-(0.35158020081654806 + 0.8300503955389886*$(x[1]) + -0.8477304064925764*$(x[2]) + 0.1840810006485709*$(x[3])))))/(1 + exp(-(0.8285261526161638 + 0.6863822914423507*(0.7833824073251248 + 0.1390648465658577*$(x[1]) + -0.6617731040629309*$(x[2]) + -0.5974574624055107*$(x[3]))/(1 + exp(-(0.7833824073251248 + 0.1390648465658577*$(x[1]) + -0.6617731040629309*$(x[2]) + -0.5974574624055107*$(x[3])))) + -0.7060930690780838*(-0.7125949557702831 + 0.9109805514170546*$(x[1]) + 0.42115237771853087*$(x[2]) + 0.7134573642319126*$(x[3]))/(1 + exp(-(-0.7125949557702831 + 0.9109805514170546*$(x[1]) + 0.42115237771853087*$(x[2]) + 0.7134573642319126*$(x[3])))) + -0.2956343589188948*(-0.36270963442709325 + 0.07407412170431016*$(x[1]) + -0.6239570881148437*$(x[2]) + -0.6542866380333012*$(x[3]))/(1 + exp(-(-0.36270963442709325 + 0.07407412170431016*$(x[1]) + -0.6239570881148437*$(x[2]) + -0.6542866380333012*$(x[3])))) + 0.32353098534638214*(-0.5451890615285553 + -0.19626106399785304*$(x[1]) + 0.5284606976702149*$(x[2]) + 0.2694706847237054*$(x[3]))/(1 + exp(-(-0.5451890615285553 + -0.19626106399785304*$(x[1]) + 0.5284606976702149*$(x[2]) + 0.2694706847237054*$(x[3])))) + 0.6963022628505415*(0.35158020081654806 + 0.8300503955389886*$(x[1]) + -0.8477304064925764*$(x[2]) + 0.1840810006485709*$(x[3]))/(1 + exp(-(0.35158020081654806 + 0.8300503955389886*$(x[1]) + -0.8477304064925764*$(x[2]) + 0.1840810006485709*$(x[3]))))))) + (-0.2246585184074279 + 0.8591320330102441*(0.7833824073251248 + 0.1390648465658577*$(x[1]) + -0.6617731040629309*$(x[2]) + -0.5974574624055107*$(x[3]))/(1 + exp(-(0.7833824073251248 + 0.1390648465658577*$(x[1]) + -0.6617731040629309*$(x[2]) + -0.5974574624055107*$(x[3])))) + -0.5037522661683527*(-0.7125949557702831 + 0.9109805514170546*$(x[1]) + 0.42115237771853087*$(x[2]) + 0.7134573642319126*$(x[3]))/(1 + exp(-(-0.7125949557702831 + 0.9109805514170546*$(x[1]) + 0.42115237771853087*$(x[2]) + 0.7134573642319126*$(x[3])))) + 0.9621050832972355*(-0.36270963442709325 + 0.07407412170431016*$(x[1]) + -0.6239570881148437*$(x[2]) + -0.6542866380333012*$(x[3]))/(1 + exp(-(-0.36270963442709325 + 0.07407412170431016*$(x[1]) + -0.6239570881148437*$(x[2]) + -0.6542866380333012*$(x[3])))) + -0.1995700954953512*(-0.5451890615285553 + -0.19626106399785304*$(x[1]) + 0.5284606976702149*$(x[2]) + 0.2694706847237054*$(x[3]))/(1 + exp(-(-0.5451890615285553 + -0.19626106399785304*$(x[1]) + 0.5284606976702149*$(x[2]) + 0.2694706847237054*$(x[3])))) + 0.442251462608958*(0.35158020081654806 + 0.8300503955389886*$(x[1]) + -0.8477304064925764*$(x[2]) + 0.1840810006485709*$(x[3]))/(1 + exp(-(0.35158020081654806 + 0.8300503955389886*$(x[1]) + -0.8477304064925764*$(x[2]) + 0.1840810006485709*$(x[3])))))/(1 + exp(-(-0.2246585184074279 + 0.8591320330102441*(0.7833824073251248 + 0.1390648465658577*$(x[1]) + -0.6617731040629309*$(x[2]) + -0.5974574624055107*$(x[3]))/(1 + exp(-(0.7833824073251248 + 0.1390648465658577*$(x[1]) + -0.6617731040629309*$(x[2]) + -0.5974574624055107*$(x[3])))) + -0.5037522661683527*(-0.7125949557702831 + 0.9109805514170546*$(x[1]) + 0.42115237771853087*$(x[2]) + 0.7134573642319126*$(x[3]))/(1 + exp(-(-0.7125949557702831 + 0.9109805514170546*$(x[1]) + 0.42115237771853087*$(x[2]) + 0.7134573642319126*$(x[3])))) + 0.9621050832972355*(-0.36270963442709325 + 0.07407412170431016*$(x[1]) + -0.6239570881148437*$(x[2]) + -0.6542866380333012*$(x[3]))/(1 + exp(-(-0.36270963442709325 + 0.07407412170431016*$(x[1]) + -0.6239570881148437*$(x[2]) + -0.6542866380333012*$(x[3])))) + -0.1995700954953512*(-0.5451890615285553 + -0.19626106399785304*$(x[1]) + 0.5284606976702149*$(x[2]) + 0.2694706847237054*$(x[3]))/(1 + exp(-(-0.5451890615285553 + -0.19626106399785304*$(x[1]) + 0.5284606976702149*$(x[2]) + 0.2694706847237054*$(x[3])))) + 0.442251462608958*(0.35158020081654806 + 0.8300503955389886*$(x[1]) + -0.8477304064925764*$(x[2]) + 0.1840810006485709*$(x[3]))/(1 + exp(-(0.35158020081654806 + 0.8300503955389886*$(x[1]) + -0.8477304064925764*$(x[2]) + 0.1840810006485709*$(x[3]))))))) + (0.5761249271362421 + -0.042724357144423575*(0.7833824073251248 + 0.1390648465658577*$(x[1]) + -0.6617731040629309*$(x[2]) + -0.5974574624055107*$(x[3]))/(1 + exp(-(0.7833824073251248 + 0.1390648465658577*$(x[1]) + -0.6617731040629309*$(x[2]) + -0.5974574624055107*$(x[3])))) + -0.4537006876708345*(-0.7125949557702831 + 0.9109805514170546*$(x[1]) + 0.42115237771853087*$(x[2]) + 0.7134573642319126*$(x[3]))/(1 + exp(-(-0.7125949557702831 + 0.9109805514170546*$(x[1]) + 0.42115237771853087*$(x[2]) + 0.7134573642319126*$(x[3])))) + -0.12519521024691826*(-0.36270963442709325 + 0.07407412170431016*$(x[1]) + -0.6239570881148437*$(x[2]) + -0.6542866380333012*$(x[3]))/(1 + exp(-(-0.36270963442709325 + 0.07407412170431016*$(x[1]) + -0.6239570881148437*$(x[2]) + -0.6542866380333012*$(x[3])))) + -0.7426009177864223*(-0.5451890615285553 + -0.19626106399785304*$(x[1]) + 0.5284606976702149*$(x[2]) + 0.2694706847237054*$(x[3]))/(1 + exp(-(-0.5451890615285553 + -0.19626106399785304*$(x[1]) + 0.5284606976702149*$(x[2]) + 0.2694706847237054*$(x[3])))) + -0.8758186284143128*(0.35158020081654806 + 0.8300503955389886*$(x[1]) + -0.8477304064925764*$(x[2]) + 0.1840810006485709*$(x[3]))/(1 + exp(-(0.35158020081654806 + 0.8300503955389886*$(x[1]) + -0.8477304064925764*$(x[2]) + 0.1840810006485709*$(x[3])))))/(1 + exp(-(0.5761249271362421 + -0.042724357144423575*(0.7833824073251248 + 0.1390648465658577*$(x[1]) + -0.6617731040629309*$(x[2]) + -0.5974574624055107*$(x[3]))/(1 + exp(-(0.7833824073251248 + 0.1390648465658577*$(x[1]) + -0.6617731040629309*$(x[2]) + -0.5974574624055107*$(x[3])))) + -0.4537006876708345*(-0.7125949557702831 + 0.9109805514170546*$(x[1]) + 0.42115237771853087*$(x[2]) + 0.7134573642319126*$(x[3]))/(1 + exp(-(-0.7125949557702831 + 0.9109805514170546*$(x[1]) + 0.42115237771853087*$(x[2]) + 0.7134573642319126*$(x[3])))) + -0.12519521024691826*(-0.36270963442709325 + 0.07407412170431016*$(x[1]) + -0.6239570881148437*$(x[2]) + -0.6542866380333012*$(x[3]))/(1 + exp(-(-0.36270963442709325 + 0.07407412170431016*$(x[1]) + -0.6239570881148437*$(x[2]) + -0.6542866380333012*$(x[3])))) + -0.7426009177864223*(-0.5451890615285553 + -0.19626106399785304*$(x[1]) + 0.5284606976702149*$(x[2]) + 0.2694706847237054*$(x[3]))/(1 + exp(-(-0.5451890615285553 + -0.19626106399785304*$(x[1]) + 0.5284606976702149*$(x[2]) + 0.2694706847237054*$(x[3])))) + -0.8758186284143128*(0.35158020081654806 + 0.8300503955389886*$(x[1]) + -0.8477304064925764*$(x[2]) + 0.1840810006485709*$(x[3]))/(1 + exp(-(0.35158020081654806 + 0.8300503955389886*$(x[1]) + -0.8477304064925764*$(x[2]) + 0.1840810006485709*$(x[3]))))))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
module escape_num
export escape_exit_num
using DifferentialEquations
using Roots
using Printf
# using MATLAB
# using ProgressMeter

include("leap_frog_definitions.jl")
max_hit_q1 = 1000
max_hit_p1=1000
barrier = 5

condition_max_hits(u, t, integrator) = u[3] > max_hit_q1 || u[6]>max_hit_p1 || maximum([abs(u[1]),abs(u[2]),abs(u[4]),abs(u[5])]) > barrier
affect_stop!(integrator) = terminate!(integrator)
function condition_hits_PSS_q1(u, t, integrator) # Event when event_f(u,t) == 0
    u[1]
end

function condition_hits_PSS_p2(u, t, integrator) # Event when event_f(u,t) == 0
    u[5]
end

function affect_update_iterator_q1!(integrator)
    integrator.u[3] = integrator.u[3] + 1
end
function affect_update_iterator_p2!(integrator)
    integrator.u[6] = integrator.u[6] + 1
end

callback_max_hits = DiscreteCallback(condition_max_hits, affect_stop!)
callback_hits_PSS_q1 = ContinuousCallback(condition_hits_PSS_q1, affect_update_iterator_q1!, nothing)
callback_hits_PSS_p2 = ContinuousCallback(condition_hits_PSS_p2, affect_update_iterator_p2!, nothing)

cb = CallbackSet(callback_hits_PSS_q1, callback_hits_PSS_p2, callback_max_hits)

function escape_exit_num(mesh_list,t_end,Energy_A)
    Q=mesh_list[2]; P=mesh_list[1];
    H=(2*Energy_A)^2
    tol_dist=1e-5
    Y=Yfind(Q,P,H)

    if ~isempty(Y)
        u0=zeros(5)
        u0[1]=0 #X
        u0[2]=P #P
        u0[3]=Q #Q
        u0[4]=Y #Y
        u0[5]=0 #Y

        prob = ODEProblem(Eq_of_M,u0,(0., t_end))
        # sol=solve(prob,RK4(),maxiters=1e20, reltol=1e-6,abstol=1e-8,callback=cb)
        sol=solve(prob, RK4(),maxiters=1e20, reltol=1e-8,abstol=1e-10,callback=cb,save_start=true,save_end=true,save_everystep=false)
        uf=zeros(5)
        uf[1]=sol[1,end] #X
        uf[2]=sol[2,end] #P
        uf[3]=sol[3,end] #Q
        uf[4]=sol[4,end] #Y
        uf[5]=0
        dH=abs(H_test(uf)-H)
        if sol.u[end][5]==max_hit
            return t_end
        else
            return sol.t[end]

        end
        # if dH<1e-5
        #     if is_it(uf,H,tol_dist)
        #         return sol.u[end][5]
        #     else
        #         return -1
        #     end
        # end
    else
        return 0.0
    end
end

end


text:
"""
connect_mesh(EToV,fv)

Initialize quad elements connectivity matrices, element to element and
element to face connectivity. Uses fv = array of indices of face vertices

```jldoctest
```
"""
function connect_mesh(EToV,fv)
        Nfaces = length(fv)
        K = size(EToV,1)

        # sort and find matches
        fnodes = [[sort(EToV[e,ids]) for ids = fv, e = 1:K]...]
        p = sortperm(fnodes) # sorts by lexicographic ordering by default
        fnodes = fnodes[p,:]

        FToF = reshape(collect(1:Nfaces*K),Nfaces,K)
        for f = 1:size(fnodes,1)-1
                if fnodes[f,:]==fnodes[f+1,:]
                        f1 = FToF[p[f]]
                        f2 = FToF[p[f+1]]
                        FToF[p[f]] = f2
                        FToF[p[f+1]] = f1
                end
        end
        return FToF
end

max_stars_repo_name: yiminllin/ESDG-Fourier

max_issues_repo_name: yiminllin/ESDG-Fourier

max_forks_repo_name: yiminllin/ESDG-Fourier


text:
# syntax: proto3
using ProtoBuf
import ProtoBuf.meta

mutable struct CostGraphDef_Node_InputInfo <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function CostGraphDef_Node_InputInfo(; kwargs...)
        obj = new(meta(CostGraphDef_Node_InputInfo), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct CostGraphDef_Node_InputInfo
const __meta_CostGraphDef_Node_InputInfo = Ref{ProtoMeta}()
function meta(::Type{CostGraphDef_Node_InputInfo})
    ProtoBuf.metalock() do
        if !isassigned(__meta_CostGraphDef_Node_InputInfo)
            __meta_CostGraphDef_Node_InputInfo[] = target = ProtoMeta(CostGraphDef_Node_InputInfo)
            allflds = Pair{Symbol,Union{Type,String}}[:preceding_node => Int32, :preceding_port => Int32]
            meta(target, CostGraphDef_Node_InputInfo, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_CostGraphDef_Node_InputInfo[]
    end
end
function Base.getproperty(obj::CostGraphDef_Node_InputInfo, name::Symbol)
    if name === :preceding_node
        return (obj.__protobuf_jl_internal_values[name])::Int32
    elseif name === :preceding_port
        return (obj.__protobuf_jl_internal_values[name])::Int32
    else
        getfield(obj, name)
    end
end

mutable struct CostGraphDef_Node_OutputInfo <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function CostGraphDef_Node_OutputInfo(; kwargs...)
        obj = new(meta(CostGraphDef_Node_OutputInfo), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct CostGraphDef_Node_OutputInfo
const __meta_CostGraphDef_Node_OutputInfo = Ref{ProtoMeta}()
function meta(::Type{CostGraphDef_Node_OutputInfo})
    ProtoBuf.metalock() do
        if !isassigned(__meta_CostGraphDef_Node_OutputInfo)
            __meta_CostGraphDef_Node_OutputInfo[] = target = ProtoMeta(CostGraphDef_Node_OutputInfo)
            allflds = Pair{Symbol,Union{Type,String}}[:size => Int64, :alias_input_port => Int64, :shape => TensorShapeProto, :dtype => Int32]
            meta(target, CostGraphDef_Node_OutputInfo, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_CostGraphDef_Node_OutputInfo[]
    end
end
function Base.getproperty(obj::CostGraphDef_Node_OutputInfo, name::Symbol)
    if name === :size
        return (obj.__protobuf_jl_internal_values[name])::Int64
    elseif name === :alias_input_port
        return (obj.__protobuf_jl_internal_values[name])::Int64
    elseif name === :shape
        return (obj.__protobuf_jl_internal_values[name])::TensorShapeProto
    elseif name === :dtype
        return (obj.__protobuf_jl_internal_values[name])::Int32
    else
        getfield(obj, name)
    end
end

mutable struct CostGraphDef_Node <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function CostGraphDef_Node(; kwargs...)
        obj = new(meta(CostGraphDef_Node), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct CostGraphDef_Node
const __meta_CostGraphDef_Node = Ref{ProtoMeta}()
function meta(::Type{CostGraphDef_Node})
    ProtoBuf.metalock() do
        if !isassigned(__meta_CostGraphDef_Node)
            __meta_CostGraphDef_Node[] = target = ProtoMeta(CostGraphDef_Node)
            fnum = Int[1,2,3,4,5,6,12,10,11,16,9,14,15,7,8,17]
            pack = Symbol[:control_input]
            allflds = Pair{Symbol,Union{Type,String}}[:name => AbstractString, :device => AbstractString, :id => Int32, :input_info => Base.Vector{CostGraphDef_Node_InputInfo}, :output_info => Base.Vector{CostGraphDef_Node_OutputInfo}, :temporary_memory_size => Int64, :persistent_memory_size => Int64, :host_temp_memory_size => Int64, :device_temp_memory_size => Int64, :device_persistent_memory_size => Int64, :compute_cost => Int64, :compute_time => Int64, :memory_time => Int64, :is_final => Bool, :control_input => Base.Vector{Int32}, :inaccurate => Bool]
            meta(target, CostGraphDef_Node, allflds, ProtoBuf.DEF_REQ, fnum, ProtoBuf.DEF_VAL, pack, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_CostGraphDef_Node[]
    end
end
function Base.getproperty(obj::CostGraphDef_Node, name::Symbol)
    if name === :name
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :device
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :id
        return (obj.__protobuf_jl_internal_values[name])::Int32
    elseif name === :input_info
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{CostGraphDef_Node_InputInfo}
    elseif name === :output_info
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{CostGraphDef_Node_OutputInfo}
    elseif name === :temporary_memory_size
        return (obj.__protobuf_jl_internal_values[name])::Int64
    elseif name === :persistent_memory_size
        return (obj.__protobuf_jl_internal_values[name])::Int64
    elseif name === :host_temp_memory_size
        return (obj.__protobuf_jl_internal_values[name])::Int64
    elseif name === :device_temp_memory_size
        return (obj.__protobuf_jl_internal_values[name])::Int64
    elseif name === :device_persistent_memory_size
        return (obj.__protobuf_jl_internal_values[name])::Int64
    elseif name === :compute_cost
        return (obj.__protobuf_jl_internal_values[name])::Int64
    elseif name === :compute_time
        return (obj.__protobuf_jl_internal_values[name])::Int64
    elseif name === :memory_time
        return (obj.__protobuf_jl_internal_values[name])::Int64
    elseif name === :is_final
        return (obj.__protobuf_jl_internal_values[name])::Bool
    elseif name === :control_input
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{Int32}
    elseif name === :inaccurate
        return (obj.__protobuf_jl_internal_values[name])::Bool
    else
        getfield(obj, name)
    end
end

mutable struct CostGraphDef_AggregatedCost <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function CostGraphDef_AggregatedCost(; kwargs...)
        obj = new(meta(CostGraphDef_AggregatedCost), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct CostGraphDef_AggregatedCost
const __meta_CostGraphDef_AggregatedCost = Ref{ProtoMeta}()
function meta(::Type{CostGraphDef_AggregatedCost})
    ProtoBuf.metalock() do
        if !isassigned(__meta_CostGraphDef_AggregatedCost)
            __meta_CostGraphDef_AggregatedCost[] = target = ProtoMeta(CostGraphDef_AggregatedCost)
            allflds = Pair{Symbol,Union{Type,String}}[:cost => Float32, :dimension => AbstractString]
            meta(target, CostGraphDef_AggregatedCost, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_CostGraphDef_AggregatedCost[]
    end
end
function Base.getproperty(obj::CostGraphDef_AggregatedCost, name::Symbol)
    if name === :cost
        return (obj.__protobuf_jl_internal_values[name])::Float32
    elseif name === :dimension
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    else
        getfield(obj, name)
    end
end

mutable struct CostGraphDef <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function CostGraphDef(; kwargs...)
        obj = new(meta(CostGraphDef), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct CostGraphDef
const __meta_CostGraphDef = Ref{ProtoMeta}()
function meta(::Type{CostGraphDef})
    ProtoBuf.metalock() do
        if !isassigned(__meta_CostGraphDef)
            __meta_CostGraphDef[] = target = ProtoMeta(CostGraphDef)
            allflds = Pair{Symbol,Union{Type,String}}[:node => Base.Vector{CostGraphDef_Node}, :cost => Base.Vector{CostGraphDef_AggregatedCost}]
            meta(target, CostGraphDef, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_CostGraphDef[]
    end
end
function Base.getproperty(obj::CostGraphDef, name::Symbol)
    if name === :node
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{CostGraphDef_Node}
    elseif name === :cost
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{CostGraphDef_AggregatedCost}
    else
        getfield(obj, name)
    end
end

export CostGraphDef_Node_InputInfo, CostGraphDef_Node_OutputInfo, CostGraphDef_Node, CostGraphDef_AggregatedCost, CostGraphDef


text:
using LinearAlgebra, StatsBase, TimerOutputs, DataFrames, StatsPlots, BenchmarkTools

function structured_P(L::Int; p::Float64 = 0.45, r::Float64 = 0.01)::Matrix{Float64}
    q = 1 - p - r
    P = diagm(fill(r,L)) + diagm(-1=>fill(q,L-1)) + diagm(1 => fill(p,L-1))
    P[1,1] = 1-p
    P[L,L] = 1-q
    return P
end

structured_π(L::Int; p::Float64 = 0.45, r::Float64 = 0.01)::Vector{Float64} = begin
    q = 1 - p - r
    [(p/q)^i  for i in 1:L] * (q-p) / p / (1-(p/q)^L) #Explicit expression (birth death)
end;

# Method 1: Linear Algebra
function method_1(input_matrix, L)
    # Transpose both sides, add a row to the bottom to represent the additional requirement
    return [I - input_matrix'; ones(1, L)] \ [zeros(L); 1]
end

# Method 2: Take a high matrix power
function method_2(input_matrix, L)
    matrix_power = 10^5 # arbitrary
    # Return the first column - the 'j' doesn't matter
    return (input_matrix^(matrix_power))[1, :]
end

# Method 3: Eigen trickery
function method_3(input_matrix, L)
    eigens = eigen(input_matrix')
    correct_vector = eigens.vectors[:, findfirst(isapprox(1), eigens.values)]
    return correct_vector * UniformScaling(1 / sum(correct_vector))
end

# Method 4: Random sample
function method_4(input_matrix, L)
    repeats = 10^5 # arbitrary
    states = 1:L
    current_value = rand(states)
    chain = zeros(L) # we don't need the order, just the counts
    for i in 1:repeats
        # Take a weighted sample based on the current state
        current_value = sample(states, Weights(view(input_matrix, current_value, :)))
        chain[current_value] += 1
    end
    # Apply the 1/n
    return map(x->x/repeats, chain)
end

# The L's given by the task
L_given = [2:5..., 10:10:50..., 100:100:500..., 1000]

# Get the benchmarks warmed up and ready to go
suite = BenchmarkGroup()
for method in [method_1, method_2, method_3, method_4]
    suite[String(Symbol(method))] = BenchmarkGroup()
end

# Prepare data structures
norms_df = DataFrame(a=L_given)
approx_methods_for_df = zip('b':'e', [method_1, method_2, method_3, method_4])

for approx_method ∈ approx_methods_for_df
    # Benchmark each L for each function, add it to suite
    for L in L_given
        suite[String(Symbol(approx_method[2]))][L] = @benchmarkable $(approx_method[2])(structured_P($L), $L)
    end

    # Put the error into the df for plotting
    norms_df[!, string(approx_method[1])] = map(L->norm(approx_method[2](structured_P(L), L) - structured_π(L)), L_given)
end

# Run our thing
results = run(suite, seconds = 2, samples = 5)

# Plot the DataFrames

# Make sure none of the values break the graph
inds = (norms_df.c .> 0) .& (norms_df.d .> 0) .& (norms_df.b .> 0) .& (norms_df.e .> 0)
# Plot the absolute errors of each method
@show @df norms_df plot(
    :a[inds], [:b[inds] :c[inds] :d[inds] :e[inds]],
    xaxis=:log10,
    xtickfont = font(5, "Courier"),
    xlabel = "L",
    ylabel = "Absolute Error",
    yscale=:log10,
    label = ["Method 1 (Linear Algebra)" "Method 2 (Matrix Power)" "Method 3 (Eigenvector)" "Method 4 (Random Seq)"],
    title = "Demonstration of Approximation Methods"
)

# Prepare the benchmarks for plotting (load their median times into a dataframe)
benchmarks_df = DataFrame(L=L_given)
for method in keys(results)
    benchmarks_df[!, String(Symbol(method))] = map(L->median(results[method][L]).time, L_given)
end

# Plot the benchmarks
@show @df benchmarks_df plot(
    :L, [:method_1 :method_2 :method_3 :method_4],
    xaxis=:log10,
    xtickfont = font(5, "Courier"),
    xlabel = "L",
    ylabel = "Time (in one unit or another)",
    yscale=:log10,
    label = ["Method 1 (Linear Algebra)" "Method 2 (Matrix Power)" "Method 3 (Eigenvector)" "Method 4 (Random Seq)"],
    title = "Demonstration of Approximation Methods"
)


text:
using PGFPlots, GridInterpolations, Colors, ColorBrewer
using AutomotivePOMDPs
using AutomotiveDrivingModels
using POMDPPolicies
using SARSOP
using LinearAlgebra
using POMDPs
using POMDPModelTools
using POMDPPolicies
using ProgressMeter

function get_XY_grid(pomdp::SingleOCPOMDP)
    env = pomdp.env
    rl = env.params.roadway_length
    cw = env.params.crosswalk_width
    lw = env.params.lane_width
    cl = env.params.crosswalk_length
    x_ped = rl/2
    Y = LinRange(-cl/4, cl/4, Int(floor(cl/2/pomdp.pos_res)) + 1)
    X = LinRange(5., rl/2 + 2*cw , Int(floor((rl/2 + 2*cw -5)/pomdp.pos_res)) + 1)
    return X, Y
end

function get_grid_data(policy::AlphaVectorPolicy, grid::RectangleGrid, pomdp::SingleOCPOMDP, v_ego::Float64, v_ped::Float64)
    values = zeros(n_actions(pomdp),length(grid))
    for i=1:length(grid)
        x, y = ind2x(grid, i)
        ego = AutomotivePOMDPs.xv_to_state(pomdp, x, v_ego)
        ped = AutomotivePOMDPs.yv_to_state(pomdp, y, v_ped)
        s = SingleOCState(collision_checker(ego, ped, pomdp.ego_type, pomdp.ped_type), ego, ped)
        # si = stateindex(pomdp, s)
        values[:,i] =  actionvalues(policy, Deterministic(s))
    end
    return values
end

function get_belief(y::Float64, sig::Float64, n_bins::Int64)
    y_noise = LinRange(y-2*sig, y+2*sig, n_bins)
    weights = zeros(n_bins)
    for i=1:n_bins
        weights[i] = 1/(sig*sqrt(2*pi))*exp(-(y_noise[i] - y)^2/(2*sig^2))
    end
    normalize!(weights, 1)
    return y_noise, weights
end

function get_belief(pomdp::SingleOCPOMDP, x::Float64, y::Float64, sig::Float64, n_bins::Int64, v_ego::Float64, v_ped::Float64)
    y_noise = LinRange(y-2*sig, y+2*sig, n_bins)
    belief = SparseCat(Vector{SingleOCState}(undef, n_bins), Vector{Float64}(undef, n_bins))
    for i=1:n_bins
        ego = AutomotivePOMDPs.xv_to_state(pomdp, x, v_ego)
        ped = AutomotivePOMDPs.yv_to_state(pomdp, y_noise[i], v_ped)
        s = SingleOCState(collision_checker(ego, ped, pomdp.ego_type, pomdp.ped_type), ego, ped)
        weight = 1/(sig*sqrt(2*pi))*exp(-(y_noise[i] - y)^2/(2*sig^2))
        belief.vals[i] = s
        belief.probs[i] = weight
    end
    normalize!(belief.probs, 1)
    return belief
end

function get_deterministic_belief(pomdp::SingleOCPOMDP, x::Float64, y::Float64, sig::Float64, n_bins::Int64, v_ego::Float64, v_ped::Float64)
    ego = AutomotivePOMDPs.xv_to_state(pomdp, x, v_ego)
    ped = AutomotivePOMDPs.yv_to_state(pomdp, y, v_ped)
    s = SingleOCState(collision_checker(ego, ped, pomdp.ego_type, pomdp.ped_type), ego, ped)
    return Deterministic(s)
end

function interpolate_belief(pomdp::SingleOCPOMDP, x::Float64, y::Float64, sig::Float64, n_bins::Int64, v_ego::Float64, v_ped::Float64)
    y_noise = LinRange(y-2*sig, y+2sig, n_bins)
    vals = SingleOCState[]
    probs = Float64[]
    X, Y = get_XY_grid(pomdp::SingleOCPOMDP)
    X = X[1:21]
    grid = RectangleGrid(X, Y)
    for i=1:n_bins
        itp_xy, itp_weights = interpolants(grid, [x, y_noise[i]])
        for (j, xy_ind) in enumerate(itp_xy)
            xg, yg = ind2x(grid, xy_ind)
            ego = AutomotivePOMDPs.xv_to_state(pomdp, xg, v_ego)
            ped = AutomotivePOMDPs.yv_to_state(pomdp, yg, v_ped)
            s = SingleOCState(collision_checker(ego, ped, pomdp.ego_type, pomdp.ped_type), ego, ped)
            weight = itp_weights[j] * 1/(sig*sqrt(2*pi))*exp(-(y_noise[i] - y)^2/(2*sig^2))
            if !(s in vals)
                push!(vals, s)
                push!(probs, weight)
            else
                si = findfirst(isequal(s), vals)
                probs[si] += weight
            end
        end
    end
    normalize!(probs, 1)
    return SparseCat(vals, probs)
end

function compute_acts(pomdp::SingleOCPOMDP, policy::AlphaVectorPolicy, sig::Float64 = 0.01, v_ego::Float64 = 5., v_ped::Float64 = 1., n_bins::Int64 = 20, n_pts::Int64 = 100)
    X, Y = get_XY_grid(pomdp::SingleOCPOMDP)
    X = X[1:21]
    grid = RectangleGrid(X, Y)
    values = get_grid_data(policy, grid, pomdp, v_ego, v_ped)

    xrange = LinRange(X[1], X[end], n_pts)
    yrange = LinRange(Y[1], Y[end], n_pts)
    action_map = zeros(n_pts, n_pts)
    n_acts = n_actions(pomdp)

    for (i,x) in enumerate(xrange)
        for (j,y) in enumerate(yrange)
            y_noise, weights = get_belief(y, sig, n_bins)
            acts_xy = zeros(n_acts)
            for k=1:n_acts
                for l=1:n_bins
                    val = interpolate(grid, values[k,:], [x,y_noise[l]])
                    acts_xy[k] += weights[l]*val
                end
            end
            max, indmax = findmax(acts_xy)
            action_map[i,j] = policy.action_map[indmax].acc
        end
    end
    acts = action_map
    acts = acts'
    acts = acts[end:-1:1,1:end]
    return acts
end

function compute_acts2(pomdp::SingleOCPOMDP, policy::AlphaVectorPolicy, sig::Float64 = 0.01, v_ego::Float64 = 5., v_ped::Float64 = 1., n_bins::Int64 = 20, n_pts::Int64 = 100)
    X, Y = get_XY_grid(pomdp)
    X = X[1:21]
    grid = RectangleGrid(X, Y)

    xrange = LinRange(X[1], X[end], n_pts)
    yrange = LinRange(Y[1], Y[end], n_pts)
    action_map = zeros(n_pts, n_pts)
    n_acts = n_actions(pomdp)

    @showprogress for (i,x) in enumerate(xrange)
        for (j,y) in enumerate(yrange)
            belief = get_deterministic_belief(pomdp, x, y, sig, n_bins, v_ego, v_ped)
            # belief = interpolate_belief(pomdp, x, y, sig, n_bins, v_ego, v_ped)
            max_val, ind_max = findmax(itp_values(policy, belief))
            action_map[i,j] = policy.action_map[ind_max].acc
            # action_map[i, j] = action(policy, belief).acc
        end
    end
    acts = action_map
    acts = acts'
    acts = acts[end:-1:1,1:end]
    return acts
end

function itp_values(policy::AlphaVectorPolicy, b)
    num_vectors = length(policy.alphas)
    max_values = -Inf * ones(length(actions(policy.pomdp)))
    for i=1:num_vectors
        ai = actionindex(policy.pomdp, policy.action_map[i])
        temp_value = itp_value(policy.alphas[i], b, policy.pomdp)
        # for (j, s) in enumerate(b.vals)
        #     s_itps, w_itps = interpolate_state(s, policy.pomdp)
        #     for (k, s_itp) in enumerate(s_itps)
        #         si = stateindex(policy.pomdp, s_itp)
        #         temp_value += policy.alphas[i][si]*b.probs[j]*w_itps[k]
        #     end
        # end
        if temp_value > max_values[ai]
            max_values[ai] = temp_value
        end
    end
    return max_values
end

function itp_value(alpha::AbstractArray, b::SparseCat, pomdp::POMDP)
    val = 0.0
    for (j, s) in enumerate(b.vals)
        s_itps, w_itps = interpolate_state(s, pomdp)
        for (k, s_itp) in enumerate(s_itps)
            si = stateindex(pomdp, s_itp)
            val += alpha[si]*b.probs[j]*w_itps[k]
        end
    end
    return val
end

function itp_value(alpha::AbstractArray, b::Deterministic, pomdp::POMDP)
    val = 0.0
    s_itps, w_itps = interpolate_state(b.val, pomdp)
    for (k, s_itp) in enumerate(s_itps)
        si = stateindex(pomdp, s_itp)
        val += alpha[si]*w_itps[k]
    end
    return val
end

function compute_acts3(pomdp::SingleOCPOMDP, policy::AlphaVectorPolicy, sig::Float64 = 0.01, v_ego::Float64 = 5., v_ped::Float64 = 1., n_bins::Int64 = 20, n_pts::Int64 = 100)
    X, Y = get_XY_grid(pomdp)
    X = X[1:21]
    action_map = zeros(length(X), length(Y))

    for (j, y) in enumerate(Y)
        for (i, x) in enumerate(X)
            ego = AutomotivePOMDPs.xv_to_state(pomdp, x, v_ego)
            ped = AutomotivePOMDPs.yv_to_state(pomdp, y, v_ped)
            s = SingleOCState(collision_checker(ego, ped, pomdp.ego_type, pomdp.ped_type), ego, ped)
            a = action(policy, Deterministic(s))
            action_map[i, j] = a.acc
        end
    end
    acts = action_map
    acts = acts'
    acts = acts[end:-1:1,1:end]
    return acts
end

function policy_plot(pomdp::SingleOCPOMDP, policy::Policy, ;
                     sig::Float64 = 0.01, v_ego::Float64 = 5., v_ped::Float64 = 1., n_bins::Int64 = 20, n_pts::Int64 = 100)
    acts = compute_acts2(pomdp, policy, sig, v_ego, v_ped, n_bins, n_pts)
    r = RGB{Float64}(1., 0., 0.) # red
    o = RGB{Float64}(1., 69/255.0, 0.) # orange
    y = RGB{Float64}(1., 1., 0.) # yellow
    g = RGB{Float64}(0., 1., 0.) # green
    colors = [r;o;y;g]
    cmap = ColorMaps.RGBArrayMap(colormap("RdBu"))
    X, Y = get_XY_grid(pomdp::SingleOCPOMDP)
    X = X[1:21]
    ax =Axis(Plots.Image(acts, (X[1], X[end]), (Y[1], Y[end]), colormap = cmap),
    xlabel = "Ego car position on the road",
    ylabel = "Pedestrian position along the crosswalk",
    title = "Policy Plot for a fixed ego velocity at $v_ego m/s and fixed pedestrian speed at $v_ped m/s")
    return ax
end

# function compute_acts(pomdp::SingleOCPOMDP, policy::NNPolicy, sig::Float64 = 0.01, v_ego::Float64 = 5., v_ped::Float64 = 1., n_bins::Int64 = 20, n_pts::Int64 = 100)
#     X, Y = get_XY_grid(pomdp::SingleOCPOMDP)
#
#     xrange = LinRange(X[1], X[end], n_pts)
#     yrange = LinRange(Y[1], Y[end], n_pts)
# action_map = zeros(n_pts, n_pts)
# n_acts = n_actions(pomdp)
#     for (i,x) in enumerate(xrange)
#         for (j,y) in enumerate(yrange)
#             acts_xy = zeros(1, n_acts)
#             o = normalize_o(x, v_ego, y, v_ped, pomdp)
#             acts_xy = policy.nn[:value](o)
#             max, indmax = findmax(acts_xy)
#             action_map[i,j] = policy.action_map[indmax].acc
#         end
#     end

max_stars_repo_name: MaximeBouton/AutomotivePOMDPs

max_issues_repo_name: MaximeBouton/AutomotivePOMDPs

max_forks_repo_name: MaximeBouton/AutomotivePOMDPs


text:
using StackOverflow
using HTTP

# the function 'geterrors' returns a 2d array containing the error_names, error ids and description in that order.

function geterrors()
    error_id = []
    description = []
    error_name = []

    r = HTTP.request("GET", "https://api.stackexchange.com/2.2/errors")

    json_obj = StackOverflow.convert_HTTP_Response_To_JSON(r)

    for (k, v) in json_obj
        if occursin("items", k)
            for item in v
                for (key, value) in item
                    if key == "error_id"
                        push!(error_id, value)
                    elseif key == "description"
                        push!(description, value)
                    elseif key == "error_name"
                        push!(error_name, value)
                    end
                end

            end
            return [error_name, error_id, description]
        end
    end

end


text:
module MyModule

using JSON
export hello

function hello()
    x = Dict("status"=>200, "message"=>"hello world")
    return JSON.json(x)
end

end


text:
using Elefridge
using PyPlot
using PyCall
xr = pyimport("xarray")
using TranscodingStreams, CodecZstd
using JLD

ZstdCompressorL3 = ZstdCompressor(level=3)
TranscodingStreams.initialize(ZstdCompressorL3)

ZstdCompressorL10 = ZstdCompressor(level=10)
TranscodingStreams.initialize(ZstdCompressorL10)

ZstdCompressorL22 = ZstdCompressor(level=22)
TranscodingStreams.initialize(ZstdCompressorL22)

path = "/Users/milan/cams/gridded"
filelist = filter(x->endswith(x,".grib"),readdir(path))

n = length(filelist)
varnames = fill("",n)
methods = ["ZstdL3","ZstdL3+bittranspose",
            "ZstdL10","ZstdL10+bittranspose",
            "ZstdL22","ZstdL22+bittranspose"]
sbits = Array(1:23)
C = fill(0.0,n,length(methods),length(sbits))

for (i,file) in enumerate(filelist)
    varname = split(split(file,"0_")[end],".")[1]
    varnames[i] = varname

    println("---")
    println("Reading $varname")
    X = xr.open_dataarray(joinpath(path,file),engine="cfgrib").data

    # transpose array
    X = permutedims(X,[3,2,1])  # unravels as longitude, latitude, level
    Xs = X[1:225,1:225,:]       # only 1/8th of the globe

    for (is,s) in enumerate(sbits)
        print(s)
        Xsr = round(Xs,s)
        Xsrt = bittranspose(Xsr)

        Xsr8 = unsafe_wrap(Array, Ptr{UInt8}(pointer(Xsr)), sizeof(Xsr))
        Xsrt8 = unsafe_wrap(Array, Ptr{UInt8}(pointer(Xsrt)), sizeof(Xsrt))

        C[i,1,is] = sizeof(Xs)/sizeof(transcode(ZstdCompressorL3,Xsr8))
        C[i,2,is] = sizeof(Xs)/sizeof(transcode(ZstdCompressorL3,Xsrt8))
        println("ZstdL3 done.")
        C[i,3,is] = sizeof(Xs)/sizeof(transcode(ZstdCompressorL10,Xsr8))
        C[i,4,is] = sizeof(Xs)/sizeof(transcode(ZstdCompressorL10,Xsrt8))
        println("ZstdL10 done.")
        C[i,5,is] = sizeof(Xs)/sizeof(transcode(ZstdCompressorL22,Xsr8))
        C[i,6,is] = sizeof(Xs)/sizeof(transcode(ZstdCompressorL22,Xsrt8))
        println("ZstdL22 done.")
    end
    @save "/Users/milan/cams/entropy/compression_transpose_all_zstd_1_8.jld" C varnames methods sbits
end


text:
module SymbolicManipulations

using SymbolicUtils
using Setfield
using DataStructures

export together, fullsimplify, integrate

include("utils.jl")
include("together.jl")
include("fullsimplify.jl")
include("integrate.jl")

end


text:
# TODO: Comment every single line

# ConicObj represents an affine function of the variables
# it is stored as a dictionary of (key, value) pairs
# keys are unique ids of variables
# values are their coefficients in the affine function
# so for example, {unique_id(x)=>5, unique_id(y)=>6} represents the function 5x + 6y
# we store the affine functions in this form for efficient manipulation of sparse affine functions
struct ConicObj
    mapping::OrderedDict{UInt64,Tuple{Value,Value}}
end
ConicObj() = ConicObj(OrderedDict{UInt64,Tuple{Value,Value}}())
Base.iterate(c::ConicObj, s...) = iterate(c.mapping, s...)
Base.keys(c::ConicObj) = keys(c.mapping)
Base.haskey(c::ConicObj, var::Integer) = haskey(c.mapping, UInt64(var))
Base.getindex(c::ConicObj, var::Integer) = c.mapping[UInt64(var)]
function Base.setindex!(c::ConicObj, val, var::Integer)
    return setindex!(c.mapping, val, UInt64(var))
end
Base.copy(c::ConicObj) = ConicObj(copy(c.mapping))

# helper function to negate conic objectives
# works by changing each (key, val) pair to (key, -val)
function Base.:-(c::ConicObj)
    new_obj = copy(c)
    for var in keys(new_obj)
        x1 = new_obj[var][1] * (-1)
        x2 = new_obj[var][2] * (-1)
        new_obj[var] = (x1, x2)
    end
    return new_obj
end

function Base.:+(c::ConicObj, d::ConicObj)
    new_obj = copy(c)
    for var in keys(d)
        if !haskey(new_obj, var)
            new_obj[var] = d[var]
        else
            # .+ does not preserve sparsity
            # need to override behavior
            if size(new_obj[var][1]) == size(d[var][1])
                x1 = new_obj[var][1] + d[var][1]
            else
                x1 = broadcast(+, new_obj[var][1], d[var][1])
            end
            if size(new_obj[var][2]) == size(d[var][2])
                x2 = new_obj[var][2] + d[var][2]
            else
                x2 = broadcast(+, new_obj[var][2], d[var][2])
            end
            new_obj[var] = (x1, x2)
        end
    end
    return new_obj
end

function get_row(c::ConicObj, row::Int)
    new_obj = ConicObj()
    for (var, coeff) in c
        x1 = coeff[1][row, :]
        x2 = coeff[2][row, :]
        new_obj[var] = (x1, x2)
    end
    return new_obj
end

function Base.:*(v::Value, c::ConicObj)
    # TODO: this part is time consuming, esp new_obj[var] = v * new_obj[var]...
    new_obj = copy(c)
    for var in keys(new_obj)
        x1 = v * new_obj[var][1]
        x2 = v * new_obj[var][2]
        new_obj[var] = (x1, x2)
    end
    return new_obj
end

function promote_size(c::ConicObj, vectorized_size::Int)
    new_obj = copy(c)
    for var in keys(new_obj)
        x1 = repeat(new_obj[var][1], vectorized_size, 1)
        x2 = repeat(new_obj[var][2], vectorized_size, 1)
        new_obj[var] = (x1, x2)
    end
    return new_obj
end

# A conic constraint is of the form [affine_expr1, affine_expr2, ..., affine_exprk] \in cone
# we represent each affine expressions as a ConicObj
# we represent the cone as a Symbol, like :SOC, :SDP, etc.
# See the details of `get_MOI_set` for the full list of symbols used.
# and we record the sizes of the affine expressions (XXX check...)
# XXX might it be better to represent objs as a single ConicObj rather than an array of them?
struct ConicConstr
    objs::Vector{ConicObj}
    cone::Symbol
    sizes::Vector{Int}
end

# in conic form, every expression e is represented by a ConicObj together with a collection of ConicConstrs
# for each expression e, UniqueExpMap maps (e.head, unique_id(e)) to that expression's ConicObj
const UniqueExpMap = OrderedDict{Tuple{Symbol,UInt64},ConicObj}
# for each expression e, UniqueExpMap maps (e.head, unique_id(e)) to the index of expression's ConicConstr in UniqueConstrList
const UniqueConstrMap = OrderedDict{Tuple{Symbol,UInt64},Int}
# records each ConicConstr created
const UniqueConstrList = Vector{ConicConstr}
# map variables' hash to the variable itself
const IdToVariables = OrderedDict{UInt64,AbstractVariable}
const ConicConstrToConstr = OrderedDict{ConicConstr,Constraint}
# UniqueConicForms caches all the conic forms of expressions we've parsed so far
struct UniqueConicForms
    exp_map::UniqueExpMap
    constr_map::UniqueConstrMap
    constr_list::UniqueConstrList
    id_to_variables::IdToVariables
    conic_constr_to_constr::ConicConstrToConstr
end

function UniqueConicForms()
    return UniqueConicForms(
        UniqueExpMap(),
        UniqueConstrMap(),
        ConicConstr[],
        IdToVariables(),
        ConicConstrToConstr(),
    )
end

function has_conic_form(conic_forms::UniqueConicForms, exp::AbstractExpr)
    return haskey(conic_forms.exp_map, (exp.head, exp.id_hash))
end

function has_conic_form(conic_forms::UniqueConicForms, constr::Constraint)
    return haskey(conic_forms.constr_map, (constr.head, constr.id_hash))
end

function get_conic_form(conic_forms::UniqueConicForms, exp::AbstractExpr)
    return conic_forms.exp_map[(exp.head, exp.id_hash)]
end

function get_conic_form(conic_forms::UniqueConicForms, constr::Constraint)
    return conic_forms.constr_map[(constr.head, constr.id_hash)]
end

function cache_conic_form!(
    conic_forms::UniqueConicForms,
    exp::AbstractExpr,
    new_conic_form::ConicObj,
)
    return conic_forms.exp_map[(exp.head, exp.id_hash)] = new_conic_form
end

function cache_conic_form!(
    conic_forms::UniqueConicForms,
    constr::Constraint,
    new_conic_form::ConicConstr,
)
    conic_forms.constr_map[(constr.head, constr.id_hash)] = 0
    return push!(conic_forms.constr_list, new_conic_form)
end

function cache_conic_form!(
    conic_forms::UniqueConicForms,
    constr::Constraint,
    new_conic_forms::UniqueConstrList,
)
    conic_forms.constr_map[(constr.head, constr.id_hash)] = 0
    return append!(conic_forms.constr_list, new_conic_forms)
end

function add_to_id_to_variables!(
    conic_forms::UniqueConicForms,
    var::AbstractVariable,
)
    return conic_forms.id_to_variables[var.id_hash] = var
end


text:
mutable struct Sound <: AbstractSound{sfSound}
    ptr::Ptr{sfSound}
    buffer::Maybe{SoundBuffer}

    function Sound()
        ptr = sfSound_create()
        ptr == C_NULL && error("Couldn't create Sound object.")
        this = new(ptr, nothing)
        finalizer(this) do obj
            try
                !isnothing(obj.buffer) && destroy!(obj.buffer)
                destroy!(obj)
            finally
                return
            end
        end
        return this
    end

    function Sound(file::AbstractString)
        this = Sound()
        setbuffer!(this, SoundBuffer(file))
        return this
    end
end

function getbuffer(sound::Sound)
    if !isnothing(sound.buffer)
        return sound.buffer
    else
        error("sound object doesn't have a SoundBuffer")
    end
end

destroy!(sound::Sound) = sfSound_destroy(sound)
"""
        setbuffer!(sound, buffer::SoundBuffer) -> Nothing

Set the source buffer containing the audio data to play.
"""
function setbuffer!(sound::Sound, buffer::SoundBuffer)
    sound.buffer = buffer
    sfSound_setBuffer(sound, buffer)
end

loop!(sound::Sound, loop::Bool) = sfSound_setLoop(sound, loop ? sfTrue : sfFalse)
islooping(sound::Sound) = sfSound_getLoop(sound) == sfTrue

timepos!(sound::Sound, offset::Integer) = sfSound_setPlayingOffset(sound, sfTime(offset))
timepos(sound::Sound) = sfSound_getPlayingOffset(sound).microseconds
samplepos(sound::Sound) = timeToSamples(timepos(sound), samplerate(sound), nchannels(sound))

volume!(sound::Sound, volume::Real) = sfSound_setVolume(sound, volume)
volume(sound::Sound) = sfSound_getVolume(sound)

pitch(sound::Sound) = sfSound_getPitch(sound)
pitch!(sound::Sound, pitch::Real) = sfSound_setPitch(sound, pitch)

getstatus(sound::Sound) = Status[ sfSound_getStatus(sound) ]

play!(sound::Sound)  = sfSound_play(sound)
stop!(sound::Sound)  = sfSound_stop(sound)
pause!(sound::Sound) = sfSound_pause(sound)

nchannels(this::Sound)    = nchannels(getbuffer(this))
duration(this::Sound)     = duration(getbuffer(this))
getsamples(this::Sound)   = getsamples(getbuffer(this))
nsamples(this::Sound)     = nsamples(getbuffer(this))
samplerate(this::Sound)   = samplerate( getbuffer(this))

isplaying(sound::Sound) = getstatus(sound) == :playing
ispaused(sound::Sound)  = getstatus(sound) == :paused
isstopped(sound::Sound) = getstatus(sound) == :stopped

getprops(sound::Sound) = (
    loop = islooping(sound),
    status = getstatus(sound),
    timepos = timepos(sound),
    volume = volume(sound),
    pitch = pitch(sound),
    channels = nchannels(sound),
    duration = duration(sound),
    nsamples = nsamples(sound),
    samplerate = samplerate(sound),
)

function Base.show(io::IO, sound::Sound)
    d = round(1e-6 * duration(sound), digits = 2)
    n = nchannels(sound)
    sr = round(1e-3 * samplerate(sound), digits = 1)
    print(io, "Sound ♫ $(d)s ($(sr)kHz, $(n)-channels)")
end

Base.getindex(sound::Sound, idx...) = getindex(getbuffer(sound), idx...)

Base.size(sound::Sound) = size(getbuffer(sound))
Base.length(sound::Sound) = length(getbuffer(sound))

Base.firstindex(sound::Sound) = firstindex(getbuffer(sound))
Base.lastindex(sound::Sound) = lastindex(getbuffer(sound))

Base.first(sound::Sound) = first(getbuffer(sound))
Base.last(sound::Sound)  = last(getbuffer(sound))

FFTW.fft(sound::Sound) = fft(getbuffer(sound)[:])
FFTW.ifft(sound::Sound) = ifft(getbuffer(sound)[:])


text:
# Breadth First Search Example

using ThreeJS
using Escher
using LightGraphs
using NetworkViz

main(window) = begin
    push!(window.assets, "widgets")
    push!(window.assets,("ThreeJS","threejs"))
    nodes = Signal(10)
    btn = Signal(Any, leftbutton)
    start = Signal(false)
    map(btn) do _
      push!(start,!start.value)
    end
    map(start,nodes,sampleon(btn,nodes)) do t,n1,n2
    vbox(
        vbox(
          h1("Breadth-First Traversal"),
          vskip(2em),
          "Number of Nodes",
          slider(10:100,disabled=!t) >>> nodes,
          vskip(2em),
          hbox(
            hskip(4em),
            button("Draw",disabled=!t) >>> btn
          )
        ) |> pad(2em),
        hbox(
          hskip(2em),
          vbox(
            if t
              vbox(
                md"""## Original Graph""",
                drawGraph(CompleteGraph(n1),z=1)
              )
            else
              vbox(
                md"""## BFS Traversal of the Graph""",
                drawGraph(bfs_tree(CompleteGraph(n2),1),z=1)
              )
            end
          )
        )
      )
  end
end


text:
using ExtremeStats
using Distributions
using Plots
using VisualRegressionTests
using Test, Pkg, Random

# workaround GR warnings
ENV["GKSwstype"] = "100"

# environment settings
islinux = Sys.islinux()
istravis = "TRAVIS" ∈ keys(ENV)
datadir = joinpath(@__DIR__,"data")
visualtests = !istravis || (istravis && islinux)
if !istravis
  Pkg.add("Gtk")
  using Gtk
end

# list of tests
testfiles = [
  "maxima.jl",
  "fitting.jl",
  "plotrecipes.jl"
]

@testset "ExtremeStats.jl" begin
  for testfile in testfiles
    include(testfile)
  end
end