text:
using CategoricalArrays

struct OrderedRuleSet <: ParticleClassifier
    name::String
    rules::Vector{Tuple{String,Function}}

    OrderedRuleSet(name::String, rules::Tuple{String,Function}...) = new(name, collect(rules))

    function OrderedRuleSet(orss::OrderedRuleSet...)
        rules = Vector{Tuple{String,Function}}()
        name = String[]
        for ors in orss
            append!(rules, ors.rules)
            push!(name, ors.name)
        end
        return new(join(name, "+"), rules)
    end
end

Base.show(io::IO, ors::OrderedRuleSet) = print(io, "JORS[$(ors.name)]")
classnames(sr::OrderedRuleSet) = collect(map(r -> r[1], sr.rules))

"""
    classify(zep::Zeppelin, catcls::CategoricalArray{String})

Assign the CLASS column to the contents of this CategoricalArray.
"""
function classify(zep::Zeppelin, catcls::CategoricalArray{String})::Zeppelin
    @assert length(catcls) == nrow(zep.data)
    # Figure out where to put the column...
    cols = collect(1:ncol(zep.data))
    clscol = findfirst(n -> n == :CLASS, propertynames(zep.data))
    if !isnothing(clscol)
        cols = deleteat!(cols, clscol)
    else
        tmp = findfirst(n -> n == :FIRSTELM, propertynames(zep.data))
        clscol = isnothing(tmp) ? ncol(zep.data) + 1 : tmp
    end
    # Copy the data and insert the column
    result = Zeppelin(zep.headerfile, copy(zep.header), copy(zep.data[:, cols]))
    insertcols!(result.data, clscol, :CLASS => catcls)
    return result
end

"""
    classify(zep::Zeppelin, )::Zeppelin

Classify the rows in `zep`, row-by-row by calling the applying the `OrderedRuleSet` to each row. Returns a new
Zeppelin object.
"""
function classify(zep::Zeppelin, sr::OrderedRuleSet)::Zeppelin
    cat = categorical(append!(classnames(sr), ["Other", "Missing", "Error"]), compress=false, ordered = true) # set up the levels
    cat = cat[1:0] # delete the data but retain the levels
    errcx = 0
    missingCols = Dict(s => 0.0 for s in filter(
        sy -> !(sy in propertynames(zep.data)),
        map(z -> convert(Symbol,elements[z]), 1:95),
    ))
    for row in eachparticle(zep)
        try
            input = Dict(n => zep.data[row, n] for n in propertynames(zep.data))
#            if all(map(v->!ismissing(v),values(input)))
                merge!(input, missingCols)
                push!(cat, classify(sr, input))
#            else
#                push!(cat, "Missing")
#            end
        catch err
            push!(cat, "Other")
            if (errcx += 1) <= 4
                @info "ERROR[$errcx]: " exception=(err, catch_backtrace())
            end
        end
    end
    res = classify(zep, cat)
    res.header["RULE_FILE"]=repr(sr)
    return res
end

function classify(sr::OrderedRuleSet, input::Dict{Symbol,Any})
    i = findfirst(rule -> rule[2](input), sr.rules)
    return isnothing(i) ? "Other" : sr.rules[i][1]
end

const NullRules = OrderedRuleSet("Null", ("Unclassified", inp -> true) )

const BaseRules = OrderedRuleSet(
    "Base",
    ("LowCounts", inp -> haskey(inp,:COUNTS) && inp[:COUNTS] < 2000),
    (
     "Maraging",
     inp -> (inp[:FIRSTELM] == n"Fe") &&
            (inp[:NI] > 5) && (inp[:CO] > 2) && (inp[:FE] + inp[:NI] + inp[:CO] + inp[:MO] > 80.0),
    ),
    ("Quartz-like", inp -> inp[:SI] > 80),
    (
     "Biotite",
     inp -> (inp[:K] > 5) &&
            (inp[:MG] + inp[:FE] > 20) &&
            (inp[:AL] > 10) &&
            (inp[:SI] > 10) && (inp[:K] + inp[:MG] + inp[:FE] + inp[:AL] + inp[:SI] > 80) && (inp[:O] > 5),
    ),
    (
     "Plagioclaise",
     inp -> (inp[:CA] + inp[:NA] > 20) &&
            (inp[:AL] > 10) && (inp[:SI] > 10) && (inp[:NA] + inp[:CA] + inp[:AL] + inp[:SI] > 80) && (inp[:O] > 5),
    ),
    (
     "Othoclaise",
     inp -> (inp[:K] > 5) &&
            (inp[:AL] > 20) && (inp[:SI] > 20) && (inp[:O] > 10) && (inp[:K] + inp[:AL] + inp[:SI] > 80),
    ),
    ("Dolomite", inp -> (inp[:CA] > 40) && (inp[:MG] > 10) && (inp[:CA] + inp[:MG] > 70)),
    ("Olivine", inp -> (inp[:MG] + inp[:FE] > 40) && (inp[:SI] > 20) && (inp[:MG] + inp[:FE] + inp[:SI] > 80)),
    (
     "Feldspar",
     inp -> (inp[:K] + inp[:NA] + inp[:CA] > 20) &&
            (inp[:AL] > 10) && (inp[:SI] > 20) && (inp[:K] + inp[:NA] + inp[:CA] + inp[:AL] + inp[:SI] > 80),
    ),
    ("Willemite", inp -> (inp[:ZN] > 20) && (inp[:SI] > 20) && (inp[:ZN] + inp[:SI] > 80)),
    (
     "Garnet",
     inp -> (inp[:MG] + inp[:FE] + inp[:MN] + inp[:CA] > 10) &&
            (inp[:AL] > 5) && (inp[:SI] > 10) && (inp[:MG] + inp[:FE] + inp[:MN] + inp[:CA] + inp[:AL] + inp[:SI] > 80),
    ),
    (
     "Al-Cu alloy (2XXX)",
     inp -> (inp[:AL] > 60) && (inp[:SECONDELM] == n"Cu") && (inp[:CU] > 3) && (inp[:AL] + inp[:CU] > 90),
    ),
    (
     "Al-Mn alloy (3XXX)",
     inp -> (inp[:AL] > 60) && (inp[:SECONDELM] == n"Mn") && (inp[:MN] > 3) && (inp[:AL] + inp[:MN] > 90),
    ),
    (
     "Al-Si alloy (4XXX)",
     inp -> (inp[:AL] > 80) && (inp[:SECONDELM] == n"Si") && (inp[:SI] > 1) && (inp[:AL] + inp[:SI] > 90),
    ),
    (
     "Al-Mg alloy (5XXX)",
     inp -> (inp[:FIRSTELM] == n"Al") && (inp[:SECONDELM] == n"Mg") && (inp[:MG] > 5) && (inp[:AL] + inp[:MG] > 80),
    ),
    (
     "Al-Si-Mg (6XXX)",
     inp -> (inp[:AL] > 90) &&
            (inp[:SI] > 0) && (inp[:MG] > 0) && (inp[:AL] + inp[:SI] + inp[:MG] + inp[:MN] + inp[:CU] + inp[:FE] > 98),
    ),
    (
     "Al-Zn alloy (7XXX)",
     inp -> (inp[:AL] > 80) && (inp[:SECONDELM] == n"Zn") && (inp[:ZN] > 2) && (inp[:AL] + inp[:ZN] > 90),
    ),
    ("Al-Fe alloy", inp -> (inp[:AL] > 80) && (inp[:SECONDELM] == n"Fe") && (inp[:FE] > 2) && (inp[:AL] + inp[:FE] > 90)),
    ("Aluminosilicate", inp -> (inp[:AL] > 30) && (inp[:SI] > 10) && (inp[:AL] + inp[:SI] > 80)),
    ("Zircon", inp -> (inp[:ZR] > 20) && (inp[:ZR] + inp[:SI] > 80)),
    (
     "Elgiloy",
     inp -> (inp[:CO] > 20) &&
            (inp[:CR] > 10) &&
            (inp[:NI] > 5) &&
            (inp[:FE] > 5) &&
            ((inp[:MO] > 2) || (inp[:S] > 2)) &&
            (inp[:CO] + inp[:CR] + inp[:NI] + inp[:FE] + inp[:MO] + inp[:MN] + inp[:S] > 80),
    ),
    ("Aluminum", inp -> inp[:AL] > 80),
    ("Fe-Al alloy", inp -> (inp[:AL] > 20) && (inp[:FE] > 20) && (inp[:AL] + inp[:FE] > 80)),
    ("Calcium silicate", inp -> (inp[:CA] > 40) && (inp[:SI] > 10) && (inp[:CA] + inp[:SI] > 80)),
    ("Fe-Cr stainless", inp -> (inp[:FE] > 50) && (inp[:CR] > 5) && (inp[:NI] < 3) && (inp[:FE] + inp[:CR] > 90)),
    (
     "Fe-Cr-Ni Stainless",
     inp -> (inp[:FIRSTELM] == n"Fe") && (inp[:CR] > 5) && (inp[:NI] > 5) && (inp[:FE] + inp[:CR] + inp[:NI] > 80),
    ),
    ("Fe-Ni stainless", inp -> (inp[:FE] > 50) && (inp[:NI] > 10) && (inp[:FE] + inp[:NI] > 90)),
    ("Iron oxide", inp -> (inp[:FE] > 80) && (inp[:O] > 10)),
    ("Iron", inp -> inp[:FE] > 80),
    ("Calcite", inp -> (inp[:CA] > 80) && (inp[:O] > 5)),
    ("Zinc", inp -> inp[:ZN] > 80),
    ("Silicate", inp -> inp[:FIRSTELM] == n"Si" ),
    (
     "Anglesite",
     inp -> (inp[:FIRSTELM] == n"Pb") && (inp[:SECONDELM] == n"S") && (inp[:PB] + inp[:S] > 80) && (inp[:O] > 5),
    ),
    ("Hashemite", inp -> (inp[:BA] > 20) && (inp[:CR] > 20) && (inp[:BA] + inp[:CR] > 80) && (inp[:O] > 5)),
    ("Barite", inp -> (inp[:BA] > 40) && (inp[:S] > 10) && (inp[:BA] + inp[:S] > 80) && (inp[:O] > 5)),
    ("Gypsum", inp -> (inp[:CA] > 20) && (inp[:S] > 20) && (inp[:CA] + inp[:S] > 80) && (inp[:O] > 5)),
    ("Lead", inp -> (inp[:PB] > 60) && (inp[:PB] + inp[:S] > 80)),
    ("Salt", inp -> (inp[:K] + inp[:NA] > 20) && (inp[:CL] > 20) && (inp[:K] + inp[:NA] + inp[:CL] > 80)),
    ("Bismuth", inp -> (inp[:BI] > 80)),
    ("Brass", inp -> (inp[:CU] > 20) && (inp[:ZN] > 20) && (inp[:CU] + inp[:ZN] > 80)),
    ("Bronze", inp -> (inp[:CU] > 20) && (inp[:SN] > 20) && (inp[:CU] + inp[:SN] > 80)),
    ("Cupronickel", inp -> (inp[:CU] > 30) && (inp[:NI] > 10) && (inp[:CU] + inp[:NI] > 80)),
    ("Sodium sulfide", inp -> (inp[:NA] > 20) && (inp[:S] > 10) && (inp[:NA] + inp[:S] > 80)),
    ("Wurtzite (ZnS)", inp -> (inp[:ZN] > 20) && (inp[:S] > 10) && (inp[:ZN] + inp[:S] > 80)),
    ("Pyrite (FeS2)", inp -> (inp[:FE] > 20) && (inp[:S] > 20) && (inp[:FE] + inp[:S] > 80)),
    ("Silver sulfide", inp -> (inp[:AG] > 20) && (inp[:S] > 10) && (inp[:AG] + inp[:S] > 80)),
    (
     "Sn62 solder",
     inp -> (inp[:SN] > 30) && (inp[:PB] > 10) && (inp[:AG] > 1) && (inp[:SN] + inp[:PB] + inp[:AG] > 80),
    ),
    (
     "SACZ solder",
     inp -> (inp[:SN] + inp[:AG] + inp[:CU] + inp[:ZN] > 80) &&
            (inp[:SN] > 10) && (inp[:CU] > 10) && (inp[:AG] > 10) && (inp[:ZN] > 2),
    ),
    ("Lead solder", inp -> (inp[:PB] > 20) && (inp[:SN] > 20) && (inp[:PB] + inp[:SN] > 80)),
    (
     "SAC solder",
     inp -> (inp[:SN] + inp[:AG] + inp[:CU] > 80) && (inp[:SN] > 10) && (inp[:CU] > 10) && (inp[:AG] > 10),
    ),
    ("Titanium", inp -> (inp[:TI] > 80)),
    ("Calcium phosphate", inp -> (inp[:CA] > 40) && (inp[:P] > 10) && (inp[:CA] + inp[:P] > 80)),
    ("Calcium fluoride", inp -> (inp[:CA] > 40) && (inp[:F] > 10) && (inp[:CA] + inp[:F] > 80)),
    ("Iron-bearing", inp -> (inp[:FE] > 50)),
    ("Ca-bearing", inp -> inp[:CA] > 50),
    ("Other Silicate", inp -> (inp[:SI] > 20) && (inp[:AL] + inp[:FE] + inp[:CA] + inp[:SI] > 80)),
    ("Gold-Silver alloy", inp -> (inp[:FIRSTELM] == n"Au") && (inp[:SECONDELM] == n"Ag") && (inp[:AU] + inp[:AG] > 80)),
    ("Gold - other", inp -> inp[:AU] > 80),
    ("Tin", inp -> inp[:SN] > 80),
    ("Celestine", inp -> (inp[:SR] > 40) && (inp[:S] > 10) && (inp[:SR] + inp[:S] > 80) && (inp[:O] > 5)),
    ("Barite-Celestine", inp -> (inp[:BA] + inp[:SR] > 40) && (inp[:BA] + inp[:SR] + inp[:S] > 80) && (inp[:O] > 5)),
    ("Cu-rich", inp -> inp[:CU] > 80),
    ("Cr-bearing", inp -> inp[:CR] > 50),
    ("Silver sulfide", inp -> (inp[:AG] > 40) && (inp[:S] > 10) && (inp[:AG] + inp[:S] > 80)),
    ("Silver", inp -> inp[:AG] > 80),
    ("MoS2", inp -> inp[:MO] > 40 && inp[:S] > 30 && (inp[:MO] + inp[:S] > 90)),
    (
     "Monazite",
     inp -> (inp[:LA] + inp[:CE] + inp[:ND] + inp[:TH] + inp[:Y] > 30) &&
            (inp[:P] > 10) && (inp[:LA] + inp[:CE] + inp[:ND] + inp[:TH] + inp[:Y] + inp[:P] + inp[:SI] > 80),
    ),
    ("Lanthanide", inp -> (inp[:LA] + inp[:CE] + inp[:ND] + inp[:Y] > 60)),
    ("Silver chloride", inp -> (inp[:AG] > 40) && (inp[:CL] > 10)),
    ("Magnesium oxide", inp -> (inp[:MG] > 80) && (inp[:O] > 10)),
    ("Nickel", inp -> inp[:NI] > 80),
    ("Ca-Fe silicate", inp -> (inp[:CA] > 20) && (inp[:FE] > 20) && (inp[:SI] > 5)),
    ("Ca-Ti silicate", inp -> (inp[:CA] > 20) && (inp[:TI] > 20) && (inp[:SI] > 10)),
    ("Al-Zr-Cl", inp -> (inp[:AL] > 20) && (inp[:ZR] > 20) && (inp[:CL] > 5)),
    ("Chlorine-rich", inp -> inp[:CL] > 80),
    ("Au-Cu alloy", inp -> (inp[:AU] > 40) && (inp[:CU] > 10) && (inp[:AU] + inp[:CU] + inp[:NI] > 80)),
    ("Uranium-bearing", inp -> inp[:U] > 20),
    ("Thorium-bearing", inp -> inp[:TH] > 20),
    ("Platinum", inp -> inp[:PT] > 80)
);


text:
function str2et( dates::Array{ASCIIString,1} )
        ets = Array{Float64}( length(dates) )
    for i = 1:length(dates)
         ets[i] = str2et( dates[i] )
    end
    return ets
end

function cnmfrm( cnames::Array{ASCIIString,1} )
    frcodes = Array{Int32,1}( length(cnames) )
    frnames = Array{String,1}( length(cnames) )
    founds = Array{Bool,1}( length(cnames) )
    for i = 1:length(cnames)
        (frcodes[i], frnames[i], founds[i]) = cnmfrm( cnames[i] )
    end
    return (frcodes, frnames, founds)
end

function bodc2n( codes::Array{Int64,1} )
    names = Array{String,2}( 1, length(codes) )
    founds = Array{Bool,2}( 1, length(codes) )
    for i = 1:length(codes)
        (names[1, i], founds[1, i]) = bodc2n( codes[i] )
    end
    return (names, founds )
end
#

function bodc2s( codes::Array{Int64,1} )
    names = Array{String,1}( length(codes) )
    for i = 1:length(codes)
        names[i] = bodc2s( codes[i] )
    end
    return names
end

function cidfrm( cents::Array{Int64,1} )
    frcodes = Array{Int32,1}( length(cents) )
    frnames = Array{String,1}( length(cents) )
    founds = Array{Bool,1}( length(cents) )
    for i = 1:length(cents)
        (frcodes[i], frnames[i], founds[i]) = cidfrm( cents[i] )
    end
    return (frcodes, frnames, founds)
end

function cnmfrm( cnames::Array{ASCIIString,1} )
    frcodes = Array{Int32,1}( length(cnames) )
    frnames = Array{String,1}( length(cnames) )
    founds = Array{Bool,1}( length(cnames) )
    for i = 1:length(cnames)
        (frcodes[i], frnames[i], founds[i]) = cnmfrm( cnames[i] )
    end
    return (frcodes, frnames, founds)
end

function et2lst(
    ets::Array{Float64}, # SpiceDouble
    body::Int64, # SpiceInt
    lon::Float64, # SpiceDouble
    _type::AbstractString) # Ptr{ConstSpiceChar}
    hrs = Array{Int64,1}( length(ets) )      # Ptr{SpiceInt}
    mns = Array{Int64,1}( length(ets) )   # Ptr{SpiceInt}
    scs = Array{Int64,1}( length(ets) )   # Ptr{SpiceInt}
    times = Array{String,1}( length(ets) )  # Ptr{SpiceChar}
    ampm = Array{String,1}( length(ets) )   # Ptr{SpiceChar}
    for i = 1;length(ets)
        (hrs[i], mns[i], scs[i], times[i], ampm[i]) = et2lst( ets[i], body, lon, _type )
    end
    return (hrs, mns, scs, times, ampms)
end

ext: jl

max_stars_repo_name: lintondf/NasaSpice

max_issues_repo_name: lintondf/NasaSpice

max_forks_repo_name: lintondf/NasaSpice


text:
import Knet: relu, atype
import Sloth: F

@testset "layers" begin
    batch_size, in_features, out_features = 2, 3, 4
    x2d = atype(randn(in_features, batch_size))
    @testset "linear" begin
        import Sloth: Linear
        layer = Linear(in_features=in_features, out_features=out_features)
        @test layer.weight * x2d .+ layer.bias ≈ layer(x2d)
        @test size(layer.weight) == (out_features, in_features)
        @test layer.in_features == in_features
        @test layer.out_features == out_features
        @test size(layer.bias) == (out_features, 1)

        layer = Linear(
            in_features=in_features,
            out_features=out_features,
            bias=false)
        @test size(layer.bias) == ()
        @test layer.bias == F(0.0)
    end

    width, height, in_channels, out_channels, kernel_size = 8, 8, 2, 3, 5
    x4d = atype(rand(width, height, in_channels, batch_size))
    @testset "conv" begin
        import Sloth: Conv, ConvTranspose

        layer = Conv(
            in_channels=in_channels,
            out_channels=out_channels,
            kernel_size=kernel_size)
        @test conv4(layer.weight, x4d) .+ layer.bias ≈ layer(x4d)

        y4d = layer(x4d)
        layer = ConvTranspose(
            in_channels=out_channels,
            out_channels=in_channels,
            kernel_size=kernel_size)
        @test (size(layer(y4d)) == size(x4d))
        @test (layer(y4d) ≈ deconv4(layer.weight ,y4d) .+ layer.bias)
    end

    @testset "pool" begin
        import Sloth: MaxPool, AvgPool, PaddedAvgPool
        layer = MaxPool()
        @test (layer(x4d) ≈ pool(x4d))
        layer = AvgPool()
        @test (layer(x4d) ≈ pool(x4d, mode=2))
        layer = PaddedAvgPool()
        @test (layer(x4d) ≈ pool(x4d, mode=1))
    end

    @testset "batchnorm" begin
        import Sloth: BatchNorm
        layer = BatchNorm(num_features=in_features)
        @test layer(x2d) ≈ batchnorm(x2d, layer.moments, layer.weight)
        layer = BatchNorm(num_features=in_channels)
        @test layer(x4d) ≈ batchnorm(x4d, layer.moments, layer.weight)
    end

    x1d = rand(1:in_features, batch_size)
    @testset "embedding" begin
        import Sloth: Embedding
        layer = Embedding(
            num_embeddings=in_features,
            embedding_dim=out_features)
        @test (layer(x1d) ≈ layer.weight[:, x1d])
    end

    @testset "dropout" begin
        import Sloth: Dropout
        layer = Dropout()
        @test layer(x2d) == x2d
        @test begin
            Knet.seed!(1); y1 = layer(x2d; drop=true)
            Knet.seed!(1); y2 = dropout(x2d, layer.p; drop=true)
            y1 ≈ y2
        end
    end

    @testset "activation" begin
        import Sloth: Relu, Tanh, Sigm, LeakyRelu
        layer = Relu(); @test layer(x2d) ≈ relu.(x2d)
        layer = Tanh(); @test layer(x2d) ≈ tanh.(x2d)
        layer = Sigm(); @test layer(x2d) ≈ sigm.(x2d)
        layer = LeakyRelu(); F = eltype(atype())
        @test layer.α > F(0)
        @test layer(x2d) ≈ max.(F(0), x2d) + layer.α .* min.(F(0), x2d)
    end

    input_size, hidden_size = 3, 5
    @testset "recurrent" begin
        import Sloth: ReluRNN, TanhRNN, LSTM, GRU, BiLSTM
        relurnn = ReluRNN(input_size=input_size, hidden_size=hidden_size)
        tanhrnn = TanhRNN(input_size=input_size, hidden_size=hidden_size)
        lstm = LSTM(input_size=input_size, hidden_size=hidden_size)
        gru = GRU(input_size=input_size, hidden_size=hidden_size)
        bilstm = BiLSTM(input_size=input_size, hidden_size=hidden_size)

        @test relurnn.mode == 0
        @test tanhrnn.mode == 1
        @test lstm.mode == 2
        @test gru.mode == 3
        @test lstm.inputSize == input_size
        @test lstm.hiddenSize == hidden_size
        @test lstm.direction == 0
        @test bilstm.mode == 2
        @test bilstm.direction == 1
    end

    @testset "sequential" begin
        import Sloth: Sequential, Linear
        net = Sequential(
            Linear(in_features=in_features, out_features=out_features),
            Linear(in_features=out_features, out_features=out_features))
        @test net(x2d) ≈ net.layers[2](net.layers[1](x2d))
    end
end


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

# Test functionality of IntSet

## IntSet

# Construction, collect
data_in = (1,5,100)
s = IntSet(data_in)
data_out = collect(s)
@test all(map(d->in(d,data_out), data_in))
@test length(data_out) == length(data_in)

# eltype, similar
@test eltype(IntSet()) === Int64
@test eltype(IntSet) === Int64
@test isequal(similar(IntSet([1,2,3])), IntSet())

# show
@test sprint(show, IntSet()) == "IntSet([])"
@test sprint(show, IntSet([1,2,3])) == "IntSet([1, 2, 3])"

s = IntSet([1,2,10,20,200,300,1000,10000,10002])
@test last(s) == 10002
@test first(s) == 1
@test length(s) == 9
@test pop!(s) == 10002
@test_throws KeyError pop!(s, -1)
@test length(s) == 8
@test shift!(s) == 1
@test length(s) == 7
@test !in(1,s)
@test !in(10002,s)
@test in(10000,s)
@test in(10000.0,s)
@test !in(10002.0,s)
@test_throws ArgumentError first(IntSet())
@test_throws ArgumentError last(IntSet())
t = copy(s)
sizehint!(t, 20000) #check that hash does not depend on size of internal Array{UInt32, 1}
@test hash(s) == hash(t)

@test setdiff(IntSet([1, 2, 3, 4]), IntSet([2, 4, 5, 6])) == IntSet([1, 3])
@test symdiff(IntSet([1, 2, 3, 4]), IntSet([2, 4, 5, 6])) == IntSet([1, 3, 5, 6])

s2 = IntSet([1, 2, 3, 4])
setdiff!(s2, IntSet([2, 4, 5, 6]))

@test s2 == IntSet([1, 3])

# == with last-bit set (groups.google.com/forum/#!topic/julia-users/vZNjiIEG_sY)
s = IntSet(255)
@test s == s

# issue #7851
@test_throws ArgumentError IntSet(-1)
@test !(-1 in IntSet(1:10))

# # issue #8570
# This requires 2^29 bytes of storage, which is too much for a simple test
# s = IntSet(2^32)
# @test length(s) == 1
# for b in s; b; end

i = IntSet([1, 2, 3])

union!(i, [1, 2])
@test length(i) == 3
union!(i, [3, 4, 5])
@test length(i) == 5

@test_throws KeyError pop!(i, 10)

empty!(i)
@test length(i) == 0

@test_throws ArgumentError symdiff!(i, -3)
@test symdiff!(i, 3) == IntSet([3])
@test symdiff!(i, 257) == IntSet([3, 257])
@test symdiff!(i, [3, 6]) == IntSet([6, 257])

i = IntSet(1:6)
@test symdiff!(i, IntSet([6, 513])) == IntSet([1:5; 513])

i = IntSet([1, 2, 3])
k = IntSet([4, 5])
copy!(k, i)
@test k == i
@test !(k === i)
copy!(k, k)
@test k == i

# unions
i = IntSet([1, 2, 3])
j = union(i)
@test j == i
@test !(j === i)

j = IntSet([4, 5, 6])
@test union(i, j) == IntSet(1:6)

k = IntSet([7, 8, 9])
@test union(i, j, k) == IntSet(1:9)

## intersections
i = IntSet([1, 2, 3])
j = IntSet([4, 5, 6])

@test intersect(i) == i
@test !(intersect(i) === i)

@test intersect(i, j) == IntSet([])
push!(j, 257)
@test intersect(i, j) == IntSet([])
push!(j, 2, 3, 17)
@test intersect(i, j) == IntSet([2, 3])

## equality
i = IntSet([1, 2, 3])
j = IntSet([1, 2, 4])
@test i != j

push!(j, 257)
pop!(j, 257)
@test i != j
@test j != i

@test issubset(IntSet([1, 2, 4]), IntSet(1:10))
@test issubset(IntSet([]), IntSet([]))
@test IntSet([1, 2, 4]) < IntSet(1:10)
@test !(IntSet([]) < IntSet([]))
@test IntSet([1, 2, 4]) <= IntSet(1:10)
@test IntSet([1, 2, 4]) <= IntSet([1, 2, 4])
@test IntSet([]) <= IntSet([])

max_stars_repo_name: habemus-papadum/julia

max_issues_repo_name: habemus-papadum/julia

max_forks_repo_name: habemus-papadum/julia


text:
# abstract type Constraint end
"""
    Dirichlet(u, ∂Ω, f)
    Dirichlet(u, ∂Ω, f, component)

Create a Dirichlet boundary condition on `u` on the `∂Ω` part of
the boundary. `f` is a function that takes two arguments, `x` and `t`
where `x` is the spatial coordinate and `t` is the current time,
and returns the prescribed value. For example, here we create a
Dirichlet condition for the `:u` field, on the faceset called
`∂Ω` and the value given by the `sin` function:

```julia
dbc = Dirichlet(:u, ∂Ω, (x, t) -> sin(t))
```

If `:u` is a vector field we can specify which component the condition
should be applied to by specifying `component`. `component` can be given
either as an integer, or as a vector, for example:

```julia
dbc = Dirichlet(:u, ∂Ω, (x, t) -> sin(t), 1)      # applied to component 1
dbc = Dirichlet(:u, ∂Ω, (x, t) -> sin(t), [1, 3]) # applied to component 1 and 3
```

`Dirichlet` boundary conditions are added to a [`ConstraintHandler`](@ref)
which applies the condition via `apply!`.
"""
struct Dirichlet # <: Constraint
    f::Function # f(x,t) -> value
    faces::Union{Set{Int},Set{Tuple{Int,Int}}}
    field_name::Symbol
    components::Vector{Int} # components of the field
    local_face_dofs::Vector{Int}
    local_face_dofs_offset::Vector{Int}
end
function Dirichlet(field_name::Symbol, faces::Union{Set{Int},Set{Tuple{Int,Int}}}, f::Function, component::Int=1)
    Dirichlet(field_name, faces, f, [component])
end
function Dirichlet(field_name::Symbol, faces::Union{Set{Int},Set{Tuple{Int,Int}}}, f::Function, components::AbstractVector{Int})
    unique(components) == components || error("components not unique: $components")
    # issorted(components) || error("components not sorted: $components")
    return Dirichlet(f, faces, field_name, Vector(components), Int[], Int[])
end

"""
    ConstraintHandler

Collection of constraints.
"""
struct ConstraintHandler{DH<:AbstractDofHandler,T}
    dbcs::Vector{Dirichlet}
    prescribed_dofs::Vector{Int}
    free_dofs::Vector{Int}
    values::Vector{T}
    dofmapping::Dict{Int,Int} # global dof -> index into dofs and values
    dh::DH
    closed::ScalarWrapper{Bool}
end

function ConstraintHandler(dh::AbstractDofHandler)
    @assert isclosed(dh)
    ConstraintHandler(Dirichlet[], Int[], Int[], Float64[], Dict{Int,Int}(), dh, ScalarWrapper(false))
end

function Base.show(io::IO, ::MIME"text/plain", ch::ConstraintHandler)
    println(io, "ConstraintHandler:")
    if !isclosed(ch)
        print(io, "  Not closed!")
    else
        print(io, "  BCs:")
        for dbc in ch.dbcs
            print(io, "\n    ", "Field: ", dbc.field_name, ", ", "Components: ", dbc.components)
        end
    end
end

isclosed(ch::ConstraintHandler) = ch.closed[]
free_dofs(ch::ConstraintHandler) = ch.free_dofs
prescribed_dofs(ch::ConstraintHandler) = ch.prescribed_dofs

"""
    close!(ch::ConstraintHandler)

Close and finalize the `ConstraintHandler`.
"""
function close!(ch::ConstraintHandler)
    fdofs = setdiff(1:ndofs(ch.dh), ch.prescribed_dofs)
    copy!!(ch.free_dofs, fdofs)
    copy!!(ch.prescribed_dofs, unique(ch.prescribed_dofs)) # for v0.7: unique!(ch.prescribed_dofs)
    sort!(ch.prescribed_dofs) # YOLO
    fill!(resize!(ch.values, length(ch.prescribed_dofs)), NaN)
    for i in 1:length(ch.prescribed_dofs)
        ch.dofmapping[ch.prescribed_dofs[i]] = i
    end
    ch.closed[] = true
    return ch
end

function dbc_check(ch::ConstraintHandler, dbc::Dirichlet)
    # check input
    dbc.field_name in getfieldnames(ch.dh) || throw(ArgumentError("field $(dbc.field_name) does not exist in DofHandler, existing fields are $(getfieldnames(ch.dh))"))
    #TODO FIX!!
    #for component in dbc.components
    #    0 < component <= ndim(ch.dh, dbc.field_name) || error("component $component is not within the range of field $field which has $(ndim(ch.dh, field)) dimensions")
    #end
    if length(dbc.faces) == 0
        @warn("added Dirichlet Boundary Condition to set containing 0 entities")
    end
end

"""
    add!(ch::ConstraintHandler, dbc::Dirichlet)

Add a `Dirichlet boundary` condition to the `ConstraintHandler`.
"""
function add!(ch::ConstraintHandler, dbc::Dirichlet)
    dbc_check(ch, dbc)
    field_idx = find_field(ch.dh, dbc.field_name)
    # Extract stuff for the field
    interpolation = ch.dh.field_interpolations[field_idx]
    field_dim = ch.dh.field_dims[field_idx] # TODO: I think we don't need to extract these here ...
    _add!(ch, dbc, dbc.faces, interpolation, field_dim, field_offset(ch.dh, dbc.field_name))
    return ch
end

function _add!(ch::ConstraintHandler, dbc::Dirichlet, bcfaces::Set{Tuple{Int,Int}}, interpolation::Interpolation, field_dim::Int, offset::Int)
    # calculate which local dof index live on each face
    # face `i` have dofs `local_face_dofs[local_face_dofs_offset[i]:local_face_dofs_offset[i+1]-1]
    local_face_dofs = Int[]
    local_face_dofs_offset = Int[1]
    for (i, face) in enumerate(faces(interpolation))
        for fdof in face, d in 1:field_dim
            if d ∈ dbc.components # skip unless this component should be constrained
                push!(local_face_dofs, (fdof-1)*field_dim + d + offset)
            end
        end
        push!(local_face_dofs_offset, length(local_face_dofs) + 1)
    end
    copy!!(dbc.local_face_dofs, local_face_dofs)
    copy!!(dbc.local_face_dofs_offset, local_face_dofs_offset)

    # loop over all the faces in the set and add the global dofs to `constrained_dofs`
    constrained_dofs = Int[]
    #_celldofs = fill(0, ndofs_per_cell(ch.dh))
    for (cellidx, faceidx) in bcfaces
        _celldofs = fill(0, ndofs_per_cell(ch.dh, cellidx))
        celldofs!(_celldofs, ch.dh, cellidx) # extract the dofs for this cell
        r = local_face_dofs_offset[faceidx]:(local_face_dofs_offset[faceidx+1]-1)
        append!(constrained_dofs, _celldofs[local_face_dofs[r]]) # TODO: for-loop over r and simply push! to ch.prescribed_dofs
        @debug println("adding dofs $(_celldofs[local_face_dofs[r]]) to dbc")
    end

    # save it to the ConstraintHandler
    push!(ch.dbcs, dbc)
    append!(ch.prescribed_dofs, constrained_dofs)
end

function _add!(ch::ConstraintHandler, dbc::Dirichlet, bcnodes::Set{Int}, interpolation::Interpolation, field_dim::Int, offset::Int)
    if interpolation !== default_interpolation(getcelltype(ch.dh.grid))
        @warn("adding constraint to nodeset is not recommended for sub/super-parametric approximations.")
    end

    ncomps = length(dbc.components)
    nnodes = getnnodes(ch.dh.grid)
    interpol_points = getnbasefunctions(interpolation)
    _celldofs = fill(0, ndofs_per_cell(ch.dh))
    node_dofs = zeros(Int, ncomps, nnodes)
    visited = falses(nnodes)
    for (cellidx, cell) in enumerate(ch.dh.grid.cells)
        celldofs!(_celldofs, ch.dh, cellidx) # update the dofs for this cell
        for idx in 1:min(interpol_points, length(cell.nodes))
            node = cell.nodes[idx]
            if !visited[node]
                noderange = (offset + (idx-1)*field_dim + 1):(offset + idx*field_dim) # the dofs in this node
                for (i,c) in enumerate(dbc.components)
                    node_dofs[i,node] = _celldofs[noderange[c]]
                    @debug println("adding dof $(_celldofs[noderange[c]]) to node_dofs")
                end
                visited[node] = true
            end
        end
    end

    constrained_dofs = Int[]
    sizehint!(constrained_dofs, ncomps*length(bcnodes))
    sizehint!(dbc.local_face_dofs, length(bcnodes))
    for node in bcnodes
        if !visited[node]
            error("Unable to add a Dirichlet boundary condition to node $node as there are no degrees of freedom on this node.")
        end
        for i in 1:ncomps
            push!(constrained_dofs, node_dofs[i,node])
        end
        push!(dbc.local_face_dofs, node) # use this field to store the node idx for each node
    end

    # save it to the ConstraintHandler
    copy!!(dbc.local_face_dofs_offset, constrained_dofs) # use this field to store the global dofs
    push!(ch.dbcs, dbc)
    append!(ch.prescribed_dofs, constrained_dofs)
end

# Updates the DBC's to the current time `time`
function update!(ch::ConstraintHandler, time::Real=0.0)
    @assert ch.closed[]
    for dbc in ch.dbcs
        field_idx = find_field(ch.dh, dbc.field_name)
        # Function barrier
        _update!(ch.values, dbc.f, dbc.faces, dbc.field_name, dbc.local_face_dofs, dbc.local_face_dofs_offset,
                 dbc.components, ch.dh, ch.dh.bc_values[field_idx], ch.dofmapping, convert(Float64, time))
    end
end

# for faces
function _update!(values::Vector{Float64}, f::Function, faces::Set{Tuple{Int,Int}}, field::Symbol, local_face_dofs::Vector{Int}, local_face_dofs_offset::Vector{Int},
                  components::Vector{Int}, dh::DofHandler{dim,N,T,M}, facevalues::BCValues,
                  dofmapping::Dict{Int,Int}, time::Float64) where {dim,N,T,M}
    grid = dh.grid

    xh = zeros(Vec{dim, T}, N) # pre-allocate
    _celldofs = fill(0, ndofs_per_cell(dh))

    for (cellidx, faceidx) in faces
        getcoordinates!(xh, grid, cellidx)
        celldofs!(_celldofs, dh, cellidx) # update global dofs for this cell

        # no need to reinit!, enough to update current_face since we only need geometric shape functions M
        facevalues.current_face[] = faceidx

        # local dof-range for this face
        r = local_face_dofs_offset[faceidx]:(local_face_dofs_offset[faceidx+1]-1)
        counter = 1

        for location in 1:getnquadpoints(facevalues)
            x = spatial_coordinate(facevalues, location, xh)
            bc_value = f(x, time)
            @assert length(bc_value) == length(components)

            for i in 1:length(components)
                # find the global dof
                globaldof = _celldofs[local_face_dofs[r[counter]]]
                counter += 1

                dbc_index = dofmapping[globaldof]
                values[dbc_index] = bc_value[i]
                @debug println("prescribing value $(bc_value[i]) on global dof $(globaldof)")
            end
        end
    end
end

# for nodes
function _update!(values::Vector{Float64}, f::Function, nodes::Set{Int}, field::Symbol, nodeidxs::Vector{Int}, globaldofs::Vector{Int},
                  components::Vector{Int}, dh::DofHandler{dim,N,T,M}, facevalues::BCValues,
                  dofmapping::Dict{Int,Int}, time::Float64) where {dim,N,T,M}
    counter = 1
    for (idx, nodenumber) in enumerate(nodeidxs)
        x = dh.grid.nodes[nodenumber].x
        bc_value = f(x, time)
        @assert length(bc_value) == length(components)
        for v in bc_value
            globaldof = globaldofs[counter]
            counter += 1
            dbc_index = dofmapping[globaldof]
            values[dbc_index] = v
            @debug println("prescribing value $(v) on global dof $(globaldof)")
        end
    end
end

# Saves the dirichlet boundary conditions to a vtkfile.
# Values will have a 1 where bcs are active and 0 otherwise
function WriteVTK.vtk_point_data(vtkfile, ch::ConstraintHandler)
    unique_fields = []
    for dbc in ch.dbcs
        push!(unique_fields, dbc.field_name)
    end
    unique_fields = unique(unique_fields) # TODO v0.7: unique!(unique_fields)

    for field in unique_fields
        nd = ndim(ch.dh, field)
        data = zeros(Float64, nd, getnnodes(ch.dh.grid))
        for dbc in ch.dbcs
            dbc.field_name != field && continue
            if eltype(dbc.faces) <: Tuple
                for (cellidx, faceidx) in dbc.faces
                    for facenode in faces(ch.dh.grid.cells[cellidx])[faceidx]
                        for component in dbc.components
                            data[component, facenode] = 1
                        end
                    end
                end
            else
                for nodeidx in dbc.faces
                    for component in dbc.components
                        data[component, nodeidx] = 1
                    end
                end
            end
        end
        vtk_point_data(vtkfile, data, string(field, "_bc"))
    end
    return vtkfile
end

function apply!(v::AbstractVector, ch::ConstraintHandler)
    @assert length(v) == ndofs(ch.dh)
    v[ch.prescribed_dofs] = ch.values # .= ??
    return v
end

function apply_zero!(v::AbstractVector, ch::ConstraintHandler)
    @assert length(v) == ndofs(ch.dh)
    v[ch.prescribed_dofs] .= 0
    return v
end

function apply!(K::Union{SparseMatrixCSC,Symmetric}, ch::ConstraintHandler)
    apply!(K, eltype(K)[], ch, true)
end

function apply_zero!(K::Union{SparseMatrixCSC,Symmetric}, f::AbstractVector, ch::ConstraintHandler)
    apply!(K, f, ch, true)
end

@enum(ApplyStrategy, APPLY_TRANSPOSE, APPLY_INPLACE)

function apply!(KK::Union{SparseMatrixCSC,Symmetric}, f::AbstractVector, ch::ConstraintHandler, applyzero::Bool=false;
                strategy::ApplyStrategy=APPLY_TRANSPOSE)
    K = isa(KK, Symmetric) ? KK.data : KK
    @assert length(f) == 0 || length(f) == size(K, 1)
    @boundscheck checkbounds(K, ch.prescribed_dofs, ch.prescribed_dofs)
    @boundscheck length(f) == 0 || checkbounds(f, ch.prescribed_dofs)

    m = meandiag(K) # Use the mean of the diagonal here to not ruin things for iterative solver
    @inbounds for i in 1:length(ch.values)
        d = ch.prescribed_dofs[i]
        v = ch.values[i]

        if !applyzero && v != 0
            for j in nzrange(K, d)
                f[K.rowval[j]] -= v * K.nzval[j]
            end
        end
    end
    zero_out_columns!(K, ch.prescribed_dofs)
    if strategy == APPLY_TRANSPOSE
        K′ = copy(K)
        transpose!(K′, K)
        zero_out_columns!(K′, ch.prescribed_dofs)
        transpose!(K, K′)
    elseif strategy == APPLY_INPLACE
        K[ch.prescribed_dofs, :] = 0
    else
        error("Unknown apply strategy")
    end
    @inbounds for i in 1:length(ch.values)
        d = ch.prescribed_dofs[i]
        v = ch.values[i]
        K[d, d] = m
        # We will only enter here with an empty f vector if we have assured that v == 0 for all dofs
        if length(f) != 0
            vz = applyzero ? zero(eltype(f)) : v
            f[d] = vz * m
        end
    end
end

# columns need to be stored entries, this is not checked
function zero_out_columns!(K, dofs::Vector{Int}) # can be removed in 0.7 with #24711 merged
    @debug @assert issorted(dofs)
    for col in dofs
        r = nzrange(K, col)
        K.nzval[r] .= 0.0
    end
end

function meandiag(K::AbstractMatrix)
    z = zero(eltype(K))
    for i in 1:size(K, 1)
        z += abs(K[i, i])
    end
    return z / size(K, 1)
end

# TODO Remove - almost identical to the previous
function _update!(values::Vector{Float64}, f::Function, faces::Set{Tuple{Int,Int}}, field::Symbol, local_face_dofs::Vector{Int}, local_face_dofs_offset::Vector{Int},
                  components::Vector{Int}, dh::MixedDofHandler{dim,C,T}, facevalues::BCValues,
                  dofmapping::Dict{Int,Int}, time::Float64) where {dim,C,T}
    grid = dh.grid
    #xh = zeros(Vec{dim, T}, N) # pre-allocate
    for (cellidx, faceidx) in faces
        _celldofs = fill(0, ndofs_per_cell(dh, cellidx))
        xh = getcoordinates(grid, cellidx)
        # zeros(Vec{dim, T}, N) # pre-allocate
        #getcoordinates!(xh, grid, cellidx)
        celldofs!(_celldofs, dh, cellidx) # update global dofs for this cell

        # no need to reinit!, enough to update current_face since we only need geometric shape functions M
        facevalues.current_face[] = faceidx

        # local dof-range for this face
        r = local_face_dofs_offset[faceidx]:(local_face_dofs_offset[faceidx+1]-1)
        counter = 1

        for location in 1:getnquadpoints(facevalues)
            x = spatial_coordinate(facevalues, location, xh)
            bc_value = f(x, time)
            @assert length(bc_value) == length(components)

            for i in 1:length(components)
                # find the global dof
                globaldof = _celldofs[local_face_dofs[r[counter]]]
                counter += 1

                dbc_index = dofmapping[globaldof]
                values[dbc_index] = bc_value[i]
                @debug println("prescribing value $(bc_value[i]) on global dof $(globaldof)")
            end
        end
    end
end

import JuAFEM.add!, JuAFEM.update!
function add!(ch::ConstraintHandler, fh::FieldHandler, dbc::Dirichlet)
    #dbc_check(ch, dbc)
    field_idx = find_field(fh, dbc.field_name)
    # Extract stuff for the field
    interpolation = getfieldinterpolations(fh)[field_idx]
    field_dim = getfielddims(fh)[field_idx]
    JuAFEM._add!(ch, dbc, dbc.faces, interpolation, field_dim, field_offset(fh, dbc.field_name))
    return ch
end

# Updates the DBC's to the current time `time`
function update!(ch::ConstraintHandler, dbcmap::Dict, time::Real=0.0)
    @assert ch.closed[]
    for dbc in ch.dbcs
        fh = dbcmap[dbc]
        field_idx = find_field(fh, dbc.field_name)
        _update!(ch.values, dbc.f, dbc.faces, dbc.field_name, dbc.local_face_dofs,
                 dbc.local_face_dofs_offset, dbc.components, ch.dh,
                 fh.bc_values[field_idx], ch.dofmapping, convert(Float64, time))
    end
end


text:
###################################################
## network/solvers/fraclp.jl
## Fractional LP (symmetric cost function) that finds new traveltimes
###################################################

"""
    fraclpTimes :
    optimize travel times to minimize L1 error from data with given paths
"""
function fraclpTimes(s::IterativeState; args...)
    g = s.data.network.graph
    paths = s.paths
    tripData = s.trips

    #Create the model (will be changed to avoid hard-coded parameters)
    # !BarConvTol needs to be changed
    m = Model(solver = GurobiSolver(TimeLimit=10000, Crossover=0, Method=2, NumericFocus=3, ScaleFlag=0))#Crossover=0, Method=3, BarConvTol=1e-6; args...))

    # DECISION VARIABLES
    # Road times
    @variable(m, t[i=vertices(g), j=outneighbors(g,i)] >= 0)
    # Absolute difference between tripData times and computed times
    @variable(m, epsilon[d=eachindex(tripData)] >= 0)
    @variable(m, y >= 0)

    # OBJECTIVE
    @objective(m, Min, sum(tripData[d].weight * epsilon[d] for d in eachindex(tripData)))

    # CONSTRAINTS
    # absolute values contraints (define epsilon), equal to time of first path
    @constraint(m, epsLower[d=eachindex(tripData)],
        epsilon[d] >=
        sum(paths[d][1][edge] * t[src(edge), dst(edge)] for edge in keys(paths[d][1])) - y * tripData[d].time
        )
    @constraint(m, epsUpper[d=eachindex(tripData)],
        epsilon[d] >=
        - sum(paths[d][1][edge] * t[src(edge), dst(edge)] for edge in keys(paths[d][1])) + y * tripData[d].time
        )

    # inequality constraints
    @constraint(m, inequalityPath[d=eachindex(tripData), p=1:(length(paths[d])-1)],
        sum(paths[d][p+1][edge] * t[src(edge), dst(edge)], edge=keys(paths[d][p+1])) >=
        sum(paths[d][1][edge] * t[src(edge), dst(edge)], edge=keys(paths[d][1]))
        )

    # fractional programming constraint
    @constraint(m, fracProgram,
        sum(paths[d][1][edge] * t[src(edge), dst(edge)], d=eachindex(tripData) for edge in keys(paths[d][1]))
        + y * sum(tripData[d].time for d in eachindex(tripData)) == 1
        )

    # new bounds on edge velocities
    @constraint(m, speedLimits[i=vertices(g), j=outneighbors(g,i)],
        t[i,j] >= s.data.minTimes[i,j] * y)

    # SOLVE LP
    status = solve(m)
    if status == :Infeasible
        buildInternalModel(m)
        print_iis_gurobi(m)
    end
    times = getvalue(t)
    yVal = getvalue(y)

    # Export result as sparse matrix
    result = spzeros(Float64, nv(g), nv(g))
    for i in vertices(g), j in outneighbors(g,i)
        result[i,j] = times[i,j] / yVal
    end

    return result
end


text:
"""
Solve any non-conjugate likelihood using Variational Inference
by making a numerical approximation (quadrature or MC integration)
of the expected log-likelihood ad its gradients
"""
abstract type NumericalVI{T<:Real} <: Inference{T} end

include("quadratureVI.jl")
include("MCVI.jl")

""" `NumericalVI(integration_technique::Symbol=:quad;ϵ::T=1e-5,nMC::Integer=1000,nGaussHermite::Integer=20,optimizer::Optimizer=Adam(α=0.1))`

General constructor for Variational Inference via numerical approximation.

**Argument**

    -`integration_technique::Symbol` : Method of approximation can be `:quad` for quadrature see [QuadratureVI](@ref) or `:mc` for MC integration see [MCIntegrationVI](@ref)

**Keyword arguments**

    - `ϵ::T` : convergence criteria, which can be user defined
    - `nMC::Int` : Number of samples per data point for the integral evaluation (for the MCIntegrationVI)
    - `nGaussHermite::Int` : Number of points for the integral estimation (for the QuadratureVI)
    - `optimizer::Optimizer` : Optimizer used for the variational updates. Should be an Optimizer object from the [GradDescent.jl]() package. Default is `Adam()`
"""
function NumericalVI(integration_technique::Symbol=:quad;ϵ::T=1e-5,nMC::Integer=1000,nGaussHermite::Integer=20,optimizer::Optimizer=Adam(α=0.1)) where {T<:Real}
    if integration_technique == :quad
        QuadratureVI{T}(ϵ,nGaussHermite,0,optimizer,false)
    elseif integration_technique == :mc
        MCIntegrationVI{T}(ϵ,nMC,0,optimizer,false)
    else
        @error "Only possible integration techniques are quadrature : :quad or mcmc integration :mcmc"
    end
end

""" `NumericalSVI(integration_technique::Symbol=:quad;ϵ::T=1e-5,nMC::Integer=1000,nGaussHermite::Integer=20,optimizer::Optimizer=Adam(α=0.1))`

General constructor for Stochastic Variational Inference via numerical approximation.

**Argument**

    -`nMinibatch::Integer` : Number of samples per mini-batches
    -`integration_technique::Symbol` : Method of approximation can be `:quad` for quadrature see [QuadratureVI](@ref) or `:mc` for MC integration see [MCIntegrationVI](@ref)

**Keyword arguments**

    - `ϵ::T` : convergence criteria, which can be user defined
    - `nMC::Int` : Number of samples per data point for the integral evaluation (for the MCIntegrationVI)
    - `nGaussHermite::Int` : Number of points for the integral estimation (for the QuadratureVI)
    - `optimizer::Optimizer` : Optimizer used for the variational updates. Should be an Optimizer object from the [GradDescent.jl]() package. Default is `Adam()`
"""
function NumericalSVI(nMinibatch::Integer,integration_technique::Symbol=:quad;ϵ::T=1e-5,nMC::Integer=200,nGaussHermite::Integer=20,optimizer::Optimizer=Adam(α=0.1)) where {T<:Real}
    if integration_technique == :quad
        QuadratureVI{T}(ϵ,nGaussHermite,0,optimizer,true,nMinibatch)
    elseif integration_technique == :mc
        MCIntegrationVI{T}(ϵ,nMC,0,optimizer,true,nMinibatch)
    else
        @error "Only possible integration techniques are quadrature : :quad or mcmc integration :mc"
    end
end

function Base.show(io::IO,inference::NumericalVI{T}) where T
    print(io,"$(inference.Stochastic ? "Stochastic numerical" : "Numerical") inference by $(isa(inference,MCIntegrationVI) ? "Monte Carlo Integration" : "Quadrature")")
end

function init_inference(inference::NumericalVI{T},nLatent::Integer,nFeatures::Integer,nSamples::Integer,nSamplesUsed::Integer) where {T<:Real}
    inference.nSamples = nSamples
    inference.nSamplesUsed = nSamplesUsed
    inference.MBIndices = 1:nSamplesUsed
    inference.ρ = nSamples/nSamplesUsed
    inference.HyperParametersUpdated = true
    inference.optimizer_η₁ = [copy(inference.optimizer_η₁[1]) for _ in 1:nLatent]
    inference.optimizer_η₂ = [copy(inference.optimizer_η₂[1]) for _ in 1:nLatent]
    inference.∇η₁ = [zeros(T,nFeatures) for _ in 1:nLatent];
    inference.∇η₂ = [Symmetric(Diagonal(ones(T,nFeatures))) for _ in 1:nLatent]
    inference.∇μE = [zeros(T,nSamplesUsed) for _ in 1:nLatent];
    inference.∇ΣE = [zeros(T,nSamplesUsed) for _ in 1:nLatent]
    return inference
end

function variational_updates!(model::VGP{<:Likelihood,<:NumericalVI})
    compute_grad_expectations!(model)
    natural_gradient!(model)
    global_update!(model)
end

function variational_updates!(model::SVGP{<:Likelihood,<:NumericalVI}) where {L<:Likelihood,T}
    compute_grad_expectations!(model)
    natural_gradient!(model)
    global_update!(model)
end

function natural_gradient!(model::VGP{<:Likelihood,<:NumericalVI})
    model.inference.∇η₁ .= model.Σ.*(model.inference.∇μE .- model.invKnn.*model.μ)
    model.inference.∇η₂ .= Symmetric.(Diagonal.(model.inference.∇ΣE).-0.5.*model.invKnn .- model.η₂)
end

function natural_gradient!(model::SVGP{<:Likelihood,<:NumericalVI})
    model.inference.∇η₁ .= model.Σ.*(model.inference.ρ.*transpose.(model.κ).*model.inference.∇μE .- model.invKmm.*model.μ)
    model.inference.∇η₂ .= Symmetric.(model.inference.ρ.*transpose.(model.κ).*Diagonal.(model.inference.∇ΣE).*model.κ.-0.5.*model.invKmm .- model.η₂)
end

function global_update!(model::AbstractGP{<:Likelihood,<:NumericalVI})
    model.η₁ .= model.η₁ .+ update.(model.inference.optimizer_η₁,model.inference.∇η₁)
    for k in 1:model.nLatent
        Δ = update(model.inference.optimizer_η₂[k],model.inference.∇η₂[k])
        α=1.0
        while true
            try
                @assert isposdef(-(model.η₂[k]+α*Δ))
                model.η₂[k] = Symmetric(model.η₂[k]+α*Δ)
                model.η₁[k] .+= update(model.inference.optimizer_η₁[k],model.inference.∇η₁[k])
                break;
            catch e
                if isa(e,AssertionError)
                    println("Error, results not pos def with α=$α")
                    α *= 0.5
                    if α < 1e-6
                        @error "α too small, stopping loop"
                        rethrow()
                    end
                else
                    rethrow()
                end
            end
        end
        if isa(model.inference.optimizer_η₂[k],Adam)
            model.inference.optimizer_η₂[k].α = min(model.inference.optimizer_η₂[k].α*α*2.0,1.0)
            # model.inference.optimizer_η₁[k].α = min(model.inference.optimizer_η₁[k].α*α*2.0,1.0)
        elseif isa(model.inference.optimizer_η₂[k],VanillaGradDescent)
            # model.inference.optimizer_η₂[k].η = min(model.inference.optimizer_η₂[k].η*α*2.0,1.0)
            # model.inference.optimizer_η₁[k].η = min(model.inference.optimizer_η₁[k].η*α*2.0,1.0)
        elseif isa(model.inference.optimizer_η₂[k],ALRSVI)
        elseif isa(model.inference.optimizer_η₂[k],InverseDecay)
        end
    end
    model.Σ .= -0.5.*inv.(model.η₂)
    model.μ .= model.η₁
    # model.μ .= model.Σ.*model.η₁
end

function ELBO(model::AbstractGP{<:Likelihood,<:NumericalVI})
    return expecLogLikelihood(model) - GaussianKL(model)
end

function expec_μ(model::AbstractGP{<:Likelihood,<:NumericalVI},index::Integer)
    return model.inference.∇μE[index]
end

function expec_μ(model::AbstractGP{<:Likelihood,<:NumericalVI})
    return model.inference.∇μE
end

function expec_Σ(model::AbstractGP{<:Likelihood,<:NumericalVI},index::Integer)
    return model.inference.∇ΣE[index]
end

function expec_Σ(model::AbstractGP{<:Likelihood,<:NumericalVI})
    return model.inference.∇ΣE
end

function global_update!(model::SVGP{L,NumericalVI{T}}) where {L<:Likelihood,T}
    if model.inference.Stochastic
    else
        model.η₁ .= model.inference.∇η₁ .+ model.η₁
        model.η₂ .= Symmetric.(model.inference.∇η₂ .+ model.η₂)
    end
    model.Σ .= inv.(model.η₂)*(-0.5)
    model.μ .= model.Σ.*model.η₁
end

function convert(::Type{T1},x::T2) where {T1<:VGP{<:Likelihood,T3} where {T3<:NumericalVI},T2<:VGP{<:Likelihood,<:AnalyticVI}}
    #TODO Check if likelihood is compatible
    inference = T3(x.inference.ϵ,x.inference.nIter,x.inference.optimizer,defaultn(T3),x.inference.Stochastic,x.inference.nSamples,x.inference.nSamplesUsed,x.inference.MBIndices,x.inference.ρ,x.inference.HyperParametersUpdated,x.inference.∇η₁,x.inference.∇η₂,copy(expec_μ(x)),copy(expec_Σ(x)))
    likelihood =isaugmented(x.likelihood) ? remove_augmentation(x.likelihood) : likelihood
    return T1(x.X,x.y,x.nSample,x.nDim,x.nFeature,x.nLatent,x.IndependentPriors,x.nPrior,x.μ,x.Σ,x.η₁,x.η₂,x.Knn,x.invKnn,x.kernel,likelihood,inference,x.verbose,x.Autotuning,x.atfrequency,x.Trained)
end


text:
batch_count = 10
max_batch_items = 3
max_ndims = 3
sizes = 20:50

rand_shape = () -> rand(sizes, rand(1:max_ndims))
batches = [[rand(Float32, rand_shape()...) for _ in 1:rand(1:max_batch_items)]
                                           for _ in 1:batch_count]
cubatches = CuIterator(batch for batch in batches) # ensure generators are accepted

previous_cubatch = missing
for (batch, cubatch) in zip(batches, cubatches)
    global previous_cubatch
    @test ismissing(previous_cubatch) || all(x -> x.storage === nothing, previous_cubatch)
    @test batch == Array.(cubatch)
    @test all(x -> x isa CuArray, cubatch)
    previous_cubatch = cubatch
end


text:
# uncomment if first time using LightXML
# Pkg.add("LightXML");
using LightXML;

xWidth = 19900.0
yWidth = 13500.0
latRoot = 46.9877911
lonRoot = 7.131336
pathString = "M"

function mappedLongitude(lon)

    if (lon < 0)
        (((180 + (180 - abs(lon))) / 360) * xWidth)
    else
        ((lon / 360) * xWidth)
    end
end

function mappedLatitude(lat)
    with_rounding(Float64,RoundDown) do
        ((180 - (lat + 90)) / 180) * yWidth
    end
end

gpxDoc = parse_file("/Users/rene/Development/gps-test/github/gpx/1k.gpx")

xroot = root(gpxDoc)
trk = get_elements_by_tagname(xroot, "trk")[1]
trkseg = get_elements_by_tagname(trk, "trkseg")[1]
trkpt = get_elements_by_tagname(trkseg, "trkpt")
nodeCounter = 1
lastMappedLongitude = 1.0
lastMappedLatitude = 1.0
for c in child_nodes(trkseg)
    if is_elementnode(c)
        #thisNode = XMLElement(c)
        ad = attributes_dict(XMLElement(c))
        if (nodeCounter == 1)
            lastMappedLongitude = mappedLongitude(parse(Float32,(ad["lon"])))
            lastMappedLatitude = mappedLatitude(parse(Float32,(ad["lat"])))
            mappedLonRoot = mappedLongitude(lonRoot)
            mappedLatRoot = mappedLatitude(latRoot)
            pathString = string(pathString,round(lastMappedLongitude - mappedLonRoot,5) , "," , round(lastMappedLatitude - mappedLatRoot,5))
        else
            currentMappedLongitude = mappedLongitude(parse(Float32,(ad["lon"])))
            currentMappedLatitude = mappedLatitude(parse(Float32,(ad["lat"])))
            pathString = string(pathString,"l",round((currentMappedLongitude - lastMappedLongitude),5) , "," , round((currentMappedLatitude - lastMappedLatitude),5))
            lastMappedLongitude = currentMappedLongitude
            lastMappedLatitude = currentMappedLatitude
        end
        nodeCounter = nodeCounter + 1
    end
end

println(pathString)

max_stars_repo_name: OpenDGPS/gpx-to-svg

max_issues_repo_name: OpenDGPS/gpx-to-svg

max_forks_repo_name: OpenDGPS/gpx-to-svg


text:
using SimpleConstants
using Test

@testset "SimpleConstants.jl" begin
    λ = rand()  # wavelength
    @test 2πc₀/λ == 2π * c₀ / λ
    @test ℎc₀/λ == ℎ * c₀ / λ

    @test exp(i2π/3) == exp(im * 2π / 3)
    @test exp(iπ/3) == exp(im * π / 3)

    @test e⁻ < 0
    @test e⁺ > 0
    @test e⁻ + e⁺ == 0

    @test 180° ≈ π

    λ_inµm = rand()
    E_ineV = meter2eV(λ_inµm * micro)
    @test eV2meter(E_ineV) ≈ λ_inµm * micro
    @test meter2eV(λ_inµm * micro) ≈ E_ineV
    @test round(λ_inµm * E_ineV, digits=2) == 1.24  # (λ in µm) times (E in eV) is roughly 1.24
end


text:
function sqlvalue(values, i)
    temp_val_ptr = unsafe_load(values, i)
    valuetype = sqlite3_value_type(temp_val_ptr)

    if valuetype == SQLITE_INTEGER
        if Sys.WORD_SIZE == 64
            return sqlite3_value_int64(temp_val_ptr)
        else
            return sqlite3_value_int(temp_val_ptr)
        end
    elseif valuetype == SQLITE_FLOAT
        return sqlite3_value_double(temp_val_ptr)
    elseif valuetype == SQLITE_TEXT
        return unsafe_string(sqlite3_value_text(temp_val_ptr))
    elseif valuetype == SQLITE_BLOB
        nbytes = sqlite3_value_bytes(temp_val_ptr)
        blob = sqlite3_value_blob(temp_val_ptr)
        buf = zeros(UInt8, nbytes)
        unsafe_copyto!(pointer(buf), convert(Ptr{UInt8}, blob), nbytes)
        return sqldeserialize(buf)
    else
        return missing
    end
end

sqlreturn(context, ::Missing)           = sqlite3_result_null(context)
sqlreturn(context, val::Int32)          = sqlite3_result_int(context, val)
sqlreturn(context, val::Int64)          = sqlite3_result_int64(context, val)
sqlreturn(context, val::Float64)        = sqlite3_result_double(context, val)
sqlreturn(context, val::AbstractString) = sqlite3_result_text(context, val)
sqlreturn(context, val::Vector{UInt8})  = sqlite3_result_blob(context, val)

sqlreturn(context, val::Bool) = sqlreturn(context, Int(val))
sqlreturn(context, val) = sqlreturn(context, sqlserialize(val))

# Internal method for generating an SQLite scalar function from
# a Julia function name
function scalarfunc(func,fsym=Symbol(string(func)))
    # check if name defined in Base so we don't clobber Base methods
    nm = isdefined(Base,fsym) ? :(Base.$fsym) : fsym
    return quote
        #nm needs to be a symbol or expr, i.e. :sin or :(Base.sin)
        function $(nm)(context::Ptr{Cvoid}, nargs::Cint, values::Ptr{Ptr{Cvoid}})
            args = [SQLite.sqlvalue(values, i) for i in 1:nargs]
            ret = $(func)(args...)
            SQLite.sqlreturn(context, ret)
            nothing
        end
        return $(nm)
    end
end
function scalarfunc(expr::Expr)
    f = eval(expr)
    return scalarfunc(f)
end

# convert a byteptr to an int, assumes little-endian
function bytestoint(ptr::Ptr{UInt8}, start::Int, len::Int)
    s = 0
    for i in start:start+len-1
        v = unsafe_load(ptr, i)
        s += v * 256^(i - start)
    end

    # swap byte-order on big-endian machines
    # TODO: this desperately needs testing on a big-endian machine!!!!!
    return htol(s)
end

function stepfunc(init, func, fsym=Symbol(string(func)*"_step"))
    nm = isdefined(Base,fsym) ? :(Base.$fsym) : fsym
    return quote
        function $(nm)(context::Ptr{Cvoid}, nargs::Cint, values::Ptr{Ptr{Cvoid}})
            args = [sqlvalue(values, i) for i in 1:nargs]

            intsize = sizeof(Int)
            ptrsize = sizeof(Ptr)
            acsize = intsize + ptrsize
            acptr = convert(Ptr{UInt8}, sqlite3_aggregate_context(context, acsize))

            # acptr will be zeroed-out if this is the first iteration
            ret = ccall(
                :memcmp, Cint, (Ptr{UInt8}, Ptr{UInt8}, Cuint),
                zeros(UInt8, acsize), acptr, acsize,
            )
            if ret == 0
                acval = $(init)
                valsize = 256
                # avoid the garbage collector using malloc
                valptr = convert(Ptr{UInt8}, Libc.malloc(valsize))
                valptr == C_NULL && throw(SQLiteException("memory error"))
            else
                # size of serialized value is first sizeof(Int) bytes
                valsize = bytestoint(acptr, 1, intsize)
                # ptr to serialized value is last sizeof(Ptr) bytes
                valptr = reinterpret(
                    Ptr{UInt8}, bytestoint(acptr, intsize+1, ptrsize)
                )
                # deserialize the value pointed to by valptr
                acvalbuf = zeros(UInt8, valsize)
                unsafe_copyto!(pointer(acvalbuf), valptr, valsize)
                acval = sqldeserialize(acvalbuf)
            end

            local funcret
            try
                funcret = sqlserialize($(func)(acval, args...))
            catch
                Libc.free(valptr)
                rethrow()
            end

            newsize = sizeof(funcret)
            if newsize > valsize
                # TODO: increase this in a cleverer way?
                tmp = convert(Ptr{UInt8}, Libc.realloc(valptr, newsize))
                if tmp == C_NULL
                    Libc.free(valptr)
                    throw(SQLiteException("memory error"))
                else
                    valptr = tmp
                end
            end
            # copy serialized return value
            unsafe_copyto!(valptr, pointer(funcret), newsize)

            # copy the size of the serialized value
            unsafe_copyto!(
                acptr,
                pointer(reinterpret(UInt8, [newsize])),
                intsize
            )
            # copy the address of the pointer to the serialized value
            valarr = reinterpret(UInt8, [valptr])
            for i in 1:length(valarr)
                unsafe_store!(acptr, valarr[i], intsize+i)
            end
            nothing
        end
        return $(nm)
    end
end

function finalfunc(init, func, fsym=Symbol(string(func)*"_final"))
    nm = isdefined(Base,fsym) ? :(Base.$fsym) : fsym
    return quote
        function $(nm)(context::Ptr{Cvoid}, nargs::Cint, values::Ptr{Ptr{Cvoid}})
            acptr = convert(Ptr{UInt8}, sqlite3_aggregate_context(context, 0))

            # step function wasn't run
            if acptr == C_NULL
                sqlreturn(context, $(init))
            else
                intsize = sizeof(Int)
                ptrsize = sizeof(Ptr)
                acsize = intsize + ptrsize

                # load size
                valsize = bytestoint(acptr, 1, intsize)
                # load ptr
                valptr = reinterpret(
                    Ptr{UInt8}, bytestoint(acptr, intsize+1, ptrsize)
                )

                # load value
                acvalbuf = zeros(UInt8, valsize)
                unsafe_copyto!(pointer(acvalbuf), valptr, valsize)
                acval = sqldeserialize(acvalbuf)

                local ret
                try
                    ret = $(func)(acval)
                finally
                    Libc.free(valptr)
                end
                sqlreturn(context, ret)
            end
            nothing
        end
        return $(nm)
    end
end

# User-facing macro for convenience in registering a simple function
# with no configurations needed
macro register(db, func)
    :(register($(esc(db)), $(esc(func))))
end

# User-facing method for registering a Julia function to be used within SQLite
function register(db, func::Function; nargs::Int=-1, name::AbstractString=string(func), isdeterm::Bool=true)
    @assert nargs <= 127 "use -1 if > 127 arguments are needed"
    # assume any negative number means a varargs function
    nargs < -1 && (nargs = -1)
    @assert sizeof(name) <= 255 "size of function name must be <= 255"

    f = eval(scalarfunc(func,Symbol(name)))

    cfunc = @cfunction($f, Cvoid, (Ptr{Cvoid}, Cint, Ptr{Ptr{Cvoid}}))
    # TODO: allow the other encodings
    enc = SQLITE_UTF8
    enc = isdeterm ? enc | SQLITE_DETERMINISTIC : enc

    @CHECK db sqlite3_create_function_v2(
        db.handle, name, nargs, enc, C_NULL, cfunc, C_NULL, C_NULL, C_NULL
    )
end

# as above but for aggregate functions
newidentity() = @eval x->x
function register(
    db, init, step::Function, final::Function= newidentity();
    nargs::Int=-1, name::AbstractString=string(step), isdeterm::Bool=true
)
    @assert nargs <= 127 "use -1 if > 127 arguments are needed"
    nargs < -1 && (nargs = -1)
    @assert sizeof(name) <= 255 "size of function name must be <= 255 chars"

    s = eval(stepfunc(init, step, Base.nameof(step)))
    cs = @cfunction($s, Cvoid, (Ptr{Cvoid}, Cint, Ptr{Ptr{Cvoid}}))
    f = eval(finalfunc(init, final, Base.nameof(final)))
    cf = @cfunction($f, Cvoid, (Ptr{Cvoid}, Cint, Ptr{Ptr{Cvoid}}))

    enc = SQLITE_UTF8
    enc = isdeterm ? enc | SQLITE_DETERMINISTIC : enc

    @CHECK db sqlite3_create_function_v2(
        db.handle, name, nargs, enc, C_NULL, C_NULL, cs, cf, C_NULL
    )
end

# annotate types because the MethodError makes more sense that way
regexp(r::AbstractString, s::AbstractString) = occursin(Regex(r), s)
# macro for preserving the special characters in a string
macro sr_str(s) s end


text:
# tests for the various forms of stParallel for the ScatteringTransform
# exit()
# using Revise
using Distributed
addprocs(min((Sys.CPU_THREADS) - 2 - nprocs(), 2))
@everywhere using Interpolations, ContinuousWavelets, LinearAlgebra
@everywhere using FFTW
@everywhere using ScatteringTransform
@everywhere using SharedArrays
using Test
using Flux, FourierFilterFlux, CUDA, AbstractFFTs, LinearAlgebra, Statistics
using Zygote

include("planTests.jl")
include("cwtTests.jl")
include("pathTests.jl")
include("fluxtests.jl")


text:
using SafeTestsets

@time begin
    @safetestset "UnivariateDiscrete" begin include("univariate_discrete.jl") end
    @safetestset "UnivariateContinuous" begin include("univariate_continuous.jl") end
    @safetestset "Mixtures" begin include("mixtures.jl") end
end


text:
using Plots;
pyplot();
using Distributions
using LinearAlgebra
using KernelFunctions
# Translational invariants kernels

default(; lw=3.0, titlefontsize=28, tickfontsize=18)

x₀ = 0.0;
l = 0.1;
n_grid = 101
fill(x₀, n_grid, 1)
xrange = reshape(collect(range(-3, 3; length=n_grid)), :, 1)

k = SqExponentialKernel() ∘ ScaleTransform(1.0)
K1 = kernelmatrix(k, xrange; obsdim=1)
p = heatmap(
    K1;
    yflip=true,
    colorbar=false,
    framestyle=:none,
    background_color=RGBA(0.0, 0.0, 0.0, 0.0),
)
savefig(joinpath(@__DIR__, "src", "assets", "heatmap_sqexp.png"))

k = @kernel Matern32Kernel FunctionTransform(x -> (sin.(x)) .^ 2)
K2 = kernelmatrix(k, xrange; obsdim=1)
p = heatmap(
    K2;
    yflip=true,
    colorbar=false,
    framestyle=:none,
    background_color=RGBA(0.0, 0.0, 0.0, 0.0),
)
savefig(joinpath(@__DIR__, "src", "assets", "heatmap_matern.png"))

k = PolynomialKernel(; c=0.0, d=2.0) ∘ LinearTransform(randn(3, 1))
K3 = kernelmatrix(k, xrange; obsdim=1)
p = heatmap(
    K3;
    yflip=true,
    colorbar=false,
    framestyle=:none,
    background_color=RGBA(0.0, 0.0, 0.0, 0.0),
)
savefig(joinpath(@__DIR__, "src", "assets", "heatmap_poly.png"))

k =
    0.5 * SqExponentialKernel() * (LinearKernel() ∘ ScaleTransform(0.5)) +
    0.4 * (@kernel Matern32Kernel() FunctionTransform(x -> sin.(x)))
K4 = kernelmatrix(k, xrange; obsdim=1)
p = heatmap(
    K4;
    yflip=true,
    colorbar=false,
    framestyle=:none,
    background_color=RGBA(0.0, 0.0, 0.0, 0.0),
)
savefig(joinpath(@__DIR__, "src", "assets", "heatmap_prodsum.png"))

plot(heatmap.([K1, K2, K3, K4], yflip=true, colorbar=false)...; layout=(2, 2))
savefig(joinpath(@__DIR__, "src", "assets", "heatmap_combination.png"))

##

for k in [SqExponentialKernel, ExponentialKernel]
    K = kernelmatrix(k(), xrange; obsdim=1)
    v = rand(MvNormal(K + 1e-7I))
    display(plot(xrange, v; lab="", title="f(x)", framestyle=:none))
    savefig(joinpath(@__DIR__, "src", "assets", "GP_sample_$(k).png"))
    display(plot(xrange, kernel.(k(), x₀, xrange); lab="", ylims=(0, 1.1), title="k(0,x)"))
    savefig(joinpath(@__DIR__, "src", "assets", "kappa_function_$(k).png"))
end

for k in [GammaExponentialKernel(1.0, 1.5)]
    sparse = 1
    while !isposdef(kernelmatrix(k, xrange * sparse; obsdim=1) + 1e-5I)
        sparse += 1
    end
    v = rand(MvNormal(kernelmatrix(k, xrange * sparse; obsdim=1) + 1e-7I))
    display(plot(xrange, v; lab="", title="f(x)", framestyle=:none))
    savefig(joinpath(@__DIR__, "src", "assets", "GP_sample_GammaExponentialKernel.png"))
    display(plot(xrange, kernel.(k, x₀, xrange); lab="", ylims=(0, 1.1), title="k(0,x)"))
    savefig(
        joinpath(@__DIR__, "src", "assets", "kappa_function_GammaExponentialKernel.png")
    )
end

for k in [MaternKernel, Matern32Kernel, Matern52Kernel]
    K = kernelmatrix(k(), xrange; obsdim=1)
    v = rand(MvNormal(K + 1e-7I))
    display(plot(xrange, v; lab="", title="f(x)", framestyle=:none))
    savefig(joinpath(@__DIR__, "src", "assets", "GP_sample_$(k).png"))
    display(plot(xrange, kernel.(k(), x₀, xrange); lab="", ylims=(0, 1.1), title="k(0,x)"))
    savefig(joinpath(@__DIR__, "src", "assets", "kappa_function_$(k).png"))
end

for k in [RationalQuadraticKernel]
    K = kernelmatrix(k(), xrange; obsdim=1)
    v = rand(MvNormal(K + 1e-7I))
    display(plot(xrange, v; lab="", title="f(x)", framestyle=:none))
    savefig(joinpath(@__DIR__, "src", "assets", "GP_sample_$(k).png"))
    display(plot(xrange, kernel.(k(), x₀, xrange); lab="", ylims=(0, 1.1), title="k(0,x)"))
    savefig(joinpath(@__DIR__, "src", "assets", "kappa_function_$(k).png"))
end


text:
using .LibClang.CEnum: name_value_pairs

cxname2clname(x::AbstractString) = "CL" * last(split(x, '_'; limit=2))
cxname2clname(x::Symbol) = cxname2clname(string(x))

# each CXCursorKind enum gets a specific Julia type wrapping
# so that we can dispatch directly on node kinds.
abstract type CLCursor end

const CXCursorMap = Dict()

for (sym, val) in name_value_pairs(CXCursorKind)
    clsym = Symbol(cxname2clname(sym))
    @eval begin
        struct $clsym <: CLCursor
            cursor::CXCursor
        end
        CXCursorMap[$sym] = $clsym
    end
end

CLCursor(c::CXCursor) = CXCursorMap[c.kind](c)
Base.convert(::Type{T}, x::CXCursor) where {T<:CLCursor} = CLCursor(x)
Base.cconvert(::Type{CXCursor}, x::CLCursor) = x
Base.unsafe_convert(::Type{CXCursor}, x::CLCursor) = x.cursor

Base.show(io::IO, x::CXCursor) = print(io, "$(kind(x)): $(spelling(x))")
Base.show(io::IO, x::CLCursor) = print(io, "CLCursor (", typeof(x), ") ", name(x))

# each CXTypeKind enum gets a specific Julia type wrapping
# so that we can dispatch directly on node kinds.
abstract type CLType end

const CLTypeMap = Dict()

for (sym, val) in name_value_pairs(CXTypeKind)
    clsym = Symbol(cxname2clname(sym))
    @eval begin
        struct $clsym <: CLType
            type::CXType
        end
        CLTypeMap[$sym] = $clsym
    end
end

CLType(c::CXType) = CLTypeMap[c.kind](c)
Base.convert(::Type{T}, x::CXType) where {T<:CLType} = CLType(x)
Base.cconvert(::Type{CXType}, x::CLType) = x
Base.unsafe_convert(::Type{CXType}, x::CLType) = x.type

Base.show(io::IO, x::CXType) = print(io, "$(kind(x)): $(spelling(x))")
Base.show(io::IO, x::CLType) = print(io, "CLType (", typeof(x), ") ")

# CLToken
abstract type CLToken end

const CLTokenMap = Dict()

for (sym, val) in name_value_pairs(CXTokenKind)
    clsym = Symbol(last(split(string(sym), '_'; limit=2)))
    @eval begin
        struct $clsym <: CLToken
            token::CXToken
            kind::CXTokenKind
            text::String
        end
        CLTokenMap[$sym] = $clsym
    end
end

Base.cconvert(::Type{CXToken}, x::CLToken) = x
Base.unsafe_convert(::Type{CXToken}, x::CLToken) = x.token

Base.show(io::IO, x::CXToken) = print(io, "$(kind(x))")
Base.show(io::IO, x::CLToken) = print(io, typeof(x), "(\"", x.text, "\")")