text:
"""
**MultinomialLogistic**

Multiclass likelihood with Multinomial logistic link : ``p(y=i|{fₖ}) = \\prod_{i\\neq k} \\Theta(f_i-f_k) ``
"""
struct MultinomialLogisticLikelihood{T<:Real} <: MultinomialLogisticLikelihood{T}
    Y::AbstractVector{BitVector} #Mapping from instances to classes
    class_mapping::AbstractVector{Any} # Classes labels mapping
    ind_mapping::Dict{Any,Int} # Mapping from label to index
    y_class::AbstractVector{Int64} # GP Index for each sample
    c::AbstractVector{AbstractVector{T}} # Second moment of fₖ
    θ::AbstractVector{AbstractVector{T}} # Variational parameter of Polya-Gamma distribution
    function MultinomialLogisticLikelihood{T}() where {T<:Real}
        new{T}()
    end
    function MultinomialLogisticLikelihood{T}(Y::AbstractVector{<:BitVector},
    class_mapping::AbstractVector, ind_mapping::Dict{<:Any,<:Int},y_class::AbstractVector{<:Int}) where {T<:Real}
        new{T}(Y,class_mapping,ind_mapping,y_class)
    end
    function MultinomialLogisticLikelihood{T}(Y::AbstractVector{<:BitVector},
    class_mapping::AbstractVector, ind_mapping::Dict{<:Any,<:Int},y_class::AbstractVector{<:Int},
    c::AbstractVector{<:AbstractVector{<:Real}}, θ::AbstractVector{<:AbstractVector}) where {T<:Real}
        new{T}(Y,class_mapping,ind_mapping,y_class,c,θ)
    end
end

function MultinomialLogisticLikelihood()
    MultinomialLogisticLikelihood{Float64}()
end

function pdf(l::MultinomialLogisticLikelihood,f::AbstractVector)
    multinomiallogistic(f)
end

function pdf(l::MultinomialLogisticLikelihood,y::Integer,f::AbstractVector)
    multinomiallogistic(f)[y]
end

function multinomiallogistic(f::AbstractVector)
    y = zero(f)
    for (i,x) in enumerate(f)
        y[i]=prod(logistic.(x.-f[1:length(f).!=i]))
    end
    return y
end

function multinomialheaviside(f::AbstractVector)
    y = zero(f)
    for (i,x) in enumerate(f)
        y[i]=prod(x.>f[1:length(f).!=i])
    end
    return y
end

function multinomialprobit(f::AbstractVector)
    y = zero(f)
    for (i,x) in enumerate(f)
        y[i]=prod(cdf.(Normal(0,1),x.-f[1:length(f).!=i]))
    end
    return y
end

function Base.show(io::IO,model::AugmentedMultinomialLogisticLikelihood{T}) where T
    print(io,"Augmented Multinomial Logistic likelihood")
end

function init_likelihood(likelihood::AugmentedMultinomialLogisticLikelihood{T},nLatent::Integer,nSamplesUsed::Integer) where T
    c = [ones(T,nSamplesUsed) for i in 1:nLatent]
    θ = [abs.(rand(T,nSamplesUsed))*2 for i in 1:nLatent]
    AugmentedMultinomialLogisticLikelihood{T}(likelihood.Y,likelihood.class_mapping,likelihood.ind_mapping,likelihood.y_class,c,θ)
end

struct MultinomialLogisticLikelihood{T<:Real} <: AbstractMultinomialLogisticLikelihood{T}
    Y::AbstractVector{BitVector} #Mapping from instances to classes
    class_mapping::AbstractVector{Any} # Classes labels mapping
    ind_mapping::Dict{Any,Int} # Mapping from label to index
    y_class::AbstractVector{Int64} #GP Index for each sample
    function MultinomialLogisticLikelihood{T}() where {T<:Real}
        new{T}()
    end
    function MultinomialLogisticLikelihood{T}(Y::AbstractVector{<:BitVector},
    class_mapping::AbstractVector, ind_mapping::Dict{<:Any,<:Int},y_class::AbstractVector{<:Int}) where {T<:Real}
        new{T}(Y,class_mapping,ind_mapping,y_class)
    end
end

function MultinomialLogisticLikelihood()
    MultinomialLogisticLikelihood{Float64}()
end

function Base.show(io::IO,model::MultinomialLogisticLikelihood{T}) where T
    print(io,"logistic-softmax likelihood")
end

function init_likelihood(likelihood::MultinomialLogisticLikelihood{T},nLatent::Integer,nSamplesUsed::Integer) where T
    return likelihood
end

function local_updates!(model::VGP{AugmentedMultinomialLogisticLikelihood{T},AnalyticVI{T},T,V}) where {T<:Real,V<:AbstractVector{T}}
    model.likelihood.c .= broadcast((Σ::V,μ::V)->Σ.+abs2.(μ),diag.(model.Σ),model.μ)
    c = (x->sqrt.(x)).(model.likelihood.Y.*model.likelihood.c .+ model.likelihood.c)
    model.likelihood.θ .= broadcast((c::V)->1.0./c.*tanh.(0.5*c),model.likelihood.c)
end

function natural_gradient!(model::VGP{AugmentedMultinomialLogisticLikelihood{T},AnalyticVI{T},T,V}) where {T<:Real,V<:AbstractVector{T}}
    up_y = 0.5*((model.nLatent-1).+2*dot(!y,μ.*θ))
    up_f = 0.5*(1+2*μ[model.likelihood.class_y[i]]*θ[i])
    upsig_y = -0.5*(Diagonal([dot(!getindexmodel(model.likelihood.Y,[i]),getindexmodel(model.likelihood.θ,[i])) for i in 1:model.nSamples]) + invKmm) - model.η₂
    upsig_g = θ + invKmm
end

function local_updates!(model::SVGP{AugmentedMultinomialLogisticLikelihood{T},AnalyticVI{T},T,V}) where {T<:Real,V<:AbstractVector{T}}
    model.likelihood.c .= broadcast((μ::V,Σ::Symmetric{T,Matrix{T}},κ::Matrix{T},K̃::V)->sqrt.(K̃+opt_diag(κ*Σ,κ)+abs2.(κ*μ)),
                                    model.μ,model.Σ,model.κ,model.K̃)
    for _ in 1:5
        model.likelihood.γ .= broadcast((c::V,κμ::V,ψα::V)->(0.5/(model.likelihood.β[1]))*exp.(ψα).*safe_expcosh.(-0.5*κμ,0.5*c),
                                    model.likelihood.c,model.κ.*model.μ,[digamma.(model.likelihood.α)])
        model.likelihood.α .= 1.0.+(model.likelihood.γ...,)
    end
    model.likelihood.θ .= broadcast((y::BitVector,γ::V,c::V)->0.5.*(y[model.inference.MBIndices]+γ)./c.*tanh.(0.5.*c),
                                    model.likelihood.Y,model.likelihood.γ,model.likelihood.c)
    return nothing
end

function sample_local!(model::VGP{<:AugmentedMultinomialLogisticLikelihood,<:GibbsSampling})
    if model.inference.nIter <= 1
        # model.likelihood.α .= 10.0.*model.likelihood.α./model.likelihood.β
    end
    model.likelihood.γ .= broadcast(μ::AbstractVector{<:Real}->rand.(Poisson.(0.5*model.likelihood.α.*safe_expcosh.(-0.5*μ,0.5*μ))), model.μ)
    model.likelihood.α .= rand.(Gamma.(1.0.+(model.likelihood.γ...,),1.0./model.likelihood.β))
    model.likelihood.θ .= broadcast((y::BitVector,γ::AbstractVector{<:Real},μ::AbstractVector{<:Real})->PolyaGammaDist().draw.(y.+γ,μ),model.likelihood.Y,model.likelihood.γ,model.μ)
    return nothing
end

function sample_global!(model::VGP{<:AugmentedMultinomialLogisticLikelihood,<:GibbsSampling})
    model.Σ .= inv.(Symmetric.(Diagonal.(model.likelihood.θ).+model.invKnn))
    model.μ .= rand.(MvNormal.(0.5.*model.Σ.*(model.likelihood.Y.-model.likelihood.γ),model.Σ))
    return nothing
end

""" Return the gradient of the expectation for latent GP `index` """
function expec_μ(model::VGP{<:AugmentedMultinomialLogisticLikelihood},index::Int)
    0.5.*(model.likelihood.Y[index]-model.likelihood.γ[index])
end

function ∇μ(model::VGP{<:AugmentedMultinomialLogisticLikelihood})
    0.5.*(model.likelihood.Y.-model.likelihood.γ)
end

""" Return the gradient of the expectation for latent GP `index` """
function expec_μ(model::SVGP{<:AugmentedMultinomialLogisticLikelihood},index::Integer)
    0.5.*(model.likelihood.Y[index][model.inference.MBIndices]-model.likelihood.γ[index])
end

function ∇μ(model::SVGP{<:AugmentedMultinomialLogisticLikelihood})
    0.5.*(getindex.(model.likelihood.Y,[model.inference.MBIndices]).-model.likelihood.γ)
end

function expec_Σ(model::AbstractGP{<:AugmentedMultinomialLogisticLikelihood},index::Integer)
    0.5.*model.likelihood.θ[index]
end

function ∇Σ(model::AbstractGP{AugmentedMultinomialLogisticLikelihood{T}}) where {T<:Real}
    Diagonal{T}.(0.5.*model.likelihood.θ)
end

function ELBO(model::AbstractGP{<:AugmentedMultinomialLogisticLikelihood})
    return expecLogLikelihood(model) - GaussianKL(model) - GammaImproperKL(model) - PoissonKL(model) - PolyaGammaKL(model)
end

function expecLogLikelihood(model::VGP{<:AugmentedMultinomialLogisticLikelihood})
    model.likelihood.c .= broadcast((Σ,μ)->sqrt.(Σ.+abs2.(μ)),diag.(model.Σ),model.μ)
    tot = -model.nSample*logtwo
    tot += -sum(sum.(model.likelihood.γ))*logtwo
    tot +=  0.5*sum(broadcast((y,μ,γ,θ,c)->sum(μ.*(y-γ)-θ.*abs2.(c)),
                    model.likelihood.Y,model.μ,model.likelihood.γ,model.likelihood.θ,model.likelihood.c))
    return tot
end

function expecLogLikelihood(model::SVGP{<:AugmentedMultinomialLogisticLikelihood})
    model.likelihood.c .= broadcast((μ::AbstractVector,Σ::AbstractMatrix,κ::AbstractMatrix,K̃::AbstractVector)->sqrt.(K̃+opt_diag(κ*Σ,κ)+abs2.(κ*μ)),
                                    model.μ,model.Σ,model.κ,model.K̃)
    tot = -model.inference.nSamplesUsed*logtwo
    tot += -sum(sum.(model.likelihood.γ))*logtwo
    tot += 0.5*sum(broadcast((y,κμ,γ,θ,c)->sum((κμ).*(y[model.inference.MBIndices]-γ)-θ.*abs2.(c)),
                    model.likelihood.Y,model.κ.*model.μ,model.likelihood.γ,model.likelihood.θ,model.likelihood.c))
    return model.inference.ρ*tot
end

function grad_samples(model::AbstractGP{<:MultinomialLogisticLikelihood,<:NumericalVI,T},samples::AbstractMatrix{T},index::Integer) where {T<:Real}
    class = model.likelihood.y_class[index]::Int64
    grad_μ = zeros(T,model.nLatent)
    grad_Σ = zeros(T,model.nLatent)
    g_μ = similar(grad_μ)
    nSamples = size(samples,1)
    @views @inbounds for i in 1:nSamples
        σ = logistic.(samples[i,:])
        samples[i,:]  .= MultinomialLogistic(samples[i,:])
        s = samples[i,class]
        g_μ .= grad_MultinomialLogistic(samples[i,:],σ,class)/s
        grad_μ += g_μ
        grad_Σ += diaghessian_MultinomialLogistic(samples[i,:],σ,class)/s - abs2.(g_μ)
    end
    for k in 1:model.nLatent
        model.inference.∇μE[k][index] = grad_μ[k]/nSamples
        model.inference.∇ΣE[k][index] = 0.5.*grad_Σ[k]/nSamples
    end
end

function log_like_samples(model::AbstractGP{<:MultinomialLogisticLikelihood,<:Inference,T},samples::AbstractMatrix,index::Integer) where {T<:Real}
    class = model.likelihood.y_class[index]
    nSamples = size(samples,1)
    loglike = zero(T)
    for i in 1:nSamples
        σ = logistic.(samples[i,:])
        loglike += log(σ[class])-log(sum(σ))
    end
    return loglike/nSamples
end

function remove_augmentation(l::Type{AugmentedMultinomialLogisticLikelihood{T}}) where T
    return MultinomialLogisticLikelihood{T}(l.Y,l.class_mapping,l.ind_mapping)
end

function grad_MultinomialLogistic(s::AbstractVector{<:Real},σ::AbstractVector{<:Real},i::Integer)
    s[i]*(δ.(i,eachindex(σ)).-s).*(1.0.-σ)
end

function diaghessian_MultinomialLogistic(s::AbstractVector{<:Real},σ::AbstractVector{<:Real},i::Integer)
    s[i]*(1.0.-σ).*(
    abs2.(δ.(i,eachindex(σ))-s).*(1.0.-σ)
    -s.*(1.0.-s).*(1.0.-σ)
    -σ.*(δ.(i,eachindex(σ))-s))
end

function hessian_MultinomialLogistic(s::AbstractVector{T},σ::AbstractVector{T},i::Integer) where {T<:Real}
    m = length(s)
    hessian = zeros(T,m,m)
    @inbounds for j in 1:m
        for k in 1:m
            hessian[j,k] = (1-σ[j])*s[i]*(
            (δ(i,k)-s[k])*(1.0-σ[k])*(δ(i,j)-s[j])
            -s[j]*(δ(j,k)-s[k])*(1.0-σ[k])
            -δ(k,j)*σ[j]*(δ(i,j)-s[j]))
        end
    end
    return hessian
end


text:
# For testing Composite
struct Foo
    x
    y::Float64
end

# For testing Composite: it is an invarient of the type that x2 = 2x
# so simple addition can not be defined
struct StructWithInvariant
    x
    x2

    StructWithInvariant(x) = new(x, 2x)
end

@testset "Composite" begin
    @testset "convert" begin
        @test convert(NamedTuple, Composite{Foo}(x=2.5)) == (; x=2.5)
        @test convert(Tuple, Composite{Tuple{Float64,}}(2.0)) == (2.0,)
    end

    @testset "indexing, iterating, and properties" begin
        @test propertynames(Composite{Foo}(x=2.5)) == (:x,)
        @test Composite{Foo}(x=2.5).x == 2.5

        @test propertynames(Composite{Tuple{Float64,}}(2.0)) == (1,)
        @test getproperty(Composite{Tuple{Float64,}}(2.0), 1) == 2.0

        @test length(Composite{Foo}(x=2.5)) == 1
        @test length(Composite{Tuple{Float64,}}(2.0)) == 1

        @test eltype(Composite{Foo}(x=2.5)) == Float64
        @test eltype(Composite{Tuple{Float64,}}(2.0)) == Float64

        # Testing iterate via collect
        @test collect(Composite{Foo}(x=2.5)) == [2.5]
        @test collect(Composite{Tuple{Float64,}}(2.0)) == [2.0]
    end

    @testset "conj" begin
        @test conj(Composite{Foo}(x=2.0+3.0im)) == Composite{Foo}(x=2.0-3.0im)
        @test ==(
            conj(Composite{Tuple{Float64,}}(2.0+3.0im)),
            Composite{Tuple{Float64,}}(2.0-3.0im)
        )
    end

    @testset "extern" begin
        @test extern(Composite{Foo}(x=2.0)) == (;x=2.0)
        @test extern(Composite{Tuple{Float64,}}(2.0)) == (2.0,)

        # with differentials on the inside
        @test extern(Composite{Foo}(x=@thunk(0+2.0))) == (;x=2.0)
        @test extern(Composite{Tuple{Float64,}}(@thunk(0+2.0))) == (2.0,)
    end

    @testset "+ with other composites" begin
        @testset "Structs" begin
            @test (
                Composite{Foo}(x=1.5) +
                Composite{Foo}(x=2.5)
            ) == Composite{Foo}(x=4.0)
            @test (
                Composite{Foo}(y=1.5) +
                Composite{Foo}(x=2.5)
            ) == Composite{Foo}(y=1.5, x=2.5)
            @test (
                Composite{Foo}(y=1.5, x=1.5) +
                Composite{Foo}(x=2.5)
                ) == Composite{Foo}(y=1.5, x=4.0)
        end

        @testset "Tuples" begin
            @test (
                Composite{Tuple{Float64, Float64}}(1.0, 2.0) +
                Composite{Tuple{Float64, Float64}}(1.0, 1.0)
            ) == Composite{Tuple{Float64, Float64}}(2.0, 3.0)
        end

        @testset "NamedTuples" begin
            nt1 = (;a=1.5, b=0.0)
            nt2 = (;a=0.0, b=2.5)
            nt_sum = (a=1.5, b=2.5)
            @test (
                Composite{typeof(nt1)}(; nt1...) +
                Composite{typeof(nt2)}(; nt2...)
            ) == Composite{typeof(nt_sum)}(; nt_sum...)
        end
    end

    @testset "+ with Primals" begin
        @testset "Structs" begin
            @test Foo(3.5, 1.5) + Composite{Foo}(x=2.5) == Foo(6.0, 1.5)
            @test Composite{Foo}(x=2.5) + Foo(3.5, 1.5) == Foo(6.0, 1.5)
        end

        @testset "Tuples" begin
            @test ((1.0, 2.0) + Composite{Tuple{Float64, Float64}}(1.0, 1.0)) == (2.0, 3.0)
            @test (Composite{Tuple{Float64, Float64}}(1.0, 1.0)) + (1.0, 2.0) == (2.0, 3.0)
        end

        @testset "NamedTuple" begin
            ntx = (; a=1.5)
            @test Composite{typeof(ntx)}(; ntx...) + ntx == (; a=3.0)

            nty = (; a=1.5, b=0.5)
            @test Composite{typeof(nty)}(; nty...) + nty == (; a=3.0, b=1.0)
        end

    end

    @testset "+ with Primals, with inner constructor" begin
        value = StructWithInvariant(10.0)
        diff = Composite{StructWithInvariant}(x=2.0, x2=6.0)
        @test_throws ChainRulesCore.PrimalAdditionFailedException (value + diff)
        @test_throws ChainRulesCore.PrimalAdditionFailedException (diff + value)

        # Now we define constuction for ChainRulesCore.jl's purposes:
        # It is going to determine the root quanity of the invarient
        function ChainRulesCore.construct(::Type{StructWithInvariant}, nt::NamedTuple)
            x = (nt.x + nt.x2/2)/2
            return StructWithInvariant(x)
        end
        @test value + diff == StructWithInvariant(12.5)
        @test diff + value == StructWithInvariant(12.5)
    end

    @testset "Scaling" begin
        @test (
            2 *  Composite{Foo}(y=1.5, x=2.5)
            == Composite{Foo}(y=3.0, x=5.0)
            == Composite{Foo}(y=1.5, x=2.5) * 2
        )
        @test (
            2 * Composite{Tuple{Float64, Float64}}(2.0, 4.0)
            == Composite{Tuple{Float64, Float64}}(4.0, 8.0)
            == Composite{Tuple{Float64, Float64}}(2.0, 4.0) * 2
        )
    end

    @testset "show" begin
        @test repr(Composite{Foo}(x=1,)) == "Composite{Foo}(x = 1,)"
        @test repr(Composite{Tuple{Int,Int}}(1, 2)) == "Composite{Tuple{Int64,Int64}}(1, 2)"
    end

    @testset "internals" begin
        @testset "Can't do backing on primative type" begin
            @test_throws Exception ChainRulesCore.backing(1.4)
        end

        @testset "Internals don't allocate a ton" begin
            bk = (; x=1.0, y=2.0)
            @test (@allocated(ChainRulesCore.construct(Foo, bk))) <= 32
            @test (@allocated ChainRulesCore.elementwise_add(bk, bk)) <= 48
        end
    end
end

ext: jl


text:
using Test
using DebugAdapter

@testset "DebugAdapter" begin
end


text:
include("history.jl")
using Test

b1 = UInt128(90000000000000000000000000000000000005)
b0 = UInt128(200000000000000000000000000000000000007)
b2 = UInt128(90000000000000000000000000000000000000)
x0 = UInt256(b0, b1)
x1 = UInt256(b0, b2)

@testset "UInt256" begin
    @test UInt256(7) == UInt256(0, 7)
    @test UInt256(10^10) < UInt256(1, 0)
    @test UInt256(7) > UInt256(0)
    @test (x0 >> 3) << 3 == x1
    @test (x0 >>> 3) << 3 == x1
    @test x0 | x1 == x0
    @test x0 | UInt256(0) == x0
    @test x0 ⊻ x0 == UInt256(0)
    @test x0 & x1 == x1
    @test peek(x0, UInt8(255))
    @test !peek(x0, UInt8(254))
    @test peek(x0, UInt8(0))
    @test !peek(x0, UInt8(1))
    @test togglebit(togglebit(x0, UInt8(8)), UInt8(8)) == x0
    @test zero(x0) == UInt256(0)
    @test zero(UInt256) == UInt256(0)
    @test one(x0) == UInt256(1)
    @test one(UInt256) == UInt256(1)
end;

@testset "Bidding" begin
    @test bid2index(2,1) == 0x16
    @test index2bid(0x16) == (2,1)
    @test index2bid(0x17) == (0, 0)
    @test parsebid(42) == (4,2)
    @test addbid(x0, 3, 7) == addbid(x0, 37)
    @test addbid(UInt256(0), 1, 0) == UInt256(0, 1)
    x = addbid(addbid(zero(UInt256), 20), 21)
    println(bidhistory(x))

end;

max_stars_repo_name: jeffreyrosenbluth/LP

max_issues_repo_name: jeffreyrosenbluth/LP

max_forks_repo_name: jeffreyrosenbluth/LP


text:
module CancellationTokens

export CancellationTokenSource, get_token, is_cancellation_requested, cancel, OperationCanceledException

include("event.jl")

@enum CancellationTokenSourceStates NotCanceledState=1 NotifyingState=2 NotifyingCompleteState=3

mutable struct CancellationTokenSource
    _state::CancellationTokenSourceStates
    _timer::Union{Nothing,Timer}
    _kernel_event::Union{Nothing,Event} # TODO Event is Julia > 1.1, make it work on 1.0

    function CancellationTokenSource()
        return new(NotCanceledState, nothing, nothing)
    end
end

function CancellationTokenSource(timespan_in_seconds::Real)
    x = CancellationTokenSource()

    x._timer = Timer(timespan_in_seconds) do _
        _internal_notify(x)
    end

    return x
end

function _internal_notify(x::CancellationTokenSource)
    if x._state==NotCanceledState
        x._state = NotifyingState

        if x._timer!==nothing
            close(x._timer)
            x._timer = nothing
        end

        if x._kernel_event!==nothing
            notify(x._kernel_event)
            x._kernel_event = nothing
        end

        x._state = NotifyingCompleteState
    end
end

function cancel(x::CancellationTokenSource)
    _internal_notify(x)

    return
end

is_cancellation_requested(x::CancellationTokenSource) = x._state > NotCanceledState

function _waithandle(x::CancellationTokenSource)
    if x._kernel_event===nothing
        x._kernel_event = Event()
    end

    return x._kernel_event
end

# CancellationToken

struct CancellationToken
    _source::CancellationTokenSource
end

get_token(x::CancellationTokenSource) = CancellationToken(x)

is_cancellation_requested(x::CancellationToken) = is_cancellation_requested(x._source)

_waithandle(x::CancellationToken) = _waithandle(x._source)

function Base.wait(x::CancellationToken)
    if is_cancellation_requested(x)
        return
    else
        wait(_waithandle(x))
    end
end

# OperationCanceledException

struct OperationCanceledException <: Exception
    _token::CancellationToken
end

get_token(x::OperationCanceledException) = x._token

function CancellationTokenSource(tokens::CancellationToken...)
    x = CancellationTokenSource()

    for t in tokens
        @async begin
            wait(t)
            _internal_notify(x)
        end
    end

    return x
end

include("augment_base.jl")

end # module


text:
# Provides a representative plot demonstrating how target weight varies with
# distance

using SubmodularMaximization
using PyPlot

close("all")

grid_size = 41
horizon = SubmodularMaximization.default_horizon

range_limits = [5.0, 10.0, 15.0, 20.0, Inf]

for range_limit in range_limits
  sensor = RangingSensor(range_limit=range_limit)

  # place the target in the center of the grid
  target_state = (21,1)

  grid = Grid(grid_size, grid_size)

  # Create some uncertainty. Two steps is probably a reasonable approximation of
  # the kind of target uncertainty that we will encounter
  histogram_filter = Filter(grid, target_state)
  process_update!(histogram_filter, transition_matrix(grid))

  range = 1:grid_size

  @time weights = thread_map(range) do x
    state = (x, 1)

    problem = MultiRobotTargetTrackingProblem([state, state],
                                              [histogram_filter],
                                              grid=grid,
                                              sensor=sensor,
                                              horizon=horizon)

    weights = compute_weight_matrix(problem)

    sum(weights) / 2
  end

  label = string("Range. lim.: ", range_limit)
  plot(range, weights, label=label)

  ratio = maximum(weights)/minimum(weights)

  println("Range limit: ", range_limit, " Weights falloff ratio: ", ratio)
end

plot_title = string("Weights Falloff")

title(plot_title)
xlabel("X position")
ylabel("Weight (bits)")
legend()

save_fig("fig", plot_title)

nothing

max_stars_repo_name: mcorah/MultiAgentSensing

max_issues_repo_name: mcorah/MultiAgentSensing

max_forks_repo_name: mcorah/MultiAgentSensing


text:
module Day20

export day20

using CombinedParsers
using CombinedParsers.Regexp

function parseline(line)
    p = join(Repeat(CharIn(".#")), "")
    line |> p
end

function parseinput(filename)
    temptiles = Dict{Int,Vector{Char}}()
    id = 0
    for line in readlines(filename)
        line = strip(line)
        if line == ""
            continue
        elseif line[1:4] == "Tile"
            id = parse(Int, line[6:end-1])
        else
            temptiles[id] = vcat(get(temptiles,id,Vector{Char}()), parseline(line))
        end
    end

    tiles = Dict{Int,Vector{String}}()

    for (id,tile) in temptiles
        tmp = reshape(tile,10,10)
        tiles[id] = map(join, [tmp[1,:], tmp[:,end], tmp[end,:], tmp[:,1]])
    end
    tiles
end

function match_against(side, tile)
    revside = reverse(side)
    for side1 in tile
        if side1 == side || side1 == revside
            return true
        end
    end
    false
end

function day20(::Val{:a}, filename)
    tiles = parseinput(filename)
    corners = Vector{Int}()
    for (id1,tile1) in tiles
        matched_sides = 0
        for side1 in tile1
            for (id2,tile2) in tiles
                if id1 != id2 && match_against(side1,tile2)
                    matched_sides += 1
                    break
                end
            end
        end
        if matched_sides <= 2
            push!(corners, id1)
        end
    end
    prod(corners)
end

end


text:
#=
This file is auto-generated. Do not edit.
=#

"""Data Structure for the procurement products for system simulations."""
mutable struct StaticReserve <: Reserve
    name::String
    contributingdevices::Vector{Device}  # devices from which the product can be procured
    timeframe::Float64  # the relative saturation timeframe
    requirement::Float64  # the required quantity of the product should be scaled by a Forecast
    _forecasts::InfrastructureSystems.Forecasts
    internal::InfrastructureSystemsInternal
end

function StaticReserve(name, contributingdevices, timeframe, requirement, _forecasts=InfrastructureSystems.Forecasts(), )
    StaticReserve(name, contributingdevices, timeframe, requirement, _forecasts, InfrastructureSystemsInternal())
end

function StaticReserve(; name, contributingdevices, timeframe, requirement, _forecasts=InfrastructureSystems.Forecasts(), )
    StaticReserve(name, contributingdevices, timeframe, requirement, _forecasts, )
end

# Constructor for demo purposes; non-functional.

function StaticReserve(::Nothing)
    StaticReserve(;
        name="init",
        contributingdevices=[ThermalStandard(nothing)],
        timeframe=0.0,
        requirement=0.0,
        _forecasts=InfrastructureSystems.Forecasts(),
    )
end

"""Get StaticReserve name."""
get_name(value::StaticReserve) = value.name
"""Get StaticReserve contributingdevices."""
get_contributingdevices(value::StaticReserve) = value.contributingdevices
"""Get StaticReserve timeframe."""
get_timeframe(value::StaticReserve) = value.timeframe
"""Get StaticReserve requirement."""
get_requirement(value::StaticReserve) = value.requirement
"""Get StaticReserve _forecasts."""
get__forecasts(value::StaticReserve) = value._forecasts
"""Get StaticReserve internal."""
get_internal(value::StaticReserve) = value.internal


text:
abstract type Model{T} end

floating_base_dim(::Model) = 3
nominal_configuration(::Model) = zeros(model.nq)


text:
module FighterJets

export FighterJet

"FighterJet is a very fast vehicle with powerful weapons."
mutable struct FighterJet

    "power status: true = on, false = off"
    power::Bool

    "current direction in radians"
    direction::Float64

    "current position coordinate (x,y)"
    position::Tuple{Float64, Float64}

end

# Import generic functions
import Vehicle: power_on!, power_off!, turn!, move!, position

# Implementation of Vehicle interface

function power_on!(fj::FighterJet)
    fj.power = true
    println("Powered on: ", fj)
    nothing
end

function power_off!(fj::FighterJet)
    fj.power = false
    println("Powered off: ", fj)
    nothing
end

function turn!(fj::FighterJet, direction)
    fj.direction = direction
    println("Changed direction to ", direction, ": ", fj)
    nothing
end

function move!(fj::FighterJet, distance)
    x, y = fj.position
    dx = round(distance * cos(fj.direction), digits = 2)
    dy = round(distance * sin(fj.direction), digits = 2)
    fj.position = (x + dx, y + dy)
    println("Moved (", dx, ",", dy, "): ", fj)
    nothing
end

function position(fj::FighterJet)
    fj.position
end

end # module


text:
#=
!NOTES: The NUTS kernel and stepsize adaption are adapted from the excellent package
  https://github.com/tpapp/DynamicHMC.jl

The DynamicHMC.jl package is licensed under the MIT "Expat" License:

> Copyright (c) 2020: Tamas K. Papp.
>
> Permission is hereby granted, free of charge, to any person obtaining a copy
> of this software and associated documentation files (the "Software"), to deal
> in the Software without restriction, including without limitation the rights
> to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
> copies of the Software, and to permit persons to whom the Software is
> furnished to do so, subject to the following conditions:
>
> The above copyright notice and this permission notice shall be included in all
> copies or substantial portions of the Software.
>
> THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
> IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
> FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
> AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
> LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
> OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
> SOFTWARE.
>
=#

############################################################################################
# Constants
"Maximum number of iterations"
const MAX_DIRECTIONS_DEPTH = 32

"Default maximum depth for trees."
const DEFAULT_MAX_TREE_DEPTH = 10

############################################################################################
#Include
include("treebuilding.jl")
include("kernel.jl")
include("config.jl")
include("diagnostics.jl")
include("postprocessing.jl")

############################################################################################
#export


text:
using Random, SnpArrays, DataFrames, GLM
using Suppressor
using LinearAlgebra
using BenchmarkTools
beta = [1.0, 5.0]

n = 10
p = 2
d = 2
m = 2

function generateSPDmatrix(n)
	A = rand(n)
	m = 0.5 * (A * A')
	PDmat = m + (n * Diagonal(ones(n)))
end

function generateRandomVCM(n::Int64, p::Int64, d::Int64, m::Int64)
	# n-by-p design matrix
	X = randn(n, p)

	# p-by-d mean component regression coefficient for each trait
	B = hcat(ones(p, 1), rand(p))

	V = ntuple(x -> zeros(n, n), m)
	for i = 1:m-1
	  copy!(V[i], generateSPDmatrix(n))
	end
	copy!(V[end], Diagonal(ones(n))) # last covarianec matrix is identity

	# a tuple of m d-by-d variance component parameters
	Σ = ntuple(x -> zeros(d, d), m)
	for i in 1:m
	  copy!(Σ[i], generateSPDmatrix(d))
	end
	return(X, B, Σ, V)
end

X, B, Σ, V = generateRandomVCM(n, p, d, m)
dist = Normal()
link = IdentityLink()
glmtraitobject = GLMTrait(X, beta, dist, link)
y_normal  = simulate(glmtraitobject)
@test size(y_normal) == size(X*beta)

@test_throws ErrorException("function not supported by $(typeof(glmtraitobject))") nvc(glmtraitobject)

@test_throws ErrorException("function not supported by $(typeof(glmtraitobject))") ntraits(glmtraitobject)

@test_throws ErrorException("function not supported by $(typeof(glmtraitobject))") noutcomecategories(glmtraitobject)

number_independent_simulations  = 5
@test length(simulate(glmtraitobject, number_independent_simulations)) == number_independent_simulations

@test glmtraitobject isa TraitSimulation.AbstractTraitModel

# make sure that ordered multinomial models works
link = LogitLink()
θ = [1.0, 1.2, 1.4]
Ordinal_Model_Test = OrderedMultinomialTrait(X, beta, θ, link)

@test Ordinal_Model_Test isa TraitSimulation.AbstractTraitModel
@test eltype(simulate(Ordinal_Model_Test)) == Int64
@test length(simulate(Ordinal_Model_Test, number_independent_simulations))  == number_independent_simulations

@test noutcomecategories(Ordinal_Model_Test) == 4

# make sure GLM simulation works
dist = Poisson()
link = IdentityLink()
beta = [21.0, 5.0]
glmtraitobject2 = GLMTrait(abs.(X), beta, dist, link)

# check if clamper worked
@test sum(glmtraitobject2.η .> 20.0) == 0

@test  eltype(simulate(glmtraitobject2)) == Int64

dist = Bernoulli()
link = LogitLink()
glmtraitobject3 = GLMTrait(X, beta, dist, link)

@test eltype(simulate(glmtraitobject3)) == Bool

@test length(simulate(glmtraitobject3, number_independent_simulations))  == number_independent_simulations

dist = NegativeBinomial()
link = LogLink()
glmtraitobject4 = GLMTrait(X, beta, dist, link)

@test eltype(simulate(glmtraitobject4)) == Int64

@test length(simulate(glmtraitobject4, number_independent_simulations)) == number_independent_simulations

dist = Gamma()
link = LogLink()
glmtraitobject5 = GLMTrait(X, beta, dist, link)

@test eltype(simulate(glmtraitobject5)) == Float64

@test length(simulate(glmtraitobject5, number_independent_simulations)) == number_independent_simulations

varcomp = @vc Σ[1] ⊗ V[1] + Σ[2] ⊗ V[2]
# make sure GLMM works
glmtraitobject6 = GLMMTrait(X, B, varcomp, dist, link)
@test glmtraitobject6 isa TraitSimulation.AbstractTraitModel

@test nsamplesize(glmtraitobject6) == n
@test neffects(glmtraitobject6) == p
@test nvc(glmtraitobject6) == m
@test ntraits(glmtraitobject6) == d

y_glmm = simulate(glmtraitobject6)
@test size(y_glmm) == (size(X,  1),  size(B, 2))

@test length(simulate(glmtraitobject6, number_independent_simulations)) == number_independent_simulations

# simulate the SnpArray
G = snparray_simulation([0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5], n)
γ = rand(7)

glmtraitobject7 =  GLMTrait(X, beta, G, γ, dist, link)

@test glmtraitobject7.dist  == Gamma
@test eltype(simulate(glmtraitobject7)) == eltype(dist)

@test length(simulate(glmtraitobject7, number_independent_simulations)) == number_independent_simulations

@test GLMTrait(X, beta, G, γ, dist, link) != empty

@test nsamplesize(glmtraitobject7) == n
@test neffects(glmtraitobject7) == p

model_info_ordinal = @capture_out show(Ordinal_Model_Test)
model_info_glmm = @capture_out show(glmtraitobject6)
model_info_glm = @capture_out show(glmtraitobject7)

@test model_info_ordinal != nothing
@test model_info_glmm != nothing
@test model_info_glm != nothing


text:
using GaussianSimulation
using GeoStatsBase
using Variography
using Plots, VisualRegressionTests
using Test, Pkg, Random

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

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

# list of tests
testfiles = [
  "lu.jl",
  "fft.jl",
  "seq.jl"
]

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


text:
export ProbabilityDistribution, ConditionalDistribution, JointProbabilityDistribution
export is_probability_distribution, is_conditional_distribution
export Probabilities, Conditionals, JointProbabilities
export outcome_probabilities

"""
    abstract type ProbabilityDistribution <: AbstractVector{Real} end

The abstract type representing a discrete probability distribution.
"""
abstract type ProbabilityDistribution{T} <: AbstractVector{T} end
Base.size(probs::ProbabilityDistribution) = size(probs.distribution)
Base.getindex(probs::ProbabilityDistribution, id::Int) = getindex(probs.distribution, id...)
Base.setindex!(probs::ProbabilityDistribution, v, id::Int) = (probs.distribution[id...] = v)

"""
    abstract type ConditionalDistribution <: AbstractMatrix{Real} end

The abstract type representing a discrete conditional probability distribution.
"""
abstract type ConditionalDistribution{T} <: AbstractMatrix{T} end
Base.size(probs::ConditionalDistribution) = size(probs.distribution)
Base.getindex(probs::ConditionalDistribution, id::Vararg{Int,2}) = getindex(probs.distribution, id...)
Base.setindex!(probs::ConditionalDistribution, v, id::Vararg{Int,2}) = (probs.distribution[id...] = v)

"""
    abstract type JointProabilityDistribution <: AbstractMatrix{Real} end

The abstract type representing a discrete joint probability distribution.
"""
abstract type JointProbabilityDistribution{T} <: AbstractMatrix{T} end
Base.size(probs::JointProbabilityDistribution) = size(probs.distribution)
Base.getindex(probs::JointProbabilityDistribution, id::Vararg{Int,2}) = getindex(probs.distribution, id...)
Base.setindex!(probs::JointProbabilityDistribution, v, id::Vararg{Int,2}) = (probs.distribution[id...] = v)

"""
    is_probability_distribution(
        probabilities :: Vector{<:Real};
        atol=ATOL :: Float64
    ) :: Bool

Returns `true` if the provided vector is a valid probability distribution:
* `sum(probabilities) ≈ 1`
* `probabilities[i] ≥ 0 ∀ i`
"""
function is_probability_distribution(probabilities :: AbstractVector{<:Real}; atol=ATOL :: Float64) :: Bool
    _is_prob(prob :: Real) :: Bool = begin
        (prob ≥ 0) || isapprox(prob, 0, atol=atol)
    end

    if findfirst(!_is_prob, probabilities) !== nothing
        return false
    end

    if !isapprox(sum(probabilities), 1, atol=atol)
        return false
    end

    return true
end
is_probability_distribution(::ProbabilityDistribution) = true
is_probability_distribution(::JointProbabilityDistribution) = true

"""
    is_conditional_distribution( probabilities :: AbstractMatrix{<:Real}; atol=ATOL :: Flaot64 ) :: Bool

Returns `true` if the provided matrix is column stochastic. That is, each column
is a valid probability distribution.
"""
function is_conditional_distribution(conditionals :: Matrix{<:Real}; atol=ATOL :: Float64) :: Bool
    all(col -> is_probability_distribution(col, atol=atol), eachcol(conditionals))
end
is_conditional_distribution(::ConditionalDistribution) = true

"""
    Probabilities( distribution :: Vector{<:Real}; atol=ATOL :: Float64 ) <: ProbabilityDistribution

A struct representing a discrete probability distribution. All elements
in the marginal distribution must be positive and their sum must be one.
"""
struct Probabilities{T} <: ProbabilityDistribution{T}
    distribution :: Vector{T}
    atol :: Float64
    Probabilities(
        distribution :: Vector{<:Real};
        atol=ATOL :: Float64
    ) = is_probability_distribution(distribution, atol=atol) ? new{eltype(distribution)}(distribution,atol) : throw(
        DomainError(distribution, "is not a valid probability distribution.")
    )
end

"""
    Conditionals( distribution :: Matrix{<:Real}; atol=ATOL :: Float64 ) <: ConditionalDistribution

A struct representing a conditional probability distribution function. `Conditionals`
are organized in a matrix with rows correpsonding to outputs and columns corresponding
to inputs. For example if there are ``M`` inputs and ``N`` outputs, the corresponding
conditionals matrix takes the form:

```math
\\begin{bmatrix}
p(1|1) & \\dots & p(1|M) \\\\
\\vdots & \\ddots & \\vdots \\\\
p(N|1) & \\dots & p(N|M) \\\\
\\end{bmatrix}
```
"""
struct Conditionals{T} <: ConditionalDistribution{T}
    distribution :: Matrix{T}
    atol :: Float64
    Conditionals(
        distribution :: Matrix{<:Real};
        atol=ATOL :: Float64
    ) = is_conditional_distribution(
            distribution, atol=atol
        ) ? new{eltype(distribution)}(distribution, atol) : throw(
            DomainError(distribution, "is not a valid conditional probability distribution.")
        )
end

"""
    JointProbabilities(
        distribution :: Matrix{<:Real};
        atol=ATOL :: Float64
    ) <: JointProbabilityDistribution

A struct representing a discrete probability distribution. All elements
in the marginal distribution must be positive and their sum must be one.
For convenience, the `JointProbabilities` can also be constructed by passing in
a `ConditionalDistribution` and a `ProbabilityDistribution`,

    JointProbabilities(
        priors :: AbstractVector{<:Real},
        conditionals :: AbstractMatrix{<:Real};
        atol=ATOL :: Float64
    )

Or, two `ProbabilityDistribution`s can be provided

    JointProbabilities(
        priors1 :: AbstractVector{<:Real},
        priors2 :: AbstractVector{<:Real};
        atol=ATOL :: Float64
    )

A `DomainError` is thrown if the joint probability distribution is invalid.
"""
struct JointProbabilities{T} <: JointProbabilityDistribution{T}
    distribution :: Matrix{T}
    atol :: Float64
    JointProbabilities(
        distribution :: Matrix{<:Real};
        atol=ATOL :: Float64
    ) = is_probability_distribution(
            distribution[:], atol=atol
        ) ? new{eltype(distribution)}(distribution, atol) : throw(
            DomainError(distribution, "is not a valid joint probability distribution.")
        )
    JointProbabilities(
        priors :: AbstractVector{<:Real},
        conditionals :: AbstractMatrix{<:Real};
        atol=ATOL :: Float64
    ) = JointProbabilities(priors' .* conditionals, atol=atol)
    JointProbabilities(
        conditionals :: AbstractMatrix{<:Real},
        priors :: AbstractVector{<:Real};
        atol=ATOL :: Float64
    ) = JointProbabilities(priors' .* conditionals, atol=atol)
    JointProbabilities(
        priors1 :: AbstractVector{<:Real},
        priors2 :: AbstractVector{<:Real};
        atol=ATOL :: Float64
    ) = JointProbabilities(priors2 * priors1', atol=atol)
end

"""
    *(
        C :: ConditionalDistribution,
        P :: ProbabilityDistribution;
        atol=ATOL :: Float64
    ) :: ProbabilityDistribution

Multiplication of a conditional with a prior distribution yields the outcome probabilites.
"""
*(
    conditional :: ConditionalDistribution,
    prior :: ProbabilityDistribution;
    atol=ATOL :: Float64
) :: ProbabilityDistribution = Probabilities(conditional.distribution * prior.distribution, atol=atol)

"""
    outcome_probabilities(
        conditionals::ConditionalDistribution,
        priors :: ProbabilityDistribution;
        atol=ATOL :: Float64
    ) :: ProbabilityDistribution

Returns the probability of each outcome given priors and conditional probabilities.
For convenience, this method can be called with an abstract vector/matrix

    outcome_probabilities(
        conditionals :: AbstractMatrix{<:Real},
        priors :: AbstractVector{<:Real};
        atol=ATOL :: Float64
    )

Furthermore, the ordering of the arguments can be reversed.
A `DomainError` is thrown if the constructed vector is not a valid probability
distribution.
"""
function outcome_probabilities(
    conditionals :: ConditionalDistribution,
    priors :: ProbabilityDistribution;
    atol=ATOL :: Float64
) :: ProbabilityDistribution
    *(conditionals,priors,atol=atol)
end
function outcome_probabilities(
    priors :: ProbabilityDistribution,
    conditionals :: ConditionalDistribution;
    atol=ATOL :: Float64
) :: ProbabilityDistribution
    *(conditionals,priors, atol=atol)
end
function outcome_probabilities(
    conditionals :: AbstractMatrix{<:Real},
    priors :: AbstractVector{<:Real};
    atol=ATOL :: Float64
) :: ProbabilityDistribution
    Probabilities(conditionals*priors, atol=atol)
end
function outcome_probabilities(
    priors :: AbstractVector{<:Real},
    conditionals :: AbstractMatrix{<:Real};
    atol=ATOL :: Float64
) :: ProbabilityDistribution
    Probabilities(conditionals*priors, atol=atol)
end


text:
"""
    LogNormal(μ,σ)

The *log normal distribution* is the distribution of the exponential of a [`Normal`](@ref) variate: if ``X \\sim \\operatorname{Normal}(\\mu, \\sigma)`` then
``\\exp(X) \\sim \\operatorname{LogNormal}(\\mu,\\sigma)``. The probability density function is

```math
f(x; \\mu, \\sigma) = \\frac{1}{x \\sqrt{2 \\pi \\sigma^2}}
\\exp \\left( - \\frac{(\\log(x) - \\mu)^2}{2 \\sigma^2} \\right),
\\quad x > 0
```

```julia
LogNormal()          # Log-normal distribution with zero log-mean and unit scale
LogNormal(mu)        # Log-normal distribution with log-mean mu and unit scale
LogNormal(mu, sig)   # Log-normal distribution with log-mean mu and scale sig

params(d)            # Get the parameters, i.e. (mu, sig)
meanlogx(d)          # Get the mean of log(X), i.e. mu
varlogx(d)           # Get the variance of log(X), i.e. sig^2
stdlogx(d)           # Get the standard deviation of log(X), i.e. sig
```

External links

* [Log normal distribution on Wikipedia](http://en.wikipedia.org/wiki/Log-normal_distribution)

"""
struct LogNormal{T<:Real} <: ContinuousUnivariateDistribution
    μ::T
    σ::T
    LogNormal{T}(μ::T, σ::T) where {T} = new{T}(μ, σ)
end

function LogNormal(μ::T, σ::T; check_args=true) where {T <: Real}
    check_args && @check_args(LogNormal, σ ≥ zero(σ))
    return LogNormal{T}(μ, σ)
end

LogNormal(μ::Real, σ::Real) = LogNormal(promote(μ, σ)...)
LogNormal(μ::Integer, σ::Integer) = LogNormal(float(μ), float(σ))
LogNormal(μ::T) where {T <: Real} = LogNormal(μ, one(T))
LogNormal() = LogNormal(0.0, 1.0, check_args=false)

@distr_support LogNormal 0.0 Inf

#### Conversions
convert(::Type{LogNormal{T}}, μ::S, σ::S) where {T <: Real, S <: Real} = LogNormal(T(μ), T(σ))
convert(::Type{LogNormal{T}}, d::LogNormal{S}) where {T <: Real, S <: Real} = LogNormal(T(d.μ), T(d.σ), check_args=false)

#### Parameters

params(d::LogNormal) = (d.μ, d.σ)
partype(::LogNormal{T}) where {T} = T

#### Statistics

meanlogx(d::LogNormal) = d.μ
varlogx(d::LogNormal) = abs2(d.σ)
stdlogx(d::LogNormal) = d.σ

mean(d::LogNormal) = ((μ, σ) = params(d); exp(μ + σ^2/2))
median(d::LogNormal) = exp(d.μ)
mode(d::LogNormal) = ((μ, σ) = params(d); exp(μ - σ^2))

function var(d::LogNormal)
    (μ, σ) = params(d)
    σ2 = σ^2
    (exp(σ2) - 1) * exp(2μ + σ2)
end

function skewness(d::LogNormal)
    σ2 = varlogx(d)
    e = exp(σ2)
    (e + 2) * sqrt(e - 1)
end

function kurtosis(d::LogNormal)
    σ2 = varlogx(d)
    e = exp(σ2)
    e2 = e * e
    e3 = e2 * e
    e4 = e3 * e
    e4 + 2*e3 + 3*e2 - 6
end

function entropy(d::LogNormal)
    (μ, σ) = params(d)
    (1 + log(twoπ * σ^2))/2 + μ
end

#### Evalution

function pdf(d::LogNormal, x::Real)
    if x ≤ zero(x)
        logx = log(zero(x))
        x = one(x)
    else
        logx = log(x)
    end
    return pdf(Normal(d.μ, d.σ), logx) / x
end

function logpdf(d::LogNormal, x::Real)
    if x ≤ zero(x)
        logx = log(zero(x))
        b = zero(logx)
    else
        logx = log(x)
        b = logx
    end
    return logpdf(Normal(d.μ, d.σ), logx) - b
end

function cdf(d::LogNormal, x::Real)
    logx = x ≤ zero(x) ? log(zero(x)) : log(x)
    return cdf(Normal(d.μ, d.σ), logx)
end

function ccdf(d::LogNormal, x::Real)
    logx = x ≤ zero(x) ? log(zero(x)) : log(x)
    return ccdf(Normal(d.μ, d.σ), logx)
    end

function logcdf(d::LogNormal, x::Real)
    logx = x ≤ zero(x) ? log(zero(x)) : log(x)
    return logcdf(Normal(d.μ, d.σ), logx)
end

function logccdf(d::LogNormal, x::Real)
    logx = x ≤ zero(x) ? log(zero(x)) : log(x)
    return logccdf(Normal(d.μ, d.σ), logx)
end

quantile(d::LogNormal, q::Real) = exp(quantile(Normal(params(d)...), q))
cquantile(d::LogNormal, q::Real) = exp(cquantile(Normal(params(d)...), q))
invlogcdf(d::LogNormal, lq::Real) = exp(invlogcdf(Normal(params(d)...), lq))
invlogccdf(d::LogNormal, lq::Real) = exp(invlogccdf(Normal(params(d)...), lq))

function gradlogpdf(d::LogNormal, x::Real)
    outofsupport = x ≤ zero(x)
    y = outofsupport ? one(x) : x
    z = (gradlogpdf(Normal(d.μ, d.σ), log(y)) - 1) / y
    return outofsupport ? zero(z) : z
end

# mgf(d::LogNormal)
# cf(d::LogNormal)

#### Sampling

rand(rng::AbstractRNG, d::LogNormal) = exp(randn(rng) * d.σ + d.μ)

## Fitting

function fit_mle(::Type{<:LogNormal}, x::AbstractArray{T}) where T<:Real
    lx = log.(x)
    μ, σ = mean_and_std(lx)
    LogNormal(μ, σ)
end

function StatsBase.fit(::Type{LogNormal}, m::AbstractMoments)
    # https://en.wikipedia.org/wiki/Log-normal_distribution
    n_moments(m) >= 2 || error("Need mean and variance to estimate lognormal")
    γ = 1+var(m)/mean(m)^2
    μ = log(mean(m)/sqrt(γ))
    σ = sqrt(log(γ))
    return LogNormal(μ,σ)
end

function StatsBase.fit(::Type{LogNormal}, lower::QuantilePoint, upper::QuantilePoint)
    #length(qset) == 2 || error("only implemented yet for exactly two quantiles.")
    #qset_log = [QuantilePoint(qp, q = log(qp.q)) for qp in qset]
    lower_log = QuantilePoint(lower, q = log(lower.q))
    upper_log = QuantilePoint(upper, q = log(upper.q))
    DN = fit(Normal, lower_log, upper_log)
    LogNormal(params(DN)...)
end

function fit_mean_quantile(::Type{LogNormal}, mean::Real, qp::QuantilePoint)
    # solution of
    # (1) mean = exp(mu + sigma^2/2)
    # (2) upper = mu + sigmaFac sigma
    # see R packaage lognorm inst/doc/coefLognorm.Rmd for derivation
    sigmaFac = quantile(Normal(),qp.p)
    m = log(mean)
    discr = sigmaFac^2 - 2*(log(qp.q) - m)
    (discr > 0) || error("Cannot fit LogNormal with quantile $(qp) and mean $(mean).")
    sigma = sigmaFac > sqrt(discr) ? (sigmaFac - sqrt(discr)) : sigmaFac + sqrt(discr)
    mu = m - sigma^2/2
    LogNormal(mu, sigma)
end

function fit_mode_quantile(::Type{LogNormal}, mode::Real, qp::QuantilePoint)
    # solution of
    # (1) mle = exp(mu - sigma^2)
    # (2) upper = mu + sigmaFac sigma
    # see R packaage lognorm inst/doc/coefLognorm.Rmd for derivation
    sigmaFac = quantile(Normal(),qp.p)
    m = log(mode)
    discr = sigmaFac^2/4 + (log(qp.q) - m)
    (discr > 0) || error("Cannot fit LogNormal with quantile $(qp) and mode $(mode).")
    root_discr = sqrt(discr)
    mfh = -sigmaFac/2
    sigma = mfh > root_discr ? (mfh - root_discr) : (mfh + root_discr)
    #sigma = mfh + root_discr
    mu = m + sigma^2
    LogNormal(mu, sigma)
end


text:
function ptrace_ref(A::Matrix, out, sysdims::NTuple)
    # Function adapted from a [MATLAB function](http://www.dr-qubit.org/Matlab_code.html) by Toby Cubitt and licensed under GPL2.
    # First, calculate systems, dimensions, etc.
    out = collect(out)
    nsub = length(sysdims) # Number of subsytems
    rdims = reverse(sysdims) # Dimensions of subsystems in reverse order
    keep = setdiff(1:nsub,out) # Indices of subsystems to keep
    dimstrace = prod(sysdims[out])::Int # Total dim of subsystems to trace out
    dimskeep = prod(sysdims[keep])::Int # Total dim of subsystems to keep
    # Reshape density matrix into tensor with one row and one column index
    # for each subsystem, permute traced subsystem indices to the end,
    perm = (nsub + 1) .- [reverse(keep); reverse(keep).-nsub; out; out.-nsub]
    B = permutedims(reshape(A,rdims...,rdims...),perm)
    # Reshape again so that first two indices are row and column
    # multi-indices for kept subsystems and third index is a flattened index
    # for traced subsystems
    C = reshape(B,dimskeep,dimskeep,dimstrace^2)
    # Sum third index over "diagonal" entries
    return dropdims(sum(C[:,:,1:dimstrace+1:dimstrace^2],3),dims=3), sysdims[keep]
end

function ptrace_ref(x::Vector, out::NTuple{no,Int}, sysdims::NTuple{ns,Int}) where {no,ns}
    # Make sure the input arguments make sense
    issorted(out) || throw(ArgumentError("subsystem indices $out must be sorted"))
    no>ns && throw(ArgumentError("more subsystem indices ($no) than number of subsystems ($ns)"))
    for i = 1:no
        out[i]>ns && throw(ArgumentError("subsystem index $(out[i]) larger than number of subsystems"))
    end
    # First, calculate systems, dimensions, etc.
    keep = sorted_setdiff(ntuple(identity,Val(ns)),out) # Subsystems to keep
    rsysdims = revtuple(sysdims) # Dimensions of subsystems in reverse order
    keepdims = gettuple(sysdims,keep) # Dimensions of subsystems to keep
    outdims  = gettuple(sysdims,out) # Dimensions of subsystems to trace out
    # Reshape state vector to "reverse" ket on traced subsystems into a bra
    perm = (ns + 1) .- [collect(revtuple(keep)); collect(out)]
    y = permutedims(reshape(x,rsysdims...),perm)
    y = reshape(y,prod(keepdims),prod(outdims))
    # Take outer product
    return y*y', keepdims
end


text:
module RadioAstroTools

using FITSIO
using WCS
using Unitful, UnitfulAstro
using StatsBase

export AbstractSpectralCube, AbstractMap
export Map, SpectralCube
export create_map, spectralcube
export MapRange, CubeRange
export xy_range, xyz_range
export read_map_from_file, read_cube_from_file
export save
export submap, subcube
export apply_mask, apply_mask!
export copy

# abstract types
abstract type AbstractMap end
abstract type AbstractSpectralCube end

include("helpers.jl")
include("ranges.jl")
include("io.jl")
include("WCShelpers.jl")
include("maps.jl")
include("cubes.jl")
include("index_funcs.jl")
include("mask_funcs.jl")

end


text:
struct GitWorker
    # Git stuff
    sys_root::String
    remote_url::String

    # System data
    dat::Dict{Any, Any}

    GitWorker(sys_root, remote_url) = new(string(sys_root), string(remote_url), Dict{Symbol, Any}())
    GitWorker(;sys_root, remote_url) = GitWorker(sys_root, remote_url)
    GitWorker(remote_url) = GitWorker(_GW_SYSTEM_DFLT_ROOT, remote_url)
end

# overwrite base
import Base.show

function Base.show(io::IO, gw::GitWorker)
    println(io, "GitWorker(;")
    println(io, "   sys_root = \"", gw.sys_root, "\",")
    println(io, "   remote_url = \"", gw.remote_url, "\"")
    println(io, ")")
end

import Base.get!
get!(f::Function, gw::GitWorker, key) = get!(f, gw.dat, key)
get!(gw::GitWorker, key, val) = get!(gw.dat, key, val)

import Base.get
get(f::Function, gw::GitWorker, key) = get(f, gw.dat, key)
get(gw::GitWorker, key, val) = get(gw.dat, key, val)

set!(gw::GitWorker, key, val) = (gw.dat[key] = val)
set!(f::Function, gw::GitWorker, key) = (gw.dat[key] = f())

remote_url(gw::GitWorker) = gw.remote_url
sys_root(gw::GitWorker) = gw.sys_root

max_stars_repo_name: josePereiro/GitWorkers

max_issues_repo_name: josePereiro/GitWorkers

max_forks_repo_name: josePereiro/GitWorkers


text:
#=
Reduce all value to logical, checking if that cell is empty
=#

logical(A::Assoc) = Assoc(A.row,A.col,promote([1.0],A.val)[1],spones(A.A))

########################################################
# D4M: Dynamic Distributed Dimensional Data Model
# Architect: Dr. Jeremy Kepner (kepner@ll.mit.edu)
# Software Engineer: Alexander Chen (alexc89@mit.edu)
########################################################


text:
module TestUtil

using ..LyceumDevTools: flattendir

using BenchmarkTools: BenchmarkTools, memory
using Distributed: Distributed
using LyceumBase
using MacroTools: MacroTools
using Random
using Shapes

using Test
using Test: AbstractTestSet, DefaultTestSet, Result, Pass, Fail, Error, Broken

export @qe
export @test_inferred, @test_noalloc
export @includetests
include("macros.jl")

include("abstractenvironment.jl")
export testenv_correctness, testenv_inferred, testenv_allocations

end # module


text:
module CoupledDipoles

using LinearAlgebra
using Distances

using DifferentialEquations
using Logging: global_logger
using TerminalLoggers: TerminalLogger
global_logger(TerminalLogger())
# using Random
# using Statistics

using Printf
using LsqFit
using Optim: minimizer, optimize, Options
using SpecialFunctions
using Folds
using ThreadPools
using LazyArrays
# using Memoize
using SharedArrays
using HCubature
using MKL

const k₀ = 1
const Γ = 1
const probability_threshold = 0.9999
const R1_threshold = 0.5
const farField_factor = 100

include("structs.jl")
include("atom_cube.jl")
include("atom_sphere.jl")
include("atom_cylinder.jl")

include("formulas.jl")
export ρ_of, b₀_of

include("IO.jl")
export Atom
export Square, Circle
export Cube, Sphere, Cylinder

export get_dimension
export Laser
export PlaneWave2D, PlaneWave3D
export Gaussian2D, Gaussian3D

export Physics, Linear
export LinearOptics, Scalar, Vectorial
export NonLinearOptics, MeanField, BBGKY

include("linear_problems.jl")
include("nonlinear_problems.jl")

include("atoms.jl")
export cube_inputs, sphere_inputs, cylinder_inputs
export get_rₘᵢₙ
export get_pairwise_matrix

include("kernels.jl")
export get_interaction_matrix
export green_scalar!

include("lasers.jl")
export apply_laser_over_atoms, apply_laser_over_oneAtom
export apply_laser_over_sensors, apply_laser_over_oneSensor

include("eigen_analysis.jl")
export get_spectrum
export get_IPRs, get_PRs
export get_localization_length
export get_spatial_profile_single_mode
export classify_modes

include("exponential_fit.jl")
export select_points

include("plot_modes.jl")
export get_spatial_profile_single_mode_and_fit
export get_decay_fit

include("scattering.jl")
export get_intensities_over_sensors
export get_intensity_over_an_angle
export get_intensity_over_an_angle_legacy

include("sensors.jl")
export get_sensors_ring
export get_sensors_sphereSurface

include("dynamics.jl")
export get_steady_state
export time_evolution, default_evolution_initial_condition, get_evolution_function

include("transmission.jl")
export get_transmission, how_far_is_FarField
export _create_sphere_sensor, _create_plane_sensor

end


text:
#=
  (Simple) True skill problem.

  Example inspired by
  Johannes Borgstrom, Andrew D. Gordon, Michael Greenberg, James Margetson, and Jurgen Van Gael:
  "Measure Transformer Semantics for Bayesian Machine Learning"
  https://www.microsoft.com/en-us/research/publication/measure-transformer-semantics-for-bayesian-machine-learning-2011/?from=http%3A%2F%2Fresearch.microsoft.com%2Fpubs%2F135344%2Fmsr-tr-2011-18.pdf

  Note that we constraint the generated worlds so they satisfies the
  observation constraints (a > b > c).

  For num_people = 5:

Summary Statistics
      parameters       mean       std   naive_se      mcse        ess      rhat   ess_per_sec
          Symbol    Float64   Float64    Float64   Float64    Float64   Float64       Float64

        skill[1]   102.2371    2.8288     0.0283    0.1671   231.3914    0.9999        5.5389
        skill[2]   100.8802    2.5754     0.0258    0.1498   281.8016    1.0053        6.7455
        skill[3]    99.7618    2.5047     0.0250    0.1523   215.7630    1.0022        5.1648
        skill[4]    99.1769    2.7811     0.0278    0.1690   206.5924    1.0058        4.9452
        skill[5]    97.8249    2.8215     0.0282    0.1792   185.3145    1.0026        4.4359
  performance[1]   105.4878    3.5954     0.0360    0.2291   198.0813    1.0005        4.7415
  performance[2]   102.3919    2.5597     0.0256    0.1493   184.1535    1.0264        4.4081
  performance[3]   100.0782    2.5635     0.0256    0.1485   201.1306    1.0179        4.8145
  performance[4]    97.8592    2.7721     0.0277    0.1637   196.4140    1.0023        4.7016
  performance[5]    94.4979    3.2780     0.0328    0.2024   158.5955    1.0044        3.7963

skill[1] > skill[2]: 0.6344
skill[1] > skill[3]: 0.7565
skill[1] > skill[4]: 0.7846
skill[1] > skill[5]: 0.8963
skill[2] > skill[3]: 0.6219
skill[2] > skill[4]: 0.7113
skill[2] > skill[5]: 0.7706
skill[3] > skill[4]: 0.5337
skill[3] > skill[5]: 0.6878
skill[4] > skill[5]: 0.6122
performance[1] > performance[2]:1.0
performance[1] > performance[3]:1.0
performance[1] > performance[4]:1.0
performance[1] > performance[5]:1.0
performance[2] > performance[3]:0.9999
performance[2] > performance[4]:0.9999
performance[2] > performance[5]:0.9995
performance[3] > performance[4]:0.9995
performance[3] > performance[5]:0.9992
performance[4] > performance[5]:0.9992

  Cf: ~/cplint/trueskill.pl
      ~/blog/true_skill.blog
      ~/psi/true_skill.psi
      ~/webppl/true_skill.wpl

   Compare with tug_of_war.jl

=#

using Turing, StatsPlots, DataFrames
include("jl_utils.jl")

@model function true_skill(num_people)
    # Each person has an unknown (latent) Skill and a
    # known performance, where the skill is
    # reflected in the performance (with uncertainties).
    skill       = tzeros(num_people)
    performance = tzeros(num_people)
    for p in 1:num_people
        skill[p] ~ Normal(100,sqrt(10))
        performance[p] ~ Normal(skill[p], sqrt(15))
    end

    # Nicer (and slower), but it's not supported by IS or PG!
    # skill .~ Normal(100,sqrt(10))
    # @. performance ~ Normal(skill, sqrt(15)) # Note The @. since we refer to skill[p]
    for p1 in 1:num_people-1, p2 in p1+1:num_people
        true ~ Dirac(performance[p1] > performance[p2])
    end

    # What are their performance and (underlying) skills?

end

# num_people = 3
num_people = 5
model = true_skill(num_people)

num_chns = 4

# chns = sample(model, Prior(), 10_000)
# chns = sample(model, MH(), 10_000)
chns = sample(model, PG(15), 10_000)
# chns = sample(model, SMC(), 10_000)
# chns = sample(model, IS(), 10_000)

display(chns)

# display(plot(chns))
# display(histogram(chns))
# display(gelmandiag(chns)) # requires number of chns >= 2

df = DataFrame(chns);

# Show the latent skills
for p1 in 1:num_people-1, p2 in p1+1:num_people
    println("skill[$p1] > skill[$p2]: ", mean(df[!,"skill[$p1]"] .>df[!,"skill[$p2]"]))
end

# Show the observed performance
for p1 in 1:num_people-1, p2 in p1+1:num_people
    println("performance[$p1] > performance[$p2]:", mean(df[!,"performance[$p1]"] .>df[!,"performance[$p2]"]))
end

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
#primary structure for portfolio math
#some auxillery structures to improve missing value efficiency
struct Portfolio
  U::Universe
  wₜ::Vector{Float64}
  R::Vector{Float64}
  pμ::Base.RefValue{Float64}

  #the below vector contains all valid weights, and does not include invalid weights
  #thus the shape of wvec is a ragged array, with each array of length equal to the
  #number of valid values
  wvec::Vector{Vector{Float64}} #store this so we don't have to reallocate
end

#update the portfolio object
#innermost loop and likely a performance bottleneck, so highly optimized
function Portfolio!(P::Portfolio, wₜ::AbstractVector; normalizeto::Float64 = 1.0,
  checkadjustment::Bool = true)::Portfolio

  P.wₜ .= wₜ

  for c ∈ 1:P.U.N
    for r ∈ P.U.missingindices[c]
      P.wvec[c][r] = wₜ[r]
    end

    adjfactor::Float64 = sum(P.wvec[c]) / normalizeto
    P.wvec[c] ./= adjfactor
    checkadjustment && (adjfactor < .05) && @warn("Extreme adjustment factor $adjfactor -> scaling =
      $(1/adjfactor). Expect unstable results. Consider shifting the test distribution more
      positive or increasing data completeness.")

  end

  #build the returns (Note: Threads.@threads provides signficant performance gains)
  P.R .= 0.0
  for c ∈ 1:P.U.N
    for r ∈ P.U.missingindices[c]
      P.R[c] += P.wvec[c][r] * P.U.returns[r,c]
    end
  end

  #update the pointer
  P.pμ[] = mean(P.R)
  return P
end

#construct a brand new portfolio object
function Portfolio(U::Universe, wₜ::AbstractVector)::Portfolio
  local wvec::Vector{Vector{Float64}} = (n->zeros(U.M)).(1:U.N)

  P::Portfolio = Portfolio( #pre-allocate and call the update function
    U,
    Vector{Float64}(undef, U.M),
    Vector{Float64}(undef, U.N),
    Ref(-1.0),
    wvec)

  return Portfolio!(P, wₜ)
end

#helper method which creates an equi-weighted portoflio
Portfolio(U::Universe)::Portfolio = Portfolio(U, fill(1.0/U.M, U.M))

#creates a new portfolio from another one, re-using some of the references
function Portfolio(sourceP::Portfolio, wₜ::AbstractVector)

  P::Portfolio = Portfolio(sourceP.U,
    Vector{Float64}(undef, sourceP.U.M),
    Vector{Float64}(undef, sourceP.U.N),
    Ref(-1.0),
    sourceP.wvec)

  Portfolio!(P, wₜ)
  return P
end

#draws a random test portfolio using the universe from G
function redrawtestportfolio!(P::Portfolio;numtestassets::Int = 100,
  randomweightdist::T = RANDOM_WEIGHT_DIST,
  weightoftestasset::Function = i::Int -> rand(randomweightdist), #easy to modify if needed
  minadjustment::Float64 = MIN_ADJUSTMENT
  )::Portfolio where T<:Union{Distribution, Vector}

  wₜ::Vector{Float64} = zeros(P.U.M)
  for i ∈ 1:numtestassets
    wₜ[i] = weightoftestasset(i)
  end

  tot::Float64 = sum(wₜ)

  #Don't want extreme weights or a flipped distribution (numerical instability)
  #Value of 0.75 effectively caps the max weight to 200%
  if tot ≤ minadjustment
    return redrawtestportfolio!(P)
  end
  wₜ ./= tot

  shuffle!(wₜ)
  try
    return Portfolio!(P, wₜ, checkadjustment=false)
  catch
    return redrawtestportfolio!(P)
  end
end

#initialize the test portfolio from G
function initializetestportfolio(G::Portfolio)::Portfolio
  P::Portfolio = Portfolio(G, G.wₜ)

  return redrawtestportfolio!(P)
end

function Statistics.cov(P1::Portfolio, P2::Portfolio)::Float64

  return dot(P1.R, P2.R) / P1.U.N - P1.pμ[]*P2.pμ[]
end

function Statistics.var(P1::Portfolio)::Float64
  return cov(P1,P1)
end

function cov!(portfolios::Vector{Portfolio},
  S::AbstractMatrix = Matrix{Float64}(undef, length(portfolios), length(portfolios)))

  #compute the covariance matrix
  for c ∈ 1:length(portfolios)
    for r ∈ 1:c
      S[r,c] = cov(portfolios[r], portfolios[c])
      if r ≠ c
        S[c,r] = S[r,c] #fill in the lower left triangle
      end
    end
  end

  return S
end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
# MIT License
#
# Copyright (c) 2018 Martin Biel
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

@everywhere module TestDecisionVariable

using StochasticPrograms
using Test

function test_decision_no_bound(Structure)
    ξ₁ = @scenario a = 1 probability = 0.5
    ξ₂ = @scenario a = 2 probability = 0.5
    sp = StochasticProgram([ξ₁,ξ₂], Structure...)
    @first_stage sp = begin
        @decision(sp, x)
        @variable(sp, w)
    end
    @second_stage sp = begin
        @variable(sp, z)
        @recourse(sp, y)
    end
    # First-stage
    x = sp[1,:x]
    @test StochasticPrograms.stage(x) == 1
    @test !JuMP.has_lower_bound(x)
    @test !JuMP.has_upper_bound(x)
    @test !JuMP.is_fixed(x)
    @test "x" == @inferred JuMP.name(x)
    @test x == decision_by_name(owner_model(x), 1, "x")
    # Scenario-dependent
    y = sp[2,:y]
    @test StochasticPrograms.stage(y) == 2
    @test_throws ErrorException JuMP.has_lower_bound(y)
    @test !JuMP.has_lower_bound(y, 1)
    @test !JuMP.has_lower_bound(y, 2)
    @test_throws ErrorException JuMP.has_upper_bound(y)
    @test !JuMP.has_upper_bound(y, 1)
    @test !JuMP.has_upper_bound(y, 2)
    @test_throws ErrorException JuMP.is_fixed(y)
    @test !JuMP.is_fixed(y, 1)
    @test !JuMP.is_fixed(y, 2)
    @test_throws ErrorException JuMP.is_integer(y)
    @test !JuMP.is_integer(y, 1)
    @test !JuMP.is_integer(y, 2)
    @test_throws ErrorException JuMP.is_binary(y)
    @test !JuMP.is_binary(y, 1)
    @test !JuMP.is_binary(y, 2)
    @test "y" == @inferred JuMP.name(y)
    @test y == decision_by_name(owner_model(y), 2, "y")
end

function test_decision_lower_bound(Structure)
    ξ₁ = @scenario a = 1 probability = 0.5
    ξ₂ = @scenario a = 2 probability = 0.5
    sp = StochasticProgram([ξ₁,ξ₂], Structure...)
    @first_stage sp = begin
        @decision(sp, x >= 0, Bin)
        @variable(sp, w)
    end
    @second_stage sp = begin
        @uncertain a
        @variable(sp, z)
        @recourse(sp, y >= a, Bin)
    end
    # First-stage
    x = sp[1,:x]
    @test JuMP.has_lower_bound(x)
    @test 0.0 == @inferred JuMP.lower_bound(x)
    @test !JuMP.has_upper_bound(x)
    @test !JuMP.is_fixed(x)
    @test JuMP.is_binary(x)
    @test !JuMP.is_integer(x)
    JuMP.delete_lower_bound(x)
    @test !JuMP.has_lower_bound(x)
    # Scenario-dependent
    y = sp[2,:y]
    @test JuMP.has_lower_bound(y, 1)
    @test JuMP.has_lower_bound(y, 2)
    @test 1.0 == @inferred JuMP.lower_bound(y, 1)
    @test 2.0 == @inferred JuMP.lower_bound(y, 2)
    @test !JuMP.has_upper_bound(y, 1)
    @test !JuMP.has_upper_bound(y, 2)
    @test !JuMP.is_fixed(y, 1)
    @test !JuMP.is_fixed(y, 2)
    @test JuMP.is_binary(y, 1)
    @test JuMP.is_binary(y, 2)
    @test !JuMP.is_integer(y, 1)
    @test !JuMP.is_integer(y, 2)
    JuMP.delete_lower_bound(y, 1)
    @test !JuMP.has_lower_bound(y, 1)
    JuMP.delete_lower_bound(y, 2)
    @test !JuMP.has_lower_bound(y, 2)
end

function test_decision_upper_bound(Structure)
    ξ₁ = @scenario a = 1 probability = 0.5
    ξ₂ = @scenario a = 2 probability = 0.5
    sp = StochasticProgram([ξ₁,ξ₂], Structure...)
    @first_stage sp = begin
        @decision(sp, x <= 1.0, Int)
        @variable(sp, w)
    end
    @second_stage sp = begin
        @uncertain a
        @variable(sp, z)
        @recourse(sp, y <= a, Int)
    end
    # First-stage
    x = sp[1,:x]
    @test !JuMP.has_lower_bound(x)
    @test JuMP.has_upper_bound(x)
    @test 1.0 == @inferred JuMP.upper_bound(x)
    @test !JuMP.is_fixed(x)
    @test !JuMP.is_binary(x)
    @test JuMP.is_integer(x)
    JuMP.delete_upper_bound(x)
    @test !JuMP.has_upper_bound(x)
    # Scenario-dependent
    y = sp[2,:y]
    @test !JuMP.has_lower_bound(y, 1)
    @test !JuMP.has_lower_bound(y, 2)
    @test JuMP.has_upper_bound(y, 1)
    @test JuMP.has_upper_bound(y, 2)
    @test 1.0 == @inferred JuMP.upper_bound(y, 1)
    @test 2.0 == @inferred JuMP.upper_bound(y, 2)
    @test !JuMP.is_fixed(y, 1)
    @test !JuMP.is_fixed(y, 2)
    @test !JuMP.is_binary(y, 1)
    @test !JuMP.is_binary(y, 2)
    @test JuMP.is_integer(y, 1)
    @test JuMP.is_integer(y, 2)
    JuMP.delete_upper_bound(y, 1)
    @test !JuMP.has_upper_bound(y, 1)
    JuMP.delete_upper_bound(y, 2)
    @test !JuMP.has_upper_bound(y, 2)
end

function test_decision_fix(Structure)
    ξ₁ = @scenario a = 1 probability = 0.5
    ξ₂ = @scenario a = 2 probability = 0.5
    sp = StochasticProgram([ξ₁,ξ₂], Structure...)
    @first_stage sp = begin
        @decision(sp, x == 1)
        @variable(sp, w)
    end
    @second_stage sp = begin
        @uncertain a
        @variable(sp, z)
        @recourse(sp, y == a)
    end
    # First-stage
    x = sp[1,:x]
    @test !JuMP.has_lower_bound(x)
    @test !JuMP.has_upper_bound(x)
    @test JuMP.is_fixed(x)
    @test 1.0 == @inferred JuMP.value(x)
    JuMP.unfix(x)
    @test !JuMP.is_fixed(x)
    JuMP.fix(x, 2.0)
    @test !JuMP.has_lower_bound(x)
    @test !JuMP.has_upper_bound(x)
    @test JuMP.is_fixed(x)
    @test 2.0 == @inferred JuMP.value(x)
    # Scenario-dependent
    y = sp[2,:y]
    @test !JuMP.has_lower_bound(y, 1)
    @test !JuMP.has_lower_bound(y, 2)
    @test !JuMP.has_upper_bound(y, 1)
    @test !JuMP.has_upper_bound(y, 2)
    @test JuMP.is_fixed(y, 1)
    @test JuMP.is_fixed(y, 2)
    @test 1.0 == @inferred JuMP.value(y, 1)
    @test 2.0 == @inferred JuMP.value(y, 2)
    JuMP.unfix(y, 1)
    @test !JuMP.is_fixed(y, 1)
    @test JuMP.is_fixed(y, 2)
    JuMP.unfix(y, 2)
    @test !JuMP.is_fixed(y, 2)
    JuMP.fix(y, 1, 2.0)
    @test !JuMP.has_lower_bound(y, 1)
    @test !JuMP.has_upper_bound(y, 1)
    @test JuMP.is_fixed(y, 1)
    @test !JuMP.is_fixed(y, 2)
    @test 2.0 == @inferred JuMP.value(y, 1)
    JuMP.fix(y, 2, 1.0)
    @test !JuMP.has_lower_bound(y, 2)
    @test !JuMP.has_upper_bound(y, 2)
    @test JuMP.is_fixed(y, 2)
    @test 1.0 == @inferred JuMP.value(y, 2)
end

function test_decision_custom_index_sets(Structure)
    ξ₁ = @scenario a = 1 probability = 0.5
    ξ₂ = @scenario a = 2 probability = 0.5
    sp = StochasticProgram([ξ₁,ξ₂], Structure...)
    @first_stage sp = begin
        @decision(sp, x[0:1, 10:20, 1:1] >= 2)
        @variable(sp, w)
    end
    @second_stage sp = begin
        @uncertain a
        @variable(sp, z)
        @recourse(sp, y[0:1, 10:20, 1:1] >= a)
    end
    # First-stage
    x = sp[1,:x]
    @test JuMP.has_lower_bound(x[0, 15, 1])
    @test 2 == @inferred JuMP.lower_bound(x[0, 15, 1])
    @test !JuMP.has_upper_bound(x[0, 15, 1])
    @first_stage sp = begin
        @decision(sp, x[i in -10:10, s in [:a,:b]] <= 5.5, Int)
        @variable(sp, w)
    end
    generate!(sp)
    x = sp[1,:x]
    @test 5.5 == @inferred JuMP.upper_bound(x[-4, :a])
    @test "x[-10,a]" == @inferred JuMP.name(x[-10,:a])
    @test x[-10,:a] == decision_by_name(sp, 1, "x[-10,a]")
    # Scenario-dependent
    y = sp[2,:y]
    @test JuMP.has_lower_bound(y[0, 15, 1], 1)
    @test JuMP.has_lower_bound(y[0, 15, 1], 2)
    @test 1 == @inferred JuMP.lower_bound(y[0, 15, 1], 1)
    @test 2 == @inferred JuMP.lower_bound(y[0, 15, 1], 2)
    @test !JuMP.has_upper_bound(y[0, 15, 1], 1)
    @test !JuMP.has_upper_bound(y[0, 15, 1], 2)
    @second_stage sp = begin
        @variable(sp, z)
        @recourse(sp, y[i in -10:10, s in [:a,:b]] <= 5.5, Int)
    end
    y = sp[2,:y]
    @test 5.5 == @inferred JuMP.upper_bound(y[-4, :a], 1)
    @test 5.5 == @inferred JuMP.upper_bound(y[-4, :a], 2)
    @test "y[-10,a]" == @inferred JuMP.name(y[-10,:a])
    @test nothing === decision_by_name(sp, 1, "y[-10,a]")
    @test y[-10,:a] == decision_by_name(sp, 2, "y[-10,a]")
end

function test_variable_is_valid_delete(Structure)
    ξ₁ = @scenario a = 1 probability = 0.5
    ξ₂ = @scenario a = 2 probability = 0.5
    sp = StochasticProgram([ξ₁,ξ₂], Structure...)
    @first_stage sp = begin
        @decision(sp, x)
        @variable(sp, w)
    end
    @second_stage sp = begin
        @variable(sp, z)
        @recourse(sp, y)
    end
    # First-stage
    x = sp[1,:x]
    @test JuMP.is_valid(sp, x)
    JuMP.delete(sp, x)
    @test !JuMP.is_valid(sp, x)
    @test_throws Exception JuMP.delete(sp, x)
    # Scenario-dependent
    y = sp[2,:y]
    @test_throws ErrorException JuMP.is_valid(sp, y)
    @test JuMP.is_valid(sp, y, 1)
    @test JuMP.is_valid(sp, y, 2)
    @test_throws ErrorException JuMP.delete(sp, y)
    JuMP.delete(sp, y, 1)
    @test !JuMP.is_valid(sp, y, 1)
    @test JuMP.is_valid(sp, y, 2)
    JuMP.delete(sp, y, 2)
    @test !JuMP.is_valid(sp, y, 2)
    @test_throws Exception JuMP.delete(sp, y, 1)
    @first_stage sp = begin
        @decision(sp, x[1:3] >= 1)
        @variable(sp, w)
    end
    @second_stage sp = begin
        @variable(sp, z)
        @recourse(sp, y[1:3] >= 1)
    end
    # First-stage
    x = sp[1,:x]
    @test all(is_valid.(sp, x))
    delete(sp, x)
    @test all((!is_valid).(sp, x))
    second_sp = StochasticProgram([Scenario()], Structure...)
    @test_throws Exception JuMP.delete(second_sp, x[1])
    @test_throws Exception JuMP.delete(second_sp, x)
    # Scenario-dependent
    y = sp[2,:y]
    @test_throws ErrorException is_valid.(sp, y)
    @test all(is_valid.(sp, y, 1))
    @test all(is_valid.(sp, y, 2))
    @test_throws ErrorException delete.(sp, y)
    delete(sp, y, 1)
    @test all((!is_valid).(sp, y, 1))
    @test all(is_valid.(sp, y, 2))
    delete(sp, y, 2)
    @test all((!is_valid).(sp, y, 2))
end

function test_variable_bounds_set_get(Structure)
    ξ₁ = @scenario a = 1 probability = 0.5
    ξ₂ = @scenario a = 2 probability = 0.5
    sp = StochasticProgram([ξ₁,ξ₂], Structure...)
    @first_stage sp = begin
        @decision(sp, 0 <= x <= 2)
        @variable(sp, w)
    end
    @second_stage sp = begin
        @uncertain a
        @variable(sp, z)
        @recourse(sp, 0 <= y <= a)
    end
    # First-stage
    x = sp[1,:x]
    @test 0 == @inferred JuMP.lower_bound(x)
    @test 2 == @inferred JuMP.upper_bound(x)
    set_lower_bound(x, 1.)
    @test 1. == @inferred JuMP.lower_bound(x)
    set_upper_bound(x, 3.)
    @test 3. == @inferred JuMP.upper_bound(x)
    # Scenario-dependent
    y = sp[2,:y]
    @test_throws ErrorException JuMP.lower_bound(y)
    @test 0 == @inferred JuMP.lower_bound(y, 1)
    @test 0 == @inferred JuMP.lower_bound(y, 2)
    @test_throws ErrorException JuMP.upper_bound(y)
    @test 1 == @inferred JuMP.upper_bound(y, 1)
    @test 2 == @inferred JuMP.upper_bound(y, 2)
    set_lower_bound(y, 1, 1.)
    @test 1. == @inferred JuMP.lower_bound(y, 1)
    set_lower_bound(y, 2, 1.)
    @test 1. == @inferred JuMP.lower_bound(y, 2)
    JuMP.set_lower_bound(y, 2.)
    @test 2. == @inferred JuMP.lower_bound(y, 1)
    @test 2. == @inferred JuMP.lower_bound(y, 2)
    set_upper_bound(y, 1, 3.)
    @test 3. == @inferred JuMP.upper_bound(y, 1)
    set_upper_bound(y, 2, 3.)
    @test 3. == @inferred JuMP.upper_bound(y, 2)
    JuMP.set_upper_bound(y, 2.)
    @test 2. == @inferred JuMP.upper_bound(y, 1)
    @test 2. == @inferred JuMP.upper_bound(y, 2)
end

function test_variable_starts_set_get(Structure)
    ξ₁ = @scenario a = 1 probability = 0.5
    ξ₂ = @scenario a = 2 probability = 0.5
    sp = StochasticProgram([ξ₁,ξ₂], Structure...)
    @first_stage sp = begin
        @decision(sp, x[1:3])
        @variable(sp, w)
    end
    @second_stage sp = begin
        @variable(sp, z)
        @recourse(sp, y[1:3])
    end
    # First-stage
    x = sp[1,:x]
    x0 = collect(1:3)
    JuMP.set_start_value.(x, x0)
    @test JuMP.start_value.(x) == x0
    @test JuMP.start_value.([x[1],x[2],x[3]]) == x0
    # First-stage
    y = sp[2,:y]
    @test_throws ErrorException JuMP.set_start_value.(y, x0)
    JuMP.set_start_value.(y, 1, x0)
    JuMP.set_start_value.(y, 2, x0)
    @test_throws ErrorException JuMP.start_value.(y)
    @test JuMP.start_value.(y, 1) == x0
    @test JuMP.start_value.(y, 2) == x0
    @test_throws ErrorException JuMP.start_value.([y[1],y[2],y[3]])
    @test JuMP.start_value.([y[1],y[2],y[3]], 1) == x0
    @test JuMP.start_value.([y[1],y[2],y[3]], 2) == x0
end

function test_variable_integrality_set_get(Structure)
    ξ₁ = @scenario a = 1 probability = 0.5
    ξ₂ = @scenario a = 2 probability = 0.5
    sp = StochasticProgram([ξ₁,ξ₂], Structure...)
    @first_stage sp = begin
        @decision(sp, x[1:3])
        @variable(sp, w)
    end
    @second_stage sp = begin
        @variable(sp, z)
        @recourse(sp, y[1:3])
    end
    # First-stage
    x = sp[1,:x]
    JuMP.set_integer(x[2])
    JuMP.set_integer(x[2])
    @test JuMP.is_integer(x[2])
    JuMP.unset_integer(x[2])
    @test !JuMP.is_integer(x[2])
    JuMP.set_binary(x[1])
    JuMP.set_binary(x[1])
    @test JuMP.is_binary(x[1])
    @test_throws Exception JuMP.set_integer(x[1])
    JuMP.unset_binary(x[1])
    @test !JuMP.is_binary(x[1])
    # Scenario-dependent
    y = sp[2,:y]
    @test_throws ErrorException JuMP.set_integer(y[2])
    JuMP.set_integer(y[2], 1)
    JuMP.set_integer(y[2], 1)
    JuMP.set_integer(y[3], 2)
    JuMP.set_integer(y[3], 2)
    @test_throws ErrorException JuMP.is_integer(y[2])
    @test JuMP.is_integer(y[2], 1)
    @test !JuMP.is_integer(y[2], 2)
    @test JuMP.is_integer(y[3], 2)
    @test !JuMP.is_integer(y[3], 1)
    @test_throws ErrorException JuMP.unset_integer(y[2])
    JuMP.unset_integer(y[2], 1)
    JuMP.unset_integer(y[3], 2)
    @test !JuMP.is_integer(y[2], 1)
    @test_throws ErrorException JuMP.set_binary(y[2])
    JuMP.set_binary(y[1], 1)
    JuMP.set_binary(y[1], 1)
    JuMP.set_binary(y[1], 2)
    JuMP.set_binary(y[1], 2)
    @test_throws ErrorException JuMP.is_binary(y[1])
    @test JuMP.is_binary(y[1], 1)
    @test JuMP.is_binary(y[1], 2)
    @test_throws Exception JuMP.set_integer(y[1], 1)
    @test_throws Exception JuMP.set_integer(y[1], 2)
    @test_throws ErrorException JuMP.unset_binary(y[1])
    JuMP.unset_binary(y[1], 1)
    @test !JuMP.is_binary(y[1], 1)
    JuMP.unset_binary(y[1], 2)
    @test !JuMP.is_binary(y[1], 2)
end

function test_variables_constrained_on_creation(Structure)
    ξ₁ = @scenario a = 1 probability = 0.5
    ξ₂ = @scenario a = 2 probability = 0.5
    sp = StochasticProgram([ξ₁,ξ₂], Structure...)
    @first_stage sp = begin
        @decision(sp, x[1:2] in SecondOrderCone())
        @variable(sp, w)
    end
    @second_stage sp = begin
        @variable(sp, r)
        @recourse(sp, y[1:2] in SecondOrderCone())
    end
    x = sp[1,:x]
    y = sp[2,:y]
    @test num_constraints(sp, 1, typeof(x), MOI.SecondOrderCone) == 1
    @test num_constraints(sp, 2, typeof(y), MOI.SecondOrderCone) == 1
    @test name(x[1]) ==  "x[1]"
    @test name(x[2]) ==  "x[2]"
    @test name(y[1]) == "y[1]"
    @test name(y[2]) == "y[2]"
    @first_stage sp = begin
        @decision(sp, x[1:2] in SecondOrderCone())
        @decision(sp, [1:2] in SecondOrderCone())
        @variable(sp, w)
    end
    @second_stage sp = begin
        @variable(sp, r)
        @recourse(sp, y[1:2] in SecondOrderCone())
        @recourse(sp, [1:2] in SecondOrderCone())
    end
    x = sp[1,:x]
    y = sp[2,:y]
    @test num_constraints(sp, 1, typeof(x), MOI.SecondOrderCone) == 2
    @test num_constraints(sp, 2, typeof(y), MOI.SecondOrderCone) == 2
    @first_stage sp = begin
        @decision(sp, x[1:2] in SecondOrderCone())
        @decision(sp, [1:2] in SecondOrderCone())
        @decision(sp, [1:3] in SecondOrderCone())
        @variable(sp, w)
    end
    @second_stage sp = begin
        @variable(sp, r)
        @recourse(sp, y[1:2] in SecondOrderCone())
        @recourse(sp, [1:2] in SecondOrderCone())
        @recourse(sp, [1:3] in SecondOrderCone())
    end
    x = sp[1,:x]
    y = sp[2,:y]
    @test num_constraints(sp, 1, typeof(x), MOI.SecondOrderCone) == 3
    @test num_constraints(sp, 2, typeof(y), MOI.SecondOrderCone) == 3
    @first_stage sp = begin
        @decision(sp, z in MOI.Semiinteger(1.0, 2.0))
        @variable(sp, w)
    end
    generate!(sp)
    z = sp[1,:z]
    @test num_constraints(sp, 1, typeof(z), MOI.Semiinteger{Float64}) == 1
    @first_stage sp = begin
        @variable(sp, w)
        @decision(sp, z in MOI.Semiinteger(1.0, 2.0))
        @decision(sp, set = MOI.Semiinteger(1.0, 2.0))
    end
    generate!(sp)
    z = sp[1,:z]
    @test num_constraints(sp, 1, typeof(z), MOI.Semiinteger{Float64}) == 2
end

function test_all_decision_variables(Structure)
    ξ₁ = @scenario a = 1 probability = 0.5
    ξ₂ = @scenario a = 2 probability = 0.5
    sp = StochasticProgram([ξ₁,ξ₂], Structure...)
    @first_stage sp = begin
        @decision(sp, x₁)
        @decision(sp, x₂)
        @variable(sp, w)
    end
    @second_stage sp = begin
        @variable(sp, z)
        @recourse(sp, y₁)
        @recourse(sp, y₂)
    end
    # First-stage
    x₁ = sp[1,:x₁]
    x₂ = sp[1,:x₂]
    @test [x₁, x₂] == @inferred all_decision_variables(sp, 1)
    # Scenario-dependent
    y₁ = sp[2,:y₁]
    y₂ = sp[2,:y₂]
    @test [y₁, y₂] == @inferred all_decision_variables(sp, 2)
end

function runtests()
    @testset "DecisionVariable" begin
        for structure in [(Deterministic(),),
                          (StageDecomposition(),),
                          (ScenarioDecomposition(),),
                          (Deterministic(), () -> MOIU.MockOptimizer(MOIU.Model{Float64}()))]
            name = length(structure) == 1 ? "$(structure[1])" : "$(structure[1]) with decision bridges"
            @testset "$name" begin
                for name in names(@__MODULE__; all = true)
                    if !startswith("$(name)", "test_")
                        continue
                    end
                    f = getfield(@__MODULE__, name)
                    @testset "$(name)" begin
                        f(structure)
                    end
                end
            end
        end
    end
end

function run_dtests()
    @testset "DecisionVariable" begin
        for structure in [(DistributedStageDecomposition(),), (DistributedScenarioDecomposition(),)]
            @testset "$(structure[1])" begin
                for name in names(@__MODULE__; all = true)
                    if !startswith("$(name)", "test_")
                        continue
                    end
                    f = getfield(@__MODULE__, name)
                    @testset "$(name)" begin
                        f(structure)
                    end
                end
            end
        end
    end
end

end


text:
using DataFrames
using MySQL
using BioServices.EUtils
using XMLDict
using LightXML
import Base.parse

#------------------ BioMedQuery -------------------
    @testset "Testing Eutils/PubMed" begin
    #testset "globals"
    narticles = 10
    ids = Array{Int64,1}()
    efetch_doc = ""
    dfs_efetch = Dict{String,DataFrame}()
    verbose = false
    articles = []

    @testset "Testing ESearch/EFetch PubMed" begin
        println("-----------------------------------------")
        println("       Testing ESearch/EFetch for PubMed     ")

        search_term="""(obstructive sleep apnea[MeSH Major Topic]) AND (journal article[Publication Type])"""

        esearch_response = esearch(db="pubmed", term = search_term,
        retstart = 0, retmax = narticles, tool ="BioJulia")

        #convert xml to dictionary
        esearch_dict = parse_xml(String(esearch_response.body))

        #examine how many ids were returned
        @test haskey(esearch_dict, "IdList")

        for id_node in esearch_dict["IdList"]["Id"]
            push!(ids, Base.parse(Int64, id_node))
        end

        @test length(ids)==narticles

        efetch_response = efetch(db = "pubmed", tool = "BioJulia", retmode = "xml", rettype = "null", id = ids)

        #convert xml to dictionary
        efetch_doc = root(parse_string(String(efetch_response.body)))
        dfs_efetch = PubMed.parse_articles(efetch_doc)

        @test name(efetch_doc) == "PubmedArticleSet"

        #articles should be an array of narticles
        @test narticles == length(collect(child_elements(efetch_doc)))

    end

    @testset "Test Save PMID MySQL" begin
        println("-----------------------------------------")
        println("       Test Save PMID MySQL     ")
        dbname = "entrez_test"
        host = "127.0.0.1";
        user = "root"
        pwd = ""

        conn = DBUtils.init_mysql_database(host, user, pwd, dbname)
        PubMed.create_pmid_table!(conn)
        PubMed.save_pmids!(conn, ids)

        #query the article table and make sure the count is correct
        all_pmids = PubMed.all_pmids(conn)
        @test length(all_pmids) == narticles

        #clean-up
        db_query(conn, "DROP DATABASE IF EXISTS $dbname;")

    end

    @testset "Testing Citations" begin
        println("-----------------------------------------")
        println("       Testing Citations      ")
        #parse 1st article
        # art = BioMedQuery.PubMed.MedlineArticle(articles[1])
        # println(art)
        citation = PubMed.CitationOutput("endnote", "./citations_temp.endnote", true)
        nsucceses = PubMed.save_efetch!(citation, efetch_doc, verbose)

        #test that citations are the same as the ones already stored
        lines=[]
        open("./citations_temp.endnote") do f
           lines = readlines(f)
        end

        nlines = 0
        for li=1:length(lines)
            if occursin("%0 Journal Article", lines[li])
              nlines+=1
            end
        end

        @test nlines == nsucceses
        rm("./citations_temp.endnote")
    end

    # save the results of an entrez fetch to a sqlite database
    @testset "Testing SQLite Saving" begin
        println("-----------------------------------------")
        println("       Testing SQLite Saving")

        db_path = "./test_db.db"

        conn = SQLite.DB(db_path)
        PubMed.create_tables!(conn)
        PubMed.save_efetch!(conn, efetch_doc,false, true)

        # query the df and db tables to make sure the count is correct
        for (table, df) in dfs_efetch
            query_res = db_query(conn, "select count(*) from $table")
            @test size(df)[1] == query_res[1][1]
        end

        #query the article table and make sure the count is correct
        all_pmids = PubMed.all_pmids(conn)
        @test length(all_pmids) == narticles

        #query the article table and make sure the count is correct
        all_abstracts = PubMed.abstracts(conn)
        @test size(all_abstracts)[1] >0 && size(all_abstracts)[1] <= narticles

        #check we can get the MESH descriptor for an article
        mesh = PubMed.get_article_mesh(conn, all_pmids[1])
        @test length(mesh) > 0

        #check that remainder of tables are not empty
        tables = ["author_ref", "mesh_desc",
        "mesh_qual", "mesh_heading", "pub_type", "abstract_structured"]

        for t in tables
            q = SQLite.query(conn, "SELECT count(*) FROM "*t*";")
            count = q[1][1]
            @test count > 0
        end

        # remove temp files
        if isfile(db_path)
            rm(db_path)
        end
    end

    @testset "Testing MySQL Saving" begin
        println("-----------------------------------------")
        println("       Testing MySQL Saving")

        dbname = "efetch_test"
        host = "127.0.0.1";
        user = "root"
        pwd = ""

        conn = DBUtils.init_mysql_database(host, user, pwd, dbname)
        PubMed.create_tables!(conn)
        @time PubMed.save_efetch!(conn, efetch_doc, false, true)

        # query the df and db tables to make sure the count is correct
        for (table, df) in dfs_efetch
            query_res = db_query(conn, "select count(*) from $table")
            @test size(df)[1] == query_res[1][1]
        end

        #query the article table and make sure the count is correct
        all_pmids = PubMed.all_pmids(conn)
        @test length(all_pmids) == narticles

        #query the article table and make sure the count is correct
        all_abstracts = PubMed.abstracts(conn)
        @test size(all_abstracts)[1] >0 && size(all_abstracts)[1] <= narticles

        #check we can get the MESH descriptor for an article
        mesh = PubMed.get_article_mesh(conn, all_pmids[1])
        # println(mesh)
        @test length(mesh) > 0

        #check that reminder of tables are not empty
        tables = ["author_ref", "mesh_desc",
        "mesh_qual", "mesh_heading", "pub_type", "abstract_structured"]

        for t in tables
            q = MySQL.query(conn, "SELECT count(*) FROM "*t*";")
            count = q[1][1]
            @test count > 0
        end

        #clean-up
        db_query(conn, "DROP DATABASE IF EXISTS $dbname;")

    end

    println("------------End Test Eutils/PubMed--------------")
    println("-----------------------------------------")

end


text:
struct TryCatch; end
register(TryCatch, Deprecation(
    "`try` without `catch` has been deprecated. Add a catch block.",
    "julia",
    v"0.1.0", v"0.7.0-DEV", typemax(VersionNumber)
))

match(TryCatch, CSTParser.Try) do x
    # Match ccall only
    dep, expr, resolutions, context = x
    # Don't match try blocks that already have a catch for finally clause
    for c in children(expr)
        if isexpr(c, CSTParser.KEYWORD, Tokens.FINALLY) ||
           isexpr(c, CSTParser.KEYWORD, Tokens.CATCH)
           return
        end
    end

    # Look for trivia between the END keyword and the preceeding node
    end_kw = children(expr)[end]
    ws_in_between = string(prev_node_ws(end_kw), leading_trivia(end_kw))
    rn = ChildReplacementNode(nothing, children(expr)[1:end-1], expr)
    if !('\n' in ws_in_between)
        # Add the `catch; ` right before the end
        push!(children(rn), ReplacementNode("catch", leading_trivia(end_kw),"; "))
        push!(children(rn), TriviaReplacementNode(rn, end_kw, "", trailing_trivia(end_kw)))
    else
        # Add `catch\n` on the line before the end with the same indent
        idx = findlast(ws_in_between, '\n')
        indent = ws_in_between[idx:end]
        push!(children(rn), ReplacementNode("catch", leading_trivia(end_kw),""))
        push!(children(rn), TriviaReplacementNode(rn, end_kw, indent, trailing_trivia(end_kw)))
    end
    buf = IOBuffer()
    print_replacement(buf, rn, false, false)
    push!(resolutions, TextReplacement(dep, expr.span, String(take!(buf))))
end


text:
module Utils
    export secondsToSamples, samplesToSeconds, ∫

    using FFTW
    function ∫(xs, ys)
        ∫ydx = 0
        for i = 1:(length(ys)-1)
            ∫ydx += (ys[i+1] + ys[i]) * (xs[i+1] - xs[i]) / 2
        end
        return ∫ydx
    end

    function spectrum(window)
        S = fft(window)
        S = S[1:Int(length(S)/2)]
        scale = 2 / length(window)
        return [sqrt(real(n)^2 + imag(n)^2) * scale for n in S]
    end

    function hannWindow(window)
        N = length(window)
        return [window[n] * 0.5 * (1 - cos(2π * n / N)) for n = 1:N]
    end

    function secondsToSamples(sec, Fs)
        return Int64(round(sec * Fs))
    end

    function samplesToSeconds(samples, Fs)
        return samples / Fs
    end

    function frequencyToSamples(freq, Fs; roundFunc = round)
        samples = Fs / freq
        return Int(roundFunc(samples))
    end

    function samplesToFrequency(samples, Fs)
        return Fs / samples
    end

    function getWindows(y, Wl; hop_size = Wl)
        num_steps = Int(ceil(length(y) / hop_size - 1))
        windows = [zeros(Wl) for i in 1:num_steps]
        for i in 1:num_steps
            x_start = (i - 1) * hop_size + 1
            copyto!(windows[i], y[x_start:min(x_start + Wl - 1, end)])
        end
        return windows
    end

    function binIdx(freq, bin_size; roundFunc::Function = round)
        return Int(roundFunc(freq / bin_size)) + 1
    end

    function idxToBin(idx, bin_size)
        return (idx - 1) * bin_size
    end

    function closestBin(freq, bins)
        i = 1
        while bins[i] < freq
            i += 1
        end
        if i > 1 && (abs(bins[i] - freq) > abs(bins[i-1] - freq))
            i -= 1
        end
        @assert i >= 1
        return i
    end

    function closestPowerOf2(n)
        i = 0
        while 2^i < n
            i += 1
        end

        if abs(2^i - n) > abs(2^(i - 1) - n)
            return 2^(i - 1)
        else
            return 2^i
        end
    end

    function calculateCents(f1, f2)
        high = max(f1, f2)
        low = min(f1, f2)
        cents = Int(floor(1200 * log(2, high / low)))
        return cents
    end

    function getNoteFromFreq(freq)
        if freq === nothing
            return nothing
        end
        noteNum = 12 * log(2, freq / 440) + 49
        trueNote = Int(floor(noteNum))
        error = calculateCents(getNote(trueNote), freq) / 100

        return trueNote + error
    end

    function getNoteName(noteNum)
        trueNote = Int(floor(noteNum))
        error = noteNum - trueNote

        octave = Int(floor((trueNote + 8) / 12))
        noteName = ((trueNote + 8) % 12) + 1
        if noteName < 1
            noteName += 12
        end
        noteName = ["C", "C#", "D", "D#", "E", "F", "F#", "G", "G#", "A", "A#", "B"][noteName]

        str = string(noteName, octave)
        if error != 0
            str = string(str, " + ", Int(round(error * 100)), "c")
        end
        return (str)
    end

    function getNote(noteNum::Integer)
        return 2^((noteNum - 49) // 12) * 440
    end

    function getRange(low_note, high_note)
        i = 0
        while getNote(i) < low_note
            i += 1
        end

        @assert getNote(i) < high_note

        j = i
        while getNote(j) < high_note
            j += 1
        end

        return getNote.(i:j)
    end

    function findpartialsort(v, k; opts...)
        inds = [partialsortperm(v, k; opts...)...]
        return v[[inds...]], inds
    end

    Base.findmax(v::AbstractArray, k::AbstractUnitRange; opts...) = findpartialsort(v, k; opts..., rev = true)
    Base.findmin(v::AbstractArray, k::AbstractUnitRange; opts...) = findpartialsort(v, k; opts...)
end


text:
"""
Register a function pointer as an R NativeSymbol.

This is completely undocumented, so may break: we technically are supposed to
use R_registerRoutines, but this is _much_ easier for just 1 function.
"""
function makeNativeSymbolRef(fptr::Ptr{Void})
    # mirror Rf_MakeNativeSymbolRef of Rdynload.c
    rexfn = ccall((:R_MakeExternalPtrFn,libR), ExtPtrSxpPtr,
                     (Ptr{Void}, Ptr{Void}, Ptr{Void}),
                     fptr, sexp(Symbol("native symbol")), sexp(Const.NilValue))
    setattrib!(rexfn, Const.ClassSymbol, sexp("NativeSymbol"))
    preserve(rexfn)
    rexfn
end

"Create an ExtPtrSxpPtr object"
makeExternalPtr(ptr::Ptr{Void}, tag=Const.NilValue, prot=Const.NilValue) =
    ccall((:R_MakeExternalPtr,libR), ExtPtrSxpPtr,
          (Ptr{Void}, UnknownSxpPtr, UnknownSxpPtr),
          ptr, sexp(tag), sexp(prot))

"""
The function called by R .External for Julia callbacks.

It receives a `ListSxpPtr` containing
 - a pointer to the function itself (`ExtPtrSxpPtr`)
 - a pointer to the Julia function (`ExtPtrSxpPtr`)
 - any arguments (as `SxpPtr`)
"""
function julia_extptr_callback(p::ListSxpPtr)
    try
        l = cdr(p) # skip callback pointer

        # julia function pointer
        f_sxp = car(l)::ExtPtrSxpPtr
        f_sxprec = unsafe_load(f_sxp)
        f = unsafe_pointer_to_objref(f_sxprec.ptr)
        l = cdr(l)

        # # extract arguments
        args = Any[]
        kwargs = Any[]
        for (k,a) in l
            # TODO: provide a mechanism for users to specify their own
            # conversion routines
            if k == sexp(Const.NilValue)
                push!(args,rcopy(a))
            else
                push!(kwargs,(rcopy(Symbol,k),rcopy(a)))
            end
        end

        # call function
        y = f(args...;kwargs...)

        # return appropriate sexp
        return p = convert(UnknownSxpPtr,sexp(y))::UnknownSxpPtr
    catch e
        ccall((:Rf_error,libR),Ptr{Void},(Ptr{Cchar},),string(e))
        return convert(UnknownSxpPtr,sexp(Const.NilValue))::UnknownSxpPtr
    end
end

"""
Julia types (typically functions) which are wrapped in `ExtPtrSxpPtr` are
stored here to prevent garbage collection by Julia.
"""
const jtypExtPtrs = Dict{ExtPtrSxpPtr, Any}()

"""
Called by the R finalizer.
"""
function decref_extptr(p::ExtPtrSxpPtr)
    delete!(jtypExtPtrs, p)
    return nothing
end

const juliaDecref = Ref{Ptr{Void}}()

"""
Register finalizer to be called by the R GC.
"""
function registerCFinalizerEx(s::ExtPtrSxpPtr)
    protect(s)
    ccall((:R_RegisterCFinalizerEx,libR),Void,
          (Ptr{ExtPtrSxp}, Ptr{Void}, Cint),
          s,juliaDecref[],0)
    unprotect(1)
end

sexp(::Type{ExtPtrSxp}, s::Ptr{ExtPtrSxp}) = s
sexp(::Type{ExtPtrSxp}, r::RObject{ExtPtrSxp}) = sexp(r)
sexp(::Type{ClosSxp}, s::Ptr{ClosSxp}) = s
sexp(::Type{ClosSxp}, r::RObject{ClosSxp}) = sexp(r)

"""
Wrap a Julia object an a R `ExtPtrSxpPtr`.

We store the pointer and the object in a const Dict to prevent it being
removed by the Julia GC.
"""
function sexp(::Type{ExtPtrSxp}, j)
    jptr = pointer_from_objref(j)
    s = makeExternalPtr(jptr)
    jtypExtPtrs[s] = j
    registerCFinalizerEx(s)
    s
end

const juliaCallback = RObject{ExtPtrSxp}()

"""
Wrap a callable Julia object `f` an a R `ClosSxpPtr`.

Constructs the following R code

    function(...) .External(juliaCallback, fExPtr, ...)

"""
function sexp(::Type{ClosSxp}, f)
    fptr = protect(sexp(ExtPtrSxp,f))
    body = protect(rlang_p(Symbol(".External"),
                           juliaCallback,
                           fptr,
                           Const.DotsSymbol))
    local clos
    try
        lang = rlang_p(:function, sexp_arglist_dots(), body)
        clos = reval_p(lang)
    finally
        unprotect(2)
    end
    clos
end

"""
Create an argument list for an R function call, with a varargs "dots" at the end.
"""
function sexp_arglist_dots(args...;kwargs...)
    rarglist = protect(allocList(length(args)+length(kwargs)+1))
    try
        rr = rarglist
        for var in args
            settag!(rr, sexp(var))
            setcar!(rr, Const.MissingArg)
            rr = cdr(rr)
        end
        for (var,val) in kwargs
            settag!(rr, sexp(var))
            setcar!(rr, sexp(val))
            rr = cdr(rr)
        end
        settag!(rr, Const.DotsSymbol)
        setcar!(rr, Const.MissingArg)
    finally
        unprotect(1)
    end
    rarglist
end

sexp(f::Function) = sexp(ClosSxp, f)