text:
# Licensed to the Apache Software Foundation (ASF) under one
# or more contributor license agreements.  See the NOTICE file
# distributed with this work for additional information
# regarding copyright ownership.  The ASF licenses this file
# to you under the Apache License, Version 2.0 (the
# "License"); you may not use this file except in compliance
# with the License.  You may obtain a copy of the License at
#
#   http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing,
# software distributed under the License is distributed on an
# "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied.  See the License for the
# specific language governing permissions and limitations
# under the License.

module TestMetric

using MXNet
using Base.Test

################################################################################
# Supporting functions
################################################################################

"""
Returns a random n x m array in which each column defines a discrete probability distribution.
Each column contains numbers between 0 and 1, and each column sums to 1.
"""
function generate_probs(n, m)
    # Init
    result = rand(n, m)

    # Normalize: ensure each column sums to 1
    for j = 1:m
        colsum = sum(result[:, j])
        for i = 1:n
            result[i, j] /= colsum
        end
    end
    result
end

function loglikelihood(labels::Vector{T}, probs::Array{T, 2}) where T <: AbstractFloat
    LL = 0.0
    eps = convert(T, 1.0e-8)
    for i = 1:size(labels, 1)
        LL += log(probs[Int(labels[i]) + 1, i] + eps)    # labels are zero-based
    end
    LL / size(labels, 1)
end

################################################################################
# Test Implementations
################################################################################

function test_ace()
    info("EvalMetric::ACE")
    n_categories   = 4
    n_observations = 100
    labels         = convert(Vector{Float32}, rand(0:(n_categories - 1), n_observations))    # MXNet uses Float32
    probs          = convert(Array{Float32}, generate_probs(n_categories, n_observations))
    LL             = loglikelihood(labels, probs)
    metric         = mx.ACE()    # For categorical variables, ACE == -LL
    mx._update_single_output(metric, labels, probs)
    LL_v2 = metric.ace_sum / metric.n_sample
    @test LL ≈ LL_v2 atol=1e-12
end

function test_nmse()
    info("EvalMetric::NMSE")

    @testset "EvalMetric::NMSE::update!" begin
        metric = mx.NMSE()
        labels = Array{mx.NDArray}(
            [mx.NDArray([100.0, 0.0]), mx.NDArray([10.0, 0.0])])
        preds = Array{mx.NDArray}(
            [mx.NDArray([20.0, 0.0]), mx.NDArray([2.0, 0.0])])

        mx.update!(metric, labels, preds)
        @test metric.nmse_sum ≈ 0.64 * 2
    end

    @testset "EvalMetric::NMSE::reset!" begin
        metric = mx.NMSE()
        metric.nmse_sum = sum(rand(10))
        metric.n_sample = 42

        mx.reset!(metric)

        @test metric.nmse_sum == 0.0
        @test metric.n_sample == 0
    end

    @testset "EvalMetric::NMSE::get" begin
        metric = mx.NMSE()
        metric.nmse_sum = 100.0
        metric.n_sample = 20

        @test mx.get(metric) == [(:NMSE, 5.0)]
    end
end

################################################################################
# Run tests
################################################################################
test_ace()
test_nmse()

end

max_stars_repo_name: Vikas-kum/incubator-mxnet

max_issues_repo_name: Vikas-kum/incubator-mxnet

max_forks_repo_name: Vikas-kum/incubator-mxnet


text:
module GitCommand

import Git_jll

export git

include("git.jl")
include("git_cmd.jl")

end # module

ext: jl


text:
module Test
@time eval(:(using Log4jl))

@time const logger = @Log4jl.logger "Log4jlExample"

function test()
error(logger, "test 1")
@error "test {}" 2
end

# println("Name:")
# @time @Log4jl.logger "AAA"
# println()

# println("Root:")
# @time @Log4jl.rootlogger
# println()

# println("Msg:")
# @time @Log4jl.logger MSG=LOG4JL_DEFAULT_STATUS_LEVEL
# println()

# println("Vaild URI:")
# @time @Log4jl.logger(URI="log4jl.yml") |> println

# println("Invalid URI:")
# @time @Log4jl.logger(URI="log4j1l.yml") |> println

# println("Bad parser:")
# @Log4jl.logger URI="log4jl.xxx"
# println()

# println("Code:")
# @Log4jl.logger begin 1 end
# println()

# println("ALL:")
# @Log4jl.logger URI="log4jl.yml" MSG=LOG4JL_DEFAULT_STATUS_LEVEL
# println()

# println("ALL:")
# @Log4jl.logger "AAA" MSG=LOG4JL_DEFAULT_STATUS_LEVEL begin 1 end
# println()
end

Test.test()


text:
using Distances
using Statistics
using SlidingDistancesBase
using SlidingDistancesBase: znorm

@test advance!(0) == 0
@test SlidingDistancesBase.setup_normalizer(Nothing, 1, 1) == (1, 1)
@test SlidingDistancesBase.normalize(Nothing, 1) == 1

@test SlidingDistancesBase.setup_normalizer(ZNormalizer, [1,0.5], [1,0.5])[1] == znorm([1,0.5])

n = 10
x = randn(100)
z = ZNormalizer(x,n)

@test length(z) == n
@test size(z) == (n,)
@test ndims(z) == 1
@test actuallastlength(z) == length(z.x)

advance!(z)  # Init
inds = 1:n
@test z == x[inds]
@test mean(z) ≈ mean(x[inds])
@test std(z) ≈ std(x[inds], corrected=false)
@test z.i == 1

@test @inferred(z[!,1]) == (x[1]-mean(z))/std(z)
@test @inferred(z[!,n]) == (x[n]-mean(z))/std(z)

advance!(z)
inds = inds .+ 1

@test z.i == 2
@test z == x[inds]
@test mean(z) ≈ mean(x[inds])
@test std(z) ≈ std(x[inds], corrected=false)

@test @inferred(z[!,1]) == (x[2]-mean(z))/std(z)
@test z[!,n] == (x[n+1]-mean(z))/std(z)

@test_throws BoundsError z[n+1]

@test normalize(ZNormalizer, z.x[2:11]) ≈ normalize(ZNormalizer, z)

for i = 1:89
    advance!(z)
end

@test z.bufi == 0
@test z[!, 1] ≈ (z[1]-mean(z))/std(z)
@test z.bufi == 1

@test normalize(ZNormalizer, z.x[91:end]) ≈ normalize(ZNormalizer, z) ≈ z.buffer
@test z.bufi == n

advance!(z)
@test_throws BoundsError z[n]

y = normalize(ZNormalizer, randn(10))
@test mean(y) ≈ 0 atol=1e-12
@test std(y, corrected=false) ≈ 1 atol=1e-12

## Diagonal ==================================================================================

n = 10
x = randn(2,100)
z = DiagonalZNormalizer(x,n)

@test length(z) == 2n
@test size(z) == (2,n)
@test ndims(z) == 2
@test actuallastlength(z) == size(z.x,2)

advance!(z)  # Init
inds = 1:n
@test Matrix(z) == x[:,inds]
@test mean(z) ≈ mean(x[:,inds], dims=2)
@test std(z) ≈ std(x[:,inds], corrected=false, dims=2)
@test z.i == 1
@test z[1,1] == x[1,1]
@test z[1:2,1] == x[1:2,1]

@test @inferred(z[!,1]) ≈ (x[:,1]-mean(z))./std(z)
@test @inferred(z[!,n]) ≈ (x[:,n]-mean(z))./std(z)

advance!(z)
inds = inds .+ 1

@test z[1,1] == x[1,2]
@test z[1:2,1] == x[1:2,2]

@test z.i == 2
@test Matrix(z) == x[:,inds]
@test mean(z) ≈ mean(x[:, inds], dims=2)
@test std(z) ≈ std(x[:, inds], dims=2, corrected=false)

@test @inferred(z[!,1]) ≈ (x[:,2]-mean(z))./std(z)
@test @inferred(z[!,n]) ≈ (x[:,n+1]-mean(z))./std(z)

@test z[1,1] == x[1,2]
@test z[1:2,1] == x[1:2,2]

@test normalize(DiagonalZNormalizer, z.x[:,2:11]) ≈ normalize(DiagonalZNormalizer, z)

for i = 1:89
    advance!(z)
end

@test z.bufi == 0
@test z[!, 1] ≈ (z[:,1]-mean(z))./std(z)
@test z.bufi == 1

@test normalize(DiagonalZNormalizer, z.x[:,91:end]) ≈ normalize(DiagonalZNormalizer, z) ≈ z.buffer
@test z.bufi == n

@test_throws BoundsError advance!(z)

y = normalize(DiagonalZNormalizer, randn(2,10))
@test mean(y) ≈ 0 atol=1e-12
@test std(y, corrected=false) ≈ 1 atol=1e-12

## Multidim ZNormalizer ==================================================================================

n = 10
x = 1 .+ randn(2,100)
z = ZNormalizer(x,n)

@test length(z) == 2n
@test size(z) == (2,n)
@test ndims(z) == 2
@test actuallastlength(z) == size(z.x,2)

advance!(z)  # Init
inds = 1:n
@test Matrix(z) == x[:,inds]
@test mean(z) ≈ mean(x[:,inds])
@test std(z) ≈ std(x[:,inds], corrected=false)
@test z.i == 1
@test z[1,1] == x[1,1]
@test z[1:2,1] == x[1:2,1]

@test @inferred(z[!,1]) ≈ (x[:,1].-mean(z))./std(z)
@test @inferred(z[!,n]) ≈ (x[:,n].-mean(z))./std(z)

advance!(z)
inds = inds .+ 1

@test z[1,1] == x[1,2]
@test z[1:2,1] == x[1:2,2]

@test z.i == 2
@test Matrix(z) == x[:,inds]
@test mean(z) ≈ mean(x[:, inds])
@test std(z) ≈ std(x[:, inds], corrected=false)

@test @inferred(z[!,1]) ≈ (x[:,2].-mean(z))./std(z)
@test @inferred(z[!,n]) ≈ (x[:,n+1].-mean(z))./std(z)

@test z[1,1] == x[1,2]
@test z[1:2,1] == x[1:2,2]

@test normalize(ZNormalizer, z.x[:,2:11]) ≈ normalize(ZNormalizer, z)

for i = 1:89
    advance!(z)
end

@test z.bufi == 0
@test z[!, 1] ≈ (z[:,1].-mean(z))./std(z)
@test z.bufi == 1

@test normalize(ZNormalizer, z.x[:,91:end]) ≈ normalize(ZNormalizer, z) ≈ z.buffer
@test z.bufi == n

y = normalize(ZNormalizer, randn(2,10))
@test mean(y) ≈ 0 atol=1e-12
@test std(y, corrected=false) ≈ 1 atol=1e-12

# NormNormalizer =============================================================================

@test SlidingDistancesBase.setup_normalizer(NormNormalizer, [1 0.5], [1 0.5])[1] ≈ [1 0.5]./norm([1 0.5])
@test SlidingDistancesBase.setup_normalizer(NormNormalizer, [1, 0.5], [1, 0.5])[1] ≈ [1, 0.5]./norm([1, 0.5])

n = 10
x = randn(2,100)
z = NormNormalizer(x,n)

@test length(z) == 2n
@test size(z) == (2,n)
@test ndims(z) == 2
@test actuallastlength(z) == size(z.x,2)

advance!(z)  # Init
inds = 1:n
@test Matrix(z) == x[:,inds]
@test z.i == 1
@test z[1,1] == x[1,1]
@test z[1:2,1] == x[1:2,1]

@test @inferred(z[!,1]) ≈ (x[:,1])./norm(z)
@test @inferred(z[!,n]) ≈ (x[:,n])./norm(z)

advance!(z)
inds = inds .+ 1

@test z[1,1] == x[1,2]
@test z[1:2,1] == x[1:2,2]

@test z.i == 2
@test Matrix(z) == x[:,inds]

@test @inferred(z[!,1]) ≈ (x[:,2])./norm(z)
@test @inferred(z[!,n]) ≈ (x[:,n+1])./norm(z)

@test z[1,1] == x[1,2]
@test z[1:2,1] == x[1:2,2]

@test normalize(NormNormalizer, z.x[:,2:11]) ≈ normalize(NormNormalizer, z)

for i = 1:89
    advance!(z)
end

@test z.bufi == 0
@test z[!, 1] ≈ (z[:,1])./norm(z)
@test z.bufi == 1

@test normalize(NormNormalizer, z.x[:,91:end]) ≈ normalize(NormNormalizer, z) ≈ z.buffer
@test z.bufi == n

@test_throws BoundsError advance!(z)

y = normalize(NormNormalizer, randn(2,10))
@test norm(y) ≈ 1 atol=1e-12

# SqNormNormalizer =============================================================================

@test SlidingDistancesBase.setup_normalizer(SqNormNormalizer, [1 0.5], [1 0.5])[1] ≈ [1 0.5]./norm([1 0.5])^2
@test SlidingDistancesBase.setup_normalizer(SqNormNormalizer, [1, 0.5], [1, 0.5])[1] ≈ [1, 0.5]./norm([1, 0.5])^2

n = 10
x = randn(2,100)
z = SqNormNormalizer(x,n)

@test length(z) == 2n
@test size(z) == (2,n)
@test ndims(z) == 2
@test actuallastlength(z) == size(z.x,2)

advance!(z)  # Init
inds = 1:n
@test Matrix(z) == x[:,inds]
@test z.i == 1
@test z[1,1] == x[1,1]
@test z[1:2,1] == x[1:2,1]

@test @inferred(z[!,1]) ≈ (x[:,1])./norm(z)^2
@test @inferred(z[!,n]) ≈ (x[:,n])./norm(z)^2

advance!(z)
inds = inds .+ 1

@test z[1,1] == x[1,2]
@test z[1:2,1] == x[1:2,2]

@test z.i == 2
@test Matrix(z) == x[:,inds]

@test @inferred(z[!,1]) ≈ (x[:,2])./norm(z)^2
@test @inferred(z[!,n]) ≈ (x[:,n+1])./norm(z)^2

@test z[1,1] == x[1,2]
@test z[1:2,1] == x[1:2,2]

@test normalize(SqNormNormalizer, z.x[:,2:11]) ≈ normalize(SqNormNormalizer, z)

for i = 1:89
    advance!(z)
end

@test z.bufi == 0
@test z[!, 1] ≈ (z[:,1])./norm(z)^2
@test z.bufi == 1

@test normalize(SqNormNormalizer, z.x[:,91:end]) ≈ normalize(SqNormNormalizer, z) ≈ z.buffer
@test z.bufi == n

@test_throws BoundsError advance!(z)

# y = normalize(SqNormNormalizer, randn(2,10))
# @test √norm(y) ≈ 1 atol=1e-12


text:
lib_path =  joinpath(dirname(@__FILE__()), "../bin/libScenGen.so")
const lib = normpath(lib_path)

function lcp_solve(M::Matrix{Float64}, q::Vector{Float64})
    n = length(q)
    w,z = Array{Float64}(undef, n), Array{Float64}(undef, n)
    return_code = ccall((:lcp_julia_solve, lib), Int, (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Int), w, z, M, q, n)
    if return_code == 1
        warn("lcp_solve reached maximum number of iterations")
    elseif return_code == 2
        warn("lcp_solve failed due to negative diagonal term")
    end
    return w,z
end


text:
using LogParser
using Base.Test
using GZip

#Read in gzipped file
jbapachecombined = readdlm(gzopen(joinpath(dirname(@__FILE__), "data", "juliabloggers-apachecombined.gz")), '\t')

#Parse file
jbparsed = parseapachecombined(vec(jbapachecombined))

#Test that array is 122,143 elements long
@test size(jbparsed)[1] == 122143

#Test that array is of type Array{ApacheLog,1}
@test typeof(jbparsed) == Array{ApacheLog,1}

#Test DataFrame method
jbparsed_df = DataFrame(jbparsed)

#Test that a DataFrame was returned
@test typeof(jbparsed_df) <: DataFrame

#Test that DataFrame is 122143x9
@test size(jbparsed_df) == (122143,9)


text:
function autofun(f,t::AbstractVector,h::AbstractVector{<:Vector{T}},maxlag;kwargs...) where T<:Number
    res = map(zip(t,h)) do (t,h)
        f(t,h,maxlag;kwargs...)
    end
    Tau = getindex.(res,1)
    Y = getindex.(res,2)
    τ = reduce(vcat, Tau)
    y = reduce(vcat, Y)

    perm = @views sortperm(τ)
    return τ[perm],y[perm]
end

StatsBase.autocov(t::AbstractVector,h::AbstractVector{<:Vector{T}},maxlag;kwargs...) where T<:Number = autofun(autocov,t,h,maxlag;kwargs...)
StatsBase.autocor(t::AbstractVector,h::AbstractVector{<:Vector{T}},maxlag;kwargs...) where T<:Number = autofun(autocor,t,h,maxlag;kwargs...)

"""
    τ, acf = autocov(t::AbstractVector, y::AbstractVector{T}, maxlag; normalize=false) where T <: Number

Calculate the autocovariance function for a signal `y` sampled at sample times `t`. This is useful if the signal is sampled non-equidistantly.

If `y` and `t` are vectors of vectors, the acf for all data is estimated and concatenated in the output, which is always two vectors of numbers.

The returned vectors are samples of the acf at lags `τ`. `τ` contains differences between entries in `t`, whereas `lags` supplied to the standard metohd of `autocov` are the differences between *indices* in `y` to compare. The result will in general be noisier that the stadard `autocov` and may require further estimation, e.g., by means of gaussian process regression or quantile regression. The following should hold `mean(acf[τ.==0]) ≈ var(y)`.

#Arguments:
- `t`: Vector of sample locations/times
- `y`: signal
- `maxlag`: Specifies the maximum difference between two values in `t` for which to calculate the ACF
- `normalize`: If false (default) the behaviour mimics the standard method from StatsBase and the ACF decreases automatically for higher lags, even if the theoretical ACF does not. If `true`, the theoretical ACF is estimated. `normalize = true` leads to sharper spectral estimates when taking the fft of the ACF. Note that the variance in the estimated ACF will be high for large lags if `true`.
"""
function StatsBase.autocov(t::AbstractVector{tT},y::AbstractVector{T},maxlag::Number; normalize=false) where {T<:Number, tT<:Number}
    isequidistant(t) || return _autocov(t,y,maxlag;normalize=normalize)
    ly = length(y)
    length(t) == ly || throw(ArgumentError("t and y must be the same length"))

    acf = zeros(T, (ly^2+ly)÷2)
    τ = zeros(tT, (ly^2+ly)÷2)
    k = 1
    @inbounds for i in eachindex(y)
        for j in 0:length(y)-1
            i+j > ly && continue
            τi = abs(t[i+j]-t[i])
            τi > maxlag && continue
            τ[k] = τi
            c = dot(y, 1:ly-j, y, 1+j:ly) / (ly-normalize*(j-1)) # the -1 is important, it corresponds to `corrected=true` for `var`. The variance increases for large lags and this becomes extra important close to lag=length(y)-1
            acf[k] = c
            k += 1
        end
    end

    perm = @views sortperm(τ[1:k-1])
    if all(x->x==y[1], y) || var(y) < eps()
        return τ[perm], zeros(T,length(perm))
    end
    return τ[perm],acf[perm]
end

"""
    τ, acf = autocor(t::AbstractVector, y::AbstractVector{T}, maxlag; normalize=false) where T <: Number

Calculate the auto correlation function for a signal `y` sampled at sample times `t`. This is useful if the signal is sampled non-equidistantly.

If `y` and `t` are vectors of vectors, the acf for all data is estimated and concatenated in the output, which is always two vectors of numbers.

The returned vectors are samples of the acf at lags `τ`. `τ` contains differences between entries in `t`, whereas `lags` supplied to the standard metohd of `autocov` are the differences between *indices* in `y` to compare. The result will in general be noisier that the stadard `autocov` and may require further estimation, e.g., by means of gaussian process regression or quantile regression. The following should hold `mean(acf[τ.==0]) ≈ var(y)`.

#Arguments:
- `t`: Vector of sample locations/times
- `y`: signal
- `maxlag`: Specifies the maximum difference between two values in `t` for which to calculate the ACF
- `normalize`: If false (default) the behaviour mimics the standard method from StatsBase and the ACF decreases automatically for higher lags, even if the theoretical ACF does not. If `true`, the theoretical ACF is estimated. `normalize = true` leads to sharper spectral estimates when taking the fft of the ACF. Note that the variance in the estimated ACF will be high for large lags if `true`.
"""
function StatsBase.autocor(t::AbstractVector{tT},y::AbstractVector{T},maxlag::Number; normalize=false) where {T<:Number, tT<:Number}
    isequidistant(t) || return _autocor(t,y,maxlag;normalize=normalize)
    ly = length(y)
    length(t) == ly || throw(ArgumentError("t and y must be the same length"))

    acf = zeros(T, (ly^2+ly)÷2)
    τ = zeros(tT, (ly^2+ly)÷2)
    k = 1
    dd = dot(y,y)
    @inbounds for i in eachindex(y)
        for j in 0:length(y)-1
            i+j > ly && continue
            τi = abs(t[i+j]-t[i])
            τi > maxlag && continue
            τ[k] = τi
            c = dot(y, 1:ly-j, y, 1+j:ly) / (dd*(ly - normalize*(j-0))/ly)
            acf[k] = c
            k += 1
        end
    end
    perm = @views sortperm(τ[1:k-1])
    if dd < eps()
        return τ[perm], ones(T,length(perm))
    end
    return τ[perm],acf[perm]
end

isequidistant(v::Union{<:UnitRange, <:StepRange, <:StepRangeLen}) = step(v) > 0
function isequidistant(v)
    d = v[2]-v[1]
    d > 0 || return false
    d = abs(d)
    for i in 3:length(v)
        abs(abs(v[i]-v[i-1])-d) < 20d*eps() || return false
    end
    true
end

function _autocov(t::AbstractVector{tT},y::AbstractVector{T},maxlag::Number; normalize=false) where {T<:Number, tT<:Number}
    ly = length(y)
    length(t) == ly || throw(ArgumentError("t and y must be the same length"))

    acf = zeros(T, (ly^2+ly)÷2)
    τ = zeros(tT, (ly^2+ly)÷2)
    k = 1
    @inbounds for i in eachindex(y)
        for j in 0:length(y)-1
            i+j > ly && continue
            τi = abs(t[i+j]-t[i])
            τi > maxlag && continue
            τ[k] = τi
            c = y[i]*y[i+j]
            acf[k] = c
            k += 1
        end
    end

    perm = @views sortperm(τ[1:k-1])
    if all(x->x==y[1], y) || var(y) < eps()
        return τ[perm], zeros(T,length(perm))
    end
    return τ[perm],acf[perm]
end

function _autocor(t::AbstractVector{tT},y::AbstractVector{T},maxlag::Number; normalize=false) where {T<:Number, tT<:Number}
    ly = length(y)
    length(t) == ly || throw(ArgumentError("t and y must be the same length"))

    acf = zeros(T, (ly^2+ly)÷2)
    τ = zeros(tT, (ly^2+ly)÷2)
    k = 1
    dd = var(y)
    @inbounds for i in eachindex(y)
        for j in 0:length(y)-1
            i+j > ly && continue
            τi = abs(t[i+j]-t[i])
            τi > maxlag && continue
            τ[k] = τi
            c = y[i]*y[i+j]/dd
            acf[k] = c
            k += 1
        end
    end
    perm = @views sortperm(τ[1:k-1])
    if dd < eps()
        return τ[perm], ones(T,length(perm))
    end
    zeroinds = τ .== 0
    acf[zeroinds] .= 1
    return τ[perm],acf[perm]
end


text:
isanonymous(x) = isanonymous(typeof(x))
isanonymous(::Type) = false

abstract type Cstruct_anonymous <: Cstruct end
isanonymous(::Type{<:Cstruct_anonymous}) = true
Base.show(io::IO, ::Type{CS}) where {CS<:Cstruct_anonymous} = print(io, "<anonymous-struct>")

abstract type Cunion_anonymous <: Cunion end
isanonymous(::Type{<:Cunion_anonymous}) = true
Base.show(io::IO, ::Type{CU}) where {CU<:Cunion_anonymous}  = print(io, "<anonymous-union>")

# <:Caggregate functions
function (::Type{CA})(cc::Cconst{CA}) where {CA<:Caggregate}
	T = concrete(CA)
	result = T(undef)
	setfield!(result, :mem, getfield(cc, :mem))
	return result
end

function (::Type{CA})(u::UndefInitializer) where {CA<:Caggregate}
	T = concrete(CA)
	return T(u)
end

function (::Type{CA})(; kwargs...) where {CA<:Caggregate}
	T = concrete(CA)
	result = T(undef)
	if isempty(kwargs)
		setfield!(result, :mem, map(zero, getfield(result, :mem)))
	else
		T <: Cunion && length(kwargs) > 1 && error("Expected only a single keyword argument when constructing Cunion's")
		foreach(kwarg -> _initproperty!(result, kwarg...), kwargs)
	end
	return result
end

Base.zero(::Type{CA}) where {CA<:Caggregate} = CA()
Base.convert(::Type{CA}, nt::NamedTuple) where {CA<:Caggregate} = CA(; nt...)
Base.isequal(x::CA, y::CA) where {CA<:Caggregate} = getfield(x, :mem) == getfield(y, :mem)
Base.:(==)(x::CA, y::CA) where {CA<:Caggregate} = isequal(x, y)

function Base.show(io::IO, ca::Union{Caggregate, Cconst{<:Caggregate}})
	if !(ca isa get(io, :typeinfo, Nothing))
		if ca isa Cconst
			print(io, typeof(ca).name, "(")
			show(io, nonconst(typeof(ca)))
		else
			show(io, typeof(ca))
		end
	end
	print(io, "(")
	for (ind, name) in enumerate(propertynames(typeof(ca)))
		print(io, ind > 1 ? ", " : "")
		print(io, name, "=")
		show(io, getproperty(ca, name))
	end
	print(io, ")")
	if !(ca isa get(io, :typeinfo, Nothing))
		ca isa Cconst && print(io, ")")
	end
end

macro cstruct(exprs...) return _caggregate(__module__, nothing, :cstruct, exprs...) end
macro cunion(exprs...) return _caggregate(__module__, nothing, :cunion, exprs...) end

function _caggregate(mod::Module, deps::Union{Vector{Pair{Symbol, Expr}}, Nothing}, kind::Symbol, name::Symbol)
	return _caggregate(mod, deps, kind, name, nothing, nothing)
end

function _caggregate(mod::Module, deps::Union{Vector{Pair{Symbol, Expr}}, Nothing}, kind::Symbol, body::Expr)
	return _caggregate(mod, deps, kind, nothing, body, nothing)
end

function _caggregate(mod::Module, deps::Union{Vector{Pair{Symbol, Expr}}, Nothing}, kind::Symbol, body::Expr, strategy::Symbol)
	return _caggregate(mod, deps, kind, nothing, body, strategy)
end

function _caggregate(mod::Module, deps::Union{Vector{Pair{Symbol, Expr}}, Nothing}, kind::Symbol, name::Union{Symbol, Nothing}, body::Union{Expr, Nothing})
	return _caggregate(mod, deps, kind, name, body, nothing)
end

todo"need to handle unknown-length aggregates with last field like `char c[]`"
function _caggregate(mod::Module, deps::Union{Vector{Pair{Symbol, Expr}}, Nothing}, kind::Symbol, name::Union{Symbol, Nothing}, body::Union{Expr, Nothing}, strategy::Union{Symbol, Nothing})
	isnothing(body) || Base.is_expr(body, :braces) || Base.is_expr(body, :bracescat) || error("Expected @$(kind) to have a `{ ... }` expression for the body of the type, but found `$(body)`")
	isnothing(body) && !isnothing(strategy) && error("Expected @$(kind) to have a body if alignment strategy is to be specified")
	isnothing(strategy) || (startswith(String(strategy), "__") && endswith(String(strategy), "__") && length(String(strategy)) > 4) || error("Expected @$(kind) to have packing specified as `__STRATEGY__`, such as `__packed__` or `__native__`")

	strategy = isnothing(strategy) ? :(CBinding.ALIGN_NATIVE) : :(Calignment{$(QuoteNode(Symbol(String(strategy)[3:end-2])))})
	isanon = isnothing(name)
	super = kind === :cunion ? (isanon ? :(Cunion_anonymous) : :(Cunion)) : (isanon ? :(Cstruct_anonymous) : :(Cstruct))
	name = isanon ? gensym("anonymous-$(kind)") : name
	escName = esc(name)
	concreteName = esc(Symbol(name, "\u200B"))  # this is so incredibly evil: appending a 0-width space to the type name

	isOuter = isnothing(deps)
	deps = isOuter ? Pair{Symbol, Expr}[] : deps
	if isnothing(body)
		push!(deps, name => quote
			abstract type $(escName) <: $(super) end
		end)
	else
		fields = []
		for arg in body.args
			arg = _expand(mod, deps, arg)
			if Base.is_expr(arg, :align, 1)
				align = arg.args[1]
				push!(fields, :(Calignment{$(align)}))
			elseif Base.is_expr(arg, :escape, 1) && !startswith(String(arg.args[1]), "##anonymous")
				# this is just a type definition, not a field
			else
				Base.is_expr(arg, :(::)) && length(arg.args) != 2 && error("Expected @$(kind) to have a `fieldName::FieldType` expression in the body of the type, but found `$(arg)`")

				argType = Base.is_expr(arg, :(::)) ? arg.args[end] : arg
				args = !Base.is_expr(arg, :(::)) ? nothing : arg.args[1]
				args = Base.is_expr(args, :tuple) ? args.args : (args,)
				for arg in args
					if isnothing(arg)
						push!(fields, :(Ctypespec($(argType))))
					elseif arg isa Symbol
						push!(fields, :(Pair{$(QuoteNode(arg)), Ctypespec($(argType))}))
					elseif Base.is_expr(arg, :escape, 1) && isnothing(arg.args[1])
						push!(fields, :(Ctypespec($(argType))))
					elseif Base.is_expr(arg, :escape, 1) && arg.args[1] isa Symbol
						push!(fields, :(Pair{$(QuoteNode(arg.args[1])), Ctypespec($(argType))}))
					elseif Base.is_expr(arg, :call, 3) && arg.args[1] === :(:) && arg.args[3] isa Integer
						push!(fields, :(Pair{$(QuoteNode(arg.args[2])), Ctypespec($(argType), Val($(arg.args[3])))}))
					elseif Base.is_expr(arg, :call, 3) && Base.is_expr(arg.args[1], :escape, 1) && arg.args[1].args[1] === :(:) && Base.is_expr(arg.args[3], :escape, 1) && arg.args[3].args[1] isa Integer
						push!(fields, :(Pair{$(QuoteNode(arg.args[2].args[1])), Ctypespec($(argType), Val($(arg.args[3].args[1])))}))
					else
						function _parseAugmentedField(a)
							a = deepcopy(a)
							if Base.is_expr(a, :(::), 2)
								(aname, atype) = _parseAugmentedField(a.args[2])
								isnothing(aname) || error("Unable to parse @$(kind) field, unexpected expression `$(a)`")
								return (Base.is_expr(a.args[1], :escape, 1) && a.args[1].args[1] isa Symbol ? a.args[1].args[1] : a.args[1], atype)
							elseif Base.is_expr(a, :curly) && length(a.args) >= 3 && a.args[1] in (:Carray, esc(:Carray), :(CBinding.Carray), Expr(:., esc(:CBinding), esc(QuoteNode(:Carray))))
								(aname, atype) = _parseAugmentedField(a.args[2])
								isnothing(aname) || error("Unable to parse @$(kind) field, unexpected expression `$(a)`")
								a.args[2] = atype
								if length(a.args) == 4 && Base.is_expr(a.args[4], :call, 2) && a.args[4].args[1] === :sizeof
									(aname, atype) = _parseAugmentedField(a.args[4].args[2])
									isnothing(aname) || error("Unable to parse @$(kind) field, unexpected expression `$(a)`")
									a.args[4].args[2] = atype
								end
								return (nothing, a)
							elseif Base.is_expr(a, :curly) && length(a.args) >= 1 && a.args[1] in (:Ptr, esc(:Ptr), :(Base.Ptr), Expr(:., esc(:Base), esc(QuoteNode(:Ptr))))
								if length(a.args) == 1
									push!(a.args, argType)
								else
									(aname, atype) = _parseAugmentedField(a.args[end])
									isnothing(aname) || error("Unable to parse @$(kind) field, unexpected expression `$(a)`")
									a.args[end] = atype
								end
								return (nothing, a)
							elseif Base.is_expr(a, :braces, 0)
								return (nothing, argType)
							else
								error("Expected @$(kind) to have a `fieldName`, `fieldName::Ptr{}`, or `fieldName::{}[N]` field name expression or some combination of them, but found `$(a)` within `$(arg)`")
							end
						end

						(aname, atype) = _parseAugmentedField(arg)
						push!(fields, :(Pair{$(QuoteNode(aname)), Ctypespec($(atype))}))
					end
				end
			end
		end

		push!(deps, name => quote
			abstract type $(escName) <: $(super) end
			mutable struct $(concreteName) <: $(escName)
				mem::NTuple{Ctypelayout(Ctypespec($(super), $(strategy), Tuple{$(fields...)})).size÷8, UInt8}

				$(concreteName)(::UndefInitializer) = new()
			end
			#=
				TypeSpec = Tuple{
					Pair{:sym1, Tuple{PrimType}},  # primitive field `sym1`
					Pair{:sym2, Tuple{PrimType, NBits}},  # bit field `sym2`
					Pair{:sym3, Ctypespec{FieldType, AggType, AggStrategy, AggTypeSpec}}},  # nested aggregate `sym3`
					Ctypespec{FieldType, AggType, AggStrategy, AggTypeSpec}},  # anonymous nested aggregate
					Calignment{align}  # alignment "field"
				}
			=#
			CBinding.concrete(::Type{$(escName)}) = $(concreteName)
			CBinding.concrete(::Type{$(concreteName)}) = $(concreteName)
			CBinding.strategy(::Type{$(concreteName)}) = $(strategy)
			CBinding.specification(::Type{$(concreteName)})  = Tuple{$(fields...)}
			Base.sizeof(::Type{$(escName)}) = sizeof(CBinding.concrete($(concreteName)))
		end)
	end

	return isOuter ? quote $(map(last, deps)...) ; $(escName) end : escName
end


text:
@enum PlatformID::UInt16 begin
    UNICODE   = 0
    MACINTOSH = 1
    ISO       = 2
    WINDOWS   = 3
    CUSTOM    = 4
end

@enum EncodingUnicode::UInt16 begin
    UNICODE_1_0                 = 0
    UNICODE_1_1                 = 1
    ISO_IEC_10646               = 2
    UNICODE_2_0_BMP             = 3
    UNICODE_2_0_FULL            = 4
    UNICODE_VARIATION_SEQUENCES = 5
    UNICODE_FULL                = 6
end

abstract type CmapSubtable end

struct SequentialMapGroup
    char_range::UnitRange{UInt32}
    start_glyph_id::UInt32
end

const ManyToOneRangeMappings = SequentialMapGroup

struct SegmentedCoverage <: CmapSubtable
    groups::Vector{SequentialMapGroup}
end

struct ByteEncodingTable <: CmapSubtable
    glyph_id_array::Vector{UInt8}
end

struct EncodingRecord
    platform_id::PlatformID
    encoding_id::UInt16
    subtable_offset::UInt32
end

struct CharToGlyph
    encoding_records::Vector{EncodingRecord}
    tables::Dict{Int,CmapSubtable}
end

function Base.parse(io::IO, table_mappings, ::Type{CharToGlyph})
    skip(io, 2)
    ntables = read(io, UInt16)
    records = map(1:ntables) do _
        EncodingRecord(
            PlatformID(read(io, UInt16)),
            [read(io, T) for T in fieldtypes(EncodingRecord)[2:3]]...
        )
    end
    d = Dict()
    foreach(records) do rec
        seek(io, table_mappings["cmap"].offset + rec.subtable_offset)
        format = read(io, UInt16)
        table = if format == 0
            skip(io, 4)
            ByteEncodingTable([read(io, UInt8) for _ in 1:256])
        elseif format == 12
            skip(io, 10)
            ngroups = read(io, UInt32)
            groups = map(1:ngroups) do i
                SequentialMapGroup(
                    read(io, UInt32):read(io, UInt32),
                    read(io, UInt32),
                )
            end
            SegmentedCoverage(groups)
        else
            return
        end
        d[format] = table
    end
    CharToGlyph(records, d)
end


text:
abstract type ALifeModel end

function init! end
function iterate! end
function set_state! end

function run_simulation(model::ALifeModel, num_iterations)
    organisms_lived = Organism[]

    init!(model)

    for _ in 1:num_iterations
        new_organisms = iterate!(model)
        append!(organisms_lived, new_organisms)
    end

    for (id, organism) in enumerate(organisms_lived)
        organism.id = id
    end

    return organisms_lived
end

max_stars_repo_name: tochsner/ALifeBenchmark

max_issues_repo_name: tochsner/ALifeBenchmark

max_forks_repo_name: tochsner/ALifeBenchmark


text:
struct SmoothBlock
	depth::Int
	i::Int
	lb::Float64
	ub::Float64
end
Base.isless(b1::SmoothBlock,b2::SmoothBlock) = isless(b1.ub-b1.lb, b2.ub-b2.lb)

adjustbounds(lb,ub,::WeightTree,::Any,::Any) = lb,ub
function adjustbounds(lb,ub,tree::MinMaxTree,depth,i)
	m = tree.mins[depth][i]
	M = tree.maxes[depth][i]
	min(m*lb,m*ub), max(M*lb,M*ub) # assumes that 0≤lb≤ub
end

function smoothexact(::Type{T},i,C,D,X,Y,xeval,leafSize) where T
	n = length(X)
	k1 = (i-1)*leafSize+1
	k2 = min(i*leafSize, n)

	# TODO: optimize (avoid allocations etc)
	# TODO: use dispatch instead
	if isnothing(Y)
		sum(ϕ.(C.*(X[k1:k2].-xeval),D))
	else
		sum(Y[k1:k2].*ϕ.(C.*(X[k1:k2].-xeval),D))
	end
end

function smoothbounds(::Type{T}, depth, i, C, D, X, tree, xeval) where T
	n = length(X)
	intervalSize = tree.leafSize * 2^(length(tree.intervalSums)-depth)

	k1 = (i-1)*intervalSize+1
	k2 = min(i*intervalSize, n)
	npoints = k2-k1+1

	a = C*(X[k1]-xeval)
	b = C*(X[k2]-xeval)

	xMean = C*(tree.intervalSums[depth][i]/npoints-xeval)

	if a<0 && b>0
		# special case, interval covers xeval
		# TODO: improve bounds? (We can use Jensen's inequality for this case too, if -1≤a≤b≤1, but ϕbounds() implementation needs to be updated)
		lb = ϕ(max(-a,b),D)
		ub = ϕ(zero(T))
		# return npoints.*(lb,ub)
		return adjustbounds((npoints.*(lb,ub))..., tree, depth, i)
	elseif a<0
		a,b,xMean = -b,-a,-xMean # mirror
	end

	# npoints.*ϕbounds(a,b,xMean)
	adjustbounds((npoints.*ϕbounds(a,b,xMean,D))..., tree, depth, i)
end

function smoothstep!(::Type{T},heap,C,D,X,Y,tree,xeval) where T
	block = pop!(heap)

	# remove existing bounds
	lb = -block.lb
	ub = -block.ub

	if block.depth >= depth(tree) # Fallback to exact sum
		# s = 0.0 # TODO: compute exact sum
		s = smoothexact(T,block.i,C,D,X,Y,xeval,tree.leafSize)
		lb += s
		ub += s
	else # create child blocks
		d = block.depth+1
		i = 2block.i-1 # first child
		nbrBlocks = length(tree.intervalSums[d])

		lb1,ub1 = smoothbounds(T,d,i,C,D,X,tree,xeval)
		lb += lb1
		ub += ub1
		lb1!=ub1 && push!(heap, SmoothBlock(d, i, lb1, ub1))

		if i<nbrBlocks
			lb2,ub2 = smoothbounds(T,d,i+1,C,D,X,tree,xeval)
			lb += lb2
			ub += ub2
			lb2!=ub2 && push!(heap, SmoothBlock(d, i+1, lb2, ub2))
		end
	end
	lb,ub
end

"""
	GaussianKernelSmoother(x,y; leafSize=10)

Create a `GaussianKernelSmoother` object of a set of data points with coordinates `x` and values `y`.
NB: References to `x` and `y` are stored in the GaussianKernelSmoother object, i.e. if you change `x` or `y` after creating the GaussianKernelSmoother, you will get incorrect results.

See also `smooth`, `density`.
"""
struct GaussianKernelSmoother{T,TX,TY}
	x::TX
	y::TY
	minMaxTree::MinMaxTree{T}
	weightTree::WeightTree{T}
end
function GaussianKernelSmoother(x, y; leafSize=10)
	@assert length(x)==length(y)
	if !issorted(x)
		perm = sortperm(x)
		x = x[perm]
		y = y[perm]
	end
	minMaxTree = MinMaxTree(x,y,leafSize)
	weightTree = WeightTree(minMaxTree)
	GaussianKernelSmoother(x,y,minMaxTree,weightTree)
end

# Treat as scalar for broadcasting, used by e.g. smooth
Base.Broadcast.broadcastable(gks::GaussianKernelSmoother) = Ref(gks)

function smoothapprox(gks::GaussianKernelSmoother{T},C,D,xeval;atol,rtol) where T
	heapn = BinaryMaxHeap{SmoothBlock}()
	heapd = BinaryMaxHeap{SmoothBlock}()

	lbn,ubn = smoothbounds(T,1,1,C,D,gks.x,gks.minMaxTree,xeval)
	push!(heapn,SmoothBlock(1,1,lbn,ubn))
	lbd,ubd = smoothbounds(T,1,1,C,D,gks.x,gks.weightTree,xeval)
	push!(heapd,SmoothBlock(1,1,lbd,ubd))

	while true
		lb = min(lbn/lbd,lbn/ubd)
		ub = max(ubn/lbd,ubn/ubd)

		isempty(heapn) && isempty(heapd) && return lb,ub
		isapprox(lb,ub;atol=atol,rtol=rtol) && return lb,ub

		if !isempty(heapn)
			lbn,ubn = (lbn,ubn) .+ smoothstep!(T,heapn,C,D,gks.x,gks.y,  gks.minMaxTree,xeval)
		end
		if !isempty(heapd)
			lbd,ubd = (lbd,ubd) .+ smoothstep!(T,heapd,C,D,gks.x,nothing,gks.weightTree,xeval)
		end
	end
end

function densityapprox(x,weightTree::WeightTree{T},C,D,xeval;atol,rtol) where T
	heap = BinaryMaxHeap{SmoothBlock}()
	lb,ub = smoothbounds(T,1,1,C,D,x,weightTree,xeval)
	push!(heap,SmoothBlock(1,1,lb,ub))
	while !isempty(heap) && !isapprox(lb,ub;atol=atol,rtol=rtol)
		lb,ub = (lb,ub) .+ smoothstep!(T,heap,C,D,x,nothing,weightTree,xeval)
	end
	lb,ub
end

# compute constant D such that the weight of the closest point is normalized to 1, which is critical for numerical precision.
function weightscale(x,C,xeval)
	r = searchsorted(x,xeval)
	D1 = (x[min(first(r),length(x))]-xeval)
	D2 = (x[max(last(r),1)]-xeval)
	C * (abs(D1)<abs(D2) ? D1 : D2)
end

"""
	smooth(gks::GaussianKernelSmoother, bandwidth, xeval; rtol=atol>0 ? 0 : 1e-3, atol=0)

Evaluate the Gaussian Kernel Smoother for the given `bandwidth` and `xeval`.
The value of the smoothed function `f` at `x₀` is given by `f(x₀) := ∑ᵢyᵢwᵢ / ∑ᵢwᵢ`, where `wᵢ := exp(-(x₀-xᵢ)²/2bandwidth²)`.

The accuracy of the result is controlled by `rtol` and `atol`. Lower and upper bounds `lb≤f(x₀)≤ub` are gradually improved until `isapprox(lb,ub;rtol=rtol,atol=atol)`.

# Examples
```julia-repl
julia> g = GaussianKernelSmoother([0.0,1.0],[2.0,4.0]);
julia> smooth.(g, 1, [0.1, 0.5, 0.9])
3-element Array{Float64,1}:
 2.802624679775096
 3.0
 3.197375320224904
```

See also `GaussianKernelSmoother`, `density`.
"""
function smooth(gks::GaussianKernelSmoother, bandwidth::Real, xeval::Real; atol=0, rtol=atol>0 ? 0 : 1e-3)
	x,y = gks.x,gks.y
	C = 1 / bandwidth
	D = weightscale(x,C,xeval)
	lb,ub = smoothapprox(gks,C,D,xeval; atol=atol, rtol=rtol)
	(lb+ub)/2
end

"""
	smooth(x, y, bandwidth, xeval; leafSize=10, rtol=atol>0 ? 0 : 1e-3, atol=0)

Evaluate the Gaussian Kernel Smoother of a set of data points with coordinates `x` and values `y`, with the specified `bandwidth`, at the coordinates in `xeval`.
The value of the smoothed function `f` at `x₀` is given by `f(x₀) := ∑ᵢyᵢwᵢ / ∑ᵢwᵢ`, where `wᵢ := exp(-(x₀-xᵢ)²/2bandwidth²)`.

The accuracy of the result is controlled by `rtol` and `atol`. Lower and upper bounds `lb≤f(x₀)≤ub` are gradually improved until `isapprox(lb,ub;rtol=rtol,atol=atol)`.

It is much more efficient to call `smooth` once with vector/matrix arguments for `xeval` and/or `bandwidth` than to call `smooth` multiple times.

# Examples
```julia-repl
julia> smooth([0.0,1.0], [2.0,4.0], 1.0, [0.1, 0.5, 0.9])
3-element Array{Float64,1}:
 2.802624679775096
 3.0
 3.197375320224904
```
See also `GaussianKernelSmoother`, `density`.
"""
function smooth(x::AbstractVector{T}, y::AbstractVector{T}, bandwidth, xeval; leafSize=10, kwargs...) where T
	gks = GaussianKernelSmoother(x,y; leafSize=leafSize)
	smooth.(gks, bandwidth, xeval; kwargs...)
end

function _density(x,weightTree::WeightTree,bandwidth,scale,xeval; atol=1e-20, rtol=atol>1e-20 ? 0 : 1e-3)
	C = 1 / bandwidth
	lb,ub = densityapprox(x,weightTree,C,0,xeval;atol=atol/scale,rtol)
	(lb+ub)/2*scale
end

"""
	density(x, bandwidth, xeval; leafSize=10, rtol=atol>1e-20 ? 0 : 1e-3, atol=1e-20)

Evaluate the density of a Gaussian Kernel Smoother of a set of data points with coordinates `x`, with the specified `bandwidth`, at the coordinates in `xeval`.
The value of the density `g` at `x₀` is given by `g(x₀) := C∑ᵢexp(-(x₀-xᵢ)²/2bandwidth²)`, where `C` is a constant ensuring that `g` is a probability density function.

The accuracy of the result is controlled by `rtol` and `atol`. Lower and upper bounds `lb≤f(x₀)≤ub` are gradually improved until `isapprox(lb,ub;rtol=rtol,atol=atol)`.

It is much more efficient to call `density` once with vector/matrix arguments for `xeval` and/or `bandwidth` than to call `density` multiple times.

# Examples
```julia-repl
julia> density([0.0,1.0], 1.0, [0.1, 0.5, 0.9])
3-element Array{Float64,1}:
 1.6619892900511566
 1.764993805169191
 1.6619892900511566
```
See also `smooth`, `GaussianKernelSmoother`.
"""
function density(x::AbstractVector, bandwidth, xeval; leafSize=10, kwargs...)
	x = issorted(x) ? x : sort(x)
	scale = 1 ./((sqrt(2π)*length(x)).*bandwidth)
	_density.(Ref(x),Ref(WeightTree(x,leafSize)),bandwidth,scale,xeval; kwargs...)
end

"""
	density(gks::GaussianKernelSmoother, bandwidth, xeval; leafSize=10, rtol=atol>1e-20 ? 0 : 1e-3, atol=1e-20)

Evaluate the density of a Gaussian Kernel Smoother, with the specified `bandwidth`, at the coordinates in `xeval`.
The value of the density `g` at `x₀` is given by `g(x₀) := C∑ᵢexp(-(x₀-xᵢ)²/2bandwidth²)`, where `C` is a constant ensuring that `g` is a probability density function.

The accuracy of the result is controlled by `rtol` and `atol`. Lower and upper bounds `lb≤f(x₀)≤ub` are gradually improved until `isapprox(lb,ub;rtol=rtol,atol=atol)`.

It is much more efficient to call `density` once with vector/matrix arguments for `xeval` and/or `bandwidth` than to call `density` multiple times.

# Examples
```julia-repl
julia> g = GaussianKernelSmoother([0.0,1.0],[2.0,4.0]);
julia> density.(g, 1.0, [0.1, 0.5, 0.9])
3-element Array{Float64,1}:
 1.6619892900511566
 1.764993805169191
 1.6619892900511566
```
See also `smooth`, `GaussianKernelSmoother`.
"""
function density(gks::GaussianKernelSmoother, bandwidth, xeval; kwargs...)
	scale = 1 ./((sqrt(2π)*length(gks.x)).*bandwidth)
	_density.(Ref(gks.x),Ref(gks.weightTree),bandwidth,scale,xeval; kwargs...)
end


text:
module Dalliance

export genomebrowser

const dalliancejs = Pkg.dir("Dalliance", "data", "dalliance-compiled.js")

import Base: writemime

using Compat, HttpServer, JSON, Bio, Bio.Intervals
import HttpCommon

const HTTP_SERVER_PORT = 8001

immutable DallianceBrowser
    genome_config::Dict
    sources::Vector{Dict}
end

function prepare_dataset(dataset::IntervalCollection{BEDMetadata})
    out = IOBuffer()
    write(out, BigBed, dataset)
    return takebuf_array(out)
end

function prepare_dataset(dataset::IntervalStream{BEDMetadata})
    ic = IntervalCollection{BEDMetadata}()
    for interval in dataset
        push!(ic, interval)
    end
    return prepare_dataset(ic)
end

function prepare_dataset{T <: Number}(dataset::IntervalStream{T})
    out = IOBuffer()
    write(out, BigWig, dataset)
    return takebuf_array(out)
end

# catchall
function prepare_dataset(dataset)
    error("Dalliance does not currently support data of type $(typeof(dataset))")
end

function genomebrowser(genome::String, datasets...)
    if !haskey(supported_genomes, genome)
        error("Genome $(genome) is not currently supported.")
    end

    sources = Dict[]
    for dataset in datasets
        id = string(Base.Random.uuid4())
        DATASETS[id] = prepare_dataset(dataset)
        push!(sources, @compat Dict(
            "name"   => "User Data", # TODO: some way to set names
            "bwgURI" => "http://localhost:$(HTTP_SERVER_PORT)/$(id)"))
    end

    push!(sources,
        @compat Dict("name"      => "Genome",
                     "twoBitURI" => supported_genomes[genome]["sequence"],
                     "tier_type" => "sequence"))

    return DallianceBrowser(supported_genomes[genome], sources)
end

const supported_genomes = @compat Dict(
    "GRCh38" => Dict(
        "speciesName" => "Human",
        "taxon"       => 9606,
        "auth"        => "GRCh",
        "version"     => "38",
        "ucscName"    => "hg38",
        "chr"         => "1",
        "viewStart"   => 1,
        "viewEnd"     => 100000,
        "sequence"    => "http://www.biodalliance.org/datasets/hg38.2bit"
    ),

    "GRCh37" => Dict(
        "speciesName" => "Human",
        "taxon"       => 9606,
        "auth"        => "GRCh",
        "version"     => "37",
        "ucscName"    => "hg19",
        "chr"         => "1",
        "viewStart"   => 1,
        "viewEnd"     => 100000,
        "sequence"    => "http://www.biodalliance.org/datasets/hg19.2bit"
    )

    # TODO: More genomes.
)

# Associate unique identifiers with datasets being served through http to dalliance.
const DATASETS = Dict{String, Vector{Uint8}}()

# HTTP Server instance.
let HTTP_SERVER = @compat Nullable{Server}()
    global start_server
    function start_server()
        if isnull(HTTP_SERVER)
            HTTP_SERVER = Nullable{Server}(Server(http_handler))
            @async run(get(HTTP_SERVER), host=IPv4(127,0,0,1), port=8001)
            #run(get(HTTP_SERVER), host=IPv4(127,0,0,1), port=8001)
            # TODO: if dalliance tries to access a track before the server is
            # ready, things break. Presumably there is a better solution than
            # this. Or does it...?
            sleep(5)
        end
    end
end

function http_handler(req::Request, res::Response)
    # debug nonsense:
    #logfile = open("/Users/dcjones/dialliance-log.txt", "a")
    #println(logfile, req)
    #close(logfile)

    h = HttpCommon.headers()
    h["Access-Control-Allow-Origin"] = "*"
    h["Access-Control-Allow-Headers"] = "Range"
    h["Access-Control-Max-Age"] = "36000"

    if req.method == "OPTIONS"
        return Response(200, h)
    elseif req.method != "GET" || !haskey(req.headers, "Range")
        return Response(501, h)
    end

    resource = match(r"^/([^?]*)", req.resource).captures[1]
    if !haskey(DATASETS, resource)
        return Response(404, h)
    end
    data = DATASETS[resource]

    m = match(r"bytes=(\d+)-(\d+)", req.headers["Range"])
    first = parse(Int, m.captures[1]) + 1
    last  = parse(Int, m.captures[2]) + 1
    range = first:min(last, length(data))

    return Response(data[range], h)
end

function writemime(io::IO, ::MIME"text/html", browser::DallianceBrowser)
    start_server()

    divid = string(Base.Random.uuid4())

    config = @compat Dict(
        "pageName" => divid,
        "chr"       => browser.genome_config["chr"],
        "viewStart" => browser.genome_config["viewStart"],
        "viewEnd"   => browser.genome_config["viewEnd"],
        "coordSystem" => Dict(
            "speciesName" => browser.genome_config["speciesName"],
            "taxon"       => browser.genome_config["taxon"],
            "auth"        => browser.genome_config["auth"],
            "version"     => browser.genome_config["version"],
            "ucscName"    => browser.genome_config["ucscName"]
        ),
        "sources" => browser.sources)

    #println(io, "<script language=\"javascript\" src=\"file://", dalliancejs, "\"></script>")
    println(io, "<script language=\"javascript\">\n", readall(dalliancejs), "</script>")
    println(io,
        """
            <script language="javascript">
                new window.Browser(
        """)

    println(io, json(config))
    println(io,
        """
                );
            </script>
        """)
    println(io, "<div id=\"", divid, "\"/>")
end

end


text:
@testcase "type destructing" begin
    @lift @data internal S{A, B} begin
        S_1{A, B} :: (a :: A, b :: B) => S{A, B}
    end

    s = S_1(1, "2")
    @test @match s begin
        ::S{A, B} where {A, B} => A == Int && B == String
    end

    @testset "curried destructing" begin
        s = S_1(S_1(1, 2), "2")
        @match s begin
            ::S{String} => false
            ::S{A} where A => A <: S{Int, Int}
            _ => false
        end
    end

    @test @match S_1(nothing, nothing) begin
        ::S{Int, Int} => false
        ::S{T} where T => T >: Nothing
        _ => false
    end
end


text:
using Test
using VSOP87

# Testing

@testset "VSOP87 Main Series" begin

    # Mercury
    pv,ierr = vsop87(2451545.0,0,1,0.0)
    @test pv ≈ [0.3870982122, 4.4026057779, 0.0446647518, 0.2007208958, 0.0406161541, 0.0456355031] atol=1E-10

    pv,ierr = vsop87(2122820.0,0,1,0.0)
    @test pv ≈ [0.3870983087, 5.5532501234, 0.0496101939, 0.1993701567, 0.0400177912, 0.0467766328] atol=1E-10

    # Venus
    pv,ierr = vsop87(2451545.0,0,2,0.0)
    @test pv ≈ [0.7233269304, 3.176135091, -0.0045086077, 0.0050312182, 0.0068248058, 0.0288221481] atol=1E-10

    pv,ierr = vsop87(2122820.0,0,2,0.0)
    @test pv ≈ [0.7233247251, 3.5192700749, -0.0047739162, 0.0053755162, 0.0055732704, 0.0291355398] atol=1E-10

    # Mars
    pv,ierr = vsop87(2451545.0,0,4,0.0)
    @test pv ≈ [1.5236789887, 6.2038757099, 0.0853133078, -0.0378067117, 0.0104705229, 0.0122862593] atol=1E-10

    pv,ierr = vsop87(2122820.0,0,4,0.0)
    @test pv ≈ [1.5236045525, 3.0145860938, 0.0818247408, -0.0434649474, 0.0102833001, 0.0132404725] atol=1E-10

    # Jupiter
    pv,ierr = vsop87(2451545.0,0,5,0.0)
    @test pv ≈ [5.2042662908, 0.5999772955, 0.0469877116, 0.0130817658, -0.0020729462, 0.0111943279] atol=1E-10

    pv,ierr = vsop87(2122820.0,0,5,0.0)
    @test pv ≈ [5.2026637693, 1.3989872239, 0.0462126899, 0.010229216, -0.0018026304, 0.0114181037] atol=1E-10

    # Saturn
    pv,ierr = vsop87(2451545.0,0,6,0.0)
    @test pv ≈ [9.5820161867, 0.872743095, 0.0003336009, 0.0557224686, -0.0086968779, 0.0198660071] atol=1E-10

    pv,ierr = vsop87(2122820.0,0,6,0.0)
    @test pv ≈ [9.53164267, 3.6837162325, 0.0034527067, 0.0581970245, -0.0093872672, 0.0193059636] atol=1E-10

    # Uranus
    pv,ierr = vsop87(2451545.0,0,7,0.0)
    @test pv ≈ [19.2294229491, 5.4713756706, -0.0438022683, 0.0073046995, 0.0018594791, 0.0064803904] atol=1E-10

    pv,ierr = vsop87(2122820.0,0,7,0.0)
    @test pv ≈ [19.2685452353, 1.0026901573, -0.0428389228, 0.0113976086, 0.0019824675, 0.0066094091] atol=1E-10

    # Neptune
    pv,ierr = vsop87(2451545.0,0,8,0.0)
    @test pv ≈ [30.1036169706, 5.3268987909, 0.0089053321, 0.0068181684, -0.0102818994, 0.0115021344] atol=1E-10

    pv,ierr = vsop87(2122820.0,0,8,0.0)
    @test pv ≈ [30.2002490334, 2.4115747749, 0.003226111, 0.0104510184, -0.010289848, 0.0114907692] atol=1E-10

    # Earth-Moon (EMB)
    pv,ierr = vsop87(2451545.0,0,9,0.0)
    @test pv ≈ [0.9999964222, 1.7534128816, -0.0037339066, 0.0162796346, -6.036e-7, 6.284e-7] atol=1E-10

    pv,ierr = vsop87(2122820.0,0,9,0.0)
    @test pv ≈ [1.0000134926, 1.8519621673, -0.0029638175, 0.0168402194, 0.0010301901, -5.30795e-5] atol=1E-10

end

@testset "VSOP87 A-Series   " begin

    # Mercury
    pv,ierr = vsop87(2451545.0,1,1,0.0)
    @test pv ≈ [-.1300934115 , -.4472876716 , -.0245983802 , .0213663982 , -.0064479797 , -.0024878668] atol=1E-10

    pv,ierr = vsop87(2122820.0,1,1,0.0)
    @test pv ≈ [.2146329139 , -.3752296250 , -.0503982597 , .0188319750 , .0155086039 , -.0005235332] atol=1E-10

    # Venus
    pv,ierr = vsop87(2451545.0,1,2,0.0)
    @test pv ≈ [-.7183022797, -.0326546017, .0410142975, .0007981261, -.0202952188, -.0003234551] atol=1E-10

    pv,ierr = vsop87(2122820.0,1,2,0.0)
    @test pv ≈ [-.6660158465, -.2753592311, .0357874175, .0075852982, -.0187913517, -.0006523229 ] atol=1E-10

    # Earth
    pv,ierr = vsop87(2451545.0,1,3,0.0)
    @test pv ≈ [-.1771354586, .9672416237, -.0000039000, -.0172076240, -.0031587881, .0000001069] atol=1E-10

    pv,ierr = vsop87(2122820.0,1,3,0.0)
    @test pv ≈ [-.2763146784, .9433985307, .0019115387, -.0168030101, -.0048929320, -.0000120286] atol=1E-10

    # Mars
    pv,ierr = vsop87(2451545.0,1,4,0.0)
    @test pv ≈ [1.3907159264, -.0134157043, -.0344677967, .0006714930, .0151872479, .0003016546] atol=1E-10

    pv,ierr = vsop87(2122820.0,1,4,0.0)
    @test pv ≈ [-1.6278485157, .3060194814, .0494215740, -.0019834234, -.0126075128, -.0002068580] atol=1E-10

    # Jupiter
    pv,ierr = vsop87(2451545.0,1,5,0.0)
    @test pv ≈ [4.0011740268, 2.9385810077, -.1017837501, -.0045683226, .0064432013, .0000755806] atol=1E-10

    pv,ierr = vsop87(2122820.0,1,5,0.0)
    @test pv ≈ [.4207861894, 5.1019591310, -.0280087663, -.0076034843, .0009707381, .0001701694] atol=1E-10

    # Saturn
    pv,ierr = vsop87(2451545.0,1,6,0.0)
    @test pv ≈ [6.4064068573, 6.5699929449, -.3690768029,  -.0042923542, .0038903162, .0001029504] atol=1E-10

    pv,ierr = vsop87(2122820.0,1,6,0.0)
    @test pv ≈ [-7.9395559173, -5.8435867016, .4165601597, .0029826848, -.0044772655, -.0000311309] atol=1E-10

    # Uranus
    pv,ierr = vsop87(2451545.0,1,7,0.0)
    @test pv ≈ [14.4318934159, -13.7343162527, -.2381421963, .0026781013, .0026726895, -.0000247716] atol=1E-10

    pv,ierr = vsop87(2122820.0,1,7,0.0)
    @test pv ≈ [11.8546461039, 15.5595370553, -.0950189343, -.0031646232, .0022092867, .0000505971] atol=1E-10

    # Neptune
    pv,ierr = vsop87(2451545.0,1,8,0.0)
    @test pv ≈ [16.8121116576, -24.9916630908, .1272190171, .0025792564, .0017769299, -.0000959082] atol=1E-10

    pv,ierr = vsop87(2122820.0,1,8,0.0)
    @test pv ≈ [-22.7985170870, 19.5994768858, .1206374034,  -.0020722808, -.0023628830, .0000962862] atol=1E-10

    # Earth-Moon (EMB)
    pv,ierr = vsop87(2451545.0,1,9,0.0)
    @test pv ≈ [-.1771591440, .9672192891, -.0000009536, -.0172031075, -.0031639188, .0000000258] atol=1E-10

    pv,ierr = vsop87(2122820.0,1,9,0.0)
    @test pv ≈ [-.2762837551, .9433889919, .0019144183, -.0168007538, -.0048864356, -.0000118518] atol=1E-10

end

@testset "VSOP87 E-Series   " begin

    # Mercury
    pv,ierr = vsop87(2451545.0,5,1,0.0)
    @test pv ≈ [-.1372349394, -.4500758422, -.0243922379, .0213717757, -.0064553869, -.0024879611] atol=1E-10

    pv,ierr = vsop87(2122820.0,5,1,0.0)
    @test pv ≈ [.2171621774, -.3801214149, -.0504926506, .0188386487, .0155090442, -.0005236921] atol=1E-10

    # Venus
    pv,ierr = vsop87(2451545.0,5,2,0.0)
    @test pv ≈ [-.7254438061, -.0354427729, .0412204390, .0008035035, -.0203026261, -.0003235494] atol=1E-10

    pv,ierr = vsop87(2122820.0,5,2,0.0)
    @test pv ≈ [-.6634865823, -.2802510190, .0356930259, .0075919721, -.0187909113, -.0006524817] atol=1E-10

    # Earth
    pv,ierr = vsop87(2451545.0,5,3,0.0)
    @test pv ≈ [-.1842769826, .9644534530, .0002022407, -.0172022466, -.0031661954, .0000000125] atol=1E-10

    pv,ierr = vsop87(2122820.0,5,3,0.0)
    @test pv ≈ [-.2737854147, .9385067424, .0018171487, -.0167963362, -.0048924916, -.0000121874] atol=1E-10

    # Mars
    pv,ierr = vsop87(2451545.0,5,4,0.0)
    @test pv ≈ [1.3835744053, -.0162038666, -.0342616574, .0006768704, .0151798404, .0003015603] atol=1E-10

    pv,ierr = vsop87(2122820.0,5,4,0.0)
    @test pv ≈ [-1.6253192611, .3011276886, .0493271851, -.0019767501, -.0126070725, -.0002070169] atol=1E-10

    # Jupiter
    pv,ierr = vsop87(2451545.0,5,5,0.0)
    @test pv ≈ [3.9940325025, 2.9357928287, -.1015776146, -.0045629454, .0064357940, .0000754864] atol=1E-10

    pv,ierr = vsop87(2122820.0,5,5,0.0)
    @test pv ≈ [.4233154540, 5.0970673238, -.0281031462, -.0075968103, .0009711785, .0001700105] atol=1E-10

    # Saturn
    pv,ierr = vsop87(2451545.0,5,6,0.0)
    @test pv ≈ [6.3992653459, 6.5672047374, -.3688706482, -.0042869773, .0038829091, .0001028562] atol=1E-10

    pv,ierr = vsop87(2122820.0,5,6,0.0)
    @test pv ≈ [-7.9370266647, -5.8484785039, .4164657927, .0029893590, -.0044768249, -.0000312897] atol=1E-10

    # Uranus
    pv,ierr = vsop87(2451545.0,5,7,0.0)
    @test pv ≈ [14.4247519568, -13.7371045087, -.2379360887, .0026834778, .0026652832, -.0000248664] atol=1E-10

    pv,ierr = vsop87(2122820.0,5,7,0.0)
    @test pv ≈ [11.8571754273, 15.5546452202, -.0951132908, -.0031579496, .0022097279, .0000504382] atol=1E-10

    # Neptune
    pv,ierr = vsop87(2451545.0,5,8,0.0)
    @test pv ≈ [16.8049701269, -24.9944513569, .1274251215, .0025846351, .0017695229, -.0000960030] atol=1E-10

    pv,ierr = vsop87(2122820.0,5,8,0.0)
    @test pv ≈ [-22.7959876638, 19.5945850298, .1205430330, -.0020656054, -.0023624426, .0000961274] atol=1E-10

    # Sun
    pv,ierr = vsop87(2451545.0,5,0,0.0)
    @test pv ≈ [-.0071415279, -.0027881715, .0002061418, .0000053774, -.0000074073, -.0000000943] atol=1E-10

    pv,ierr = vsop87(2122820.0,5,0,0.0)
    @test pv ≈ [.0025292639, -.0048917910, -.0000943913, .0000066739, .0000004404, -.0000001588] atol=1E-10

end


text:
@reexport baremodule LibKROME
  using CBinding: 𝐣𝐥

  𝐣𝐥.Base.include((𝐣𝐥.@__MODULE__), 𝐣𝐥.joinpath(𝐣𝐥.dirname(𝐣𝐥.@__DIR__), "deps", "libkrome.jl"))
end


text:
using POMDPs
#using RLESUtils, Observers
import Base: getindex, setindex!, haskey

# State node in the search tree
type StateNode
  n::Array{Int64,1} # number of visits at the node for each action
  Q::Array{Reward,1} # estimated value for each action
  action_map::Vector{Action} # for converting action idxs to action types
  nextstates::Vector{Nullable{UInt64}} #deterministic state transition means 1-1 mapping between actions and sp's
  explored::Bool #true if subtree is fully explored
end
function StateNode(nA, space::AbstractSpace, q0::Float64=0.0, explored::Bool=false)
  n = zeros(Int32, nA)
  Q = fill(q0, nA) #zeros(Reward, nA)
  # creates the action map
  am = Action[]
  for a in iterator(space)
    push!(am, a)
  end
  nextstates = [Nullable{UInt64}() for i = 1:nA]
  StateNode(n, Q, am, nextstates, explored)
end

# MCTS solver type
type MCTSSolver <: POMDPs.Solver
  n_iterations::Int64	# number of iterations during each action() call
  depth::Int64 # the max depth of the tree
  exploration_constant::Float64 # constant balancing exploration and exploitation
  maxmod::Bool
  rng::AbstractRNG # random number generator
  tree::Dict{UInt64, StateNode} # the search tree
end
# solver constructor
function MCTSSolver(;n_iterations::Int64=100,
                    depth::Int64=10,
                    exploration_constant::Float64=1.0,
                    maxmod::Bool=false, #use the max modification when updating q
                    rng = MersenneTwister(1))
  tree = Dict{UInt64, StateNode}() #key is hash(State)
  return MCTSSolver(n_iterations, depth, exploration_constant, maxmod, rng, tree)
end

# MCTS policy type
type MCTSPolicy <: POMDPs.Policy
  mcts::MCTSSolver # containts the solver parameters
  mdp::POMDP # model
  action_space::AbstractSpace # pre-allocated for rollout
  state::State # pre-allocated for sampling
  action::Action # pre-allocated for sampling
  distribution::AbstractDistribution # pre-allocated for memory efficiency
  rollout_policy::Policy # used in the rollout evaluation
  q0::Float64
  best_reward::Reward
  best_state::State
  best_at_eval::Int64
  totalevals::Int64
end
# policy constructor
function MCTSPolicy(mcts::MCTSSolver, mdp::POMDP,
                    rollout_policy=RandomPolicy(mdp, mcts.rng); # random policy is default
                    q0::Float64=0.0,
                    best_at_eval::Int64=0,
                    totalevals::Int64=0)
  # pre-allocate action space
  as = actions(mdp)
  # pre-allocate the state distrbution
  d = create_transition_distribution(mdp)
  s = create_state(mdp)
  a = create_action(mdp)
  best_reward = q0
  best_state = create_state(mdp)
  return MCTSPolicy(mcts, mdp, as, s, a, d, rollout_policy, q0, best_reward, best_state, best_at_eval, totalevals)
end

# for convenience - no computation is done in solve
function POMDPs.solve(solver::MCTSSolver, mdp::POMDP, policy::MCTSPolicy)
  policy
end

# retuns an approximately optimal action
function POMDPs.action(policy::MCTSPolicy, state::State)
  n_iterations = policy.mcts.n_iterations
  depth = policy.mcts.depth
  tree = policy.mcts.tree
  # build the tree
  for n = 1:n_iterations
    isexplored(tree, state) && break #can't do better, exit early

    #@notify_observer(policy.observer, "iteration", [n])
    simulate(policy, state, depth)
  end
  snode = tree[state]
  # find the index of action with highest q val
  i = indmax(snode.Q)
  #@notify_observer(policy.observer, "Qs", [snode.Q])
  #@notify_observer(policy.observer, "best_i", [i])
  #@notify_observer(policy.observer, "tree", [tree, state])
  # use map to conver index to mdp action
  return snode.action_map[i]
end

# runs a simulation from the passed in state to the specified depth
function POMDPs.simulate(policy::MCTSPolicy, state::State, depth::Int64)
  # model parameters
  mdp = policy.mdp
  na = n_actions(mdp, state)
  discount_factor = discount(mdp)
  sp = policy.state
  rng = policy.mcts.rng

  # solver parameters
  n_iterations = policy.mcts.n_iterations
  tree = policy.mcts.tree
  exploration_constant = policy.mcts.exploration_constant

  # if unexplored state add to the tree
  isnew = false
  if !haskey(tree, state)
    as = actions(mdp, state, policy.action_space)
    tree[state] = StateNode(na, as, policy.q0)
    isnew = true
  end

  snode = tree[state]

  if isterminal(mdp, state) #terminal states don't have corresponding snodes?
    snode.explored = true
    r = reward(mdp, state)
    policy.totalevals += 1 #expr eval call
    #@notify_observer(policy.observer, "terminal_reward", [r, state])
    if r > policy.best_reward
      policy.best_reward = r
      copy!(policy.best_state, state)
      policy.best_at_eval = policy.totalevals
    end
    return r
  end

  if isnew #new but not terminal
    return rollout(policy, depth, state)
  end

  # don't choose fully-explored subtrees
  best_i = -1
  best_val = -realmax(Float64)
  for i = 1:na
    if !isexplored(tree, snode.nextstates[i])
      #selection criteria
      val =  snode.Q[i] + exploration_constant * real(sqrt(complex(log(sum(snode.n)) / snode.n[i])))

      if best_i < 0 || val > best_val #accept the first unexplored one or accept if multiple ties for best
        best_val = val
        best_i = i
      end
    end
  end
  i = best_i
  a = snode.action_map[i]
  r = reward(mdp, state, a)
  # transition to a new state
  d = policy.distribution
  d = transition(mdp, state, a, d)
  sp = rand!(rng, sp, d)

  #add this transition to snode
  if !isnull(snode.nextstates[i])
    @assert get(snode.nextstates[i]) == hash(sp)
  else
    snode.nextstates[i] = Nullable{UInt64}(hash(sp))
  end

  # update the Q and n values
  q = r + discount_factor * simulate(policy, sp, depth - 1)
  snode.n[i] += 1 # increase number of node visits by one

  if policy.mcts.maxmod #maxmod
    snode.Q[i] = max(snode.Q[i], q)
  else
    snode.Q[i] += ((q - snode.Q[i]) / (snode.n[i])) # moving average of Q value
  end

  #if all children are fully-explored, mark this node as fully-explored
  snode.explored = isexplored(tree, snode.nextstates)

  return q
end

# recursive rollout to specified depth, returns the accumulated discounted reward
function rollout(policy::MCTSPolicy, depth::Int64, state::State)
  mdp = policy.mdp
  # finish when depth is zero or reach terminal state
  if isterminal(mdp, state)
    r = reward(mdp, state)
    policy.totalevals += 1 #expr eval call
    if r > policy.best_reward
      policy.best_reward = r
      copy!(policy.best_state, state)
      policy.best_at_eval = policy.totalevals
    end
    return r
  end
  d = policy.distribution
  discount_factor = discount(mdp)
  rng = policy.mcts.rng
  sp = policy.state
  # follow the rollout policy
  a = action(policy.rollout_policy, state)
  # compute reward
  r = reward(mdp, state, a)
  # sample the next state
  d = transition(mdp, state, a, d)
  sp = rand!(rng, sp, d)
  return r + discount_factor * rollout(policy, depth - 1, sp)
end

#override to use the hash instead
Base.getindex{V}(D::Dict{UInt64,V}, key::State) = D[hash(key)]
Base.getindex{V}(D::Dict{UInt64,V}, key::Nullable{UInt64}) = D[get(key)]
Base.setindex!{V}(D::Dict{UInt64,V}, value::V, key::State) = D[hash(key)] = value
Base.setindex!{V}(D::Dict{UInt64,V}, value::V, key::Nullable{UInt64}) = D[get(key)] = value
Base.haskey{V}(D::Dict{UInt64,V}, key::State) = haskey(D, hash(key))
Base.haskey{V}(D::Dict{UInt64,V}, key::Nullable{UInt64}) = haskey(D, get(key))

isexplored(tree::Dict{UInt64,StateNode}, s::State) = haskey(tree, s) && tree[s].explored
isexplored(tree::Dict{UInt64,StateNode}, h::Nullable{UInt64}) = !isnull(h) && tree[get(h)].explored

function isexplored(tree::Dict{UInt64,StateNode}, nextstates::Vector{Nullable{UInt64}})
  explored = true
  for h in nextstates
    if !isexplored(tree, h)
      explored = false
      break
    end
  end
  return explored
end

max_stars_repo_name: sisl/ExprSearch

max_issues_repo_name: sisl/ExprSearch


text:
module Resource
using HTTP, JSON3, Sockets
using ..Model, ..Service

const ROUTER = HTTP.Router()

health(req) = Dict("status"=>"ok")
HTTP.@register(ROUTER, "GET", "/", health)

createPaper(req) = Service.createPaper(JSON3.read(req.body))::Paper
HTTP.@register(ROUTER, "POST", "/paper", createPaper)

getPaper(req) = Service.getPaper(
                                 parse(Int, HTTP.URIs.splitpath(req.target)[2])
                                )::Paper
HTTP.@register(ROUTER, "GET", "/paper/*", getPaper)

updatePaper(req) = Service.updatePaper(
                                       parse(Int, HTTP.URIs.splitpath(req.target)[2]),
                                       JSON3.read(req.body, Paper)
                                      )::Paper
HTTP.@register(ROUTER, "PUT", "/paper/*", updatePaper)

deletePaper(req) = Service.deletePaper(
                                  parse(Int, HTTP.URIs.splitpath(req.target)[2])
                                 )
HTTP.@register(ROUTER, "DELETE", "/paper/*", deletePaper)

pickPaperToRead(req) = Service.pickPaperToRead()::Paper
HTTP.@register(ROUTER, "GET", "/paper", pickPaperToRead)

function requetHandler(req)
    obj = HTTP.handle(ROUTER, req)
    resp = HTTP.Response(200, JSON3.write(obj))
    return resp
end

function run()
    HTTP.serve(requetHandler, "0.0.0.0", 8880)
end

end


text:
##### Beginning of file

function _is_filesystem_root(path::AbstractString)::Bool
    path::String = abspath(strip(path))
    if path == dirname(path)
        return true
    else
        return false
    end
end

function _is_package_directory(path::AbstractString)::Bool
    path::String = abspath(strip(path))
    if isfile(joinpath(path, "Project.toml"))
        return true
    else
        return false
    end
end

function _find_package_directory(path::AbstractString)::String
    path::String = abspath(strip(path))
    if _is_package_directory(path)
        return path
    elseif _is_filesystem_root(path)
        error(string("Could not find the Project.toml file"))
    else
        result = _find_package_directory(dirname(path))
        return result
    end
end

"""
    package_directory()::String

Return the Snapshots package directory.
"""
function package_directory()::String
    result::String = _find_package_directory(abspath(strip(@__FILE__)))
    return result
end

function _location(m::Method)::String
    result::String = abspath(first(functionloc(m)))
    return result
end

function _location(f::Function)::String
    result::String = abspath(first(functionloc(f)))
    return result
end

function _location(f::Function, types::Tuple)::String
    result::String = abspath(first(functionloc(f, types)))
    return result
end

function _location(m::Module)::String
    result::String = abspath(_location(getfield(m, :eval)))
    return result
end

"""
    package_directory(parts...)::String

Equivalent to `abspath(joinpath(abspath(package_directory()), parts...))`.
"""
function package_directory(parts...)::String
    result::String = abspath(joinpath(abspath(package_directory()), parts...))
    return result
end

"""
    package_directory(m::Method)::String

If method `m`
is part of a Julia package, returns the package root directory.

If method `m`
is not part of a Julia package, throws an error.
"""
function package_directory(m::Method)::String
    m_module_directory::String = abspath(_location(m))
    m_package_directory::String = abspath(
        _find_package_directory(m_module_directory)
        )
    return m_package_directory
end

"""
    package_directory(m::Method, parts...)::String

Equivalent to
`result = abspath(joinpath(abspath(package_directory(m)), parts...))`.
"""
function package_directory(m::Method, parts...)::String
    result::String = abspath(joinpath(abspath(package_directory(m)), parts...))
    return result
end

"""
    package_directory(f::Function)::String

If function `f`
is part of a Julia package, returns the package root directory.

If function `f`
is not part of a Julia package, throws an error.
"""
function package_directory(f::Function)::String
    m_module_directory::String = abspath(_location(f))
    m_package_directory::String = abspath(
        _find_package_directory(m_module_directory)
        )
    return m_package_directory
end

"""
    package_directory(f::Function, parts...)::String

Equivalent to
`result = abspath(joinpath(abspath(package_directory(f)), parts...))`.
"""
function package_directory(f::Function, parts...)::String
    result::String = abspath(joinpath(abspath(package_directory(f)), parts...))
    return result
end

"""
    package_directory(f::Function, types::Tuple)::String

If function `f` with type signature `types`
is part of a Julia package, returns the package root directory.

If function `f` with type signature `types`
is not part of a Julia package, throws an error.
"""
function package_directory(f::Function, types::Tuple)::String
    m_module_directory::String = abspath(_location(f, types))
    m_package_directory::String = abspath(
        _find_package_directory(m_module_directory)
        )
    return m_package_directory
end

"""
    package_directory(f::Function, types::Tuple, parts...)::String

Equivalent to
`result = abspath(joinpath(abspath(package_directory(f, types)), parts...))`.
"""
function package_directory(f::Function, types::Tuple, parts...)::String
    result::String = abspath(joinpath(abspath(package_directory(f, types)), parts...))
    return result
end

"""
    package_directory(m::Module)::String

If module `m`
is part of a Julia package, returns the package root directory.

If module `m`
is not part of a Julia package, throws an error.
"""
function package_directory(m::Module)::String
    m_module_directory::String = abspath(_location(m))
    m_package_directory::String = abspath(
        _find_package_directory(m_module_directory)
        )
    return m_package_directory
end

"""
    package_directory(m::Module, parts...)::String

Equivalent to
`result = abspath(joinpath(abspath(package_directory(m)), parts...))`.
"""
function package_directory(m::Module, parts...)::String
    result::String = abspath(joinpath(abspath(package_directory(m)), parts...))
    return result
end

##### End of file


text:
immutable InputLayerState <: AbstractScatteredLayerState
    blobs::Vector{Mocha.Blob}
    layer::InputLayer
end

function InputLayerState(backend::Mocha.Backend, layer::InputLayer)
    ranges = ntuple(k -> kthrange(layer, k), ndims(layer.data))
    blob = ScatteredBlob(
        DataStructures.SortedDict((Path() => Node(layer.data, ranges),)))
    blobs = Mocha.Blob[blob]
    return InputLayerState(blobs, layer)
end

function Mocha.setup(
        backend::Mocha.Backend,
        layer::InputLayer,
        inputs::Vector{Mocha.Blob},
        diffs::Vector{Mocha.Blob})
    @assert length(diffs) == 0
    @assert length(inputs) == 0
    return InputLayerState(backend, layer)
end

kthrange(layer::InputLayer, k::Int) =
    (PathKey(layer.symbols[k]) => (1:1:size(layer.data, k)))


text:
######################################################################
# third pass: order collated dataset by start dates
######################################################################

using AMDB: data_path
using JLD2
using DocStringExtensions: SIGNATURES
using FileIO: load, save
using LargeColumns: meta_path, MmappedColumns
using ProgressMeter: @showprogress
import Base.Mmap: sync!

collated = MmappedColumns(data_path("collated"));

meta = load(meta_path(collated, "meta.jld2"))

const date_col_index = findfirst(meta[AMDB.META_COLUMN_NAMES], :STICHTAG)

@assert date_col_index > 0

"""
    $SIGNATURES

Iterating through indices in `ix`, sort spells for each individual by the start
date in `date_col_index`.
"""
function order_date_start!(ix, collated)
    @showprogress for r in ix
        collated[r] .= sort(collated[r], by = x -> x[date_col_index])
    end
end

order_date_start!(meta[AMDB.META_IX], collated) # sort by dates
sync!(collated)                 # make sure it is synced to disk


text:
# core plotting functions

plotbeta(α, β)
plotbeta() = plotBeta(1, 1)

plotunif(a, b)
plotunif() = plotUnif(0.0, 1.0)

plotbinom(n, p)
plotbinom() = plotBinom(1, 0.5)

plotbernoulli(p)
plotbernoulli() = plotBernolli(0.5)

plottest(t)


text:
using pigpiod_if2_jll
include("constants.jl")
# I2C
#  i2c_open                     Opens an I2C device
#  int i2c_open(int pi, unsigned i2c_bus, unsigned i2c_addr, unsigned i2c_flags)
#  Returns handle for the device at address addr on bus "bus"
function i2c_open(p::Int32, bus::Int, addr::Int, flags::Int)
    return ccall((:i2c_open, libpigpiod_if2), Int32, (Int32, UInt32, UInt32, UInt32), p, bus, addr, flags)
end

#  i2c_close                    Closes an I2C device
#  int i2c_close(int pi, unsigned handle)
#  Closes i2c device with associated handle
function i2c_close(p::Int32, handle::Int32)
    return ccall((:i2c_close, libpigpiod_if2), Int32, (Int32, UInt32), p, handle)
end

#  i2c_write_quick              smbus write quick
#  int i2c_write_quick(int pi, unsigned handle, unsigned bit)
#  Sends a single bit (0 ot 1) to the device associated with "handle"
function i2c_write_quick(p::Int32, handle::Int32, bit::Int)
    return ccall((:i2c_write_quick, libpigpiod_if2), Int32, (Int32, UInt32, UInt32), p, handle, bit)
end

#  i2c_read_byte                smbus read byte
#  int i2c_read_byte(int pi, unsigned handle)
#  Reads single byte from device with handle
function i2c_read_byte(p::Int32, handle::Int32)
    return ccall((:i2c_read_byte, libpigpiod_if2), Int32, (Int32, UInt32), p, handle)
end

#  i2c_write_byte               smbus write byte
#  int i2c_write_byte(int pi, unsigned handle, unsigned bVal)
#  Sends a single byte to the device with handle
function i2c_write_byte(p::Int32, handle::Int32, bVal::Int)
    return ccall((:i2c_write_byte, libpigpiod_if2), Int32, (Int32, UInt32, UInt32), p, handle, bVal)
end

#  i2c_read_byte_data           smbus read byte data
#  int i2c_read_byte_data(int pi, unsigned handle, unsigned i2c_reg)
#  Reads single byte from specified register of device with associated handle
function i2c_read_byte_data(p::Int32, handle::Int32, reg::Int)
    return ccall((:i2c_read_byte_data, libpigpiod_if2), Int32, (Int32, UInt32, UInt32), p, handle, reg)
end

#  i2c_write_byte_data          smbus write byte data
#  int i2c_write_byte_data(int pi, unsigned handle, unsigned i2c_reg, unsigned bVal)
#  Writes single byte to the specified register of the device associated with handle
function i2c_write_byte_data(p::Int32, handle::Int32, reg::Int, bVal::Int)
    return ccall((:i2c_write_byte_data, libpigpiod_if2), Int32, (Int32, UInt32, UInt32, UInt32), p, handle, reg, bVal)
end

#  i2c_read_word_data           smbus read word data
#  int i2c_read_word(int pi, unsigned handle, unsigned i2c_reg)
#  Reads single 16 bit word from specified register of device with associated handle
function i2c_read_word_data(p::Int32, handle::Int32, reg::Int)
    return ccall((:i2c_read_word_data, libpigpiod_if2), Int32, (Int32, UInt32, UInt32), p, handle, reg)
end

#  i2c_write_word_data          smbus write word data
#  int i2c_write_word_data(int pi, unsigned handle, unsigned i2c_reg, unsigned wVal)
#  Writes a single 16 bit word to specified register of device associated with handle
function i2c_write_word_data(p::Int32, handle::Int32, reg::Int, wVal::Int)
    return ccall((:i2c_write_word_data, libpigpiod_if2), Int32, (Int32, UInt32, UInt32, UInt32), p, handle, reg, wVal)
end

#  i2c_read_block_data          smbus read block data
#  int i2c_read_block_data(int pi, unsigned handle, unsigned i2c_reg, char *buf)
#  Reads block of up to 32 bytes from the specified register of the device associated with handle.
function i2c_read_block_data(p::Int32, handle::Int32, reg::Int, buf::Vector)
    return ccall((:i2c_read_block_data, libpigpiod_if2), Int32, (Int32, UInt32, UInt32, Ref{UInt8}), p, handle, reg, buf)
end

#  i2c_write_block_data         smbus write block data
#  int i2c_write_block_data(int pi, unsigned handle, unsigned i2c_reg, char *buf, unsigned count)
#  Writes up to 32 bytes to the specified register of the device associated with handle
function i2c_write_block_data(p::Int32, handle::Int32, reg::Int, buf::Vector, count::Int)
    return ccall((:i2c_write_block_data, libpigpiod_if2), Int32, (Int32, UInt32, UInt32, Ref{UInt8}, UInt32), p, handle, reg, buf, count)
end

#  i2c_read_i2c_block_data      smbus read I2C block data
#  int i2c_read_i2c_block_data(int pi, unsigned handle, unsigned i2c_reg, char *buf, unsigned count)
#  Reads count bytes from the specified register of the device associated with handle
function i2c_read_i2c_block_data(p::Int32, handle::Int32, reg::Int, buf::Vector, count::Int)
    return ccall((:i2c_read_i2c_block_data, libpigpiod_if2), Int32, (Int32, UInt32, UInt32, Ref{UInt8}, UInt32), p, handle, reg, buf, count)
end

#  i2c_write_i2c_block_data     smbus write I2C block data
#  int i2c_write_i2c_block_data(int pi, unsigned handle, unsigned i2c_reg, char *buf, unsigned count)
#  Writes count bytes from the specified register of the device associated with handle
function i2c_write_i2c_block_data(p::Int32, handle::Int32, reg::Int, buf::Vector, count::Int)
    return ccall((:i2c_write_i2c_block_data, libpigpiod_if2), Int32, (Int32, UInt32, UInt32, Ref{UInt8}, UInt32), p, handle, reg, buf, count)
end

#  i2c_read_device              Reads the raw I2C device
#  int i2c_read_device(int pi, unsigned handle, char *buf, unsigned count)
#  Reads count bytes from ray device into buf
function i2c_read_device(p::Int32, handle::Int32, buf::Vector, count::Int)
    return ccall((:i2c_read_device, libpigpiod_if2), Int32, (Int32, UInt32, Ref{UInt8}, UInt32), p, handle, buf, count)
end

#  i2c_write_device             Writes the raw I2C device
#  int i2c_write_device(int pi, unsigned handle, char *buf, unsigned count)
#  Writes count bytes from buf to raw device
function i2c_write_device(p::Int32, handle::Int32, buf::Vector, count::Int)
    return ccall((:i2c_write_device, libpigpiod_if2), Int32, (Int32, UInt32, Ref{UInt8}, UInt32), p, handle, buf, count)
end

#  i2c_process_call             smbus process call
#  int i2c_process_call(int pi, unsigned handle, unsigned i2c_reg, unsigned wVal)
#  Writes 2 byte word data to specified register of device with handle and reads 2 bytes of data in return
function i2c_process_call(p::Int32, handle::Int32, reg::Int, wVal::Int)
    return ccall((:i2c_process_call, libpigpiod_if2), Int32, (Int32, UInt32, UInt32, UInt32), p, handle, reg, wVal)
end

#  i2c_block_process_call       smbus block process call
#  int i2c_block_process_call(int pi, unsigned handle, unsigned i2c_reg, char *buf, unsigned count)
#  Writes data bytes to specific register of of device associated with handle and reads a device specified number of bytes in return
function i2c_block_process_call(p::Int32, handle::Int32, reg::Int, buf::Vector, count::Int)
    return ccall((:i2c_block_process_call, libpigpiod_if2), Int32, (Int32, UInt32, UInt32, Ref{UInt8}, UInt32), p, handle, reg, buf, count)
end

#  i2c_zip                      Performs multiple I2C transactions
#  int i2c_zip(int pi, unsigned handle, char *inBuf, unsigned inLen, char *outBuf, unsigned outLen)
#  This function executes a sequence of I2C operations. The operations to be performed are specified by the contents of inBuf which contains the concatenated command codes and associated data.
function i2c_zip(p::Int32, handle::Int32, inBuf::Vector, inLen::Int, outBuf::Vector, outLen::Int)
    return ccall((:i2c_zip, libpigpiod_if2), Int32, (Int32, UInt32, Ref{UInt8}, UInt32, Ref{UInt8}, UInt32), p, handle, inBuf, inLen, outBuf, outLen)
end

# I2C BIT BANG
#  bb_i2c_open          Opens GPIO for bit banging I2C
#  Sets a pair of GPIO for bit banging I2C at a specific baud rate
function bb_i2c_open(p::Int32, sda::Int, scl::Int, baud::Int)
    return ccall((:bb_i2c_open, libpigpiod_if2), Int32, (Int32, UInt32, UInt32, UInt32), p, sda, scl, baud)
end

#  bb_i2c_close         Closes GPIO for bit banging I2C
#  int bb_i2c_close(int pi, unsigned SDA)
#  Stops bit banging I2C on a pair of GPIO previously opened
function bb_i2c_close(p::Int32, sda::Int)
    return ccall((:bb_i2c_close, libpigpiod_if2), Int32, (Int32, UInt32), p, sda)
end

#  bb_i2c_zip           Performs bit banged I2C transactions
#  int bb_i2c_zip(int pi, unsigned SDA, char *inBuf, unsigned inLen, char *outBuf, unsigned outLen)
#  This function executes a sequence of bit banged I2C operations. The operations to be performed are specified by the contents of inBuf which contains the concatenated command codes and associated data.
function bb_i2c_zip(p::Int32, sda::Int, scl::Int, inBuf::Vector, inLen::Int, outBuf::Vector, outLen::Int)
    return ccall((:bb_i2c_zip, libpigpiod_if2), Int32, (Int32, UInt32, UInt32, Ref{UInt8}, UInt32, Ref{UInt8}, UInt32), p, sda, scl, inBuf, inLen, outBuf, outLen)
end

# I2C/SPI SLAVE
mutable struct bsc_xfer_t
    control::Int
    rxCnt::Int
    rxBuf::Vector{Int}
    txCnt::Int
    txBuf::Vector{Int}
end

#  bsc_xfer             I2C/SPI as slave transfer
#  int bsc_xfer(int pi, bsc_xfer_t *bscxfer)
#  This function provides a low-level interface to the SPI/I2C Slave peripheral on the BCM chip.
function bsc_xfer(p::Int32, bscxfer::bsc_xfer_t)
    return ccall((:bsc_xfer, libpigpiod_if2), Int32, (Int32, Ref{bsc_xfer_t}), p, bscxfer)
end

#  bsc_i2c              I2C as slave transfer
#  int bsc_i2c(int pi, int i2c_addr, bsc_xfer_t *bscxfer)
#  Allows Pi to act as a slace I2C device
function bsc_i2c(p::Int32, addr::Int32, bscxfer::bsc_xfer_t)
    return ccall((:bsc_i2c, libpigpiod_if2), Int32, (Int32, Int32, Ref{bsc_xfer_t}), p, addr, bscxfer)
end


text:
using Random

rand_config = bitrand(20)

for i in rand_config
    print(Int(i))
    sleep(0.05)
end

sleep(1)


text:
# using Pkg
# Pkg.activate("./")

using JLD2, FileIO
using Plots
using Printf
gr()

Ts = 1.2:0.1:3.4
n = size(Ts, 1)

ising_magn = Vector{Float64}(undef, n)
rbm_magn = Vector{Float64}(undef, n)
JLD2.@load joinpath("results", "ising_magnetization.jld2") ising_magn
JLD2.@load joinpath("results", "rbm_magnetization.jld2") rbm_magn

deviations = abs.(ising_magn .- rbm_magn)
display(deviations)

plot(Ts, rbm_magn, label = "RBM")
plot!(Ts, ising_magn, label = "Ising")
gui()

max_stars_repo_name: edwinb-ai/IsingRBM

max_issues_repo_name: edwinb-ai/IsingRBM

max_forks_repo_name: edwinb-ai/IsingRBM


text:
using Adapt
using CUDA: CuArray
using OffsetArrays: OffsetArray
using Oceananigans.Fields: fill_halo_regions!
using Oceananigans.Architectures: arch_array
using Oceananigans.BoundaryConditions: FBC

import Oceananigans.TurbulenceClosures: ivd_upper_diagonal,
                                        ivd_lower_diagonal

import Oceananigans.TurbulenceClosures: immersed_∂ⱼ_τ₁ⱼ,
                                        immersed_∂ⱼ_τ₂ⱼ,
                                        immersed_∂ⱼ_τ₃ⱼ,
                                        immersed_∇_dot_qᶜ

#####
##### Some conveniences for grid fitted boundaries
#####

abstract type AbstractGridFittedBoundary <: AbstractImmersedBoundary end
const GFIBG = ImmersedBoundaryGrid{<:Any, <:Any, <:Any, <:Any, <:Any, <:AbstractGridFittedBoundary}

#####
##### GridFittedBottom (2.5D immersed boundary with modified bottom height)
#####

"""
    GridFittedBottom(bottom)

Return an immersed boundary with an irregular bottom fit to the underlying grid.
"""
struct GridFittedBottom{B} <: AbstractGridFittedBoundary
    bottom_height :: B
end

"""
    ImmersedBoundaryGrid(grid, ib::GridFittedBottom)

Return a grid with `GridFittedBottom` immersed boundary.

Computes ib.bottom_height and wraps in an array.
"""
function ImmersedBoundaryGrid(grid, ib::GridFittedBottom)
    arch = grid.architecture
    bottom_field = Field{Center, Center, Nothing}(grid)
    set!(bottom_field, ib.bottom_height)
    fill_halo_regions!(bottom_field)
    offset_bottom_array = dropdims(bottom_field.data, dims=3)
    new_ib = GridFittedBottom(offset_bottom_array)
    return ImmersedBoundaryGrid(grid, new_ib)
end

function ImmersedBoundaryGrid(grid, ib::GridFittedBottom{<:OffsetArray})
    TX, TY, TZ = topology(grid)
    # TODO: check size
    return ImmersedBoundaryGrid{TX, TY, TZ}(grid, ib)
end

@inline function immersed_cell(i, j, k, underlying_grid, ib::GridFittedBottom)
    z = znode(c, c, c, i, j, k, underlying_grid)
    return @inbounds z < ib.bottom_height[i, j]
end

on_architecture(arch, ib::GridFittedBottom) = GridFittedBottom(arch_array(arch, ib.bottom_height))
Adapt.adapt_structure(to, ib::GridFittedBottom) = GridFittedBottom(adapt(to, ib.bottom_height))

#####
##### Implicit vertical diffusion
#####

####
#### For a center solver we have to check the interface "solidity" at faces k+1 in both the
#### Upper diagonal and the Lower diagonal
#### (because of tridiagonal convention where lower_diagonal on row k is found at k-1)
#### Same goes for the face solver, where we check at centers k in both Upper and lower diagonal
####

@inline immersed_ivd_peripheral_node(LX, LY, ::Center, i, j, k, ibg) = immersed_peripheral_node(LX, LY, Face(), i, j, k+1, ibg)
@inline immersed_ivd_peripheral_node(LX, LY, ::Face, i, j, k, ibg)   = immersed_peripheral_node(LX, LY, Center(), i, j, k, ibg)

# Extend the upper and lower diagonal functions of the batched tridiagonal solver

for location in (:upper_, :lower_)
    immersed_func = Symbol(:immersed_ivd_, location, :diagonal)
    ordinary_func = Symbol(:ivd_ ,         location, :diagonal)
    @eval begin
        # Disambiguation
        @inline $ordinary_func(i, j, k, ibg::GFIBG, closure, K, id, ℓx, ℓy, ℓz::Face, clock, Δt, κz) =
                $immersed_func(i, j, k, ibg::GFIBG, closure, K, id, ℓx, ℓy, ℓz, clock, Δt, κz)

        @inline $ordinary_func(i, j, k, ibg::GFIBG, closure, K, id, ℓx, ℓy, ℓz::Center, clock, Δt, κz) =
                $immersed_func(i, j, k, ibg::GFIBG, closure, K, id, ℓx, ℓy, ℓz, clock, Δt, κz)

        @inline function $immersed_func(i, j, k, ibg::GFIBG, closure, K, id, ℓx, ℓy, ℓz, clock, Δt, κz)
            return ifelse(immersed_ivd_peripheral_node(ℓx, ℓy, ℓz, i, j, k, ibg),
                          zero(eltype(ibg.underlying_grid)),
                          $ordinary_func(i, j, k, ibg.underlying_grid, closure, K, id, ℓx, ℓy, ℓz, clock, Δt, κz))
        end
    end
end

#####
##### GridFittedBoundary (experimental 3D immersed boundary)
#####

struct GridFittedBoundary{M} <: AbstractGridFittedBoundary
    mask :: M
end

@inline immersed_cell(i, j, k, underlying_grid, ib::GridFittedBoundary{<:AbstractArray}) = @inbounds ib.mask[i, j, k]

@inline function immersed_cell(i, j, k, underlying_grid, ib::GridFittedBoundary)
    x, y, z = node(c, c, c, i, j, k, underlying_grid)
    return ib.mask(x, y, z)
end

function compute_mask(grid, ib)
    arch = architecture(grid)
    mask_field = Field{Center, Center, Center}(grid, Bool)
    set!(mask_field, ib.mask)
    fill_halo_regions!(mask_field)
    return mask_field
end

function ImmersedBoundaryGrid(grid, ib::GridFittedBoundary; precompute_mask=true)
    TX, TY, TZ = topology(grid)

    if precompute_mask
        mask_field = compute_mask(grid, ib)
        new_ib = GridFittedBoundary(mask_field)
        return ImmersedBoundaryGrid{TX, TY, TZ}(grid, new_ib)
    else
        return ImmersedBoundaryGrid{TX, TY, TZ}(grid, ib)
    end
end

function ImmersedBoundaryGrid(grid, ib::GridFittedBoundary{<:OffsetArray}; kw...)
    TX, TY, TZ = topology(grid)
    return ImmersedBoundaryGrid{TX, TY, TZ}(grid, ib)
end

on_architecture(arch, ib::GridFittedBoundary{<:AbstractArray}) = GridFittedBoundary(arch_array(arch, ib.mask))
on_architecture(arch, ib::GridFittedBoundary{<:Field}) = GridFittedBoundary(compute_mask(on_architecture(arch, ib.mask.grid), ib))
on_architecture(arch, ib::GridFittedBoundary) = ib # need a workaround...

Adapt.adapt_structure(to, ib::GridFittedBoundary) = GridFittedBoundary(adapt(to, ib.mask))


text:
"""
    DiscreteVoronoi <: NeutralLandscapeMaker

    DiscreteVoronoi(; n=3)
    DiscreteVoronoi(n)

This type provides a rasterization of a Voronoi-like diagram.
Assigns a value to each patch using a 1-NN algorithmm with `n` initial clusters.
It is a `NearestNeighborElement` algorithmm with `k` neighbors set to 1.
The default is to use three clusters.
"""
@kwdef struct DiscreteVoronoi <: NeutralLandscapeMaker
    n::Int64 = 3
    function DiscreteVoronoi(n::Int64)
        @assert n > 0
        new(n)
    end
end

_landscape!(mat, alg::DiscreteVoronoi) = _landscape!(mat, NearestNeighborElement(alg.n, 1))


text:
"""
```
plot_scenario(m, var, class, scen; title = "", kwargs...)

plot_scenario(m, vars, class, scen; untrans = false, fourquarter = false,
    plotroot = figurespath(m, \"scenarios\"), titles = [], tick_size = 1,
    kwargs...)
```

Plot `var` or `vars` *in deviations from baseline* for the alternative scenario
specified by `key` and `vint`.

### Inputs

- `var::Symbol` or `vars::Vector{Symbol}`: variable(s) to be plotted,
  e.g. `:obs_gdp` or `[:obs_gdp, :obs_nominalrate]`
- `class::Symbol`
- `scen::AbstractScenario`: scenario

### Keyword Arguments

- `untrans::Bool`: whether to plot untransformed (model units) forecast
- `fourquarter::Bool`: whether to plot four-quarter forecast
- `plotroot::String`: if nonempty, plots will be saved in that directory
- `title::String` or `titles::Vector{String}`
- `tick_size::Int`: x-axis (time) tick size in units of years
- `legend`

See `?histforecast` for additional keyword arguments, all of which can be passed
into `plot_scenario`.

### Output

- `p::Plot` or `plots::OrderedDict{Symbol, Plot}`
"""
function plot_scenario(m::AbstractModel, var::Symbol, class::Symbol, scen::AbstractScenario;
                       title::String = "", kwargs...)

    plots = plot_scenario(m, [var], class, scen;
                          titles = isempty(title) ? String[] : [title],
                          kwargs...)
    return plots[var]
end

function plot_scenario(m::AbstractModel, vars::Vector{Symbol}, class::Symbol,
                       scen::AbstractScenario; untrans::Bool = false, fourquarter::Bool = false,
                       plotroot::String = figurespath(m, "scenarios"),
                       titles::Vector{String} = String[],
                       tick_size::Int = 1, legend = :none,
                       verbose::Symbol = :low,
                       kwargs...)
    # Determine output_var
    fcast_prod = if untrans && fourquarter
        error("Only one of untrans or fourquarter can be true")
    elseif untrans
        :forecastut
    elseif fourquarter
        :forecast4q
    else
        :forecast
    end

    # Read in MeansBands
    hist  = MeansBands()
    fcast = read_scenario_mb(m, scen, Symbol(fcast_prod, class))

    # Get titles if not provided
    if isempty(titles)
        detexify_title = typeof(Plots.backend()) == Plots.GRBackend
        titles = map(var -> describe_series(m, var, class, detexify = detexify_title), vars)
    end

    # Loop through variables
    plots = OrderedDict{Symbol, Plots.Plot}()
    for (var, title) in zip(vars, titles)
        ylabel = series_ylabel(m, var, class, untrans = untrans, fourquarter = fourquarter)
        ylabel = ylabel * " (deviations from baseline)"

        plots[var] = histforecast(var, hist, fcast;
                                  start_date = date_forecast_start(m),
                                  title = title, legend = legend, kwargs...)

        # Save if output_file provided
        output_file = if isempty(plotroot)
            ""
        else
            get_scenario_filename(m, scen, Symbol(fcast_prod, "_", detexify(var)),
                                  pathfcn = figurespath,
                                  fileformat = plot_extension())
        end
        save_plot(plots[var], output_file, verbose = verbose)
    end
    return plots
end


text:
# test AUROC compute functions
tprs = zeros(101)
fprs = zeros(101)

scores = collect(1:-0.01:0)
classes = [:b,:a,:a,:a,:b,:a,:a,:b,:b,:b,:a,:a,:b,:b,:b,:b,:a,:b,:a,:a,
    :a,:b,:a,:a,:b,:a,:b,:a,:a,:a,:b,:b,:a,:a,:a,:a,:a,:b,:b,:b,:a,:b,
    :a,:a,:b,:a,:a,:a,:a,:b,:a,:a,:b,:a,:a,:a,:a,:a,:b,:a,:b,:b,:a,:b,
    :a,:a,:a,:b,:b,:a,:a,:b,:b,:a,:b,:b,:b,:a,:b,:a,:b,:b,:a,:a,:b,:b,
    :b,:a,:b,:b,:b,:a,:a,:a,:a,:a,:a,:b,:a,:a,:b]
@test isapprox(auroc(scores, classes, Set([:a]))[1], 0.513556, atol=1e-6)
@test venn(scores, classes, Set([:a])) == 0.5357142857142857

scores = collect(1:-0.1:0)
classes = [:b, :b, :b, :b, :a, :a, :a, :a, :a, :a, :a]
@test auroc(scores, classes, Set([:a]))[1] == 0.0
@test venn(scores, classes, Set([:a])) == 0.0

scores = collect(1:-0.1:0)
classes = [:a, :a, :a, :a, :a, :a, :a, :a, :b, :b, :b]
@test auroc(scores, classes, Set([:a]))[1] == 1.0
@test venn(scores, classes, Set([:a])) == 1.0

# this algorithm is a direct implementation from the Boyd paper
function compute_auprc_posthoc(p, r)
    unqs = unique(r)
    pmins = p[[findlast(r, unq) for unq in unqs]]
    pmaxes = p[[findfirst(r, unq) for unq in unqs]]
    n = length(unqs)
    auprc = 0.0
    for i in n-1:-1:1
        auprc += (pmins[i] + pmaxes[i+1])/2 * (unqs[i+1] - unqs[i])
    end
    auprc
end

# Test multiple positive labels
classes = [:c,:a,:a,:b,:b,:b,:c,:c,:a,:c,:a,:a]
scores = Float64[80,5,157,169,158,166,115,18,6,78,31,5]
@test venn(scores, classes, Set([:a])) == 0.0
@test venn(scores, classes, Set([:b])) == 1.0
@test venn(scores, classes, Set([:a, :b])) == 1.0

@test auroc(scores, classes, Set([:a]))[1] == 0.14285714285714288
@test auroc(scores, classes, Set([:b]))[1] == 1.0
@test auroc(scores, classes, Set([:a, :b]))[1] == 0.53125

@test auprc(scores, classes, Set([:b]))[1] == 1.0
a, p, r = auprc(scores, classes, Set([:a, :b]))
@test a == compute_auprc_posthoc(p, r)
a, p, r = auprc(scores, classes, Set([:a]))
@test a == compute_auprc_posthoc(p, r)

# test for #30
classes = [:a, :a, :a, :b, :c, :a, :c, :c, :c, :b, :b]
scores = [-10., -8, -6, -4, -2, -1, 1, 1, 2, 3, 7]

@test venn(scores, classes, Set([:a]), rev=true) == 0.0
@test venn(scores, classes, Set([:a]), rev=false) == 1.0
@test venn(abs.(scores), classes, Set([:a, :b]), rev=true) == 1.0

@test auprc(scores, classes, Set([:a]), rev=true)[1] ≈ 0.21246843434343435
@test auprc(scores, classes, Set([:a]), rev=false)[1] ≈ 0.9083333333333333
@test auprc(abs.(scores), classes, Set([:a, :b]), rev=true)[1] ≈ 0.9662698412698412

function compute_bias(recalls, precisions, X, Y)
    true_auprc = 0.0
    for i in 2:length(recalls)
        true_auprc += (precisions[i-1] + precisions[i])/2*(recalls[i - 1] - recalls[i])
    end

    num_runs = 1000
    auprc_sum = 0.0
    cat_dist = Categorical([1-π, π])
    @inbounds for run in 1:num_runs
        cat = rand(cat_dist, 10^3)
        x, y = StatsBase.counts(cat, 1:2)
        scores = [rand(X, x); rand(Y, y)]
        auprc_sum += auprc(scores, [repmat(classes[1:1], x); repmat(classes[2:2], y)], Set([:b]))[1]
    end
    isapprox(auprc_sum/num_runs, true_auprc, atol=0.025)
end

π = 0.1
test_dists = Array[
     [Normal(0, 1), Normal(1,1)],
     [Beta(2, 5), Beta(5, 2)],
     [Uniform(0, 1), Uniform(0.5, 1.5)]
]
# true precision, recall functions
recall(xs, Y) = 1-cdf.(Y, xs)
function precision(xs, π, X, Y)
    res = π*recall(xs, Y)./(π*recall(xs, Y) + (1-π)*(1-cdf.(X, xs)))
    res[isnan.(res)] = 1
    res
end
xs = linspace(-10, 10, 1000)

for dists in test_dists
    X, Y = dists
    recalls = recall(xs, Y)
    precisions = precision(xs, π, X, Y)
    classes = [:a, :b]
    @test compute_bias(recalls, precisions, X, Y)
end


text:
@inline unorm(U) = sum(map((x->x.^2),U))

@inline function unpackfields(eqn::Euler{d},U) where {d}
    return first(U),ntuple(i->U[i+1],d),last(U)
end

"""
    prim_to_cons(eqn::Euler{d},rho,U,p) where {d}

convert primitive variables (ρ,U...,p) to conservative vars (ρ,ρU,E).
n-dimensional version where U = tuple(u1,...,u_d)
"""
function prim_to_cons(eqn::Euler{d},Q) where {d}
    rho,U,p = unpackfields(eqn,Q)
    rhoU = map(x->rho.*x,U)
    Unorm = unorm(U)
    E = @. p /(eqn.γ-1) + .5*rho*Unorm
    return SVector{d+2}(rho,rhoU...,E)
end

"""
    cons_to_prim_beta(eqn::Euler{d},U) where {d}

converts conservative variables to `primitive' variables which make evaluating EC
fluxes cheaper.
"""
@inline function cons_to_prim_beta(eqn::Euler{d},U) where {d}
    rho,rhoU,E = unpackfields(eqn,U)
    return SVector{d+2}(rho, map(x->x./rho,rhoU)..., betafun(eqn,U))
end

"""
    pfun(eqn::Euler{d},U) where {d}

Computes pressure (assuming ideal gas law) given array/tuple of conservative variables `U`
"""
@inline function pfun(eqn::Euler{d},U) where {d}
    rho,rhoU,E = unpackfields(eqn,U)
    rhoUnorm = unorm(rhoU)
    return @. (eqn.γ - 1)*(E - .5*rhoUnorm / rho)
end

"""
    betafun(eqn::Euler{d},U) where {d}

Converts "inverse temperature" β = ρ/(2p) given array/tuple of conservative variables `U`
"""
@inline function betafun(eqn::Euler{d},U) where {d}
    rho = first(U)
    p = pfun(eqn,U)
    return @. rho/(2*p)
end

function sfun(eqn::Euler{d},U) where {d}
    rho = first(U)
    p = pfun(eqn,U)
    return @. log(p/(rho^eqn.γ))
end

"""
    Sfun(eqn::Euler{d},U) where {d}

Converts mathematical entropy -ρs(U) given array/tuple of conservative variables `U`
"""
function Sfun(eqn::Euler{d},U) where {d}
    rho = first(U)
    return -rho.*sfun(eqn,U)
end

"""
    v_ufun(eqn::Euler{d}, U) where {d}

Returns entropy variables given tuple/array of conservative variables `U`.
"""
function v_ufun(eqn::Euler{d}, U) where {d}
    @unpack γ = eqn

    s = sfun(eqn,U)
    p = pfun(eqn,U)

    rho,rhoU,E = unpackfields(eqn,U)
    v1 = @. (γ+1-s) - (γ-1)*E/p
    vU = (x->@. x*(γ-1)/p).(rhoU)
    vE = (x->@. x*(γ-1)/p)(-rho)

    return SVector(v1,vU...,vE)
end

function s_vfun(eqn::Euler{d},V) where {d}
    @unpack γ = eqn
    v1,vU,vE = unpackfields(eqn,V)
    vUnorm = unorm(vU)
    return @. γ - v1 + vUnorm/(2*vE)
end

function rhoe_vfun(eqn::Euler{d},V) where {d}
    @unpack γ = eqn
    s = s_vfun(eqn,V)
    v1,vU,vE = unpackfields(eqn,V)
    return @. ((γ-1)/((-vE)^γ))^(1/(γ-1)) * exp(-s/(γ-1))
end

"""
    u_vfun(eqn::Euler{d},V) where {d}

Returns conservative variables given tuple/array of entropy variables `V`.
"""
function u_vfun(eqn::Euler{d},V) where {d}
    v1,vU,vE  = unpackfields(eqn,V)
    rhoeV     = rhoe_vfun(eqn,V)
    vUnorm    = unorm(vU)
    rho       = (@. rhoeV*(-vE))
    rhoU      = (x->rhoeV.*x).(vU)
    E         = (@. rhoeV*(1-vUnorm/(2*vE)))
    return SVector{d+2}(rho,rhoU...,E)
end

"""
    wavespeed(eqn::Euler{1},U)

Returns 1D wavespeed (used for flux penalization terms) given conservative variables U
"""
function wavespeed(eqn::Euler{1},U)
    rho,rhou,_ = unpackfields(eqn,U)
    p = pfun(eqn,U)
    cvel = @. sqrt(eqn.γ*p/rho)
    return @. abs(rhou/rho) + cvel
end

# function dUdV_explicit(v1,vU1,vU2,vE)
#     rho,rhou,rhov,E = u_vfun(v1,vU1,vU2,vE)
#     u,v = (x->x./rho).((rhou,rhov))
#     p = pfun(rho,rhou,rhov,E)
#     a2 = γ*p/rho
#     H = a2/(γ-1) + (u^2+v^2)/2
#
#     dUdV = @SMatrix [rho  rhou        rhov        E;
#                      rhou rhou*u + p  rhou*v      rhou*H;
#                      rhov rhov*u      rhov*v + p  rhov*H;
#                      E    rhou*H      rhov*H      rho*H^2-a2*p/(γ-1)]
#
#     return dUdV*(1/(γ-1))
# end
#
# function dVdU_explicit(rho,rhou,rhov,E)
#     rhoe = rhoe_ufun(rho,rhou,rhov,E)
#     V = v_ufun(rho,rhou,rhov,E)
#     k = .5*(V[2]^2+V[3]^2)/V[4]
#
#     dVdU = @SMatrix [γ+k^2      k*V[2]          k*V[3]         V[4]*(k+1);
#                     k*V[2]      V[2]^2-V[4]     V[2]*V[3]      V[2]*V[4];
#                     k*V[3]      V[2]*V[3]       V[3]^2-V[4]    V[3]*V[4]
#                     V[4]*(k+1)  V[2]*V[4]       V[3]*V[4]      V[4]^2]
#     return -dVdU/(rhoe*V[4])
# end


text:
using Revise
using Test

wait_for_key(prompt) = (println(stdout, prompt); read(stdin, 1); nothing)

names = [
    "test_operators_dense.jl",
    "test_superoperators.jl",
    "test_metrics.jl",
    "test_molecules.jl",
    "test_aggregate.jl",
    "test_evolution.jl",
    "test_schroedinger.jl",
    "test_liouville.jl",
    "test_master.jl",
]

names = ["test_master.jl"]

while true
    for name in names
        if startswith(name, "test_") && endswith(name, ".jl")
            try
                include(name)
            catch
                nothing
            end
        end
    end

    wait_for_key("press any key to run tests")
end

max_stars_repo_name: detrin/OpenQuantumSystems

max_issues_repo_name: detrin/OpenQuantumSystems