text:
module Cowsay

using TextWrap

export cowsay
export cowthink
export cowsaid
export cowthunk

include("cows/blowfish.cow.jl")
include("cows/bunny.cow.jl")
include("cows/cower.cow.jl")
include("cows/default.cow.jl")
include("cows/dragon-and-cow.cow.jl")
include("cows/dragon.cow.jl")
include("cows/elephant-in-snake.cow.jl")
include("cows/elephant.cow.jl")
include("cows/eyes.cow.jl")
include("cows/flaming-sheep.cow.jl")
include("cows/fox.cow.jl")
include("cows/kitty.cow.jl")
include("cows/koala.cow.jl")
include("cows/mech-and-cow.cow.jl")
include("cows/meow.cow.jl")
include("cows/moofasa.cow.jl")
include("cows/moose.cow.jl")
include("cows/mutilated.cow.jl")
include("cows/sheep.cow.jl")
include("cows/skeleton.cow.jl")
include("cows/small.cow.jl")
include("cows/stegosaurus.cow.jl")
include("cows/supermilker.cow.jl")
include("cows/surgery.cow.jl")
include("cows/three-eyes.cow.jl")
include("cows/turkey.cow.jl")
include("cows/turtle.cow.jl")
include("cows/tux.cow.jl")
include("cows/udder.cow.jl")
include("cows/vader-koala.cow.jl")
include("cows/vader.cow.jl")
include("cows/www.cow.jl")

"""
    cowsay([io,] message::AbstractString; kwargs...)
    cowthink([io,] message::AbstractString; kwargs...)

Print an ASCII picture of a cow saying/thinking `message`

# Arguments
- `message::AbstractString`: Tell the cow what to say

# Keywords
- `cow=default`: Specify a particular function handle to print the ASCII art. See below for
    more details on what constitutes a valid function handle.
- `eyes::AbstractString="oo"`: A two-character string to be drawn in for the eyes. Not all
    cowfiles support this, though.
- `tongue::AbstractString="  "`: A two-character string to be drawn in for the tongue. Not
    all cowfiles support this.
- `wrap::Int=40`: The number of characters at which to wrap `message` to a new line
- `nowrap::Bool=false`: Don't perform text wrapping on `message`

## Cow appearance Keywords
- `borg::Bool=false`: Initiates Borg mode
- `dead::Bool=false`: Causes the cow to appear dead
- `greedy::Bool=false`: Invokes greedy mode
- `paranoid::Bool=false`: Causes a state of paranoia to come over the cow
- `stoned::Bool=false`: Makes the cow appear thoroughly stoned
- `tired::Bool=false`: Yields a tired cow
- `wired::Bool=false`: Somewhat the opposite of `tired`, and initiates wired mode
- `young::Bool=false`: Brings on the cow's youthful appearance

# Example
```jldoctest
julia> cowsay("Have you mooed today?")
 _______________________
< Have you mooed today? >
 -----------------------
        \\   ^__^
         \\  (oo)\\_______
            (__)\\       )\\/\\
                ||----w |
                ||     ||

julia> cowthink("Have I mooed today?")
 _____________________
( Have I mooed today? )
 ---------------------
        o   ^__^
         o  (oo)\\_______
            (__)\\       )\\/\\
                ||----w |
                ||     ||
```
"""
function cowsay(message::AbstractString; kwargs...)
    cowsay(stdout, message; kwargs...)
end

function cowsay(io::IO, message::AbstractString; kwargs...)
    println(io, cowsaid(message; kwargs...))
end

"""
    cowthink([io,] message::AbstractString; kwargs...)

Print an ASCII picture of a cow thinking `message`

See [`cowsay`](@ref)
"""
function cowthink(message::AbstractString; kwargs...)
    cowthink(stdout, message; kwargs...)
end

function cowthink(io::IO, message::AbstractString; kwargs...)
    println(io, cowthunk(message; kwargs...))
end

"""
    cowsaid(message::AbstractString; kwargs...)
    cowthunk(message::AbstractString; kwargs...)

Constructs cow art and returns it as a string instead of printing it. See [`cowsay`](@ref)
for implementation details.

# Example

```jldoctest
julia> cowsaid("Have you mooed today?")
" _______________________\\n< Have you mooed today? >\\n -----------------------\\n        \\\\   ^__^\\n         \\\\  (oo)\\\\_______\\n            (__)\\\\       )\\\\/\\\\\\n                ||----w |\\n                ||     ||\\n"
```

```jldoctest
julia> cowthunk("Have I mooed today?")
" _____________________\\n( Have I mooed today? )\\n ---------------------\\n        o   ^__^\\n         o  (oo)\\\\_______\\n            (__)\\\\       )\\\\/\\\\\\n                ||----w |\\n                ||     ||\\n"
```
"""
function cowsaid(message::AbstractString; kwargs...)
    return cowmoo(message, :say; kwargs...)
end

"""
    cowthunk(message::AbstractString; kwargs...)

Construct art of a cow thinking `message`.

See [`cowsaid`](@ref) and [`cowsay`](@ref)
"""
function cowthunk(message::AbstractString; kwargs...)
    return cowmoo(message, :think; kwargs...)
end

"""
    cowmoo(message::AbstractString, mode, kwargs...)

Converts `message` into either saying or thinking art based on the value of `mode`. All of
the default values for `kwargs` for pretty much every other function live here. See
[`cowsay`](@ref) for their descriptions.
"""
function cowmoo(message::AbstractString, mode; kwargs...)
    cow = dict_or_default(kwargs, :cow, Cowsay.default)
    eyes = dict_or_default(kwargs, :eyes, "oo")
    tongue = dict_or_default(kwargs, :tongue, "  ")
    wrap = dict_or_default(kwargs, :wrap, 40)
    nowrap = dict_or_default(kwargs, :nowrap, false)

    eyes, tongue = construct_face!(eyes, tongue; kwargs...)

    # Default to 'say' mode
    if mode == :think
        balloon = thinkballoon
        thoughts = "o"
    else
        balloon = sayballoon
        thoughts = "\\"
    end

    wrapped_message = !nowrap ? TextWrap.wrap(message, width=wrap, replace_whitespace=false) : message
    speechbubble = balloon(wrapped_message)

    return string(speechbubble, cow(eyes=eyes, tongue=tongue, thoughts=thoughts))
end

function construct_face!(eyes, tongue; kwargs...)
    borg = dict_or_default(kwargs, :borg, false)
    dead = dict_or_default(kwargs, :dead, false)
    greedy = dict_or_default(kwargs, :greedy, false)
    paranoid = dict_or_default(kwargs, :paranoid, false)
    stoned = dict_or_default(kwargs, :stoned, false)
    tired = dict_or_default(kwargs, :tired, false)
    wired = dict_or_default(kwargs, :wired, false)
    young = dict_or_default(kwargs, :young, false)

    if borg
        eyes = "=="
    end

    if dead
        eyes = "xx"
        tongue = "U "
    end

    if greedy
        eyes = "\$\$"
    end

    if paranoid
        eyes = "@@"
    end

    if stoned
        eyes = "**"
        tongue = "U "
    end

    if tired
        eyes = "--"
    end

    if wired
        eyes = "OO"
    end

    if young
        eyes = ".."
    end

    return eyes, tongue
end

"""
    dict_or_default(dict, key, value)

If `key` is present in `dict`, then return `dict[key]`, otherwise return `value`

# Example
```jldoctest
julia> dict = Dict(:holstein => "lotsa milk", :jersey => "butterfat", :shorthorn => "worthless");

julia> Cowsay.dict_or_default(dict, :holstein, "fluid")
"lotsa milk"

julia> Cowsay.dict_or_default(dict, :brown_swiss, "cheese")
"cheese"
```
"""
function dict_or_default(dict, key, value)
    if haskey(dict, key)
        return dict[key]
    else
        return value
    end
end

"""
    function sayballoon(message::AbstractString)

Creates the balloon portion of the art for [`cowsay`](@ref)

See also [`cowsay`](@ref), [`thinkballoon`](@ref)
"""
function sayballoon(message::AbstractString)
    messagelines = split(message, "\n")
    nlines = length(messagelines)

    linelength = max(length.(messagelines)...)
    paddinglength = linelength + 2

    if nlines == 1
        balloon = string(
            " ",
            repeat("_", paddinglength),
            "\n< ",
            message,
            " >\n ",
            repeat("-", paddinglength),
            "\n"
        )
    else
        balloon = string(
            " ",
            repeat("_", paddinglength),
            "\n/ ",
            rpad(messagelines[1], linelength),
            " \\\n"
        )
        if nlines > 2
            for i = 2:nlines-1
                balloon = string(
                    balloon,
                    "| ",
                    rpad(messagelines[i], linelength),
                    " |\n"
                )
            end
        end
        balloon = string(
            balloon,
            "\\ ",
            rpad(messagelines[end], linelength),
            " /\n ",
            repeat("-", paddinglength),
            "\n"
        )
    end

    return balloon

end #function

"""
    function thinkballoon(message::AbstractString)

Creates the balloon portion of the art for [`cowthink`](@ref)

See also [`cowthink`](@ref), [`sayballoon`](@ref)
"""
function thinkballoon(message::AbstractString)
    messagelines = split(message, "\n")
    nlines = length(messagelines)

    linelength = max(length.(messagelines)...)
    paddinglength = linelength + 2

    balloon = string(
        " ",
        repeat("_", paddinglength),
        "\n"
    )
    for i in 1:nlines
        balloon = string(
            balloon,
            "( ",
            rpad(messagelines[i], linelength),
            " )\n"
        )
    end #for
    balloon = string(
        balloon,
        " ",
        repeat("-", paddinglength),
        "\n"
    )

    return balloon
end #function

end


text:
macro def(genre::Symbol, ex::Union{Expr,Symbol}, vars::Symbol...)
    o = genre == :prop ? :(::) : genre == :vars ? :(=) : error("@def(genre = $genre, ...) is invalid.")
    n = length(vars)
    e = Vector{Expr}(undef, n)
    @inbounds for i in 1:n
        e[i] = Expr(o, vars[i], ex)
    end
    return Expr(:escape, Expr(:block, e...))
end

macro get(obj::Symbol, vars::Symbol...)
    n = length(vars)
    e = Vector{Expr}(undef, n)
    @inbounds for i in 1:n
        vari = vars[i]
        e[i] = :($vari = $obj.$vari)
    end
    return Expr(:escape, Expr(:block, e...))
end

macro copy!(ex::Expr...)
    n = length(ex)
    e = Vector{Expr}(undef, n)
    @inbounds for i in 1:n
        ei = ex[i]
        ai = ei.args
        if ai[1] == :(=>)
            ai[1] = :unsafe_copy!
            ai[2], ai[3] = ai[3], ai[2]
        else
            error("$e[$i] has invalid symbol.")
        end
        e[i] = ei
    end
    return Expr(:escape, Expr(:block, e...))
end

ext: jl


text:
# Test multihypo computation assembly

using Test

using IncrementalInference

@testset "test IncrementalInference.assembleHypothesesElements! without multihypothesis..." begin

s2_1_gt1 = 1:2
s2_1_gt2 = (1:20,Int[])
s2_1_gt3 = [(1,1:2); (2,Int[])]
s2_1_gt4 = Int[]
s2_1 = IncrementalInference.assembleHypothesesElements!(nothing, 20, 1, 2 )
@test sum([s2_1_gt1;] .- [s2_1[1];]) == 0
@test sum([s2_1_gt2[1];] .- [s2_1[2][1];]) == 0
@test sum(s2_1_gt2[2] .- s2_1[2][2]) == 0
@test s2_1_gt3[1][1] == s2_1[3][1][1]
@test sum([s2_1_gt3[1][2];] .- [s2_1[3][1][2];]) == 0
@test s2_1_gt3[2][1] == s2_1[3][2][1]
@test sum(s2_1_gt3[2][2] .- s2_1[3][2][2]) == 0
@test sum(s2_1_gt4 .- s2_1[4])  == 0

s2_2_gt1 = 1:2
s2_2_gt2 = (1:20,Int[])
s2_2_gt3 = [(1,1:2); (2,Int[])]
s2_2_gt4 = Int[]

s2_2 = IncrementalInference.assembleHypothesesElements!(nothing, 20, 2, 2 )

@test sum([s2_2_gt1;] .- [s2_2[1];]) == 0
@test sum([s2_2_gt2[1];] .- [s2_2[2][1];]) == 0
@test sum(s2_2_gt2[2] .- s2_2[2][2]) == 0
@test s2_2_gt3[1][1] == s2_2[3][1][1]
@test sum([s2_2_gt3[1][2];] .- [s2_2[3][1][2];]) == 0
@test s2_2_gt3[2][1] == s2_2[3][2][1]
@test sum(s2_2_gt3[2][2] .- s2_2[3][2][2]) == 0
@test sum(s2_2_gt4 .- s2_2[4])  == 0

end

@testset "test IncrementalInference.assembleHypothesesElements! with bi-modality..." begin

# certainidx = 1
# sfidx=1, mhidx=1:  ah = []
# sfidx=1, mhidx=2:  ah = [1;2]
# sfidx=1, mhidx=3:  ah = [1;3]
s3_1_gt1 = [1]
s3_1_gt2 = (0,4,4,40)
s3_1_gt3 = [(1,Int[]); (2,Int[1;2]); (3,Int[1;3])]
s3_1_gt4 = 40

s3_1 = IncrementalInference.assembleHypothesesElements!(Categorical([0.0;0.5;0.5]), 40, 1, 3)

@test sum([s3_1_gt1;] - [s3_1[1];]) == 0
@test sum([s3_1_gt2[1];] .- [s3_1[2][1];]) == 0
@test length(s3_1[2][2]) > s3_1_gt2[2]
@test length(s3_1[2][3]) > s3_1_gt2[3]
@test length(s3_1[2][2]) + length(s3_1[2][3]) == s3_1_gt2[4]
@test s3_1_gt3[1][1] == s3_1[3][1][1]
@test s3_1_gt3[2][1] == s3_1[3][2][1]
@test s3_1_gt3[3][1] == s3_1[3][3][1]
@test sum(s3_1_gt3[1][2] .- s3_1[3][1][2]) == 0
@test sum(s3_1_gt3[2][2] .- s3_1[3][2][2]) == 0
@test sum(s3_1_gt3[3][2] .- s3_1[3][3][2]) == 0

@test sum(s3_1[4] .== 2) > s3_1_gt2[2]
@test sum(s3_1[4] .== 3) > s3_1_gt2[3]

@test sum( [1:40;][s3_1[4] .== 2] .== s3_1[2][2] ) == length(s3_1[2][2])
@test sum( [1:40;][s3_1[4] .== 3] .== s3_1[2][3] ) == length(s3_1[2][3])
@test length(s3_1[4]) == s3_1_gt4

# certainidx = 1
# sfidx=2, mhidx=1:  ah = [1;2]
# sfidx=2, mhidx=2:  ah = [2;3]
# sfidx=2, mhidx=3:  [2;3], 2 should take a value from 3
s3_2_gt1 = [1]
s3_2_gt2 = (0,4,4,40)
s3_2_gt3 = [(1,Int[1;2]); (2,Int[1;2]); (3,Int[2;3])]
s3_2_gt4 = 40

s3_2 = IncrementalInference.assembleHypothesesElements!(Categorical([0.0;0.5;0.5]), 40, 2, 3 )

@test sum(s3_2_gt1 - s3_2[1]) == 0
@test sum(s3_2_gt2[1] .- s3_2[2][1]) == 0
@test length(s3_2[2][2]) > s3_2_gt2[2]
@test length(s3_2[2][3]) > s3_2_gt2[3]
@test length(s3_2[2][2]) + length(s3_2[2][3]) == s3_2_gt2[4]
@test s3_2_gt3[1][1] == s3_2[3][1][1]
@test s3_2_gt3[2][1] == s3_2[3][2][1]
@test s3_2_gt3[3][1] == s3_2[3][3][1]
@test sum(s3_2_gt3[1][2] .- s3_2[3][1][2]) == 0
@test sum(s3_2_gt3[2][2] .- s3_2[3][2][2]) == 0
@test sum(s3_2_gt3[3][2] .- s3_2[3][3][2]) == 0

@test sum(s3_2[4] .== 2) > s3_2_gt2[2]
@test sum(s3_2[4] .== 3) > s3_2_gt2[3]

@test sum( [1:40;][s3_2[4] .== 2] .== s3_2[2][2] ) == length(s3_2[2][2])
@test sum( [1:40;][s3_2[4] .== 3] .== s3_2[2][3] ) == length(s3_2[2][3])
@test length(s3_2[4]) == s3_2_gt4

# certainidx = 1
# sfidx=3, mhidx=1:  ah = [1;3]
# sfidx=3, mhidx=2:  [2:3], 3 should take a value from 2
# sfidx=3, mhidx=3:  ah = [1;3]
s3_3_gt1 = [1]
s3_3_gt2 = (0,4,4,40)
s3_3_gt3 = [(1,Int[1;3]); (2,Int[2;3]); (3,Int[1;3])]
s3_3_gt4 = 40

s3_3 = IncrementalInference.assembleHypothesesElements!(Categorical([0.0;0.5;0.5]), 40, 3, 3 )

@test sum(s3_3_gt1 - s3_3[1]) == 0
@test sum(s3_3_gt2[1] .- s3_3[2][1]) == 0
@test length(s3_3[2][2]) > s3_3_gt2[2]
@test length(s3_3[2][3]) > s3_3_gt2[3]
@test length(s3_3[2][2]) + length(s3_3[2][3]) == s3_3_gt2[4]
@test s3_3_gt3[1][1] == s3_3[3][1][1]
@test s3_3_gt3[2][1] == s3_3[3][2][1]
@test s3_3_gt3[3][1] == s3_3[3][3][1]
@test sum(s3_3_gt3[1][2] .- s3_3[3][1][2]) == 0
@test sum(s3_3_gt3[2][2] .- s3_3[3][2][2]) == 0
@test sum(s3_3_gt3[3][2] .- s3_3[3][3][2]) == 0

@test sum(s3_3[4] .== 2) > s3_3_gt2[2]
@test sum(s3_3[4] .== 3) > s3_3_gt2[3]

@test sum( [1:40;][s3_3[4] .== 2] .== s3_3[2][2] ) == length(s3_3[2][2])
@test sum( [1:40;][s3_3[4] .== 3] .== s3_3[2][3] ) == length(s3_3[2][3])
@test length(s3_3[4]) == s3_3_gt4

end

# @testset "test IncrementalInference.assembleHypothesesElements! with bi-modality backwards permutation..." begin

# certainidx = 1
# sfidx=1, mhidx=1:  ah = []
# sfidx=1, mhidx=2:  ah = [1;2]
# sfidx=1, mhidx=3:  ah = [1;3]
# s3_1_gt1 = [1]
# s3_1_gt2 = (0,4,4,20)
# s3_1_gt3 = [(1,Int[]); (2,Int[1;2]); (3,Int[1;3])]
# s3_1_gt4 = 20
#
# s3_1 = IncrementalInference.assembleHypothesesElements!(Categorical([0.5;0.5;0.0]), 20, 1, 3)

# @test sum(s3_1_gt1 - s3_1[1]) == 0
# @test sum(s3_1_gt2[1] .- s3_1[2][1]) == 0
# @test length(s3_1[2][2]) > s3_1_gt2[2]
# @test length(s3_1[2][3]) > s3_1_gt2[3]
# @test length(s3_1[2][2]) + length(s3_1[2][3]) == s3_1_gt2[4]
# @test s3_1_gt3[1][1] == s3_1[3][1][1]
# @test s3_1_gt3[2][1] == s3_1[3][2][1]
# @test s3_1_gt3[3][1] == s3_1[3][3][1]
# @test sum(s3_1_gt3[1][2] .- s3_1[3][1][2]) == 0
# @test sum(s3_1_gt3[2][2] .- s3_1[3][2][2]) == 0
# @test sum(s3_1_gt3[3][2] .- s3_1[3][3][2]) == 0
#
# @test sum(s3_1[4] .== 2) > s3_1_gt2[2]
# @test sum(s3_1[4] .== 3) > s3_1_gt2[3]
#
# @test sum( [1:20;][s3_1[4] .== 2] .== s3_1[2][2] ) == length(s3_1[2][2])
# @test sum( [1:20;][s3_1[4] .== 3] .== s3_1[2][3] ) == length(s3_1[2][3])
# @test length(s3_1[4]) == s3_1_gt4
#

# end

@testset "test IncrementalInference.assembleHypothesesElements! with tri-modality..." begin

N = 50
s4_1_gt1 = [1]
s4_1_gt2 = (0,4,4,4,N)
s4_1_gt3 = [(1,Int[]); (2,Int[1;2]); (3,Int[1;3]); (4,Int[1;4])]
s4_1_gt4 = N

s4_1 = IncrementalInference.assembleHypothesesElements!(Categorical([0.0;0.33;0.33;0.34]), N, 1, 4 )

@test sum(s4_1_gt1 - s4_1[1]) == 0
@test sum(s4_1_gt2[1] .- s4_1[2][1]) == 0
@test length(s4_1[2][2]) > s4_1_gt2[2]
@test length(s4_1[2][3]) > s4_1_gt2[3]
@test length(s4_1[2][4]) > s4_1_gt2[4]
@test length(s4_1[2][2]) + length(s4_1[2][3]) + length(s4_1[2][4]) == s4_1_gt2[5]

@test s4_1_gt3[1][1] == s4_1[3][1][1]
@test s4_1_gt3[2][1] == s4_1[3][2][1]
@test s4_1_gt3[3][1] == s4_1[3][3][1]
@test s4_1_gt3[4][1] == s4_1[3][4][1]
@test sum(s4_1_gt3[1][2] .- s4_1[3][1][2]) == 0
@test sum(s4_1_gt3[2][2] .- s4_1[3][2][2]) == 0
@test sum(s4_1_gt3[3][2] .- s4_1[3][3][2]) == 0
@test sum(s4_1_gt3[4][2] .- s4_1[3][4][2]) == 0

@test sum(s4_1[4] .== 2) > s4_1_gt2[2]
@test sum(s4_1[4] .== 3) > s4_1_gt2[3]
@test sum(s4_1[4] .== 4) > s4_1_gt2[4]

@test sum( [1:N;][s4_1[4] .== 2] .== s4_1[2][2] ) == length(s4_1[2][2])
@test sum( [1:N;][s4_1[4] .== 3] .== s4_1[2][3] ) == length(s4_1[2][3])
@test sum( [1:N;][s4_1[4] .== 4] .== s4_1[2][4] ) == length(s4_1[2][4])
@test length(s4_1[4]) == s4_1_gt4

N = 50
s4_2_gt1 = [1]
s4_2_gt2 = (0,4,4,4,N)
s4_2_gt3 = [(1,Int[1;2]); (2,Int[1;2]); (3,Int[2;3;4]); (4,Int[2;3;4])]
s4_2_gt4 = N

s4_2 = IncrementalInference.assembleHypothesesElements!(Categorical([0.0;0.33;0.33;0.34]), N, 2, 4 )

@test sum(s4_2_gt1 - s4_2[1]) == 0
@test sum(s4_2_gt2[1] .- s4_2[2][1]) == 0
@test length(s4_2[2][2]) > s4_2_gt2[2]
@test length(s4_2[2][3]) > s4_2_gt2[3]
@test length(s4_2[2][4]) > s4_2_gt2[4]
@test length(s4_2[2][2]) + length(s4_2[2][3]) + length(s4_2[2][4]) == s4_2_gt2[5]

@test s4_2_gt3[1][1] == s4_2[3][1][1]
@test s4_2_gt3[2][1] == s4_2[3][2][1]
@test s4_2_gt3[3][1] == s4_2[3][3][1]
@test s4_2_gt3[4][1] == s4_2[3][4][1]
@test sum(s4_2_gt3[1][2] .- s4_2[3][1][2]) == 0
@test sum(s4_2_gt3[2][2] .- s4_2[3][2][2]) == 0
@test sum(s4_2_gt3[3][2] .- s4_2[3][3][2]) == 0
@test sum(s4_2_gt3[4][2] .- s4_2[3][4][2]) == 0

@test sum(s4_2[4] .== 2) > s4_2_gt2[2]
@test sum(s4_2[4] .== 3) > s4_2_gt2[3]
@test sum(s4_2[4] .== 4) > s4_2_gt2[4]

@test sum( [1:N;][s4_2[4] .== 2] .== s4_2[2][2] ) == length(s4_2[2][2])
@test sum( [1:N;][s4_2[4] .== 3] .== s4_2[2][3] ) == length(s4_2[2][3])
@test sum( [1:N;][s4_2[4] .== 4] .== s4_2[2][4] ) == length(s4_2[2][4])
@test length(s4_2[4]) == s4_2_gt4

N = 50
s4_3_gt1 = [1]
s4_3_gt2 = (0,4,4,4,N)
s4_3_gt3 = [(1,Int[1;3]); (2,Int[2;3;4]); (3,Int[1;3]); (4,Int[2;3;4])]
s4_3_gt4 = N

s4_3 = IncrementalInference.assembleHypothesesElements!(Categorical([0.0;0.33;0.33;0.34]), N, 3, 4 )

@test sum(s4_3_gt1 - s4_3[1]) == 0
@test sum(s4_3_gt2[1] .- s4_3[2][1]) == 0
@test length(s4_3[2][2]) > s4_3_gt2[2]
@test length(s4_3[2][3]) > s4_3_gt2[3]
@test length(s4_3[2][4]) > s4_3_gt2[4]
@test length(s4_3[2][2]) + length(s4_3[2][3]) + length(s4_3[2][4]) == s4_3_gt2[5]

@test s4_3_gt3[1][1] == s4_3[3][1][1]
@test s4_3_gt3[2][1] == s4_3[3][2][1]
@test s4_3_gt3[3][1] == s4_3[3][3][1]
@test s4_3_gt3[4][1] == s4_3[3][4][1]
@test sum(s4_3_gt3[1][2] .- s4_3[3][1][2]) == 0
@test sum(s4_3_gt3[2][2] .- s4_3[3][2][2]) == 0
@test sum(s4_3_gt3[3][2] .- s4_3[3][3][2]) == 0
@test sum(s4_3_gt3[4][2] .- s4_3[3][4][2]) == 0

@test sum(s4_3[4] .== 2) > s4_3_gt2[2]
@test sum(s4_3[4] .== 3) > s4_3_gt2[3]
@test sum(s4_3[4] .== 4) > s4_3_gt2[4]

@test sum( [1:N;][s4_3[4] .== 2] .== s4_3[2][2] ) == length(s4_3[2][2])
@test sum( [1:N;][s4_3[4] .== 3] .== s4_3[2][3] ) == length(s4_3[2][3])
@test sum( [1:N;][s4_3[4] .== 4] .== s4_3[2][4] ) == length(s4_3[2][4])
@test length(s4_3[4]) == s4_3_gt4

N = 50
s4_4_gt1 = [1]
s4_4_gt2 = (0,4,4,4,N)
s4_4_gt3 = [(1,Int[1;4]); (2,Int[2;3;4]); (3,Int[2;3;4]); (4,Int[1;4])]
s4_4_gt4 = N

s4_4 = IncrementalInference.assembleHypothesesElements!(Categorical([0.0;0.33;0.33;0.34]), N, 4, 4 )

@test sum(s4_4_gt1 - s4_4[1]) == 0
@test sum(s4_4_gt2[1] .- s4_4[2][1]) == 0
@test length(s4_4[2][2]) > s4_4_gt2[2]
@test length(s4_4[2][3]) > s4_4_gt2[3]
@test length(s4_4[2][4]) > s4_4_gt2[4]
@test length(s4_4[2][2]) + length(s4_4[2][3]) + length(s4_4[2][4]) == s4_4_gt2[5]

@test s4_4_gt3[1][1] == s4_4[3][1][1]
@test s4_4_gt3[2][1] == s4_4[3][2][1]
@test s4_4_gt3[3][1] == s4_4[3][3][1]
@test s4_4_gt3[4][1] == s4_4[3][4][1]
@test sum(s4_4_gt3[1][2] .- s4_4[3][1][2]) == 0
@test sum(s4_4_gt3[2][2] .- s4_4[3][2][2]) == 0
@test sum(s4_4_gt3[3][2] .- s4_4[3][3][2]) == 0
@test sum(s4_4_gt3[4][2] .- s4_4[3][4][2]) == 0

@test sum(s4_4[4] .== 2) > s4_4_gt2[2]
@test sum(s4_4[4] .== 3) > s4_4_gt2[3]
@test sum(s4_4[4] .== 4) > s4_4_gt2[4]

@test sum( [1:N;][s4_4[4] .== 2] .== s4_4[2][2] ) == length(s4_4[2][2])
@test sum( [1:N;][s4_4[4] .== 3] .== s4_4[2][3] ) == length(s4_4[2][3])
@test sum( [1:N;][s4_4[4] .== 4] .== s4_4[2][4] ) == length(s4_4[2][4])
@test length(s4_4[4]) == s4_4_gt4

@warn "only partially testing tri-modality"

end


text:
"""
    barplot(x, y; kwargs...)

Plots a barplot; `y` defines the height.  `x` and `y` should be 1 dimensional.

## Attributes
$(ATTRIBUTES)
"""
@recipe(BarPlot, x, y) do scene
    Attributes(;
        fillto = automatic,
        color = theme(scene, :patchcolor),
        colormap = theme(scene, :colormap),
        colorrange = automatic,
        dodge = automatic,
        n_dodge = automatic,
        x_gap = 0.2,
        dodge_gap = 0.03,
        marker = Rect,
        stack = automatic,
        strokewidth = theme(scene, :patchstrokewidth),
        strokecolor = theme(scene, :patchstrokecolor),
        width = automatic,
        direction = :y,
        visible = theme(scene, :visible),
        inspectable = theme(scene, :inspectable),
        cycle = [:color => :patchcolor],
    )
end

conversion_trait(::Type{<: BarPlot}) = PointBased()

function bar_rectangle(x, y, width, fillto)
    # y could be smaller than fillto...
    ymin = min(fillto, y)
    ymax = max(fillto, y)
    w = abs(width)
    return FRect(x - (w / 2f0), ymin, w, ymax - ymin)
end

flip(r::Rect2D) = Rect2D(reverse(origin(r)), reverse(widths(r)))

function xw_from_dodge(x, width, minimum_distance, x_gap, dodge, n_dodge, dodge_gap)
    width === automatic && (width = (1 - x_gap) * minimum_distance)
    if dodge === automatic
        i_dodge = 1
    elseif eltype(dodge) <: Integer
        i_dodge = dodge
    else
        ArgumentError("The keyword argument `dodge` currently supports only `AbstractVector{<: Integer}`") |> throw
    end
    n_dodge === automatic && (n_dodge = maximum(i_dodge))
    dodge_width = scale_width(dodge_gap, n_dodge)
    shifts = shift_dodge.(i_dodge, dodge_width, dodge_gap)
    return x .+ width .* shifts, width * dodge_width
end

function AbstractPlotting.plot!(p::BarPlot)

    in_y_direction = lift(p.direction) do dir
        if dir == :y
            true
        elseif dir == :x
            false
        else
            error("Invalid direction $dir. Options are :x and :y.")
        end
    end

    bars = lift(p[1], p.fillto, p.width, p.dodge, p.n_dodge, p.x_gap, p.dodge_gap, p.stack, in_y_direction) do xy, fillto, width, dodge, n_dodge, x_gap, dodge_gap, stack, in_y_direction

        x = first.(xy)
        y = last.(xy)

        minimum_distance = nothing
        # only really compute `minimum_distance` if `width` is `automatic`
        if width === automatic
            x_unique = unique(filter(isfinite, x))
            x_diffs = diff(sort(x_unique))
            minimum_distance = isempty(x_diffs) ? 1.0 : minimum(x_diffs)
        end

        # compute width of bars and x̂ (horizontal position after dodging)
        x̂, barwidth = xw_from_dodge(x, width, minimum_distance, x_gap, dodge, n_dodge, dodge_gap)

        # --------------------------------
        # ----------- Stacking -----------
        # --------------------------------

        if stack === automatic
            if fillto === automatic
                fillto = 0.0
            end
        elseif eltype(stack) <: Integer
            fillto === automatic || @warn "Ignore keyword fillto when keyword stack is provided"
            i_stack = stack

            from, to = stack_grouped_from_to(i_stack, y, (x = x̂,))
            y, fillto = to, from
        else
            ArgumentError("The keyword argument `stack` currently supports only `AbstractVector{<: Integer}`") |> throw
        end

        rects = @. bar_rectangle(x̂, y, barwidth, fillto)
        return in_y_direction ? rects : flip.(rects)
    end

    poly!(
        p, bars, color = p.color, colormap = p.colormap, colorrange = p.colorrange,
        strokewidth = p.strokewidth, strokecolor = p.strokecolor, visible = p.visible,
        inspectable = p.inspectable
    )
end

scale_width(dodge_gap, n_dodge) = (1 - (n_dodge - 1) * dodge_gap) / n_dodge

function shift_dodge(i, dodge_width, dodge_gap)
    (dodge_width - 1) / 2 + (i - 1) * (dodge_width + dodge_gap)
end

function stack_grouped_from_to(i_stack, y, grp)

	from = Array{Float64}(undef, length(y))
	to   = Array{Float64}(undef, length(y))

	groupby = StructArray((; grp..., is_pos = y .> 0))

	grps = StructArrays.finduniquesorted(groupby)

	for (grp, inds) in grps

		fromto = stack_from_to(i_stack[inds], y[inds])

		from[inds] .= fromto.from
		to[inds] .= fromto.to

	end

	(from = from, to = to)
end

function stack_from_to(i_stack, y)
	# save current order
	order = 1:length(y)
	# sort by i_stack
	perm = sortperm(i_stack)
	# restore original order
	inv_perm = sortperm(order[perm])

	from, to = stack_from_to_sorted(view(y, perm))

	(from = view(from, inv_perm), to = view(to, inv_perm))
end

function stack_from_to_sorted(y)
	to = cumsum(y)
	from = [0.0; to[firstindex(to):end-1]]

	(from = from, to = to)
end


text:
module PriorityChannels
export PriorityChannel
using DataStructures
import Base: notify_error, register_taskdone_hook, check_channel_state

const PriorityElement{T,I<:Real} = Tuple{T,I}
const ordering = Base.By(x->x[2])

"""
    PriorityChannel{T}(sz::Int)

Constructs a `PriorityChannel` with an internal buffer that can hold a maximum of `sz` objects
of type `T`, each assigned an real-valued priority (low number = higher priority).
[`put!`](@ref) calls on a full channel block until an object is removed with [`take!`](@ref).

`PriorityChannel(0)` constructs an unbuffered channel. `put!` blocks until a matching `take!` is called.
And vice-versa.

Other constructors:

* `PriorityChannel(Inf)`: equivalent to `PriorityChannel{Any,Int}(typemax(Int))`
* `PriorityChannel(sz)`: equivalent to `PriorityChannel{Any,Int}(sz)`
"""
mutable struct PriorityChannel{T,I} <: AbstractChannel{T}
    cond_take::Condition                 # waiting for data to become available
    cond_put::Condition                  # waiting for a writeable slot
    state::Symbol
    excp::Union{Exception, Nothing}         # exception to be thrown when state != :open

    data::Vector{PriorityElement{T,I}}
    sz_max::Int                          # maximum size of channel

    # Used when sz_max == 0, i.e., an unbuffered channel.
    waiters::Int
    takers::Vector{Task}
    putters::Vector{Task}

    function PriorityChannel{T,I}(sz::Float64) where {T,I<:Real}
        if sz == Inf
            PriorityChannel{T,I}(typemax(Int))
        else
            PriorityChannel{T,I}(convert(Int, sz))
        end
    end
    function PriorityChannel{T,I}(sz::Integer) where {T,I<:Real}
        if sz <= 0
            throw(ArgumentError("Channel size must be a positive integer or Inf"))
        end
        ch = new(Condition(), Condition(), :open, nothing, Vector{PriorityElement{T,I}}(), sz, 0)
        return ch
    end
end

PriorityChannel(sz) = PriorityChannel{Any,Int}(sz)

# special constructors
"""
    PriorityChannel(func::Function; ctype=Any, csize=1, taskref=nothing)

Create a new task from `func`, bind it to a new channel of type
`ctype` and size `csize`, and schedule the task, all in a single call.

`func` must accept the bound channel as its only argument.

If you need a reference to the created task, pass a `Ref{Task}` object via
keyword argument `taskref`.

Return a `PriorityChannel`.

# Examples
```jldoctest
julia> chnl = PriorityChannel(c->foreach(i->put!(c,i), 1:4));

julia> typeof(chnl)
PriorityChannel{Any,Int64}

julia> for i in chnl
           @show i
       end;
i = 1
i = 2
i = 3
i = 4
```

Referencing the created task:

```jldoctest
julia> taskref = Ref{Task}();

julia> chnl = PriorityChannel(c->(@show take!(c)); taskref=taskref);

julia> istaskdone(taskref[])
false

julia> put!(chnl, "Hello");
take!(c) = "Hello"

julia> istaskdone(taskref[])
true
```
"""
function PriorityChannel(func::Function; ctype=Any, csize=1, taskref=nothing)
    csize > 0 || throw(ArgumentError("Channel size must be a positive integer or Inf"))
    chnl = PriorityChannel{ctype,Int}(csize)
    task = Task(() -> func(chnl))
    bind(chnl, task)
    yield(task) # immediately start it

    isa(taskref, Ref{Task}) && (taskref[] = task)
    return chnl
end

closed_exception() = InvalidStateException("Channel is closed.", :closed)

isbuffered(c::PriorityChannel) = true

function Base.check_channel_state(c::PriorityChannel)
    if !isopen(c)
        c.excp !== nothing && throw(c.excp)
        throw(closed_exception())
    end
end

function Base.close(c::PriorityChannel)
    c.state = :closed
    c.excp = closed_exception()
    notify_error(c)
    nothing
end
Base.isopen(c::PriorityChannel) = (c.state == :open)

function Base.bind(c::PriorityChannel, task::Task)
    ref = WeakRef(c)
    register_taskdone_hook(task, tsk->Base.close_chnl_on_taskdone(tsk, ref))
    c
end

"""
    channeled_tasks(n::Int, funcs...; ctypes=fill(Any,n), csizes=fill(0,n))

A convenience method to create `n` channels and bind them to tasks started
from the provided functions in a single call. Each `func` must accept `n` arguments
which are the created channels. PriorityChannel types and sizes may be specified via
keyword arguments `ctypes` and `csizes` respectively. If unspecified, all channels are
of type `Channel{Any}(0)`.

Returns a tuple, `(Array{Channel}, Array{Task})`, of the created channels and tasks.
"""
function channeled_tasks(n::Int, funcs...; ctypes=fill(Any,n), csizes=fill(0,n))
    @assert length(csizes) == n
    @assert length(ctypes) == n

    chnls = map(i -> PriorityChannel{ctypes[i],Int}(csizes[i]), 1:n)
    tasks = Task[ Task(() -> f(chnls...)) for f in funcs ]

    # bind all tasks to all channels and schedule them
    foreach(t -> foreach(c -> bind(c, t), chnls), tasks)
    foreach(schedule, tasks)
    yield() # Allow scheduled tasks to run

    return (chnls, tasks)
end

"""
    put!(c::PriorityChannel, v, p)

Append an item `v` to the channel `c` with priority `p`. Blocks if the channel is full.
"""
function Base.put!(c::PriorityChannel{T,I}, v,i::I = 0) where {T,I<:Real}
    check_channel_state(c)
    while length(c.data) == c.sz_max
        wait(c.cond_put)
    end
    heappush!(c.data, PriorityElement((v,i)), ordering)

    # notify all, since some of the waiters may be on a "fetch" call.
    notify(c.cond_take, nothing, true, false)
    v
end

Base.push!(c::PriorityChannel, v, i=0) = put!(c, v, i)

"""
    fetch(c::PriorityChannel)

Wait for and get the highest priority item from the channel. Does not
remove the item.
"""
function Base.fetch(c::PriorityChannel)
    wait(c)
    c.data[1][1]
end

"""
    take!(c::PriorityChannel)

Remove and return the highest priority value from a [`PriorityChannel`](@ref). Blocks until data is available.
"""
function Base.take!(c::PriorityChannel)
    wait(c)
    v = heappop!(c.data, ordering)[1]
    notify(c.cond_put, nothing, false, false) # notify only one, since only one slot has become available for a put!.
    v
end

Base.popfirst!(c::PriorityChannel) = take!(c)

Base.isready(c::PriorityChannel) = n_avail(c) > 0
n_avail(c::PriorityChannel) = length(c.data)

function Base.wait(c::PriorityChannel)
    while !isready(c)
        check_channel_state(c)
        wait(c.cond_take)
    end
    nothing
end

function Base.notify_error(c::PriorityChannel, err)
    Base.notify_error(c.cond_take, err)
    Base.notify_error(c.cond_put, err)
end
Base.notify_error(c::PriorityChannel) = Base.notify_error(c, c.excp)

Base.eltype(::Type{PriorityChannel{T,I}}) where {T,I} = T

Base.show(io::IO, c::PriorityChannel) = print(io, "$(typeof(c))(sz_max:$(c.sz_max),sz_curr:$(n_avail(c)))")

function Base.iterate(c::PriorityChannel, state=nothing)
    try
        return (take!(c), nothing)
    catch e
        if isa(e, InvalidStateException) && e.state==:closed
            return nothing
        else
            rethrow()
        end
    end
end

Base.IteratorSize(::Type{<:PriorityChannel}) = SizeUnknown()

end # module


text:
const mc_sampling = :mc_sampling
const gauss_hermite = :gauss_hermite
const gauss_legendre = :gauss_legendre
const gauss_laguerre = :gauss_laguerre
const trapezoid = :trapezoid

"""
    generate_measure_data(params::Union{InfiniteOpt.ParameterRef,
                          AbstractArray{<:InfiniteOpt.ParameterRef}},
                          num_supports::Int,
                          lb::Union{Number, JuMPC.SparseAxisArray, Nothing} = nothing,
                          ub::Union{Number, JuMPC.SparseAxisArray, Nothing} = nothing;
                          eval_method::Symbol = mc_sampling,
                          name::String = "measure",
                          weight_func::Function = InfiniteOpt._w,
                          kwargs...)::InfiniteOpt.AbstractMeasureData

Generate an [`AbstractMeasureData`](@ref) object that automatically generate
supports based on a set of given information. The information is required to
include parameters and number of supports. Other optional information to input
includes lower and upper bounds, measure name, support generation method, and
weight functions. The users could supply extra keyword arguments if necessary
for their custom support generation methods.

**Example**
```jldoctest; setup = :(using InfiniteOpt, JuMP; model = InfiniteModel())
julia> @infinite_parameter(m, x in [0., 1.]);

julia> measure_data = generate_measure_data(x, 3, 0.3, 0.7, method = gauss_legendre)
DiscreteMeasureData(x, [0.1111111111111111, 0.17777777777777776, 0.1111111111111111], [0.3450806661517033, 0.5, 0.6549193338482967], "", InfiniteOpt._w)
```
"""
function generate_measure_data(params::Union{InfiniteOpt.ParameterRef,
                               AbstractArray{<:InfiniteOpt.ParameterRef}},
                               num_supports::Int,
                               lb::Union{Number, JuMPC.SparseAxisArray, Nothing} = nothing,
                               ub::Union{Number, JuMPC.SparseAxisArray, Nothing} = nothing;
                               eval_method::Symbol = mc_sampling,
                               name::String = "measure",
                               weight_func::Function = InfiniteOpt._w,
                               kwargs...)::InfiniteOpt.AbstractMeasureData
    if isa(params, InfiniteOpt.ParameterRef)
        set = InfiniteOpt._parameter_set(params)
    else
        params = convert(JuMPC.SparseAxisArray, params)
        set = InfiniteOpt._parameter_set(first(params))
        sets = InfiniteOpt._parameter_set.(params)
        if any(typeof(set) != typeof(s) for s in sets)
            error("Automatic measure data generation for multi-dimensional " *
                  "parameters with mixed set types is not supported.")
        end
    end
    (supports, coeffs) = generate_supports_and_coeffs(set, params, num_supports, lb, ub, Val(eval_method); kwargs...)
    return InfiniteOpt.DiscreteMeasureData(params, coeffs, supports,
                                           name = name, weight_function = weight_func)
end

"""
    generate_supports_and_coeffs(set::InfiniteOpt.AbstractInfiniteSet,
                                 params::Union{InfiniteOpt.ParameterRef,
                                 AbstractArray{<:InfiniteOpt.ParameterRef}},
                                 num_supports::Int,
                                 lb::Union{Number, JuMPC.SparseAxisArray, Nothing},
                                 ub::Union{Number, JuMPC.SparseAxisArray, Nothing},
                                 method::Val{eval_method}; [kwargs...])::Tuple

Generate supports and coefficients using the arguments and
`kwargs` as appropriate. This will dispatch to `eval_method` in accordance with the
type of `set`. This is intended as an internal method for
[`generate_measure_data`](@ref InfiniteOpt.MeasureEvalMethods.generate_measure_data)
and will need to be extended for user-defined
infinite set types. The output tuple should consists of two data objects, the
first being the supports and the second being the coefficients. Supports must
be a `Vector` of `Number` or of `AbstractArray{<:Number}`, and coefficients
must be a `Vector` of `Number`. Extensions will be needed to implement their
own version of function `generate_supports_and_coeffs(set::AbstractInfiniteSet, params, num_supports, lb, ub, method::Val{my_method};kwargs...)`.
Refer to the other dispatches of `generate_supports_and_coeffs` for details of
specific evaluation method implementations.
"""
function generate_supports_and_coeffs end

# fallback
function generate_supports_and_coeffs(set::InfiniteOpt.AbstractInfiniteSet,
                                      params::Union{InfiniteOpt.ParameterRef,
                                      AbstractArray{<:InfiniteOpt.ParameterRef}},
                                      num_supports::Int,
                                      lb::Union{Number, JuMPC.SparseAxisArray, Nothing},
                                      ub::Union{Number, JuMPC.SparseAxisArray, Nothing},
                                      method; kwargs...)::Tuple
    error("`generate_supports_and_coeffs` is not extended for parameters in sets " *
          "of type $(typeof(set)) with method " * string(method)[6:length(string(method))-3] * ".")
end

# MC Sampling for IntervalSet
"""
    generate_supports_and_coeffs(set::InfiniteOpt.IntervalSet,
                                 params::Union{InfiniteOpt.ParameterRef,
                                 AbstractArray{<:InfiniteOpt.ParameterRef}},
                                 num_supports::Int,
                                 lb::Union{Number, JuMPC.SparseAxisArray, Nothing},
                                 ub::Union{Number, JuMPC.SparseAxisArray, Nothing},
                                 method::Val{mc_sampling})::Tuple

Return a tuple that contains supports and coefficients generated by Monte Carlo
sampling from a uniform distribution between the lower and upper bounds provided.

**Example**
```jldoctest; setup = :(using InfiniteOpt, JuMP; model = InfiniteModel(seed = true))
julia> @infinite_parameter(model, t in [0., 1.]);

julia> (supps, coeffs) = generate_supports_and_coeffs(IntervalSet(0,1), t, 5, 0., 1., Val(mc_sampling))
([0.8236475079774124, 0.9103565379264364, 0.16456579813368521, 0.17732884646626457, 0.278880109331201], [0.2, 0.2, 0.2, 0.2, 0.2])

julia> supps
5-element Array{Float64,1}:
 0.8236475079774124
 0.9103565379264364
 0.16456579813368521
 0.17732884646626457
 0.278880109331201
```
"""
function generate_supports_and_coeffs(set::InfiniteOpt.IntervalSet,
                                      params::Union{InfiniteOpt.ParameterRef,
                                      AbstractArray{<:InfiniteOpt.ParameterRef}},
                                      num_supports::Int,
                                      lb::Union{Number, JuMPC.SparseAxisArray, Nothing},
                                      ub::Union{Number, JuMPC.SparseAxisArray, Nothing},
                                      method::Val{mc_sampling})::Tuple
    _nothing_test(method, lb, ub)
    if isa(params, InfiniteOpt.ParameterRef)
        return _mc_sampling(set, params, num_supports, lb, ub)
    else
        samples_dict = Dict()
        for i in eachindex(lb)
            (samples_dict[i], _) = _mc_sampling(set, params, num_supports, lb[i], ub[i])
        end
        samples = Array{JuMPC.SparseAxisArray, 1}(undef, num_supports)
        for j in 1:num_supports
            samples[j] = JuMPC.SparseAxisArray(Dict(k => samples_dict[k][j]
                                                    for k in eachindex(lb)))
        end
        return (samples, ones(num_supports) / num_supports * prod(ub .- lb))
    end
end

"""
    generate_supports_and_coeffs(set::InfiniteOpt.DistributionSet,
                                 params::Union{InfiniteOpt.ParameterRef,
                                 AbstractArray{<:InfiniteOpt.ParameterRef}},
                                 num_supports::Int,
                                 lb::Union{Number, JuMPC.SparseAxisArray, Nothing},
                                 ub::Union{Number, JuMPC.SparseAxisArray, Nothing},
                                 method::Val{mc_sampling})::Tuple

Return a tuple that contains supports and coefficients generated by Monte Carlo
sampling from a given distribution.

**Example**
```jldoctest; setup = :(using InfiniteOpt, Distributions, Random; model = InfiniteModel(seed = true))
julia> dist = Normal(0., 1.)
Normal{Float64}(μ=0.0, σ=1.0)

julia> @infinite_parameter(model, x in dist)
x

julia> (supps, coeffs) = generate_supports_and_coeffs(DistributionSet(dist), x, 10, nothing, nothing, Val(mc_sampling))
([0.6791074260357777, 0.8284134829000359, -0.3530074003005963, -0.13485387193052173, 0.5866170746331097, 0.29733585084941616, 0.06494754854834232, -0.10901738508171745, -0.514210390833322, 1.5743302021369892], [0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1])
```
"""
function generate_supports_and_coeffs(set::InfiniteOpt.DistributionSet,
                                      params::Union{InfiniteOpt.ParameterRef,
                                      AbstractArray{<:InfiniteOpt.ParameterRef}},
                                      num_supports::Int,
                                      lb::Union{Number, JuMPC.SparseAxisArray, Nothing},
                                      ub::Union{Number, JuMPC.SparseAxisArray, Nothing},
                                      method::Val{mc_sampling})::Tuple
    dist = set.distribution
    # create truncated distribution if necessary
    if !isa(lb, Nothing) || !isa(ub, Nothing)
        if isa(dist, Distributions.MultivariateDistribution)
            @warn("Truncated distribution for multivariate distribution is " *
                  "not supported. Lower bounds and upper bounds are ignored.")
        else
            isa(lb, Number) ? temp_lb = lb : temp_lb = -Inf
            isa(ub, Number) ? temp_ub = ub : temp_ub = Inf
            dist = Distributions.Truncated(dist, temp_lb, temp_ub)
        end
    end
    return _mc_sampling(dist, params, num_supports)
end

"""
    generate_supports_and_coeffs(set::InfiniteOpt.IntervalSet,
                                 params::InfiniteOpt.ParameterRef,
                                 num_supports::Int,
                                 lb::Union{Number, JuMPC.SparseAxisArray, Nothing},
                                 ub::Union{Number, JuMPC.SparseAxisArray, Nothing},
                                 method::Val{gauss_legendre})::Tuple

Return a tuple that contains supports and coefficients generated using
Gauss-Legendre quadrature method. This is useful for univariate parameter in a
finite interval.

**Example**
```jldoctest; setup = :(using InfiniteOpt, JuMP; model = InfiniteModel())
julia> @infinite_parameter(model, t in [0, 1])
t

julia> (supps, coeffs) = generate_supports_and_coeffs(IntervalSet(0,1), t, 5, 0., 1., Val(gauss_legendre))
([0.04691007703066802, 0.23076534494715845, 0.5, 0.7692346550528415, 0.9530899229693319], [0.11846344252809454, 0.23931433524968324, 0.28444444444444444, 0.23931433524968324, 0.11846344252809454])

julia> supps
5-element Array{Float64,1}:
 0.04691007703066802
 0.23076534494715845
 0.5
 0.7692346550528415
 0.9530899229693319
```
"""
function generate_supports_and_coeffs(set::InfiniteOpt.IntervalSet,
                                      params::InfiniteOpt.ParameterRef,
                                      num_supports::Int,
                                      lb::Union{Number, JuMPC.SparseAxisArray, Nothing},
                                      ub::Union{Number, JuMPC.SparseAxisArray, Nothing},
                                      method::Val{gauss_legendre})::Tuple
    _nothing_test(method, lb, ub)
    _univariate_bounds(method, lb, ub)
    (supports, coeffs) = FastGaussQuadrature.gausslegendre(num_supports)
    supports = (ub - lb) / 2 * supports .+ (ub + lb) / 2
    coeffs = (ub - lb) / 2 * coeffs
    return (supports, coeffs)
end

"""
    generate_supports_and_coeffs(set::InfiniteOpt.IntervalSet,
                                 params::InfiniteOpt.ParameterRef,
                                 num_supports::Int,
                                 lb::Union{Number, JuMPC.SparseAxisArray, Nothing},
                                 ub::Union{Number, JuMPC.SparseAxisArray, Nothing},
                                 method::Val{gauss_hermite})::Tuple

Return a tuple that contains supports and coefficients generated using
Gauss-Hermite quadrature method. This is useful for univariate parameter in an
infinite interval.

**Example**
```jldoctest; setup = :(using InfiniteOpt, JuMP; model = InfiniteModel())
julia> @infinite_parameter(model, t in [-Inf, Inf])
t

julia> (supps, coeffs) = generate_supports_and_coeffs(IntervalSet(-Inf,Inf), t, 5, -Inf, Inf, Val(gauss_hermite))
([-2.0201828704560856, -0.9585724646138196, -8.881784197001252e-16, 0.9585724646138196, 2.0201828704560856], [1.1814886255359844, 0.986580996751429, 0.9453087204829428, 0.986580996751429, 1.1814886255359844])

julia> supps
5-element Array{Float64,1}:
 -2.0201828704560856
 -0.9585724646138196
 -8.881784197001252e-16
  0.9585724646138196
  2.0201828704560856
```
"""
function generate_supports_and_coeffs(set::InfiniteOpt.IntervalSet,
                                      params::InfiniteOpt.ParameterRef,
                                      num_supports::Int,
                                      lb::Union{Number, JuMPC.SparseAxisArray, Nothing},
                                      ub::Union{Number, JuMPC.SparseAxisArray, Nothing},
                                      method::Val{gauss_hermite})::Tuple
    _nothing_test(method, lb, ub)
    _univariate_bounds(method, lb, ub)
    if lb != -Inf || ub != Inf
        error("Lower/upper bound is not infinity. Use other measure evaluation " *
              "methods.")
    end
    (supports, coeffs) = FastGaussQuadrature.gausshermite(num_supports)
    coeffs = coeffs .* exp.(supports.^2)
    return (supports, coeffs)
end

"""
    generate_supports_and_coeffs(set::InfiniteOpt.IntervalSet,
                                 params::InfiniteOpt.ParameterRef,
                                 num_supports::Int,
                                 lb::Union{Number, JuMPC.SparseAxisArray, Nothing},
                                 ub::Union{Number, JuMPC.SparseAxisArray, Nothing},
                                 method::Val{gauss_laguerre})::Tuple

Return a tuple that contains supports and coefficients generated using
Gauss-Laguerre quadrature method. This is useful for univariate parameter in a
semi-infinite interval.

**Example**
```jldoctest; setup = :(using InfiniteOpt, JuMP; model = InfiniteModel())
julia> @infinite_parameter(model, t in [-Inf, 0])
t

julia> (supps, coeffs) = generate_supports_and_coeffs(IntervalSet(-Inf,0), t, 5, -Inf, 0, Val(gauss_laguerre))
([-0.2635603197181408, -1.413403059106515, -3.596425771040715, -7.08581000585883, -12.640800844275773], [0.6790940422077494, 1.638487873602747, 2.7694432423708255, 4.3156569009208585, 7.219186354354335])

julia> supps
5-element Array{Float64,1}:
  -0.2635603197181408
  -1.413403059106515
  -3.596425771040715
  -7.08581000585883
 -12.640800844275773
```
"""
function generate_supports_and_coeffs(set::InfiniteOpt.IntervalSet,
                                      params::InfiniteOpt.ParameterRef,
                                      num_supports::Int,
                                      lb::Union{Number, JuMPC.SparseAxisArray, Nothing},
                                      ub::Union{Number, JuMPC.SparseAxisArray, Nothing},
                                      method::Val{gauss_laguerre})::Tuple
    _nothing_test(method, lb, ub)
    _univariate_bounds(method, lb, ub)
    # default Gaussian quadrature method for semi-infinite domain
    if ub == Inf
        if lb == -Inf
            error("The range is infinite. Use other measure evaluation methods.")
        end
        (supports, coeffs) = FastGaussQuadrature.gausslaguerre(num_supports)
        coeffs = coeffs .* exp.(supports)
        supports = supports .+ lb
        return (supports, coeffs)
    elseif lb == -Inf
        (supports, coeffs) = FastGaussQuadrature.gausslaguerre(num_supports)
        coeffs = coeffs .* exp.(supports)
        supports = -supports .+ ub
        return (supports, coeffs)
    else
        error("The range is bounded. Use other measure evaluation methods.")
    end
end

"""
    generate_supports_and_coeffs(set::InfiniteOpt.IntervalSet,
                                 params::InfiniteOpt.ParameterRef,
                                 num_supports::Int,
                                 lb::Union{Number, JuMPC.SparseAxisArray, Nothing},
                                 ub::Union{Number, JuMPC.SparseAxisArray, Nothing},
                                 method::Val{trapezoid})::Tuple

Return a tuple that contains supports and coefficients generated a uniform
trapezoid rule. This is useful for both univariate parameter in a bounded
interval. For multivariate parameters, `num_supports` will be the number of
grid points on one dimension.

**Example**
```jldoctest; setup = :(using InfiniteOpt, JuMP; model = InfiniteModel())
julia> @infinite_parameter(model, t in [0, 1])
t

julia> (supps, coeffs) = generate_supports_and_coeffs(IntervalSet(0,1), t, 3, 0, 1, Val(trapezoid))
([0.0, 0.5, 1.0], [0.25, 0.5, 0.25])

julia> supps
3-element Array{Float64,1}:
 0.0
 0.5
 1.0
```
"""
function generate_supports_and_coeffs(set::InfiniteOpt.IntervalSet,
                                      params::InfiniteOpt.ParameterRef,
                                      num_supports::Int,
                                      lb::Union{Number, JuMPC.SparseAxisArray, Nothing},
                                      ub::Union{Number, JuMPC.SparseAxisArray, Nothing},
                                      method::Val{trapezoid})::Tuple
    _nothing_test(method, lb, ub)
    _univariate_bounds(method, lb, ub)
    return _trapezoid_univ(lb, ub, num_supports)
end

function _trapezoid_univ(lb::Number, ub::Number, num_supports::Int)::Tuple
    increment = (ub - lb) / (num_supports - 1)
    supps = [(i - 1) * increment + lb for i in 1:num_supports]
    coeffs = increment * ones(num_supports)
    coeffs[1] = coeffs[1] / 2
    coeffs[num_supports] = coeffs[num_supports] / 2
    return (supps, coeffs)
end

# MC sampling for univariate IntervalSet
function _mc_sampling(set::InfiniteOpt.IntervalSet,
                      params::Union{InfiniteOpt.ParameterRef,
                      AbstractArray{<:InfiniteOpt.ParameterRef}},
                      num_supports::Int,
                      lb::Number, ub::Number)::Tuple
    # MC sampling from uniform distribution over the interval [lb, ub]
    if lb == -Inf || ub == Inf
        return infinite_transform(set, params, num_supports, lb, ub, Val(mc_sampling))
    else
        samples = rand(num_supports) .* (ub - lb) .+ lb
        return (samples, ones(num_supports) / num_supports * (ub - lb))
    end
end

# MC sampling for univariate DistributionSet
function _mc_sampling(dist::Distributions.UnivariateDistribution,
                     param::InfiniteOpt.ParameterRef,
                     num_supports::Int)::Tuple
    # MC sampling - Distribution Set
    samples = rand(dist, num_supports)
    return (samples, ones(num_supports) / num_supports)
end

# MC sampling for multivariate DistributionSet
function _mc_sampling(dist::Distributions.MultivariateDistribution,
                     params::AbstractArray{<:InfiniteOpt.ParameterRef},
                     num_supports::Int)::Tuple
    samples_matrix = rand(dist, num_supports)
    ordered_pairs = sort(collect(params.data), by=x->x.second.index)
    samples_dict = Dict()
    for i in eachindex(ordered_pairs)
        samples_dict[ordered_pairs[i].first] = samples_matrix[i, :]
    end
    samples = Array{JuMPC.SparseAxisArray, 1}()
    for j in 1:num_supports
        append!(samples, [JuMP.Containers.SparseAxisArray(
                          Dict(k.first => samples_dict[k.first][j]
                          for k in ordered_pairs))])
    end
    return (samples, ones(num_supports) / num_supports)
end

function _nothing_test(method::Val,
                       lb::Union{Number, JuMPC.SparseAxisArray, Nothing},
                       ub::Union{Number, JuMPC.SparseAxisArray, Nothing})
    if isa(lb, Nothing) || isa(ub, Nothing)
        error("Method " * string(method)[6:length(string(method))-3] *
              " cannot be nothing.")
    end
end

# test if the lower/upper bound is univariate
function _univariate_bounds(method::Val,
                            lb::Union{Number, JuMPC.SparseAxisArray, Nothing},
                            ub::Union{Number, JuMPC.SparseAxisArray, Nothing})
    if isa(lb, JuMPC.SparseAxisArray) || isa(ub, JuMPC.SparseAxisArray)
        error("Method " * string(method)[6:length(string(method))-3] *
              " is not applicable to multivariate infinite parameters.")
    end
end

"""
    infinite_transform(lb::Number, ub::Number, num_supports::Int;
                       [sub_method::Function = mc_sampling,
                       transform_x::Function = _default_x,
                       transform_dx::Function = _default_dx,
                       t_lb::Number = -convert(Number, lb == -Inf && ub == Inf),
                       t_ub::Number = 1., kwargs...])::Tuple

Returns a tuple that contains supports and coefficients generated for a
parameter in an infinite or semi-infinite interval. It works by transforming
the original unbounded interval to a finite interval, on which a support
generation method for finite intervals is applied. Then, the generated supports
are transformed back to the original interval. The user is allowed to specify
the support generation method for finite intevals to use, as well as the
transform function. The default transform function is
``t \\in [-\\infty, \\infty] \\rightarrow x \\in [-1, 1]: t(x) = \\frac{t}{1-t^2}``
``t \\in [a, \\infty] \\rightarrow x \\in [0, 1]: t(x) = a + \\frac{t}{1-t}``
``t \\in [-\\infty, a] \\rightarrow x \\in [0, 1]: t(x) = a - \\frac{1-t}{t}``

**Example**
```jldoctest; setup = :(using InfiniteOpt)
julia> (supps, coeffs) = infinite_transform(-Inf, Inf, 5, sub_method = gauss_legendre)
([-5.06704059565454, -0.7583532171678754, 0.0, 0.7583532171678754, 5.06704059565454], [13.490960583398396, 1.2245949721571516, 0.5688888888888889, 1.2245949721571516, 13.490960583398396])
```
"""
function infinite_transform(set::InfiniteOpt.IntervalSet,
                            params::Union{InfiniteOpt.ParameterRef,
                            AbstractArray{<:InfiniteOpt.ParameterRef}},
                            num_supports::Int,
                            lb::Union{Number, JuMPC.SparseAxisArray, Nothing},
                            ub::Union{Number, JuMPC.SparseAxisArray, Nothing},
                            sub_method::Val = Val(mc_sampling);
                            transform_x::Function = _default_x,
                            transform_dx::Function = _default_dx,
                            t_lb::Number = -convert(Number, lb == -Inf && ub == Inf),
                            t_ub::Number = 1.)::Tuple
    # transform (semi-)infinite domain to finite domain
    if lb != -Inf && ub != Inf
        error("The range is not (semi-)infinite. Use evaluation methods for " *
              "bounded domains.")
    end
    (t_supports, t_coeffs) = generate_supports_and_coeffs(set, params, num_supports, t_lb, t_ub, sub_method)
    supports = transform_x.(t_supports, lb, ub)
    coeffs = t_coeffs .* transform_dx.(t_supports, lb, ub)
    return (supports, coeffs)
end

function _default_x(t::Number, lb::Number, ub::Number)::Number
    if lb > -Inf
        return lb + t / (1 - t)
    elseif ub < Inf
        return ub - (1 - t) / t
    else
        return t / (1 - t^2)
    end
end

function _default_dx(t::Number, lb::Number, ub::Number)::Number
    if lb > -Inf
        return 1 / (1 - t)^2
    elseif ub < Inf
        return 1 / t^2
    else
        return (1 + t^2) / (1 - t^2)^2
    end
end

# TODO: consider truncated distribution
# TODO: consider adding uniform grids


text:
# Basic library functionality

# TODO: put error here to avoid need for forward decl

#
# API versioning
#

const mapping = Dict{Symbol,Symbol}()
const minreq = Dict{Symbol,VersionNumber}()

if libcuda_version >= v"3.2"
    mapping[:cuDeviceTotalMem]           = :cuDeviceTotalMem_v2
    mapping[:cuCtxCreate]                = :cuCtxCreate_v2
    mapping[:cuModuleGetGlobal]          = :cuModuleGetGlobal_v2
    mapping[:cuMemGetInfo]               = :cuMemGetInfo_v2
    mapping[:cuMemAlloc]                 = :cuMemAlloc_v2
    mapping[:cuMemAllocPitch]            = :cuMemAllocPitch_v2
    mapping[:cuMemFree]                  = :cuMemFree_v2
    mapping[:cuMemGetAddressRange]       = :cuMemGetAddressRange_v2
    mapping[:cuMemAllocHost]             = :cuMemAllocHost_v2
    mapping[:cuMemHostGetDevicePointer]  = :cuMemHostGetDevicePointer_v2
    mapping[:cuMemcpyHtoD]               = :cuMemcpyHtoD_v2
    mapping[:cuMemcpyDtoH]               = :cuMemcpyDtoH_v2
    mapping[:cuMemcpyDtoD]               = :cuMemcpyDtoD_v2
    mapping[:cuMemcpyDtoA]               = :cuMemcpyDtoA_v2
    mapping[:cuMemcpyAtoD]               = :cuMemcpyAtoD_v2
    mapping[:cuMemcpyHtoA]               = :cuMemcpyHtoA_v2
    mapping[:cuMemcpyAtoH]               = :cuMemcpyAtoH_v2
    mapping[:cuMemcpyAtoA]               = :cuMemcpyAtoA_v2
    mapping[:cuMemcpyHtoAAsync]          = :cuMemcpyHtoAAsync_v2
    mapping[:cuMemcpyAtoHAsync]          = :cuMemcpyAtoHAsync_v2
    mapping[:cuMemcpy2D]                 = :cuMemcpy2D_v2
    mapping[:cuMemcpy2DUnaligned]        = :cuMemcpy2DUnaligned_v2
    mapping[:cuMemcpy3D]                 = :cuMemcpy3D_v2
    mapping[:cuMemcpyHtoDAsync]          = :cuMemcpyHtoDAsync_v2
    mapping[:cuMemcpyDtoHAsync]          = :cuMemcpyDtoHAsync_v2
    mapping[:cuMemcpyDtoDAsync]          = :cuMemcpyDtoDAsync_v2
    mapping[:cuMemcpy2DAsync]            = :cuMemcpy2DAsync_v2
    mapping[:cuMemcpy3DAsync]            = :cuMemcpy3DAsync_v2
    mapping[:cuMemsetD8]                 = :cuMemsetD8_v2
    mapping[:cuMemsetD16]                = :cuMemsetD16_v2
    mapping[:cuMemsetD32]                = :cuMemsetD32_v2
    mapping[:cuMemsetD2D8]               = :cuMemsetD2D8_v2
    mapping[:cuMemsetD2D16]              = :cuMemsetD2D16_v2
    mapping[:cuMemsetD2D32]              = :cuMemsetD2D32_v2
    mapping[:cuArrayCreate]              = :cuArrayCreate_v2
    mapping[:cuArrayGetDescriptor]       = :cuArrayGetDescriptor_v2
    mapping[:cuArray3DCreate]            = :cuArray3DCreate_v2
    mapping[:cuArray3DGetDescriptor]     = :cuArray3DGetDescriptor_v2
    mapping[:cuTexRefSetAddress]         = :cuTexRefSetAddress_v2
    mapping[:cuTexRefGetAddress]         = :cuTexRefGetAddress_v2
    mapping[:cuGraphicsResourceGetMappedPointer] = :cuGraphicsResourceGetMappedPointer_v2
end

if libcuda_version >= v"4.0"
    mapping[:cuCtxDestroy]               = :cuCtxDestroy_v2
    mapping[:cuCtxPopCurrent]            = :cuCtxPopCurrent_v2
    mapping[:cuCtxPushCurrent]           = :cuCtxPushCurrent_v2
    mapping[:cuStreamDestroy]            = :cuStreamDestroy_v2
    mapping[:cuEventDestroy]             = :cuEventDestroy_v2
end

if libcuda_version >= v"4.1"
    mapping[:cuTexRefSetAddress2D]       = :cuTexRefSetAddress2D_v3
end

if libcuda_version >= v"6.5"
    mapping[:cuLinkCreate]              = :cuLinkCreate_v2
    mapping[:cuLinkAddData]             = :cuLinkAddData_v2
    mapping[:cuLinkAddFile]             = :cuLinkAddFile_v2
end

if libcuda_version >= v"6.5"
    mapping[:cuMemHostRegister]         = :cuMemHostRegister_v2
    mapping[:cuGraphicsResourceSetMapFlags] = :cuGraphicsResourceSetMapFlags_v2
end

if v"3.2" <= libcuda_version < v"4.1"
    mapping[:cuTexRefSetAddress2D]      = :cuTexRefSetAddress2D_v2
end

## Version-dependent features

minreq[:cuLinkCreate]       = v"5.5"
minreq[:cuLinkDestroy]      = v"5.5"
minreq[:cuLinkComplete]     = v"5.5"
minreq[:cuLinkAddFile]      = v"5.5"
minreq[:cuLinkAddData]      = v"5.5"

minreq[:cuDummyAvailable]   = v"0"      # non-existing functions
minreq[:cuDummyUnavailable] = v"999"    # for testing purposes

# explicitly mark unavailable symbols, signaling `@apicall` to error out
for (api_function, minimum_version) in minreq
    if libcuda_version < minimum_version
        mapping[api_function] = :unavailable
    end
end

#
# API call wrapper
#

# forward definition of CuError
const CuError_t = Cint
"""
    CuError(code::Integer)
    CuError(code::Integer, info::String)

Create a CUDA error object with error code `code`. The optional `info` parameter indicates
whether extra information, such as error logs, is known.
"""
struct CuError <: Exception
    code::CuError_t
    meta::Any

    CuError(code, meta=nothing) = new(code, meta)
end

# ccall wrapper for calling functions in NVIDIA libraries
const apicall_hook = Ref{Union{Nothing,Function}}(nothing)
macro apicall(funspec, argtypes, args...)
    isa(funspec, QuoteNode) || error("first argument to @apicall should be a symbol")
    fun = funspec.value

    # resolve the function
    global mapping, version_requirements
    api_fun = get(mapping, fun, fun)
    if api_fun == :unavailable
        return :(throw(CuVersionError($(QuoteNode(fun)), $(minreq[fun]))))
    end

    configured || return :(error("CUDAdrv.jl has not been configured. Please run Pkg.build(\"CUDAdrv\") and restart Julia."))

    return quote
        # NOTE: this hook is used by CUDAnative.jl to initialize upon the first API call
        apicall_hook[] !== nothing && apicall_hook[]($funspec)

        status = ccall(($(QuoteNode(api_fun)), libcuda), CuError_t,
                       $(esc(argtypes)), $(map(esc, args)...))

        if status != SUCCESS.code
            err = CuError(status)
            throw(err)
        end
    end
end

#
# Basic functionality
#

"""
Returns a string identifying the vendor of your CUDA driver.
"""
function vendor()
    return libcuda_vendor
end


text:
using LoopVectorization, StrideArraysCore, Test

function issue238_noavx!(output, matrix, input)
  for i in axes(output, 2), v in axes(output, 1)
    res = zero(eltype(output))
    for ii in axes(matrix, 2)
      res += matrix[i, ii] * input[v, ii]
    end
    output[v, i] = res
  end
  return nothing
end
function issue238!(output, matrix, input)
  @turbo inline = true for i in axes(output, 2), v in axes(output, 1)
    res = zero(eltype(output))
    for ii in axes(matrix, 2)
      res += matrix[i, ii] * input[v, ii]
    end
    output[v, i] += res
  end
  return nothing
end
function issue238_v2!(output, matrix, input)
  @turbo inline = true for i in axes(output, 2), v in axes(output, 1)
    res = output[v, i]
    for ii in axes(matrix, 2)
      res += matrix[i, ii] * input[v, ii]
    end
    output[v, i] = res
  end
  return nothing
end
const MAXTESTSIZE = 8 # N^3 loop iterations, meaning N^3 different functions compiled
function n2testloop(output1, output2, output3, output_nonstatic0, output_nonstatic1)
  n1, n3 = size(output1)
  for n2 ∈ 1:MAXTESTSIZE
    # @show n1, n2, n3
    input = randn(n1, n2)
    matrix = randn(n3, n2)
    smatrix = StrideArray(matrix, (StaticInt(n1), StaticInt(n2)))
    issue238_noavx!(output_nonstatic0, matrix, input)
    fill!(output_nonstatic1, 0)
    issue238!(output_nonstatic1, matrix, input)
    fill!(output1, 0)
    issue238!(output1, matrix, input)
    fill!(output2, 0)
    issue238!(output2, smatrix, input)
    fill!(output3, 0)
    issue238_v2!(output3, smatrix, input)
    if !(output1 ≈ output2 ≈ output3 ≈ output_nonstatic0 ≈ output_nonstatic1)
      @show n1, n2, n3
    end
    @test output_nonstatic0 ≈ output1 # static test
    @test output_nonstatic0 ≈ output2 # static test
    @test output_nonstatic0 ≈ output3 # static test
    @test output_nonstatic0 ≈ output_nonstatic1 # nonstatic test
  end
end
function update!(B⁻¹yₖ, B⁻¹, yₖ, sₖᵀyₖ⁻¹)
  yₖᵀB⁻¹yₖ = zero(eltype(B⁻¹))
  @inbounds @fastmath for c ∈ axes(B⁻¹, 2)
    t = zero(yₖᵀB⁻¹yₖ)
    for r ∈ axes(B⁻¹, 1)
      t += yₖ[r] * B⁻¹[r, c]
    end
    B⁻¹yₖ[c] = t * sₖᵀyₖ⁻¹
    yₖᵀB⁻¹yₖ += t * yₖ[c]
  end
  yₖᵀB⁻¹yₖ
end
function update_turbo!(B⁻¹yₖ, B⁻¹, yₖ, sₖᵀyₖ⁻¹)
  yₖᵀB⁻¹yₖ = zero(eltype(B⁻¹))
  @turbo for c ∈ axes(B⁻¹, 2)
    t = zero(yₖᵀB⁻¹yₖ)
    for r ∈ axes(B⁻¹, 1)
      t += yₖ[r] * B⁻¹[r, c]
    end
    B⁻¹yₖ[c] = t * sₖᵀyₖ⁻¹
    yₖᵀB⁻¹yₖ += t * yₖ[c]
  end
  yₖᵀB⁻¹yₖ
end

function maxabs(x)
  s = -Inf
  @turbo for i ∈ eachindex(x)
    s = max(s, abs(x[i]))
  end
  s
end
function sum_turbo(x)
  s = zero(eltype(x))
  @turbo for i ∈ eachindex(x)
    s += x[i]
  end
  s
end

@testset "Statically Sized Arrays" begin
  @show @__LINE__
  for n1 ∈ 1:MAXTESTSIZE, n3 ∈ 1:MAXTESTSIZE
    # @show n1, n3
    output1 = StrideArray(undef, StaticInt(n1), StaticInt(n3))
    output2 = StrideArray(undef, StaticInt(n1), StaticInt(n3))
    output3 = StrideArray(undef, StaticInt(n1), StaticInt(n3))
    output_nonstatic0 = StrideArray(undef, n1, n3)
    output_nonstatic1 = StrideArray(undef, n1, n3)
    n2testloop(output1, output2, output3, output_nonstatic0, output_nonstatic1)

    y = StrideArray(undef, StaticInt(max(n1, n3)))
    y .= rand.()
    By0 = StrideArray(undef, StaticInt(n3))
    By1 = StrideArray(undef, StaticInt(n3))
    GC.@preserve By0 By1 output1 y begin
      @test update_turbo!(By0, output1, y, 0.124) ≈ update!(By1, output1, y, 0.124)
      @test By0 ≈ By1
    end
  end
  for i = 1:65
    # @show i
    x = StrideArray(undef, StaticInt(i))
    x .= randn.()
    @test maxabs(x) == maximum(abs, x)
    @test sum_turbo(x) ≈ sum(x)
  end
end


text:
#!/usr/bin/env julia

using SEM
using Plots,LinearAlgebra

n = 8
x,w = SEM.FastGaussQuadrature.gausslobatto(n)
xo  = SEM.linspace(-1,1,10*n)

J = SEM.interpMat(xo,x)
D = SEM.derivMat(x)

# neumann operators
B = diagm(w)
A = D'*B *D
C(v) = B * diagm(v) * D

## CASE SETUP
function utrue(x,v,t)
    xx = x - v*t
    return cos(2pi*xx)
#   return exp(-5*x*x)
end
f = @. 0+0*x
vl= 1.0
v = @. vl+0*x
ν = 0e-0
u = utrue.(x,vl,0.0)
id= Matrix(I,n,n)

#R = id[2:end-1,:] # double dirichlet
R = id[1:end-0,:] # periodic
Q = semq(1,n,true)

CFL = 0.1
dx = minimum(diff(x))
dt = CFL * dx / vl
println("CFL=$CFL, dt=$dt")

## global system (full rank)
xx = R*x
AA = R*A*R' .* ν
BB = R*B
CC = R*C(v)*R'
bb = BB*f
AC = AA + CC

## local system (rank deficient)
b  = R'*bb
Al = R'*AA*R
Bl = R'*BB
Cl = R'*CC*R
Hl = Al + Bl

#println("er: ",norm(u - ut,Inf))
p=plot(xo,J*u,width=3,ylims=(-1.5,1.5))
display(p)

uh = [zero(u) for i in 1:3] # histories
gh = [zero(u) for i in 1:3] # explicit term

T = 1.0
t = zeros(4,1)
nstep = Int(floor(T/dt))

for istep=1:nstep
    global u,rhs,Hl,bdfA,bdfB

    updateHist!(t)
    updateHist!(u,uh)

    t[1] += dt
    bdfA,bdfB = bdfExtK(t)

    gh[1] = -Cl*uh[1] # convection

    rhs = Bl * f
    for i=1:3
        gh[i] = -Cl*uh[i]                 # convection
        rhs  += -bdfB[1+i] .* (Bl *uh[i]) # histories
        rhs  +=  bdfA[i]   .* gh[i]       # explicit term
    end

    Hl = Al + bdfB[1]*Bl # Helmholtz op

    Hl  = Q*Q'*Hl
    rhs = Q*Q'*rhs
    u = pcg(rhs,Hl)

    ut = utrue.(x,vl,t[1])
    er = norm(u .- ut,Inf)
    println("Pointwise error: $er")

    plt = plot(xo,J*u,width=3)
    plt = plot!(title="Step $istep, Time $(t[1])")
    plt = plot!(ylims=(-1.5,1.5))
    display(plt)
end
println("CFL=$CFL, dt=$dt")
#----------------------------------------------------------------------#

#prob = ODEProblem(dudt,u0,tspn)
#sol  = solve(prob,SSPRK33(),dt=0.01,saveat=0.1,cb=callback)
#----------------------------------------------------------------------#
return


text:
"""
    IndicatorActiveOnFalseBridge{T}

The `IndicatorActiveOnFalseBridge` replaces an indicator constraint activated
on 0 with a variable ``z_0`` with the constraint activated on 1, with a variable ``z_1``.
It stores the added `variable_index` and added constraints:
- ``z_1 \\in \\mathbb{B}`` in `zero_one_cons`
- ``z_0 + z_1 == 1`` in `` in `disjunction_cons`
- The added `ACTIVATE_ON_ONE` indicator constraint in `indicator_cons_index`.
"""
struct IndicatorActiveOnFalseBridge{T, F <: MOI.AbstractVectorFunction, S <: MOI.AbstractScalarSet} <: AbstractBridge
    variable_index::MOI.VariableIndex
    zero_one_cons::MOI.ConstraintIndex{MOI.SingleVariable, MOI.ZeroOne}
    disjunction_cons::MOI.ConstraintIndex{MOI.ScalarAffineFunction{T}, MOI.EqualTo{T}}
    indicator_cons_index::MOI.ConstraintIndex{F, MOI.IndicatorSet{MOI.ACTIVATE_ON_ONE, S}}
end

function bridge_constraint(::Type{<:IndicatorActiveOnFalseBridge}, model::MOI.ModelLike, f::MOI.VectorAffineFunction{T}, s::IS) where {S <: MOI.AbstractScalarSet, T <: Real, F, IS <: MOI.IndicatorSet{MOI.ACTIVATE_ON_ZERO, S}}
    z1 = f.terms[1].scalar_term.variable_index
    z2 = MOI.add_variable(model)
    zo_cons = MOI.add_constraint(model, MOI.SingleVariable(z2), MOI.ZeroOne())

    # z1 + z2 == 1
    dcons = MOI.add_constraint(model,
        MOI.ScalarAffineFunction(
            [MOI.ScalarAffineTerm(one(T), z1), MOI.ScalarAffineTerm(one(T), z2)], zero(T),
        ),
        MOI.EqualTo(one(T)),
    )
    vec_terms2 = [t for t in f.terms if t.output_index == 2]
    f2 = MOI.VectorAffineFunction(
        vcat(
            [MOI.VectorAffineTerm(1, MOI.ScalarAffineTerm(1.0, z2))],
            vec_terms2,
        ),
        f.constants
    )
    ci = MOI.add_constraint(model, f2, MOI.IndicatorSet{MOI.ACTIVATE_ON_ONE}(s.set))

    return IndicatorActiveOnFalseBridge{T, MOI.VectorAffineFunction{T}, S}(z2, zo_cons, dcons, ci)
end

function MOI.supports_constraint(::Type{<:IndicatorActiveOnFalseBridge},
                                 ::Type{<:MOI.VectorAffineFunction},
                                 ::Type{<:MOI.IndicatorSet{MOI.ACTIVATE_ON_ZERO}})
    return true
end

function added_constraint_types(::Type{IndicatorActiveOnFalseBridge{T, F, S}}) where {T, F, S}
    return [(F, MOI.IndicatorSet{MOI.ACTIVATE_ON_ONE, S})]
end

function concrete_bridge_type(::Type{<:IndicatorActiveOnFalseBridge{T}},
                              ::Type{F},
                              ::Type{MOI.IndicatorSet{MOI.ACTIVATE_ON_ZERO, S}}) where {T, F, S}
    return IndicatorActiveOnFalseBridge{T, F, S}
end


text:
using Random
include("Params.jl")

struct Solution
    params::Params               # Access to the problem and dataset parameters
    permutation::Array{Int64}    # Permutation
    cost::Int64                  # Value of objective function for the solution
end

# Constructor of struct Solution (random permutation)
function Solution(params::Params)::Solution
    # p = convert(Array{Int64}, transpose(randperm(params.datasetSize)))
    p = randperm(params.datasetSize)
    sol = Solution(params, p)
    return sol
end

# Constructor of struct Solution given a permutation
function Solution(params::Params, p::Array{Int64})::Solution
    cost = evaluate(params, p)
    sol = Solution(params, p, cost)
    return sol
end

# Evaluate solution (calculate value of objective function)
function evaluate(params::Params, p::Array{Int64})::Float64
    n = params.datasetSize
    A = params.A
    B = params.B
    cost = sum(A[i,j]*B[p[i],p[j]] for i = 1:n, j = 1:n)
    if cost < 0
        throw("Cost should be non negative!")
    end
    return cost
end

# Print solution permutation and value of objective function (cost)
function printSolution(sol::Solution)
    println("Permutation=",sol.permutation)
    println("Cost=",sol.cost)
end

# Generate neighborhood of solutions (array of Solutions) obtained from
# all pairwise exchanges of "facilities"
function generateSwapNeighbors(sol::Solution)::Array{Solution}
    n = sol.params.datasetSize
    neighborhood = Array{Solution}(undef, Int((n-1)*n/2))
    k = 1
    for i = 1:n-1
        for j = i+1:n
            newP = copy(sol.permutation)
            # Exchange facilities
            tmp = newP[i]
            newP[i] = newP[j]
            newP[j] = tmp
            # Add new solution to neighborhood
            neighborhood[k] = Solution(sol.params,newP)
            k += 1
        end
    end
    return neighborhood
end

function generateCyclicNeighbors(sol::Solution)::Array{Solution}
    n = sol.params.datasetSize
    neighborhood = Array{Solution}(undef, n)
    for k=1:n
        newP = circshift(sol.permutation, k)
        neighborhood[k] = Solution(sol.params,newP)
    end
    return neighborhood
end

import Base: isless

isless(a::Solution, b::Solution) = isless(a.cost, b.cost)

max_stars_repo_name: BSAraujo/QAP-Metaheuristics

max_issues_repo_name: BSAraujo/QAP-Metaheuristics

max_forks_repo_name: BSAraujo/QAP-Metaheuristics


text:
# Please follow the longer description in the documentation:
# https://csertegt3.github.io/RANSAC.jl/stable/example/

using FileIO
m = load("fandisk_input.obj");

using RANSACVisualizer
showgeometry(m.position, m.normals)

using RANSAC
pc = RANSACCloud(m.position, m.normals, 8)
p = ransacparameters()
newparams = (ϵ=0.05, α=deg2rad(10),)
p = ransacparameters(p, sphere=newparams, cone=newparams,
    plane=newparams, cylinder=newparams, iteration=(τ=50, itermax=100_000,))

extr, _ = ransac(pc, p, true)

showbytype(pc, extr; show_axis = false)


text:
@testset "building blocks" begin
    A = rand(5, 5)
    x = rand(5)
    R = [1:2, 3:5]
    @test _basis_block(A, R, (1, 2), (2, 1)) ==
        kron(A[R[2], R[1]], A[R[1], R[2]])
    @test @inferred(_basis_block(A, R, (1, 2, 2), (2, 1, 1))) ≈
        kron(kron(A[R[2], R[1]],
                  A[R[2], R[1]]),
             A[R[1], R[2]])
    @test @inferred(_coordinate_block(x, R, (1, 2, 2))) ==
        tuple.(repeat(x[R[1]], outer = 9),
               repeat(x[R[2]], inner = 2, outer = 3),
               repeat(x[R[2]], inner = 6))
end

@testset "basis and coordinates" begin
    B = HomogeneousBasis(Chebyshev(), (4, 4), 2)
    A, v = basis_matrix_and_coordinates(B)

    # consistency checks
    d = degrees_of_freedom(B)
    @test size(A) == (d, d)
    @test size(v) == (d, )
    @test eltype(v) == SVector{2, Float64}

    # compare to manually constructed matrix
    R = set_range.(Chebyshev(), 1:3)
    x0 = all_gridpoints(Chebyshev(), 4)
    A0 = evaluate(9, Chebyshev(), x0)

    iter = CappedCartesianIndices(4, (4, 4))
    let row_sum = 0             # cumulative row indx
        for i_row in iter
            col_sum = 0         # cumulative column index
            row_s = 0           # when ≠0, (row) size of blocks
            for i_col in iter
                a = fill(1.0, 1, 1)
                for (i, j) in zip(i_row, i_col)
                    a = kron(A0[R[i], R[j]], a)
                end
                s1, s2 = size(a)
                if row_s == 0
                    row_s = s1      # first, just save
                else
                    @test row_s == s1
                end
                @test a ≈ A[row_sum .+ (1:s1), col_sum .+ (1:s2)]
                col_sum += s2
            end
            row_sum += row_s
            @test col_sum == d
        end
        @test row_sum == d
    end
    let p = 0
        for i in iter
            for j in CartesianIndices(map(i -> R[i], i))
                p += 1
                @test v[p] ≈ SVector(map(j -> x0[j], Tuple(j)))
            end
        end
        @test p == d
    end

    # interpolation
    f(x) = sin(x[1]) + cos(x[2] * 1.7)
    fx = f.(v)
    c = A \ fx

    # test at exactly interpolated points
    for i in 1:d
        @test interpolate(B, c, v[i]) ≈ fx[i]
        a = interpolated_basis(B, v[i])
        @test length(a) == d
        @test dot(a, c) ≈ fx[i]
    end

    # test at other points
    xs = [rand(2) .* 2 .- 1 for _ in 1:50]
    for x in xs
        f̃x = interpolate(B, c, x)
        @test f̃x ≈ f(x) atol = 1e-2
        a = interpolated_basis(B, x)
        @test length(a) == d
        @test dot(a, c) ≈ f̃x
    end
end


text:
using Poirot, Test, Random
using Poirot: Singleton

Random.seed!(0)

@test infer(() -> 1) == Singleton(1)
@test infer(() -> "foo") == Singleton("foo")

@test infer(() -> rand()) == Uniform(0, 1)
@test infer(() -> randn()) == Normal(0, 1)
@test infer(() -> rand(Bool)) == Bernoulli(1//2)

coin() = rand(Bool)

@test infer(() -> coin() & coin()) == Bernoulli(0.25)
@test infer(() -> coin() | coin()) == Bernoulli(0.75)

@test infer() do
  a = coin()
  b = coin()
  observe(a | b)
  a & b
end == Bernoulli(1/3)

infer() do
  a = coin()
  b = coin()
  observe(a | b)
  observe(!b)
  a
end == Singleton(true)

@test infer() do
  x = randn()
  observe(x > 0)
  x
end isa Poirot.Empirical

d = infer() do
  x = 10*rand()
  x^2 > 50
end

@test d isa Bernoulli
@test mean(d) ≈ (10-sqrt(50))/10

d = infer() do
  x = 10*rand()
  observe(x > 5)
  x^2 > 50
end

@test d isa Bernoulli
@test mean(d) ≈ (10-sqrt(50))/5


text:
using PDMesh
using Test

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

max_stars_repo_name: sajidmannan/PDMesh

max_issues_repo_name: sajidmannan/PDMesh

max_forks_repo_name: sajidmannan/PDMesh


text:
function DNDR(ALPHA,BETA)
      A=[-0.0180   -0.0280   -0.0370   -0.0480   -0.0430   -0.0520   -0.0620
         -0.0520   -0.0510   -0.0410   -0.0450   -0.0440   -0.0340   -0.0340
         -0.0520   -0.0430   -0.0380   -0.0450   -0.0410   -0.0360   -0.0270
         -0.0520   -0.0460   -0.0400   -0.0450   -0.0410   -0.0360   -0.0280
         -0.0540   -0.0450   -0.0400   -0.0440   -0.0400   -0.0350   -0.0270
         -0.0490   -0.0490   -0.0380   -0.0450   -0.0380   -0.0280   -0.0270
         -0.0590   -0.0570   -0.0370   -0.0470   -0.0340   -0.0240   -0.0230
         -0.0510   -0.0520   -0.0300   -0.0480   -0.0350   -0.0230   -0.0230
         -0.0300   -0.0300   -0.0270   -0.0490   -0.0350   -0.0200   -0.0190
         -0.0370   -0.0330   -0.0240   -0.0450   -0.0290   -0.0160   -0.0090
         -0.0260   -0.0300   -0.0190   -0.0330   -0.0220   -0.0100   -0.0250
         -0.0130   -0.0080   -0.0130   -0.0160   -0.0090   -0.0140   -0.0100]
      #
      S= 0.2 * ALPHA
      K= floor(S)+3
      if(K <= 1)
            K= 2
      end
      if(K >=  12)
            K=  11
      end
      DA= S - K +3
      L = K + fix( 1.1*sign(DA) )
      #
      S= 0.1 * BETA
      M= floor(S)+4
      if(M <= 1)
            M= 2
      end
      if(M >=  7)
            M=  6
      end
      DB= S - M + 4
      N= M + fix( 1.1*sign(DB) )
      #
      K=convert(Int64, K)
      L=convert(Int64, L)
      M=convert(Int64, M)
      N=convert(Int64, N)
      T= A[K,M]
      U= A[K,N]
      V= T + abs(DA) * (A[L,M] - T)
      W= U + abs(DA) * (A[L,N] - U)
      return dndrv = V + (W-V)  * abs(DB)
end


text:
"""

-T paranoid|sneaky|polite|normal|aggressive|insane

See also: [Nmap Timing and Performance](https://nmap.org/book/man-performance.html)
"""
@enum Timing begin
    paranoid    = 0
    sneaky      = 1
    polite      = 2
    normal      = 3
    aggressive  = 4
    insane      = 5
end

const AllowedTimings = Union{Timing, Int}

function timingtemplate(timing::AllowedTimings) :: Option
    @assert Int(timing) in Int.(instances(Timing)) "Timing template values range from 0 to 5"
    return Option(
        function(scan::Scanner)
            push!(scan.args, string("-T", Int(timing)))
        end
    )
end


text:
@testset "weights ui" begin
    x = randn()

    um1 = ui(-1)
    woff = Float64[]
    compute_weights(um1, 2.0, woff)
    @test woff ==Float64[]

    u4 = ui(0, zeros(0), ones(0), [0.0])

    woff = zeros(1)

    compute_weights(u4, 2.0, woff)
    @test woff[1]==2.0

    u4 = ui(2, zeros(2), ones(2), [0.0; 0.0; 0.0])

    woff = zeros(3)

    compute_weights(u4, x, woff)
    @test isapprox(woff[1], x, atol = 1e-10)
    @test isapprox(woff[2], rbf(0.0, 1.0)(x), atol = 1e-10)
    @test isapprox(woff[3], rbf(0.0, 1.0)(x), atol = 1e-10)

    u4.ξi[1] = -5.0
    u4.σi[1] = 5.0

    u4.ξi[2] = 2.0
    u4.σi[2] = 2.0

    compute_weights(u4, x, woff)
    @test isapprox(woff[1], x, atol = 1e-10)
    @test isapprox(woff[2], rbf(-5.0, 5.0)(x), atol = 1e-10)
    @test isapprox(woff[3], rbf(2.0, 2.0)(x), atol = 1e-10)
end

@testset "weights uk" begin
    x = randn()
    um1 = uk(-1)
    wdiag = Float64[]
    w∂k = Float64[]
    compute_weights(um1, 2.0, wdiag, w∂k)
    @test wdiag ==Float64[]
    @test w∂k ==Float64[]

    u1 = uk(0)
    wdiag = zeros(2)
    w∂k = zeros(1)
    compute_weights(u1, 2.0, wdiag, w∂k)
    @test wdiag == [1.0; 2.0]
    @test w∂k[1] === 1.0

    u4 = uk(0)
    wdiag = zeros(2)
    w∂k = zeros(1)

    compute_weights(u4, 2.0, wdiag, w∂k)

    @test wdiag == [1.0; 2.0]
    @test w∂k[1] === 1.0

    u4 = uk(3)
    wdiag = zeros(3+3)
    w∂k = zeros(3+2)

    compute_weights(u4, x, wdiag, w∂k)
    @test size(wdiag,1)== 3+3
    @test size(w∂k,1)== 3+2
    @test norm(wdiag - [1.0;
                        ψ₀(0.0, 1.0)(x);
                        ψj(0.0, 1.0)(x);
                        ψj(0.0, 1.0)(x);
                        ψj(0.0, 1.0)(x);
                        ψpp1(0.0, 1.0)(x)])<1e-10

    @test norm(w∂k -[ψ₀′(0.0, 1.0)(x);
                     rbf(0.0, 1.0)(x);
                     rbf(0.0, 1.0)(x);
                     rbf(0.0, 1.0)(x);
                   ψpp1′(0.0, 1.0)(x)])<1e-10

    u4 = uk(3)
    wdiag = zeros(3+3)
    w∂k = zeros(3+2)

    u4.ξk .= [-1.0; 2.0; -3.0; 4.0; -5.0]
    u4.σk .= [1.0; 2.0; 3.0; 4.0; 5.0]

    compute_weights(u4, x, wdiag, w∂k)
    @test size(wdiag,1)== 3+3
    @test size(w∂k,1)== 3+2

    @test norm(wdiag - [1.0;
                        ψ₀(-1.0, 1.0)(x);
                         ψj(2.0, 2.0)(x);
                        ψj(-3.0, 3.0)(x);
                         ψj(4.0, 4.0)(x);
                      ψpp1(-5.0, 5.0)(x)])<1e-10
end

@testset "Weights of a RadialMapComponent map" begin
    M = RadialMap(1, 0)
    x = randn()
    y = randn(3)

    woff, wdiag, w∂k = create_weights(M)
    # wdiag = zeros(2)
    # w∂k = zeros(1)
    # woff = zeros(0)

    # w, w∂k = compute_weights(TransportMap.component(C.U[1],1), [2.0], w, w∂k)
    compute_weights(TransportBasedInference.component(M.C[1],1), x, wdiag, w∂k)

    @test size(wdiag,1) == 2
    @test size(w∂k, 1) == 1
    @test size(woff, 1) == 0

    @test wdiag == [1.0; x]
    @test w∂k   == [1.0]

    #
    M = RadialMap(1, 1)

    woff, wdiag, w∂k = create_weights(M)

    compute_weights(M, [x], woff,  wdiag, w∂k)
    #
    @test size(woff,1)==0
    @test size(wdiag,1) == 1+3
    @test size(w∂k, 1) == 1+2
    #
    @test woff == Float64[]
    @test isapprox(wdiag, [1.0; ψ₀(0.0, 1.0)(x); ψj(0.0, 1.0)(x); ψpp1(0.0, 1.0)(x)], atol = 1e-10)
    @test isapprox(w∂k, [ψ₀′(0.0, 1.0)(x); rbf(0.0, 1.0)(x); ψpp1′(0.0, 1.0)(x)], atol = 1e-10)
    #
    #
    M = RadialMap(2, 0)

    woff, wdiag, w∂k = create_weights(M)

    #
    compute_weights(M, y[1:2], woff, wdiag, w∂k)
    #
    @test size(woff, 1) == 1
    @test size(wdiag,1) == 2+2
    @test size(w∂k, 1) == 1+1
    #
    #
    @test isapprox(wdiag, [1.0; y[1]; 1.0; y[2]], atol = 1e-10)
    @test isapprox(w∂k, [1.0; 1.0], atol = 1e-10)
    @test isapprox(woff, [y[1]], atol = 1e-10)
    #
    #
    M = RadialMap(3, 2)

    woff, wdiag, w∂k = create_weights(M)

    compute_weights(M, y, woff, wdiag, w∂k)
    #
    @test size(woff,1) == 2*(2+1)
    @test size(wdiag,1) == 3*(2+3)
    @test size(w∂k, 1) == 3*(2+2)
    #
    #
    @test isapprox(woff, [y[1]; rbf(0.0, 1.0)(y[1]); rbf(0.0, 1.0)(y[1]);
                          y[2]; rbf(0.0, 1.0)(y[2]); rbf(0.0, 1.0)(y[2])], atol = 1e-10)

    @test isapprox(wdiag, [1.0; ψ₀(0.0, 1.0)(y[1]); ψj(0.0, 1.0)(y[1]);
                           ψj(0.0, 1.0)(y[1]); ψpp1(0.0, 1.0)(y[1]);
                           1.0; ψ₀(0.0, 1.0)(y[2]); ψj(0.0, 1.0)(y[2]); ψj(0.0, 1.0)(y[2]);ψpp1(0.0, 1.0)(y[2]);
                           1.0; ψ₀(0.0, 1.0)(y[3]); ψj(0.0, 1.0)(y[3]); ψj(0.0, 1.0)(y[3]); ψpp1(0.0, 1.0)(y[3])],
                           atol = 1e-10)

    @test isapprox(w∂k, [ψ₀′(0.0, 1.0)(y[1]); rbf(0.0, 1.0)(y[1]) ; rbf(0.0, 1.0)(y[1]); ψpp1′(0.0, 1.0)(y[1]);
                        ψ₀′(0.0, 1.0)(y[2]); rbf(0.0, 1.0)(y[2]) ; rbf(0.0, 1.0)(y[2]); ψpp1′(0.0, 1.0)(y[2]);
                        ψ₀′(0.0, 1.0)(y[3]); rbf(0.0, 1.0)(y[3]) ; rbf(0.0, 1.0)(y[3]); ψpp1′(0.0, 1.0)(y[3])],
                        atol = 1e-10)
    Nx = 5
    Ne = 50
    p = 2

    X = randn(Nx, Ne) .* randn(Nx, Ne)
    M = RadialMap(Nx, p)

    W = create_weights(M, X)

    woff  = zero(W.woff)
    wdiag = zero(W.wdiag)
    w∂k   = zero(W.w∂k)
    compute_weights(M, X, woff, wdiag, w∂k)

    compute_weights(M, X, W)

    @test isapprox(woff, W.woff, atol = 1e-10)
    @test isapprox(wdiag, W.wdiag, atol = 1e-10)
    @test isapprox(w∂k, W.w∂k, atol = 1e-10)

    woval, wval, w∂kval = create_weights(M)
    compute_weights(M, zeros(Nx), woval, wval, w∂kval)
    compute_weights(M, zeros(Nx,Ne), woff, wdiag, w∂k)

    @test woff[1,:] == woval
    @test wdiag[1,:] == wval
    @test w∂k[1,:] == w∂kval
end

@testset "Weights SparseRadialMapComponentMap I" begin
    x = randn()

    ## k = 1 and p = -1
    C = SparseRadialMapComponent(1, [-1])

    no, nd, n∂ = ncoeff(C.Nx, C.p)

    @test no ==0
    @test nd ==0
    @test n∂ ==0

    wo, wd, w∂ = compute_weights(C, [x])
    @test wo == Float64[]
    @test wd == Float64[]
    @test w∂ == Float64[]

    ## Nx = 1 and p = 0
    C = SparseRadialMapComponent(1, [0])

    no, nd, n∂ = ncoeff(C.Nx, C.p)

    @test no ==0
    @test nd ==1
    @test n∂ ==1

    wo, wd, w∂ = compute_weights(C, [x])
    @test wo == Float64[]
    @test isapprox(wd, [x], atol = 1e-10)
    @test w∂ == [1.0]

    z = randn(1,5)
    wo, wd, w∂ = compute_weights(C, z)
    @test wo == zeros(0,0)
    @test isapprox(wd, z', atol = 1e-10)
    @test w∂ == ones(5,1)

    ## Nx = 1 and p = 2
    C = SparseRadialMapComponent(1, [2])

    no, nd, n∂ = ncoeff(C.Nx, C.p)

    @test no ==0
    @test nd ==2+2
    @test n∂ ==2+2

    wo, wd, w∂ = compute_weights(C, [x])
    @test wo == Float64[]
    @test isapprox(wd, [ψ₀(0.0, 1.0)(x); ψj(0.0, 1.0)(x);ψj(0.0, 1.0)(x);ψpp1(0.0, 1.0)(x)], atol = 1e-10)
    @test isapprox(w∂, [ψ₀′(0.0, 1.0)(x); rbf(0.0, 1.0)(x);rbf(0.0, 1.0)(x);ψpp1′(0.0, 1.0)(x)], atol = 1e-10)

    z = randn(1,5)

    wo, wd, w∂ = compute_weights(C, z)
    @test wo == zeros(0,0)

    for i=1:5
        zi  = z[1,i]
        @test wd[i,:] == [ψ₀(0.0, 1.0)(zi); ψj(0.0, 1.0)(zi);ψj(0.0, 1.0)(zi);ψpp1(0.0, 1.0)(zi)]
        @test w∂[i,:] == [ψ₀′(0.0, 1.0)(zi); rbf(0.0, 1.0)(zi);rbf(0.0, 1.0)(zi);ψpp1′(0.0, 1.0)(zi)]
    end

end

@testset "Weights SparseRadialMapComponent II" begin
    ## k=2 and p = [-1 -1]
    y = randn(2)
    C = SparseRadialMapComponent(2, [-1; -1])

    no, nd, n∂ = ncoeff(C.Nx, C.p)

    @test no ==0
    @test nd ==0
    @test n∂ ==0

    wo, wd, w∂ = compute_weights(C, [1.0; 2.0])
    @test wo == Float64[]
    @test wd == Float64[]
    @test w∂ == Float64[]

    wo, wd, w∂ = compute_weights(C, randn(2,5))
    @test wo == zeros(0, 0)
    @test wd == zeros(0, 0)
    @test w∂ == zeros(0, 0)

    ## k = 2 and p = [-1 0]
    C = SparseRadialMapComponent(2, [-1; 0])
    no, nd, n∂ = ncoeff(C.Nx, C.p)

    @test no ==0
    @test nd ==1
    @test n∂ ==1

    wo, wd, w∂ = compute_weights(C, y)
    @test wo == Float64[]
    @test wd == [y[2]]
    @test w∂ == [1.0]

    z = randn(2, 5)
    wo, wd, w∂ = compute_weights(C, deepcopy(z))
    @test wo == zeros(0,0)
    @test isapprox(wd[:,1], z[2,:], atol = 1e-10)
    @test isapprox(w∂[:,1], ones(5), atol = 1e-10)

    ## k = 2 and p = [-1 1]
    C = SparseRadialMapComponent(2, [-1; 1])
    no, nd, n∂ = ncoeff(C.Nx, C.p)

    @test no ==0
    @test nd ==3
    @test n∂ ==3

    wo, wd, w∂ = compute_weights(C, y)
    @test wo == Float64[]
    @test isapprox(wd, [ψ₀(0.0, 1.0)(y[2]); ψj(0.0, 1.0)(y[2]);ψpp1(0.0, 1.0)(y[2])], atol = 1e-10)
    @test isapprox(w∂, [ψ₀′(0.0, 1.0)(y[2]); rbf(0.0, 1.0)(y[2]);ψpp1′(0.0, 1.0)(y[2])], atol = 1e-10)

    z = randn(2,5)

    wo, wd, w∂ = compute_weights(C, z)
    @test wo == zeros(0,0)

    for i=1:5
        zi  = z[:,i]
        @test wd[i,:] == [ψ₀(0.0, 1.0)(zi[2]); ψj(0.0, 1.0)(zi[2]);ψpp1(0.0, 1.0)(zi[2])]
        @test w∂[i,:] == [ψ₀′(0.0, 1.0)(zi[2]); rbf(0.0, 1.0)(zi[2]);ψpp1′(0.0, 1.0)(zi[2])]
    end

    ## k = 2 and p = [0 1]

    C = SparseRadialMapComponent(2, [0; 1])
    no, nd, n∂ = ncoeff(C.Nx, C.p)

    @test no ==1
    @test nd ==3
    @test n∂ ==3

    wo, wd, w∂ = compute_weights(C, y)
    @test isapprox(wo, [y[1]], atol = 1e-10)
    @test isapprox(wd, [ψ₀(0.0, 1.0)(y[2]); ψj(0.0, 1.0)(y[2]);ψpp1(0.0, 1.0)(y[2])], atol = 1e-10)
    @test isapprox(w∂, [ψ₀′(0.0, 1.0)(y[2]); rbf(0.0, 1.0)(y[2]);ψpp1′(0.0, 1.0)(y[2])], atol = 1e-10)

    z = randn(2,5)

    wo, wd, w∂ = compute_weights(C, deepcopy(z))

    for i=1:5
        zi  = z[:,i]
        @test norm(wo[i,:] - [zi[1]])<1e-10

        @test norm(wd[i,:]  - [ψ₀(0.0, 1.0)(zi[2]); ψj(0.0, 1.0)(zi[2]);ψpp1(0.0, 1.0)(zi[2])])<1e-10
        @test norm(w∂[i,:]  - [ψ₀′(0.0, 1.0)(zi[2]); rbf(0.0, 1.0)(zi[2]);ψpp1′(0.0, 1.0)(zi[2])])<1e-10
    end

    # k = 2 and p = [2 1]

    C = SparseRadialMapComponent(2, [2; 1])
    C.ξ[1] .= randn(2)
    C.ξ[2] .= randn(3)
    C.σ[1] .= rand(2)
    C.σ[2] .= rand(3)

    no, nd, n∂ = ncoeff(C.Nx, C.p)

    @test no ==3
    @test nd ==3
    @test n∂ ==3

    wo, wd, w∂ = compute_weights(C, y)
    @test isapprox(wo, [y[1]; rbf(C.ξ[1][1], C.σ[1][1])(y[1]); rbf(C.ξ[1][2], C.σ[1][2])(y[1])], atol = 1e-10)
    @test isapprox(wd, [ψ₀(C.ξ[2][1], C.σ[2][1])(y[2]); ψj(C.ξ[2][2], C.σ[2][2])(y[2]);
                 ψpp1(C.ξ[2][3], C.σ[2][3])(y[2])], atol = 1e-10)
    @test isapprox(w∂, [ψ₀′(C.ξ[2][1], C.σ[2][1])(y[2]); rbf(C.ξ[2][2], C.σ[2][2])(y[2]);
                       ψpp1′(C.ξ[2][3], C.σ[2][3])(y[2])], atol = 1e-10)

    z = randn(2,5)

    wo, wd, w∂ = compute_weights(C, deepcopy(z))

    for i=1:5
        zi  = z[:,i]
        @test norm(wo[i,:]  - [zi[1]; rbf(C.ξ[1][1], C.σ[1][1])(zi[1]); rbf(C.ξ[1][2], C.σ[1][2])(zi[1])])<1e-14

        @test norm(wd[i,:]  - [ψ₀(C.ξ[2][1], C.σ[2][1])(zi[2]); ψj(C.ξ[2][2], C.σ[2][2])(zi[2]);ψpp1(C.ξ[2][3], C.σ[2][3])(zi[2])])<1e-14
        @test norm(w∂[i,:]  - [ψ₀′(C.ξ[2][1], C.σ[2][1])(zi[2]); rbf(C.ξ[2][2], C.σ[2][2])(zi[2]);ψpp1′(C.ξ[2][3], C.σ[2][3])(zi[2])])<1e-14
    end

    ## k = 4 p = [-1 2 0 2]

    C = SparseRadialMapComponent(4, [-1; 2; 0; 2])

    no, nd, n∂ = ncoeff(C.Nx, C.p)

    @test no ==3+1
    @test nd ==4
    @test n∂ ==4

    y = randn(4)

    wo, wd, w∂ = compute_weights(C, deepcopy(y))

    @test norm(wo - [y[2]; rbf(0.0,1.0)(y[2]); rbf(0.0,1.0)(y[2]); y[3]])<1e-14
    @test norm(wd - [ψ₀(0.0, 1.0)(y[4]); ψj(0.0, 1.0)(y[4]); ψj(0.0, 1.0)(y[4]);ψpp1(0.0, 1.0)(y[4])])<1e-14
    @test norm(w∂ - [ψ₀′(0.0, 1.0)(y[4]); rbf(0.0, 1.0)(y[4]); rbf(0.0, 1.0)(y[4]);ψpp1′(0.0, 1.0)(y[4])])<1e-14

    z = randn(4, 5)

    wo, wd, w∂ = compute_weights(C, deepcopy(z))
    for i=1:5
        zi  = z[:,i]
        @test norm(wo[i,:]  - [zi[2]; rbf(0.0,1.0)(zi[2]); rbf(0.0,1.0)(zi[2]); zi[3]])<1e-14

        @test norm(wd[i,:]  - [ψ₀(0.0, 1.0)(zi[4]); ψj(0.0, 1.0)(zi[4]); ψj(0.0, 1.0)(zi[4]);ψpp1(0.0, 1.0)(zi[4])])<1e-14
        @test norm(w∂[i,:]  - [ψ₀′(0.0, 1.0)(zi[4]); rbf(0.0, 1.0)(zi[4]); rbf(0.0, 1.0)(zi[4]);ψpp1′(0.0, 1.0)(zi[4])])<1e-14
    end

    ## k = 5 and p = [2 -1 0 -1 0]

    C = SparseRadialMapComponent(5, [2; -1; 0; -1; 0])

    no, nd, n∂ = ncoeff(C.Nx, C.p)

    @test no ==3+1
    @test nd ==1
    @test n∂ ==1

    y = randn(5)
    wo, wd, w∂ = compute_weights(C, y)

    @test norm(wo - [y[1]; rbf(0.0,1.0)(y[1]); rbf(0.0,1.0)(y[1]); y[3]])<1e-14
    @test norm(wd - [y[5]])<1e-14
    @test norm(w∂ - [1.0])<1e-14

    z = randn(5, 10)

    wo, wd, w∂ = compute_weights(C, z)
    for i=1:10
        zi  = z[:,i]
        @test norm(wo[i,:]  - [zi[1]; rbf(0.0,1.0)(zi[1]); rbf(0.0,1.0)(zi[1]); zi[3]])<1e-14

        @test norm(wd[i,:]  - [zi[5]])<1e-14
        @test norm(w∂[i,:]  - [1.0])<1e-14
    end

end


text:
# This file is a part of TypeDBClient.  License is MIT: https://github.com/Humans-of-Julia/TypeDBClient.jl/blob/main/LICENSE

#
# package typedb.client.api.query;
#
# import typedb.client.api.TypeDBOptions;
# import typedb.client.api.answer.ConceptMap;
# import typedb.client.api.answer.ConceptMapGroup;
# import typedb.client.api.answer.Numeric;
# import typedb.client.api.answer.NumericGroup;
# import graql.lang.query.GraqlDefine;
# import graql.lang.query.GraqlDelete;
# import graql.lang.query.GraqlInsert;
# import graql.lang.query.GraqlMatch;
# import graql.lang.query.GraqlUndefine;
# import graql.lang.query.GraqlUpdate;
#
# import javax.annotation.CheckReturnValue;
# import java.util.stream.Stream;
#
# public interface QueryManager {
#
#     @CheckReturnValue
#     Stream<ConceptMap> match(GraqlMatch query);
#
#     @CheckReturnValue
#     Stream<ConceptMap> match(GraqlMatch query, TypeDBOptions options);
#
#     @CheckReturnValue
#     QueryFuture<Numeric> match(GraqlMatch.Aggregate query);
#
#     @CheckReturnValue
#     QueryFuture<Numeric> match(GraqlMatch.Aggregate query, TypeDBOptions options);
#
#     @CheckReturnValue
#     Stream<ConceptMapGroup> match(GraqlMatch.Group query);
#
#     @CheckReturnValue
#     Stream<ConceptMapGroup> match(GraqlMatch.Group query, TypeDBOptions options);
#
#     @CheckReturnValue
#     Stream<NumericGroup> match(GraqlMatch.Group.Aggregate query);
#
#     @CheckReturnValue
#     Stream<NumericGroup> match(GraqlMatch.Group.Aggregate query, TypeDBOptions options);
#
#     Stream<ConceptMap> insert(GraqlInsert query);
#
#     Stream<ConceptMap> insert(GraqlInsert query, TypeDBOptions options);
#
#     QueryFuture<Void> delete(GraqlDelete query);
#
#     QueryFuture<Void> delete(GraqlDelete query, TypeDBOptions options);
#
#     Stream<ConceptMap> update(GraqlUpdate query);
#
#     Stream<ConceptMap> update(GraqlUpdate query, TypeDBOptions options);
#
#     QueryFuture<Void> define(GraqlDefine query);
#
#     QueryFuture<Void> define(GraqlDefine query, TypeDBOptions options);
#
#     QueryFuture<Void> undefine(GraqlUndefine query);
#
#     QueryFuture<Void> undefine(GraqlUndefine query, TypeDBOptions options);
# }


text:
using RNGPool
using CoupledConditionalSMC
import Statistics: mean, var

include("lgModel.jl")

function runDemo(n::Int64, N::Int64, b::Int64, m::Int64, maxit::Int64)
  model, lM, ko = setupLGModel(n, 1.0)

  function h(path::Vector{Float64Particle})
    return path[1].x
  end

  println("\n-----\nTrue value = ", ko.smoothingMeans[1])

  println("\nBackward sampling, dependent initialization")
  resultsBS = CoupledConditionalSMC.unbiasedEstimates(model, lM, h, N, n,
    b, m, :BS, false, true, maxit)
  println("mean of estimates = ", mean(resultsBS[2]))
  println("mean no. iterations = ", mean(resultsBS[1]))
  println("estimated variance of estimator = ", var(resultsBS[2]))
  println("should be approx. standard normal: ",
    (mean(resultsBS[2]) - ko.smoothingMeans[1])/sqrt(var(resultsBS[2])/m))

  println("\nBackward sampling, independent initialization")
  resultsBS = CoupledConditionalSMC.unbiasedEstimates(model, lM, h, N, n,
    b, m, :BS, true, true, maxit)
  println("mean of estimates = ", mean(resultsBS[2]))
  println("mean no. iterations = ", mean(resultsBS[1]))
  println("estimated variance of estimator = ", var(resultsBS[2]))
  println("should be approx. standard normal: ",
    (mean(resultsBS[2]) - ko.smoothingMeans[1])/sqrt(var(resultsBS[2])/m))

  println("\nAncestral tracing, dependent initialization")
  results = CoupledConditionalSMC.unbiasedEstimates(model, h, N, n, b, m,
    false, true, maxit)
  println("mean of estimates = ", mean(results[2]))
  println("mean no. iterations = ", mean(results[1]))
  println("estimated variance of estimator = ", var(results[2]))
  println("should be approx. standard normal: ",
    (mean(results[2]) - ko.smoothingMeans[1])/sqrt(var(results[2])/m))

  println("\nAncestral tracing, independent initialization")
  results = CoupledConditionalSMC.unbiasedEstimates(model, h, N, n, b, m,
    true, true, maxit)
  println("mean of estimates = ", mean(results[2]))
  println("mean no. iterations = ", mean(results[1]))
  println("estimated variance of estimator = ", var(results[2]))
  println("should be approx. standard normal: ",
    (mean(results[2]) - ko.smoothingMeans[1])/sqrt(var(results[2])/m))

  println("\nAncestor sampling, dependent initialization")
  results = CoupledConditionalSMC.unbiasedEstimates(model, lM, h, N, n, b, m,
    :AS, false, true, maxit)
  println("mean of estimates = ", mean(results[2]))
  println("mean no. iterations = ", mean(results[1]))
  println("estimated variance of estimator = ", var(results[2]))
  println("should be approx. standard normal: ",
    (mean(results[2]) - ko.smoothingMeans[1])/sqrt(var(results[2])/m))

  println("\nAncestor sampling, independent initialization")
  results = CoupledConditionalSMC.unbiasedEstimates(model, lM, h, N, n, b, m,
    :AS, true, true, maxit)
  println("mean of estimates = ", mean(results[2]))
  println("mean no. iterations = ", mean(results[1]))
  println("estimated variance of estimator = ", var(results[2]))
  println("should be approx. standard normal: ",
    (mean(results[2]) - ko.smoothingMeans[1])/sqrt(var(results[2])/m))
  println("-----\n")
end

setRNGs(12345)

runDemo(10, 16, 1, 1000, 100000)
runDemo(20, 16, 1, 1000, 100000)
runDemo(30, 32, 1, 1000, 100000)