text:
# Specify model parameters and/or initial values to optimize
function get_search_index()::Tuple{Array{Int64,1},Array{Int64,1}}
    # parameters
    search_idx_params::Vector{Int} = [
        C.kassaikk
        C.kassaikknfkb
        C.kassanfkbikk
        C.kdisaikk
        C.kdisaikknfkb
        C.kdisanfkbikk
        C.kassa
        C.kdisa
        C.kdegbounda
        C.kdegboundaIKK
        C.kdegfreea
        C.kdegfreeaIKK
        C.kshutboundikbain
        C.kshutboundikbaout
        C.kshutfreeikbain
        C.kshutfreeikbaout
        C.k0mrnaikba
        C.kdegmrnaikba
        C.khillprodmrnaikba
        C.kpikba
        C.kprodmrnaikba

        C.kassbikk
        C.kassbikknfkb
        C.kassbnfkbikk
        C.kdisbikk
        C.kdisbikknfkb
        C.kdisbnfkbikk
        C.kassb
        C.kdegboundb
        C.kdegboundbIKK
        C.kdegfreeb
        C.kdegfreebIKK
        C.kdisb
        C.kshutboundikbbin
        C.kshutboundikbbout
        C.kshutfreeikbbin
        C.kshutfreeikbbout
        C.k0mrnaikbb
        C.kdegmrnaikbb
        # C.khillprodmrnaikbb
        C.kpikbb
        # C.kprodmrnaikbb
        C.kasseikk
        C.kasseikknfkb
        C.kassenfkbikk
        C.kdiseikk
        C.kdiseikknfkb
        C.kdisenfkbikk

        C.kasse
        C.kdegbounde
        C.kdegboundeIKK
        C.kdegfreee
        C.kdegfreeeIKK
        C.kdise
        C.kshutboundikbein
        C.kshutboundikbeout
        C.kshutfreeikbein
        C.kshutfreeikbeout
        C.k0mrnaikbe
        C.kdegmrnaikbe
        C.khillprodmrnaikbe
        C.kpikbe
        C.kprodmrnaikbe

        C.kshutnfkbin
        C.kshutnfkbout

        C.k0mrnaa20
        C.kdegmrnaa20
        C.kprodmrnaa20
        C.kpa20
        C.kdega20
        C.khillprodmrnaa20
    ]

    # initial values
    search_idx_initials::Vector{Int} = [
        # V.(variableName)
    ]

    return search_idx_params, search_idx_initials
end

function get_search_region()::Matrix{Float64}
    p::Vector{Float64} = param_values()
    u0::Vector{Float64} = initial_values()

    search_idx::Tuple{Array{Int64,1},Array{Int64,1}} = get_search_index()
    search_param::Vector{Float64} = initialize_search_param(search_idx, p, u0)

    search_rgn::Matrix{Float64} = zeros(2, length(p) + length(u0))

    # Default: 0.1 ~ 10x
    for (i, j) in enumerate(search_idx[1])
        search_rgn[1,j] = search_param[i] * 0.1  # lower bound
        search_rgn[2,j] = search_param[i] * 10.0  # upper bound
    end

    # Default: 0.5 ~ 2x
    for (i, j) in enumerate(search_idx[2])
        search_rgn[1,j + length(p)] = search_param[i + length(search_idx[1])] * 0.5  # lower bound
        search_rgn[2,j + length(p)] = search_param[i + length(search_idx[1])] * 2.0  # upper bound
    end

    # search_rgn[:, C.param_name] = [lower_bound, upper_bound]
    # search_rgn[:, V.var_name+length(p)] = [lower_bound, upper_bound]
    search_rgn[:, C.khillprodmrnaikba] = [1.0, 2.0]
    search_rgn[:, C.khillprodmrnaikbe] = [1.0, 2.0]
    search_rgn[:, C.khillprodmrnaa20] = [1.0, 2.0]

    search_rgn = convert_scale!(search_rgn, search_idx)

    return search_rgn
end

function update_param(indiv::Vector{Float64})::Tuple{Array{Float64,1},Array{Float64,1}}
    p::Vector{Float64} = param_values()
    u0::Vector{Float64} = initial_values()

    search_idx::Tuple{Array{Int64,1},Array{Int64,1}} = get_search_index()

    for (i, j) in enumerate(search_idx[1])
        @inbounds p[j] = indiv[i]
    end
    for (i, j) in enumerate(search_idx[2])
        @inbounds u0[j] = indiv[i + length(search_idx[1])]
    end

    return p, u0
end

function decode_gene2val(indiv_gene)::Vector{Float64}
    search_rgn::Matrix{Float64} = get_search_region()
    indiv::Vector{Float64} = zeros(length(indiv_gene))

    for i in eachindex(indiv_gene)
        indiv[i] = 10^(
            indiv_gene[i] * (
                search_rgn[2,i] - search_rgn[1,i]
            ) + search_rgn[1,i]
        )
    end

    return round.(indiv, sigdigits=7)
end

function encode_val2gene(indiv::Vector{Float64})
    search_rgn::Matrix{Float64} = get_search_region()
    indiv_gene::Vector{Float64} = zeros(length(indiv))

    for i in eachindex(indiv)
        indiv_gene[i] = (
            log10(indiv[i]) - search_rgn[1,i]
        ) / (
            search_rgn[2,i] - search_rgn[1,i]
        )
    end

    return indiv_gene
end

function encode_bestIndivVal2randGene(
        gene_idx::Int64,
        best_indiv::Vector{Float64},
        p0_bounds::Vector{Float64})::Float64
    search_rgn::Matrix{Float64} = get_search_region()
    rand_gene::Float64 = (
        log10(
            best_indiv[gene_idx] * 10^(
                rand() * log10(p0_bounds[2] / p0_bounds[1]) + log10(p0_bounds[1])
            )
        ) - search_rgn[1,gene_idx]
    ) / (
        search_rgn[2,gene_idx] - search_rgn[1,gene_idx]
    )
    return rand_gene
end

function initialize_search_param(
        search_idx::Tuple{Array{Int64,1},Array{Int64,1}},
        p::Vector{Float64},
        u0::Vector{Float64})::Vector{Float64}
    duplicate::Vector{String} = []
    if length(search_idx[1]) != length(unique(search_idx[1]))
        for idx in findall(
            [count(x -> x == i, search_idx[1]) for i in unique(search_idx[1])] .!= 1
        )
            push!(duplicate, C.NAMES[search_idx[1][idx]])
        end
        error(
            "Duplicate parameters (C.): $duplicate"
        )
    elseif length(search_idx[2]) != length(unique(search_idx[2]))
        for idx in findall(
            [count(x -> x == i, search_idx[2]) for i in unique(search_idx[2])] .!= 1
        )
            push!(duplicate, V.NAMES[search_idx[2][idx]])
        end
        error(
            "Duplicate initial conditions (V.): $duplicate"
        )
    end
    search_param = zeros(
        length(search_idx[1]) + length(search_idx[2])
    )
    for (i, j) in enumerate(search_idx[1])
        @inbounds search_param[i] = p[j]
    end
    for (i, j) in enumerate(search_idx[2])
        @inbounds search_param[i + length(search_idx[1])] = u0[j]
    end

    if any(x -> x == 0.0, search_param)
        msg::String = "search_param must not contain zero."
        for idx in search_idx[1]
            if p[idx] == 0.0
                error(
                    @sprintf(
                        "`C.%s` in search_idx_params: ", C.NAMES[idx]
                    ) * msg
                )
            end
        end
        for idx in search_idx[2]
            if u0[idx] == 0.0
                error(
                    @sprintf(
                        "`V.%s` in search_idx_initials: ", V.NAMES[idx]
                    ) * msg
                )
            end
        end
    end

    return search_param
end

function convert_scale!(
        search_rgn::Matrix{Float64},
        search_idx::Tuple{Array{Int64,1},Array{Int64,1}})::Matrix{Float64}
    for i = 1:size(search_rgn, 2)
        if minimum(search_rgn[:,i]) < 0.0
            msg = "search_rgn[lower_bound,upper_bound] must be positive.\n"
            if i <= C.NUM
                error(@sprintf("`C.%s` ", C.NAMES[i]) * msg)
            else
                error(@sprintf("`V.%s` ", V.NAMES[i - C.NUM]) * msg)
            end
        elseif minimum(search_rgn[:,i]) == 0.0 && maximum(search_rgn[:,i]) != 0.0
            msg = "lower_bound must be larger than 0.\n"
            if i <= C.NUM
                error(@sprintf("`C.%s` ", C.NAMES[i]) * msg)
            else
                error(@sprintf("`V.%s` ", V.NAMES[i - C.NUM]) * msg)
            end
        elseif search_rgn[2,i] - search_rgn[1,i] < 0.0
            msg = "lower_bound must be smaller than upper_bound.\n"
            if i <= C.NUM
                error(@sprintf("`C.%s` ", C.NAMES[i]) * msg)
            else
                error(@sprintf("`V.%s` ", V.NAMES[i - C.NUM]) * msg)
            end
        end
    end

    nonzero_idx::Vector{Int} = []
    for i = 1:size(search_rgn, 2)
        if search_rgn[:,i] != [0.0,0.0]
            push!(nonzero_idx, i)
        end
    end
    difference::Vector{Int} = collect(
        symdiff(
            Set(nonzero_idx),
            Set(append!(search_idx[1], C.NUM .+ search_idx[2]))
        )
    )
    if length(difference) > 0
        for idx in difference
            if j <= C.NUM
                println(@sprintf("`C.%s`", C.NAMES[Int(idx)]))
            else
                println(@sprintf("`V.%s`", V.NAMES[Int(idx) - C.NUM]))
            end
        end
        error(
            "Set these search_params in both search_idx and search_rgn."
        )
    end

    search_rgn = search_rgn[:,nonzero_idx]

    return log10.(search_rgn)
end


text:
# Simulate an Levy driven SDE

using Bridge, PyPlot
Random.seed!(5.)

# Define Gamma process
γ = 10.
α = 5.
P = GammaProcess(γ, α)

# Define Grid

T = 5.0
tt = 0.0:T/n:T
m = 1000
sigma = (t, x = 0.0) -> exp(-cos(t))

# Driving Levy process is a difference of Gamma processes
L = sample(tt, P) .- sample(tt, P)

# Solve Levy driven SDE
X = solve(EulerMaruyama(), 0.0, L, ((t,x)-> -0.*x, sigma))

# Plot result
figure();plot(X.tt, X.yy, label="X")
plot(X.tt, L.yy, label="L")
plot(X.tt, sigma.(X.tt), label=L"\sigma")
legend()
PyPlot.savefig("output/levysde.png")

ext: jl


text:
"""
    mlm_perms(data, nPerms; permFun, hasXIntercept, hasZIntercept,
              weights, targetType, isMainEff)

Obtains permutation p-values for MLM t-statistics.

# Arguments

- data = RawData object
- nPerms = number of permutations. Defaults to `1000`.

# Keyword arguments

- permFun = function used to permute `Y`. Defaults to `shuffle_rows`
  (shuffles rows of `Y`).
- hasXIntercept = boolean flag indicating whether or not to include an `X`
  intercept (row main effects). Defaults to `true`.
- hasZIntercept = boolean flag indicating whether or not to include a `Z`
  intercept (column main effects). Defaults to `true`.
- weights = 1d array of floats to use as column weights for `Y`, or `nothing`.
  If the former, must be the same length as the number of columns of `Y`.
  Defaults to `nothing`.
- targetType = string indicating the target type toward which to shrink the
  error variance, or `nothing`. If the former, acceptable inputs are "A", "B",
  "C", and "D". Defaults to `nothing`.
    - "A": Target is identity matrix
    - "B": Target is diagonal matrix with constant diagonal
    - "C": Target is has same diagonal element, and same off-diagonal element
    - "D": Target is diagonal matrix with unequal entries
- isMainEff = boolean flag indicating whether or not to include p-values for
  the main effects

# Value

Tuple
- tStats: 2d array of floats; t-statistics
- pvals: 2d array of floats; permutation p-values

# Some notes

Permutations are computed in parallel when possible.

"""
function mlm_perms(data::RawData, nPerms::Int64=1000;
                   permFun::Function=shuffle_rows,
                   hasXIntercept::Bool=true, hasZIntercept::Bool=true,
                   weights=nothing, targetType=nothing, isMainEff::Bool=false)

    # Wrapper function that performs MLM and gets t-statistics
    function mlm_t_stat(data::RawData)
        return t_stat(mlm(data; hasXIntercept=hasXIntercept,
                      hasZIntercept=hasZIntercept, weights=weights,
                      targetType=targetType), isMainEff)
    end

    # Get MLM t-statistcs and permutation p-values
    return perm_pvals(mlm_t_stat, data, nPerms; permFun=permFun)
end


text:
using HWfuncapp
using Test

@testset "HWfuncapp.jl" begin
	# test that q1 with 15 nodes has an error smaller than 1e-9

    # test that the integral of function h [-10,0] ≈ 1.46039789878568

    @test false  # by default fail
end


text:
__precompile__()

module ImageCore

using Colors, FixedPointNumbers, MappedArrays, Graphics
using Colors: Fractional

using Base: tail, @pure, Indices

import FixedPointNumbers: ufixed8, ufixed10, ufixed12, ufixed14, ufixed16
import Graphics: width, height

export
    ## Types
    ChannelView,
    ColorView,
    ## functions
    # views
    channelview,
    colorview,
    permuteddimsview,
    rawview,
    ufixedview,
    # conversions
#    float16,
    float32,
    float64,
    u8,
    ufixed8,
    ufixed10,
    ufixed12,
    ufixed14,
    ufixed16,
    u16,
    # traits
    assert_timedim_last,
    coords_spatial,
    height,
    indices_spatial,
    nimages,
    pixelspacing,
    sdims,
    size_spatial,
    spacedirections,
    width,
    widthheight

include("colorchannels.jl")
include("convert_reinterpret.jl")
include("traits.jl")
include("deprecated.jl")

"""
    rawview(img::AbstractArray{FixedPoint})

returns a "view" of `img` where the values are interpreted in terms of
their raw underlying storage. For example, if `img` is an `Array{U8}`,
the view will act like an `Array{UInt8}`.
"""
rawview{T<:FixedPoint}(a::AbstractArray{T}) = mappedarray((x->x.i, y->T(y,0)), a)
rawview{T<:Real}(a::AbstractArray{T}) = a

"""
    ufixedview([T], img::AbstractArray{Unsigned})

returns a "view" of `img` where the values are interpreted in terms of
`UFixed` number types. For example, if `img` is an `Array{UInt8}`, the
view will act like an `Array{UFixed8}`.  Supply `T` if the element
type of `img` is `UInt16`, to specify whether you want a `UFixed10`,
`UFixed12`, `UFixed14`, or `UFixed16` result.
"""
ufixedview{T<:FixedPoint,S<:Unsigned}(::Type{T}, a::AbstractArray{S}) = mappedarray((y->T(y,0),x->x.i), a)
ufixedview(a::AbstractArray{UInt8}) = ufixedview(U8, a)
ufixedview{T<:UFixed}(::Type{T}, a::AbstractArray{T}) = a

"""
    permuteddimsview(A, perm)

returns a "view" of `A` with its dimensions permuted as specified by
`perm`. This is like `permutedims`, except that it produces a view
rather than a copy of `A`; consequently, any manipulations you make to
the output will be mirrored in `A`. Compared to the copy, the view is
much faster to create, but generally slower to use.
"""
permuteddimsview(A, perm) = Base.PermutedDimsArrays.PermutedDimsArray(A, perm)

end ## module


text:
macro OK(func)
    return :($func == SQLITE_OK)
end

macro CHECK(db,ex)
    return quote
        if !(@OK $ex)
            sqliteerror($db)
        end
        SQLITE_OK
    end
end

macro NULLCHECK(ptr)
    return quote
        if $ptr == C_NULL
            error()
        end
    end
end

regexp(r::AbstractString, s::AbstractString) = ismatch(Regex(r), s)

sqliteopen(file,handle) = sqlite3_open(file, handle)
sqliteopen(file::UTF16String,handle) = sqlite3_open16(file, handle)

type SQLiteException <: Exception
    msg::AbstractString
end
sqliteerror() = throw(SQLiteException(bytestring(sqlite3_errmsg())))
sqliteerror(db) = throw(SQLiteException(bytestring(sqlite3_errmsg(db.handle))))


text:
#using CuArrays
using Statistics: mean
using Tables
using CategoricalArrays
using MappedArrays: mappedarray

export find_best_split

"""
    find_best_split(loss, df::DataFrameLike, feature, target, warmstart, lambda, gamma)

Find the best (binary) split point by optimizing ∑ loss(warmstart + δx, target) using order-2
Taylor-series expexpansion.

Does not assume that Feature, target, and warmstart are sorted and will sort them for you.
"""

function find_best_split(loss, df, feature::Symbol, target::Symbol, warmstart::AbstractVector, lambda, gamma; verbose = false, kwargs...)
	 @assert Tables.istable(df)

	 if verbose
	     println("Choosing a split on ", feature)
	 end

	 dfc = Tables.columns(df)

	 x = getproperty(dfc, feature)
	 target_vec = getproperty(dfc, target);

	 split_res = find_best_split(loss, x, target_vec, warmstart, lambda, gamma; verbose = verbose, kwargs...)
	 (feature = feature, split_res...)
end

"""
    find_best_split(loss, feature, target, warmstart, lambda, gamma)

Find the best (binary) split point by optimizing ∑ loss(warmstart + δx, target) using order-2 Taylor series expexpansion.

Does not assume that Feature, target, and warmstart sorted and will sort them for you.
"""
function find_best_split(loss, features::AbstractVector, target::AbstractVector, warmstart::AbstractVector, lambda::Number, gamma::Number; kwargs...)
	@assert length(features) == length(target)
	@assert length(features) == length(warmstart)

    if issorted(features)
        res = _find_best_split(loss, features, target, warmstart, lambda, gamma; kwargs...)
    else
        s = fsortperm(features)
        res = _find_best_split(loss, @view(features[s]), @view(target[s]), @view(warmstart[s]), lambda, gamma; kwargs...)
    end
end

"""
	_find_best_split(fn, f, t, p, lambda, gamma, verbose)

Assume that f, t, p are iterable and that they are sorted. Intended for advanced users only
"""
function _find_best_split(loss::LogitLogLoss, feature::AbstractVector, target::CategoricalVector, warmstart::AbstractVector, lambda::Number, gamma::Number; kwargs...)
	@assert length(levels(target)) == 2

	find_best_split(loss, feature, 2 .- target.refs, warmstart, lambda, gamma; kwargs...)
end

function _find_best_split(loss::LogitLogLoss, feature::AbstractVector, target::SubArray{A, B, C, D, E}, warmstart::AbstractVector, lambda::Number, gamma::Number; kwargs...) where {A, B, C<:CategoricalArray, D, E}
	@assert length(levels(target)) == 2

	find_best_split(loss, feature, mappedarray(x->2 - x.level, target), warmstart, lambda, gamma; kwargs...)
end

function _find_best_split(loss, feature, target, warmstart, lambda::Number, gamma::Number; min_child_weight = 1, verbose = false)
    @assert length(feature) >= 2
	@assert length(target) == length(feature)
	@assert length(warmstart) == length(feature)

	cg = cumsum(g.(loss, target, warmstart))
    ch = cumsum(h.(loss, target, warmstart))

    max_cg = cg[end]
    max_ch = ch[end]

    last_feature = feature[1]
    cutpt = zero(Int)
    lweight = 0.0
    rweight = 0.0
    best_gain = typemin(Float64)

	for (i, (f, cg, ch)) in enumerate(zip(drop(feature,1) , @view(cg[1:end-1]), @view(ch[1:end-1])))
		if f != last_feature
			left_split = cg^2 /(ch + lambda)
			right_split = (max_cg-cg)^(2) / ((max_ch - ch) + lambda)
			no_split = max_cg^2 /(max_ch + lambda)
			gain = left_split +  right_split - no_split - gamma
			if gain > best_gain
				best_gain = gain
				cutpt = i
				lweight = -cg/(ch+lambda)
				rweight = -(max_cg - cg)/(max_ch - ch + lambda)
			end
			last_feature = f
		end
	end
	# end

	# set the split at the point at the end
    split_at = feature[end]
    further_split = false

	# the child weight is the hessian
    if cutpt >= 1
    	split_at = feature[cutpt]
		if ch[cutpt] < min_child_weight
            # the weight is less than child weight; do not split further
            # TODO assess if this is appropriate
			split_at = feature[end]
			cutpt = 0
			no_split = max_cg^2 /(max_ch + lambda)
	    	gain = no_split - gamma
	    	lweight = -cg[end]/(ch[end]+lambda)
            rweight = -cg[end]/(ch[end]+lambda)
            further_split = false
        else
            further_split = true
		end
    end

    (split_at = split_at, cutpt = cutpt, gain = best_gain, lweight = lweight, rweight = rweight, further_split = further_split)
end

# TODO more reseach into GPU friendliness
# function _find_best_split(loss, feature, target::CuArray, warmstart::CuArray, lambda::Number, gamma::Number; verbose = false)
# 	g1 = g.(loss, target, warmstart)
# 	h1 = h.(loss, target, warmstart)
#
# 	cg = cumsum(g1)
# 	ch = cumsum(h1)
#
# 	max_cg = sum(g1)
# 	max_ch = sum(h1)
#
# 	lambda = 0.0
# 	gamma = 0.0
# 	left_split = cumsum(g1).^2 ./ (cumsum(h1) .+ lambda)
# 	right_split = (max_cg .- cumsum(g1)).^(2) ./ ((max_ch .- cumsum(h1)) .+ lambda)
# 	no_split = max_cg^2 / (max_ch + lambda)
# 	gain = left_split .+  right_split .- no_split .- gamma
#
# 	# find the positions where values change because it's only meaningful to cut
# 	# at where the values change
# 	i = findall(!=(0), diff(feature))
#
# 	split_at, cutpt_i = findmax(gain[i])
# 	cutpt = i[cutpt_i]
# 	best_gain = gain[cutpt]
# 	lweight = -cg[cutpt] / (ch[cutpt] + lambda)
# 	rweight = -(max_cg - cg[cutpt])/(max_ch - ch[cutpt] + lambda)
#
# 	(split_at = split_at, cutpt = cutpt, gain = best_gain, lweight = lweight, rweight = rweight)
# end


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

mutable struct IoK8sApiCoreV1ServiceSpec <: SwaggerModel
    clusterIP::Any # spec type: Union{ Nothing, String } # spec name: clusterIP
    externalIPs::Any # spec type: Union{ Nothing, Vector{String} } # spec name: externalIPs
    externalName::Any # spec type: Union{ Nothing, String } # spec name: externalName
    externalTrafficPolicy::Any # spec type: Union{ Nothing, String } # spec name: externalTrafficPolicy
    healthCheckNodePort::Any # spec type: Union{ Nothing, Int32 } # spec name: healthCheckNodePort
    ipFamily::Any # spec type: Union{ Nothing, String } # spec name: ipFamily
    loadBalancerIP::Any # spec type: Union{ Nothing, String } # spec name: loadBalancerIP
    loadBalancerSourceRanges::Any # spec type: Union{ Nothing, Vector{String} } # spec name: loadBalancerSourceRanges
    ports::Any # spec type: Union{ Nothing, Vector{IoK8sApiCoreV1ServicePort} } # spec name: ports
    publishNotReadyAddresses::Any # spec type: Union{ Nothing, Bool } # spec name: publishNotReadyAddresses
    selector::Any # spec type: Union{ Nothing, Dict{String, String} } # spec name: selector
    sessionAffinity::Any # spec type: Union{ Nothing, String } # spec name: sessionAffinity
    sessionAffinityConfig::Any # spec type: Union{ Nothing, IoK8sApiCoreV1SessionAffinityConfig } # spec name: sessionAffinityConfig
    topologyKeys::Any # spec type: Union{ Nothing, Vector{String} } # spec name: topologyKeys
    type::Any # spec type: Union{ Nothing, String } # spec name: type

    function IoK8sApiCoreV1ServiceSpec(;clusterIP=nothing, externalIPs=nothing, externalName=nothing, externalTrafficPolicy=nothing, healthCheckNodePort=nothing, ipFamily=nothing, loadBalancerIP=nothing, loadBalancerSourceRanges=nothing, ports=nothing, publishNotReadyAddresses=nothing, selector=nothing, sessionAffinity=nothing, sessionAffinityConfig=nothing, topologyKeys=nothing, type=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1ServiceSpec, Symbol("clusterIP"), clusterIP)
        setfield!(o, Symbol("clusterIP"), clusterIP)
        validate_property(IoK8sApiCoreV1ServiceSpec, Symbol("externalIPs"), externalIPs)
        setfield!(o, Symbol("externalIPs"), externalIPs)
        validate_property(IoK8sApiCoreV1ServiceSpec, Symbol("externalName"), externalName)
        setfield!(o, Symbol("externalName"), externalName)
        validate_property(IoK8sApiCoreV1ServiceSpec, Symbol("externalTrafficPolicy"), externalTrafficPolicy)
        setfield!(o, Symbol("externalTrafficPolicy"), externalTrafficPolicy)
        validate_property(IoK8sApiCoreV1ServiceSpec, Symbol("healthCheckNodePort"), healthCheckNodePort)
        setfield!(o, Symbol("healthCheckNodePort"), healthCheckNodePort)
        validate_property(IoK8sApiCoreV1ServiceSpec, Symbol("ipFamily"), ipFamily)
        setfield!(o, Symbol("ipFamily"), ipFamily)
        validate_property(IoK8sApiCoreV1ServiceSpec, Symbol("loadBalancerIP"), loadBalancerIP)
        setfield!(o, Symbol("loadBalancerIP"), loadBalancerIP)
        validate_property(IoK8sApiCoreV1ServiceSpec, Symbol("loadBalancerSourceRanges"), loadBalancerSourceRanges)
        setfield!(o, Symbol("loadBalancerSourceRanges"), loadBalancerSourceRanges)
        validate_property(IoK8sApiCoreV1ServiceSpec, Symbol("ports"), ports)
        setfield!(o, Symbol("ports"), ports)
        validate_property(IoK8sApiCoreV1ServiceSpec, Symbol("publishNotReadyAddresses"), publishNotReadyAddresses)
        setfield!(o, Symbol("publishNotReadyAddresses"), publishNotReadyAddresses)
        validate_property(IoK8sApiCoreV1ServiceSpec, Symbol("selector"), selector)
        setfield!(o, Symbol("selector"), selector)
        validate_property(IoK8sApiCoreV1ServiceSpec, Symbol("sessionAffinity"), sessionAffinity)
        setfield!(o, Symbol("sessionAffinity"), sessionAffinity)
        validate_property(IoK8sApiCoreV1ServiceSpec, Symbol("sessionAffinityConfig"), sessionAffinityConfig)
        setfield!(o, Symbol("sessionAffinityConfig"), sessionAffinityConfig)
        validate_property(IoK8sApiCoreV1ServiceSpec, Symbol("topologyKeys"), topologyKeys)
        setfield!(o, Symbol("topologyKeys"), topologyKeys)
        validate_property(IoK8sApiCoreV1ServiceSpec, Symbol("type"), type)
        setfield!(o, Symbol("type"), type)
        o
    end
end # type IoK8sApiCoreV1ServiceSpec

const _property_map_IoK8sApiCoreV1ServiceSpec = Dict{Symbol,Symbol}(Symbol("clusterIP")=>Symbol("clusterIP"), Symbol("externalIPs")=>Symbol("externalIPs"), Symbol("externalName")=>Symbol("externalName"), Symbol("externalTrafficPolicy")=>Symbol("externalTrafficPolicy"), Symbol("healthCheckNodePort")=>Symbol("healthCheckNodePort"), Symbol("ipFamily")=>Symbol("ipFamily"), Symbol("loadBalancerIP")=>Symbol("loadBalancerIP"), Symbol("loadBalancerSourceRanges")=>Symbol("loadBalancerSourceRanges"), Symbol("ports")=>Symbol("ports"), Symbol("publishNotReadyAddresses")=>Symbol("publishNotReadyAddresses"), Symbol("selector")=>Symbol("selector"), Symbol("sessionAffinity")=>Symbol("sessionAffinity"), Symbol("sessionAffinityConfig")=>Symbol("sessionAffinityConfig"), Symbol("topologyKeys")=>Symbol("topologyKeys"), Symbol("type")=>Symbol("type"))
const _property_types_IoK8sApiCoreV1ServiceSpec = Dict{Symbol,String}(Symbol("clusterIP")=>"String", Symbol("externalIPs")=>"Vector{String}", Symbol("externalName")=>"String", Symbol("externalTrafficPolicy")=>"String", Symbol("healthCheckNodePort")=>"Int32", Symbol("ipFamily")=>"String", Symbol("loadBalancerIP")=>"String", Symbol("loadBalancerSourceRanges")=>"Vector{String}", Symbol("ports")=>"Vector{IoK8sApiCoreV1ServicePort}", Symbol("publishNotReadyAddresses")=>"Bool", Symbol("selector")=>"Dict{String, String}", Symbol("sessionAffinity")=>"String", Symbol("sessionAffinityConfig")=>"IoK8sApiCoreV1SessionAffinityConfig", Symbol("topologyKeys")=>"Vector{String}", Symbol("type")=>"String")
Base.propertynames(::Type{ IoK8sApiCoreV1ServiceSpec }) = collect(keys(_property_map_IoK8sApiCoreV1ServiceSpec))
Swagger.property_type(::Type{ IoK8sApiCoreV1ServiceSpec }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiCoreV1ServiceSpec[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1ServiceSpec }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1ServiceSpec[property_name]

function check_required(o::IoK8sApiCoreV1ServiceSpec)
    true
end

function validate_property(::Type{ IoK8sApiCoreV1ServiceSpec }, name::Symbol, val)
end


text:
function poi2videofile(folder::String, poi::POI, i::Int)
    name = joinpath(folder,"POI_videos", "$i.mp4")
    isfile(name) && return name
    if poi.start.file ≠ poi.stop.file
        file = tempname()
        open(file, "w") do o
            for f in (poi.start.file, poi.stop.file)
                fullname = joinpath(folder, f)
                println(o, "file $fullname")
            end
        end
        Δ = duration(poi, folder)
        run(`ffmpeg -f concat -safe 0 -i $file -c copy -ss $(Dates.value(poi.start.time)) -to $Δ $name`)
    else
        fullname = joinpath(folder, poi.start.file)
        if poi.start.time == poi.stop.time
            run(`ffmpeg -i $fullname -c copy -ss $(Dates.value(poi.start.time)) -to $(Dates.value(poi.stop.time) + 1) $name`)
        else
            run(`ffmpeg -i $fullname -c copy -ss $(Dates.value(poi.start.time)) -to $(Dates.value(poi.stop.time)) $name`)
        end
    end
    return name
end

function fragment(folder::String)
    a = loadAssociation(folder)
    poi_folder = joinpath(folder,"POI_videos")
    isdir(poi_folder) || mkdir(poi_folder)
    for (i,p) in enumerate(a.pois)
        poi2videofile(folder, p, i)
    end
end

function coordinates_gui(folder::String)

    a = loadAssociation(folder)
    poi2do = [p for (i, p) in enumerate(a.pois) if isfile(joinpath(folder, "log", "$i.h5"))]
    isempty(poi2do) && return
    builder = Builder(filename=joinpath(@__DIR__, "choose2track.glade"))

    # data
    n = length(poi2do)
    # widgets
    done = button(widget=builder["done"])
    previous = button(widget=builder["previous"])
    next = button(widget=builder["next"])
    track = button(widget=builder["track"])
    # functions
    down = map(_ -> -1, previous)
    up = map(_ -> +1, next)
    step = merge(down, up)
    _state = foldp(1, step) do x,y
        clamp(x + y, 1, n)
    end
    state = droprepeats(_state)
    pb = progressbar(n, widget=builder["progressbar"], signal=state)
    poi = map(state) do i
        poi2do[i]
    end
    poilabel = map(p -> "POI: $(p.name), $(p.label)", poi)
    label(value(poilabel), widget=builder["poi"], signal=poilabel)
    runlabel = map(poi) do p
        r = String[]
        for (pp, rr) in a.associations
            if p == pp
                push!(r, """
                      Run
                      Repetition: $(rr.repetition)""")
                for (k, v) in rr.run.metadata
                    push!(r, "$k: $v")
                end
                push!(r, "Comment: $(rr.run.comment)")
            end
        end
        join(r, "\n")
    end
    label(value(runlabel), widget=builder["run"], signal=runlabel)

    poi_folder = joinpath(folder,"POI_videos")
    isdir(poi_folder) || mkdir(poi_folder)
    foreach(track, init=nothing) do _
        p = value(poi)
        i = first(find(p == x for x in a.pois))
        name = poi2videofile(folder, p, i)
        #=h5open(f5name, "w") do o
            @write o xyt
        end=#
        nothing
    end
    foreach(done,  init = nothing) do _
        destroy(builder["window"])
    end
    showall(builder["window"])

end

#=
using Images
import ImageView
import VideoIO
video_file = "/home/yakir/datasturgeon/projects/marie/afterLog/lana/POI_videos/5.mp4"
=##=io = VideoIO.open(video_file)
f = VideoIO.openvideo(io)=##=
# As a shortcut for just video, you can upen the file directly
# with openvideo
f = VideoIO.openvideo(video_file)
# One can seek to an arbitrary position in the video
seek(f,2.5)  ## The second parameter is the time in seconds and must be Float64
img = read(f)

canvas, _ = ImageView.view(img)

while !eof(f)
    read!(f, img)
    ImageView.view(canvas, img)
    #sleep(1/30)
end

=##=using GtkReactive, Gtk.ShortNames
c = canvas(UserUnit)
win = Window(c)
img = Signal(testimage("lighthouse"))
imgsig = map(img) do i
    i
end
redraw = draw(c, imgsig) do cnvs, image
    copy!(cnvs, image)
end
showall(win)=##=

import ImageView
import VideoIO
f = VideoIO.testvideo("annie_oakley")  # downloaded if not available
VideoIO.playvideo(f)  # no sound
=#


text:
const braille_signs = ['⠁' '⠂' '⠄' '⡀';
                       '⠈' '⠐' '⠠' '⢀']

"""
The type of canvas with the highest resolution
for Unicode-based plotting.
It uses the Unicode characters for
the Braille symbols to represent individual pixel.
This effectively turns every character into 8 pixels
that can individually be manipulated using binary operations.
"""
type BrailleCanvas <: Canvas
    grid::Array{Char,2}
    colors::Array{UInt8,2}
    pixel_width::Int
    pixel_height::Int
    origin_x::Float64
    origin_y::Float64
    width::Float64
    height::Float64
end

@inline pixel_width(c::BrailleCanvas) = c.pixel_width
@inline pixel_height(c::BrailleCanvas) = c.pixel_height
@inline origin_x(c::BrailleCanvas) = c.origin_x
@inline origin_y(c::BrailleCanvas) = c.origin_y
@inline width(c::BrailleCanvas) = c.width
@inline height(c::BrailleCanvas) = c.height
@inline nrows(c::BrailleCanvas) = size(c.grid, 2)
@inline ncols(c::BrailleCanvas) = size(c.grid, 1)

@inline x_pixel_per_char(::Type{BrailleCanvas}) = 2
@inline y_pixel_per_char(::Type{BrailleCanvas}) = 4

function BrailleCanvas(char_width::Int, char_height::Int;
                       origin_x::Real = 0.,
                       origin_y::Real = 0.,
                       width::Real  = 1.,
                       height::Real = 1.)
    char_width = max(char_width, 5)
    char_height = max(char_height, 2)
    pixel_width = char_width * x_pixel_per_char(BrailleCanvas)
    pixel_height = char_height * y_pixel_per_char(BrailleCanvas)
    width > 0 || throw(ArgumentError("Width has to be positive"))
    height > 0 || throw(ArgumentError("Height has to be positive"))
    grid = fill(Char(0x2800), char_width, char_height)
    colors = fill(0x00, char_width, char_height)
    BrailleCanvas(grid, colors,
                  pixel_width, pixel_height,
                  Float64(origin_x), Float64(origin_y),
                  Float64(width), Float64(height))
end

function pixel!(c::BrailleCanvas, pixel_x::Int, pixel_y::Int, color::Symbol)
    0 <= pixel_x <= c.pixel_width || return c
    0 <= pixel_y <= c.pixel_height || return c
    pixel_x = pixel_x < c.pixel_width ? pixel_x : pixel_x - 1
    pixel_y = pixel_y < c.pixel_height ? pixel_y : pixel_y - 1
    cw, ch = size(c.grid)
    tmp = pixel_x / c.pixel_width * cw
    char_x = floor(Int, tmp) + 1
    char_x_off = (pixel_x % 2) + 1
    if char_x < round(Int, tmp, RoundNearestTiesUp) + 1 && char_x_off == 1
        char_x = char_x + 1
    end
    char_y = floor(Int, pixel_y / c.pixel_height * ch) + 1
    char_y_off = (pixel_y % 4) + 1
    c.grid[char_x,char_y] = Char(UInt64(c.grid[char_x,char_y]) | UInt64(braille_signs[char_x_off, char_y_off]))
    c.colors[char_x,char_y] = c.colors[char_x,char_y] | color_encode[color]
    c
end

function printrow(io::IO, c::BrailleCanvas, row::Int)
    nunrows = nrows(c)
    0 < row <= nunrows || throw(ArgumentError("Argument row out of bounds: $row"))
    y = row
    for x in 1:ncols(c)
        print_color(c.colors[x,y], io, c.grid[x,y])
    end
end


text:
# This file includes code originally part of Julia.n
# Licensed under MIT License, see LICENSE.md

# The plan is to rewrite all of the functionality to not use the utf8proc library,
# and to use tables loaded up on initialization, as with StringLiterals.jl

# Currently, for characters outside the Latin1 subset, this depends on the following C calls:
# utf8proc_decompose
# utf8proc_reencode
# utf8proc_charwidth
# utf8proc_category
# utf8proc_category_string
# utf8proc_tolower
# utf8proc_toupper
# utf8proc_totitle

# For grapheme segmentation, this currently depends on the following 2 C calls:
# utf8proc_grapheme_break
# utf8proc_grapheme_break_stateful

# Unicode category constants
module Uni
const CN = 0
const LU = 1
const LL = 2
const LT = 3
const LM = 4
const LO = 5
const MN = 6
const MC = 7
const ME = 8
const ND = 9
const NL = 10
const NO = 11
const PC = 12
const PD = 13
const PS = 14
const PE = 15
const PI = 16
const PF = 17
const PO = 18
const SM = 19
const SC = 20
const SK = 21
const SO = 22
const ZS = 23
const ZL = 24
const ZP = 25
const CC = 26
const CF = 27
const CS = 28
const CO = 29

# strings corresponding to the category constants
const catstrings = [
    "Other, not assigned",
    "Letter, uppercase",
    "Letter, lowercase",
    "Letter, titlecase",
    "Letter, modifier",
    "Letter, other",
    "Mark, nonspacing",
    "Mark, spacing combining",
    "Mark, enclosing",
    "Number, decimal digit",
    "Number, letter",
    "Number, other",
    "Punctuation, connector",
    "Punctuation, dash",
    "Punctuation, open",
    "Punctuation, close",
    "Punctuation, initial quote",
    "Punctuation, final quote",
    "Punctuation, other",
    "Symbol, math",
    "Symbol, currency",
    "Symbol, modifier",
    "Symbol, other",
    "Separator, space",
    "Separator, line",
    "Separator, paragraph",
    "Other, control",
    "Other, format",
    "Other, surrogate",
    "Other, private use",
    "Invalid, too high",
    "Malformed, bad data",
]

const STABLE    = 1<<1
const COMPAT    = 1<<2
const COMPOSE   = 1<<3
const DECOMPOSE = 1<<4
const IGNORE    = 1<<5
const REJECTNA  = 1<<6
const NLF2LS    = 1<<7
const NLF2PS    = 1<<8
const NLF2LF    = NLF2LS | NLF2PS
const STRIPCC   = 1<<9
const CASEFOLD  = 1<<10
const CHARBOUND = 1<<11
const LUMP      = 1<<12
const STRIPMARK = 1<<13
end

############################################################################

utf8proc_error(result) =
    error(unsafe_string(ccall(:utf8proc_errmsg, Cstring, (Cssize_t,), result)))

function utf8proc_map(::Type{T}, str::MaybeSub{T}, options::Integer) where {T<:Str}
    nwords = ccall(:utf8proc_decompose, Int, (Ptr{UInt8}, Int, Ptr{UInt8}, Int, Cint),
                   str, sizeof(str), C_NULL, 0, options)
    nwords < 0 && utf8proc_error(nwords)
    buffer = Base.StringVector(nwords*4)
    nwords = ccall(:utf8proc_decompose, Int, (Ptr{UInt8}, Int, Ptr{UInt8}, Int, Cint),
                   str, sizeof(str), buffer, nwords, options)
    nwords < 0 && utf8proc_error(nwords)
    nbytes = ccall(:utf8proc_reencode, Int, (Ptr{UInt8}, Int, Cint), buffer, nwords, options)
    nbytes < 0 && utf8proc_error(nbytes)
    Str(cse(T), String(resize!(buffer, nbytes)))
end

utf8proc_charwidth(ch) = Int(ccall(:utf8proc_charwidth, Cint, (UInt32,), ch))

utf8proc_map(str::MaybeSub{T}, options::Integer) where {T<:Str} =
    utf8proc_map(T, convert(UTF8Str, str), options)
utf8proc_map(str::MaybeSub{<:Str{UTF8CSE}}, options::Integer) =
    utf8proc_map(UTF8Str, str, options)

@inline utf8proc_cat(ch) = ccall(:utf8proc_category, Cint, (UInt32,), ch)
@inline utf8proc_cat_abbr(ch) =
    unsafe_string(ccall(:utf8proc_category_string, Cstring, (UInt32,), ch))

@inline _lowercase_u(ch) = ccall(:utf8proc_tolower, UInt32, (UInt32,), ch)
@inline _uppercase_u(ch) = ccall(:utf8proc_toupper, UInt32, (UInt32,), ch)
@inline _titlecase_u(ch) = ccall(:utf8proc_totitle, UInt32, (UInt32,), ch)

############################################################################

function _normalize(::Type{T}, str::AbstractString;
                   stable::Bool      = false,
                   compat::Bool      = false,
                   compose::Bool     = true,
                   decompose::Bool   = false,
                   stripignore::Bool = false,
                   rejectna::Bool    = false,
                   newline2ls::Bool  = false,
                   newline2ps::Bool  = false,
                   newline2lf::Bool  = false,
                   stripcc::Bool     = false,
                   casefold::Bool    = false,
                   lump::Bool        = false,
                   stripmark::Bool   = false,
                   ) where {T<:Str}
    (compose & decompose) && unierror(UTF_ERR_DECOMPOSE_COMPOSE)
    (!(compat | stripmark) & (compat | stripmark)) && unierror(UTF_ERR_COMPAT_STRIPMARK)
    newline2ls + newline2ps + newline2lf > 1 && unierror(UTF_ERR_NL_CONVERSION)
    flags =
        ifelse(stable,      Uni.STABLE, 0) |
        ifelse(compat,      Uni.COMPAT, 0) |
        ifelse(decompose,   Uni.DECOMPOSE, 0) |
        ifelse(compose,     Uni.COMPOSE, 0) |
        ifelse(stripignore, Uni.IGNORE, 0) |
        ifelse(rejectna,    Uni.REJECTNA, 0) |
        ifelse(newline2ls,  Uni.NLF2LS, 0) |
        ifelse(newline2ps,  Uni.NLF2PS, 0) |
        ifelse(newline2lf,  Uni.NLF2LF, 0) |
        ifelse(stripcc,     Uni.STRIPCC, 0) |
        ifelse(casefold,    Uni.CASEFOLD, 0) |
        ifelse(lump,        Uni.LUMP, 0) |
        ifelse(stripmark,   Uni.STRIPMARK, 0)
    T(utf8proc_map(str, flags))
end

normalize(str::T, options::Integer) where {T<:Str} = _normalize(T, UTF8Str(str), options)
normalize(str::Str{UTF8CSE}, options::Integer) = _normalize(UTF8Str, str, options)

normalize(str::Str, nf::Symbol) =
    utf8proc_map(str, nf == :NFC ? (Uni.STABLE | Uni.COMPOSE) :
                 nf == :NFD ? (Uni.STABLE | Uni.DECOMPOSE) :
                 nf == :NFKC ? (Uni.STABLE | Uni.COMPOSE | Uni.COMPAT) :
                 nf == :NFKD ? (Uni.STABLE | Uni.DECOMPOSE | Uni.COMPAT) :
                 unierror(UTF_ERR_NORMALIZE, nf))

############################################################################

# iterators for grapheme segmentation

is_grapheme_break(c1::CodeUnitTypes, c2::CodeUnitTypes) =
    !(c1 <= 0x10ffff && c2 <= 0x10ffff) ||
    ccall(:utf8proc_grapheme_break, Bool, (UInt32, UInt32), c1, c2)

# Stateful grapheme break required by Unicode-9 rules: the string
# must be processed in sequence, with state initialized to Ref{Int32}(0).
# Requires utf8proc v2.0 or later.
is_grapheme_break!(state::Ref{Int32}, c1::CodeUnitTypes, c2::CodeUnitTypes) =
    ((c1 <= 0x10ffff && c2 <= 0x10ffff)
     ? ccall(:utf8proc_grapheme_break_stateful, Bool, (UInt32, UInt32, Ref{Int32}), c1, c2, state)
     : (state[] = 0; true))

is_grapheme_break(c1::Chr, c2::Chr) = is_grapheme_break(codepoint(c1), codepoint(c2))
is_grapheme_break(state::Ref{UInt32}, c1::Chr, c2::Chr) =
    is_grapheme_break(state, codepoint(c1), codepoint(c2))


text:
using Test, Statistics, LinearAlgebra, DSP, DSP.Util
using Plots, InteractiveViz
using SignalAnalysis
using SignalAnalysis.Units

@testset "signals" begin

  x = signal(randn(8000), 1000)
  @test x isa AbstractArray
  @test length(x) == 8000
  @test nframes(x) == 8000
  @test framerate(x) == 1000
  @test nchannels(x) == 1
  @test !isanalytic(x)
  @test domain(x) ≈ 0:1/1000:7999/1000

  x = signal(randn((8000,2)), 1000)
  @test x isa AbstractArray
  @test size(x) == (8000, 2)
  @test nframes(x) == 8000
  @test framerate(x) == 1000
  @test nchannels(x) == 2
  @test !isanalytic(x)
  @test domain(x) ≈ 0:1/1000:7999/1000

  x1 = signal(x, 1000)
  @test x === x1

  x1 = signal(x, 8000)
  @test framerate(x) == 1000
  @test framerate(x1) == 8000

  x = signal(randn(8000), 1000)
  @test x isa AbstractArray
  @test length(x) == 8000
  @test nframes(x) == 8000
  @test framerate(x) == 1000
  @test nchannels(x) == 1
  @test !isanalytic(x)
  @test domain(x) ≈ 0:1/1000:7999/1000

  x1 = signal(randn(8000), framerate(x))
  @test x1 isa AbstractArray
  @test nframes(x1) == 8000
  @test nchannels(x1) == 1
  @test framerate(x1) == 1000

  x1 = signal(randn(8000), framerate(x)/2)
  @test x1 isa AbstractArray
  @test nframes(x1) == 8000
  @test nchannels(x1) == 1
  @test framerate(x1) == 500

  x1 = signal(randn(8000), 1000)
  x = analytic(x1)
  @test nframes(x1) == nframes(x)
  @test nchannels(x1) == nchannels(x)
  @test framerate(x1) == framerate(x)
  @test isanalytic(x)
  @test √2*real(x) ≈ x1
  @test rms(x) ≈ rms(x1) atol=0.001
  x2 = analytic(x)
  @test isanalytic(x2)
  @test x === x2

  x1 = signal(randn((8000,2)), 1000)
  x = analytic(x1)
  @test nframes(x1) == nframes(x)
  @test nchannels(x1) == nchannels(x)
  @test framerate(x1) == framerate(x)
  @test isanalytic(x)
  @test √2*real(x) ≈ x1
  @test rms(x) ≈ rms(x1) atol=0.001
  x2 = analytic(x)
  @test isanalytic(x2)
  @test x === x2

  x = randn((8000,2))
  x1 = signal(x, 1000)
  @test x1 !== x
  @test samples(x1) === x
  @test samples(x) === x

  x = randn(8000)
  x1 = padded(x, (10, 5); delay=1)
  @test x1[-9] == 0
  @test x1[8005] == 0
  @test x1[1] == 0
  @test x1[2] == x[1]
  @test x1[8001] == x[8000]

  x1 = padded(x, (10, 5); delay=-2)
  @test x1[-9] == 0
  @test x1[8005] == 0
  @test x1[-1] == x[1]
  @test x1[7998] == x[8000]

  x = signal(collect(1:10), 1.0)
  @test length(collect(partition(x, 5))) == 2
  @test length(collect(partition(x, 4))) == 3
  @test length(collect(partition(x, 5; flush=false))) == 2
  @test length(collect(partition(x, 4; flush=false))) == 2
  @test length(collect(partition(x, 5; step=2))) == 5
  @test length(collect(partition(x, 4; step=2))) == 5
  @test length(collect(partition(x, 5; step=2, flush=false))) == 3
  @test length(collect(partition(x, 4; step=2, flush=false))) == 4

  x = signal(ones(1000), 8kHz)
  x2 = map(enumerate(partition(x, 250))) do (blknum, x1)
    @test size(x1) == (250,)
    sum(x1)*blknum
  end
  @test x2 == [250.0, 500.0, 750.0, 1000.0]
  x2 = map(enumerate(partition(x, 250))) do (blknum, x1)
    @test size(x1) == (250,)
    (blknum, sum(x1)*blknum)
  end
  @test x2 == [(1, 250.0), (2, 500.0), (3, 750.0), (4, 1000.0)]
  x2 = hcat(map(partition(x, 250)) do x1
    @test size(x1) == (250,)
    [sum(x1), prod(x1)]
  end...)'
  @test x2 == [250.0 1.0; 250.0 1.0; 250.0 1.0; 250.0 1.0;]
  for (blknum, x1) ∈ enumerate(partition(x, 250))
    x1 .= blknum
  end
  @test (x[1], x[251], x[501], x[751]) == (1.0, 2.0, 3.0, 4.0)

  x = signal(ones(1000, 2), 8kHz)
  x2 = map(enumerate(partition(x, 250))) do (blknum, x1)
    @test size(x1) == (250, 2)
    sum(x1)*blknum
  end
  @test x2 == [2*250.0, 2*500.0, 2*750.0, 2*1000.0]
  for x1 ∈ partition(x, 250)
    @test size(x1) == (250, 2)
  end

  x = signal(randn(2000), 8kHz)
  @test toframe(0.2s, x) == 1601
  @test toframe([0.2s, 201ms], x) == [1601, 1609]
  @test toframe(0.2:0.01:0.3, x)[1:3] == [1601, 1681, 1761]

  x = signal(randn(8000), 1000)
  t = toframe(0:0.1:1, x)
  @test t == 1:100:1001

end

@testset "generate" begin

  x = cw(5000, 0.2, 44100)
  @test nframes(x) == 8820
  @test nchannels(x) == 1
  @test framerate(x) == 44100
  @test isanalytic(x)
  @test abs.(x) ≈ ones(8820)
  @test mean(ifrequency(x)) ≈ 5000

  x1 = cw(-5000, 0.2, 44100)
  @test samples(x) .* samples(x1) ≈ ones(8820)
  @test x .* x1 ≈ ones(8820)

  x1 = cw(5kHz, 200ms, 44.1kHz)
  @test samples(x1) ≈ samples(x)
  @test framerate(x1) == framerate(x)

  x1 = cw(5kHz, 200ms, 44.1kHz; window=(tukey, 0.05))
  @test nframes(x1) == 8820
  @test nchannels(x1) == 1
  @test framerate(x1) == 44100
  @test isanalytic(x1)
  @test energy(x1) < energy(x)
  @test abs(x1[1]) < abs(x1[2]) < 1
  @test abs(x1[end]) < abs(x1[end-1]) < 1

  x2 = cw(5kHz, 200ms, 44.1kHz; window=hamming)
  @test nframes(x2) == 8820
  @test nchannels(x2) == 1
  @test framerate(x2) == 44100
  @test isanalytic(x2)
  @test energy(x2) < energy(x1)
  @test abs(x2[1]) < abs(x2[2]) < 1
  @test abs(x2[end]) < abs(x2[end-1]) < 1

  x1 = cw(5kHz, 200ms, 44.1kHz; phase=90°)
  @test nframes(x1) == 8820
  @test nchannels(x1) == 1
  @test framerate(x1) == 44100
  @test isanalytic(x1)
  @test energy(x1) ≈ energy(x)
  @test samples(x1) .* -1im ≈ samples(x)

  x = chirp(5kHz, 7kHz, 100ms, 44.1kHz)
  @test nframes(x) == 4410
  @test nchannels(x) == 1
  @test framerate(x) == 44100
  @test isanalytic(x)
  @test abs.(x) ≈ ones(4410)

  x3 = ifrequency(x)
  @test x3[2] ≈ 5000
  @test x3[end-1] ≈ 7000

  x1 = chirp(7kHz, 5kHz, 100ms, 44.1kHz)
  @test nframes(x) == 4410
  @test nchannels(x) == 1
  @test framerate(x) == 44100
  @test isanalytic(x)
  @test abs.(x) ≈ ones(4410)

  x1 = ifrequency(x1)
  @test x1[3] ≈ 7000 atol=10
  @test x1[end-2] ≈ 5000 atol=10

  x1 = chirp(5kHz, 7kHz, 100ms, 44.1kHz; shape=:hyperbolic, window=(tukey,0.05))
  @test nframes(x1) == 4410
  @test nchannels(x1) == 1
  @test framerate(x1) == 44100
  @test isanalytic(x1)
  @test energy(x1) < energy(x)
  @test abs(x1[1]) < abs(x1[2]) < 1
  @test abs(x1[end]) < abs(x1[end-1]) < 1

  x2 = ifrequency(x1)
  @test x2[3] ≈ 5000 atol=10
  @test x2[2205] < x3[2205]
  @test x2[end-2] ≈ 7000 atol=10

  x2 = chirp(5kHz, 7kHz, 100ms, 44.1kHz; phase=45°, window=hamming)
  @test nframes(x2) == 4410
  @test nchannels(x2) == 1
  @test framerate(x2) == 44100
  @test isanalytic(x2)
  @test energy(x2) < energy(x1)
  @test abs(x2[1]) < abs(x2[2]) < 1
  @test abs(x2[end]) < abs(x2[end-1]) < 1

  x2 = chirp(5kHz, 7kHz, 100ms, 44.1kHz; phase=90°)
  @test nframes(x2) == 4410
  @test nchannels(x2) == 1
  @test framerate(x2) == 44100
  @test isanalytic(x2)
  @test energy(x2) ≈ energy(x)
  @test samples(x2) .* -1im ≈ samples(x)

  x1 = chirp(7kHz, 5kHz, 100ms, 44.1kHz; shape=:hyperbolic)
  @test nframes(x1) == 4410
  @test nchannels(x1) == 1
  @test framerate(x1) == 44100
  @test isanalytic(x1)

  x2 = ifrequency(x1)
  @test x2[2] ≈ 7000
  @test x2[end-1] ≈ 5000

end

@testset "basic" begin

  x1 = signal(randn(8000), 1000)
  x2 = signal(randn((8000,2)), 1000)

  @test energy(x1) isa Real
  @test energy(x1) ≈ sum(abs2, x1)/1000
  @test energy(samples(x1)) isa Real
  @test energy(samples(x1)) ≈ 1000*energy(x1)
  @test energy(samples(x1); fs=1000) ≈ energy(x1)

  @test size(energy(x2)) == (2,)
  @test energy(x2) ≈ [energy(x2[:,1]), energy(x2[:,2])]
  @test energy(samples(x2)) ≈ 1000*energy(x2)
  @test energy(samples(x2); fs=1000) ≈ energy(x2)

  x1 = signal(vcat(ones(8000), zeros(4000)), 1000)
  x2 = signal(ones((8000,2)), 1000)

  @test meantime(x1) ≈ 4 atol=0.01
  @test meantime(samples(x1); fs=1000) ≈ 4 atol=0.01
  @test meantime(x2) ≈ [4, 4] atol=0.01
  @test meantime(samples(x2); fs=1000) ≈ [4, 4] atol=0.01

  @test rmsduration(x1) ≈ 2.3094 atol=0.01
  @test rmsduration(samples(x1); fs=1000) ≈ 2.3094 atol=0.01
  @test rmsduration(x2) ≈ [2.3094, 2.3094] atol=0.01
  @test rmsduration(samples(x2); fs=1000) ≈ [2.3094, 2.3094] atol=0.01

  x1 = cw(5000, 0.2, 44100)
  x2 = [x1 x1]
  @test mean(ifrequency(x1)) ≈ 5000
  @test meanfrequency(x1) ≈ 5000 atol=100
  @test dropdims(mean(ifrequency(x2); dims=1); dims=1) ≈ [5000, 5000]
  @test meanfrequency(x2) ≈ [5000, 5000] atol=100

  x1 = chirp(5kHz, 7kHz, 100ms, 44.1kHz)
  x2 = [x1 chirp(7kHz, 5kHz, 100ms, 44.1kHz)]
  @test meanfrequency(x1) ≈ 6000 atol=100
  @test meanfrequency(x2) ≈ [6000, 6000] atol=100*√2
  @test rmsbandwidth(x1) ≈ 577.35 atol=100
  @test rmsbandwidth(x2) ≈ [577.35, 577.35]  atol=100

end

@testset "dsp" begin

  x = 1 .+ randn(100000)

  f = fir(127, 0, 5kHz; fs=44.1kHz)
  @test length(f) == 127
  y = filt(f, 1, x)
  @test energy(y) < energy(x)*2/3
  @test mean(y) ≈ mean(x) atol=0.1

  f = fir(127, 15kHz; fs=44.1kHz)
  @test length(f) == 127
  y = filt(f, 1, x)
  @test energy(y) < energy(x)/2
  @test mean(y) ≈ 0 atol=0.01

  f = fir(127, 5kHz, 10kHz; fs=44.1kHz)
  @test length(f) == 127
  y = filt(f, 1, x)
  @test energy(y) < energy(x)/2
  @test mean(y) ≈ 0 atol=0.01

  @test mean(removedc(x)) ≈ 0 atol=0.01
  @test mean(x) ≈ 1 atol=0.01
  removedc!(x)
  @test mean(x) ≈ 0 atol=0.01

  x = randn(100000) .* (1 .+ sin.(2π*50 .* (1:100000)./44100))
  x = signal(x, 44.1kHz)
  y = demon(x)
  @test length(y) == 400
  @test meanfrequency(y) ≈ 50 atol=10
  y = demon(x; method=:mean)
  @test length(y) == 400
  @test meanfrequency(y) ≈ 50 atol=10
  y = demon(x; method=:fir)
  @test length(y) == 400
  @test meanfrequency(y) ≈ 50 atol=10
  y = demon(samples(x); fs=framerate(x))
  @test length(y) == 400
  @test meanfrequency(y) ≈ 50 atol=10

  x = randn(100000,2) .* (1 .+ sin.(2π*50 .* (1:100000)./44100))
  x = signal(x, 44.1kHz)
  y = demon(x)
  @test size(y) == (400, 2)
  @test meanfrequency(y) ≈ [50, 50] atol=10*√2
  y = demon(x; method=:mean)
  @test size(y) == (400, 2)
  @test meanfrequency(y) ≈ [50, 50] atol=10*√2
  y = demon(x; method=:fir)
  @test size(y) == (400, 2)
  @test meanfrequency(y) ≈ [50, 50] atol=10*√2
  y = demon(samples(x); fs=framerate(x))
  @test size(y) == (400, 2)
  @test meanfrequency(y) ≈ [50, 50] atol=10*√2

  x = signal(2*round.(Int,rand(1000)).-1 + 1im*(2*round.(Int,rand(1000)).-1), 100)/√2
  y = upconvert(x, 10, 100)
  @test framerate(y) == 1000
  @test rms(y[100:end-100])*sqrt(10) ≈ 1.0 atol=0.01
  @test meanfrequency(y) ≈ 100.0 atol=5.0
  z = downconvert(y, 10, 100)[12:end-11]
  @test framerate(z) == 100
  @test amp2db(rms(x-z)/rms(z)) < -40.0

  x = signal(2*round.(Int,rand(1000,2)).-1 + 1im*(2*round.(Int,rand(1000,2)).-1), 100)/√2
  y = upconvert(x, 10, 100)
  @test nchannels(y) == 2
  @test framerate(y) == 1000
  @test rms(y[100:end-100,1])*sqrt(10) ≈ 1.0 atol=0.01
  @test rms(y[100:end-100,2])*sqrt(10) ≈ 1.0 atol=0.01
  @test meanfrequency(y) ≈ [100.0, 100.0] atol=5.0
  z = downconvert(y, 10, 100)[12:end-11,:]
  @test nchannels(z) == 2
  @test framerate(z) == 100
  @test amp2db(rms(x-z)/rms(z)) < -40.0

  x = rcosfir(0.25, 10)
  @test length(x) == 221
  @test x[1:10:end] ≈ vcat(zeros(11), 0.326598866403003, zeros(11))
  @test sum(x.^2) ≈ 1.0

  x = rrcosfir(0.25, 10)
  @test length(x) == 221
  y = conv(x, x)[1:10:end]
  @test y ≈ vcat(zeros(22), 1.0, zeros(22)) atol=0.01
  @test sum(x.^2) ≈ 1.0

  for j ∈ 2:10
    x = mseq(j)
    @test x isa Array{Float64}
    @test length(x) == 2^j-1
    @test all(abs.(x) .== 1)
    y = circconv(x)
    @test all(y .== circconv(x, x))
    @test y[1] ≈ length(x)
    @test all(y[2:end] .≈ -1)
    x = gmseq(j)
    @test x isa Array{ComplexF64}
    @test length(x) == 2^j-1
    @test all(abs.(x) .== 1)
    y = circconv(x)
    @test y[1] ≈ length(x)
    @test rms(y[2:end]) ≈ 0 atol=1e-10
  end

  @test mseq((1,3)) == mseq(3)
  @test all(gmseq(3, 0) .== 1.0)

  x = cw(10kHz, 0.1s, 80kHz)
  @test abs(goertzel(x, 10kHz)) ≈ length(x) atol=1e-6
  @test abs(goertzel(x, 9kHz)) ≈ 0 atol=1e-6
  @test abs(goertzel(x, 11kHz)) ≈ 0 atol=1e-6
  @test abs.(goertzel(x, 10kHz, 512))[1:end-1] ≈ repeat([512.0], 15) atol=1e-6
  @test abs.(goertzel(x, 9375, 512))[1:end-1] ≈ repeat([0.0], 15) atol=1e-6
  x = hcat(cw(10kHz, 0.1s, 80kHz), cw(11kHz, 0.1s, 80kHz))
  @test abs.(goertzel(x, 10kHz)) ≈ [size(x,1), 0] atol=1e-6
  @test abs.(goertzel(x, 9kHz)) ≈ [0, 0] atol=1e-6
  @test abs.(goertzel(x, 11kHz)) ≈ [0, size(x,1)] atol=1e-6

  x = cw(1kHz, 1s, 80kHz)
  x2 = x .+ 0.25*rand(size(x)...) .+ 0.25im*rand(size(x)...)
  e2 = √mean(abs2.(x[800:end] .- x2[800:end]))
  x3 = pll(x2)
  @test meanfrequency(x3[800:end]) ≈ meanfrequency(x[800:end]) atol=0.1
  e3 = √mean(abs2.(x[800:end] .- x3[800:end]))
  @test e2/e3 > 3

  x = cw(7kHz, 1s, 44.1kHz)
  f = fir(127, 5kHz, 10kHz; fs=44.1kHz)
  x1 = sfilt(f, x)
  @test framerate(x1) == framerate(x)
  x1 = sfilt(f, 1, x)
  @test framerate(x1) == framerate(x)
  x1 = sfiltfilt(f, x)
  @test framerate(x1) == framerate(x)
  x1 = sresample(x, 3//2)
  @test framerate(x1) == 3 * framerate(x) / 2

  x = signal(randn(100), 10kHz)
  x1 = signal(vcat(zeros(1000), x/2, zeros(1000)), 10kHz)
  x2 = mfilter(x, x1)
  @test argmax(abs.(x2)) == 1001

  x = analytic(signal(randn(100), 10kHz))
  x1 = signal(vcat(zeros(1000), x/2, zeros(1000)), 10kHz)
  x2 = mfilter(x, x1)
  @test argmax(abs.(x2)) == 1001

  onesideds = [true, false]
  nfft = 256
  windows = [nothing, rect, tukey(nfft, 0.5), hanning, hamming]
  noverlaps = [0, 0, (15*nfft)÷16, (5*nfft)÷6, (5*nfft)÷6]
  for onesided in onesideds
    if onesided === true
      x = randn(96000)
    else
      x = randn(96000) + im .* randn(96000)
    end
    for (window, noverlap) in zip(windows, noverlaps)
      xstft = stft(x, nfft, noverlap; window=window, onesided=onesided)
      outputtype = isreal(x) ? Real : Complex
      x̂ = istft(outputtype, xstft; nfft=nfft, noverlap=noverlap, window=window)
      outputtype === Complex && (@test x̂ == istft(xstft; nfft=nfft, noverlap=noverlap, window=window))
      @test rms(x[nfft:length(x̂)-nfft]-x̂[nfft:end-nfft]) / rms(x[nfft:length(x̂)-nfft]) ≈ 0. atol=0.001
    end
  end

  fs = 9600
  hpf = fir(9, 1000; fs=fs)
  x = filtfilt(hpf, randn(96000))
  y = 0.1 .* real(chirp(500, 1000, 1.0, fs))
  x[9600+1:2*9600] .+= y
  nfft = 256
  noverlap = 0
  x̃ = whiten(x; nfft=nfft, noverlap=noverlap, window=nothing)
  ỹ = x̃[9600+1:2*9600]
  @test cor(y, ỹ) > 0.3
  x̃ = whiten(x; nfft=nfft, noverlap=noverlap, window=nothing, γ=0.)
  @test rms(x[1:length(x̃)] - x̃) / rms(x[1:length(x̃)]) ≈ 0. atol=0.0001

end

@testset "rand" begin

  x = rand(RedGaussian(100000))
  @test length(x) == 100000
  @test mean(x) ≈ 0 atol=1e-1
  @test std(x) ≈ 1 atol=1e-1
  x1 = 0
  x2 = 0
  for j ∈ 1:100
    x1 += abs2(sum(x .* cis.(2π*(100+j) .* (1:100000)./10000)))
    x2 += abs2(sum(x .* cis.(2π*(1000+j) .* (1:100000)./10000)))
  end
  @test 1 < log10(x1/x2) < 3

  x = rand(PinkGaussian(100000))
  @test length(x) == 100000
  @test mean(x) ≈ 0 atol=1e-1
  @test std(x) ≈ 1 atol=1e-1
  x1 = 0
  x2 = 0
  for j ∈ 1:100
    x1 += abs2(sum(x .* cis.(2π*(100+j) .* (1:100000)./10000)))
    x2 += abs2(sum(x .* cis.(2π*(1000+j) .* (1:100000)./10000)))
  end
  @test 0 < log10(x1/x2) < 2

  types = [Float16, Float32, Float64]
  for t ∈ types
    @test eltype(rand(RedGaussian(1000, one(t)))) == t
    @test eltype(rand(PinkGaussian(1000, one(t)))) == t
  end

end

@testset "plots" begin

  x = chirp(5kHz, 7kHz, 500ms, 44.1kHz)

  # plots may perhaps be tested with VisualRegressionTests.jl
  # for now we just test that the recipes don't die
  p = plot(x)
  @test p isa Plots.Plot
  p = plot(real(x))
  @test p isa Plots.Plot
  p = plot([x -x])
  @test p isa Plots.Plot
  p = plot(real([x -x]))
  @test p isa Plots.Plot
  p = plot(signal(randn(61040, 2), 1000))
  @test p isa Plots.Plot
  p = psd(x)
  @test p isa Plots.Plot
  p = psd(real(x))
  @test p isa Plots.Plot
  p = psd(real(x); xscale=:log10)
  @test p isa Plots.Plot
  p = psd([x -x])
  @test p isa Plots.Plot
  p = psd(real([x -x]))
  @test p isa Plots.Plot
  p = psd(real([x -x]); xscale=:log10)
  @test p isa Plots.Plot
  p = psd(samples(x); fs=framerate(x))
  @test p isa Plots.Plot
  p = psd(samples([x -x]); fs=framerate(x))
  @test p isa Plots.Plot
  p = specgram(x)
  @test p isa Plots.Plot
  p = specgram(x; pooling=nothing)
  @test p isa Plots.Plot
  p = specgram(samples(x); fs=framerate(x))
  @test p isa Plots.Plot
  f = fir(127, 15kHz; fs=44.1kHz)
  p = freqresp(f)
  @test p isa Plots.Plot
  p = freqresp(f; xscale=:log10)
  @test p isa Plots.Plot
  p = freqresp(f, [1])
  @test p isa Plots.Plot

  # interactive plots are hard to test
  # we just test that the calls don't crash

  iplot(x)
  @test true
  iplot(real(x))
  @test true
  iplot([x -x])
  @test true
  ispecgram(x)
  @test true
  ispecgram(real(x))
  @test true
  ispecgram(real(samples(x)))
  @test true
  ispecgram(real(samples(x)); fs=44100)
  @test true
  ispecgram(real(samples(x)); fs=44.1kHz)
  @test true

end

@testset "array" begin

  c = 1500.0
  θ = range(0.0, π; length=181)
  sd1 = steering(0.0:1.0:5.0, c, θ)
  @test size(sd1) == (6, 181)
  @test sd1[:,1] ≈ (2.5:-1.0:-2.5)/c atol=1e-6
  @test sd1[:,91] ≈ zeros(6) atol=1e-6
  @test sd1[:,181] ≈ (-2.5:1.0:2.5)/c atol=1e-6
  sd2 = steering(-2.0:1.0:3.0, c, θ)
  @test sd1 == sd2
  sd2 = steering(vcat((0.0:1.0:5.0)', zeros(1,6)), c, θ)
  @test sd1 == sd2
  sd2 = steering(vcat((0.0:1.0:5.0)', zeros(2,6)), c, θ)
  @test sd1 == sd2

  θ = range(0.0, 2π; length=9)
  rxpos = [-1.0 0.0 1.0 0.0; 0.0 1.0 0.0 -1.0]
  sd1 = steering(rxpos, c, θ)
  @test size(sd1) == (4, 9)
  @test sd1[:,1] ≈ [1.0, 0.0, -1.0, 0.0]/c atol=1e-6
  @test sd1[:,3] ≈ [0.0, -1.0, 0.0, 1.0]/c atol=1e-6
  @test sd1[:,5] ≈ -[1.0, 0.0, -1.0, 0.0]/c atol=1e-6
  @test sd1[:,7] ≈ -[0.0, -1.0, 0.0, 1.0]/c atol=1e-6
  @test sd1[1,4] ≈ sd1[2,4] atol=1e-6
  @test sd1[3,4] ≈ sd1[4,4] atol=1e-6
  @test sd1[3,4]-sd1[2,4] ≈ √2/c atol=1e-6

  θ = deg2rad.(1:360)
  fc = 500.0
  fs = 44100.0
  sd1 = steering(rxpos, c, θ)
  @test size(sd1) == (4, 360)
  x = cw(fc, 1.0, fs)
  x̄ = padded(x, 0; delay=round(Int, √2/c*fs))
  x4 = [x x x̄ x̄]
  bfo = beamform(x4, fc, 4096, sd1)
  @test argmax(bfo) == 135
  bfo = beamform(samples(x4), fc, 4096, sd1; fs=framerate(x4))
  @test argmax(bfo) == 135
  bfo = beamform(x4, fc, 4096, sd1; method=Bartlett())
  @test argmax(bfo) == 135
  bfo = beamform(x4 .+ 0.001*randn(size(x4)), fc, 4096, sd1; method=Capon())
  @test argmax(bfo) == 135
  bfo = beamform(x4, fc, 4096, sd1; method=Capon(0.1))
  @test argmax(bfo) == 135
  bfo = beamform(x4, fc, 4096, sd1; method=Music())
  @test argmax(bfo) == 135
  bfo = beamform(x4, fc, 4096, sd1; method=Music(1))
  @test argmax(bfo) == 135
  y4 = goertzel(x4, fc, 4096)
  bfo = beamform(y4, fc, sd1)
  @test argmax(bfo) == 135
  bfo = beamform(x4, sd1)
  @test nchannels(bfo) == 360
  e = energy(bfo)
  @test argmax(e) == 135
  @test maximum(e) ≈ 16.0 atol=0.1
  @test meanfrequency(bfo[:,135]) ≈ fc atol=10.0

end

@testset "tfa" begin

  x = real.(chirp(5kHz, 10kHz, 1s, 44.1kHz))
  y = tfd(x, Spectrogram())
  @test y isa SignalAnalysis.TFD
  @test time(y) ≈ 0.0029024943310657597:0.005804988662131519:0.9955555555555555
  @test freq(y) ≈ range(0.0, 22050.0; length=129)
  @test power(y) isa AbstractMatrix
  @test size(power(y)) == (129, 172)
  @test freq(y)[argmax(power(y)[:,86])] ≈ 7500.0 atol=100.0
  y = tfd(x, Spectrogram(nfft=512, noverlap=256, window=hamming))
  @test y isa SignalAnalysis.TFD
  @test time(y) ≈ 0.005804988662131519:0.005804988662131519:0.9926530612244898
  @test freq(y) ≈ range(0.0, 22050.0; length=257)
  @test power(y) isa AbstractMatrix
  @test size(power(y)) == (257, 171)
  @test freq(y)[argmax(power(y)[:,86])] ≈ 7500.0 atol=10.0

  x = real.(chirp(5kHz, 10kHz, 0.01s, 44.1kHz))
  y = tfd(x, Wigner())
  @test y isa SignalAnalysis.TFD
  @test time(y) ≈ range(0.0, 0.00998; length=441)
  @test freq(y) ≈ range(0.0, 22050.0, length=442)
  @test power(y) isa AbstractMatrix
  @test size(power(y)) == (442, 441)
  @test freq(y)[argmax(power(y)[:,220])] ≈ 7500.0 atol=1.0
  y = tfd(x, Wigner(nfft=512, smooth=10, method=:CM1980, window=hamming))
  @test y isa SignalAnalysis.TFD
  @test time(y) ≈ range(0.0, 0.00998; length=441)
  @test freq(y) ≈ range(0.0, 22050.0, length=257)
  @test power(y) isa AbstractMatrix
  @test size(power(y)) == (257, 441)
  @test freq(y)[argmax(power(y)[:,220])] ≈ 7500.0 atol=100.0

end


text:
"""
	WaylandCore

Provides basic Wayland functionality shared between client- and server-side. Connecting, passing messages, message-level types.
"""
module WaylandCore

export WlVersion, WlInt, WlUInt, WlFixed, WlString, WlID, WlNewID, WlObjID, WlArray, WlFD, WlMsgType, TypeofWlMsgType

import Sockets, Base.read, Base.write
using FixedPointNumbers

# Utility
"""
    typewrap(u)

Get a Type Union that matches all types in `u`.
"""
typewrap(u) = u isa Union ? Union{Type{u.a}, typewrap(u.b)} : Type{u}

const ByteBuffer = Vector{UInt8}

# Native Wayland types
# Opaque types
"""
	WlFD

A file descriptor, corresponds to "fd" wayland protocol type.
"""
const WlFD = RawFD

# Primitives:
"""
	WlUInt

Generic wayland message 32-bit word. Corresponds to "uint" wayland protocol type.
"""
const WlUInt = UInt32
"""
	WlInt

The "int" wayland message argument type.
"""
const WlInt = Int32
"""
	WlID

Either [`WlObjID`](@ref) or [`WlNewID`](@ref) type, as their differences are purely semantic.
"""
const WlID = WlUInt
"""
	WlObjID

Describes a [`WlID`](@ref) argument which represents an existing object. (The "object" wayland protocol type.)
"""
const WlObjID = WlID
"""
	WlNewID

Describes a [`WlID`](@ref) argument which will represent a new object. (The "new_id" wayland protocol type.)
"""
const WlNewID = WlID
"""
    WlFixed

24.8 signed fixed-point number based on Int32
1 sign bit, 23 bits integer prexision, 8 bits decimal precision
Corresponds to "fixed" wayland protocol type.

See also: [`Fixed`](@ref)
"""
const WlFixed = Fixed{Int32,8}

# Transparent types:
"""
	WlArray{T}

Direct representation of the "array" wayland protocol type.
"""
struct WlArray{T}
	size::WlUInt # size of (only, excluding the size) content in bytes
	content::Vector{T} # the content including padding to 32-bit boundary
end
"""
	WlByteArray

A [`WlArray`](@ref) containing generic bytes.
"""
const WlByteArray = WlArray{UInt8}
"""
	WlString

Direct representation of the "string" wayland protocol type.
"""
struct WlString
	length::WlUInt
	content::Vector{Cchar} # the content of the string including NUL terminator and padding to 32-bit boundary
end

# Meta types
"""
	WlVersion

The type used to denote interfaces' and their members' versions.
"""
const WlVersion = Int
"""
	WlMsgType

Represents any message argument type.
"""
const WlMsgType = Union{WlInt, WlUInt, WlFixed, WlString, WlNewID, WlObjID, WlArray, WlFD}
"""
    TypeofWlMsgType

The type that matches all (and only) types in WlMsgType
"""
const TypeofWlMsgType = typewrap(WlMsgType)
"""
	WaylandConnection

A connection that is able to receive all Wayland messages. Currently only a Unix domain socket.
"""
const WaylandConnection = Sockets.PipeEndpoint

# sendmsg- & recvmsg-based message passing:
const Csocklen_t = Cint
const default_outer_flags = 0
const default_inner_flags = 0
"""
	Ciovec

Corresponds to `struct iovec` from `uio.h`.
"""
struct Ciovec
	iov_base::ByteBuffer # Gets passed as a pointer to the buffer.
	iov_len::Csize_t # Byte length of `iov_base`.
end
"Constructor from a DenseVector{UInt8}, which knows its length"
Ciovec(vec::DenseVector{UInt8}) = Ciovec(vec, length(vec))
"""
	Cmsghdr

Corresponds to `struct msghdr` from `socket.h`.
"""
struct Cmsghdr
	msg_name::Ref{Cvoid} # The socket name if given socket is not connected. Never used, will be `C_NULL`
	msg_namelen::Csocklen_t # Length of `msg_name`. Never used, will be `0`
	msg_iov::DenseVector{Ciovec} # Vector of `Ciovec`s, which each represent a normal data buffer.
	msg_iovlen::Csize_t # Count of msg_iov elements.
	msg_control::ByteBuffer # Buffer for control messages
	msg_controllen::Csize_t # Byte length of msg_control
	msg_flags::Cint # Flags of the returned message
end
"Constructor without socket name and with optional flags."
function Cmsghdr(iov::DenseVector{Ciovec}, iovlen::Csize_t, control::DenseVector{UInt8}, controllen::Csize_t, flags::Cint = default_inner_flags)
	Cmsghdr(C_NULL, 0, iov, iovlen, control, controllen, flags)
end
"Constructor from `DenseVector`s, which know their length."
function Cmsghdr(iov::DenseVector{Ciovec}, control::DenseVector{UInt8}, flags::Cint = default_inner_flags)
	Cmsghdr(iov, length(iov), control, length(control), flags)
end

"""
    send(io::WaylandConnection, msg::Cmsghdr, flags::Cint = default_outer_flags)

The lowest level message sending mechanism, thin wrapper around `sendmsg` syscall.
"""
function send(io::WaylandConnection, msg::Cmsghdr, flags::Cint = default_outer_flags)
	ret = ccall(:sendmsg, Cssize_t,
	(RawFD, Ref{Cmsghdr}, Cint),
	fd(io), Ref(msg), flags)
	if ret == -1
		error("An error occured.") # TODO: what error?
	end
end
"""
    receive(io::WaylandConnection, msg::Cmsghdr, flags::Cint = default_outer_flags)

The lowest level message receiving mechanism, thin wrapper around `recvmsg` syscall.
"""
function receive(io::WaylandConnection, msg::Cmsghdr, flags::Cint = default_outer_flags)
	ret = ccall(:recvmsg, Cssize_t,
	(RawFD, Ref{Cmsghdr}, Cint),
	fd(io), Ref(msg), flags)
	if ret == 0
		# Server has shut down.
	elseif ret == -1
		error("An error occured.") # TODO: what error?
	end
	return ret
end

# Core type methods
"""
    read(io::IO, ::Type{WlByteArray})

Read a WlByteArray.
"""
function read(io::IO, ::Type{WlByteArray})
	size = read(io, UInt32)
	if size == 0
		return WlByteArray(size, [])
	end
	vec = Vector{UInt8}()
	for i in 1:size
		push!(vec, read(io, UInt8))
	end
	WlByteArray(size, vec)
end
"""
    read(io::IO, ::WlArray)

Read a WlArray of type T.
"""
function read(io::IO, ::Type{WlArray{T}}) where T
	size = read(io, UInt32)
	if size == 0
		WlArray(size, [])
	else
		vec = Vector{T}()
		tsize = sizeof(T)
		for i in tsize:tsize:size
			push!(vec, read(io, UInt8))
		end
		WlArray(size, vec)
	end
end
"""
    read(io::IO, ::Type{WlString})

Read a WlString.
"""
function read(io::IO, ::Type{WlString})
	length = read(io, WlUInt)
	contents = Array{Cchar}()
	for i in 1:length - 1 # length includes the terminating NUL, which we don't care about
		push!(contents, read(io, Cchar))
	end
	# Pull the terminating NUL and assert it's actually NUL
	@assert read(io, UInt8) == 0 "Last byte of string not actually NUL!"
	# Now pull the padding
	for i in 1:length % 4
		read(io, UInt8) # We don't care about padding contents, they are undefined anyway, we just want them out of the stream.
	end
end

# Library core types
"""
	WaylandMessage

Supertype for all messages.
"""
abstract type WaylandMessage end
"""
	BufferedMessage

Medium level, between high level messages and low level wrapped syscalls.
"""
struct BufferedMessage <: WaylandMessage
	from::UInt32
	opcode::UInt16
	size::UInt16
	payload::ByteBuffer
	ancillary::ByteBuffer
end
function receive(io::WaylandConnection, ::Type{BufferedMessage})
	body
end
"""
	LookupTable

An interface for looking up argument types. Contains a vector of message argument types in the order they should be received in. Anything implementing it should at least define getindex(x, id::UInt32, opcode::UInt16).
"""
abstract type LookupTable end
#= GenericMessage seems to not be compatible with passing FDs.
"""
	(::Type{<: WaylandMessage})(from::WlID, size::UInt16, opcode::UInt16)

Constructor with an empty payload for any message, which can accept `nothing` as payload.
"""
(type::Type{<: WaylandMessage})(from::WlID, size::UInt16, opcode::UInt16) = type(from, size, opcode, nothing)
"""
	GenericMessage

A generic, low-level message. Directly corresponds to the wayland wire format spec, using IOBuffer as storage. Inner constructor enforces padding to 32-bit word boundary.

**Do not** use for messages including file descriptors as they are passed with special calls on the socket and will be lost and introduce binary corruption if processed as just data.
"""
struct GenericMessage <: WaylandMessage
	from::WlID
	size::UInt16 # Overall (including the `from, size, opcode` header) message size in bytes. This means that a message without any payload has the size of 8. (4 bytes, 2 bytes, 2 bytes)
	opcode::UInt16 # Request/event opcode
	payload::Union{IOBuffer, Nothing} # Aligned to 32-bit words payload.
	GenericMessage(from::WlID, size::UInt16, opcode::UInt16, n::Nothing) = new(from, size, opcode, n)
	function GenericMessage(from::WlID, size::UInt16, opcode::UInt16, buf::IOBuffer)
		over_boundary = bytesavailable(buf) % 4
		if over_boundary == 0
			new(from, size, opcode, buf)
		elseif over_boundary == 1
			write(buf, 0x000000)
			new(from, size, opcode, buf)
		elseif over_boundary == 2
			write(buf, 0x0000)
			new(from, size, opcode, buf)
		elseif over_boundary == 3
			write(buf, 0x00)
			new(from, size, opcode, buf)
		end
	end
end
"""
	GenericMessage(from::WlID, size::UInt16, opcode::UInt16)

Constructor for messages with no payload.
"""
GenericMessage(from::WlID, size::UInt16, opcode::UInt16) = GenericMessage(from, size, opcode, nothing)
"""
	GenericMessage(from::WlID, size::UInt16, opcode::UInt16, iterable)

Constructor using an iterable collection as source of arguments.
"""
function GenericMessage(from::WlID, size::UInt16, opcode::UInt16, iterable)
	buf = IOBuffer()
	for value in iterable
		write(buf, value)
	end
	GenericMessage(from, size, opcode, buf)
end
"""
	GenericMessage(msg::WaylandMessage)

Constructor from any compliant WaylandMessage.
"""
function GenericMessage(msg::WaylandMessage)
	if WlObjID in msg.payload || WlNewID in msg.payload || WlID in msg.payload
		throw(ArgumentError("GenericMessage cannot be used with file descriptors!"))
	end
	GenericMessage(msg.from, msg.size, msg.opcode, msg.payload)
end=#
"""
	VectorMessage

A generic, low-level message. Uses Vector{WlMsgType} to store arguments.
"""
struct VectorMessage <: WaylandMessage
	from::WlID
	size::UInt16 # Overall (including the `from, size, opcode` header) message size in bytes. This means that a message without any payload has the size of 8. (4 bytes, 2 bytes, 2 bytes)
	opcode::UInt16 # Request/event opcode
	payload::Union{Vector{WlMsgType}, Nothing} # Stored arguments. Need to be aligned to 32-bit word boundary by write()
end
VectorMessage(from::WlID, size::UInt16, opcode::UInt16) = VectorMessage(from, size, opcode, nothing)

abstract type AbstractQueue end
"""
	GenericQueue

A message queue. This is a parametric queue, that can be specialised to only hold a stricter subset of messages for optimization.
"""
struct GenericQueue{T <: WaylandMessage} <: AbstractQueue
	queue::Vector{T}
end
"""
	MessageQueue

A generic message queue able to hold any messages.
"""
const MessageQueue = GenericQueue{WaylandMessage}

# Core library-side functions.
#= Commented out because they're currently incompatible with passing FDs.
# Message I/O methods:
"""
    write(io::IO, msg::WaylandMessage)

Write any `WaylandMessage` to the io. This requires a `GenericMessage(msg)` compatible constructor to exist.
"""
function write(io::IO, msg::WaylandMessage)
	write(io, GenericMessage(msg))
end
"""
	write(io::IO, msg::GenericMessage)

Write a `GenericMessage` to the io.
"""
function write(io::IO, msg::GenericMessage)
	if msg.payload == nothing
		write(io, msg.from, msg.size, msg.opcode)
	else
		write(io, msg.from, msg.size, msg.opcode, msg.payload)
	end
end
"""
    read(io::IO, ::Type{GenericMessage})

Read a `GenericMessage`. **Do not** use this for messages including file descriptors as they are passed with special calls on the socket and will not be read.
"""
function read(io::IO, ::Type{GenericMessage})
	from = read(io, WlID)
	size = read(io, UInt16)
	opcode = read(io, UInt16)
	if size == 8
		payload = nothing
	elseif size > 8
		payload = IOBuffer()
		remaining = size - 8
		while remaining > 0
			write(payload, read(io, UInt8))
			remaining -= 1
		end
	else
		error("Wrong size when reading a message!")
	end
	GenericMessage(from, size, opcode, payload)
end
"""
	read(io::IO, ::Type{VectorMessage}, lookup::Dict{UInt32, Dict{UInt16, AbstractVector{DataType}}})

Read a `VectorMessage` using a lookup table mapping object IDs to `opcode=>typevec` dictionaries.
"""
function read(io::IO, ::Type{VectorMessage}, lookup::LookupTable)
	from = read(io, UInt32)
	size = read(io, UInt16)
	@assert size % 4 == 0 "Wrong message size! Not multiple of 4."
	opcode = read(io, UInt16)
	pulled_bytes = 8
	payload = Vector{WlMsgType}()
	types = lookup[from, opcode]
	for type in types
		if size - pulled_bytes <= 0
			# If this gets executed, we're about to try reading something that exceeds message size.
			error("Wrong type vector, exceeded message size!")
		end
		arg = read(io, type)
		push!(payload, arg)
		if type isa WlString
			pulled_bytes += bytesize(arg)
		elseif type isa WlArray
			pulled_bytes += bytesize(arg)
		else
			pulled_bytes += sizeof(type)
		end
	end
	VectorMessage(from, size, opcode, payload)
end=#
# Connection and I/O
"""
    connect(name::AbstractString)

Connect to the named display. Returns a connection IO stream.
"""
function connect(path::AbstractString)
	path_isabsolute = path[1] == '/'
	runtimedir = get(ENV, "XDG_RUNTIME_DIR", nothing)
	if runtimedir == nothing && !path_isabsolute
		throw(ArgumentError("Cannot connect - XDG_RUNTIME_DIR is unset and non-absolute path was provided!"))
	elseif !path_isabsolute
		path = runtimedir * '/' * path
	end
	Sockets.connect(path)
end
"""
    connect()

Connect to default display, which is named in "WAYLAND_DISPLAY" env variable. If that is missing, connect to the hard default, "wayland-0".
"""
connect() = connect(get(ENV, "WAYLAND_DISPLAY", "wayland-0"))
"""
    disconnect(connection::IO)

Disconnect from a connection.

A high-level API should instead create and use a method for its Display object.
"""
function disconnect(connection::IO)
	Sockets.close(connection)
end

end # module


text:
#=
  German tank problem

  https://en.wikipedia.org/wiki/German_tank_problem

  This is a real (float) version of German tank problem

  See  ~/stan/german_tank.py
       ~/blog/german_tank_int.blog
       ~/webppl/german_tank.blog

  Cf german_tank_int.jl

  The theoretical answer is 382.5
  The advantage of this model using y as Uniform() is that we can use
  NN = 10000 and it's still fast (the full integer version is much slower).

  Interestingly, with NN = 10000 and N = randomInteger(NN), then
  enumerate calculates the value to exactly 383.0076434324275 (in 1.3s)
  which is almost exactly the same as the full integer version in german_tank_int2.wppl
  (which takes about 2min51s).

=#
# using Memoization
using Turing, StatsPlots, DataFrames
using ReverseDiff, Zygote, Tracker
# Turing.setadbackend(:reversediff)
# Turing.setadbackend(:zygote)
# Turing.setadbackend(:tracker)

include("jl_utils.jl")

# "Theoretical" Bayesian mean
# https://en.wikipedia.org/wiki/German_tank_problem
function theoretical(a)
    m = maximum(a)
    len = length(a)
    return (m-1)*(len-1)/(len-2)
end

# Note: Compared to  german_tank_in.jl both
# n and y are Uniform
@model function german_tank(y,nn=10_000)
    # n ~ DiscreteUniform(1,nn)
    n ~ Uniform(1,nn)
    # n ~ truncated(Normal(nn/2,nn/2),1,nn)
    for i in 1:length(y)
        y[i] ~ Uniform(1,n+1) # using Uniform instead of DiscreteUniform
    end

end

y = [10,256,202,97]
# y = [60] # Mosteller's version
nn = 10_000
println("theoretical: $(theoretical(y))")
model = german_tank(y,nn)

num_chains = 4

# n: 5005.3523   2889.8752
# chains = sample(model, Prior(), MCMCThreads(), 40_000, num_chains)

#
# chains = sample(model, MH(), 40_000)
# 382.1528 +/- 218.3875
chains = sample(model, MH(), MCMCThreads(), 40_000, num_chains)

# n:398.0311 +/- 321.7577  (it took 63s!)
# chains = sample(model, PG(20), MCMCThreads(), 40_000, num_chains)
# chains = sample(model, PG(20), 1_000)

# n: 4999.0273+/-2889.7176
# chains = sample(model, IS(), MCMCThreads(), 40_000, num_chains)
# chains = sample(model, IS(), 10_000)

# n: 956.9097 +/- 1271.8082
# chains = sample(model, SMC(), MCMCThreads(), 40_000, num_chains)
# chains = sample(model, SMC(), 10_000)

# n: 384.4288 +/- 203.1368
# chains = sample(model, NUTS(1000,0.65), 40_000)
# n: 384.5504 +/- 232.5935  (12.7s)
# chains = sample(model, NUTS(1000,0.65), MCMCThreads(), 40_000, num_chains)

# Error:
# LoadError: MethodError: no method matching phasepoint(::Random._GLOBAL_RNG, ::Array{Any,1},
# chains = sample(model, Gibbs(MH(:n),NUTS(1000,0.65,:y)), MCMCThreads(), 1_000, num_chains)
# chains = sample(model, Gibbs(MH(:n),NUTS(10,0.65,:y)), 1_000)

display(chains)

show_var_dist_pct(chains, :n,20)

max_stars_repo_name: Wikunia/hakank

max_issues_repo_name: Wikunia/hakank

max_forks_repo_name: Wikunia/hakank


text:
"""
```julia
    SoftMaxLikelihood()
```
Multiclass likelihood with [Softmax transformation](https://en.wikipedia.org/wiki/Softmax_function):

```math
p(y=i|{fₖ}) = exp(fᵢ)/ ∑ₖexp(fₖ)
```

There is no possible augmentation for this likelihood
"""
mutable struct SoftMaxLikelihood{T<:Real,A<:AbstractVector{T}} <: MultiClassLikelihood{T}
    nClasses::Int
    class_mapping::Vector{Any} # Classes labels mapping
    ind_mapping::Dict{Any,Int} # Mapping from label to index
    Y::Vector{BitVector} #Mapping from instances to classes
    y_class::Vector{Int64} #GP Index for each sample
    function SoftMaxLikelihood{T}(nClasses::Int) where {T<:Real}
        return new{T,Vector{T}}(nClasses)
    end
    function SoftMaxLikelihood{T}(
        nClasses::Int, labels::AbstractVector, ind_mapping::Dict
    ) where {T<:Real}
        return new{T,Vector{T}}(nClasses, labels, ind_mapping)
    end
    function SoftMaxLikelihood{T}(
        Y::AbstractVector{<:BitVector},
        class_mapping::AbstractVector,
        ind_mapping::Dict{<:Any,<:Int},
        y_class::AbstractVector{<:Int},
    ) where {T<:Real}
        return new{T,Vector{T}}(Y, class_mapping, ind_mapping, y_class)
    end
end

SoftMaxLikelihood(nClasses::Int) = SoftMaxLikelihood{Float64}(nClasses)
function SoftMaxLikelihood(ylabels::AbstractVector)
    return SoftMaxLikelihood{Float64}(
        length(ylabels), ylabels, Dict(value => key for (key, value) in enumerate(ylabels))
    )
end

implemented(::SoftMaxLikelihood, ::MCIntegrationVI) = true

function (::SoftMaxLikelihood)(f::AbstractVector)
    return StatsFuns.softmax(f)
end

function (l::SoftMaxLikelihood)(y::Int, f::AbstractVector{<:Real})
    return StatsFuns.softmax(f)[y]
end

function Base.show(io::IO, model::SoftMaxLikelihood{T}) where {T}
    return print(io, "Softmax likelihood")
end

function init_likelihood(
    likelihood::SoftMaxLikelihood{T},
    inference::AbstractInference{T},
    nLatent::Int,
    nSamplesUsed::Int,
) where {T}
    if inference isa GibbsSampling
        θ = [abs.(rand(T, nSamplesUsed)) * 2 for i in 1:nLatent]
        LogisticSoftMaxLikelihood{T}(
            likelihood.Y,
            likelihood.class_mapping,
            likelihood.ind_mapping,
            likelihood.y_class,
            θ,
        )
    else
        return likelihood
    end
end

function sample_local!(model::VGP{T,<:SoftMaxLikelihood,<:GibbsSampling}) where {T}
    model.likelihood.θ .= broadcast(
        (y::BitVector, γ::AbstractVector{<:Real}, μ::AbstractVector{<:Real}, i::Int64) ->
            rand.(PolyaGamma.(1.0, μ - logsumexp())),
        model.likelihood.Y,
        model.likelihood.γ,
        model.μ,
        1:(model.nLatent),
    )
    return nothing #TODO FINISH AT SOME POINT
end

function grad_samples(
    model::AbstractGP{T,<:SoftMaxLikelihood}, samples::AbstractMatrix{T}, index::Integer
) where {T}
    class = model.likelihood.y_class[index]
    grad_μ = zeros(T, nLatent(model))
    grad_Σ = zeros(T, nLatent(model))
    nSamples = size(samples, 1)
    samples .= mapslices(StatsFuns.softmax, samples; dims=2)
    t = 0.0
    @inbounds for i in 1:nSamples
        s = samples[i, class]
        @views g_μ = grad_softmax(samples[i, :], class) / s
        grad_μ += g_μ
        @views h = diaghessian_softmax(samples[i, :], class) / s
        grad_Σ += h - abs2.(g_μ)
    end
    for k in 1:nLatent(model)
        get_opt(inference(model), k).ν[index] = -grad_μ[k] / nSamples
        get_opt(inference(model), k).λ[index] = grad_Σ[k] / nSamples
    end
end

function log_like_samples(
    model::AbstractGP{T,<:SoftMaxLikelihood}, samples::AbstractMatrix, index::Integer
) where {T}
    class = model.likelihood.y_class[index]
    nSamples = size(samples, 1)
    return loglike = mapslices(logsumexp, samples; dims=2) / nSamples
end

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

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

function hessian_softmax(s::AbstractVector{<:Real}, i::Integer)
    m = length(s)
    hessian = zeros(m, m)
    for j in 1:m
        for k in 1:m
            hessian[j, k] =
                s[i] * ((δ(i, k) - s[k]) * (δ(i, j) - s[j]) - s[j] * (δ(j, k) - s[k]))
        end
    end
    return hessian
end


text:
using Polyopt
using Base.Test

approxzero{T<:Number}(p::Polyopt.Poly{T}, threshold=1e-6) = (norm(p.c, Inf) < threshold)

# example of unconstrained global optimization
let
    println("Unconstrained global optimization")

    x, z = variables(["x", "z"])
    f = 4.0*x^2 + x*z - 4*z^2 - 21//10*x^4 + 4*z^4 + 1//3*x^6

    # perturb problem to find global optimizer
    f = f + 1e-3*(x+z)
    prob = momentprob(3, f)

    X, Z, t, y, solsta = solve_mosek(prob)

    # test that (x, z) = (y[8], y[2]) is globally optimal
    xo = Polyopt.vectorize([x,z],6)*y
    @test abs(t - Polyopt.evalpoly(f, xo)) < 1e-6
end

# gloptipoly example
let
    println("Gloptipoly example")

    x1, x2, x3 = variables(["x1", "x2", "x3"])

    # This is the Gloptipoly example.
    f = -2*x1 + x2 - x3
    g = [ 24 - 20*x1 + 9*x2 - 13*x3 + 4*x1^2 - 4*x1*x2 + 4*x1*x3 + 2*x2^2 - 2*x2*x3 + 2*x3^2,
          4 - (x1 + x2 + x3),
          6 - (3*x2 + x3),
          x1, 2-x1,
          x2,
          x3, 3-x3 ]

    # perturb problem to find global optimizer
    f = f + 1e-3*(x1+x2+x3)

    prob = momentprob(4, f, g)
    X, Z, t, y, solsta = solve_mosek(prob)

    # test that (x1, x2, x3) extracted from y is globally optimal
    xo = Polyopt.vectorize([x1,x2,x3], 8)*y
    @test abs(t - Polyopt.evalpoly(f, xo)) < 1e-6
    @test all( [ Polyopt.evalpoly(gi, xo)  for gi=g ] .> -1e-4 )
end

# determine if f(x,z) is can be written as a SOS
let
    println("SOS example 1")

    x, z = variables(["x", "z"])
    f = 2*x^4 + 2*x^3*z - x^2*z^2 + 5*z^4

    prob = momentprob(2, f)
    X, Z, t, y, solsta = solve_mosek(prob)
    @test t > -1e-6

    v = monomials(2, [x,z])
    @test approxzero( f - dot(v,X[1]*v) )
end

# in this case f(x,z) is not SOS, but f(x,z)-t is SOS
let
    println("SOS example 2")

    x, z = variables(["x", "z"])

    f = (x + z + 1)^2 - 1
    prob = momentprob(1, f)

    X, Z, t, y, solsta = solve_mosek(prob)
    @test abs(t + 1) < 1e-6

    v = monomials(1, [x,z])

    # test that (f(x,z)-t) is SOS
    @test approxzero( f - t - dot(v,X[1]*v) )
end

# test duality for problem with inequality constraints
let
    println("Duality example 1")

    x, z = variables(["x", "z"])

    f = x^4 + z^2 + 1
    g = x^2 - z^2 - 1
    prob = momentprob(2, f, [g])

    X, Z, t, y, solsta = solve_mosek(prob)

    v1 = monomials(2,[x,z])
    v2 = monomials(1,[x,z])

    # test that f(x,z) - t - g(x,z)*s(x,z) is SOS,  where s(x,z) = v2'*X[2]*v2 is SOS
    @test approxzero( f - t - dot(v1,X[1]*v1) - g*dot(v2,X[2]*v2) )
end

# Test duality for problem with both inequality and equality constraints
let
    println("Duality example 2")

    x1, x2 = variables(["x1", "x2"])
    f = x1 + x2
    g = x1 - x2
    h = x1^2 + x2^2 - 1

    prob = momentprob(2, f, [g], [h])
    X, Z, t, y, solsta = solve_mosek(prob);

    u = monomials(2, [x1,x2])
    v = monomials(1, [x1,x2])

    # test optimality
    xo = Polyopt.vectorize([x1,x2], 4)*y
    @test abs(t - Polyopt.evalpoly(f, xo)) < 1e-4
    @test abs(Polyopt.evalpoly(h, xo)) < 1e-4
    @test Polyopt.evalpoly(g, xo) > -1e-4

    # test that f(x1,x2) - t - g(x1,x2)*s(x1,x2) - h(x1,x2)*w(x1,x2) is SOS,  where s1(x1,x2) = v'*X[2]*v is SOS, but w(x1,x2) = v'*Z[1]*z is not
    @test approxzero(f - t - (dot(u, X[1]*u) + g*dot(v, X[2]*v) + h*dot(v, Z[1]*v)))
end

# Example 6.23 "Sums of Squares, Moment Matrices and Polynomial Optimization", M. Laurent
let
    println("Duality example 3")

    x1, x2 = variables(["x1", "x2"])

    f = -x1 - x2
    g = [ 2*x1^4 - 8*x1^3 + 8*x1^2 + 2 - x2,
          4*x1^4 - 32*x1^3 + 88*x1^2 - 96*x1 + 36 - x2,
          x1, 3-x1,
          x2, 4-x2 ]

    prob = momentprob(4, f, g)

    X, Z, t, y, solsta = solve_mosek(prob)
    xo = Polyopt.vectorize([x1,x2], 8)*y
    @test abs(t - Polyopt.evalpoly(f, xo)) < 1e-6
    @test all( [ Polyopt.evalpoly(gi, xo) for gi=g ] .> -1e-4 )
end

# Example 6.25 "Sums of Squares, Moment Matrices and Polynomial Optimization", M. Laurent
let
    println("Constrained optimization, Laurent 1")

    x1, x2, x3 = variables(["x1", "x2", "x3"])

    f = x1^2*x2^2*(x1^2 + x2^2 - 3*x3^2) + x3^6 + 1e-2*(x1^6 + x2^6 + x3^6)
    g = 1 - (x1^2 + x2^2 + x3^2)

    prob = momentprob(4, f, [g])

    X, Z, t, y, solsta = solve_mosek(prob)
    xo = Polyopt.vectorize([x1,x2,x3], 8)*y
    @test abs(t - Polyopt.evalpoly(f, xo)) < 1e-6
    @test Polyopt.evalpoly(g, xo) > -1e-4
end

# Example 6.26 "Sums of Squares, Moment Matrices and Polynomial Optimization", M. Laurent
let
    println("Constrained optimization, Laurent 2")

    x1, x2, x3 = variables(["x1", "x2", "x3"])

    f = 1 + 1e-3*(x1+x2+x3)
    h = [ 5*x1^9 - 6*x1^5*x2 + x1*x2^4 + 2*x1*x3,
          -2*x1^6*x2 + 2*x1^2*x2^3 + 2*x2*x3,
          x1^2 + x2^2 - 0.265625 ]

    prob1 = momentprob(6, f, Polyopt.Poly{Int}[], h)

    X, Z, t, y, solsta = solve_mosek(prob1)
    xo = Polyopt.vectorize([x1,x2,x3],12)*y
    @test abs(t - Polyopt.evalpoly(f, xo)) < 1e-4
    @test norm([ Polyopt.evalpoly(hi, xo) for hi=h ], Inf) < 1e-4
end

# Example 8.8 "Sums of Squares, Moment Matrices and Polynomial Optimization", M. Laurent
let
    println("Constrained optimization, Laurent 3")

    x1, x2, x3 = variables(["x1", "x2", "x3"]);

    g1, g2 = x1^4 + (x1*x2-1)^2, x2^2*x3^2 + (x3^2-1)^2;
    f = g1+g2;
    g = [g1, g2, 1-x1^2, 1-x2^2, 1-x3^2];

    prob = momentprob(3, f, g);
    X, Z, t, y, solsta = solve_mosek(prob);

    probc = momentprob_chordalembedding(3, f, g);
    Xc, Zc, tc, yc, solstac = solve_mosek(probc);

    @test abs(t-tc) < 1e-6
end

# http://gamsworld.org/global/globallib/ex2_1_1.htm
let
    println("Globallib ex2_1_1")

    x1, x2, x3, x4, x5 = variables(["x1", "x2", "x3", "x4", "x5"])

    f =  -0.5*(100*x1^2 + 100*x2^2 + 100*x3^2 + 100*x4^2 + 100*x5^2) + 42*x1 + 44*x2 + 45*x3 + 47*x4 + 47.5*x5
    g = [40 - (20*x1 + 12*x2 + 11*x3 + 7*x4 + 4*x5), x1, 1-x1, x2, 1-x2, x3, 1-x3, x4, 1-x4, x5, 1-x5]

    prob = momentprob(3, f, g)
    X, Z, t, y, solsta = solve_mosek(prob);
    xo = Polyopt.vectorize([x1,x2,x3,x4,x5],6)*y
    @test all( [ Polyopt.evalpoly(gi, xo) for gi=g ] .> -1e-4 )
end

# http://gamsworld.org/global/globallib/ex2_1_2.htm
let
    println("Globallib ex2_1_2")

    x1, x2, x3, x4, x5, x6 = variables(["x1", "x2", "x3", "x4", "x5", "x6"])

    f = -0.5*(x1*x1 + x2*x2 + x3*x3 + x4*x4 + x5*x5) - 10.5*x1 - 7.5*x2 - 3.5*x3 - 2.5*x4 - 1.5*x5 - 10*x6
    g = [ 6.5 - (6*x1 + 3*x2 + 3*x3 + 2*x4 + x5), 20 - (10*x1 + 10*x3 + x6), x1, 1-x1, x2, 1-x2, x3, 1-x3, x4, 1-x4, x5, 1-x5, x6 ]

    prob = momentprob(2, f, g)
    X, Z, t, y, solsta = solve_mosek(prob);
    xo = Polyopt.vectorize([x1,x2,x3,x4,x5,x6],4)*y
    @test all( [ Polyopt.evalpoly(gi, xo) for gi=g ] .> -1e-4 )
end

# http://gamsworld.org/global/globallib/ex2_1_4.htm
let
    println("Globallib ex2_1_4")

    x1, x2, x3, x4, x5, x6 = variables(["x1", "x2", "x3", "x4", "x5", "x6"])

    f = (6.5*x1 - 0.5*x1*x1) - x2 - 2*x3 - 3*x4 - 2*x5 - x6
    g = [ 16-(x1 + 2*x2 + 8*x3 + x4 + 3*x5 + 5*x6),
          -1-(-8*x1 - 4*x2 - 2*x3 + 2*x4 + 4*x5 - x6),
          24-(2*x1 + 0.5*x2 + 0.2*x3 - 3*x4 - x5 - 4*x6),
          12-(0.2*x1 + 2*x2 + 0.1*x3 - 4*x4 + 2*x5 + 2*x6),
           3-(-0.1*x1 - 0.5*x2 + 2*x3 + 5*x4 - 5*x5 + 3*x6),
           x1, 1-x1, x2, x3, x4, 1-x4, x5, 1-x5, x6, 2-x6 ]

    prob = momentprob(2, f, g)
    X, Z, t, y, solsta = solve_mosek(prob);
    xo = Polyopt.vectorize([x1,x2,x3,x4,x5,x6],4)*y
    @test all( [ Polyopt.evalpoly(gi, xo) for gi=g ] .> -1e-4 )
end

# http://gamsworld.org/global/globallib/ex3_1_2.htm
let
    println("Globallib ex3_1_2")

    x1, x2, x3, x4, x5 = variables(["x1", "x2", "x3", "x4", "x5"])

    f = -40792.141 + 0.8356891*x1*x5 + 37.293239*x1 + 5.3578547*x3*x3
    g = [ 6.665593-(0.0056858*x2*x5 - 0.0022053*x3*x5 + 0.0006262*x1*x4),
          85.334407-(0.0022053*x3*x5 - 0.0056858*x2*x5 - 0.0006262*x1*x4),
          29.48751-(0.0071317*x2*x5 + 0.0021813*x3*x3 + 0.0029955*x1*x2),
          -9.48751-(-0.0071317*x2*x5 - 0.0021813*x3*x3 - 0.0029955*x1*x2),
          15.599039-(0.0047026*x3*x5 + 0.0019085*x3*x4 + 0.0012547*x1*x3),
          -10.699039-((-0.0047026*x3*x5) - 0.0019085*x3*x4 - 0.0012547*x1*x3),
          x1-78, 102-x1, x2-33, 45-x2, x3-27, 45-x3, x4-27, 45-x4, x5-27, 45-x5 ]

    prob = momentprob(2, f, g)
    X, Z, t, y, solsta = solve_mosek(prob);
    xo = Polyopt.vectorize([x1,x2,x3,x4,x5],4)*y
    @test all( [ Polyopt.evalpoly(gi, xo) for gi=g ] .> -1e-4 )
end

# http://gamsworld.org/global/globallib/ex3_1_1.htm
let
    println("Globallib ex3_1_1")

    x1, x2, x3, x4, x5, x6, x7, x8 = variables(["x1", "x2", "x3", "x4", "x5", "x6", "x7", "x8"])
    #f = x1 + x2 + x3
    #g = [ 1 - (0.0025*x4 + 0.0025*x6),
    #      1 - (-0.0025*x4 + 0.0025*x5 + 0.0025*x7),
    #      1 - (-0.01*x5 + 0.01*x8),
    #      83333.333 - (100*x1 - x1*x6 + 833.33252*x4),
    #      -(x2*x4 - x2*x7 - 1250*x4 + 1250*x5),
    #      -1250000 - (x3*x5 - x3*x8 - 2500*x5),
    #      x1-100,  10000-x1,
    #      x2-1000, 10000-x2,
    #      x3-1000, 10000-x3,
    #      x4-10,   1000-x4,
    #      x5-10,   1000-x5,
    #      x6-10,   1000-x6,
    #      x7-10,   1000-x7,
    #      x8-10,   1000-x8 ]

    # rescale model to make it easier to solve
    f = x1 + x2 + x3
    g = [ 1 - (2.5*x4 + 2.5*x6),
          1 - (-2.5*x4 + 2.5*x5 + 2.5*x7),
          1 - (-10*x5 + 10*x8),
          0.08333 - (x1 - 10*x1*x6 + 0.83333*x4),
          -(x2*x4 - x2*x7 - 0.1250*x4 + 0.125*x5),
          -0.125 - (x3*x5 - x3*x8 - 0.25*x5),
          x1-0.01, 1-x1,
          x2-0.1,  1-x2,
          x3-0.1,  1-x3,
          x4-0.01, 1-x4,
          x5-0.01, 1-x5,
          x6-0.01, 1-x6,
          x7-0.01, 1-x7,
          x8-0.01, 1-x8 ]

    probc = momentprob_chordalembedding(3, f, g)
    Xc, Zc, tc, yc, solstac = solve_mosek(probc);
    xo = Polyopt.vectorize([x1,x2,x3,x4,x5,x6,x7,x8],6)*yc
    @test all( [ Polyopt.evalpoly(gi, xo) for gi=g ] .> -1e-4 )
end

#http://gamsworld.org/global/globallib/ex5_4_2.htm
let
    println("Globallib ex5_4_2")

    x1,x2,x3,x4,x5,x6,x7,x8 = variables(["x1","x2","x3","x4","x5","x6","x7","x8"])

    #f = x1 + x2 + x3
    #g = [ 400 - (x4 + x6),
    #      300 - (-x4 + x5 + x7),
    #      100 - (-x5 + x8),
    #      83333.3333333333 - (x1 - x1*x6 + 833.333333333333*x4),
    #      -(x2*x4 - x2*x7 - 1250*x4 + 1250*x5),
    #      -1250000 - (x3*x5 - x3*x8 - 2500*x5),
    #      x1-100, 10000-x1,
    #      x2-1000, 10000-x2,
    #      x3-1000, 10000-x3,
    #      x4-10, 1000-x4,
    #      x5-10, 1000-x5,
    #      x6-10, 1000-x6,
    #      x7-10, 1000-x7,
    #      x8-10, 1000-x8 ]

    # rescale model to make it easier to solve
    f = x1 + x2 + x3
    g = [ 400/1000 - (x4 + x6),
          300/1000 - (-x4 + x5 + x7),
          100/1000 - (-x5 + x8),
          .083333 - (0.01*x1 - 10*x1*x6 + 0.83333*x4),
          -(x2*x4 - x2*x7 - 0.125*x4 + 0.125*x5),
          -0.125 - (x3*x5 - x3*x8 - 0.25*x5),
          x1-100/10000, 1-x1,
          x2-1000/10000, 1-x2,
          x3-1000/10000, 1-x3,
          x4-10/1000, 1-x4,
          x5-10/1000, 1-x5,
          x6-10/1000, 1-x6,
          x7-10/1000, 1-x7,
          x8-10/1000, 1-x8 ]

    probc = momentprob_chordalembedding(3, f, g)
    Xc, Zc, tc, yc, solstac = solve_mosek(probc);
    xo = Polyopt.vectorize([x1,x2,x3,x4,x5,x6,x7,x8],6)*yc
    @test all( [ Polyopt.evalpoly(gi, xo) for gi=g ] .> -1e-4 )
end

let
    println("Chordal relaxation")

    x1,x2,x3 = variables(["x1","x2","x3"])

    f = x1-x2+x3
    g = [1-x1*x2,
         1-x2*x3,
        1-x1^2,
        1-x2^2,
        1-x3^2 ]

    order = 2
    prob = momentprob(order, f, g)
    X, Z, t, y, solsta = solve_mosek(prob);
    xo = Polyopt.vectorize([x1,x2,x3],2*order)*y
    @test all( [ Polyopt.evalpoly(gi, xo) for gi=g ] .> -1e-4 )

    K = Array{Int,1}[ [1,2], [2,3] ]
    probc = momentprob_chordal(order, K, f, g)
    Xc, Zc, tc, yc, solstac = solve_mosek(probc);
    xo = Polyopt.vectorize([x1,x2,x3],2*order)*yc
    @test all( [ Polyopt.evalpoly(gi, xo) for gi=g ] .> -1e-4 )

    u1 = monomials(order, [x1,x2,x3][K[1]])
    u2 = monomials(order, [x1,x2,x3][K[2]])
    v1 = monomials(order-1, [x1,x2,x3][K[1]])
    v2 = monomials(order-1, [x1,x2,x3][K[2]])

    r = dot(u1, Xc[1]*u1) + dot(u2, Xc[2]*u2) +
            g[1]*dot(v1, Xc[3]*v1) +
            g[2]*dot(v2, Xc[4]*v2) +
            g[3]*dot(v1, Xc[5]*v1) +
            g[4]*(dot(v1, Xc[6]*v1) + dot(v2, Xc[7]*v2)) +
            g[5]*dot(v2, Xc[8]*v2)

    @test approxzero( f - t - r )
end


text:
using Slacker
using Test
using HTTP
using JSON

import Slacker.getConfigFile

function backup()
    src = getConfigFile()
    target = joinpath(tempdir(),"slackerbackup.json" )
    if isfile(src)
        cp(src, target, force = true)
    end
end

function restore()
    src = getConfigFile()
    target = joinpath(tempdir(),"slackerbackup.json" )
    if isfile(target)
        mv(target, src, force = true)
    end
end

@testset "SlackConfig Constructors" begin

    # does blank object instantiation work?
    @test isa(SlackConfig(), SlackConfig)

end

@testset "Creating and writing Config Files" begin

    # ensure existing configuration is not overridden
    backup()
    addConfig(SlackConfig(), "testConfig")
    @test isfile(getConfigFile())

    conf = loadConfig("testConfig")

    @test conf.user == "JuliaBot"
    @test conf.channel == "#general"
    @test conf.icon_emoji == ":ghost:"

    addConfig(SlackConfig("webhook", "name", "channel", "icon"), "testConfig2")
    conf2 = loadConfig("testConfig2")

    @test conf2.webhook_url == "webhook"
    @test conf2.user == "name"
    @test conf2.channel == "channel"
    @test conf2.icon_emoji == "icon"

    #@test_warn "Overwriting exisiting configuration." addConfig(SlackConfig("webhook", "name2", "channel2", "icon2"), "testConfig2")
    conf2 = loadConfig("testConfig2")

    removeConfigFile()
    @test !isfile(getConfigFile())
    restore()

end

@testset "Removing Configs" begin
    addConfig("DUMMY_URL", "testConfig3")
    @test hasConfig("testConfig3") == true
    removeConfig("testConfig3")
    @test hasConfig("testConfig3") == false

end

@testset "Sending Messages" begin
    message = "This is a complicated Test Message. :smile:"
    addConfig("https://postman-echo.com/post", "echo_config")
    res = sendMessage(message, "echo_config")
    @test res.status==200
    @test res.headers[1][1]=="Content-Type"
    @test res.headers[1][2]=="application/json; charset=utf-8"
    body = String(res.body)
    dic = JSON.parse(body)
    content = JSON.parse(dic["json"]["payload"])
    @test content["channel"] == "#random"
    @test content["username"] == "Juliabot"
    @test content["icon_emoji"] == ":ghost:"
    @test content["text"] == message
    removeConfig("echo_config")
    @test hasConfig("echo_config") == false

end

#TODO Test overwriting of file?
# Loadconfig not found error
# sending... figure out how :D