text:
"""
    Quaternions

A package for quaternion arithmetic.

"""
module Quaternions
import Base: convert, promote_rule, show, real, imag, conj, abs, abs2, inv, rand, randn
import Base: +, -, /, *, &, ⊻, |
import Base: inv, float, isreal, isinteger, isfinite, isnan, isinf, iszero, isone
import Base: widen, big, round

import Random
import LinearAlgebra: pinv, norm

export Quaternion, QuaternionF64, QuaternionF32, QuaternionF16
export quaternion, imagi, imagj, imagk, jm, km, ispure
export wedge, antiwedge, perp, para

"""
    Quaternion{T<:Real} <: Number

Quaternion number type with real and imaginary parts of type `T`.

`QuaternionF16`, `QuaternionF32` and `QuaternionF64` are aliases for
`Quaternion{Float16}`, `Quaternion{Float32}` and `Quaternion{Float64}` respectively.
"""
struct Quaternion{T<:Real} <: Number
    q0::T
    q1::T
    q2::T
    q3::T
end

include("qmath.jl")

Quaternion(q0::Real,q1::Real,q2::Real,q3::Real) = Quaternion(promote(q0,q1,q2,q3)...)
Quaternion(x::Real) = Quaternion(x,zero(x),zero(x),zero(x))
Quaternion(x::Irrational) = Quaternion(float(x))
Quaternion(z::Complex) = Quaternion(real(z),imag(z),zero(real(z)),zero(real(z)))

const QuaternionF64 = Quaternion{Float64}
const QuaternionF32 = Quaternion{Float32}
const QuaternionF16 = Quaternion{Float16}

convert(::Type{Quaternion}, x::Real) = Quaternion(x)
convert(::Type{Quaternion}, z::Complex) = Quaternion(z)
convert(::Type{Quaternion{T}}, x::Real) where T<:Real = Quaternion(convert(T,x))
convert(::Type{Quaternion{T}}, z::Complex) where T<:Real = Quaternion(convert(Complex{T},z))
convert(::Type{Quaternion{T}}, q::Quaternion) where T<:Real  = Quaternion{T}(convert(T,q.q0), convert(T,q.q1), convert(T,q.q2), convert(T,q.q3))
convert(::Type{Quaternion{T}}, q::Quaternion{T}) where T<:Real = q
convert(::Type{T}, q::Quaternion) where T<:Real = (iszero(q.q1) && iszero(q.q2) && iszero(q.q3)) ? convert(T,q.q0) : throw(InexactError())

promote_rule(::Type{Quaternion{T}}, ::Type{T}) where T<:Real = Quaternion{T}
promote_rule(::Type{Quaternion}, ::Type{T}) where T<:Real = Quaternion
promote_rule(::Type{Quaternion{T}}, ::Type{Complex{T}}) where T<:Real = Quaternion{T}
promote_rule(::Type{Quaternion}, ::Type{Complex{T}}) where T<:Real = Quaternion
promote_rule(::Type{Quaternion{T}}, ::Type{S}) where {T<:Real, S<:Real} = Quaternion{promote_type(T,S)}
promote_rule(::Type{Quaternion{T}}, ::Type{Complex{S}}) where {T<:Real, S<:Real} = Quaternion{promote_type(T,S)}
promote_rule(::Type{Quaternion{T}}, ::Type{Quaternion{S}}) where {T<:Real, S<:Real} = Quaternion{promote_type(T,S)}

quaternion(q0,q1,q2,q3) = Quaternion(q0,q1,q2,q3)
quaternion(x) = Quaternion(x)
quaternion(z::Complex) = Quaternion(z)
quaternion(q::Quaternion) = q

widen(::Type{Quaternion{T}}) where {T} = Quaternion{widen(T)}
float(::Type{Quaternion{T}}) where {T<:AbstractFloat} = Quaternion{T}
float(::Type{Quaternion{T}}) where {T} = Quaternion{float(T)}
float(z::Quaternion{<:AbstractFloat}) = z
float(z::Quaternion) = Quaternion(float(real(z)), float(imagi(z)), float(imagj(z)), float(imagk(z)))
big(::Type{Quaternion{T}}) where {T<:Real} = Quaternion{big(T)}
big(z::Quaternion{T}) where {T<:Real} = Quaternion{big(T)}(z)

function show(io::IO, q::Quaternion)
    compact = get(io, :compact, false)
    show(io, q.q0)
    p(i, unit) = begin
        if signbit(i) && !isnan(i)
            print(io, compact ? "-" : " - ")
            if isa(i,Signed) && !isa(i,BigInt) && i == typemin(typeof(i))
                show(io, -widen(i))
            else
                show(io, -i)
            end
        else
            print(io, compact ? "+" : " + ")
            show(io, i)
        end
        if !(isa(i,Integer) && !isa(i,Bool) || isa(i,AbstractFloat) && isfinite(i))
            print(io, "*")
        end
        print(io, unit)
    end
    p(q.q1, "im")
    p(q.q2, "jm")
    p(q.q3, "km")
end

show(io::IO, q::Quaternion{Bool}) =
    print(io, q == im ? "im" :
              q == jm ? "jm" :
              q == km ? "km" :
          "Quaternion($(q.q0),$(q.q1),$(q.q2),$(q.q3))")

function quaternion(A::Array{S}, B::Array{T}, C::Array{U}, D::Array{V}) where {S<:Real, T<:Real, U<:Real, V<:Real}
    if !(size(A)==size(B)==size(C)==size(D)); error("argument dimensions must match"); end
    F = similar(A, typeof(quaternion(zero(S),zero(T),zero(U),zero(V))))
    for i=1:length(A)
        @inbounds F[i] = quaternion(A[i], B[i], C[i], D[i])
    end
    return F
end

quaternion(x::AbstractArray{T}) where T<:Quaternion = x
function quaternion(A::AbstractArray{T}) where T
    if !isconcretetype(T)
        error("`quaternion` not defined on abstractly-typed arrays; please convert to a more specific type")
    end
    convert(AbstractArray{typeof(quaternion(zero(T)))}, A)
end

real(z::Quaternion) = z.q0
imagi(z::Quaternion) = z.q1
imagj(z::Quaternion) = z.q2
imagk(z::Quaternion) = z.q3
imagi(z::Real) = zero(z)
imagj(z::Real) = zero(z)
imagk(z::Real) = zero(z)
imag(z::Quaternion) = Quaternion(zero(z.q0),z.q1,z.q2,z.q3)

imagi(x::AbstractVector{T}) where T<:Real = zero(x)
imagj(x::AbstractVector{T}) where T<:Real = zero(x)
imagk(x::AbstractVector{T}) where T<:Real = zero(x)

ispure(z::Quaternion)::Bool = iszero(z.q0)

for fn in (:imagi,:imagj,:imagk)
    @eval begin
        ($fn)(A::AbstractArray) = map(($fn),A)
    end
end

conj(z::Quaternion) = Quaternion(z.q0, -z.q1, -z.q2, -z.q3)
abs(z::Quaternion) = sqrt(z.q0*z.q0 + z.q1*z.q1 + z.q2*z.q2 + z.q3*z.q3)
abs2(z::Quaternion) = z.q0*z.q0 + z.q1*z.q1 + z.q2*z.q2 + z.q3*z.q3
inv(z::Quaternion) = conj(z)/abs2(z)

(+)(z::Quaternion) = Quaternion(+z.q0, +z.q1, +z.q2, +z.q3)
(-)(z::Quaternion) = Quaternion(-z.q0, -z.q1, -z.q2, -z.q3)
(/)(z::Quaternion, x::Real) = Quaternion(z.q0/x, z.q1/x, z.q2/x, z.q3/x)
(+)(z::Quaternion, w::Quaternion) = Quaternion(z.q0 + w.q0, z.q1 + w.q1,
                                                           z.q2 + w.q2, z.q3 + w.q3)
(-)(z::Quaternion, w::Quaternion) = Quaternion(z.q0 - w.q0, z.q1 - w.q1,
                                                           z.q2 - w.q2, z.q3 - w.q3)
(*)(z::Quaternion, w::Quaternion) = Quaternion(z.q0*w.q0 - z.q1*w.q1 - z.q2*w.q2 - z.q3*w.q3,
                                                           z.q0*w.q1 + z.q1*w.q0 + z.q2*w.q3 - z.q3*w.q2,
                                                           z.q0*w.q2 - z.q1*w.q3 + z.q2*w.q0 + z.q3*w.q1,
                                                           z.q0*w.q3 + z.q1*w.q2 - z.q2*w.q1 + z.q3*w.q0)
(/)(z::Quaternion, w::Quaternion) = z*inv(w)

# element wise multiplication
mul_ew(z1::Quaternion, z2::Quaternion)::Quaternion = Quaternion(z1.q0*z2.q0,z1.q1*z2.q1,z1.q2*z2.q2,z1.q3*z2.q3)
(&)(z1::Quaternion, z2::Quaternion)::Quaternion = mul_ew(z1,z2)
(&)(z::Quaternion, x::Real) = x*z
(&)(z::Real, x::Quaternion) = x*z

wedge(p::Quaternion,q::Quaternion) = (p*q-q*p)/2
antiwedge(p::Quaternion,q::Quaternion) = (p*q+q*p)/2
(⊻)(z1::Quaternion, z2::Quaternion) = wedge(z1,z2)
(|)(z1::Quaternion, z2::Quaternion) = antiwedge(z1,z2)

para(p, q::Quaternion) = (p - q*p*q)/2
perp(p, q::Quaternion) = (p + q*p*q)/2

rand(r::Random.AbstractRNG, ::Random.SamplerType{Quaternion{T}}) where {T<:Integer} = Quaternion{T}(rand(r, T), rand(r, T), rand(r,T), rand(r,T))
#randn(r::Random.AbstractRNG, ::Type{Quaternion{T}}) where {T<:AbstractFloat} = Quaternion(0.5*randn(r, T), 0.5*randn(r, T), 0.5*randn(r,T), 0.5*randn(r,T))
function Base.rand(r::Random.AbstractRNG, ::Random.SamplerType{Quaternion{T}}) where {T <: AbstractFloat}
    # http://planning.cs.uiuc.edu/node198.html
    # https://github.com/KieranWynn/pyquaternion/blob/master/pyquaternion/quaternion.py#L260
    tau = 2.0 * pi
    r1 = rand(r, T)
    r2 = rand(r, T)
    r3 = rand(r, T)
    q1 = sqrt(one(T) - r1) * (sin(tau * r2))
    q2 = sqrt(one(T) - r1) * (cos(tau * r2))
    q3 = sqrt(r1)       * (sin(tau * r3))
    q4 = sqrt(r1)       * (cos(tau * r3))
    Quaternion{T}(q1, q2, q3, q4)
end

function Base.randn(r::Random.AbstractRNG, ::Type{Quaternion{T}}) where {T <: AbstractFloat}
    Quaternion{T}(randn(r, T), randn(r, T), randn(r, T), randn(r, T))*T(0.5)
end

real(::Type{Quaternion{T}}) where T<:Real = T
quaternion(::Type{T}) where T<:Real = Quaternion{T}
quaternion(::Type{Quaternion{T}}) where T<:Real = Quaternion{T}

"""
    jm

The quaternion imaginary unit `j`.

# Examples
```jldoctest
julia> jm * jm
-1 + 0im + 0jm + 0km
```
"""
const jm = Quaternion(false,false,true,false)

"""
    km

The quaternion imaginary unit `k`.

# Examples
```jldoctest
julia> km * km
-1 + 0im + 0jm + 0km
```
"""
const km = Quaternion(false,false,false,true)

function equiv!(Y::AbstractArray{Complex{T}}, X::AbstractArray{Quaternion{T}}) where T<:Real
    @assert size(X,1)<<1 == size(Y,1)
    @assert size(X,2)<<1 == size(Y,2)

    dtype = real(eltype(X))
    M = size(X,1)
    N = size(X,2)
    Xr = reshape(reinterpret(dtype, X), (M<<2,N))
    Yr = reshape(reinterpret(dtype, Y), (M<<2,N<<1))

    @inbounds Yr[1:2:end>>1,1:end>>1] = @view(Xr[1:4:end,:])
    @inbounds Yr[2:2:end>>1,1:end>>1] = @view(Xr[2:4:end,:])
    @inbounds Yr[1+end>>1:2:end,1+end>>1:end] = @view(Xr[1:4:end,:])
    @inbounds Yr[2+end>>1:2:end,1+end>>1:end] = -@view(Xr[2:4:end,:])
    @inbounds Yr[1:2:end>>1,1+end>>1:end] = @view(Xr[3:4:end,:])
    @inbounds Yr[2:2:end>>1,1+end>>1:end] = @view(Xr[4:4:end,:])
    @inbounds Yr[1+end>>1:2:end,1:end>>1] = -@view(Xr[3:4:end,:])
    @inbounds Yr[2+end>>1:2:end,1:end>>1] = @view(Xr[4:4:end,:])

    return Y
end

function equiv(X::AbstractArray{T}) where T<:Quaternion
    Y = Array{Complex{real(eltype(X))}}(undef, size(X,1)<<1, size(X,2)<<1)
    return equiv!(Y, X)
end

function equiv!(Y::AbstractArray{Quaternion{T}}, X::AbstractArray{Complex{T}}) where T<:Real
    @assert size(X,1)>>1 == size(Y,1)
    @assert size(X,2)>>1 == size(Y,2)

    dtype = real(eltype(X))
    M = size(X,1)
    N = size(X,2)
    Xr = reshape(reinterpret(dtype, X), (M<<1,N))
    Yr = reshape(reinterpret(dtype, Y), (M<<1,N>>1))

    @inbounds Yr[1:4:end,:] = @view(Xr[1:2:end>>1,1:end>>1])
    @inbounds Yr[2:4:end,:] = @view(Xr[2:2:end>>1,1:end>>1])
    @inbounds Yr[3:4:end,:] = @view(Xr[1:2:end>>1,1+end>>1:end])
    @inbounds Yr[4:4:end,:] = @view(Xr[2:2:end>>1,1+end>>1:end])

    return Y
end

function equiv(X::AbstractArray{T}) where T<:Complex
    Y = Array{Quaternion{real(eltype(X))}}(undef, size(X,1)>>1, size(X,2)>>1)
    return equiv!(Y, X)
end

for fn in (:pinv, :inv)
    @eval begin
        function ($fn)(A::Array{Quaternion{T},2}) where T<:Real
            cmat = equiv(A)
            ret = ($fn)(cmat)
            return equiv(ret)
        end
    end
end

isreal(q::Quaternion) = iszero(imagi(q)) & iszero(imagj(q)) & iszero(imagk(q))
isinteger(q::Quaternion) = isreal(q) & isinteger(real(q))
isfinite(q::Quaternion) = isfinite(real(q)) & isfinite(imagi(q)) & isfinite(imagj(q)) & isfinite(imagk(q))
isnan(q::Quaternion) = isnan(real(q)) | isnan(imagi(q)) | isnan(imagj(q)) | isnan(imagk(q))
isinf(q::Quaternion) = isinf(real(q)) | isinf(imagi(q)) | isinf(imagj(q)) | isinf(imagk(q))
iszero(q::Quaternion) = iszero(real(q)) & iszero(imagi(q)) & iszero(imagj(q)) & iszero(imagk(q))
isone(q::Quaternion) = isone(real(q)) & iszero(imagi(q)) & iszero(imagj(q)) & iszero(imagk(q))

function round(z::Quaternion, rr::RoundingMode=RoundNearest,
               ri::RoundingMode=rr, rj::RoundingMode=rr, rk::RoundingMode=rr; kwargs...)
    Quaternion(round(real(z), rr; kwargs...),
               round(imagi(z), ri; kwargs...),
               round(imagj(z), rj; kwargs...),
               round(imagk(z), rk; kwargs...))
end

end


text:
struct BB{T1} <: QuasiNewton{T1}
 approx::T1
end
BB(;inverse=false) = BB(inverse ? Inverse() : Direct())
hasprecon(::BB) = NoPrecon()
update!(scheme::BB{<:Direct}, B, s, y) = _bb(scheme, B, s, y)
update(scheme::BB{<:Direct}, B, s, y) = _bb(scheme, B, s, y)
_bb(scheme::BB{<:Direct}, M, s, y) = dot(s, y)/dot(s, s)

ext: jl


text:
function passMasterSolution(m::Model)
    for ch in m.children
        ch.ext[:Stochastic].sol = m.colVal
    end
end

# Create sparse matrix A with all parent variables. Removes them from linear constraints afterwards.
# Once created, should be able to do A*x, where x is the solution to the master problem.
function createMasterMat(m::Model)
    stoch = getStructure(m)
    parent = stoch.parent
    for (nrow,con) in enumerate(m.linconstr)
        aff = con.terms
        for (id,var) in zip(reverse(aff.vars), length(aff.vars):-1:1)
            if var.m == parent
                push!(I, nrow)
                push!(J, var.col)
                push!(V, aff.coeff[id])
                splice!(aff.vars, id)
                splice!(aff.coeffs, id)
            end
        end
    end
    stoch.parentMat = sparse(I,J,V,length(m.linconstr),parent.numCols)
end

function prepProblemBounds(m::Model)

    objaff::AffExpr = m.obj.aff

    # We already have dense column lower and upper bounds

    # Create dense objective vector
    f = zeros(m.numCols)
    for ind in 1:length(objaff.vars)
        f[objaff.vars[ind].col] += objaff.coeffs[ind]
    end

    # Create row bounds
    numRows = length(m.linconstr)
    rowlb = fill(-Inf, numRows)
    rowub = fill(+Inf, numRows)

    stoch = getStructure(m)
    if stoch.parent != nothing && !isempty(stoch.parent.colVal)
        offset = stoch.parentMat * parent.colVal
        for c in 1:numRows
            rowlb[c] = m.linconstr[c].lb - offset[c]
            rowub[c] = m.linconstr[c].ub - offset[c]
        end
    else
        rowlb[c] = m.linconstr[c].lb
        rowub[c] = m.linconstr[c].ub
    end

    return f, rowlb, rowub
end

# Solve stochastic LP with Bender's decomposition
function solveStochastic(m::Model)
    stoch = getStructure(m)
    @assert stoch.parent == nothing # make sure we're at the master problem

    @variable(m, θ) # hopefully this name doesn't conflict...

    master_stat solve(m)
    @assert master_stat == :Optimal
    passMasterSolution(m)

    m.obj += θ

    for bl in stoch.children
        createMasterMat(bl)
    end

    while true
        for bl in stoch.children
            solve(bl)
            if bl.objVal > tol
            𝛑 = getconstrsolution(bl.internalModel)
            lhs = 𝛑 * getconstrmatrix(bl.internalModel)
            lb  = 𝛑 * getconstrLB(bl.internalModel) # what to do here...
            ub  = 𝛑 * getconstrUB(bl.internalModel) # and here...
            addconstr!(m.internalModel, Int64[1:length(lhs)], lhs, lb, ub)
        end
    end

end


text:
function kpp_psis(xi)
if xi>=0
  y=1+5*xi;
elseif (xi<0) & (xi>=-1.)
  y=(1-16*xi)^(-1/2);
else
  y=(-28.86-98.96*xi)^(-1/3);
end
return y
end


text:
### A Pluto.jl notebook ###
# v0.12.16

using Markdown
using InteractiveUtils

# ╔═╡ e077ab16-22e5-11eb-15ed-5f0654bf7328
using Pkg, DrWatson

# ╔═╡ e75a9a06-22e5-11eb-2e36-8d4f62b830ed
begin
	@quickactivate "StatisticsWithJuliaPlutoNotebooks"
	using Random, Distributions, Plots, LaTeXStrings
	Random.seed!(0)
end

# ╔═╡ ca5fbdc8-22e5-11eb-0a60-cb0a127a3ca5
md"## Listing5.09"

# ╔═╡ fa68607e-22e5-11eb-0558-c9a4d9f77426
begin
	actualAlpha, actualLambda = 2,3
	gammaDist = Gamma(actualAlpha,1/actualLambda)
	n = 10^2
	sample = rand(gammaDist, n)
end

# ╔═╡ 68b12c46-3a9c-11eb-03b3-bf9862c5371c
begin
	alphaGrid = 1:0.02:3
	lambdaGrid = 2:0.02:5
end

# ╔═╡ 68b16da0-3a9c-11eb-21e2-53c689f1599c
likelihood = [prod([pdf.(Gamma(a,1/l),v) for v in sample])
                        for l in lambdaGrid, a in alphaGrid]

# ╔═╡ 68b1f8ce-3a9c-11eb-2ad9-a5af4bd66c0e
surface(alphaGrid, lambdaGrid, likelihood, lw=0.1,
	c=cgrad([:blue, :red]), legend=:none, camera = (135,20),
	xlabel=L"\alpha", ylabel=L"\lambda", zlabel="Likelihood")

# ╔═╡ 475ff888-22e6-11eb-2354-09f8f40a8e12
md"## End of listing5.09"

# ╔═╡ Cell order:
# ╟─ca5fbdc8-22e5-11eb-0a60-cb0a127a3ca5
# ╠═e077ab16-22e5-11eb-15ed-5f0654bf7328
# ╠═e75a9a06-22e5-11eb-2e36-8d4f62b830ed
# ╠═fa68607e-22e5-11eb-0558-c9a4d9f77426
# ╠═68b12c46-3a9c-11eb-03b3-bf9862c5371c
# ╠═68b16da0-3a9c-11eb-21e2-53c689f1599c
# ╠═68b1f8ce-3a9c-11eb-2ad9-a5af4bd66c0e
# ╟─475ff888-22e6-11eb-2354-09f8f40a8e12


text:
### Forward-Backward Algorithm
function forward_backward(hmm::HMM,o::Vector; scaling=true)
	if scaling
		alpha, log_p_obs, coeff = forward(hmm,o; scaling=true)
		beta = backward(hmm,o; scale_coeff=coeff)
	else
		alpha, p_obs = forward(hmm,o; scaling=false)
		log_p_obs = log(p_obs)
	    beta = backward(hmm,o)
	end
	return alpha, beta, log_p_obs
end

function forward(hmm::HMM, o::Vector; scaling=true)
	n_obs = length(o)

	# alpha[t,i] = probability of being in state 'i' given o[1:t]
	alpha = zeros(n_obs, hmm.n)

	# base case (initialize at start)
	for i = 1:hmm.n
		alpha[1,i] = hmm.p[i] * pdf(hmm.B[i],o[1])
	end

	if scaling
		c = (Float64)[] # scaling coefficients
		push!(c,1./sum(alpha[1,:]))
		alpha[1,:] *= c[end]
	end

	# induction step
	for t = 2:n_obs
		for j = 1:hmm.n
			for i = 1:hmm.n
				alpha[t,j] += hmm.A[i,j] * alpha[t-1,i]
			end
			alpha[t,j] *= pdf(hmm.B[j],o[t])
		end
		if scaling
			push!(c,1./sum(alpha[t,:]))
			alpha[t,:] *= c[end]
		end
	end

	# Calculate likelihood (or log-likelihood) of observed sequence
	if scaling
		log_p_obs = -sum(log(c)) # see Rabiner (1989), eqn 103
		return (alpha,log_p_obs,c)
	else
		p_obs = sum(alpha[end,:])
		return (alpha,p_obs)
	end
end

function backward(hmm::HMM, o::Vector; scale_coeff=nothing)
	# scale_coeff are 1/sum(alpha[t,:]) calculated by forward algorithm
	n_obs = length(o)

	# beta[t,i] = probability of being in state 'i' and then obseverving o[t+1:end]
	beta = zeros(n_obs, hmm.n)

	# base case (initialize at end)
	if scale_coeff == nothing
		beta[end,:] += 1
	else
		if length(scale_coeff) != n_obs
			error("scale_coeff is improperly defined (wrong length)")
		end
		beta[end,:] += scale_coeff[end]
	end

	# induction step
	for t = reverse(1:n_obs-1)
		for i = 1:hmm.n
			for j = 1:hmm.n
				beta[t,i] += hmm.A[i,j] * pdf(hmm.B[j],o[t+1]) * beta[t+1,j]
			end
		end
		if scale_coeff != nothing
			beta[t,:] *= scale_coeff[t]
		end
	end

	return beta
end

### Viterbi

function viterbi(hmm::HMM, o::Vector)
	n_obs = length(o)

	# delta[i,j] = highest probability of state sequence ending in state j on step i
	# psi[i,j] = most likely state on step i-1 given state j on step i (argmax of deltas)
	delta = zeros(n_obs, hmm.n)
	psi = ones(Int, n_obs, hmm.n)

	# base case, psi[:,1] is ignored so don't initialize
	for i = 1:hmm.n
		delta[1,i] = hmm.p[i] .* pdf(hmm.B[i],o[1])
	end

	# induction step
	for t = 2:n_obs
		for j = 1:hmm.n
			delta[t,j],psi[t,j] = findmax(hmm.A[:,j].*delta[t-1,:]')
			delta[t,j] *= pdf(hmm.B[j],o[t])
		end
	end

	# backtrack to uncover the most likely path / state sequence
	q = zeros(Int,n_obs) # vector holding state sequence
	q[end] = indmax(delta[end,:])

	# backtrack recursively
	for t = reverse(1:n_obs-1)
		q[t] = psi[t+1,q[t+1]]
	end
	return q
end

### Smoothed state sequence based on forward-backward algorithm
function smoothed_forward_backward(hmm::HMM, o::Vector; scaling=true)
	log_p_obs,alpha,beta,x,g = calc_stats(hmm,o;scaling=scaling)
	return g,log_p_obs
end


text:
function test14()
  filename = "test14.asc"
  # exectuablepath = null string will not run LTspice.exe.  Test parsing only.
  sim = LTspiceSimulation(filename,executablepath="")
  @test LTspice.does_circuitfilearray_file_match(sim)
  show(IOBuffer(),sim)
  @test measurementnames(sim) == ("m1",)
  @test parameternames(sim) == ("a","b","d", "e", "φ", "Ω", "this_is_a_long_name")
end
test14()


text:
global prev_theme=nothing

if haskey(ENV, "GADFLY_THEME")
    prev_theme = ENV["GADFLY_THEME"]
    pop!(ENV, "GADFLY_THEME")
end

using Base.Test, Gadfly, Compat

repo = LibGit2.GitRepo(dirname(@__DIR__))
branch = LibGit2.headname(repo)
outputdir = joinpath(@__DIR__, mapreduce(x->startswith(branch,x), |, ["master","(detac"]) ?
        "cachedoutput" : "gennedoutput")

options = LibGit2.StatusOptions(flags=LibGit2.Consts.STATUS_OPT_INCLUDE_IGNORED |
                                      LibGit2.Consts.STATUS_OPT_RECURSE_IGNORED_DIRS)
status = LibGit2.GitStatus(repo, status_opts=options)
for i in 1:length(status)
    entry = status[i]
    index_to_workdir = unsafe_load(entry.index_to_workdir)
    if index_to_workdir.status == Int(LibGit2.Consts.DELTA_IGNORED)
        filepath = unsafe_string(index_to_workdir.new_file.path)
        startswith(filepath,joinpath("test",outputdir)) || continue
        rm(joinpath(dirname(@__DIR__),filepath))
    end
end
function mimic_git_log_n1(io::IO, head)
    hash = LibGit2.GitHash(head)
    println(io, "commit ",string(hash))
    commit = LibGit2.GitCommit(repo, hash)
    author = LibGit2.author(commit)
    println(io, "Author: ",author.name," <",author.email,">")
    datetime = Dates.unix2datetime(author.time + 60*author.time_offset)
    println(io, "Date:   ",Dates.format(datetime, "e u d HH:MM:SS YYYY"))
    println(io, "    ",LibGit2.message(commit))
end
function mimic_git_status(io::IO, head)
    println(io, "On branch ",LibGit2.shortname(head))
    status = LibGit2.GitStatus(repo)
    println(io, "Changes not staged for commit:")
    for i in 1:length(status)
        entry = status[i]
        index_to_workdir = unsafe_load(entry.index_to_workdir)
        if index_to_workdir.status == Int(LibGit2.Consts.DELTA_MODIFIED)
            println(io, "    ", unsafe_string(index_to_workdir.new_file.path))
        end
    end
    println(io, "Untracked files:")
    for i in 1:length(status)
        entry = status[i]
        index_to_workdir = unsafe_load(entry.index_to_workdir)
        if index_to_workdir.status == Int(LibGit2.Consts.DELTA_UNTRACKED)
            println(io, "    ", unsafe_string(index_to_workdir.new_file.path))
        end
    end
end
head = LibGit2.head(repo)
open(io->mimic_git_log_n1(io,head), joinpath(outputdir,"git.log"), "w")
open(io->mimic_git_status(io,head), joinpath(outputdir,"git.status"), "w")

backends = Dict{AbstractString, Function}(
    "svg" => name -> SVG(joinpath(outputdir,"$(name).svg")),
    "svgjs" => name -> SVGJS(joinpath(outputdir,"$(name).js.svg"), jsmode=:linkabs),
    "png" => name -> PNG(joinpath(outputdir,"$(name).png")),
    "ps"  => name -> PS(joinpath(outputdir,"$(name).ps")),
    "pdf" => name -> PDF(joinpath(outputdir,"$(name).pdf")),
    "pgf" => name -> PGF(joinpath(outputdir,"$(name).tex"))
)

testdir = joinpath((@__DIR__),"testscripts")
testfiles = isempty(ARGS) ?
        [splitext(filename)[1] for filename in readdir(testdir) if filename[1]!='.'] :
        ARGS

@testset "Gadfly" begin
    for filename in testfiles, (backend_name, backend) in backends
        info(filename,'.',backend_name)
        srand(1)
        p = evalfile(joinpath(testdir, "$(filename).jl"))
        @test typeof(p) in [Plot,Compose.Context]
        r = draw(backend(filename), p)
        @test typeof(r) in [Bool,Void]
    end
end

output = open(string(outputdir,".html"), "w")
print(output,
    """
    <!DOCTYPE html>

    <meta charset="utf-8" />

        <title>Gadfly Test Plots</title>

    """)

for filename in testfiles
    println(output, "
", filename, "
")
    print(output, """
<object type="image/svg+xml" data="$(outputdir)/$(filename).js.svg"></object>
""")
    print(output, """""")
    print(output, """\n""")
end

print(output,
    """

    """)

close(output)

if prev_theme !== nothing
    ENV["GADFLY_THEME"] = prev_theme
end


text:
# GFF3 File Format
# ================

module GFF3

import Automa
import Automa.RegExp: @re_str
import BufferedStreams
import Bio.Exceptions: missingerror
import URIParser
importall Bio

include("record.jl")
include("reader.jl")
include("writer.jl")

end


text:
module QuantumCumulants

import SymbolicUtils
import SymbolicUtils: substitute

import Symbolics

import SciMLBase

import ModelingToolkit
const MTK = ModelingToolkit
import ModelingToolkit: ⊗ # just to avoid conflicts

using Combinatorics: partitions, combinations

export HilbertSpace, ProductSpace, ⊗, tensor,
        QSym, QTerm, @qnumbers,
        FockSpace, Destroy, Create,
        NLevelSpace, Transition,
        MeanfieldEquations,
        meanfield, commutator, acts_on,
        CNumber, Parameter, @cnumbers, cnumbers,
        Average, average, cumulant_expansion, get_order, cumulant,
        find_missing, complete, complete!, find_operators, fundamental_operators,
            unique_ops, unique_ops!,
        CorrelationFunction, Spectrum, correlation_u0, correlation_p0,
        ClusterSpace,
        scale,
        transition_superscript

include("hilbertspace.jl")
include("qnumber.jl")
include("cnumber.jl")
include("fock.jl")
include("nlevel.jl")
include("equations.jl")
include("meanfield.jl")
include("average.jl")
include("utils.jl")
include("diffeq.jl")
include("correlation.jl")
include("cluster.jl")
include("scale.jl")
include("latexify_recipes.jl")
include("printing.jl")

@deprecate heisenberg(args...; kwargs...) meanfield(args...; kwargs...)

end # module


text:
"""
function fix_events()

Perform some fixes on the events such that the spadl action dataframe can be built.
"""
function fix_events(vector_wyscout_data::Vector{WyscoutData})::Vector{WyscoutData}
    vector_wyscout_data = create_shot_coordinates(vector_wyscout_data)
    vector_wyscout_data = convert_duels(vector_wyscout_data)
    vector_wyscout_data = insert_interception_passes(vector_wyscout_data)
    vector_wyscout_data = add_offside_variable(vector_wyscout_data)
    vector_wyscout_data = convert_touches(vector_wyscout_data)

    return vector_wyscout_data
end

"""
function create_shot_coordinates()

Create short coordinates (estimates) from Wyscout tags
# """
function create_shot_coordinates(vector_wyscout_data::Vector{WyscoutData})::Vector{WyscoutData}

    for i in eachindex(vector_wyscout_data)

        if vector_wyscout_data[i].tags.position_goal_low_center || vector_wyscout_data[i].tags.position_goal_mid_center || vector_wyscout_data[i].tags.position_goal_high_center
            vector_wyscout_data[i].event_fixed.end_x = 100
            vector_wyscout_data[i].event_fixed.end_y = 50
        elseif vector_wyscout_data[i].tags.position_goal_low_right || vector_wyscout_data[i].tags.position_goal_mid_right || vector_wyscout_data[i].tags.position_goal_high_right
            vector_wyscout_data[i].event_fixed.end_x = 100
            vector_wyscout_data[i].event_fixed.end_y = 55
        elseif vector_wyscout_data[i].tags.position_goal_mid_left || vector_wyscout_data[i].tags.position_goal_low_left || vector_wyscout_data[i].tags.position_goal_high_left
            vector_wyscout_data[i].event_fixed.end_x = 100
            vector_wyscout_data[i].event_fixed.end_y = 45
        elseif vector_wyscout_data[i].tags.position_out_high_center || vector_wyscout_data[i].tags.position_post_high_center
            vector_wyscout_data[i].event_fixed.end_x = 100
            vector_wyscout_data[i].event_fixed.end_y = 50
        elseif vector_wyscout_data[i].tags.position_out_low_right || vector_wyscout_data[i].tags.position_out_mid_right || vector_wyscout_data[i].tags.position_out_high_right
            vector_wyscout_data[i].event_fixed.end_x = 100
            vector_wyscout_data[i].event_fixed.end_y = 60
        elseif vector_wyscout_data[i].tags.position_out_mid_left || vector_wyscout_data[i].tags.position_out_low_left || vector_wyscout_data[i].tags.position_out_high_left
            vector_wyscout_data[i].event_fixed.end_x = 100
            vector_wyscout_data[i].event_fixed.end_y = 40
        elseif vector_wyscout_data[i].tags.position_post_mid_left || vector_wyscout_data[i].tags.position_post_low_left || vector_wyscout_data[i].tags.position_post_high_left
            vector_wyscout_data[i].event_fixed.end_x = 100
            vector_wyscout_data[i].event_fixed.end_y = 55.38
        elseif vector_wyscout_data[i].tags.position_post_low_right || vector_wyscout_data[i].tags.position_post_mid_right || vector_wyscout_data[i].tags.position_post_high_right
            vector_wyscout_data[i].event_fixed.end_x = 100
            vector_wyscout_data[i].event_fixed.end_y = 44.62
        elseif vector_wyscout_data[i].tags.blocked
            vector_wyscout_data[i].event_fixed.end_x = vector_wyscout_data[i].event_fixed.start_x
            vector_wyscout_data[i].event_fixed.end_y = vector_wyscout_data[i].event_fixed.start_y
        end
    end

    return vector_wyscout_data
end

"""
function convert_duels()

This function converts Wyscout duels that end with the ball out of field
(subtype_id 50) into a pass for the player winning the duel to the location
of where the ball went out of field. The remaining duels are removed as
they are not on-the-ball actions.
"""
function convert_duels(vector_wyscout_data::Vector{WyscoutData})

    size_condition = length(eachindex(vector_wyscout_data)) - 2

    for i in eachindex(vector_wyscout_data)

        if i < size_condition
            # define salector for same period id
            selector_same_period = (vector_wyscout_data[i].event.period_id == vector_wyscout_data[i+2].event.period_id)

            # Define selector for duels that are followed by an 'out of field' event
            selector_duel_out_of_field = (vector_wyscout_data[i].event.type_id == 1 &&
                                            vector_wyscout_data[i+1].event.type_id == 1 &&
                                            vector_wyscout_data[i+2].event.subtype_id == 50 &&
                                            selector_same_period)

            # Define selector for current time step
            selector_0_duel_won = selector_duel_out_of_field && (vector_wyscout_data[i].event.team_id != vector_wyscout_data[i+2].event.team_id)
            selector_0_duel_won_air = selector_0_duel_won && (vector_wyscout_data[i].event.subtype_id == 10)
            selector_0_duel_won_not_air = selector_0_duel_won && (vector_wyscout_data[i].event.subtype_id != 10)

            # Define selector for next time step
            selector_1_duel_won = selector_duel_out_of_field && (vector_wyscout_data[i+1].event.team_id != vector_wyscout_data[i+2].event.team_id)
            selector_1_duel_won_air = selector_1_duel_won && (vector_wyscout_data[i+1].event.subtype_id == 10)
            selector_1_duel_won_not_air = selector_1_duel_won && (vector_wyscout_data[i+1].event.subtype_id != 10)

            # Aggregate selectors
            selector_duel_won = selector_0_duel_won || selector_1_duel_won
            selector_duel_won_air = selector_0_duel_won_air || selector_1_duel_won_air
            selector_duel_won_not_air = selector_0_duel_won_not_air || selector_1_duel_won_not_air

            # Set types and subtypes
            if selector_duel_won vector_wyscout_data[i].event_fixed.type_id = 8 end
            if selector_duel_won_air vector_wyscout_data[i].event_fixed.subtype_id = 82 end
            if selector_duel_won_not_air vector_wyscout_data[i].event_fixed.subtype_id = 85 end

            # set end location equal to ball out of field location
            if selector_duel_won vector_wyscout_data[i].tags.accurate = false end
            if selector_duel_won vector_wyscout_data[i].tags.not_accurate = true end
            if selector_duel_won vector_wyscout_data[i].event_fixed.end_x = 100 - vector_wyscout_data[i+2].event_fixed.start_x end
            if selector_duel_won vector_wyscout_data[i].event_fixed.end_y = 100 - vector_wyscout_data[i+2].event_fixed.start_y end

            # Define selector for ground attacking duels with take on
            selector_attacking_duel = (vector_wyscout_data[i].event.subtype_id == 11)
            selector_take_on = (vector_wyscout_data[i].tags.take_on_left || vector_wyscout_data[i].tags.take_on_right)
            selector_att_duel_take_on = selector_attacking_duel && selector_take_on

            # set take on type to 0
            if selector_att_duel_take_on vector_wyscout_data[i].event_fixed.type_id = 0 end
            if vector_wyscout_data[i].tags.sliding_tackle vector_wyscout_data[i].event_fixed.type_id = 0 end

        end
    end
    # drop the remaining duels
    vector_wyscout_data = filter(el -> el.event_fixed.type_id != 1, vector_wyscout_data)
    return vector_wyscout_data
end

"""
    function insert_interception_passes()
This function converts passes (type_id 8) that are also interceptions
(tag interception) in the Wyscout event data into two separate events,
first an interception and then a pass.
"""
function insert_interception_passes(vector_wyscout_data::Vector{WyscoutData})::Vector{WyscoutData}

    new_vector_wyscout_data = Vector{WyscoutData}()

    for i in eachindex(vector_wyscout_data)

        push!(new_vector_wyscout_data, vector_wyscout_data[i])
        interception = vector_wyscout_data[i].tags.interception && vector_wyscout_data[i].event_fixed.type_id == 8

        if interception
            new_wyscout_data = vector_wyscout_data[i]
            new_wyscout_data.tags.interception = true
            new_wyscout_data.event_fixed.type_id = 0
            new_wyscout_data.event_fixed.type_id = 0
            new_wyscout_data.event_fixed.end_x, new_wyscout_data.event_fixed.end_y = new_wyscout_data.event_fixed.start_x, new_wyscout_data.event_fixed.start_y
            push!(new_vector_wyscout_data, new_wyscout_data)
        end
    end

    return vector_wyscout_data
end

"""
    function add_offside_variable
This function removes the offside events in the Wyscout event data and adds
sets offside to 1 for the previous event (if this was a passing event)
"""
function add_offside_variable(vector_wyscout_data::Vector{WyscoutData})::Vector{WyscoutData}

    for i in eachindex(vector_wyscout_data)
        if i < length(vector_wyscout_data) - 1
        # Selector_offside_pass
            selector_offside = (vector_wyscout_data[i+1].event_fixed.type_id == 6) && (vector_wyscout_data[i].event_fixed.type_id == 8)
            # set boolean offside to true if offside pass
            if selector_offside
                vector_wyscout_data[i].tags.offside = true
            end
        end
    end

    # remove offside events
    vector_wyscout_data = filter(x -> x.event_fixed.type_id != 6, vector_wyscout_data)

    return vector_wyscout_data
end

"""
    function convert_touches()

This function converts the Wyscout 'touch' event (sub_type_id 72) into either
a dribble or a pass (accurate or not depending on receiver)
"""
function convert_touches(vector_wyscout_data::Vector{WyscoutData})::Vector{WyscoutData}
    min_dribble_length = 3.0

    for i in eachindex(vector_wyscout_data)
        if i < length(vector_wyscout_data) - 1

            # select touch data
            selector_touch = (vector_wyscout_data[i].event_fixed.subtype_id == 72) & !(vector_wyscout_data[i].tags.interception)

            selector_same_player = (vector_wyscout_data[i].event_fixed.player_id == vector_wyscout_data[i+1].event_fixed.player_id)
            selector_same_team = (vector_wyscout_data[i].event_fixed.team_id == vector_wyscout_data[i+1].event_fixed.team_id)

            selector_touch_same_team = (selector_touch && !selector_same_player && selector_same_team)
            selector_touch_other = (selector_touch && !selector_same_player && !selector_same_team)

            same_x = abs(vector_wyscout_data[i].event_fixed.end_x - vector_wyscout_data[i+1].event_fixed.start_x) < min_dribble_length
            same_y = abs(vector_wyscout_data[i].event_fixed.end_y - vector_wyscout_data[i+1].event_fixed.start_y) < min_dribble_length
            same_loc = same_x && same_y

            if selector_touch_same_team && same_loc
                vector_wyscout_data[i].event_fixed.type_id = 8
                vector_wyscout_data[i].event_fixed.subtype_id = 85
                vector_wyscout_data[i].tags.accurate = true
                vector_wyscout_data[i].tags.not_accurate = false
            end

            if selector_touch_other && same_loc
                vector_wyscout_data[i].event_fixed.type_id = 8
                vector_wyscout_data[i].event_fixed.subtype_id = 85
                vector_wyscout_data[i].tags.accurate = false
                vector_wyscout_data[i].tags.not_accurate = true
            end

        end

    end
    return vector_wyscout_data

end


text:
using Base.Enums
using MLStyle.Render
using DataStructures

struct Field
    name      :: Any  # an expr to represent the field name from IN_FIELDS.
    make      :: Any
    var       :: Symbol
    typ       :: Any
end

function query_routine(assigns            :: OrderedDict{Symbol, Any},
                       fn_in_fields       :: Vector{Field},
                       fn_returns         :: Any,
                       result; infer_type = true)
    @assert haskey(assigns, FN_OUT_FIELDS)

    fn_arguments = map(x -> x.var, fn_in_fields)
    fn_arg_types = Expr(:vect, map(x -> x.typ, fn_in_fields)...)

    function (inner_expr)
        let_seq = [
            Expr(:(=), Expr(:tuple, IN_FIELDS, IN_TYPES, IN_SOURCE), inner_expr),
            (:($name = $value) for (name, value) in assigns)...,
            :(@inline $FN($(fn_arguments...)) =  $fn_returns),
        ]
        if infer_type
            let type_infer = :($FN_RETURN_TYPES = $type_unpack($length($FN_OUT_FIELDS, ), $return_type($FN, $fn_arg_types)))
                push!(let_seq, type_infer)
            end
        end
        Expr(:let,
            Expr(
                :block,
                let_seq...
            ),
            result
        )
    end
end

# visitor to process the pattern `_.x, _,"x", _.(1)` inside an expression
function mk_visit(fields :: Dict{Any, Field}, assigns :: OrderedDict{Symbol, Any})
    visit = expr ->
    @match expr begin
        Expr(:. , :_, q :: QuoteNode) && Do(a = q.value) ||
        Expr(:., :_, Expr(:tuple, a)) =>
            @match a begin
                a :: Int =>
                    let field = get!(fields, a) do
                            var_sym = gen_sym(a)
                            Field(
                                a,
                                Expr(:ref, RECORD, a),
                                var_sym,
                                Expr(:ref, IN_TYPES, a)
                            )
                        end
                        field.var
                    end

                ::String && Do(b = Symbol(a)) ||
                b::Symbol =>
                    let field = get!(fields, b) do
                            idx_sym = gen_sym()
                            var_sym = gen_sym(b)
                            assigns[idx_sym] = Expr(:call, findfirst, x -> x === b, IN_FIELDS)
                            Field(
                                b,
                                Expr(:ref, RECORD, idx_sym),
                                var_sym,
                                Expr(:ref, IN_TYPES, idx_sym)
                            )
                        end
                        field.var
                    end
            end
        Expr(head, args...) => Expr(head, map(visit, args)...)
        a => a
    end
end

function fn_apply(fields :: Vector{Field})
    Expr(
        :let,
        let let_seq = (:($(field.var) = $(field.make)) for field in fields)
            Expr(:block, let_seq...)
        end,
        let args = (field.var for field in fields)
            :($FN($(args...)))
        end
    )
end

function generate_where(args :: AbstractArray)

    pred_in_fields = Dict{Any, Field}()
    assigns        = OrderedDict{Symbol, Any}(FN_OUT_FIELDS => IN_FIELDS)
    visit = mk_visit(pred_in_fields, assigns)
    pred = foldl(args, init=true) do last, arg
        boolean = visit(arg)
        if last === true
            boolean
        else
            Expr(:&&, last, boolean)
        end
    end

    fields :: Vector{Field} = pred_in_fields |> values |> collect
    # where expression generation
    query_routine(
        assigns,
        fields,
        pred,
        Expr(:tuple,
             IN_FIELDS,
             IN_TYPES,
             :($RECORD for $RECORD in $IN_SOURCE if $(fn_apply(fields)))
        );
        infer_type = false
    )
end

function generate_select(args :: AbstractArray)
    map_in_fields = Dict{Any, Field}()
    assigns = OrderedDict{Symbol, Any}()
    fn_return_elts   :: Vector{Any} = []
    fn_return_fields :: Vector{Any} = []
    visit = mk_visit(map_in_fields, assigns)
    # process selectors
    predicate_process(arg) =
        @match arg begin
        :(!$pred($ (args...) )) && Do(ab=true) ||
        :($pred($ (args...) )) && Do(ab=false) ||
        :(!$pred) && Do(ab=true, args=[])      ||
        :($pred) && Do(ab=false, args=[])      =>
            let idx_sym = gen_sym()
                assigns[idx_sym] =
                    Expr(
                        :call,
                        findall,
                        ab ?
                            :(@inline function ($ARG,) !$pred($string($ARG,), $(args...)) end) :
                            :(@inline function ($ARG,) $pred($string($ARG,), $(args...)) end)
                        , IN_FIELDS
                    )
                idx_sym
            end
        end
    foreach(args) do arg
        @match arg begin
            :_ =>
                let field = get!(map_in_fields, all) do
                        var_sym = gen_sym()
                        push!(fn_return_elts, Expr(:..., var_sym))
                        push!(fn_return_fields, Expr(:..., IN_FIELDS))
                        Field(
                            all,
                            RECORD,
                            var_sym,
                            :($Tuple{$IN_TYPES...})
                        )
                    end
                    nothing
                end

            :(_.($(args...))) =>
                let indices = map(predicate_process, args)
                    if haskey(map_in_fields, arg)
                        throw("The columns `$(string(arg))` are selected twice!")
                    elseif !isempty(indices)
                        idx_sym = gen_sym()
                        var_sym = gen_sym()
                        field = begin
                            assigns[idx_sym] =
                                length(indices) === 1 ?
                                indices[1] :
                                Expr(:call, intersect, indices...)
                            push!(fn_return_elts, Expr(:..., var_sym))
                            push!(fn_return_fields, Expr(:..., Expr(:ref, IN_FIELDS, idx_sym)))
                            Field(
                                arg,
                                Expr(:ref, RECORD, idx_sym),
                                var_sym,
                                Expr(:curly, Tuple, Expr(:..., Expr(:ref, IN_TYPES, idx_sym)))
                            )
                        end
                        map_in_fields[arg] = field
                        nothing
                    end
                end
           :($a => $new_field) || a && Do(new_field = Symbol(string(a))) =>
                let new_value = visit(a)
                    push!(fn_return_fields, QuoteNode(new_field))
                    push!(fn_return_elts, new_value)
                    nothing
                end
        end
    end
    fields = map_in_fields |> values |> collect
    assigns[FN_OUT_FIELDS] = Expr(:vect, fn_return_fields...)
    # select expression generation
    query_routine(
        assigns,
        fields,
        Expr(:tuple, fn_return_elts...),
        Expr(
            :tuple,
            FN_OUT_FIELDS,
            FN_RETURN_TYPES,
            :($(fn_apply(fields)) for $RECORD in $IN_SOURCE)
        ); infer_type = true
    )
end

function generate_groupby(args :: AbstractArray, having_args :: Union{Nothing, AbstractArray} = nothing)
    group_in_fields = Dict{Any, Field}()
    having_in_fields = Dict{Any, Field}()
    assigns = OrderedDict{Symbol, Any}()
    visit_group_fn = mk_visit(group_in_fields, assigns)
    visit_having = mk_visit(having_in_fields, assigns)

    group_fn_return_vars = []
    group_fn_return_elts = []

    foreach(args) do arg
        @match arg begin
            :($b => $a) || b && Do(a = Symbol(string(b))) =>
                let field = a
                    push!(group_fn_return_vars, a)
                    push!(group_fn_return_elts, visit_group_fn(b))
                    nothing
                end
        end
    end

    cond_expr = having_args === nothing ? nothing :
        let cond =
            foldl(having_args) do last, arg
                Expr(:&&, last, visit_having(arg))
            end

            quote
                if !($cond)
                    continue
                end
            end
        end

    group_fn_return_fields = map(QuoteNode, group_fn_return_vars)
    ngroup_key = length(group_fn_return_vars)

    group_fields = group_in_fields |> values |> collect
    bindings_inside_for_loop = [
        (:($(field.var) = $(field.make)) for field in group_fields)...,
        (:($(field.var) = $(field.make)) for (key, field) in having_in_fields if !haskey(group_in_fields, key))...
    ]

    assigns[FN_OUT_FIELDS] = Expr(:vect, group_fn_return_fields...)
    group_fn_required_vars = [field.var for field in values(group_in_fields)]

    group_key_lhs = Expr(:tuple, group_fn_return_vars...)
    group_key_rhs = Expr(:call, FN, group_fn_required_vars...)

    # groupby expression generation
    query_routine(
        assigns,
        group_fields,
        Expr(:tuple, group_fn_return_elts...),
        let out_fields = Expr(:vcat, FN_OUT_FIELDS, IN_FIELDS),
            out_types = Expr(:vcat, FN_RETURN_TYPES, AGG_TYPES)
            quote
                $AGG_TYPES = [$Vector{each} for each in $IN_TYPES]
                $GROUPS = $Dict{$Tuple{$FN_RETURN_TYPES...}, $Tuple{$AGG_TYPES...}}()
                $N = $length($IN_FIELDS,)
                for $RECORD in $IN_SOURCE
                    $(bindings_inside_for_loop...)
                    $GROUP_KEY = $group_key_lhs = $group_key_rhs
                    $cond_expr
                    $AGG =
                        $get!($GROUPS, $GROUP_KEY) do
                            Tuple(typ[] for typ in $IN_TYPES)
                        end
                    ($push!).($AGG, $RECORD,)
                end
                (
                    $out_fields,
                    $out_types,
                    ((k..., v...) for (k, v) in $GROUPS)
                )
            end
        end
    )
end


text:
using Documenter, MixedModelsBLB

# ENV["DOCUMENTER_DEBUG"] = "true"

# makedocs(
#     format = Documenter.HTML(),
#     sitename = "MixedModelsBLB",
#     modules = [MixedModelsBLB]
# )

# deploydocs(
#     repo   = "https://github.com/xinkai-zhou/MixedModelsBLB.jl.git",
#     target = "build"
# )

# using Documenter: Documenter, makedocs, deploydocs
# using PkgTemplates: PkgTemplates

makedocs(
    modules=[MixedModelsBLB],
    authors="Xinkai Zhou, Hua Zhou",
    repo="https://github.com/xinkai-zhou/MixedModelsBLB.jl.git/blob/{commit}{path}#{line}",
    sitename="MixedModelsBLB.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://xinkai-zhou.github.io/MixedModelsBLB.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "Installation" => "man/installation.md",
        "Quick Tutorial" => "man/quick-tutorial.md",
        "Detailed Usage" => "man/detailed-usage.md",
        "Contributing" => "man/contributing.md",
        "API" => "man/api.md",
    ],
)

deploydocs(;
    repo="github.com/xinkai-zhou/MixedModelsBLB.jl.git",
    target = "build"
)


text:
################################################################################
# Julia C wrapper generator using libclang from the LLVM project               #
################################################################################

module wrap_c
    version = v"0.0.0"

using Clang.cindex
using DataStructures
using Compat

export wrap_c_headers
export WrapContext

# Reserved Julia identifiers will be prepended with "_"
reserved_words = [ "abstract", "baremodule", "begin", "bitstype", "break", "catch", "ccall",
                   "const", "continue", "do", "else", "elseif", "end", "export", "finally",
                   "for", "function", "global", "if", "immutable", "import", "importall", "in",
                   "let", "local", "macro", "module", "quote", "return", "try", "type",
                   "typealias", "using", "while"]

# These argument types are unsupported
reserved_argtypes = ["va_list"]

# These argument types will be untyped in the Julia signature
untyped_argtypes = [IncompleteArray]

### InternalOptions
type InternalOptions
    wrap_structs::Bool
    immutable_structs::Bool
end
InternalOptions() = InternalOptions(true, false)

type ExprUnit
    items::Array{Any}
    deps::OrderedSet{Symbol}
    state::Symbol
end

immutable Poisoned
end

ExprUnit() = ExprUnit(Any[], OrderedSet{Symbol}(), :new)
ExprUnit(e::Union(Expr,Symbol,ASCIIString,Poisoned), deps=Any[]; state::Symbol=:new) = ExprUnit(Any[e], OrderedSet{Symbol}([target_type(dep) for dep in deps]...), state)
ExprUnit(a::Array, deps=Any[]; state::Symbol=:new) = ExprUnit(a, OrderedSet{Symbol}([target_type(dep) for dep in deps]...), state)

### WrapContext
# stores shared information about the wrapping session
type WrapContext
    index::cindex.CXIndex
    headers::Array{ASCIIString,1}
    output_file::ASCIIString
    common_file::ASCIIString
    clang_includes::Array{ASCIIString,1}         # clang include paths
    clang_args::Array{ASCIIString,1}             # additional {"-Arg", "value"} pairs for clang
    header_wrapped::Function                     # called to determine header inclusion status
                                                 #   (top_header, cursor_header) -> Bool
    header_library::Function                     # called to determine shared library for given header
                                                 #   (header_name) -> library_name::String
    header_outputfile::Function                  # called to determine output file group for given header
                                                 #   (header_name) -> output_file::String
    cursor_wrapped::Function                     # called to determine cursor inclusion status
                                                 #   (cursor_name, cursor) -> Bool
    common_buf::OrderedDict{Symbol, ExprUnit}    # output buffer for common items: typedefs, enums, etc.
    empty_structs::Set{ASCIIString}
    output_bufs::DefaultOrderedDict{ASCIIString, Array{Any}}
    options::InternalOptions
    anon_count::Int
    rewriter::Function
end

### Convenience function to initialize wrapping context with defaults
init(;args...) = init(ASCIIString[]; args...)
function init(;
            headers                         = ASCIIString[],
            index                           = None,
            output_file::ByteString        = "",
            common_file::ByteString        = "",
            output_dir::ByteString         = "",
            clang_args::Array{ASCIIString,1}
                                            = ASCIIString[],
            clang_includes::Array{ASCIIString,1}
                                            = ASCIIString[],
            clang_diagnostics::Bool         = true,
            header_wrapped                  = (header, cursorname) -> true,
            header_library                  = None,
            header_outputfile               = None,
            cursor_wrapped                  = (cursorname, cursor) -> true,
            options                         = InternalOptions(),
            rewriter                        = x -> x)

    # Set up some optional args if they are not explicitly passed.

    (index == None)         && ( index = cindex.idx_create(0, (clang_diagnostics ? 1 : 0)) )

    if (output_file == "" && header_outputfile == None)
        header_outputfile = x->joinpath(output_dir, strip(splitext(basename(x))[1]) * ".jl")
    end

    (common_file == "")    && ( common_file = output_file )
    common_file = joinpath(output_dir, common_file)

    if (header_library == None)
        header_library = x->strip(splitext(basename(x))[1])
    elseif isa(header_library, ASCIIString)
        libname = copy(header_library)
        header_library = x->libname
    end
    if (header_outputfile == None)
        header_outputfile = x->joinpath(output_dir, output_file)
    end

    # Instantiate and return the WrapContext
    global context = WrapContext(index,
                                 headers,
                                 output_file,
                                 common_file,
                                 clang_includes,
                                 clang_args,
                                 header_wrapped,
                                 header_library,
                                 header_outputfile,
                                 cursor_wrapped,
                                 OrderedDict{Symbol,ExprUnit}(),
                                 Set{ASCIIString}(),
                                 DefaultOrderedDict(ASCIIString, Array{Any}, ()->Any[]),
                                 options,
                                 0,
                                 rewriter)
    return context
end

###############################################################################
#
# Mapping from libclang types to Julia types
#
#   This is primarily used with CXTypeKind enums (aka: cindex.TypeKind)
#   but can also be used to map arbitrary type names to corresponding
#   Julia entities, for example "size_t" -> :Csize_t
#
##############################################################################
cl_to_jl = @compat Dict{Any,Any}(
    cindex.VoidType         => :Void,
    cindex.BoolType         => :Bool,
    cindex.Char_U           => :Uint8,
    cindex.UChar            => :Cuchar,
    cindex.Char16           => :Uint16,
    cindex.Char32           => :Uint32,
    cindex.UShort           => :Uint16,
    cindex.UInt             => :Uint32,
    cindex.ULong            => :Culong,
    cindex.ULongLong        => :Culonglong,
    cindex.Char_S           => :Uint8,
    cindex.SChar            => :Uint8,
    cindex.WChar            => :Char,
    cindex.Short            => :Int16,
    cindex.IntType          => :Cint,
    cindex.Long             => :Clong,
    cindex.LongLong         => :Clonglong,
    cindex.Float            => :Cfloat,
    cindex.Double           => :Cdouble,
    cindex.LongDouble       => :Float64,
    cindex.Enum             => :Cint,
    cindex.NullPtr          => :C_NULL,
    cindex.Int128           => :Int128,
    cindex.UInt128          => :Uint128,
    cindex.FirstBuiltin     => :Void,
    :size_t                 => :Csize_t,
    :ptrdiff_t              => :Cptrdiff_t,
    :uint64_t               => :Uint64,
    :uint32_t               => :Uint32,
    :uint16_t               => :Uint16,
    :uint8_t                => :Uint8,
    :int64_t                => :Int64,
    :int32_t                => :Int32,
    :int16_t                => :Int16,
    :int8_t                 => :Int8,
    :uintptr_t              => :Csize_t
    )

int_conversion = @compat Dict{Any,Any}(
    :Cint   => int32,
    :Cuint  => uint32,
    :Uint64 => uint64,
    :Uint32 => uint32,
    :Uint16 => uint16,
    :Uint8  => uint8,
    :Int64  => int64,
    :Int32  => int32,
    :Int16  => int16,
    :Int8   => int8
    )

################################################################################
#
# libclang objects to Julia representation
#
# each repr_jl function takes one or more CLCursor or CLType objects,
# and returns the appropriate representation.
#
################################################################################

function repr_jl(t::Union(cindex.Record, cindex.Typedef))
    tname = symbol(spelling(cindex.getTypeDeclaration(t)))
    return get(cl_to_jl, tname, tname)
end

function repr_jl(t::TypeRef)
    reftype = cindex.getCursorReferenced(t)
    refdef = cindex.getCursorDefinition(reftype)
    if isa(refdef, cindex.InvalidFile) ||
       isa(refdef, cindex.FirstInvalid) ||
       isa(refdef, cindex.Invalid)
        return :Void
    else
        return symbol(spelling(reftype))
    end
end

function repr_jl(ptr::cindex.Pointer)
    ptee = pointee_type(ptr)
    Expr(:curly, :Ptr, repr_jl(ptee))
end

function repr_jl(parm::cindex.ParmDecl)
    return repr_jl(cu_type(parm))
end

function repr_jl(unxp::Unexposed)
    r = spelling(cindex.getTypeDeclaration(unxp))
    r == "" ? :Void : symbol(r)
end

function repr_jl(t::ConstantArray)
    # For ConstantArray declarations, we make an immutable
    # array with that many members of the appropriate type.

    # Override pointer representation so pointee type is not written
    repr_short(t::CLType) = isa(t, Pointer) ? "Ptr" : repr_jl(t)

    # Grab the WrapContext in order to add bitstype
    global context::WrapContext

    arrsize = cindex.getArraySize(t)
    eltype = cindex.getArrayElementType(t)
    typename = string("Array_", arrsize, "_", repr_short(eltype))
    typesym = symbol(typename)

    if !(typesym in keys(context.common_buf))
        b = Expr(:block)
        e = Expr(:type, false, symbol(typename), b)
        repr = repr_jl(eltype)
        for i = 1:arrsize
            push!(b.args, Expr(:(::), symbol("d$i"), repr))
        end

        # Create a zero function for this type
        b = :($(symbol(typename))(fill(zero($repr), $arrsize)...))
        zero_call = :(zero(::Type{$(symbol(typename))}) = $b)

        context.common_buf[typesym] = ExprUnit(Any[e, zero_call])
    end

    return typesym
end

function repr_jl(t::IncompleteArray)
    eltype = cindex.getArrayElementType(t)
    Expr(:curly, :Ptr, repr_jl(eltype))
end

function repr_jl(t::UnionType)
    tdecl = cindex.getTypeDeclaration(t)
    maxelem = largestfield(tdecl)
    return repr_jl(maxelem)
end

function repr_jl(arg::cindex.CLType)
    rep = get(cl_to_jl, typeof(arg), nothing)
    rep == nothing && error("No CLType translation available for: ", arg)
    return rep
end

target_type(s::Symbol) = s
function target_type(e::Expr)
    if e.head == :curly && e.args[1] == :Ptr
        return target_type(e.args[2])
    else
        error("target_type: don't know how to handle $e")
    end
end
target_type(q) = error("target_type: don't know how to handle $q")

###############################################################################
# Get field decl sizes
#   - used for hacky union inclusion: we find the largest union field and
#     declare a block of bytes to match.
###############################################################################

typesize(t::CLType) = sizeof(eval(cl_to_jl[typeof(t)]))
typesize(t::Record) = begin warn("  incorrect typesize for Record field"); 0 end
typesize(t::Unexposed) = begin warn("  incorrect typesize for Unexposed field"); 0 end
typesize(t::ConstantArray) = cindex.getArraySize(t)
typesize(t::TypedefDecl) = typesize(cindex.getTypedefDeclUnderlyingType(t))
typesize(t::Typedef) = typesize(cindex.getTypeDeclaration(t))
typesize(t::Invalid) = begin warn("  incorrect typesize for Invalid field"); 0 end

function largestfield(cu::UnionDecl)
    maxsize,maxelem = 0,0
    fields = children(cu)
    for i in 1:length(fields)
        field_size = typesize(cu_type(fields[i]))
        if field_size > maxsize
            maxsize = field_size
            maxelem = i
        end
    end
    fields[maxelem]
end

function fieldsize(cu::FieldDecl)
    fieldsize(children(cu)[1])
end

################################################################################
# Handle function declarations
################################################################################

function wrap(context::WrapContext, buf::Array, func_decl::FunctionDecl, libname)
    func_type = cindex.cu_type(func_decl)
    if isa(func_type, FunctionNoProto)
        warn("No Prototype for $(func_decl) - assuming no arguments")
	elseif cindex.isFunctionTypeVariadic(func_type) == 1
        warn("Skipping VarArg Function $(func_decl)")
        return
    end

    funcname = symbol(spelling(func_decl))
    ret_type = repr_jl(return_type(func_decl))

    args = cindex.function_args(func_decl)

    arg_types = [cindex.getArgType(func_type, UInt32(i)) for i in 0:length(args)-1]
    arg_reps = [repr_jl(x) for x in arg_types]

    # check whether any argument types are blocked
    for arg_t in arg_types
        if spelling(arg_t) in reserved_argtypes
            warning("Skipping $(name(func_decl)) due to unsupported argument: $(name(arg_t))")
            return
        end
    end

    # Handle unnamed args and convert names to symbols
    arg_count = 0
    arg_names = convert(Vector{Symbol},
                        map(x-> symbol(begin
                            nm = name_safe(cindex.name(x))
                            nm != "" ? nm : "arg"*string(arg_count+=1)
                        end), args))

    sig = efunsig(funcname, arg_names, arg_reps)
    body = eccall(funcname, symbol(libname), ret_type, arg_names, arg_reps)
    e = Expr(:function, sig, Expr(:block, body))

    push!(buf, e)
end

################################################################################
# Handle all other declarations
################################################################################

function wrap(context::WrapContext, expr_buf::OrderedDict, cursor::EnumDecl; usename="")
    if (usename == "" && (usename = name(cursor)) == "")
        usename = name_anon()
    end
    enumname = usename

    buf = Any[]
    enum_exprs = ExprUnit(buf)
    expr_buf[symbol_safe(enumname)] = enum_exprs

    push!(buf, "# begin enum $enumname")
    enumtype = repr_jl(cindex.getEnumDeclIntegerType(cursor))
    _int = int_conversion[enumtype]
    push!(buf, :(typealias $(symbol(enumname)) $enumtype))
    for enumitem in children(cursor)
        cur_name = cindex.spelling(enumitem)
        if (length(cur_name) < 1) continue end
        cur_sym = symbol_safe(cur_name)
        push!(buf, :(const $cur_sym = $_int($(value(enumitem)))))
        expr_buf[cur_sym] = enum_exprs
    end
    push!(buf, "# end enum $enumname")

end

function wrap(context::WrapContext, expr_buf::OrderedDict, sd::StructDecl; usename = "")
    !context.options.wrap_structs && return

    if (usename == "" && (usename = name(sd)) == "")
        warn("Skipping unnamed StructDecl")
        return
    end
    usesym = symbol(usename)

    struct_fields = children(sd)

    # Generate type declaration
    b = Expr(:block)
    e = Expr(:type, !context.options.immutable_structs, usesym, b)
    deps = OrderedSet{Symbol}()
    for cu in struct_fields
        cur_name = spelling(cu)
        if (isa(cu, StructDecl) || isa(cu, UnionDecl))
            continue
        elseif !(isa(cu, FieldDecl) || isa(cu, TypeRef))
            expr_buf[usesym] = ExprUnit(Poisoned())
            warn("Skipping struct: \"$usename\" due to unsupported field: $cur_name")
            return
        elseif (length(cur_name) < 1)
            error("Unnamed struct member in: $usename ... cursor: ", string(cu))
        end

        repr = repr_jl(cu_type(cu))
        push!(b.args, Expr(:(::), symbol_safe(cur_name), repr))
        push!(deps, target_type(repr))
    end

    # Check for a previous forward ordering
    if !(usesym in keys(expr_buf)) || (expr_buf[usesym].state == :empty)
        if length(struct_fields) > 0
            expr_buf[usesym] = ExprUnit(e, deps)
        else
            # Possible forward definition
            expr_buf[usesym] = ExprUnit(e, deps, state=:empty)
        end
    end
    return
end

function wrap(context::WrapContext, expr_buf::OrderedDict, ud::UnionDecl; usename = "")
    if (usename == "" && (usename = name(ud)) == "")
        warn("Skipping unnamed UnionDecl")
        return
    end

    b = Expr(:block)
    e = Expr(:type, !context.options.immutable_structs, symbol(usename), b)
    max_cu = largestfield(ud)
    cur_sym = symbol("_"*usename)
    target = repr_jl(cu_type(max_cu))

	if string(target) == ""
        warn("Skipping UnionDecl $(usename) because largest field '$(name(max_cu))' could not be typed", string(target))
        return
	end

    push!(b.args, Expr(:(::), cur_sym, target))

    # TODO: add other dependencies
    expr_buf[cur_sym] = ExprUnit(e, Any[target])

    return
end

function efunsig(name::Symbol, args::Vector{Symbol}, types)
    x = Any[ Expr(:(::), a, t) for (a,t) in zip(args,types) ]
    Expr(:call, name, x...)
end

function eccall(funcname::Symbol, libname::Symbol, rtype, args, types)
    Expr(:ccall,
         Expr(:tuple, QuoteNode(funcname), libname),
         rtype,
         Expr(:tuple, types...),
         args...)
end

function wrap(context::WrapContext, expr_buf::OrderedDict, tdecl::TypedefDecl; usename="")
    cursor_type = cindex.cu_type(tdecl)
    td_type = cindex.getTypedefDeclUnderlyingType(tdecl)

    if isa(td_type, Unexposed)
        decl_parts = children(tdecl)

        local tdunxp

        for part in decl_parts
            # skip any leading non-type cursors
            # Attributes, ...
            if !isa(part, FirstAttr)
                tdunxp = part
                break
            end
        end

        if isa(tdunxp, TypeRef)
            td_type = tdunxp
        else
            wrap(context, expr_buf, tdunxp; usename=name(tdecl))
            return # TODO.. ugly flow
        end
    elseif isa(td_type, FunctionProto)
        return string("# Skipping Typedef: FunctionProto ", spelling(tdecl))
    end

    td_sym = symbol(spelling(tdecl))
    td_target = repr_jl(td_type)
    if !haskey(expr_buf, td_sym)
        expr_buf[td_sym] = ExprUnit(Expr(:typealias, td_sym, td_target), Any[td_target])
    end
end

################################################################################
# Handler for macro definitions
################################################################################
#
# For handling of #define'd constants, allows basic expressions
# but bails out quickly.

function handle_macro_exprn(tokens::TokenList, pos::Int)
    function trans(tok)
        ops = ["+" "-" "*" ">>" "<<" "/" "\\" "%" "|" "||" "^" "&" "&&"]
        if (isa(tok, cindex.Literal) ||
            (isa(tok,cindex.Identifier))) return 0
        elseif (isa(tok, cindex.Punctuation) && tok.text in ops) return 1
        else return -1
        end
    end

    # normalize literal with a size suffix
    function literally(tok)
        # note: put multi-character first, or it will break out too soon for those!
        literalsuffixes = ["ULL", "Ull", "uLL", "ull", "LLU", "LLu", "llU", "llu",
                           "LL", "ll", "UL", "Ul", "uL", "ul", "LU", "Lu", "lU", "lu",
                           "U", "u", "L", "l", "F", "f"]
        txt = tok.text
        if isa(tok,cindex.Identifier) || isa(tok,cindex.Punctuation)
            # pass
        elseif isa(tok,cindex.Literal)
            txt = strip(tok.text)
            for sfx in literalsuffixes
                if endswith(txt, sfx)
                    txt = txt[1:end-length(sfx)]
                    break
                end
            end
        end
        return txt
    end

    # check whether identifiers and literals alternate
    # with punctuation
    exprn = ""
    prev = 1 >> trans(tokens[pos])
    for pos = pos:length(tokens)
        tok = tokens[pos]
        state = trans(tok)
        if ( state $ prev  == 1)
            prev = state
        else
            break
        end
        exprn = exprn * literally(tok)
    end
    return (exprn,pos)
end

# TODO: This really returns many more symbols than we want,
# Functionally, it shouldn't matter, but eventually, we
# might want something more sophisticated.
# (Check: Does this functionality already exist elsewhere?)
get_symbols(s) = Any[]
get_symbols(s::Symbol) = Any[s]
get_symbols(e::Expr) = vcat(get_symbols(e.head), get_symbols(e.args))
get_symbols(xs::Array) = reduce(vcat, [get_symbols(x) for x in xs])

function wrap(context::WrapContext, expr_buf::OrderedDict, md::cindex.MacroDefinition)
    tokens = tokenize(md)
    # Skip any empty definitions
    tokens.size < 2 && return
    startswith(name(md), "_") && return

    pos = 1; exprn = ""
    if tokens[2].text == "("
        exprn,pos = handle_macro_exprn(tokens, 3)
        if (pos != endof(tokens) || tokens[pos].text != ")")
            mdef_str = join([c.text for c in tokens], " ")
            expr_buf[symbol(mdef_str)] = ExprUnit(string("# Skipping MacroDefinition: ", mdef_str))
            return
        end
        exprn = "(" * exprn * ")"
    else
        (exprn,pos) = handle_macro_exprn(tokens, 2)
        if pos != endof(tokens)
            mdef_str = join([c.text for c in tokens], " ")
            expr_buf[symbol(mdef_str)] = ExprUnit(string("# Skipping MacroDefinition: ", mdef_str))
            return
        end
    end

    # Occasionally, skipped definitions slip through
    (exprn == "" || exprn == "()") && return

    use_sym = symbol_safe(tokens[1].text)
    target = parse(exprn)
    e = Expr(:const, Expr(:(=), use_sym, target))

    deps = get_symbols(target)
    expr_buf[use_sym] = ExprUnit(e, deps)
end

# Does this actually occur in header files???
function wrap(context::WrapContext, expr_buf::OrderedDict, cursor::TypeRef; usename="")
    usename == "" && (usename = name(cursor))
    error("Found typeref: ", cursor)
    println("Printing typeref: ", cursor)
    usesym = symbol(usename)
    expr_buf[usesym] = ExprUnit(usesym)
end

function wrap(context::WrapContext, buf, cursor; usename="")
    warn("Not wrapping $(typeof(cursor))  $usename $(name(cursor))")
end

################################################################################
# Wrapping driver
################################################################################

# Any failed cursor in the loop below is put in this global for debugging
debug_cursors = CLCursor[]

function wrap_header(wc::WrapContext, topcu::CLCursor, top_hdr, obuf::Array)
    println("WRAPPING HEADER: $top_hdr")

    topcl = children(topcu)

    # Loop over all of the child cursors and wrap them, if appropriate.
    for i=1:topcl.size
        cursor = topcl[i]
        cursor_hdr = cu_file(cursor)
        cursor_name = name(cursor)

        # what should be wrapped:
        #   1. wrap if context.header_wrapped(top_header, cursor_header) == True
        #   2. wrap if context.cursor_wrapped(cursor_name, cursor) == True
        #   3. skip compiler defs and cursors already wrapped
        if !wc.header_wrapped(top_hdr, cursor_hdr) ||
           !wc.cursor_wrapped(cursor_name, cursor) ||       # client callbacks
           startswith(cursor_name, "__")           ||       # skip compiler definitions
           ((cursor_name in keys(wc.common_buf)) &&         # already wrapped
                    !(cursor_name in keys(wc.empty_structs)))
            continue
        end

        try
            if (isa(cursor, FunctionDecl))
                wrap(wc, obuf, cursor, wc.header_library(cu_file(cursor)))
            elseif !isa(cursor, TypeRef)
                wrap(wc, wc.common_buf, cursor)
            else
                continue
            end
        catch err
            push!(debug_cursors::Array{CLCursor,1}, cursor)
            info("Error thrown. Last cursor available in Clang.wrap_c.debug_cursors")
            rethrow(err)
        end
    end

    cindex.cl_dispose(topcl)
end

function parse_c_headers(wc::WrapContext)
    parsed = Dict{ASCIIString, CLCursor}()

    # Parse the headers
    for header in unique(wc.headers)
        topcu = cindex.parse_header(header;
                                    index = wc.index,
                                    args  = wc.clang_args,
                                    includes = wc.clang_includes,
                                    flags = TranslationUnit_Flags.DetailedPreprocessingRecord |
                                    TranslationUnit_Flags.SkipFunctionBodies)
        parsed[header] = topcu
    end
    return parsed
end

function sort_includes(wc::WrapContext, parsed)
    includes = mapreduce(x->search(parsed[x], InclusionDirective), append!, keys(parsed))
    header_paths = unique(map(x->cindex.getIncludedFile(x), includes))
    return unique([filter(x -> x in header_paths, wc.headers), wc.headers])
end

################################################################################
# Wrapping driver
################################################################################

# Pretty-print a buffer of expressions (and comments) to an output stream
# Adds blank lines at appropriate places for readability
function print_buffer(ostrm, obuf)

    state = :comment
    in_enum = false

    for e in obuf
        isa(e, Poisoned) && continue
        prev_state = state
        if state != :enum
            state = (isa(e, String) ? :string :
                     isa(e, Expr)   ? e.head :
                     error("output: can't handle type $(typeof(e))"))
        end

        if state == :string
            if startswith(e, "# begin enum")
                state = :enum
                println(ostrm)
            elseif startswith(e, "# Skipping")
                state = :skipping
            end
        end

        if state != :enum
            if ((state != prev_state && prev_state != :string) ||
                (state == prev_state && (state == :function ||
                                         state == :type)))
                println(ostrm)
            end
        end

        println(ostrm, e)

        if state == :enum && isa(e, String) && startswith(e, "# end enum")
            state = :end_enum
        end
    end
end

function dump_to_buf!(buf::Array, expr_buf::OrderedDict{Symbol, ExprUnit}, item::ExprUnit)
    (item.state == :done || item.state == :processing) && return
    item.state = :processing

    for dep in item.deps
        (dep_expr = get(expr_buf, dep, nothing)) == nothing && continue
        dump_to_buf!(buf, expr_buf, dep_expr)
    end

    append!(buf, item.items)
    item.state = :done
end

function dump_to_buf(expr_buf::OrderedDict{Symbol, ExprUnit})
    buf = Any[]
    for item in values(expr_buf)
        item.state == :done && continue
        dump_to_buf!(buf, expr_buf, item)
    end
    buf
end

function Base.run(wc::WrapContext)
    # Parse headers
    parsed = parse_c_headers(wc)
    # Sort includes by requirement order
    wc.headers = sort_includes(wc, parsed)

    # Helper to store file handles
    filehandles = Dict{ASCIIString,IOStream}()
    getfile(f) = (f in keys(filehandles)) ? filehandles[f] : (filehandles[f] = open(f, "w"))

    for hfile in wc.headers
        outfile = wc.header_outputfile(hfile)
        obuf = wc.output_bufs[hfile]

        # Extract header to Expr[] array
        wrap_header(wc, parsed[hfile], hfile, obuf)

        # Apply user-supplied transformation
        wc.output_bufs[hfile] = wc.rewriter(obuf)

        # Debug
        println("writing $(outfile)")

        # Write output
        ostrm = getfile(outfile)
        println(ostrm, "# Julia wrapper for header: $hfile")
        println(ostrm, "# Automatically generated using Clang.jl wrap_c, version $version\n")

        print_buffer(ostrm, wc.output_bufs[hfile])
    end

    common_buf = dump_to_buf(wc.common_buf)

    # Apply user-supplied transformation
    common_buf = wc.rewriter(common_buf)

    # Write "common" definitions: types, typealiases, etc.
    print_buffer(open(wc.common_file, "w"), common_buf)

    map(close, values(filehandles))
end

# Deprecated interface
@deprecate wrap_c_headers(wc::WrapContext, headers)   (wc.headers = headers; run(wc))

###############################################################################
# Utilities
###############################################################################

function name_anon()
    global context::WrapContext
    "ANONYMOUS_"*string(context.anon_count += 1)
end

function name_safe(cursor_name::String)
    return (cursor_name in reserved_words) ? "_"*cursor_name : cursor_name
end
symbol_safe(cursor_name::String) = symbol(name_safe(cursor_name))

###############################################################################

end # module wrap_c


text:
version = v"2.26.0"
include("../common.jl")

build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.7")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
using MinLSH
using Base.Test

numband = 5
bandwidth = 20
lsh_vars = LSH_vars(numband, bandwidth)
ret_prob_table(numband, bandwidth)

@testset "member" begin
@test lsh_vars.numband   == numband
@test lsh_vars.bandwidth == bandwidth
@test lsh_vars.numhash   == numband * bandwidth
end

# generate_signature_matrixで作成したsignature_matrixが正しくJaccard係数近似になっているかどうかのテスト
a = zeros(Bool, 20, 10)
# a[1:5, :] = 1
a[2:7, 1] = 1
a[2:8, 2] = 1
a[2:5, 3] = 1
a[7:8, 4] = 1
a[:, 6] = 1
a[:, 7] = 1
jaccard(set1::Set{UInt}, set2::Set{UInt}) = length(set1 ∩ set2) / length(set1 ∪ set2)
# sbinvec: sparse (CSC) binary vector
function jaccard(sbinvec1::SparseVector{Bool, Int64}, sbinvec2::SparseVector{Bool, Int64})
    set1 = Set{UInt}(findn(sbinvec1 .== 1))
    set2 = Set{UInt}(findn(sbinvec2 .== 1))
    return(jaccard(set1, set2))
end

b = sparse(a)
sigmat = generate_signature_matrix(lsh_vars, b, window_num = 1, slidewidth = 1)
ethreshould = 3e-2
print_with_color(:red, "estimation error within "*string(ethreshould)*" are allowed \n")
testcols = [(1, 2), (2, 3), (2, 4)]
for (col1, col2) in testcols
    count = 0
    for row in 1:size(sigmat)[1]
        if sigmat[row, col1] == sigmat[row, col2]
            count += 1
        end
    end
    estimated = count / size(sigmat)[1]
    @show jaccard(b[:, col1], b[:, col2]) - estimated
    @test_approx_eq_eps jaccard(b[:, col1], b[:, col2]) estimated ethreshould
end

max_stars_repo_name: KeitaW/MinLSH

max_issues_repo_name: KeitaW/MinLSH


text:
"""
    heatmap(A; kw...)

# Description

Draws a heatmap for the given points.
It uses the `Matrix` `A` as the values of the heatmap, with the column and row indices of the matrix
as x and y coordinates respectively.

# Usage

    heatmap(A::AbstractMatrix; $(keywords((width = 0, height = 0, xfact = 0, yfact = 0); add = (Z_DESCRIPTION..., :fix_ar)))

# Arguments

$(arguments(
    (
        A = "input matrix (color values)",
        xfact = "scale for the `x` coordinate (defaults to 0 - i.e. each column in `A` maps to one unit, `x` origin starting at 1). If set to anything else, the x origin will start at 0",
        yfact = "scale for the `y` coordinate labels (defaults to 0 - i.e. each row in `A` maps to one unit, `y` origin starting at 1). If set to anything else, the y origin will start at 0",
        xoffset = "plotting offset for the `x` coordinate (after scaling)",
        yoffset = "plotting offset for the `y` coordinate (after scaling)",
    ); add = (Z_DESCRIPTION..., :fix_ar)
))

# Returns

A plot object of type `Plot{HeatmapCanvas}`.

# Author(s)

- Rowan Katekar (github.com/rjkat)

# See also

`Plot`, `scatterplot`, `HeatmapCanvas`
"""
function heatmap(
    A::AbstractMatrix;
    xlim = KEYWORDS.xlim,
    ylim = KEYWORDS.ylim,
    zlim = KEYWORDS.zlim,
    xoffset = 0.0,
    yoffset = 0.0,
    out_stream::Union{Nothing,IO} = nothing,
    width::Int = 0,
    height::Int = 0,
    margin::Int = KEYWORDS.margin,
    padding::Int = KEYWORDS.padding,
    colormap = KEYWORDS.colormap,
    xfact = 0,
    yfact = 0,
    labels = true,
    fix_ar::Bool = false,
    kw...,
)
    nrows = size(A, 1)
    ncols = size(A, 2)

    # if scale is auto, use the matrix indices as axis labels
    # otherwise, start axis labels at zero
    X = xfact == 0 ? collect(1:ncols) : collect(0:(ncols - 1)) .* xfact
    X .+= xoffset
    xfact == 0 && (xfact = 1)
    Y = yfact == 0 ? collect(1:nrows) : collect(0:(nrows - 1)) .* yfact
    Y .+= yoffset
    yfact == 0 && (yfact = 1)

    # set the axis limits automatically
    if xlim == (0, 0) && length(X) > 0
        xlim = extrema(X)
    end
    if ylim == (0, 0) && length(Y) > 0
        ylim = extrema(Y)
    end

    # select a subset of A based on the supplied limits
    firstx = findfirst(x -> x >= xlim[1], X)
    lastx = findlast(x -> x <= xlim[2], X)
    xrange = (firstx == nothing || lastx == nothing) ? (1:0) : (firstx:lastx)
    firsty = findfirst(y -> y >= ylim[1], Y)
    lasty = findlast(y -> y <= ylim[2], Y)
    yrange = (firsty == nothing || lasty == nothing) ? (1:0) : (firsty:lasty)
    A = A[yrange, xrange]
    X = X[xrange]
    Y = Y[yrange]

    # allow A to be an array over which min and max is not defined,
    # e.g. an array of RGB color values
    minz, maxz = (0, 0)
    noextrema = true
    try
        minz, maxz = extrema(A)
        noextrema = false
    catch
    end
    if zlim != (0, 0)
        noextrema &&
            throw(ArgumentError("zlim cannot be set when the element type is $(eltype(A))"))
        minz, maxz = zlim
    end

    # if A is an rgb image, translate the colors directly to the terminal
    callback = if length(A) > 0 && all(x -> x in fieldnames(eltype(A)), [:r, :g, :b])
        (A, minz, maxz) -> rgbimgcolor(A)
    else
        colormap_callback(colormap)
    end

    nrows = ceil(Int, (ylim[2] - ylim[1]) / yfact) + 1
    ncols = ceil(Int, (xlim[2] - xlim[1]) / xfact) + 1
    data_ar = ncols / nrows  # data aspect ratio

    max_width = width == 0 ? (height == 0 ? 0 : ceil(Int, height * data_ar)) : width
    max_height = height == 0 ? (width == 0 ? 0 : ceil(Int, width / data_ar)) : height

    # 2nrows: compensate nrows(c::HeatmapCanvas) = div(size(grid(c), 2) + 1, 2)
    width, height, max_width, max_height = get_canvas_dimensions_for_matrix(
        HeatmapCanvas,
        2nrows,
        ncols,
        max_width,
        max_height,
        width,
        height,
        margin,
        padding,
        out_stream,
        fix_ar,
    )

    colorbar = if height < 7
        # for small plots, don't show colorbar by default
        !noextrema && get(kw, :colorbar, false)
    else
        # show colorbar by default, unless set to false, or labels == false
        !noextrema && get(kw, :colorbar, labels)
    end
    kw = (; kw..., colorbar = colorbar)

    xs = length(X) > 0 ? [X[1], X[end]] : Float64[0, 0]
    ys = length(Y) > 0 ? [Y[1], Y[end]] : Float64[0, 0]
    plot = Plot(
        xs,
        ys,
        HeatmapCanvas;
        grid = false,
        colormap = callback,
        colorbar = colorbar,
        colorbar_lim = (minz, maxz),
        ylim = ylim,
        xlim = xlim,
        labels = labels,
        width = width,
        height = height,
        min_width = 1,
        min_height = 1,
        kw...,
    )

    for row in 1:length(Y)
        points!(
            plot,
            X,
            fill(Y[row], length(X)),
            UserColorType[callback(v, minz, maxz) for v in A[row, :]],
        )
    end
    plot
end


text:
include("debug.jl")
include("fpump.jl")

function create_root_model!(optimizer::MOI.AbstractOptimizer, jp::JuniperProblem; fix_start=false)
    ps = jp.options.log_levels

    jp.model = Model()
    lb = jp.l_var
    ub = jp.u_var
    # all continuous we solve relaxation first
    @variable(jp.model, lb[i] <= x[i=1:jp.num_var] <= ub[i])

    for i=1:jp.num_var
        JuMP.set_start_value(x[i], jp.primal_start[i])
    end

    register_functions!(jp.model, jp.options.registered_functions)

    # TODO check whether it is supported
    if optimizer.nlp_data.has_objective
        obj_expr = MOI.objective_expr(optimizer.nlp_data.evaluator)
        obj_expr = expr_dereferencing(obj_expr, jp.model)
        try
            JuMP.set_NL_objective(jp.model, optimizer.sense, obj_expr)
        catch
            error("Have you registered a function? Then please register the function also for Juniper see: \n
            https://lanl-ansi.github.io/Juniper.jl/stable/options/#registered_functions%3A%3AUnion%7BNothing%2CVector%7BRegisteredFunction%7D%7D-%5Bnothing%5D-1")
        end
    elseif optimizer.objective !== nothing
        MOI.set(jp.model, MOI.ObjectiveFunction{typeof(optimizer.objective)}(), optimizer.objective)
        MOI.set(jp.model, MOI.ObjectiveSense(), optimizer.sense)
    end

    backend = JuMP.backend(jp.model);
    llc = optimizer.linear_le_constraints
    lgc = optimizer.linear_ge_constraints
    lec = optimizer.linear_eq_constraints
    qlc = optimizer.quadratic_le_constraints
    qgc = optimizer.quadratic_ge_constraints
    qec = optimizer.quadratic_eq_constraints
    for constr_type in [llc, lgc, lec, qlc, qgc, qec]
        for constr in constr_type
            MOI.add_constraint(backend, constr[1], constr[2])
        end
    end
    for i in 1:jp.num_nl_constr
        constr_expr = MOI.constraint_expr(optimizer.nlp_data.evaluator, i)
        constr_expr = expr_dereferencing(constr_expr, jp.model)
        try
            JuMP.add_NL_constraint(jp.model, constr_expr)
        catch
            error("Have you registered a function? Then please register the function also for Juniper see: \n
            https://lanl-ansi.github.io/Juniper.jl/stable/options/#registered_functions%3A%3AUnion%7BNothing%2CVector%7BRegisteredFunction%7D%7D-%5Bnothing%5D-1")
        end
    end

    jp.x = x
end

function fix_primal_start!(jp::JuniperProblem)
    lb = jp.l_var
    ub = jp.u_var
    x = jp.x
    for i=1:jp.num_var
        if jp.var_type[i] != :Cont && lb[i] <= jp.primal_start[i] <= ub[i]
            JuMP.set_lower_bound(jp.x[i], jp.primal_start[i])
            JuMP.set_upper_bound(jp.x[i], jp.primal_start[i])
        else
            JuMP.set_start_value(x[i], jp.primal_start[i])
        end
    end
end

function unfix_primal_start!(jp::JuniperProblem)
    lb = jp.l_var
    ub = jp.u_var
    x = jp.x
    for i=1:jp.num_var
        JuMP.set_lower_bound(jp.x[i], lb[i])
        JuMP.set_upper_bound(jp.x[i], ub[i])
        JuMP.set_start_value(x[i], jp.primal_start[i])
    end
end

function solve_root_incumbent_model(jp::JuniperProblem)
    status, backend = optimize_get_status_backend(jp.model; solver=jp.nl_solver)
    incumbent = nothing
    ps = jp.options.log_levels
    if state_is_optimal(status; allow_almost=jp.options.allow_almost_solved)
        # set incumbent
        objval = MOI.get(backend, MOI.ObjectiveValue())
        solution = JuMP.value.(jp.x)
        if are_type_correct(solution, jp.var_type, jp.disc2var_idx, jp.options.atol)
            if only_almost_solved(status) && jp.options.allow_almost_solved_integral
                @warn "Start value incumbent only almost locally solved. Disable with `allow_almost_solved_integral=false`"
            end

            incumbent = Incumbent(objval, solution, only_almost_solved(status))
            (:All in ps || :Info in ps) && println("Incumbent using start values: $objval")
        end
    else
        (:All in ps || :Info in ps) && println("Start values are not feasible.")
    end
    return incumbent
end

function solve_root_model!(jp::JuniperProblem)
    status, backend = optimize_get_status_backend(jp.model; solver=jp.nl_solver)
    jp.relaxation_status = status
    restarts = 0
    max_restarts = jp.options.num_resolve_root_relaxation
    jp.options.debug && debug_init(jp.debugDict)
    while !state_is_optimal(jp.relaxation_status; allow_almost=jp.options.allow_almost_solved) &&
        restarts < max_restarts && time()-jp.start_time < jp.options.time_limit

        # TODO freemodel for Knitro
        restart_values = generate_random_restart(jp)
        jp.options.debug && debug_restart_values(jp.debugDict,restart_values)
        for i=1:jp.num_var
            JuMP.set_start_value(jp.x[i], restart_values[i])
        end
        JuMP.optimize!(jp.model)
        jp.relaxation_status = MOI.get(backend, MOI.TerminationStatus())
        restarts += 1
    end
    if only_almost_solved(jp.relaxation_status)
        @warn "The relaxation is only almost solved."
    end
    return restarts
end


text:
@testset "requests.jl" begin
    @test FHIRClient._add_trailing_slash(HTTP.URI("https://example.com")) == HTTP.URI("https://example.com/")
    @test FHIRClient._add_trailing_slash(HTTP.URI("https://example.com/")) == HTTP.URI("https://example.com/")
end


text:
export Monomials,
    Monomial

"A basis of the monomials `x^i`."
struct Monomials{T} <: PolynomialBasis{T}
    n   ::  Int     # the degrees go from 0 to n-1
end

Monomials(n) = Monomials{Float64}(n)

show(io::IO, b::Monomials{Float64}) = print(io, "Monomials($(length(b)))")
show(io::IO, b::Monomials{T}) where T = print(io, "Monomials{$(T)}($(length(b)))")

support(dict::Monomials{T}) where {T} = DomainSets.FullSpace{T}()

size(dict::Monomials) = (dict.n,)

unsafe_eval_element(b::Monomials, idxn::PolynomialDegree, x) = x^degree(idxn)

function unsafe_eval_element_derivative(b::Monomials{T}, idxn::PolynomialDegree, x, order) where {T}
    @assert order > 0
    i = degree(idxn)
    if order > i
        zero(T)
    elseif order == 1
        i*x^(i-1)
    else
        factorial(i) / factorial(i-order) * x^(i-order)
    end
end

similar(b::Monomials, ::Type{T}, n::Int) where {T} = Monomials{T}(n)

extension(::Type{T}, src::Monomials, dest::Monomials; options...) where {T} =
    IndexExtension{T}(src, dest, 1:length(src))

restriction(::Type{T}, src::Monomials, dest::Monomials; options...) where {T} =
    IndexRestriction{T}(src, dest, 1:length(dest))

"A monomial basis function `x^i`."
struct Monomial{T} <: Polynomial{T}
    degree  ::  Int
end

Monomial{T}(p::Monomial) where {T} = Monomial{T}(p.degree)

show(io::IO, p::Monomial) = print(io, "x^$(degree(p)) (monomial)")

convert(::Type{TypedFunction{T,T}}, p::Monomial) where {T} = Monomial{T}(p.degree)

support(::Monomial{T}) where {T} = DomainSets.FullSpace{T}()

(m::Monomial)(x) = x^degree(m)

(*)(p1::Monomial, p2::Monomial) = (*)(promote(p1,p2)...)

(*)(p1::Monomial{T}, p2::Monomial{T}) where {T} = Monomial{T}(degree(p1)+degree(p2))

basisfunction(dict::Monomials, idx) = basisfunction(dict, native_index(dict, idx))
basisfunction(dict::Monomials{T}, idx::PolynomialDegree) where {T} = Monomial{T}(degree(idx))

dictionary(p::Monomial{T}) where {T} = Monomials{T}(degree(p)+1)


text:
"""
    join(r1::Relation, r2::Relation)

Performs the natural join between two `Relation`s. If two columns share the same name, they
must match on the output.
"""
function join(r1::Relation, r2::Relation)
    joinlabels = labels(r1) ∩ labels(r2)
    _join(r1, r2, indexes(r1, joinlabels), indexes(r2, joinlabels))
end

#const ⨝ = join


text:
using Test

tests = [ "test_configuration.jl", "test_PlaneWaves.jl", "test_UniformElectronGas.jl", "test_Updates.jl", "test_CPIMC.jl" ]

for t in tests
    @testset "$t" begin
        include(t)
    end
end


text:
# Truncated normal distribution
"""
    TruncatedNormal(mu, sigma, l, u)

The *truncated normal distribution* is a particularly important one in the family of truncated distributions.
We provide additional support for this type with `TruncatedNormal` which calls `Truncated(Normal(mu, sigma), l, u)`.
Unlike the general case, truncated normal distributions support `mean`, `mode`, `modes`, `var`, `std`, and `entropy`.
"""
TruncatedNormal(mu::Real, sigma::Real, a::Real, b::Real) =
    Truncated(Normal(mu, sigma), a, b)

### statistics

minimum(d::Truncated{Normal{T},Continuous}) where {T <: Real} = d.lower
maximum(d::Truncated{Normal{T},Continuous}) where {T <: Real} = d.upper

function mode(d::Truncated{Normal{T},Continuous}) where T <: Real
    μ = mean(d.untruncated)
    d.upper < μ ? d.upper :
    d.lower > μ ? d.lower : μ
end

modes(d::Truncated{Normal{T},Continuous}) where {T <: Real} = [mode(d)]

# do not export. Used in mean, var
function _F1(x::Real, y::Real; thresh=1e-7)
    @assert 0 < thresh < Inf
    -Inf < x < Inf && -Inf < y < Inf || throw(DomainError())
    ϵ = exp(x^2 - y^2)
    if abs(x) > abs(y)
        _F1(y,x)
    elseif abs(x - y) ≤ thresh
        δ = y - x
        √π*x + (√π/2 + (-√π*x/6 + (-√π/12 + x*(√π/90 + (√π*x^2)/90)δ)δ)δ)δ
    elseif max(x,y) < 0
        (1 - ϵ) / (ϵ * erfcx(-y) - erfcx(-x))
    elseif min(x,y) > 0
        (1 - ϵ) / (erfcx(x) - ϵ * erfcx(y))
    else
        exp(-x^2) * (1 - ϵ) / (erf(y) - erf(x))
    end
end

# do not export. Used in mean, var
function _F2(x::Real, y::Real; thresh=1e-7)
    @assert 0 < thresh < Inf
    -Inf < x < Inf && -Inf < y < Inf || throw(DomainError())
    ϵ = exp(x^2 - y^2)
    if abs(x) > abs(y)
        _F2(y,x)
    elseif abs(x - y) ≤ thresh
        δ = y - x
        √π*x^2 - √π/2 + (√π*x + (√π/3 - √π*x^2/3 + (((√π/30 + √π*x^2/45)x^2 - 4*√π/45)δ - √π*x/3)δ)δ)δ
    elseif max(x,y) < 0
        (x - ϵ * y) / (ϵ * erfcx(-y) - erfcx(-x))
    elseif min(x,y) > 0
        (x - ϵ * y) / (erfcx(x) - ϵ * erfcx(y))
    else
        exp(-x^2) * (x - ϵ * y) / (erf(y) - erf(x))
    end
end

# do not export. Used in mean
function _tnmean(a, b)
    -Inf < a ≤ b < Inf || throw(DomainError())
    √(2/π) * _F1(a/√2, b/√2)
end

# do not export. Used in mean
function _tnmean(a, b, μ, σ)
    -Inf < a ≤ b < Inf || throw(DomainError())
    -Inf < μ < Inf && 0 < σ < Inf || throw(DomainError())
    α = (a - μ)/σ; β = (b - μ)/σ
    μ + _tnmean(α, β) * σ
end

# do not export. Used in var
function _tnvar(a, b)
    -Inf < a ≤ b < Inf || throw(DomainError())
    1 + 2/√π * _F2(a/√2, b/√2) - 2/π * _F1(a/√2, b/√2)^2
end

# do not export. Used in var
function _tnvar(a, b, μ, σ)
    -Inf < a ≤ b < Inf || throw(DomainError())
    -Inf < μ < Inf && 0 < σ < Inf || throw(DomainError())
    α = (a - μ)/σ; β = (b - μ)/σ
    _tnvar(α, β) * σ^2
end

function mean(d::Truncated{Normal{T},Continuous}) where T <: Real
    d0 = d.untruncated
    μ = mean(d0)
    σ = std(d0)
    if isfinite(d.lower) && isfinite(d.upper) && isfinite(μ) && isfinite(σ)
        # avoids loss of significance when truncation is far from μ.
        # See https://github.com/cossio/TruncatedNormal.jl/blob/master/notes/normal.pdf.
        _tnmean(d.lower, d.upper, μ, σ)
    else
        a = (d.lower - μ) / σ
        b = (d.upper - μ) / σ
        μ + ((normpdf(a) - normpdf(b)) / d.tp) * σ
    end
end

function var(d::Truncated{Normal{T},Continuous}) where T <: Real
    d0 = d.untruncated
    μ = mean(d0)
    σ = std(d0)
    if isfinite(d.lower) && isfinite(d.upper) && isfinite(μ) && isfinite(σ)
        # avoids loss of significance when truncation is far from μ.
        # See https://github.com/cossio/TruncatedNormal.jl/blob/master/notes/normal.pdf.
        _tnvar(d.lower, d.upper, μ, σ)
    else
        a = (d.lower - μ) / σ
        b = (d.upper - μ) / σ
        z = d.tp
        φa = normpdf(a)
        φb = normpdf(b)
        aφa = isinf(a) ? 0.0 : a * φa
        bφb = isinf(b) ? 0.0 : b * φb
        t1 = (aφa - bφb) / z
        t2 = abs2((φa - φb) / z)
        abs2(σ) * (1 + t1 - t2)
    end
end

function entropy(d::Truncated{Normal{T},Continuous}) where T <: Real
    d0 = d.untruncated
    z = d.tp
    μ = mean(d0)
    σ = std(d0)
    a = (d.lower - μ) / σ
    b = (d.upper - μ) / σ
    aφa = isinf(a) ? 0.0 : a * normpdf(a)
    bφb = isinf(b) ? 0.0 : b * normpdf(b)
    0.5 * (log2π + 1.) + log(σ * z) + (aφa - bφb) / (2.0 * z)
end

### sampling

## Use specialized sampler, as quantile-based method is inaccurate in
## tail regions of the Normal, issue #343

function rand(rng::AbstractRNG, d::Truncated{Normal{T},Continuous}) where T <: Real
    d0 = d.untruncated
    μ = mean(d0)
    σ = std(d0)
    a = (d.lower - μ) / σ
    b = (d.upper - μ) / σ
    z = randnt(rng, a, b, d.tp)
    return μ + σ * z
end

# Rejection sampler based on algorithm from Robert (1995)
#
#  - Available at http://arxiv.org/abs/0907.4010

function randnt(rng::AbstractRNG, lb::Float64, ub::Float64, tp::Float64)
    local r::Float64
    if tp > 0.3   # has considerable chance of falling in [lb, ub]
        r = randn(rng)
        while r < lb || r > ub
            r = randn(rng)
        end
        return r

    else
        span = ub - lb
        if lb > 0 && span > 2.0 / (lb + sqrt(lb^2 + 4.0)) * exp((lb^2 - lb * sqrt(lb^2 + 4.0)) / 4.0)
            a = (lb + sqrt(lb^2 + 4.0))/2.0
            while true
                r = rand(rng, Exponential(1.0 / a)) + lb
                u = rand(rng)
                if u < exp(-0.5 * (r - a)^2) && r < ub
                    return r
                end
            end
        elseif ub < 0 && ub - lb > 2.0 / (-ub + sqrt(ub^2 + 4.0)) * exp((ub^2 + ub * sqrt(ub^2 + 4.0)) / 4.0)
            a = (-ub + sqrt(ub^2 + 4.0)) / 2.0
            while true
                r = rand(rng, Exponential(1.0 / a)) - ub
                u = rand(rng)
                if u < exp(-0.5 * (r - a)^2) && r < -lb
                    return -r
                end
            end
        else
            while true
                r = lb + rand(rng) * (ub - lb)
                u = rand(rng)
                if lb > 0
                    rho = exp((lb^2 - r^2) * 0.5)
                elseif ub < 0
                    rho = exp((ub^2 - r^2) * 0.5)
                else
                    rho = exp(-r^2 * 0.5)
                end
                if u < rho
                    return r
                end
            end
        end
    end
end


text:
using PyCall
import Base.size

global const gym = PyCall.pywrap(PyCall.pyimport("gym"))

mutable struct GymState <: AbsState; data; done; end
mutable struct GymAction <: AbsAction; action; end

size(s::GymState) = size(s.data)

==(lhs::GymAction, rhs::GymAction) = lhs.action == rhs.action
isequal(lhs::GymAction, rhs::GymAction) = lhs.action == rhs.action
hash(a::GymAction) = hash(a.action)

struct Spec
	id
	nondeterministic
	reward_threshold
	tags
	timestep_limit
	trials
end

struct GymEnv <: AbsEnvironment
	env
	actions
	spec
end

function GymEnv(name::AbstractString)
	env = gym.make(name)
	actions = nothing

	if :n in keys(env[:action_space])
		actions = map(GymAction, 0:(env[:action_space][:n]-1))
	elseif :low in keys(env[:action_space])
		actions = (env[:action_space][:low], env[:action_space][:high])
	else#algorithmic
		t = env[:action_space][:shape]
        actions = []
        for i=0:(t[1][1]-1)
            for j=0:(t[2][1]-1)
                for k=0:(t[3][1]-1)
                    push!(actions, GymAction((i, j, k)))
                end
            end
        end
	end

	s = Spec(env[:spec][:id], env[:spec][:nondeterministic], env[:spec][:reward_threshold],
		env[:spec][:tags], env[:spec][:timestep_limit], env[:spec][:trials])
	GymEnv(env, actions, s)
end

getActions(s::GymState, env::GymEnv) = env.actions
getInitialState(env::GymEnv) = GymState(env.env[:reset](), false)
isTerminal(state::GymState, env::GymEnv) = state.done

function transfer(env::GymEnv, state::GymState, action::GymAction)
	obs, reward, done, info = env.env[:step](action.action)
	return (GymState(obs, done), reward)
end

monitor_start(env::GymEnv, fname::AbstractString) = env.env[:monitor][:start](fname)
monitor_close(env::GymEnv) = env.env[:monitor][:close]()
render_env(env::GymEnv) = env.env[:render]()
sample(env::GymEnv) = GymAction(env.env[:action_space][:sample]())

max_stars_repo_name: ozanarkancan/ReinforcementLearning

max_issues_repo_name: ozanarkancan/ReinforcementLearning

max_forks_repo_name: ozanarkancan/ReinforcementLearning