text:
header("SimpleNet")

@testset "SIMPLENET" begin
    str = """
    <?xml version="1.0"?>
    <pnml xmlns="http://www.pnml.org/version-2009/grammar/pnml">
        <net id="net0" type="continuous">
            <page id="page0">
            <place id="p1"> <initialMarking> <text>1.0</text> </initialMarking> </place>
            <place id="p2"> <initialMarking> <text>2.0</text> </initialMarking> </place>
            <place id="p3">
                <structure att1="doo"/>
                <frog name="hoppy" />
            </place>
            <transition id ="t1"> </transition>
            <transition id ="t2"> </transition>
            <transition id ="t3"> </transition>
            <arc id="a1" source="p1" target="t1"> <inscription/> </arc>
            <arc id="a2" source="p2" target="t1"> <inscription/> </arc>
            <arc id="a3" source="t1" target="p3"> <inscription/> </arc>
            <arc id="a4" source="p3" target="t2"> <inscription/> </arc>
            <arc id="a5" source="t2" target="p1"> <inscription/> </arc>
            <arc id="a6" source="t2" target="p2"> <inscription/> </arc>
            </page>
        </net>
    </pnml>
    """
    @test_call parse_str(str)
    model = @inferred parse_str(str)
    printnode(nets(model))

    @test_call PNML.find_nets(model, :continuous)
    v = @inferred PNML.find_nets(model, :continuous)

    @test_call PNML.first_net(model)
    @test v[begin] == @inferred PnmlNet PNML.first_net(model)
    @test first(v) == @inferred PnmlNet PNML.first_net(model)

    @test_call PNML.SimpleNet(v[begin])
    @test_call PNML.SimpleNet(model)
    @test_call PNML.SimpleNet(PNML.first_net(model))

    net  = @inferred SimpleNet PNML.SimpleNet(v[begin])
    net1 = @inferred SimpleNet PNML.SimpleNet(model)
    net2 = @inferred SimpleNet PNML.SimpleNet(PNML.first_net(model))

    for accessor in [pid, place_ids, transition_ids, arc_ids,
                     reftransition_ids, refplace_ids]
        @test accessor(net1) == accessor(net)
        @test accessor(net2) == accessor(net1)
        @test accessor(net2) == accessor(net)
    end

    for accessor in [places, transitions, arcs]
        for (a,b) in zip(accessor(net1), accessor(net))
            @test pid(a) == pid(b)
        end
        for (a,b) in zip(accessor(net2), accessor(net1))
            @test pid(a) == pid(b)
        end
        for (a,b) in zip(accessor(net2), accessor(net))
            @test pid(a) == pid(b)
        end
    end

    for top in [net, net.net, first(pages(net.net))]
        @test_call places(top)
        for p in @inferred places(top)
            @test @inferred has_place(top, pid(p))
            @test p == @inferred Maybe{Place} place(top, pid(p))
            @test pid(p) ===  p.id
            @test place(top, :bogus) === nothing
            @test typeof(marking(p)) <: typeof(default_marking(nettype(p))())
            @test @inferred(marking(p)) isa typeof(default_marking(nettype(p))())
        end
        @test_call transitions(top)
        for t in @inferred transitions(top)
            @test @inferred has_transition(top, pid(t))
            @test t == @inferred Maybe{Transition} transition(top, pid(t))
            @test pid(t) ===  t.id
            @test transition(top, :bogus) === nothing
            @test condition(t) !== nothing
            @test @inferred condition(t)
        end
        @test_call arcs(top)
        for a in @inferred arcs(top)
            @test @inferred has_arc(top, pid(a))
            @test a == @inferred Maybe{Arc} arc(top, pid(a))
            @test pid(a) ===  a.id
            @test arc(net, :bogus) === nothing
            @test @inferred(PNML.source(a)) !== nothing
            @test @inferred(PNML.target(a)) !== nothing
            @test @inferred(inscription(a)) !== nothing
        end
    end
end

header("RATE")
@testset "rate" begin
    str = """<?xml version="1.0"?>
    <pnml xmlns="http://www.pnml.org/version-2009/grammar/pnml">
        <net id="net0" type="stochastic">
        <page id="page0">
            <transition id ="birth"><rate> <text>0.3</text> </rate> </transition>
        </page>
        </net>
    </pnml>
    """
    model = parse_str(str)
    net = PNML.first_net(model)
    #printnode(net, label="rate net")

    snet = PNML.SimpleNet(net)
    @show snet
    β = PNML.rates(snet)
    @show β
    @test β == LVector(birth=0.3)
end

header("LOTKA-VOLTERRA")
@testset "lotka-volterra" begin
    str = """<?xml version="1.0"?>
    <pnml xmlns="http://www.pnml.org/version-2009/grammar/pnml">
        <net id="net0" type="continuous">
        <page id="page0">
            <place id="wolves">  <initialMarking> <text>10.0</text> </initialMarking> </place>
            <place id="rabbits"> <initialMarking> <text>100.0</text> </initialMarking> </place>
            <transition id ="birth">     <rate> <text>0.3</text> </rate> </transition>
            <transition id ="predation"> <rate> <text>0.015</text> </rate> </transition>
            <transition id ="death">     <rate> <text>0.7</text> </rate> </transition>
            <arc id="a1" source="rabbits"   target="birth"> <inscription><text>1</text> </inscription> </arc>
            <arc id="a2" source="birth"     target="rabbits"> <inscription><text>2</text> </inscription> </arc>
            <arc id="a3" source="wolves"    target="predation"> <inscription><text>1</text> </inscription> </arc>
            <arc id="a4" source="rabbits"   target="predation"> <inscription><text>1</text> </inscription> </arc>
            <arc id="a5" source="predation" target="wolves"> <inscription><text>2</text> </inscription> </arc>
            <arc id="a6" source="wolves"    target="death"> <inscription><text>1</text> </inscription> </arc>
        </page>
        </net>
    </pnml>
    """
    model = parse_str(str)
    net1 = PNML.first_net(model)
    printnode(net1)
    snet = PNML.SimpleNet(net1)

    S = PNML.place_ids(snet) # [:rabbits, :wolves]
    T = PNML.transition_ids(snet)
    PRINT_PNML && @show S, T
    for t in T
        PRINT_PNML && @show PNML.in_out(snet, t)
    end

    # keys are transition ids
    # values are input, output vectors of "tuples" place id -> inscription (integer?)
    Δ = PNML.transition_function(snet)#,T)
    tfun = LVector(
        birth=(LVector(rabbits=1), LVector(rabbits=2)),
        predation=(LVector(wolves=1, rabbits=1), LVector(wolves=2)),
        death=(LVector(wolves=1), LVector()),
    )
    @test Δ.birth     == tfun.birth
    @test Δ.predation == tfun.predation
    @test Δ.death     == tfun.death

    uX = LVector(wolves=10.0, rabbits=100.0) # initialMarking
    u0 = PNML.initialMarking(snet)
    @test u0 == uX
    βx = LVector(birth=0.3, predation=0.015, death=0.7); # transition rate
    β = PNML.rates(snet)
    if PRINT_PNML
        @show Δ
        @show u0
        @show uX
        @show βx
        @show β
    end
    @test β == βx
end


text:
function step(cur_player::AbstractPlayer)
  isnull(cur_player.piece) ||
    return soft_drop(cur_player, true)

  summon_piece!(cur_player) ||
    return end_round(cur_player)

  soft_drop(cur_player, true)
end

ext: jl


text:
const IM2COL_FLOAT_HANDLE = Libdl.dlsym(Native.library, :im2col_float)
const IM2COL_DOUBLE_HANDLE = Libdl.dlsym(Native.library, :im2col_double)
const COL2IM_FLOAT_HANDLE = Libdl.dlsym(Native.library, :col2im_float)
const COL2IM_DOUBLE_HANDLE = Libdl.dlsym(Native.library, :col2im_double)

function im2col{T}(img::Array{T}, n::Int, col::Array{T}, width::Int, height::Int, channels::Int,
    kernel::NTuple{2,Int}, pad::NTuple{2,Int}, stride::NTuple{2,Int})

  kernel_w, kernel_h = kernel
  pad_w, pad_h = pad
  stride_w, stride_h = stride

  img_ptr = pointer(img) + (width*height*channels) * (n-1) * sizeof(T)
  col_ptr = pointer(col)
  im2col_impl(img_ptr, col_ptr, width, height, channels, kernel_w, kernel_h, pad_w, pad_h, stride_w, stride_h)
end

function im2col_impl(img::Ptr{Float32}, col::Ptr{Float32}, width::Int, height::Int, channels::Int,
    kernel_w::Int, kernel_h::Int, pad_w::Int, pad_h::Int, stride_w::Int, stride_h::Int)

  ccall(IM2COL_FLOAT_HANDLE, Void, (Ptr{Float32}, Ptr{Float32}, Cint, Cint, Cint,
        Cint, Cint, # kernel
        Cint, Cint, # pad
        Cint, Cint, # stride
      ), img, col, width, height, channels, kernel_w, kernel_h, pad_w, pad_h, stride_w, stride_h)
end
function im2col_impl(img::Ptr{Float64}, col::Ptr{Float64}, width::Int, height::Int, channels::Int,
    kernel_w::Int, kernel_h::Int, pad_w::Int, pad_h::Int, stride_w::Int, stride_h::Int)

  ccall(IM2COL_DOUBLE_HANDLE, Void, (Ptr{Float64}, Ptr{Float64}, Cint, Cint, Cint,
        Cint, Cint, # kernel
        Cint, Cint, # pad
        Cint, Cint, # stride
      ), img, col, width, height, channels, kernel_w, kernel_h, pad_w, pad_h, stride_w, stride_h)
end

function col2im{T}(col::Array{T}, img::Array{T}, n::Int, img_buf::Array{T}, width::Int, height::Int, channels::Int,
    kernel::NTuple{2,Int}, pad::NTuple{2,Int}, stride::NTuple{2,Int})

  kernel_w, kernel_h = kernel
  pad_w, pad_h = pad
  stride_w, stride_h = stride

  # ignore img_buf
  img_ptr = pointer(img) + (width*height*channels) * (n-1) * sizeof(T)
  col_ptr = pointer(col)
  col2im_impl(col_ptr, img_ptr, width, height, channels, kernel_w, kernel_h, pad_w, pad_h, stride_w, stride_h)
end

function col2im_impl(col::Ptr{Float32}, img::Ptr{Float32}, width::Int, height::Int, channels::Int,
    kernel_w::Int, kernel_h::Int, pad_w::Int, pad_h::Int, stride_w::Int, stride_h::Int)

  ccall(COL2IM_FLOAT_HANDLE, Void, (Ptr{Float32}, Ptr{Float32}, Cint, Cint, Cint,
        Cint, Cint, # kernel
        Cint, Cint, # pad
        Cint, Cint, # stride
      ), col, img, width, height, channels, kernel_w, kernel_h, pad_w, pad_h, stride_w, stride_h)
end
function col2im_impl(col::Ptr{Float64}, img::Ptr{Float64}, width::Int, height::Int, channels::Int,
    kernel_w::Int, kernel_h::Int, pad_w::Int, pad_h::Int, stride_w::Int, stride_h::Int)

  ccall(COL2IM_DOUBLE_HANDLE, Void, (Ptr{Float64}, Ptr{Float64}, Cint, Cint, Cint,
        Cint, Cint, # kernel
        Cint, Cint, # pad
        Cint, Cint, # stride
      ), col, img, width, height, channels, kernel_w, kernel_h, pad_w, pad_h, stride_w, stride_h)
end


text:
using Test
using MixedDistributions
import Distributions
const Dst = Distributions
using Statistics: mean, var

const u1 = Dst.Uniform()
const u2 = Dst.Uniform(2.0,2.5)
const m = MixedDistribution([0.1,0.05,0.05],[1.,7.,3.],[0.2,0.6], [u1, u2])
const m1 = MixedDistribution(Float64[],Float64[],[0.4,0.6], [u1, u2])

@testset "Base interface" begin
    @test maximum(m1) ≈ maximum(u2)
    @test minimum(m1) ≈ minimum(u1)
    @test mean(m1)    ≈ 0.4 * mean(u1) + 0.6 * mean(u2)
    @test maximum(m)  ≈ 7.
    @test minimum(m)  ≈ 0.
end

@testset "Mixed properties" begin
    @test issorted(m.mass_points)
    @test mean(m) ≈ 0.2 * mean(u1) + 0.6 * mean(u2) + sum([0.1,0.05,0.05].*[1.,7.,3.])
end

@testset "Random generation" begin
	# each range is suppose to yield the corresponding mass point
	r1 = [0.0, 0.05,0.099, 0.1]
	r2 = [0.1025, 0.105,0.11,0.14]
	r3 = [0.151,0.16,0.19,0.2]
	# each range yields corresponding points on the uniform domains
	r4 = 0.2001:0.01:0.4
	r5 = 0.4001:0.01:1.0
	@test all(MixedDistributions.ordered_rand(m,r) ≈ 1.0 for r in r1)
	@test all(MixedDistributions.ordered_rand(m,r) ≈ 3.0 for r in r2)
	@test all(MixedDistributions.ordered_rand(m,r) ≈ 7.0 for r in r3)
	@test all(MixedDistributions.ordered_rand(m,r) ≈ ((r-0.2)/0.2) for r in r4)
	@test all(MixedDistributions.ordered_rand(m,r) ≈ ((r-0.4)/0.6)*0.5+2.0 for r in r5)
end

@testset "Homogeneous continuous behave like MixtureModel" begin
	mixture = Dst.MixtureModel([u1,u2],[0.4,0.6])
	@test mean(mixture) ≈ mean(m1)
	@test var(mixture)  ≈  var(m1)
end

@testset "Exception behaves like one" begin
	ex = MixedDistributions.CDFException(Dst.Normal())
	t = typeof(ex)
	@test t <: Exception
	@test t <: MixedDistributions.CDFException{<:Dst.Normal}
end

@testset "Simple partial expectancy" begin
	unif = MixedDistribution(Float64[],Float64[],[1.0], [u1])
	@test MixedDistributions.partial_expectancy(unif; dx = 0.0001) ≈ mean(unif)
	@test abs(MixedDistributions.partial_expectancy(m) - mean(m)) <= 0.02 * mean(m)
	@test abs(MixedDistributions.partial_expectancy(m1) - mean(m1)) <= 0.02 * mean(m1)
	ndst = Dst.Normal(1.5, 0.6)
	normal_mass = MixedDistribution([0.25], [1.5], [0.75], [ndst])
	@test mean(normal_mass) ≈ mean(ndst)
end

@testset "Interval partial expectancy" begin
	md = MixedDistribution([0.5],[0.5],[0.5], [u1])
	@test MixedDistributions.partial_expectancy(md, dx = 0.0001) ≈ mean(md)
	for δ in [0.01,0.05,0.06,0.1,0.45]
		(intervals, ys) = MixedDistributions.partial_expectancy(md, [0.5-δ, 0.5+δ], dx = 0.0001)
		@test length(ys) == 1
		@test length(intervals) == 1
		(low, high) = intervals[1]
		@test low ≈ 0.5 - δ
		@test high ≈ 0.5 + δ
		@test ys[1] ≈ (δ + 0.5) * 0.5
	end
end


text:
#= This file is a part of URIValidator. License is MIT:
https://github.com/rafpyprog/URIValidator.jl/blob/master/LICENSE.md =#

include("components.jl")

function splitnetloc(url, start=0)
    component_delims = "/?#"
    delim = length(url) # position of end of domain part of url, default is end
    has_delimiter = false
    for c in component_delims  # look for delimiters; the order is NOT important
        wdelim = searchindex(url, c, start)  # find first of this delim
        if wdelim > 0  # if found
            delim = min(delim, wdelim) # use earliest delim position
            has_delimiter = true
        end
    end

    if (has_delimiter == true) & (delim == length(url))
        rest = url[end]
        domain = url[start:end-1]
    elseif has_delimiter == true
        rest = url[delim:end]
        domain = url[start:delim-1]
    else
        domain = url[start:end]
        rest = ""
    end
    return domain, rest
end

function has_delimiter(url, delimiter)
    if isa(url, AbstractString)
        check = contains(url, delimiter)
    elseif isa(url, Char)
        check = delimiter in url
    else
        println("Invalid url type.")
        throw(DomainError())
    end
    return check
end

function URISplit(url, allow_fragments=true)
    #= Parse a URL into 5 components:
    <scheme>://<netloc>/<path>?<query>#<fragment> =#
    netloc = query = fragment = ""

    i = searchindex(url, ":")
    if i > 0
        if url[1:i-1] == "http" # optimize the common case
            scheme = lowercase(url[1:i-1])
            url = url[i+1:end]
            if url[1:2] == "//"
                netloc, url = splitnetloc(url, 3)
                if ((contains(netloc, "[")  & (!contains(netloc, "]"))) | (contains(netloc, "]")  & (!contains(netloc, "["))))
                    println("Invalid IPv6 URL")
                    throw(DomainError())
                end
            end

            if allow_fragments & has_delimiter(url, "#")
                url, fragment = split(url, "#", limit=2)
            end

            if has_delimiter(url, "?")
                url, query = split(url, "?", limit=2)
            end
            splited_URI = Dict("scheme"=>scheme, "authority"=>netloc, "path"=>url, "query"=>query, "fragment"=>fragment)
            components = URIComponents(splited_URI)
            return components
        else
            #= # make sure "url" is not actually a port number (in which case
            # "scheme" is really part of the path) =#
            rest = url[i + 1:end]
            if !isdefined(:rest) | any(!(c in "0123456789") for c in rest)
                # not a port number
                scheme, url = lowercase(url[1:i-1]), rest
            end
        end
    end

    if url[1:2] == "//"
        netloc, url = splitnetloc(url, 3)
        if ((contains(netloc, "[")  & (!contains(netloc, "]"))) | (contains(netloc, "]")  & (!contains(netloc, "["))))
            println("Invalid IPv6 URL")
            throw(DomainError())
        end
    end

    if allow_fragments & has_delimiter(url, "#")
        url, fragment = split(url, "#", limit=2)
    end

    if has_delimiter(url, "?")
        url, query = split(url, "?", limit=2)
    end
    splited_URI = Dict("scheme"=>scheme, "authority"=>netloc, "path"=>url, "query"=>query, "fragment"=>fragment)
    components = URIComponents(splited_URI)
    return components
end


text:
module BenchArray

using Benchmarkplotting
using Statistics
using Gadfly
using MLStyle
using DataFrames
import Match

import Base.getindex
getindex(asoc_lst :: Vector{Pair{Symbol, T}}, key ::Symbol) where T =
    for (search_key, value) in asoc_lst
        if search_key === key
            return value
        end
    end

"""
3 cases could succeed in validation:

[_, _, frag..., 10], sum(frag) > 10
[1, 2, 3, _...]
[_, _, 1, _, 2, 3, z], z > 10

"""

data = [
    Symbol("ok_case1") => [42, 42, 3, 3, 5, 10],
    Symbol("ok_case2") => [1, 2, 3, 42, 42, 42],
    Symbol("ok_case3") => [42, 42, 1, 42, 2, 3, 11],
    Symbol("early_fail") => [42, 42, 3, 3, 3, 10],
    Symbol("late_fail1") => [42, 42, 3, 3, 5, 42],
    Symbol("late_fail2") => [1, 2, 42],
    Symbol("late_fail3") => [42, 42, 1, 42, 2, 3, -42]
]

implementations = [
    :MLStyle => (@λ begin
        [_, _, (frag && if sum(frag) > 10 end)..., 10] -> 1
        [1, 2, 3, _...]  -> 2
        [_, _, 1, _, 2, 3, z && if z > 10 end]  -> 3
        _ -> 4
    end),
    Symbol("Match.jl") => function (x)
        Match.@match x begin
            [_, _, frag..., 10], if sum(frag) > 10 end => 1
            [1, 2, 3, _...]  => 2
            [_, _, 1, _, 2, 3, z], if z > 10 end => 3
            _ => 4
        end
    end,
    :HandWritten => function(vec)
        !(vec isa Vector)  ? 4 :
        let n = length(vec)
            if n > 1 && vec[end] === 10 && sum(vec[3:end-1]) > 10
                1
            elseif n > 3 && vec[1] === 1 && vec[2] === 2 && vec[3] === 3
                2
            elseif n === 7 && vec[3] === 1 && vec[5] === 2 && vec[6] === 3 &&
                    vec[end] > 10
                3
            else
                4
            end
        end
    end
]

criterion(x) = (meantime = mean(x.times), allocs = 1 + x.allocs)
df = bcompare(criterion, data, implementations, repeat=1)

@info df

theme = Theme(
    guide_title_position = :left,
    colorkey_swatch_shape = :circle,
    minor_label_font = "Consolas",
    major_label_font = "Consolas",
    point_size=5px
)
report_meantime, df_time = report(:meantime, df, Scale.y_log10, theme)
report_allocs, df_allocs = report(:allocs, df, theme)

open("stats/vs-hw(array).txt", "w") do f
    write(f, string(df))
end

draw(SVG("stats/vs-hw(array)-on-time.svg", 10inch, 4inch), report_meantime);
draw(SVG("stats/vs-hw(array)-on-allocs.svg", 10inch, 4inch), report_allocs);

end


text:
# Jackknife resampling.

"""
    jackknife(f, xss...)

Perform jackknife resampling for the function `f` with the data `xss`.

This uses the formulas from Peter Young's notes
(https://arxiv.org/abs/1210.3781, https://arxiv.org/pdf/1210.3781v3.pdf).
"""
function jackknife(f, xss...)
    length(xss) >= 1 || throw(DomainError(length(xss), "At least one data set."))

    N = length(xss[1])
    lengths = [length(xs) for xs in xss]
    all(lengths .== N) || throw(DomainError(lengths, "Uniform length."))

    # Regular averages.
    xss_mean = Float64[]
    for xs in xss
        push!(xss_mean, mean(xs))
    end

    # "All-but-one" averages.
    xss_J = Array{Float64,1}[]
    for xs in xss
        push!(xss_J, (sum(xs) .- xs) / (N - 1))
    end

    f_plain = f(xss_mean...)
    f_Js = f(xss_J...)

    # Unbiased mean.
    f_J = N * f_plain - (N - 1) * mean(f_Js)
    # Standard error.
    f_err = sqrt(N - 1) * std(f_Js; corrected=false)

    f_J, f_err
end


text:
using LinearAlgebra

@testset "PointForce" begin

  Re = 200
  xlim, ylim = (-2,2), (-2,4)
  Δx,Δt = setstepsizes(Re,gridRe=4)

  sys = NavierStokes(Re,Δx,xlim,ylim,Δt)
  w₀ = Nodes(Dual,size(sys))

  x0 = [1.5,0.0]
  f0 = 10.0
  t0 = 4
  σt = 1
  pulse = PointForce(w₀,x0,f0,t0,σt,sys)

  @test isapprox(abs(sum(pulse(t0))-f0),0,atol=1e-2)

  x0 = (1.5,0.0)
  pulse = PointForce(w₀,x0,f0,t0,σt,sys)

  @test isapprox(abs(sum(pulse(t0))-f0),0,atol=1e-2)

end

@testset "Spatial pulse" begin

σx = 0.5
σy = 0.5
x0 = 0
y0 = 0
A = 1
dgaussx = SpatialGaussian(σx,σy,x0,y0,A,deriv=1)
dgaussy = SpatialGaussian(σx,σy,x0,y0,A,deriv=2)
x = 0.1
y = 0.2
@test dgaussx(0.1,0.2) == dgaussy(0.2,0.1) ≈ -2*A*x/π/σx^3/σy*exp(-x^2/σx^2)*exp(-y^2/σy^2)

end


text:
@testset "arrays can be diffed" begin
    a1 = [1, 2, 3, 4]

    # one number changed
    a2 = [1, 7, 3, 4]
    d = deepdiff(a1, a2)
    @test added(d) == [2]
    @test removed(d) == [2]
    @test changed(d) == []
    @test before(d) == a1
    @test after(d) == a2

    # removed from middle
    d = deepdiff(a1, [1, 3, 4])
    @test removed(d) == [2]
    @test added(d) == []

    # removed from beginning
    d = deepdiff(a1, [2, 3, 4])
    @test removed(d) == [1]
    @test added(d) == []

    # removed from end
    d = deepdiff(a1, [1, 2, 3])
    @test removed(d) == [4]
    @test added(d) == []

    # added to end
    d = deepdiff(a1, [1, 2, 3, 4, 5])
    @test removed(d) == []
    @test added(d) == [5]

    # added to beginning
    d = deepdiff(a1, [0, 1, 2, 3, 4])
    @test removed(d) == []
    @test added(d) == [1]

    # two additions
    d = deepdiff(a1, [1, 4, 2, 5, 3, 4])
    @test removed(d) == []
    @test added(d) == [2, 4]

    # two removals
    d = deepdiff(a1, [2, 4])
    @test removed(d) == [1, 3]
    @test added(d) == []

    d = deepdiff([NaN], [NaN])
    @test removed(d) == added(d) == []

    d = deepdiff([missing], [missing])
    @test removed(d) == added(d) == []

    d = deepdiff([NaN], [])
    @test d == d

    d = deepdiff([0.0], [-0.0])
    @test removed(d) == [1]
    @test added(d) == [1]
end

max_stars_repo_name: ssfrr/StructureDiffs

max_issues_repo_name: ssfrr/StructureDiffs

max_forks_repo_name: ssfrr/StructureDiffs


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct IoK8sKubeAggregatorPkgApisApiregistrationV1beta1APIServiceSpec <: SwaggerModel
    caBundle::Any # spec type: Union{ Nothing, Vector{UInt8} } # spec name: caBundle
    group::Any # spec type: Union{ Nothing, String } # spec name: group
    groupPriorityMinimum::Any # spec type: Union{ Nothing, Int32 } # spec name: groupPriorityMinimum
    insecureSkipTLSVerify::Any # spec type: Union{ Nothing, Bool } # spec name: insecureSkipTLSVerify
    service::Any # spec type: Union{ Nothing, IoK8sKubeAggregatorPkgApisApiregistrationV1beta1ServiceReference } # spec name: service
    version::Any # spec type: Union{ Nothing, String } # spec name: version
    versionPriority::Any # spec type: Union{ Nothing, Int32 } # spec name: versionPriority

    function IoK8sKubeAggregatorPkgApisApiregistrationV1beta1APIServiceSpec(;caBundle=nothing, group=nothing, groupPriorityMinimum=nothing, insecureSkipTLSVerify=nothing, service=nothing, version=nothing, versionPriority=nothing)
        o = new()
        validate_property(IoK8sKubeAggregatorPkgApisApiregistrationV1beta1APIServiceSpec, Symbol("caBundle"), caBundle)
        setfield!(o, Symbol("caBundle"), caBundle)
        validate_property(IoK8sKubeAggregatorPkgApisApiregistrationV1beta1APIServiceSpec, Symbol("group"), group)
        setfield!(o, Symbol("group"), group)
        validate_property(IoK8sKubeAggregatorPkgApisApiregistrationV1beta1APIServiceSpec, Symbol("groupPriorityMinimum"), groupPriorityMinimum)
        setfield!(o, Symbol("groupPriorityMinimum"), groupPriorityMinimum)
        validate_property(IoK8sKubeAggregatorPkgApisApiregistrationV1beta1APIServiceSpec, Symbol("insecureSkipTLSVerify"), insecureSkipTLSVerify)
        setfield!(o, Symbol("insecureSkipTLSVerify"), insecureSkipTLSVerify)
        validate_property(IoK8sKubeAggregatorPkgApisApiregistrationV1beta1APIServiceSpec, Symbol("service"), service)
        setfield!(o, Symbol("service"), service)
        validate_property(IoK8sKubeAggregatorPkgApisApiregistrationV1beta1APIServiceSpec, Symbol("version"), version)
        setfield!(o, Symbol("version"), version)
        validate_property(IoK8sKubeAggregatorPkgApisApiregistrationV1beta1APIServiceSpec, Symbol("versionPriority"), versionPriority)
        setfield!(o, Symbol("versionPriority"), versionPriority)
        o
    end
end # type IoK8sKubeAggregatorPkgApisApiregistrationV1beta1APIServiceSpec

const _property_map_IoK8sKubeAggregatorPkgApisApiregistrationV1beta1APIServiceSpec = Dict{Symbol,Symbol}(Symbol("caBundle")=>Symbol("caBundle"), Symbol("group")=>Symbol("group"), Symbol("groupPriorityMinimum")=>Symbol("groupPriorityMinimum"), Symbol("insecureSkipTLSVerify")=>Symbol("insecureSkipTLSVerify"), Symbol("service")=>Symbol("service"), Symbol("version")=>Symbol("version"), Symbol("versionPriority")=>Symbol("versionPriority"))
const _property_types_IoK8sKubeAggregatorPkgApisApiregistrationV1beta1APIServiceSpec = Dict{Symbol,String}(Symbol("caBundle")=>"Vector{UInt8}", Symbol("group")=>"String", Symbol("groupPriorityMinimum")=>"Int32", Symbol("insecureSkipTLSVerify")=>"Bool", Symbol("service")=>"IoK8sKubeAggregatorPkgApisApiregistrationV1beta1ServiceReference", Symbol("version")=>"String", Symbol("versionPriority")=>"Int32")
Base.propertynames(::Type{ IoK8sKubeAggregatorPkgApisApiregistrationV1beta1APIServiceSpec }) = collect(keys(_property_map_IoK8sKubeAggregatorPkgApisApiregistrationV1beta1APIServiceSpec))
Swagger.property_type(::Type{ IoK8sKubeAggregatorPkgApisApiregistrationV1beta1APIServiceSpec }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_IoK8sKubeAggregatorPkgApisApiregistrationV1beta1APIServiceSpec[name]))}
Swagger.field_name(::Type{ IoK8sKubeAggregatorPkgApisApiregistrationV1beta1APIServiceSpec }, property_name::Symbol) =  _property_map_IoK8sKubeAggregatorPkgApisApiregistrationV1beta1APIServiceSpec[property_name]

function check_required(o::IoK8sKubeAggregatorPkgApisApiregistrationV1beta1APIServiceSpec)
    (getproperty(o, Symbol("groupPriorityMinimum")) === nothing) && (return false)
    (getproperty(o, Symbol("service")) === nothing) && (return false)
    (getproperty(o, Symbol("versionPriority")) === nothing) && (return false)
    true
end

function validate_property(::Type{ IoK8sKubeAggregatorPkgApisApiregistrationV1beta1APIServiceSpec }, name::Symbol, val)
    if name === Symbol("caBundle")
    end
end


text:
using CSV
using DataFrames

csv_reader = CSV.File("Datasets/monk1/monks-1.csv")
println(typeof(csv_reader))

df = DataFrame(CSV.File("Datasets/monk1/monks-1.csv",header=0))
select!(df, Not(:Column1))
select!(df, Not(:Column9))
classes=Matrix(select(df,:Column2))
select!(df, Not(:Column2))
df=Matrix(df)
unique(df) .== permutedims(df)


text:
include("../src/main.jl")
include("test_save.jl")


text:
"""
Analytic ITER-like Solov'ev equilibrium with X-point.
"""
module TokamakIterCylindrical

    import ElectromagneticFields.AxisymmetricTokamakCylindrical

    export initial_conditions_barely_passing, initial_conditions_barely_trapped,
           initial_conditions_deeply_passing, initial_conditions_deeply_trapped,
           initial_conditions_trapped

    export hamiltonian, toroidal_momentum

    AxisymmetricTokamakCylindrical.@code_iter() # inject magnetic field code

    const qᵢ = [7.0-1.4, 0.0, 0.0, 2.8166280889939737]
    const parameters = (μ = 4.607782183567846,)

    include("guiding_center_4d_common.jl")
    include("guiding_center_4d_equations.jl")

    function toroidal_momentum(t,q)
        R(t,q) * ϑ₃(t,q)
    end

    include("guiding_center_4d_diagnostics.jl")

    initial_conditions_barely_passing() = ([2.5, 0., 0., 3.425E-1], (μ = 1E-2,))
    initial_conditions_barely_trapped() = ([2.5, 0., 0., 3.375E-1], (μ = 1E-2,))
    initial_conditions_deeply_passing() = ([2.5, 0., 0.,  5E-1],    (μ = 1E-2,))
    initial_conditions_deeply_trapped() = ([2.5, 0., 0.,  1E-1],    (μ = 1E-2,))
    initial_conditions_trapped()        = ([7.0, 0., 0., -2E-3],    (μ = 1.88E-7,))

end


text:
include("main.jl")

REWRITE_DATA = false

BUSES, LINES, GENERATORS, rPTDF = load_case_data()
ERROR_VARIANCES_TRUE, ERROR_DIST_TRUE, ERROR_VARIANCES_WC, ERROR_VARIANCES_SAMPLE = init_stochasticity()

# CASE: Numerical Experiments on the difference between SOC and LinDist
# 1) Voltage Deviations
# 2) Sum of losses and respective change in objective
# => With different demands

casename = "lindist_vs_soc"
case_path = "results/$(casename)"

if REWRITE_DATA && isdir(case_path)
    rm(case_path, recursive=true)
end

if !(isdir(case_path)) || REWRITE_DATA
    # Only run if there are no saved results
    mkdir(case_path)

    demand_cases = collect(0.5:0.1:1.3)

    voltage_result_soc = DataFrame()
    voltage_result_ld = DataFrame()
    voltage_deltas = DataFrame()
    objective_result = DataFrame(case=["SOC", "LinDist", "delta"])

    for α in demand_cases
        buses_part_load = change_load_same_pf(α)

        soc_objective, soc_bus_results, soc_line_results = Solve_Dist_SOC(buses_part_load, LINES, GENERATORS)
        ld_objective, ld_bus_results, ld_line_results = Solve_LinDist(buses_part_load, LINES, GENERATORS)

        voltage_result_soc[Symbol(α)] = sqrt.(soc_bus_results[:v_squared])
        voltage_result_ld[Symbol(α)] = sqrt.(ld_bus_results[:v_squared])
        voltage_deltas[Symbol(α)] = (voltage_result_ld[Symbol(α)] -  voltage_result_soc[Symbol(α)]) ./ voltage_result_soc[Symbol(α)]

        objective_result[Symbol(α)] = [soc_objective, ld_objective, (ld_objective - soc_objective) / soc_objective]
    end

    bus_summary = DataFrame(bus=Any[], quant_10=Any[], quant_90=Any[], median=Any[])
    for i in 1:nrow(voltage_deltas)
        arr = [voltage_deltas[i,a] for a in names(voltage_deltas)]
        b_sum = [i, exp_arr_quantile(arr, 0.1), exp_arr_quantile(arr, 0.9), median(arr)]
        push!(bus_summary, b_sum)
    end

    writetable("$(case_path)/voltage_result_soc.csv", voltage_result_soc)
    writetable("$(case_path)/voltage_result_ld.csv", voltage_result_ld)
    writetable("$(case_path)/voltage_deltas.csv", voltage_deltas)
    writetable("$(case_path)/bus_summary.csv", bus_summary)
    writetable("$(case_path)/objective_result.csv", objective_result)

end


text:
@testset "MaxFuncEvals" begin

# It was reported in an GitHub Issue that MaxFuncEvals doesn't have an effect
# https://github.com/robertfeldt/BlackBoxOptim.jl/issues/157
# so we added a test to check this.

function rosenbrock2d(x)
    return (1.0 - x[1])^2 + 100.0 * (x[2] - x[1]^2)^2
end

mutable struct MyCounter
    count::Int
end

function myfit(c::MyCounter, fitness::Function)
    f(x) = begin
        c.count += 1
        fitness(x)
    end
    return f
end

MFE = 10
c = MyCounter(0)
res = bboptimize(myfit(c, rosenbrock2d); SearchRange = (-5.0, 5.0), NumDimensions = 2, MaxFuncEvals = MFE);
@test BlackBoxOptim.stop_reason(res) == "Max number of function evaluations ($MFE) reached"
@test MFE <= c.count <= (MFE+20) # Must be at least 10 but can be somewhat higher since might be several evaluations per step of the optimizer

MFE = 1042
c = MyCounter(0)
res = bboptimize(myfit(c, rosenbrock2d); SearchRange = (-5.0, 5.0), NumDimensions = 2, MaxFuncEvals = MFE);
@test BlackBoxOptim.stop_reason(res) == "Max number of function evaluations ($MFE) reached"
@test MFE <= c.count <= (MFE+20) # Must be at least 10 but can be somewhat higher since might be several evaluations per step of the optimizer

end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

module DataLakeStoreFileSystemManagementClient

using Random
using Dates
using Swagger
import Swagger: field_name, property_type, hasproperty, validate_property, SwaggerApi, SwaggerModel
import Base: convert, propertynames

include("modelincludes.jl")

include("api_FileSystemApi.jl")

# export models
export convert, AclStatus
export convert, AclStatusResult
export convert, AdlsAccessControlException
export convert, AdlsBadOffsetException
export convert, AdlsError
export convert, AdlsFileAlreadyExistsException
export convert, AdlsFileNotFoundException
export convert, AdlsIOException
export convert, AdlsIllegalArgumentException
export convert, AdlsRemoteException
export convert, AdlsRuntimeException
export convert, AdlsSecurityException
export convert, AdlsThrottledException
export convert, AdlsUnsupportedOperationException
export convert, ContentSummary
export convert, ContentSummaryResult
export convert, FileOperationResult
export convert, FileStatusProperties
export convert, FileStatusResult
export convert, FileStatuses
export convert, FileStatusesResult

# export operations
export convert, FileSystemApi

export check_required, field_name, property_type, hasproperty, propertynames, validate_property, convert

end


text:
"""
white_kernel function. Requires 0 hyperparameters.
Use with include("src/kernels/white_kernel.jl").
hyperparameters == []
"""
function white_kernel(
    hyperparameters::Vector{<:Real},
    δ::Real;
    dorder::Vector{<:Integer}=zeros(Int64, length(hyperparameters) + 2))

    if all(dorder .== 0)
        return Int(δ == 0)
    end

    return 0

end

return white_kernel, 0  # the function handle and the number of kernel hyperparameters

max_stars_repo_name: christiangil/GPLinearODEMaker

max_issues_repo_name: christiangil/GPLinearODEMaker

max_forks_repo_name: christiangil/GPLinearODEMaker


text:
include("path.jl")
include("model.jl")
using TrajectoryOptimization
using TrajOptPlots
import TrajectoryOptimization: set_state!
using Plots

function gen_car_prob(U0, path, N; s0=0.0, L=10.0, δ0=0.0, integration=RK3)

	# Discretization
	sf = L + s0
	ds = sf/(N-1)

	s_traj = range(0,sf,step=ds)
	@assert N == length(s_traj)
	# k_traj = k_itp.(s_traj)
	# d_path = Dict(s_traj.=>k_traj)

	# Model
	car = BicycleCar(path=path)
	n,m = size(car)
	p_c = 0.1 # probability of contingency plan

	x,u = zeros(car)
	car.μ = 0.9
	fiala(a) = fiala_tire_model(car.μ, car.Cαf, a, 0.0, 10000)
	plot(fiala.(range(-0.2,0.2, length=100)))

	# Objective
	Ux_des = 15.
	x0 = @SVector [δ0,     # δ
	               0,      # fx
	               0,      # r
	               0,      # Uy
	               Ux_des, # Ux
	               0,      # dpsi
	               0,      # e
	               0]      # t

	xd = @SVector [0,      # δ
	               0,      # fx
	               0,      # r
	               0,      # Uy
	               Ux_des, # Ux
	               0,      # dpsi
	               0,      # e
	               0]      # t

	Qd = @SVector [0.001,      # δ
	               0.001,      # fx
	               0.00,      # r
	               0.00,      # Uy
	               1.0,      # Ux
	               100.0,      # dpsi
	               100.0,      # e
	               0.0]      # t

	Qf = @SVector [0.000,      # δ
	               0.000,      # fx
	               0.000,      # r
	               0.000,      # Uy
	               10.0,      # Ux
	               100.0,      # dpsi
	               100.0,      # e
	               0.00]      # t

	Rd   = @SVector [0.1,  # δ_dot
	                 0.1]    # fx_dot
	Rd_c = @SVector [0.3282,  # δ_dot
	                 4e-8]    # fx_dot

	# Q = Diagonal([(1-p_c)*Qd; p_c*Qd])
	# R = Diagonal([Rd; Rd_c])

	# Qd = @SVector fill(0.1,n)
	# Rd = @SVector fill(0.01,m)
	Q = Diagonal(Qd)
	R = Diagonal(Rd)

	obj = LQRObjective(Q,R,Q,xd,N)

	# Bound Constraints
	δ_dot_bound = deg2rad(90)  # deg/s
	δ_bound = deg2rad(27)  # deg
	Fx_max = car.μ*car.mass*car.g
	Ux_min = 1.0  # m/s
	e_bnd = 1.0 # m

	x_min = fill(-Inf,n)
	x_min[1] = -δ_bound
	x_min[5] = Ux_min
	x_min[7] = -e_bnd

	x_max = fill(Inf,n)
	x_max[1] = δ_bound
	x_max[2] = Fx_max
	x_max[7] = e_bnd

	u_min = [-δ_dot_bound; -Inf]
	u_max = [ δ_dot_bound;  Inf]

	bnd = BoundConstraint(n,m, x_max=x_max, x_min=x_min, u_min=u_min, u_max=u_max)

	# Brake constraint
	brake = BrakeForceConstraint(car)

	# Constraint Set
	conSet = ConstraintSet(n,m,N)
	add_constraint!(conSet, bnd, 1:N-1)
	add_constraint!(conSet, brake, 1:N)

	# Problem
	prob = Problem(car, obj, xd, sf, constraints=conSet, x0=x0, t0=s0,
		integration=integration)
	initial_controls!(prob, U0)

	return prob
end


text:
module SCIDACRoutines

import QJuliaInterface
import QJuliaEnums

################# Interface functions, found in quda.h #######################
#be sure that a directory with libscidacwrap.so is in LD_LIBRARY_PATH

#QIO routines (based on QUDA implementation):
   read_gauge_field_qj(filename, gauge, prec, X, argc, argv) = ccall((:read_gauge_field, "libscidacwrap"), Cvoid, (Ptr{Cchar}, Ptr{Cvoid}, Cint, Ptr{Cint}, Cint, Ptr{Ptr{Cchar}}), filename, gauge, prec, X, argc, argv)

#QMP routines:
   QMPInitComms_qj(argc, argv, commDims) = ccall((:QMPInitComms, "libscidacwrap"), Cvoid, (Cint, Ptr{Ptr{Cchar}}, Ptr{Cint}), argc, argv, commDims)

   QMPFinalizeComms_qj() = ccall((:QMPFinalizeComms, "libscidacwrap"), Cvoid, (), )

end #QUDARoutines

max_stars_repo_name: alexstrel/QJulia

max_issues_repo_name: alexstrel/QJulia

max_forks_repo_name: alexstrel/QJulia


text:
"""
    ZeroInflatedDistribution{Te, Tp} <: Distributions.Sampleable{Univariate,Continuous}

Inner constructor:
    - encdist::Bernoulli : Bernoulli distribution with probability of encounter
    - posdist::Distributions.ContinuousUnivariateDistribution : Continuous
            distribution with distribution of observations conditional on an
            encounter

Outer constructor:
    - zil::ZeroInflatedLink : Link function
    - posdist::Type{ContinuousUnivariateDistribution} : Type of the distribution
        of observations conditional on an encounter
    - p1 : First linear predictor
    - p2 : Second linear predictor
    - disp : Positive observation dispersion parameter (standard deviation)
    - biascorrect::Bool : For log-normal positive observations, parameterize by
        the mean rather than the median

A distribution with positive probability of a zero observation as well as a
continuous response. Generally determined by two processes and a relevant link
function. Outer constructors are currently implemented for the following
distributions:

    - LogNormal
    - Gamma
    - InverseGamma
    - InverseGaussian

The dispersion parameter is the usual standard deviation on the log scale for
log-normal, standard deviation for the gamma and inverse gamma, and the shape
parameter for the inverse Gaussian distribution.
"""
struct ZeroInflatedDistribution{Te, Tp} <: Distributions.Sampleable{Univariate,Continuous}
    encdist::Te
    posdist::Tp

    function ZeroInflatedDistribution(encdist::Te, posdist::Tp) where
            {Te<:Distributions.Bernoulli,
             Tp<:Distributions.ContinuousUnivariateDistribution}
        new{Te, Tp}(encdist, posdist)
    end
end

function ZeroInflatedDistribution(
    zil::AbstractZeroInflatedLink, posdist::Type{LogNormal}, p1, p2, disp;
    biascorrect = true)
    p = encprob(zil, p1, p2)
    if biascorrect
        bias = disp^2 / 2
    else
        bias = zero(disp)
    end
    r = posrate(zil, p1, p2)
    ZeroInflatedDistribution(Bernoulli(p), posdist(log(r) - bias, disp))
end

function ZeroInflatedDistribution(
    zil::AbstractZeroInflatedLink, posdist::Type{Gamma}, p1, p2, disp)
    p = encprob(zil, p1, p2)
    mu = posrate(zil, p1, p2)
    alpha = (mu / disp)^2
    beta = disp^2 / mu
    ZeroInflatedDistribution(Bernoulli(p), posdist(alpha, beta))
end

function ZeroInflatedDistribution(
    zil::AbstractZeroInflatedLink, posdist::Type{InverseGamma}, p1, p2, disp)
    p = encprob(zil, p1, p2)
    # From mean and standard deviation to alpha and beta:
    # μ = β / (α - 1) for α > 1
    # σ² = β² / [(α - 1)²(α - 2)] for α > 2
    # β = μ(α - 1)
    # σ² = μ² / (α - 2)
    # α = μ² / σ² + 2
    # β = μ(μ² / σ² + 1)
    # β = μ³ / σ² + μ
    mu = posrate(zil, p1, p2)
    alpha = (mu / disp)^2 + 2
    beta = mu * (alpha - 1)
    ZeroInflatedDistribution(Bernoulli(p), posdist(alpha, beta))
end

function ZeroInflatedDistribution(
    zil::AbstractZeroInflatedLink, posdist::Type{InverseGaussian}, p1, p2, disp)
    p = encprob(zil, p1, p2)
    mu = posrate(zil, p1, p2)
    ZeroInflatedDistribution(Bernoulli(p), posdist(mu, disp))
end

# Density functions: pdf, logpdf
function Distributions.pdf(zil::ZeroInflatedDistribution, obs)
    dens = pdf(zil.encdist, obs ≠ 0)
    if obs ≠ 0
        dens *= pdf(zil.posdist, obs)
    end
    dens
end
function Distributions.logpdf(zil::ZeroInflatedDistribution, obs)
    logdens = logpdf(zil.encdist, obs ≠ 0)
    if obs ≠ 0
        logdens += logpdf(zil.posdist, obs)
    end
    logdens
end

# CDF functions: cdf
function Distributions.cdf(zil::ZeroInflatedDistribution, obs)
    # Zero has positive probability
    cum = failprob(zil.encdist)
    if obs > 0
        # Then add CDF of positive distribution rescaled to probability of
        # encounter
        cum += cdf(zil.posdist, obs) * succprob(zil.encdist)
    end
    cum
end

# Quantile functions: quantile
function Distributions.quantile(zil::ZeroInflatedDistribution, p)
    p0 = failprob(zil.encdist)
    if p ≤ p0
        q = zero(p)
    else
        p2 = (p - p0) / succprob(zil.encdist)
        q = quantile(zil.posdist, p2)
    end
    q
end

# FIXME Consider handline Bernoulli(0), Bernoulli(1) corner cases
# Support functions: minimum, maximum, insupport
function Base.minimum(zil::ZeroInflatedDistribution)
    min(minimum(zil.encdist), minimum(zil.posdist))
end
function Base.maximum(zil::ZeroInflatedDistribution)
    max(maximum(zil.encdist), maximum(zil.posdist))
end
function Distributions.insupport(zil::ZeroInflatedDistribution, x)
    minimum(zil) ≤ x ≤ maximum(zil)
end

# Statistics: mean, var, std, modes
function Statistics.mean(zil::ZeroInflatedDistribution)
    succprob(zil.encdist) * mean(zil.posdist)
end
function Statistics.var(zil::ZeroInflatedDistribution)
    ppos = succprob(zil.encdist)
    varpos = var(zil.posdist)
    meanpos = mean(zil.posdist)
    ppos * (varpos + meanpos^2) - (ppos * meanpos)^2
end
Statistics.std(zil::ZeroInflatedDistribution) = sqrt(var(zil))
Distributions.modes(zil::ZeroInflatedDistribution) = (zero(eltype(zil.encdist)), mode(zil.posdist))

# Facilities for random number generation
function Base.rand(rng::AbstractRNG, zil::ZeroInflatedDistribution)
    rand(rng, zil.encdist) * rand(rng, zil.posdist)
end
Base.rand(zil::ZeroInflatedDistribution) = rand(Random.GLOBAL_RNG, zil)


text:
#!/usr/bin/env bash
    #=
    exec julia --project="$(realpath $(dirname $0))/" "${BASH_SOURCE[0]}" "$@" -e "include(popfirst!(ARGS))" \
    "${BASH_SOURCE[0]}" "$@"
    =#

# we assume that `smart_choose_feats = true`
main_path = dirname(dirname(@__FILE__))
data_path = joinpath(main_path, "data")
main_image_path = joinpath(main_path, "data", "main")
alt_image_path = joinpath(main_path, "data", "alt")

num_classifiers = 10


text:
module NeuroIO
using Reexport

include("Struct.jl")
include("Common.jl")
include("NCS.jl")
include("PLX.jl")
include("NEV.jl")
include("Klusters.jl")

@reexport using .NEV
@reexport using .PLX
@reexport using .Klusters
@reexport using .NCS
using .Common
export validchannels, samplerate, bitspersample, voltagemultiplier, unitnumbers

end # module


text:
"""
encode(input)

Program to implement atbash cipher for the given sentence.A full description of the algorithm can be found on [wikipedia](https://en.wikipedia.org/wiki/Atbash)

# Arguments:
- `input` : The sentence needed to rotate

# Examples/Tests
```julia
julia> encode("test")
gvhg

julia> encode("abcdefghijklmnopqrstuvwxyz")
zyxwvutsrqponmlkjihgfedcba

julia> encode("hello")
svool

```

# Algorithm:

```julia

for r in input
    part *= xform(r)
    if length(part) >= 5
      push!(parts, part)
      part = ""
    end
  end
  if part != ""
    push!(parts, part)
  end
  return join(parts, " ")

```

# References:
https://en.wikipedia.org/wiki/Atbash

```

# Contributed by:- [Ihjass Thasbekha](https://github.com/Ihjass)
"""
function encode(input)
  input = replace(lowercase(input), reject_re => "")
  parts = []
  part = ""
  for r in input
    part *= xform(r)
    if length(part) >= 5
      push!(parts, part)
      part = ""
    end
  end
  if !isempty(part)

    push!(parts, part)
  end
  return join(parts, " ")
end
reject_re = r"[^a-z\d]+"
xform(r) = (r >= 'a' && r <= 'z') ? r = 25 - (r - 'a') + 'a' : r


text:
function param_fieldnames(::Type{Param}; include_subject_specific=false)
    a = [:G, :G0, :xi, :beta, :mu, :S, :L, :sigma, :bl_shape, :bl_scale, :delta]
    include_subject_specific ? vcat(a, [:b, :cind]) : a
end

function param_fieldnames(::Type{ParamSpline}; include_subject_specific=false)
    a = [:G, :G0, :xi, :beta, :mu, :S, :L, :sigma, :w, :delta]
    include_subject_specific ? vcat(a, [:b, :cind]) : a
end

function variables(param::T; include_subject_specific=false) where T <: AbstractParam
    names_str = Vector{String}(undef, 0)

    for par in param_fieldnames(typeof(param); include_subject_specific=include_subject_specific)
        prop = getproperty(param, par)
        prop_str = String(par)
        dims = size(prop)

        if length(dims) == 0
            names_str = vcat(names_str, "$(prop_str)")
        elseif length(dims) == 1
            names_str = vcat(names_str, ["$(prop_str)$(g)" for g in 1:dims[1]]...)
        elseif length(dims) == 2
            names_str = vcat(names_str, ["$(prop_str)$(g)$(k)" for k in 1:dims[2] for g in 1:dims[1]]...)
        else
            names_str = vcat(names_str, ["$(prop_str)$(g)$(k)$(l)" for l in 1:dims[3] for k in 1:dims[2] for g in 1:dims[1]]...)
        end
    end

    names_str ##Symbol.(names_str)
end

vectorize(a::T) where T <: Real = a
vectorize(a::T) where T <: AbstractArray = vec(a)

function vectorize(param::T; include_subject_specific=true) where T <: AbstractParam
    vcat([vectorize(getproperty(param, par)) for par in param_fieldnames(typeof(param); include_subject_specific=include_subject_specific)]...)
end

function summary(x::T, data::Data) where T <: Fit
    relabel!(x, MAP(x, data))
    summary(x)
end

summary(x::T) where T <: Fit = summary(x.chains)
summary(x::AbstractVector{T}) where T <: Chain = get_mcmc_chains(x)


text:
using GLMakie.Makie: Record
using GLMakie.GLFW
using GLMakie.ModernGL
using GLMakie.ShaderAbstractions
using GLMakie.ShaderAbstractions: Sampler
using GLMakie.StaticArrays
using GLMakie.GeometryBasics
using ReferenceTests.RNG

# A test case for wide lines and mitering at joints
@cell "Miter Joints for line rendering" begin
    scene = Scene()
    cam2d!(scene)
    r = 4
    sep = 4*r
    scatter!(scene, (sep+2*r)*[-1,-1,1,1], (sep+2*r)*[-1,1,-1,1])

    for i=-1:1
        for j=-1:1
            angle = pi/2 + pi/4*i
            x = r*[-cos(angle/2),0,-cos(angle/2)]
            y = r*[-sin(angle/2),0,sin(angle/2)]

            linewidth = 40 * 2.0^j
            lines!(scene, x .+ sep*i, y .+ sep*j, color=RGBAf(0,0,0,0.5), linewidth=linewidth)
            lines!(scene, x .+ sep*i, y .+ sep*j, color=:red)
        end
    end
    center!(scene)
    scene
end

@cell "Sampler type" begin
    # Directly access texture parameters:
    x = Sampler(fill(to_color(:yellow), 100, 100), minfilter=:nearest)
    scene = image(x)
    # indexing will go straight to the GPU, while only transfering the changes
    st = Stepper(scene)
    x[1:10, 1:50] .= to_color(:red)
    Makie.step!(st)
    x[1:10, end] .= to_color(:green)
    Makie.step!(st)
    x[end, end] = to_color(:blue)
    Makie.step!(st)
    st
end

# Test for resizing of TextureBuffer
@cell "Dynamically adjusting number of particles in a meshscatter" begin
    pos = Observable(RNG.rand(Point3f, 2))
    rot = Observable(RNG.rand(Vec3f, 2))
    color = Observable(RNG.rand(RGBf, 2))
    size = Observable(0.1*RNG.rand(2))

    makenew = Observable(1)
    on(makenew) do i
        pos[] = RNG.rand(Point3f, i)
        rot[] = RNG.rand(Vec3f, i)
        color[] = RNG.rand(RGBf, i)
        size[] = 0.1*RNG.rand(i)
    end

    fig, ax, p = meshscatter(pos,
        rotations=rot,
        color=color,
        markersize=size,
        limits=Rect3f(Point3(0), Point3(1))
    )
    Record(fig, [10, 5, 100, 60, 177]) do i
        makenew[] = i
    end
end

@cell "Explicit frame rendering" begin
    set_window_config!(renderloop=(screen) -> nothing)
    function update_loop(m, buff, screen)
        for i = 1:20
            GLFW.PollEvents()
            buff .= RNG.rand.(Point3f) .* 20f0
            m[1] = buff
            GLMakie.render_frame(screen)
            GLFW.SwapBuffers(GLMakie.to_native(screen))
            glFinish()
        end
    end
    fig, ax, meshplot = meshscatter(RNG.rand(Point3f, 10^4) .* 20f0)
    screen = Makie.backend_display(GLMakie.GLBackend(), fig.scene)
    buff = RNG.rand(Point3f, 10^4) .* 20f0;
    update_loop(meshplot, buff, screen)
    set_window_config!(renderloop=GLMakie.renderloop)
    fig
end

@cell "Contour and isosurface with correct depth" begin
    # Make sure shaders can recompile
    close(GLMakie.global_gl_screen())

    fig = Figure()
    left = LScene(fig[1, 1])
    contour!(left, [sin(i+j) * sin(j+k) * sin(i+k) for i in 1:10, j in 1:10, k in 1:10], enable_depth = true)
    mesh!(left, Sphere(Point3f(5), 6f0))
    right = LScene(fig[1, 2])
    volume!(right, [sin(2i) * sin(2j) * sin(2k) for i in 1:10, j in 1:10, k in 1:10], algorithm = :iso, enable_depth = true)
    mesh!(right, Sphere(Point3f(5), 6f0))
    fig
end


text:
import Base: unsigned, Unsigned, UInt8, UInt16, UInt32, UInt64, UInt128, show, typemax, typemin, <, <=, +, *, one, zero, convert, promote_rule, promote_typeof, iseven
@doc raw"""
    AInteger (or Avus Integer)

Integers that follow Avus definition for addition and multiplication. AIntegers also have
a distinct integer, ``\color{red}1``, where:
* ``0 < {\color{red}1} < 1``
* iseven(``{\color{red}1}``) == true

"""
abstract type AInteger <: Unsigned end
primitive type AInt8 <: AInteger 8 end
primitive type AInt16 <: AInteger 16 end
primitive type AInt32 <: AInteger 32 end
primitive type AInt64 <: AInteger 64 end
primitive type AInt128 <: AInteger 128 end

AInt8(x::Integer) = reinterpret(AInt8, convert(UInt8, x))
AInt16(x::Integer) = reinterpret(AInt16, convert(UInt16, x))
AInt32(x::Integer) = reinterpret(AInt32, convert(UInt32, x))
AInt64(x::Integer) = reinterpret(AInt64, convert(UInt64, x))
AInt128(x::Integer) = reinterpret(AInt128, convert(UInt128, x))

UInt8(x::AInt8) = reinterpret(UInt8, x)
UInt16(x::AInt8) = UInt16(reinterpret(UInt8, x))
UInt16(x::AInt16) = reinterpret(UInt16, x)
UInt32(x::AInt8) = UInt32(reinterpret(UInt8, x))
UInt32(x::AInt16) = UInt32(reinterpret(UInt16, x))
UInt32(x::AInt32) = reinterpret(UInt32, x)
UInt64(x::AInt8) = UInt64(reinterpret(UInt8, x))
UInt64(x::AInt16) = UInt64(reinterpret(UInt16, x))
UInt64(x::AInt32) = UInt64(reinterpret(UInt32, x))
UInt64(x::AInt64) = reinterpret(UInt64, x)
UInt128(x::AInt8) = UInt128(reinterpret(UInt8, x))
UInt128(x::AInt16) = UInt128(reinterpret(UInt16, x))
UInt128(x::AInt32) = UInt128(reinterpret(UInt32, x))
UInt128(x::AInt64) = UInt128(reinterpret(UInt64, x))
UInt128(x::AInt128) = reinterpret(UInt128, x)

function _print(io::IO, x::AInteger, y::Unsigned)
    if y == 0
        printstyled(io, "0"; color = :black)
    elseif x == red_one(x)
        printstyled(io, "1"; color = :red)
    else
        printstyled(io, Base.dec(y, 0, false); color = :black)
    end
end

Base.show(io::IO, x::AInt8) = _print(io, x, reinterpret(UInt8, x))
Base.show(io::IO, x::AInt16) = _print(io, x, reinterpret(UInt16, x))
Base.show(io::IO, x::AInt32) = _print(io, x, reinterpret(UInt32, x))
Base.show(io::IO, x::AInt64) = _print(io, x, reinterpret(UInt64, x))
Base.show(io::IO, x::AInt128) = _print(io, x, reinterpret(UInt128, x))

## integer promotions ##

# with different sizes, promote to larger type
promote_rule(::Type{AInt16}, ::Union{Type{Int8}, Type{UInt8}, Type{AInt8}}) = AInt16
promote_rule(::Type{AInt32}, ::Union{Type{Int16}, Type{Int8}, Type{UInt16}, Type{UInt8}, Type{AInt16}, Type{AInt8}}) = AInt32
promote_rule(::Type{AInt64}, ::Union{Type{Int16}, Type{Int32}, Type{Int8}, Type{UInt16}, Type{UInt32}, Type{UInt8}, Type{AInt16}, Type{AInt32}, Type{AInt8}}) = AInt64
promote_rule(::Type{AInt128}, ::Union{Type{Int16}, Type{Int32}, Type{Int64}, Type{Int8}, Type{UInt16}, Type{UInt32}, Type{UInt64}, Type{UInt8}, Type{AInt16}, Type{AInt32}, Type{AInt64}, Type{AInt8}}) = AInt128

promote_rule(::Type{Int16}, ::Type{AInt8}) = AInt16
promote_rule(::Type{Int32}, ::Union{Type{AInt16}, Type{AInt8}}) = AInt32
promote_rule(::Type{Int64}, ::Union{Type{AInt16}, Type{AInt32}, Type{AInt8}}) = AInt64
promote_rule(::Type{Int128}, ::Union{Type{AInt16}, Type{AInt32}, Type{AInt64}, Type{AInt8}}) = AInt128

# with mixed signed/unsigned/Avus and same size, Avus wins
promote_rule(::Type{AInt8}, ::Union{Type{Int8}, Type{UInt8}}) = AInt8
promote_rule(::Type{AInt16}, ::Union{Type{Int16}, Type{UInt16}}) = AInt16
promote_rule(::Type{AInt32}, ::Union{Type{Int32}, Type{UInt32}}) = AInt32
promote_rule(::Type{AInt64}, ::Union{Type{Int64}, Type{UInt64}}) = AInt64
promote_rule(::Type{AInt128}, ::Union{Type{Int128}, Type{UInt128}}) = AInt128

@doc raw"""
    red_one(x)
    red_one(T::type)

Return a ``\color{red}1`` for `x`: a value such that
``{\color{red}1} + 1 == {\color{red}1}``.
"""
function red_one(::AInteger) end

# now define min, max and red_one for each width of AInteger
typemin(::Type{AInt8}) = AInt8(0)
typemax(::Type{AInt8}) = AInt8(254)
red_one(::AInt8) = AInt8(255)
red_one(::Type{AInt8}) = AInt8(255)
red_one(x::Int8) = red_one(AInt8(x))
red_one(x::UInt8) = red_one(AInt8(x))
typemin(::Type{AInt16}) = AInt16(0)
typemax(::Type{AInt16}) = AInt16(65534)
red_one(::AInt16) = AInt16(65535)
red_one(::Type{AInt16}) = AInt16(65535)
red_one(x::Int16) = red_one(AInt16(x))
red_one(x::UInt16) = red_one(AInt16(x))
typemin(::Type{AInt32}) = AInt32(0)
typemax(::Type{AInt32}) = AInt32(4294967294)
red_one(::AInt32) = AInt32(4294967295)
red_one(::Type{AInt32}) = AInt32(4294967295)
red_one(x::Int32) = red_one(AInt32(x))
red_one(x::UInt32) = red_one(AInt32(x))
typemin(::Type{AInt64}) = AInt64(0)
typemax(::Type{AInt64}) = AInt64(0xfffffffffffffffe)
red_one(::AInt64) = AInt64(0xffffffffffffffff)
red_one(::Type{AInt64}) = AInt64(0xffffffffffffffff)
red_one(x::Int64) = red_one(AInt64(x))
red_one(x::UInt64) = red_one(AInt64(x))
typemin(::Type{AInt128}) = reinterpret(AInt128, UInt128(0))
typemax(::Type{AInt128}) = reinterpret(AInt128, Int128(-2))
red_one(::AInt128) = reinterpret(AInt128, Int128(-1))
red_one(::Type{AInt128}) = reinterpret(AInt128, Int128(-1))
red_one(x::Int128) = red_one(AInt128(x))
red_one(x::UInt128) = red_one(AInt128(x))

unsigned(x::AInt8) = x % UInt8
unsigned(x::AInt16) = x % UInt16
unsigned(x::AInt32) = x % UInt32
unsigned(x::AInt64) = x % UInt64
unsigned(x::AInt128) = x % UInt128

function string(x::AInteger)
    s = IOBuffer()
    _print(s, x, unsigned(x))
    String(take!(s))
end

@doc raw"""
    <(x::AInteger, y::AInteger)

for comparison purpose, 0 compares as ∞ and ``\color{red}1`` is the least value of any
AInteger
"""
function <(x::AInteger, y::AInteger)
    T = promote_typeof(x, y)
    xT, yT = x % T, y % T
    xu, yu = unsigned(xT), unsigned(yT)

    if xu == 0 || xu == yu
        return false
    elseif yu == 0 || xu == typemax(xu)
        return true
    else
        return xu < yu
    end
end

function <=(x::AInteger, y::AInteger)
    T = promote_typeof(x, y)
    xT, yT = x % T, y % T
    xu, yu = unsigned(xT), unsigned(yT)

    if xu == yu
        return true
    elseif yu == 0 || xu == typemax(xu)
        return true
    else
        return xu <= yu
    end
end

@doc raw"""
    iseven(x::AInteger) -> Bool

Return `true` if `x` is even (that is, divisible by 2) or ``\color{red}1``, and `false` otherwise.

# Examples
"""
iseven(n::AInteger) = n == red_one(n) || iseven(unsigned(n))

@doc raw"""
    +(x::Integer, y::Integer)

Avus Sum or min(x, y) where min(0) == ∞ and ``{\color{red}1} < 1``.
"""
function +(x::AInteger, y::AInteger)
    x < y ? x : y
end

MSB(n::AInteger) = MSB(unsigned(n))
MSB(n::Integer) = typeof(n)(sizeof(n)<<3 - leading_zeros(n) - (1))

"""
    *(x::Integer, y::Integer)

Avus Product or transitive relationship function.

Consider 3 vertices: u, v and w. Furthermore, assume that there is a path from u to v,
represented as x, and from v to w, represented as y.

The Avus Product provides the path from u to w, represented as z, e.g. z = x ⨰ y
"""
function *(x::AInteger, y::AInteger)
    # the implementation of the Avus product simply replaces the left most significant bit
    # of arg2 with the arg1

    T = promote_typeof(x, y)
    xT, yT = x % T, y % T
    xu, yu = unsigned(xT), unsigned(yT)

    # result is 0 by default (e.g. no transitive relationship)
    product = zero(T)

    # Handle zero case
    if xu == 0 || yu == 0
        return product
    end

    # get constants that conform to type needed
    max = typemax(xu)
    eins = one(xu)
    two = eins + eins

    # Handle red one and one case
    if xu == max
        if yu == eins
            product = xT
        else
            product = yT
        end
    elseif yu == max
        if xu == eins
            product = yT
        else
            product = xT
        end
    elseif xu == eins
        product = yT
    elseif yu == eins
        product = xT
    else
        s = sizeof(xu)
        if leading_zeros(xu) + leading_zeros(yu) <= s - 2
            # check to see if operation would overflow input type
            throw(OverflowError("product overflow(2): $x, $y, $(T)"))
        end

        msb = MSB(yu)

        product = T((yu & (two^msb - eins)) | (xu << msb))
        if product == max
            # this is a special case of overflow, when product == maxvalue for type
            throw(OverflowError("product overflow(1): $x, $y, $(T)"))
        end
    end

    product
end


text:
using Test, RecursiveArrayTools, StaticArrays

using DiffEqBase: NAN_CHECK

@test !NAN_CHECK(3.0+4.0im)
@test NAN_CHECK(NaN)

u1 = ones(3)
@test !NAN_CHECK(u1)
u1′ = copy(u1)
u1′[2] = NaN
@test NAN_CHECK(u1′)

u2 = [SA[1.0 1.0; 1.0 1.0] for i = 1:3]
@test !NAN_CHECK(u2)
u2′ = copy(u2)
u2′[2] = SA[1.0 NaN; 1.0 1.0]
@test NAN_CHECK(u2′)

u3 = VectorOfArray([ones(5), ones(5)])
@test !NAN_CHECK(u3)
u3′ = recursivecopy(u3)
u3′[2][3] = NaN
@test NAN_CHECK(u3′)

u4 = ArrayPartition(u1, u2, u3)
@test !NAN_CHECK(u4)
u4_1 = ArrayPartition(u1′, u2, u3)
@test NAN_CHECK(u4_1)
u4_2 = ArrayPartition(u1, u2′, u3)
@test NAN_CHECK(u4_2)
u4_3 = ArrayPartition(u1, u2, u3′)
@test NAN_CHECK(u4_3)

@test !NAN_CHECK(ArrayPartition(u4, u4))
@test NAN_CHECK(ArrayPartition(u4, u4_1))
@test NAN_CHECK(ArrayPartition(u4, u4_2))
@test NAN_CHECK(ArrayPartition(u4, u4_3))


text:
module GreekModule
   σ(x) = 1 ./ (1 + exp.(-x))
   logσ(x) = - log1p.(exp.(-x))
end


text:
using Documenter
using EasyHamiltonian

makedocs(sitename="EasyHamiltonian")


text:
export SparseClassSamples

function SparseClassSamples(label::AbstractArray{T, 4},active_channels,n_pix_per_class::Int) where T
  #get a few randomly selected pixels per class per image as the point-annotations
  cw = zeros(Float32,size(label)[1:2])
  for i in active_channels
    class_i_ind = findall(label[:,:,i,1].==1)
    class_i_ind_select = shuffle(class_i_ind)[1:n_pix_per_class]
    cw[class_i_ind_select] .= 1.0f0
  end
  return cw
end


text:
"""
    cmpr=compare(A,B)

Compare two simplices following lexicographic rules.
If B has lower rank, then A `cmpr ==-1`. If B has higher rank than A,
`cmpr==1`. If A and B are identical, `cmpr==0`.
"""
function compare(A,B)
    A=sort(A,rev=true)
    B=sort(B,rev=true)
    if length(B)!=length(A)
        println("Error: lists do not share same length!")
        println(A,B)
        return
    end

    cmpr=0
    for (a,b) in zip(A,B)#zip(A[end:-1:1],B[end:-1:1])
        if a!=b
            if a<b
                cmpr=1
                break
            else
                cmpr=-1
                break
            end
        end
    end
    return cmpr
end
function compare(a::Real,b::Real)
    cmpr=0
    if a!=b
        if a<b
            cmpr=1
        else
            cmpr=-1
        end
    end
    return cmpr
end

"""
    sort_smplx(list,smplx)

Helper function for sorting simplices in lexicographic manner.
Utilize divide an conquer strategy to find a simplex in ordered list of simpleces.
"""
function sort_smplx(list,smplx)
    idx=0
    list_length=length(list)
    if list_length>=2

        low=1
        high=list_length

        #check whether smplx is in list range
        # lower bound
        cmpr=compare(list[low],smplx)
        if cmpr<0
            idx=low
            return idx,true
        elseif cmpr==0
            idx=low
            return idx,false
        end

        #upper bound
        cmpr=compare(smplx,list[high])
        if cmpr<0
            idx=high+1
            return idx,true
        elseif cmpr==0
            idx=high
            return idx,false
        end

        if idx==0
            #divide et impera
            while high-low>1
                mid=Int(round((high+low)/2))

                cmpr=compare(list[mid],smplx)
                if cmpr>0
                    low=mid
                elseif cmpr<0
                    high=mid
                else
                    idx=mid
                    return idx, false
                end
            end
        end

        if idx==0
            #find location from last two possibilities
            cmpr_low=compare(list[low],smplx)
            cmpr_high=compare(smplx,list[high])
            if cmpr_low==0
                idx=low
                return idx,false
            elseif cmpr_high==0
                idx=high
                return idx,false
            else
                idx=high
                return idx,true
            end
        end

    elseif list_length==1
        cmpr=compare(list[1],smplx)
        if cmpr<0
            idx=1
            return idx, true
        elseif cmpr>0
            idx=2
            return idx,true
        else
            idx=1
            return idx,false
        end

    else list_length==0
        return 1,true
    end
    return nothing
end

"""
    insert_smplx!(list,smplx)

Insert simplex in ordered list of simplices.
"""
function insert_smplx!(list,smplx)
    #if !isa(smplx,Real)
    #    smplx=sort(smplx)
    #end
    idx,ins=sort_smplx(list,smplx)
    if ins
        insert!(list,idx,smplx)
    end
    return
end

"""
    idx = find_smplx(list,smplx)

Find index of simplex in ordered list of simplices.
If the simplex is not contained in the list, the returned index is `nothing`
"""
function find_smplx(list,smplx)
    #if !isa(smplx,Real)
    #    smplx=sort(smplx)
    #end
    idx,ins=sort_smplx(list,smplx)
    ins=!ins
    if ins
        return idx
    else
        return 0
    end
end