text:
using Printf
using Statistics
using LaTeXStrings

"""
    UncertainValue

Represents a floating point numerical value and an associated uncertainty (1σ).
"""
struct UncertainValue <: AbstractFloat
    value::Float64
    sigma::Float64

    function UncertainValue(val::Real, sigma::Real)
        @assert sigma >= 0.0
        return new(convert(Float64, val), convert(Float64, sigma))
    end
end

Base.convert(::Type{UncertainValue}, n::Real) = UncertainValue(convert(Float64, n), 0.0)
Base.convert(::Type{UncertainValue}, n::UncertainValue) = n

Base.zero(::Type{UncertainValue}) = UncertainValue(0.0, 0.0)
Base.one(::Type{UncertainValue}) = UncertainValue(1.0, 0.0)
Base.add_sum(x::UncertainValue, y::UncertainValue)::Real = x.value + y.value
Base.add_sum(x::AbstractFloat, y::UncertainValue)::Real = x + y.value
Base.hash(uv::UncertainValue, h::UInt) = Base.hash(uv.value, Base.hash(uv.sigma, h))

"""
    Base.round(uva::UncertainValue; defdigits=4)

Provides a (somewhat) intelligent formatter for `UncertainValue` objects that uses the
uncertainty to suggest the number of digits to display.
"""
function Base.round(uva::UncertainValue; defdigits = 4)
    r =
        (σ(uva) >= 1.0e-10 * value(uva)) && (abs(value(uva)) > 1.0e-100) ?
        convert(Int, ceil(max(1, log10(max(σ(uva), abs(value(uva))) / σ(uva))))) + 2 :
        defdigits
    return uv(
        round(value(uva), RoundNearestTiesUp, sigdigits = r),
        round(σ(uva), RoundNearestTiesUp, sigdigits = 2),
    )
end

"""
    poisson(val::Int)

Creates an UncertainValue from an integer which is assumed to have σ = sqrt(val).
"""
poisson(val::Int) = UncertainValue(convert(Float64, val), sqrt(convert(Float64, val)))

"""
    uv(val::Real, σ::Real)

Create an UncertainValue from a real value and 1σ uncertainty.
"""
uv(val::Real, σ::Real=0.0) = UncertainValue(val, σ)
uv(uuvv::UncertainValue) = uuvv

"""
    max(uv1::UncertainValue, uv2::UncertainValue)

Determines the maximum of two values by comparing first the value and second the
uncertainty. (Equal value defer to the larger uncertainty being considered the 'max'.)

If `value(uv1)==value(uv2)` then it is undefined which is larger but to ensure a constant
ordering I define the one with a larger uncertainty to be larger since it is more probable
that it could be larger.
"""
Base.max(uv1::UncertainValue, uv2::UncertainValue) =
    uv1.value ≠ uv2.value ? (uv1.value > uv2.value ? uv1 : uv2) :
    (uv1.sigma > uv2.sigma ? uv1 : uv2)

"""
    min(uv1::UncertainValue, uv2::UncertainValue)

Determines the maximum of two values by comparing first the value and second the
uncertainty. (Equal value defer to the larger uncertainty being considered the 'min'.)

If `value(uv1)==value(uv2)` then it is undefined which is smaller but to ensure a constant
ordering I define the one with a larger uncertainty to be smaller since it is more probable
that it could be smaller.
"""
Base.min(uv1::UncertainValue, uv2::UncertainValue) =
    uv1.value ≠ uv2.value ? (uv1.value > uv2.value ? uv2 : uv1) :
    (uv1.sigma > uv2.sigma ? uv1 : uv2)

Base.minimum(uvs::AbstractVector{UncertainValue}) = reduce(min, uvs)

Base.minimum(xs::Tuple{UncertainValue,Vararg{UncertainValue}}) = reduce(min, xs)

Base.maximum(uvs::AbstractVector{UncertainValue}) = reduce(max, uvs)

Base.maximum(xs::Tuple{UncertainValue,Vararg{UncertainValue}}) = reduce(max, xs)

Base.sum(xs::UncertainValue...) =
    UncertainValue(sum(value(x) for x in xs), sqrt(sum(variance(x) for x in xs)))

Base.sum(uvs::AbstractVector{UncertainValue}) =
    UncertainValue(sum(value(x) for x in uvs), sqrt(sum(variance(x) for x in uvs)))

"""
    Statistics.mean(uvs::AbstractVector{UncertainValue}, minvar=0.0e-10)

The variance weighted mean of a collection of UncertainValue items.
(see https://en.wikipedia.org/wiki/Weighted_arithmetic_mean#Variance_weights).
Be careful! Adding a single value with value!=0 and σ=0.0 will cause NaN.
"""
function Statistics.mean(uvs::AbstractVector{UncertainValue})
    return UncertainValue(
        sum(value(x) / variance(x) for x in uvs) / sum(1.0 / variance(x) for x in uvs),
        sqrt(1.0 / sum(1.0 / variance(x) for x in uvs)),
    )
end
function Statistics.mean(uvs::UncertainValue...)
    return UncertainValue(
        sum(value(x) / variance(x) for x in uvs) / sum(1.0 / variance(x) for x in uvs),
        sqrt(1.0 / sum(1.0 / variance(x) for x in uvs)),
    )
end

"""
    Statistics.std(uvs::UncertainValue...)

The standard deviation of the value portion of the collection of UncertainValue items.
"""
Statistics.std(uvs::UncertainValue...) = std(value(x) for x in uvs)
Statistics.std(uvs::AbstractVector{UncertainValue}) = std(value(x) for x in uvs)

"""
    variance(uv::UncertainValue)

Returns σ².
"""
variance(uv::UncertainValue) = uv.sigma^2

"""
    variance(f::Real)

Returns 0.
"""
variance(f::Real) = zero(typeof(f))

Base.isless(uv1::UncertainValue, uv2::UncertainValue) =
    isequal(uv1.value, uv2.value) ? !isless(uv1.sigma, uv2.sigma) :
    isless(uv1.value, uv2.value)

"""
    σ(uv::UncertainValue)

Returns the 1-σ uncertainty)
"""
σ(uv::UncertainValue) = uv.sigma

"""
    σ(r::Real)

Returns 0.
"""
σ(::Real) = zero(Float64)
σ(m::Missing) = m

"""
    uncertainty(uv::UncertainValue, k::Real=1.0)

Returns the k-σ uncertainty (defaults to k=1.0)
"""
uncertainty(uv::UncertainValue, k::Float64 = 1.0) = k * uv.sigma

"""
    uncertainty(f::Real, k::Real=1.0)

Returns 0.0.
"""
uncertainty(::Real, ::Real) = zero(Float64)

"""
    fractional(uv::UncertainValue)

Computes the fractional uncertainty.
"""
fractional(uv::UncertainValue)::Float64 =
    abs(uv.value) > 1.0e-15 * abs(uv.sigma) ? abs(uv.sigma / uv.value) : Inf64

"""
    fractional(f::Real)

Returns 0
"""
fractional(f::Real) = zero(typeof(f))

"""
    value(uv::UncertainValue)

Returns the value portion. (uv.value)
"""
value(uv::UncertainValue) = uv.value

"""
    value(f::Real)

Returns f
"""
value(f::Real) = Float64(f)
value(m::Missing) = m

"""
    pearson(uv1::UncertainValue, uv2::UncertainValue, covar::Real)

Computes the Pearson correlation coefficient given the covariance between two UncertainValue.
"""
function pearson(uv1::UncertainValue, uv2::UncertainValue, covar::Real)
    res = covar / (σ(uv1) * σ(uv2))
    @assert (res >= -1.0) & (res <= 1.0) "The Pearson correlation coefficient must be on [-1.0, 1.0] - $res"
    return res
end

"""
    covariance(uv1::UncertainValue, uv2::UncertainValue, correlation::Real)

Computes the covariance given the correlation coefficient between two UncertainValue.
"""
function covariance(uv1::UncertainValue, uv2::UncertainValue, correlation::Real)
    @assert (correlation >= -1.0) & (correlation <= 1.0) "The Pearson correlation coefficient must be on [-1.0, 1.0] - $correlation"
    return correlation * σ(uv1) * σ(uv2)
end

"""
    divide(n::UncertainValue, d::UncertainValue, cc::AbstractFloat)

Computes `n/d` where `n` and `d` are UncertainValue and `cc` is the correlation coefficient
defined as `cc = covar(n,d)/( σ(n), σ(d) )`
"""
function divide(a::UncertainValue, b::UncertainValue, cc::Float64=0.0)
    @assert (cc >= -1.0) & (cc <= 1.0) "The Pearson correlation coefficient must be on [-1.0, 1.0] - $cc"
    f, sab = a.value / b.value, covariance(a, b, cc)
    return UncertainValue(
        f, abs(1.0/b.value)*sqrt(a.sigma^2 + ((b.sigma*a.value) / b.value)^2 - 2.0 * sab * a.value / b.value)
    )
end

"""
    multiply(a::UncertainValue, b::UncertainValue, cc::AbstractFloat)

Computes `a*b` where `a` and `b` are UncertainValue and `cc` is the correlation coefficient
defined as `cc = covar(a,b)/( σ(a), σ(b) )`
"""
function multiply(a::UncertainValue, b::UncertainValue, cc::Float64=0.0)
    @assert (cc >= -1.0) & (cc <= 1.0) "The Pearson correlation coefficient must be on [-1.0, 1.0] - $cc"
    f, sab = a.value * b.value, covariance(a, b, cc)
    x, y = a.sigma * b.value, b.sigma*a.value
    return UncertainValue(f, sqrt(x^2 + y^2 + 2.0 * cc * x * y))
end

"""
    add(ka::Real, a::UncertainValue, kb::Real, b::UncertainValue, cc::AbstractFloat)

Computes `ka*a + kb*b` where `a` and `b` are UncertainValue and `cc` is the correlation
coefficient defined as `cc = covar(a,b)/( σ(a), σ(b) )`
"""
add(ka::Real, a::UncertainValue, kb::Real, b::UncertainValue, cc::AbstractFloat) =
    UncertainValue(
        convert(Float64, ka) * a.value + convert(Float64, kb) * b.value,
        sqrt((ka * a.sigma)^2 + (kb * b.sigma)^2 + 2.0 * ka * kb * covariance(a, b, cc)),
    )

Base.:*(a::Real, B::UncertainValue) =
    UncertainValue(convert(Float64, a) * B.value, abs(convert(Float64, a) * B.sigma))
Base.:*(A::UncertainValue, b::Real) =
    UncertainValue(convert(Float64, b) * A.value, abs(convert(Float64, b) * A.sigma))

Base.:-(A::UncertainValue) = UncertainValue(-A.value, A.sigma)
Base.:+(A::UncertainValue) = UncertainValue(A.value, A.sigma)

Base.abs(A::UncertainValue) = A.value >= 0.0 ? A : UncertainValue(abs(A.value), A.sigma)

Base.:/(a::Real, B::UncertainValue) = UncertainValue(
    convert(Float64, a) / B.value,
    abs((convert(Float64, a) * B.sigma) / (B.value * B.value)),
)
Base.:/(A::UncertainValue, b::Real) =
    UncertainValue(A.value / convert(Float64, b), abs(A.sigma / convert(Float64, b)))

Base.inv(B::UncertainValue) = UncertainValue(1.0 / B.value, abs(B.sigma / (B.value^2)))

Base.:^(A::UncertainValue, b::Real) =
    UncertainValue(A.value^b, abs(A.value^(b - one(b)) * b * A.sigma))

Base.log(A::UncertainValue) = UncertainValue(log(A.value), abs(A.sigma / A.value))

Base.exp(A::UncertainValue) = UncertainValue(exp(A.value), exp(A.value) * A.sigma)

Base.sin(A::UncertainValue) = UncertainValue(sin(A.value), abs(cos(A.value) * A.sigma))

Base.cos(A::UncertainValue) = UncertainValue(cos(A.value), abs(sin(A.value) * A.sigma))

Base.tan(A::UncertainValue) = UncertainValue(tan(A.value), sec(A.value)^2 * A.sigma)

Base.sqrt(A::UncertainValue) =
    UncertainValue(sqrt(A.value), abs(A.value^(-0.5) * 0.5 * A.sigma))

Base.isequal(A::UncertainValue, B::UncertainValue)::Bool =
    isequal(A.value, B.value) && isequal(A.sigma, B.sigma)

Base.isapprox(
    A::UncertainValue,
    B::UncertainValue;
    atol::Float64 = 0.0,
    rtol::Float64 = atol == 0.0 ? 0.0 : √eps(Float64),
)::Bool =
    isapprox(A.value, B.value, rtol = rtol, atol = atol) &&
    isapprox(A.sigma, B.sigma, rtol = rtol, atol = atol)

"""
    equivalent(A::UncertainValue, B::UncertainValue, k=1.0)

A and B are equivalent if |A-B| <= k σ(A-B)
"""
equivalent(A::UncertainValue, B::UncertainValue, k = 1.0) =
    abs(A.value - B.value) <= k * sqrt(A.sigma^2 + B.sigma^2)

function Base.parse(::Type{UncertainValue}, str::AbstractString)::UncertainValue
    sp = split(str, r"(?:\+\-|\-\+|±)")
    val = parse(Float64, sp[1])
    sigma = length(sp) == 2 ? parse(Float64, sp[2]) : 0.0
    UncertainValue(val, sigma)
end

function _showstr(uv::UncertainValue, pm = "±")
    smin = max(uv.sigma, 1.0e-6 * abs(uv.value), 1.0e-10)
    ls, lr = map(v -> floor(Int, log10(abs(v))), ( smin, max(1.0e-10, uv.value / smin)))
    return if (ls < -4) || (ls > 6)
        if lr < 1
            @sprintf("%0.1e %s %0.1e", uv.value, pm, uv.sigma)
        elseif lr < 2
            @sprintf("%0.2e %s %0.1e", uv.value, pm, uv.sigma)
        elseif lr < 3
            @sprintf("%0.3e %s %0.1e", uv.value, pm, uv.sigma)
        elseif lr < 4
            @sprintf("%0.4e %s %0.1e", uv.value, pm, uv.sigma)
        elseif lr < 5
            @sprintf("%0.5e %s %0.1e", uv.value, pm, uv.sigma)
        else
            @sprintf("%0.6e %s %0.1e", uv.value, pm, uv.sigma)
        end
    else # lv < 0 && lv > -4
        if lr > 6
            @sprintf("%e %s %0.1e", uv.value, pm, uv.sigma)
        elseif ls >= 0
            @sprintf("%0.0f %s %0.0f", uv.value, pm, uv.sigma)
        elseif ls == -1
            @sprintf("%0.2f %s %0.2f", uv.value, pm, uv.sigma)
        elseif ls == -2
            @sprintf("%0.3f %s %0.3f", uv.value, pm, uv.sigma)
        else# if ls==-3
            @sprintf("%0.4f %s %0.4f", uv.value, pm, uv.sigma)
        end
    end
end

Base.show(io::IO, uv::UncertainValue) = print(io, _showstr(uv))

"""
    LaTeXStrings.latexstring(uv::UncertainValue; fmt=nothing, mode=:normal[|:siunitx])::LaTeXString

Converts an `UncertainValue` to a `LaTeXString` in a reasonable manner.
`mode=:siunitx" requires `\\usepackage{siunitx}` which defines `\\num{}`.
`fmt` is a C-style format string like "%0.2f" or nothing for an "intelligent" default.
"""
function LaTeXStrings.latexstring(
    uv::UncertainValue;
    fmt = nothing,
    mode = :normal,
)::LaTeXString
    pre, post = (mode == :siunitx ? ("\\num{", "}") : ("", ""))
    if isnothing(fmt)
        return latexstring(pre * _showstr(uv, raw"\pm") * post)
    else
        fmtrfunc = generate_formatter(fmt)
        return latexstring(
            pre * fmtrfunc(value(uv)) * " \\pm " * fmtrfunc(uncertainty(uv)) * post,
        )
    end
end


text:
#
# ICU - International Components for Unicode
#

module ICU

using Compat

export set_locale

export U_FAILURE,
       U_SUCCESS,

       # ubrk
       UBRK_CHARACTER,
       UBRK_LINE,
       UBRK_SENTENCE,
       UBRK_TITLE,
       UBRK_WORD,
       UBreakIterator,
       ubrk_close,
       ubrk_next,
       ubrk_open,
       ubrk_setUText,

       # ucasemap
       ucasemap_utf8FoldCase,
       ucasemap_utf8ToLower,
       ucasemap_utf8ToTitle,
       ucasemap_utf8ToUpper,

       # ucnv
       UConverter,
       ucnv_close,
       ucnv_convertEx,
       ucnv_open,
       ucnv_toUChars,

       # ucol
       UCollator,
       ucol_close,
       ucol_open,
       ucol_strcoll,

       # ustring
       u_strFoldCase,
       u_strToLower,
       u_strToUpper,
       u_strToTitle,

       # utext
       UText,
       utext_close,
       utext_open

include("../deps/deps.jl")
include("../deps/versions.jl")

@windows_only begin
    # make sure versions match
    v1 = last(matchall(r"\d{2}", iculib))
    v2 = last(matchall(r"\d{2}", iculibi18n))
    if v1 != v2
        error("ICU library version mismatch -- please correct $(realpath("../deps/deps.jl"))")
    end
end

dliculib = Libdl.dlopen(iculib)
for (suffix,version) in [("",0);
                         [("_$i",i) for i in versions];
                         [("_$(string(i)[1])_$(string(i)[2])",i) for i in versions]]
    if Libdl.dlsym_e(dliculib, "u_strToUpper"*suffix) != C_NULL
        @eval const version = $version
        for f in (:u_strFoldCase,
                  :u_strToLower,
                  :u_strToTitle,
                  :u_strToUpper,
                  :ubrk_close,
                  :ubrk_next,
                  :ubrk_open,
                  :ubrk_setUText,
                  :ucal_add,
                  :ucal_clear,
                  :ucal_close,
                  :ucal_get,
                  :ucal_getDefaultTimeZone,
                  :ucal_getTimeZoneDisplayName,
                  :ucal_getMillis,
                  :ucal_getNow,
                  :ucal_open,
                  :ucal_set,
                  :ucal_setDate,
                  :ucal_setDateTime,
                  :ucal_setMillis,
                  :ucasemap_open,
                  :ucasemap_close,
                  :ucasemap_getBreakIterator,
                  :ucasemap_utf8FoldCase,
                  :ucasemap_utf8ToLower,
                  :ucasemap_utf8ToTitle,
                  :ucasemap_utf8ToUpper,
                  :ucnv_close,
                  :ucnv_convertEx,
                  :ucnv_open,
                  :ucnv_toUChars,
                  :ucol_close,
                  :ucol_open,
                  :ucol_strcoll,
                  :ucol_strcollUTF8,
                  :udat_close,
                  :udat_format,
                  :udat_open,
                  :udat_parse,
                  :utext_close,
                  :utext_openUChars,
                  :utext_openUTF8)
            @eval const $(symbol(string('_',f))) = $(string(f,suffix))
        end
        break
    end
end

typealias UBool Int8
typealias UChar UInt16
typealias UErrorCode Int32

U_FAILURE(x::Int32) = x > 0
U_SUCCESS(x::Int32) = x <= 0
U_BUFFER_OVERFLOW_ERROR = 15

locale = C_NULL
casemap = C_NULL
collator = C_NULL

typealias LocaleString @compat Union{ASCIIString,Ptr{Void}}

function set_locale(loc::LocaleString)
    global locale, casemap, collator
    if casemap != C_NULL
        ccall((_ucasemap_close, iculib), Void, (Ptr{Void},), casemap)
        casemap = C_NULL
    end
    if collator != C_NULL
        ucol_close(collator)
        collator = C_NULL
    end
    err = UErrorCode[0]
    casemap = ccall((_ucasemap_open, iculib), Ptr{Void},
                    (Ptr{UInt8}, Int32, Ptr{UErrorCode}),
                    loc, 0, err)
    U_FAILURE(err[1]) && error("could not set casemap")
    collator = ucol_open(loc)
    locale = loc
end

## utext ##

immutable UText
    p::Ptr{Void}
end

function utext_open(str::UTF8String)
    err = UErrorCode[0]
    p = ccall((_utext_openUTF8, iculib), Ptr{Void},
              (Ptr{Void}, Ptr{UInt8}, Int64, Ptr{UErrorCode}),
              C_NULL, str.data, sizeof(str.data), err)
    @assert U_SUCCESS(err[1])
    UText(p)
end

function utext_open(str::UTF16String)
    err = UErrorCode[0]
    p = ccall((_utext_openUChars, iculib), Ptr{Void},
              (Ptr{Void}, Ptr{UChar}, Int64, Ptr{UErrorCode}),
              C_NULL, str.data, length(str.data)-1, err)
    @assert U_SUCCESS(err[1])
    UText(p)
end

function utext_close(t::UText)
    if (t.p != C_NULL)
        ccall((_utext_close, iculib), Void, (Ptr{Void},), t.p)
        t.p = C_NULL
    end
end

## ustring ##

for f in [:u_strToLower, :u_strToUpper]
    @eval begin
        function ($f)(s::UTF16String)
            src = s.data
            destsiz = @compat Int32(length(src))
            dest = zeros(UInt16, destsiz)
            err = UErrorCode[0]
            n = ccall(($(symbol(string('_',f))), iculib), Int32,
                      (Ptr{UChar}, Int32, Ptr{UChar}, Int32, Ptr{UInt8}, Ptr{UErrorCode}),
                      dest, destsiz, src, length(src)-1, locale, err)
            # Retry with large enough buffer if got buffer overflow
            if err[1] == U_BUFFER_OVERFLOW_ERROR
                err[1] = 0
                destsiz = n + 1
                dest = zeros(UInt16, destsiz)
                n = ccall(($(symbol(string('_',f))), iculib), Int32,
                          (Ptr{UChar}, Int32, Ptr{UChar}, Int32, Ptr{UInt8}, Ptr{UErrorCode}),
                          dest, destsiz, src, length(src)-1, locale, err)
            end
            U_FAILURE(err[1]) && error("failed to map case")
            return UTF16String(dest[1:n+1])
        end
    end
end

function u_strFoldCase(s::UTF16String)
    src = s.data
    destsiz = @compat Int32(length(src))
    dest = zeros(UInt16, destsiz)
    err = UErrorCode[0]
    n = ccall((_u_strFoldCase, iculib), Int32,
              (Ptr{UChar}, Int32, Ptr{UChar}, Int32, UInt32, Ptr{UErrorCode}),
              dest, destsiz, src, length(src)-1, 0, err)
    # Retry with large enough buffer if got buffer overflow
    if err[1] == U_BUFFER_OVERFLOW_ERROR
        err[1] = 0
        destsiz = n + 1
        dest = zeros(UInt16, destsiz)
        n = ccall((_u_strFoldCase, iculib), Int32,
                  (Ptr{UChar}, Int32, Ptr{UChar}, Int32, UInt32, Ptr{UErrorCode}),
                  dest, destsiz, src, length(src)-1, 0, err)
    end
    U_FAILURE(err[1]) && error("failed to map case")
    return UTF16String(dest[1:n+1])
end

function u_strToTitle(s::UTF16String)
    src = s.data
    destsiz = @compat Int32(length(src))
    dest = zeros(UInt16, destsiz)
    err = UErrorCode[0]
    breakiter = ccall((_ucasemap_getBreakIterator, iculib), Ptr{Void}, (Ptr{Void},), casemap)
    n = ccall((_u_strToTitle, iculib), Int32,
              (Ptr{UChar}, Int32, Ptr{UChar}, Int32, Ptr{Void}, Ptr{UInt8}, Ptr{UErrorCode}),
              dest, destsiz, src, length(src)-1, breakiter, locale, err)
    # Retry with large enough buffer if got buffer overflow
    if err[1] == U_BUFFER_OVERFLOW_ERROR
        err[1] = 0
        destsiz = n + 1
        dest = zeros(UInt16, destsiz)
        n = ccall((_u_strToTitle, iculib), Int32,
                  (Ptr{UChar}, Int32, Ptr{UChar}, Int32, Ptr{Void}, Ptr{UInt8}, Ptr{UErrorCode}),
                  dest, destsiz, src, length(src)-1, breakiter, locale, err)
    end
    U_FAILURE(err[1]) && error("failed to map case")
    return UTF16String(dest[1:n+1])
end

## ucasemap ##

for f in [:ucasemap_utf8FoldCase,
          :ucasemap_utf8ToLower,
          :ucasemap_utf8ToTitle,
          :ucasemap_utf8ToUpper]
    @eval begin
        function ($f)(src::UTF8String)
            destsiz = @compat Int32(sizeof(src.data)+1)
            dest = zeros(Cchar, destsiz)
            err = UErrorCode[0]
            n = ccall(($(symbol(string('_',f))), iculib), Int32,
                      (Ptr{Void}, Ptr{Cchar}, Int32, Ptr{Cchar}, Int32, Ptr{UErrorCode}),
                      casemap, dest, destsiz, src.data, sizeof(src.data), err)
            # Retry with large enough buffer if got buffer overflow
            if err[1] == U_BUFFER_OVERFLOW_ERROR
                err[1] = 0
                destsiz = n + 1
                dest = zeros(Cchar, destsiz)
                n = ccall(($(symbol(string('_',f))), iculib), Int32,
                          (Ptr{Void}, Ptr{Cchar}, Int32, Ptr{Cchar}, Int32, Ptr{UErrorCode}),
                          casemap, dest, destsiz, src.data, sizeof(src.data), err)
            end
            U_FAILURE(err[1]) && error("failed to map case")
            return utf8(dest[1:n])
        end
    end
end

## ubrk ##

typealias UBreakIteratorType Int32
const UBRK_CHARACTER = @compat Int32(0)
const UBRK_WORD = @compat Int32(1)
const UBRK_LINE = @compat Int32(2)
const UBRK_SENTENCE = @compat Int32(3)
const UBRK_TITLE = @compat Int32(4)

type UBreakIterator
    p::Ptr{Void}
    UBreakIterator(p::Ptr) = (self = new(p); finalizer(self, close); self)
end

function ubrk_open(kind::Integer, loc::LocaleString, s::Vector{UInt16})
    err = UErrorCode[0]
    p = ccall((_ubrk_open, iculib), Ptr{Void},
              (UBreakIteratorType, Ptr{UInt8}, Ptr{UChar}, Int32, Ptr{UErrorCode}),
              kind, loc, s, length(s), err)
    @assert U_SUCCESS(err[1])
    UBreakIterator(p)
end

function ubrk_open(kind::Integer, loc::LocaleString)
    err = UErrorCode[0]
    p = ccall((_ubrk_open, iculib), Ptr{Void},
              (UBreakIteratorType, Ptr{UInt8}, Ptr{UChar}, Int32, Ptr{UErrorCode}),
              kind, loc, C_NULL, 0, err)
    @assert U_SUCCESS(err[1])
    UBreakIterator(p)
end

function ubrk_close(bi::UBreakIterator)
    if bi.p != C_NULL
        ccall((_ubrk_close, iculib), Void, (Ptr{Void},), bi.p)
        bi.p = C_NULL
    end
end

ubrk_next(bi::UBreakIterator) =
    ccall((_ubrk_next, iculib), Int32, (Ptr{Void},), bi.p)

function ubrk_setUText(bi::UBreakIterator, t::UText)
    err = UErrorCode[0]
    ccall((_ubrk_setUText, iculib), Void,
          (Ptr{Void}, Ptr{Void}, Ptr{UErrorCode}),
          bi.p, t.p, err)
    @assert U_SUCCESS(err[1])
    nothing
end

## ucnv ##

type UConverter
    p::Ptr{Void}
    UConverter(p::Ptr) = (self = new(p); finalizer(self, close); self)
end

ucnv_close(c::UConverter) =
    ccall((_ucnv_close, iculibi18n), Void, (Ptr{Void},), c.p)

type UConverterPivot
    buf::Vector{UChar}
    pos::Vector{Ptr{UChar}}

    function UConverterPivot(n::Int)
        buf = Array(UChar, n)
        p = pointer(buf)
        new(buf, [p,p])
    end
end

function ucnv_convertEx(dstcnv::UConverter, srccnv::UConverter,
                        dst::IOBuffer, src::IOBuffer, pivot::UConverterPivot,
                        reset::Bool=false, flush::Bool=true)
    p = Ptr{UInt8}[pointer(dst.data, position(dst)+1),
                   pointer(src.data, position(src)+1)]
    p0 = copy(p)
    err = UErrorCode[0]
    ccall((_ucnv_convertEx, iculibi18n), Void,
          (Ptr{Void}, Ptr{Void},
           Ptr{Ptr{UInt8}}, Ptr{UInt8}, Ptr{Ptr{UInt8}}, Ptr{UInt8},
           Ptr{UChar}, Ptr{Ptr{UChar}}, Ptr{Ptr{UChar}}, Ptr{UChar},
           UBool, UBool, Ptr{UErrorCode}),
          dstcnv.p, srccnv.p,
          pointer(p, 1), pointer(dst.data, length(dst.data)+1),
          pointer(p, 2), pointer(src.data, src.size+1),
          pointer(pivot.buf, 1),
          pointer(pivot.pos, 1),
          pointer(pivot.pos, 2),
          pointer(pivot.buf, length(pivot.buf)+1),
          reset, flush, err)
    dst.size += p[1] - p0[1]
    dst.ptr += p[1] - p0[1]
    src.ptr += p[2] - p0[2]
    err[1] == U_BUFFER_OVERFLOW_ERROR && return true
    @assert U_SUCCESS(err[1])
    false
end

function ucnv_open(name::ASCIIString)
    err = UErrorCode[0]
    p = ccall((_ucnv_open, iculibi18n), Ptr{Void},
              (Cstring, Ptr{UErrorCode}),
              name, err)
    U_SUCCESS(err[1]) || error("ICU: could not open converter ", name)
    UConverter(p)
end

function ucnv_toUChars(cnv::UConverter, b::Vector{UInt8})
    u = zeros(UInt16, 2*length(b))
    err = UErrorCode[0]
    n = ccall((_ucnv_toUChars, iculibi18n), Int32,
              (Ptr{Void}, Ptr{UChar}, Int32, Ptr{UInt8}, Int32, Ptr{UErrorCode}),
              cnv.p, u, length(u), b, length(b), err)
    U_SUCCESS(err[1]) || error("ICU: could not open converter ", name)
    UTF16String(u[1:n+1])
end

## ucol ##

type UCollator
    p::Ptr{Void}
    UCollator(p::Ptr) = (self = new(p); finalizer(self, close); self)
end

function ucol_close(c::UCollator)
    if c.p != C_NULL
        ccall((_ucol_close, iculibi18n), Void, (Ptr{Void},), c.p)
        c.p = C_NULL
    end
end

function ucol_open(loc::LocaleString)
    err = UErrorCode[0]
    p = ccall((_ucol_open, iculibi18n), Ptr{Void},
              (Ptr{UInt8}, Ptr{UErrorCode}),
              loc, err)
    U_SUCCESS(err[1]) || error("ICU: could not open collator for locale ", locale)
    UCollator(p)
end

function ucol_strcoll(c::UCollator, a::UTF8String, b::UTF8String)
    err = UErrorCode[0]
    o = ccall((_ucol_strcollUTF8, iculibi18n), Int32,
              (Ptr{Void}, Ptr{UInt8}, Int32, Ptr{UInt8}, Int32, Ptr{UErrorCode}),
              c.p, a.data, sizeof(a.data), b.data, sizeof(b.data), err)
    @assert U_SUCCESS(err[1])
    o
end

function ucol_strcoll(c::UCollator, a::UTF16String, b::UTF16String)
    err = UErrorCode[0]
    o = ccall((_ucol_strcoll, iculibi18n), Int32,
              (Ptr{Void}, Ptr{UChar}, Int32, Ptr{UChar}, Int32, Ptr{UErrorCode}),
              c.p, a.data, length(a.data)-1, b.data, length(b.data)-1, err)
    @assert U_SUCCESS(err[1])
    o
end

## calendar ##

export ICUCalendar,
       add,
       clear,
       get,
       getDefaultTimeZone,
       getMillis,
       getNow,
       getTimeZoneDisplayName,
       set,
       setDate,
       setDateTime,
       setMillis

export UCAL_ERA,
       UCAL_YEAR,
       UCAL_MONTH,
       UCAL_WEEK_OF_YEAR,
       UCAL_WEEK_OF_MONTH,
       UCAL_DATE,
       UCAL_DAY_OF_YEAR,
       UCAL_DAY_OF_WEEK,
       UCAL_DAY_OF_WEEK_IN_MONTH,
       UCAL_AM_PM,
       UCAL_HOUR,
       UCAL_HOUR_OF_DAY,
       UCAL_MINUTE,
       UCAL_SECOND,
       UCAL_MILLISECOND,
       UCAL_ZONE_OFFSET,
       UCAL_DST_OFFSET,
       UCAL_YEAR_WOY,
       UCAL_DOW_LOCAL,
       UCAL_EXTENDED_YEAR,
       UCAL_JULIAN_DAY,
       UCAL_MILLISECONDS_IN_DAY,
       UCAL_IS_LEAP_MONTH

for (i,a) in enumerate([
        :UCAL_ERA,
        :UCAL_YEAR,
        :UCAL_MONTH,
        :UCAL_WEEK_OF_YEAR,
        :UCAL_WEEK_OF_MONTH,
        :UCAL_DATE,
        :UCAL_DAY_OF_YEAR,
        :UCAL_DAY_OF_WEEK,
        :UCAL_DAY_OF_WEEK_IN_MONTH,
        :UCAL_AM_PM,
        :UCAL_HOUR,
        :UCAL_HOUR_OF_DAY,
        :UCAL_MINUTE,
        :UCAL_SECOND,
        :UCAL_MILLISECOND,
        :UCAL_ZONE_OFFSET,
        :UCAL_DST_OFFSET,
        :UCAL_YEAR_WOY,
        :UCAL_DOW_LOCAL,
        :UCAL_EXTENDED_YEAR,
        :UCAL_JULIAN_DAY,
        :UCAL_MILLISECONDS_IN_DAY,
        :UCAL_IS_LEAP_MONTH
    ])
    @eval const $a = @compat Int32($i - 1)
end

typealias UDate Float64

type ICUCalendar
    ptr::Ptr{Void}
    ICUCalendar(p::Ptr) = (self = new(p); finalizer(self, close); self)
end

ICUCalendar(timezone::AbstractString) = ICUCalendar(utf16(timezone))
function ICUCalendar(tz::UTF16String)
    err = UErrorCode[0]
    p = ccall((_ucal_open, iculibi18n), Ptr{Void},
              (Ptr{UChar}, Int32, Ptr{UInt8}, Int32, Ptr{UErrorCode}),
              tz.data, length(tz.data)-1, locale, 0, err)
    ICUCalendar(p)
end
function ICUCalendar()
    err = UErrorCode[0]
    p = ccall((_ucal_open, iculibi18n), Ptr{Void},
              (Ptr{UChar}, Int32, Ptr{UInt8}, Int32, Ptr{UErrorCode}),
              C_NULL, 0, locale, 0, err)
    ICUCalendar(p)
end

function close(cal::ICUCalendar)
    if cal.ptr != C_NULL
        ccall((_ucal_close, iculibi18n), Void, (Ptr{Void},), cal.ptr)
        cal.ptr = C_NULL
    end
end

getNow() = ccall((_ucal_getNow, iculibi18n), UDate, ())

function getMillis(cal::ICUCalendar)
    err = UErrorCode[0]
    ccall((_ucal_getMillis, iculibi18n), UDate, (Ptr{Void}, Ptr{UErrorCode}), cal.ptr, err)
end

function setMillis(cal::ICUCalendar, millis::UDate)
    err = UErrorCode[0]
    ccall((_ucal_setMillis, iculibi18n), Void,
          (Ptr{Void}, UDate, Ptr{UErrorCode}),
          cal.ptr, millis, err)
end

function setDate(cal::ICUCalendar, y::Integer, m::Integer, d::Integer)
    err = UErrorCode[0]
    ccall((_ucal_setDate, iculibi18n), Void,
          (Ptr{Void}, Int32, Int32, Int32, Ptr{UErrorCode}),
          cal.ptr, y, m-1, d, err)
end

function setDateTime(cal::ICUCalendar, y::Integer, mo::Integer, d::Integer, h::Integer, mi::Integer, s::Integer)
    err = UErrorCode[0]
    ccall((_ucal_setDateTime, iculibi18n), Void,
          (Ptr{Void}, Int32, Int32, Int32, Int32, Int32, Int32, Ptr{UErrorCode}),
          cal.ptr, y, mo-1, d, h, mi, s, err)
end

function clear(cal::ICUCalendar)
    err = UErrorCode[0]
    ccall((_ucal_clear, iculibi18n), Void,
          (Ptr{Void}, Ptr{UErrorCode}),
          cal.ptr, err)
end

function get(cal::ICUCalendar, field::Int32)
    err = UErrorCode[0]
    ccall((_ucal_get, iculibi18n), Int32,
          (Ptr{Void},Int32,Ptr{UErrorCode}),
          cal.ptr, field, err)
end
get(cal::ICUCalendar, fields::Vector{Int32}) = [get(cal,f) for f in fields]

function add(cal::ICUCalendar, field::Int32, amount::Integer)
    err = UErrorCode[0]
    ccall((_ucal_add, iculibi18n), Int32,
          (Ptr{Void},Int32,Int32,Ptr{UErrorCode}),
          cal.ptr, field, amount, err)
end

function set(cal::ICUCalendar, field::Int32, val::Integer)
             ccall((_ucal_set, iculibi18n), Void,
             (Ptr{Void}, Int32, Int32),
             cal.ptr, field, val)
end

function getTimeZoneDisplayName(cal::ICUCalendar)
    bufsz = 64
    buf = zeros(UInt16, bufsz)
    err = UErrorCode[0]
    len = ccall((_ucal_getTimeZoneDisplayName, iculibi18n), Int32,
                (Ptr{Void}, Int32, Ptr{UInt8}, Ptr{UChar}, Int32, Ptr{UErrorCode}),
                cal.ptr, 1, locale, buf, bufsz, err)
    UTF16String(buf[1:len+1])
end

function getDefaultTimeZone()
    bufsz = 64
    buf = zeros(UInt16, bufsz)
    err = UErrorCode[0]
    len = ccall((_ucal_getDefaultTimeZone, iculibi18n), Int32,
                (Ptr{UChar}, Int32, Ptr{UErrorCode}),
                buf, bufsz, err)
    UTF16String(buf[1:len+1])
end

export ICUDateFormat,
       format,
       parse

export UDAT_NONE,
       UDAT_FULL,
       UDAT_LONG,
       UDAT_MEDIUM,
       UDAT_SHORT,
       UDAT_RELATIVE

const UDAT_NONE     = @compat Int32(-1)
const UDAT_FULL     = @compat Int32(0)
const UDAT_LONG     = @compat Int32(1)
const UDAT_MEDIUM   = @compat Int32(2)
const UDAT_SHORT    = @compat Int32(3)
const UDAT_RELATIVE = @compat Int32(1<<7)

type ICUDateFormat
    ptr::Ptr{Void}
    ICUDateFormat(p::Ptr) = (self = new(p); finalizer(self, close); self)
end

ICUDateFormat(pattern::AbstractString, tz::AbstractString) = ICUDateFormat(utf16(pattern), utf16(tz))
function ICUDateFormat(pattern::UTF16String, tz::UTF16String)
    err = UErrorCode[0]
    p = ccall((_udat_open, iculibi18n), Ptr{Void},
              (Int32, Int32, Ptr{UInt8}, Ptr{UChar}, Int32, Ptr{UChar}, Int32, Ptr{UErrorCode}),
              -2, -2, locale, tz.data, length(tz.data)-1,
              pattern.data, length(pattern.data)-1, err)
    U_FAILURE(err[1]) && error("bad date format")
    ICUDateFormat(p)
end

ICUDateFormat(tstyle::Integer, dstyle::Integer, tz::AbstractString) = ICUDateFormat(tstyle, dstyle, utf16(tz))
function ICUDateFormat(tstyle::Integer, dstyle::Integer, tz::UTF16String)
    err = UErrorCode[0]
    p = ccall((_udat_open, iculibi18n), Ptr{Void},
              (Int32, Int32, Ptr{UInt8}, Ptr{UChar}, Int32, Ptr{UChar}, Int32, Ptr{UErrorCode}),
              tstyle, dstyle, locale, tz.data, length(tz.data)-1, C_NULL, 0, err)
    U_FAILURE(err[1]) && error("bad date format")
    ICUDateFormat(p)
end

close(df::ICUDateFormat) =
    ccall((_udat_close, iculibi18n), Void, (Ptr{Void},), df.ptr)

function format(df::ICUDateFormat, millis::Float64)
    err = UErrorCode[0]
    buflen = 64
    buf = zeros(UChar, buflen)
    len = ccall((_udat_format, iculibi18n), Int32,
                (Ptr{Void}, Float64, Ptr{UChar}, Int32, Ptr{Void}, Ptr{UErrorCode}),
                df.ptr, millis, buf, buflen, C_NULL, err)
    U_FAILURE(err[1]) && error("failed to format time")
    UTF16String(buf[1:len+1])
end

parse(df::ICUDateFormat, s::AbstractString) = parse(df, utf16(s))
function parse(df::ICUDateFormat, s16::UTF16String)
    err = UErrorCode[0]
    ret = ccall((_udat_parse, iculibi18n), Float64,
                (Ptr{Void}, Ptr{UChar}, Int32, Ptr{Int32}, Ptr{UErrorCode}),
                df.ptr, s16.data, length(s16.data)-1, C_NULL, err)
    U_FAILURE(err[1]) && error("failed to parse string")
    ret
end

function test_icucalendar()
    cal = ICUCalendar("America/Los_Angeles")
    setMillis(cal, getNow())
    fields = [UCAL_YEAR, UCAL_MONTH, UCAL_DATE,
              UCAL_HOUR_OF_DAY, UCAL_MINUTE, UCAL_SECOND]
    println(get(cal,fields))
    clear(cal)
    df = ICUDateFormat()
    s = format(df, getNow())
    show(s)
end

## init ##

set_locale(locale)

end # module

ext: jl


text:
using Blosc2
import Blosc2.Lib as lb
@testset "init destroy" begin
    ctx = lb.blosc2_create_cctx()
    lb.blosc2_free_ctx(ctx)
    @test 1 == 1

    ctx = lb.blosc2_create_dctx()
    lb.blosc2_free_ctx(ctx)
    @test 1 == 1
end

@testset "blosc2_compress_ctx / decompress" begin
    ctx = lb.blosc2_create_cctx(lb.blosc2_cparams(typesize = sizeof(Int64)))
    n = 10000
    data = rand(1:1000, n)
    buff = Vector{UInt8}(undef, sizeof(data) + lb.BLOSC_MAX_OVERHEAD)
    sz = lb.blosc2_compress_ctx(ctx, data, sizeof(data), buff, sizeof(buff))
    @test sz > 0
    lb.blosc2_free_ctx(ctx)

    dctx = lb.blosc2_create_dctx()

    res = Vector{Int64}(undef, n + 500)
    sz = lb.blosc2_decompress_ctx(dctx, buff, sizeof(buff), res, sizeof(res))
    @test sz == sizeof(data)
    @test data == res[1:div(sz, sizeof(Int64))]

    lb.blosc2_free_ctx(dctx)
end

@testset "blosc2_getitem_ctx" begin
    ctx = lb.blosc2_create_cctx(lb.blosc2_cparams(typesize = sizeof(Int64)))
    n = 10000
    data = rand(1:1000, n)
    buff = Vector{UInt8}(undef, sizeof(data) + lb.BLOSC_MAX_OVERHEAD)
    sz = lb.blosc2_compress_ctx(ctx, data, sizeof(data), buff, sizeof(buff))
    @test sz > 0
    lb.blosc2_free_ctx(ctx)

    dctx = lb.blosc2_create_dctx()

    res = Int64[1, 2]

    for i = 1 : div(n, 2)
        sz = lb.blosc2_getitem_ctx(dctx, buff, sizeof(buff), (i - 1) * 2, 2, res, sizeof(res))
        @test sz == sizeof(res)
        @test res == data[(i - 1)*2 + 1:i*2]
    end
    lb.blosc2_free_ctx(dctx)
end

@testset "blosc2_set_maskout" begin
    ctx = lb.blosc2_create_cctx(lb.blosc2_cparams(typesize = sizeof(Int64), blocksize = 80000))
    n = 100000
    data = rand(1:1000, n)
    buff = Vector{UInt8}(undef, sizeof(data) + lb.BLOSC_MAX_OVERHEAD)
    sz = lb.blosc2_compress_ctx(ctx, data, sizeof(data), buff, sizeof(buff))
    @test sz > 0
    lb.blosc2_free_ctx(ctx)

    dctx = lb.blosc2_create_dctx()

    mask = Vector{Bool}(undef, 10)
    fill!(mask, false)
    mask[5:end] .= true
    r = lb.blosc2_set_maskout(dctx, mask, sizeof(mask))

    res = Vector{Int64}(undef, n + 500)
    fill!(res, 0)
    sz = lb.blosc2_decompress_ctx(dctx, buff, sizeof(buff), res, sizeof(res))
    @test sz == sizeof(data)
    int_block = div(80000,sizeof(Int))
    @test data[1:int_block*4] == res[1:int_block*4]
    @test sum(res[int_block*4 + 1 : end]) == 0

    lb.blosc2_free_ctx(dctx)
end


text:
using Test, SparseArrays
using Recommenders: tfidf, bm25, compute_similarity, predict_u2i

@testset "TF-IDF" begin
    inter = sparse([1, 1, 1, 2], [1, 2, 3, 3], [1.0, 3.0, 5.0, 2.0])
    idf = [3 / (3 + 1), 3 / (1 + 1)]
    idf = log.(idf)
    expected =
        sparse([1, 1, 1, 2], [1, 2, 3, 3], [1 * idf[1], 3 * idf[1], 5 * idf[1], 2 * idf[2]])
    evaluated = tfidf(inter)
    @test evaluated == expected
end

@testset "BM25" begin
    k1 = 1.2
    b = 0.75
    inter = sparse([1, 1, 1, 2], [1, 2, 3, 3], [1.0, 3.0, 5.0, 2.0])
    idf = [(3 - 3 + 0.5) / (3 + 0.5), (3 - 1 + 0.5) / (1 + 0.5)]
    idf = log.(idf)
    dl = [1, 1, 2]
    avgdl = 4 / 3
    dl_factor = [
        (k1 + 1) / (1 + k1 * (1 - b + b * dl[1] / avgdl)),
        (k1 + 1) / (1 + k1 * (1 - b + b * dl[2] / avgdl)),
        (k1 + 1) / (1 + k1 * (1 - b + b * dl[3] / avgdl)),
    ]
    expected = sparse(
        [1, 1, 1, 2],
        [1, 2, 3, 3],
        [
            1 * idf[1] * dl_factor[1],
            3 * idf[1] * dl_factor[2],
            5 * idf[1] * dl_factor[3],
            2 * idf[2] * dl_factor[3],
        ],
    )
    evaluated = bm25(inter, k1, b)
    @test evaluated ≈ expected
end

@testset "Compute similarity matrix" begin
    X = sparse([1, 1, 2, 2], [1, 2, 1, 2], [1.0, 2.0, 2.0, 1.0])
    expected = sparse([2, 1], [1, 2], [4 / (5 + 1e-6), 4 / (5 + 1e-6)])
    evaluated = compute_similarity(X, 1, 0.0, true, false)
    @test evaluated ≈ expected

    evaluated = compute_similarity(X, 1, 0.0, true, true)
    for c in eachcol(evaluated)
        @test sum(c) ≈ 1
    end
end

@testset "Predict user to item" begin
    inter = (userid = [1, 2, 2], itemid = [1, 1, 2], rating = [1.0, 2.0, 1.0])
    similarity = sparse([2, 1], [1, 2], [4 / (5 + 1e-6), 4 / (5 + 1e-6)])
    user_history = sparse([1, 0])
    @test predict_u2i(similarity, user_history, 1) == [2]
    @test predict_u2i(similarity, user_history, 2) == [2] # score=0 なので一個しか返らない
    @test predict_u2i(similarity, user_history, 1, drop_history = true) == [2]
    @test predict_u2i(similarity, user_history, 2, drop_history = true) == [2]
    # check dispatch
    @test predict_u2i(similarity, [1, 0], 1) == [2]
end

@testset "Predict item to item" begin
    similarity = sparse([2, 1], [1, 2], [4 / (5 + 1e-6), 4 / (5 + 1e-6)])
    @test predict_i2i(similarity, 1, 1) == [2]
    @test predict_i2i(similarity, 2, 1) == [1]
    @test predict_i2i(similarity, 1, 2) == [2]
    @test predict_i2i(similarity, 2, 2) == [1]
end


text:
const Particle = Trace

"""
Data structure for particle filters
- effectiveSampleSize(pc :: ParticleContainer)
- normalise!(pc::ParticleContainer)
- consume(pc::ParticleContainer): return incremental likelihood
"""
mutable struct ParticleContainer{T<:Particle, F}
  model :: F
  num_particles :: Int
  vals  :: Array{T}
  logWs :: Array{Float64}  # Log weights (Trace) or incremental likelihoods (ParticleContainer)
  logE  :: Float64           # Log model evidence
  # conditional :: Union{Nothing,Conditional} # storing parameters, helpful for implementing rejuvenation steps
  conditional :: Nothing # storing parameters, helpful for implementing rejuvenation steps
  n_consume :: Int # helpful for rejuvenation steps, e.g. in SMC2
end
ParticleContainer{T}(m) where T = ParticleContainer{T}(m, 0)
function ParticleContainer{T}(m::F, n::Int) where {T, F}
  ParticleContainer{T, F}(m, n, Vector{T}(), Vector{Float64}(), 0.0, nothing, 0)
end

Base.collect(pc :: ParticleContainer) = pc.vals # prev: Dict, now: Array
Base.length(pc :: ParticleContainer)  = length(pc.vals)
Base.similar(pc :: ParticleContainer{T}) where T = ParticleContainer{T}(pc.model, 0)
# pc[i] returns the i'th particle
Base.getindex(pc :: ParticleContainer, i :: Real) = pc.vals[i]

# registers a new x-particle in the container
function Base.push!(pc :: ParticleContainer, p :: Particle)
  pc.num_particles += 1
  push!(pc.vals, p)
  push!(pc.logWs, 0)
  pc
end
Base.push!(pc :: ParticleContainer) = Base.push!(pc, eltype(pc.vals)(pc.model))

function Base.push!(pc :: ParticleContainer, n :: Int, spl :: Sampler, varInfo :: VarInfo)
  vals  = Vector{eltype(pc.vals)}(undef,n)
  logWs = zeros(eltype(pc.logWs), n)
  for i=1:n
    vals[i]  = eltype(pc.vals)(pc.model, spl, varInfo)
  end
  append!(pc.vals, vals)
  append!(pc.logWs, logWs)
  pc.num_particles += n
  pc
end

# clears the container but keep params, logweight etc.
function Base.empty!(pc :: ParticleContainer)
  pc.num_particles = 0
  pc.vals  = Vector{Particle}()
  pc.logWs = Vector{Float64}()
  pc
end

# clones a theta-particle
function Base.copy(pc :: ParticleContainer)
  particles = collect(pc)
  newpc     = similar(pc)
  for p in particles
    newp = fork(p)
    push!(newpc, newp)
  end
  newpc.logE        = pc.logE
  newpc.logWs       = deepcopy(pc.logWs)
  newpc.conditional = deepcopy(pc.conditional)
  newpc.n_consume   = pc.n_consume
  newpc
end

# run particle filter for one step, return incremental likelihood
function Turing.consume(pc :: ParticleContainer)
  @assert pc.num_particles == length(pc)
  # normalisation factor: 1/N
  _, z1      = weights(pc)
  n = length(pc.vals)

  particles = collect(pc)
  num_done = 0
  for i=1:n
    p = pc.vals[i]
    score = Turing.consume(p)
    if score isa Real
      score += getlogp(p.vi)
      resetlogp!(p.vi)
      increase_logweight(pc, i, Float64(score))
    elseif score == Val{:done}
      num_done += 1
    else
      error("[consume]: error in running particle filter.")
    end
  end

  if num_done == length(pc)
    res = Val{:done}
  elseif num_done != 0
    error("[consume]: mis-aligned execution traces, num_particles= $(n), num_done=$(num_done).")
  else
    # update incremental likelihoods
    _, z2      = weights(pc)
    res = increase_logevidence(pc, z2 - z1)
    pc.n_consume += 1
    # res = increase_loglikelihood(pc, z2 - z1)
  end

  res
end

function weights(pc :: ParticleContainer)
  @assert pc.num_particles == length(pc)
  logWs = pc.logWs
  Ws = exp.(logWs .- maximum(logWs))
  logZ = log(sum(Ws)) + maximum(logWs)
  Ws = Ws ./ sum(Ws)
  return Ws, logZ
end

function effectiveSampleSize(pc :: ParticleContainer)
  Ws, _ = weights(pc)
  ess = 1.0 / sum(Ws .^ 2) # sum(Ws) ^ 2 = 1.0, because weights are normalised
end

increase_logweight(pc :: ParticleContainer, t :: Int, logw :: Float64) =
  (pc.logWs[t]  += logw)

increase_logevidence(pc :: ParticleContainer, logw :: Float64) =
  (pc.logE += logw)

function resample!( pc :: ParticleContainer,
                   randcat :: Function = Turing.resample_systematic,
                   ref :: Union{Particle, Nothing} = nothing)
  n1, particles = pc.num_particles, collect(pc)
  @assert n1 == length(particles)

  # resample
  Ws, _ = weights(pc)

  # check that weights are not NaN
  @assert !any(isnan.(Ws))

  n2    = isa(ref, Nothing) ? n1 : n1-1
  indx  = randcat(Ws, n2)

  # fork particles
  empty!(pc)
  num_children = zeros(Int,n1)
  map(i->num_children[i]+=1, indx)
  for i = 1:n1
    is_ref = particles[i] == ref
    p = is_ref ? fork(particles[i], is_ref) : particles[i]
    num_children[i] > 0 && push!(pc, p)
    for k=1:num_children[i]-1
      newp = fork(p, is_ref)
      push!(pc, newp)
    end
  end

  if isa(ref, Particle)
    # Insert the retained particle. This is based on the replaying trick for efficiency
    #  reasons. If we implement PG using task copying, we need to store Nx * T particles!
    push!(pc, ref)
  end

  pc
end

########### Auxilary Functions ###################

# ParticleContainer: particles ==> (weight, results)
function getsample(pc :: ParticleContainer, i :: Int, w :: Float64 = 0.)
  p = pc.vals[i]
  predicts = Sample(p.vi).value
  predicts[:le] = pc.logE
  return Sample(w, predicts)
end

getsample(pc :: ParticleContainer) = begin
  w = pc.logE
  Ws, z = weights(pc)
  s = map((i)->getsample(pc, i, Ws[i]), 1:length(pc))
  return exp.(w), s
end


text:
const IntPair = Tuple{Int,Int}

"""
Edge with a kind index.

Each instance of ``Edge`` has three fields:

- ``s``:  index of the source vertex.
- ``t``:  index of the target vertex.
- ``k``:  the kind index, which can be used to access external data.
"""
struct Edge
    s::Int      # source index
    t::Int      # target index
    k::Int      # kind index

    Edge(s::Int, t::Int, k::Int) = new(s, t, k)
    Edge(e::IntPair, k::Int) = new(e[1], e[2], k)
end

"""
Undirected graph.
"""
mutable struct UGraph
    nv::Int     # number of vertices
    g_elst::GroupedList{Edge}   # edge list (grouped by kind)
    g_ies::GroupedList{Int}     # incoming edge indices (grouped by vertices)
    g_oes::GroupedList{Int}     # outgoing edge indices (grouped by vertices)

    function UGraph(nv::Int, nk::Int, g_elst::GroupedList)
        ne = length(g_elst)
        g_ies = grlist_by(1:ne, nv, i->g_elst[i].t)
        g_oes = grlist_by(1:ne, nv, i->g_elst[i].s)
        new(nv, g_elst, g_ies, g_oes)
    end
end

nvertices(g::UGraph) = g.nv
nedgekinds(g::UGraph) = ngroups(g.g_elst)
nedges(g::UGraph) = length(g.g_elst)

vertices(g::UGraph) = 1:g.nv
edges(g::UGraph) = g.g_elst.elems
edges_of(g::UGraph, k::Int) = group(g.g_elst, k)

indegree_at(g::UGraph, v::Int) = group_length(g.g_ies, v)
outdegree_at(g::UGraph, v::Int) = group_length(g.g_oes, v)
degree_at(g::UGraph, v::Int) = indegree_at(g, v) + outdegree_at(g, v)

inedges_at(g::UGraph, v::Int) = view(edges(g), group(g.g_ies, v))
outedges_at(g::UGraph, v::Int) = view(edges(g), group(g.g_oes, v))

"""
Construct an undirected graph with a list of edges (in vertex pairs).

# Arguments

- ``nv``:       The number of vertices
- ``edges``:    The list of edges. Each element in a vertex pair.

# Keyword arguments

- ``shared_kind``:   Whether all edges share the same kind (default = `false`).
                     When `true`, all edges share the same kind index `1`;
                     otherwise, they are assumed to have different kinds,
                     *i.e.*, the kind index of the  `i`-th edge is `i`.
"""
function UGraph(nv::Int, edges::AbstractVector{IntPair}; shared_kind::Bool=false)
    nv >= 0 || throw(ArgumentError("nv must be non-negative."))

    # make edge list
    ne = length(edges)
    nk = (shared_kind ? 1 : ne)

    elst = Vector{Edge}(ne)
    elst_cnts = Vector{Int}(nk)
    if shared_kind
        for (i, e) in enumerate(edges)
            elst[i] = Edge(e, 1)
        end
        elst_cnts[1] = ne
    else
        for (i, e) in enumerate(edges)
            elst[i] = Edge(e, i)
            elst_cnts[i] = 1
        end
    end
    g_elst = grlist_from_counts(elst, elst_cnts)

    # construct
    UGraph(nv, nk, g_elst)
end

"""
Construct an undirected graph with edges grouped by their kinds.

# Arguments

- nv:      The number of vertices.
- egrps:   A list of edge groups. `gedges[k]` is the list of edges for
           the `k`-th kind.
"""
function UGraph{G<:AbstractVector{IntPair}}(nv::Int, egrps::AbstractVector{G})
    nv >= 0 || throw(ArgumentError("nv must be non-negative."))

    # make edge list
    nk = length(egrps)
    ne = sum(length(g) for g in egrps)
    elst = Vector{Edge}(ne)
    elst_cnts = Vector{Int}(nk)

    i = 0
    for (k, g) in enumerate(egrps)
        for e in g
            elst[i += 1] = Edge(e, k)
        end
        elst_cnts[k] = length(g)
    end
    @assert i == ne
    g_elst = grlist_from_counts(elst, elst_cnts)

    # construct
    UGraph(nv, nk, g_elst)
end

# """
# Convert a tree-structured undirected graph to an instance of Tree.
#
# # Arguments
# - g:  The input graph.
# - r:  The index of the root.
#
# It throws an ArgumentError if `g` is not a connected tree.
# """
# function to_tree(g::UGraph, r::Int)
#     nv = nvertices(g)
#     0 < r <= nv || throw(BoundsError("Root index out of range."))
#     nedges(g) == nv - 1 || throw(
#         BoundsError("Unexpected number of edges for a tree."))
#
#     parents = collect(1:nv)
#     parents_lt = zeros(Int, nv)
#     visited = fill(false, nv)
#
#     # initialize the queue of active vertices
#     que = zeros(Int, nv)
#     que[1] = r
#     qf = 0
#     qr = 1
#     visited[r] = true
#
#     # main loop
#     while qf < qr       # non-empty queue
#         u = que[qf += 1]
#         p = parents[u]
#         for (v, t) in zip(neighbors(g, u), neighbor_linktypes(g, u))
#             if v != p
#                 if visited[v]
#                     throw(ArgumentError("A cycle through $v is detected."))
#                 end
#                 # tag v as visited, and store its parent
#                 visited[v] = true
#                 parents[v] = u
#                 parents_lt[v] = -t
#                 # enqueue v as an active vertex
#                 que[qr+=1] = v
#             end
#         end
#     end
#
#     if qr < nv
#         throw(ArgumentError("The input graph has disjoint components."))
#     end
#     Tree(nv, nedgetypes(g), parents, parents_lt)
# end


text:
"""
Custom backtracking line-search.
Lnsrch returns a step-length α such that the decreasing criterion,
    f(x-αd) < f(x) - γ*|g|
is satisfied.
This is implemented based on other sophisticated line-search methods in LineSearches package.
"""

immutable Lnsrch{TF, TI}
    beta::TF
    gamma::TF
    iterations::TI
    func
end

(ls::Lnsrch)(df, x, s, x_scratch, lsr, alpha, mayterminate) =
    _lnsrch!(df, x, s, x_scratch, lsr, ls.func, alpha, mayterminate,
             ls.beta, ls.gamma, ls.iterations)

function _lnsrch!{T}(df,
                          x::Array{T},
                          s::Array{T},
                          x_scratch::Array{T},
                          lsr::LineSearchResults,
                          lsfunc,
                          alpha::Real = 1.0,
                          mayterminate::Bool = false,
                          beta::Real = 0.8,
                          gamma::Real = 0.01,
                          iterations::Integer = 1_000,
                          maxstep::Real = Inf)
    # Count the total number of iterations
    iteration = 0

    # Count number of parameters
    n = length(x)

    # read f_x and slope from LineSearchResults
    f_x = lsr.value[end]
    gxp = lsr.slope[end]

    # Tentatively move a distance of alpha in the direction of s
    x_scratch .= x .+ alpha.*s

    # Backtrack until we satisfy sufficient decrease condition
    f_x_scratch = lsfunc(x_scratch)
    while isinf(f_x_scratch) || isnan(f_x_scratch)
        alpha *= beta
        x_scratch .= x .+ alpha.*s
        f_x_scratch = lsfunc(x_scratch)
    end
    f_x_scratch = lsfunc(x_scratch)

    while f_x_scratch > f_x + gamma * gxp
        # Increment the number of steps we've had to perform
        iteration += 1

        # Ensure termination
        if iteration > iterations
            break
        end

        alpha *= beta

        # Update proposed position
        x_scratch .= x .+ alpha.*s
        # Evaluate f(x) at proposed position
        f_x_scratch = lsfunc(x_scratch)
        while f_x_scratch == Inf
            alpha *= beta
            x_scratch .= x .+ alpha.*s
            f_x_scratch = lsfunc(x_scratch)
        end
    end

    return alpha
end


text:
convert(::Type{Ptr{Nothing}}, obj::JavaObject{T}) where T = Ptr(obj)

convert(::Type{JString}, str::AbstractString) = JString(str)
convert(::Type{JObject}, str::AbstractString) = convert(JObject, JString(str))

# From JProxy
convert(::Type{AbstractString}, str::JString) = unsafe_string(str)
convert(::Type{JavaObject{Symbol("java.lang.Double")}}, n::Real) = jnew(Symbol("java.lang.Double"), (jdouble,), Float64(n))
convert(::Type{JavaObject{Symbol("java.lang.Float")}}, n::Real) = jnew(Symbol("java.lang.Float"), (jfloat,), Float32(n))
convert(::Type{JavaObject{Symbol("java.lang.Long")}}, n::Real) = jnew(Symbol("java.lang.Long"), (jlong,), Int64(n))
convert(::Type{JavaObject{Symbol("java.lang.Integer")}}, n::Real) = jnew(Symbol("java.lang.Integer"), (jint,), Int32(n))
convert(::Type{JavaObject{Symbol("java.lang.Short")}}, n::Real) = jnew(Symbol("java.lang.Short"), (jshort,), Int16(n))
convert(::Type{JavaObject{Symbol("java.lang.Byte")}}, n::Real) = jnew(Symbol("java.lang.Byte"), (jbyte,), Int8(n))
convert(::Type{JavaObject{Symbol("java.lang.Character")}}, n::Real) = jnew(Symbol("java.lang.Character"), (jchar,), Char(n))
convert(::Type{JavaObject{:int}}, n) = convert(jint, n)
convert(::Type{JavaObject{:long}}, n) = convert(jlong, n)
convert(::Type{JavaObject{:byte}}, n) = convert(jbyte, n)
convert(::Type{JavaObject{:boolean}}, n) = convert(jboolean, n)
convert(::Type{JavaObject{:char}}, n) = convert(jchar, n)
convert(::Type{JavaObject{:short}}, n) = convert(jshort, n)
convert(::Type{JavaObject{:float}}, n) = convert(jfloat, n)
convert(::Type{JavaObject{:double}}, n) = convert(jdouble, n)
convert(::Type{JavaObject{:void}}, n) = convert(jvoid, n)
convert(o::Type{JavaObject{T}}, ::Nothing) where T = o(J_NULL)

#Cast java object from S to T . Needed for polymorphic calls
function convert(::Type{JavaObject{T}}, obj::JavaObject{S}) where {T,S}
    if isConvertible(T, S)   #Safe static cast
        ptr = JNI.NewLocalRef(Ptr(obj))
        ptr === C_NULL && geterror()
        return JavaObject{T}(ptr)
    end
    isnull(obj) && throw(ArgumentError("Cannot convert NULL"))
    realClass = JNI.GetObjectClass(Ptr(obj))
    if isConvertible(T, realClass)  #dynamic cast
        ptr = JNI.NewLocalRef(Ptr(obj))
        ptr === C_NULL && geterror()
        return JavaObject{T}(ptr)
    end
    throw(JavaCallError("Cannot cast java object from $S to $T"))
end

#Is java type convertible from S to T.
isConvertible(T, S) = JNI.IsAssignableFrom(Ptr(metaclass(S)), Ptr(metaclass(T))) == JNI_TRUE
isConvertible(T, S::Ptr{Nothing} ) = JNI.IsAssignableFrom(S, Ptr(metaclass(T))) == JNI_TRUE

unsafe_convert(::Type{Ptr{Nothing}}, cls::JavaMetaClass) = Ptr(cls)
unsafe_convert(::Type{Ptr{Nothing}}, obj::JavaObject) = Ptr(obj)
unsafe_convert(::Type{Ptr{Nothing}}, ref::JavaRef) = Ptr(ref)

# Get the JNI/C type for a particular Java type
function real_jtype(rettype)
    if rettype <: JavaObject || rettype <: Array || rettype <: JavaMetaClass
        jnitype = Ptr{Nothing}
    else
        jnitype = rettype
    end
    return jnitype
end

function convert_args(argtypes::Tuple, args...)
    # Previously we converted to jvalue here
    # Instead convert to jvalue in _jcall so that we can delete local refs
    convertedArgs = Array{Any}(undef, length(args))
    savedArgs = Array{Any}(undef, length(args))
    for i in 1:length(args)
        r = convert_arg(argtypes[i], args[i])
        savedArgs[i] = r[1]
        convertedArgs[i] = r[2]
    end
    return savedArgs, convertedArgs
end

# We may be able to consolidate these convert_arg functions
function convert_arg(argtype::Type{JString}, arg)
    x = convert(JString, arg)
    return x,x
end

function convert_arg(argtype::Type, arg)
    x = convert(argtype, arg)
    return x,x
end
function convert_arg(argtype::Type{T}, arg) where T<:JavaObject
    x = convert(T, arg)::T
    return x,x
end

for (x, y, z) in [(:jboolean, :(JNI.NewBooleanArray), :(JNI.SetBooleanArrayRegion)),
                  (:jchar,    :(JNI.NewCharArray),    :(JNI.SetCharArrayRegion))   ,
                  (:jbyte,    :(JNI.NewByteArray),    :(JNI.SetByteArrayRegion))   ,
                  (:jshort,   :(JNI.NewShortArray),   :(JNI.SetShortArrayRegion))  ,
                  (:jint,     :(JNI.NewIntArray),     :(JNI.SetIntArrayRegion))    ,
                  (:jlong,    :(JNI.NewLongArray),    :(JNI.SetLongArrayRegion))   ,
                  (:jfloat,   :(JNI.NewFloatArray),   :(JNI.SetFloatArrayRegion))  ,
                  (:jdouble,  :(JNI.NewDoubleArray),  :(JNI.SetDoubleArrayRegion)) ]
    m = quote
        function convert_arg(argtype::Type{Array{$x,1}}, arg)
            carg = convert(argtype, arg)
            sz=length(carg)
            arrayptr = $y(sz)
            arrayptr === C_NULL && geterror()
            $z(arrayptr, 0, sz, carg)
            # Wrap into JavaObjects so we can delete the local refs
            return carg, JavaObject{argtype}(arrayptr)
        end
    end
    eval( m)
end

function convert_arg(argtype::Type{Array{T,1}}, arg) where T<:JavaObject
    carg = convert(argtype, arg)
    sz = length(carg)
    init = sz == 0 ? C_NULL : Ptr(carg[1])
    arrayptr = JNI.NewObjectArray(sz, Ptr(metaclass(T)), init)
    arrayptr === C_NULL && geterror()
    for i=2:sz
        JNI.SetObjectArrayElement(arrayptr, i-1, Ptr(carg[i]))
    end
    return carg, JavaObject{argtype}(arrayptr)
end

convert_result(rettype::Type{T}, result) where {T<:JString} = unsafe_string(JString(result))
convert_result(rettype::Type{T}, result) where {T<:JavaObject} = T(result)
convert_result(rettype, result) = result

for (x, y, z) in [(:jboolean, :(JNI.GetBooleanArrayElements), :(JNI.ReleaseBooleanArrayElements)),
                  (:jchar,    :(JNI.GetCharArrayElements),    :(JNI.ReleaseCharArrayElements))   ,
                  (:jbyte,    :(JNI.GetByteArrayElements),    :(JNI.ReleaseByteArrayElements))   ,
                  (:jshort,   :(JNI.GetShortArrayElements),   :(JNI.ReleaseShortArrayElements))  ,
                  (:jint,     :(JNI.GetIntArrayElements),     :(JNI.ReleaseIntArrayElements))    ,
                  (:jlong,    :(JNI.GetLongArrayElements),    :(JNI.ReleaseLongArrayElements))   ,
                  (:jfloat,   :(JNI.GetFloatArrayElements),   :(JNI.ReleaseFloatArrayElements))  ,
                  (:jdouble,  :(JNI.GetDoubleArrayElements),  :(JNI.ReleaseDoubleArrayElements)) ]
    m = quote
        function convert_result(rettype::Type{Array{$(x),1}}, result)
            sz = JNI.GetArrayLength(result)
            arr = $y(result, Ptr{jboolean}(C_NULL))
            jl_arr::Array = unsafe_wrap(Array, arr, Int(sz))
            jl_arr = deepcopy(jl_arr)
            $z(result, arr, Int32(0))
            return jl_arr
        end
    end
    eval(m)
end

function convert_result(rettype::Type{Array{JavaObject{T},1}}, result) where T
    sz = JNI.GetArrayLength(result)

    ret = Array{JavaObject{T}}(undef, sz)

    for i=1:sz
        a=JNI.GetObjectArrayElement(result, i-1)
        ret[i] = JavaObject{T}(a)
    end
    return ret
end

# covers return types like Vector{Vector{T}}
function convert_result(rettype::Type{Array{T,1}}, result) where T
    sz = JNI.GetArrayLength(result)

    ret = Array{T}(undef, sz)

    for i=1:sz
        a=JNI.GetObjectArrayElement(result, i-1)
        ret[i] = convert_result(T, a)
    end
    return ret
end

function convert_result(rettype::Type{Array{JavaObject{T},2}}, result) where T
    sz = JNI.GetArrayLength(result)
    if sz == 0
        return Array{T}(undef, 0,0)
    end
    a_1 = JNI.GetObjectArrayElement(result, 0)
    sz_1 = JNI.GetArrayLength(a_1)
    ret = Array{JavaObject{T}}(undef, sz, sz_1)
    for i=1:sz
        a = JNI.GetObjectArrayElement(result, i-1)
        # check that size of the current subarray is the same as for the first one
        sz_a = JNI.GetArrayLength(a)
        @assert(sz_a == sz_1, "Size of $(i)th subrarray is $sz_a, but size of the 1st subarray was $sz_1")
        for j=1:sz_1
            x = JNI.GetObjectArrayElement(a, j-1)
            ret[i, j] = JavaObject{T}(x)
        end
    end
    return ret
end

# matrices of primitive types and other arrays
function convert_result(rettype::Type{Array{T,2}}, result) where T
    sz = JNI.GetArrayLength(result)
    if sz == 0
        return Array{T}(undef, 0,0)
    end
    a_1 = JNI.GetObjectArrayElement(result, 0)
    sz_1 = JNI.GetArrayLength(a_1)
    ret = Array{T}(undef, sz, sz_1)
    for i=1:sz
        a = JNI.GetObjectArrayElement(result, i-1)
        # check that size of the current subarray is the same as for the first one
        sz_a = JNI.GetArrayLength(a)
        @assert(sz_a == sz_1, "Size of $(i)th subrarray is $sz_a, but size of the 1st subarray was $sz_1")
        ret[i, :] = convert_result(Vector{T}, a)
    end
    return ret
end

convert(::Type{jlong}, obj::JavaObject{Symbol("java.lang.Long")}) = jcall(obj, "longValue", jlong, ())
convert(::Type{jint}, obj::JavaObject{Symbol("java.lang.Integer")}) = jcall(obj, "intValue", jint, ())
convert(::Type{jdouble}, obj::JavaObject{Symbol("java.lang.Double")}) = jcall(obj, "doubleValue", jdouble, ())
convert(::Type{jfloat}, obj::JavaObject{Symbol("java.lang.Float")}) = jcall(obj, "floatValue", jfloat, ())
convert(::Type{jboolean}, obj::JavaObject{Symbol("java.lang.Boolean")}) = jcall(obj, "booleanValue", jboolean, ())

#The second term in this addition is due to the fact that Java converts all times to local time
function convert(::Type{DateTime}, x::@jimport(java.util.Date))
    if isnull(x)
        Dates.DateTime(1970,1,1,0,0,0)
    else
        Dates.unix2datetime(jcall(x, "getTime", jlong, ())/1000) +
            Second(round(div(Dates.value(now() - now(Dates.UTC)),1000)/900)*(900))
    end
end

function convert(::Type{DateTime}, x::JavaObject)
    isnull(x) && return Dates.DateTime(1970,1,1,0,0,0)
    JDate = @jimport(java.util.Date)
    if isConvertible(JDate, x)
        return convert(DateTime, convert(JDate, x))
    elseif isConvertible(@jimport(java.util.Calendar), x)
        return convert(DateTime, jcall(x, "getTime", JDate, ()))
    end
end

function convert(::Type{@jimport(java.util.Properties)}, x::Dict)
    Properties = @jimport(java.util.Properties)
    p = Properties(())
    for (n,v) in x
        jcall(p, "setProperty", @jimport(java.lang.Object), (JString, JString), n, v)
    end
    return p
end

function convert(::Type{@jimport(java.util.HashMap)}, K::Type{JavaObject{X}}, V::Type{JavaObject{Y}},
                 x::Dict) where {X,Y}
    Hashmap = @jimport(java.util.HashMap)
    p = Hashmap(())
    for (n,v) in x
        jcall(p, "put", @jimport(java.lang.Object), (JObject, JObject), n, v)
    end
    return p
end

function convert(::Type{@jimport(java.util.Map)}, K::Type{JavaObject{X}}, V::Type{JavaObject{Y}},
                 x::Dict) where {X,Y}
    convert(@jimport(java.util.Map), convert(@jimport(java.util.HashMap), K, V, x))
end

function convert(::Type{@jimport(java.util.ArrayList)}, x::Vector, V::Type{JavaObject{X}}=JObject) where X
    ArrayList = @jimport(java.util.ArrayList)
    a = ArrayList(())
    for v in x
        jcall(a, "add", jboolean, (JObject,), convert(V, v))
    end
    return a
end

function convert(::Type{@jimport(java.util.List)}, x::Vector, V::Type{JavaObject{X}}=JObject) where X
    convert(@jimport(java.util.ArrayList), x, V)
end

# Convert a reference to a java.lang.String into a Julia string. Copies the underlying byte buffer
unsafe_string(jstr::JString) = unsafe_string(Ptr(jstr))  #jstr must be a jstring obtained via a JNI call

function unsafe_string(jstr::Ptr{Nothing})  #jstr must be a jstring obtained via a JNI call
    if jstr == C_NULL; return ""; end #Return empty string to keep type stability. But this is questionable
    pIsCopy = Array{jboolean}(undef, 1)
    #buf::Ptr{UInt8} = JNI.GetStringUTFChars(Ptr(jstr), pIsCopy)
    buf = JNI.GetStringUTFChars(jstr, pIsCopy)
    s = unsafe_string(buf)
    JNI.ReleaseStringUTFChars(jstr, buf)
    return s
end

for (x, y, z) in [(:jboolean, :(JNI.GetBooleanArrayElements), :(JNI.ReleaseBooleanArrayElements)),
                  (:jchar,    :(JNI.GetCharArrayElements),    :(JNI.ReleaseCharArrayElements))   ,
                  (:jbyte,    :(JNI.GetByteArrayElements),    :(JNI.ReleaseByteArrayElements))   ,
                  (:jshort,   :(JNI.GetShortArrayElements),   :(JNI.ReleaseShortArrayElements))  ,
                  (:jint,     :(JNI.GetIntArrayElements),     :(JNI.ReleaseIntArrayElements))    ,
                  (:jlong,    :(JNI.GetLongArrayElements),    :(JNI.ReleaseLongArrayElements))   ,
                  (:jfloat,   :(JNI.GetFloatArrayElements),   :(JNI.ReleaseFloatArrayElements))  ,
                  (:jdouble,  :(JNI.GetDoubleArrayElements),  :(JNI.ReleaseDoubleArrayElements)) ]
    m = quote
        function convert(::Type{Array{$(x),1}}, obj::JObject)
            sz = JNI.GetArrayLength(Ptr(obj))
            arr = $y(Ptr(obj), Ptr{jboolean}(C_NULL))
            jl_arr::Array = unsafe_wrap(Array, arr, Int(sz))
            jl_arr = deepcopy(jl_arr)
            $z(Ptr(obj), arr, Int32(0))
            return jl_arr
        end
    end
    eval(m)
end

function convert(::Type{Array{T, 1}}, obj::JObject) where T
    sz = JNI.GetArrayLength(Ptr(obj))
    ret = Array{T}(undef, sz)
    for i=1:sz
        ptr = JNI.GetObjectArrayElement(Ptr(obj), i-1)
        ret[i] = convert(T, JObject(ptr))
    end
    return ret
end

##Iterator
iterator(obj::JavaObject) = jcall(obj, "iterator", @jimport(java.util.Iterator), ())

"""
Given a `JavaObject{T}` narrows down `T` to a real class of the underlying object.
For example, `JavaObject{:java.lang.Object}` pointing to `java.lang.String`
will be narrowed down to `JavaObject{:java.lang.String}`
"""
function narrow(obj::JavaObject)
    c = jcall(obj,"getClass", @jimport(java.lang.Class), ())
    t = jcall(c, "getName", JString, ())
    return convert(JavaObject{Symbol(t)}, obj)
end

has_next(itr::JavaObject) = (jcall(itr, "hasNext", jboolean, ()) == JNI_TRUE)

function Base.iterate(itr::JavaObject, state=nothing)
    if has_next(itr)
        o = jcall(itr, "next", @jimport(java.lang.Object), ())
        return (narrow(o), state)
    else
        return nothing
    end
end


text:
#===============================================================================

        Macros for generating a collection of helper functions that
        are required to perform conjugate updates in `DiffusionMCMC.jl`.
        If you don't use `DiffusionMCMC.jl`, then don't worry about
        this file, as the macros below and the autogenerated functions
        that they produce will most likely be of no use to you.

===============================================================================#

#=
            These functions are overwritten when calling the
            @conjugate_gaussian macro
                                                                              =#
@doc raw"""
    phi(::Val{T}, t, x, P) where T

If the drift can be written in a form:
```math
b_θ(t, x) = θ^Tφ(t,x)+ϕ(t,x)
```
for parameters $θ$ that are being updated, then `phi` is a row of the matrix
$φ(t,x)$ that get's multiplied by the coordinate of a θ vector with a name `T`.
"""
phi(::Val, t, x, P) = error("not implemented")

"""
    num_non_hypo(Ptype::Type{T}) where T <: DiffusionProcess

Return a total number of coordinates of the diffusion process that have
non-degenerate noise structure. I.e. it is equal to the total number of
coordinates of the process minus the number of coordinates that have no direct
Wiener contribution.
"""
num_non_hypo(Ptype) = error("not implemented")

@doc raw"""
    hypo_a_inv(t, x, P)

Similar to $a^{-1}:=(σσ^T)^{-1}$, with the only exception that all of the zero
rows of $σ$ are first removed to yield $\hat{σ}$, and then,
$(\hat{σ} \hat{σ}^T)^{-1}$ is computed.
"""
hypo_a_inv(t, x, P) = error("not implemented")

"""
    nonhypo(x, P)

Return the diffusion's coordinates that have direct contribution from some
non-degenerate Wiener terms, i.e. leave out coordinates whose Wiener term is
zero.
"""
nonhypo(x, P) = x

"""
    nonhypo_σ(t::Float64, x, P)

Return a sub-matrix of the full volatility matrix `σ` that consists of non-zero
rows of `σ`.
"""
nonhypo_σ(t::Float64, x, P) = σ(t, x, P)

"""
    ignore_for_cu(::Val{T}, P) where T

A flag to indicate that a parameter with name `T` has no contribution to the
computation of conjugate updates.
"""
ignore_for_cu(::Val, P) = false

#=
                Main macro for Conjugate Gaussian updates
                                                                              =#
@doc raw"""
    @conjugate_gaussian DIFFUSION_NAME begin
        BODY
    end

Define helper functions for the conjugate gaussian updates of diffusion
DIFFUSION_NAME. In BODY information can be passed as a list with elements of the
form:

    :entry-name --> entry-value

The following entry names have special meaning: `:inplace`, `:nonhypo`,
`:hypo_a_inv`. Any other name is interpreted as a parameter name for which
function `phi` is defined. More precisely, the following can be defined:

# `:inplace`

    :inplace --> true

to indicate that compmutations are to be done in-place. Otherwise, by default
they are assumed to be done out-of-place.

# `:nonhypo`

    :nonhypo --> list-of-indices-with-non-degenerate-noise

e.g.

    :nonhypo --> [1,3,5]

to indicate which coordinates have non-degenerate (i.e. non-zero) Wiener term
Wiener term. It automatically defines `num_non_hypo` for the user.

!!! warning
    If no flag `:inplace --> true` is passed, then it will be assumed that
    computations are done **out-of-place** by default and then no matter what
    subset you passed it will be internally saved as an `SVector` so as to work
    as a static way of accessing a subset of parameters. Otherwise, no change
    will be made to `list-of-indices-with-non-degenerate-noise`.

# `hypo_a_inv`

    :hypo_a_inv --> f(t,x,P)

where `f` is a function that uses time variable `t`, state variable `x` and a
struct with the target diffusion law `P`.

# `phi` and `ignore_for_cu`

    :name-of-parameter --> (φ₁(t,x,P), φ₂(t,x,P), ...)

where `φ₁(t,x,P)` corresponds to all terms from the drift that are in the
**first** non-degenerate coordinate **and** that are multiplied by the parameter
with the name `name-of-parameter`, etc. Any parameters of `P` for which `φ` is
not defined will have a function `ignore_for_cu` defined for them automatically
and returning `true`.
"""
macro conjugate_gaussian(diff_name, ex::Expr)
    parse_conjug_gsn(diff_name, MacroTools.striplines(ex))
end

function parse_conjug_gsn(diff_name, ex::Expr)
    fns = []
    declared_inplace = findfirst(line->line.args[1]==:inplace, ex.args)
    inplace = ( declared_inplace == nothing ? false : declared_inplace )
    #NOTE currently only out of place supported
    @assert !inplace
    pnames = Symbol[]

    for line in ex.args
        @assert line.head == :-->
        name, args = eval(line.args[1]), line.args[2]
        if name == :inplace
            continue
        elseif name == :nonhypo
            push!(fns, define_nonhypo(diff_name, Val(inplace), args))
            push!(fns, define_num_non_hypo(diff_name, length(eval(args))))
        elseif (name == :hypo_a_inv)
            push!(fns, define_a_inv(diff_name, args))
        else
            push!(fns, define_phi(diff_name, name, args))
            push!(pnames, name)
        end
    end

    additional_params = [
        p for p in parameter_names(eval(diff_name)) if !(p in pnames)
    ]
    for p in additional_params
        push!(fns, define_phi_omission(diff_name, p))
    end

    for fn in fns
        eval(fn)
    end

    #[println(f) for f in fns]
    nothing
end

function define_num_non_hypo(diff_name, args)
    fn_def = Expr(
        :call,
        :num_non_hypo,
        Expr(
            :(::),
            Expr(
                :curly,
                :Type,
                Expr(
                    :(<:),
                    diff_name,
                )
            )
        )
    )
    @assert typeof(args) <: Integer

    fn_body = Expr(
        :block,
        Expr(
            :call,
            :Val,
            args
        )
    )

    Expr(
        :macrocall,
        Symbol("@inline"),
        :(#= none:1 =#),
        Expr(:(=), fn_def, fn_body)
    )
end

function define_nonhypo(diff_name, ::Val{false}, args)
    coords = eval(args)
    N = length(coords)

    fn_def = Expr(
        :call,
        :nonhypo,
        :x,
        Expr(
            :(::),
            :P,
            diff_name,
        ),
    )

    fn_body = Expr(
        :block,
        Expr(
            :ref,
            :x,
            Expr(
                :call,
                Expr(
                    :curly,
                    :SVector,
                    N,
                    Int64
                ),
                args,
            )
        )
    )

    Expr(
        :macrocall,
        Symbol("@inline"),
        :(#= none:1 =#),
        Expr(:(=), fn_def, fn_body),
    )
    #Expr(:(=), fn_def, fn_body)
end

function define_a_inv(diff_name, args)
    fn_def = Expr(
        :call,
        :hypo_a_inv,
        :t,
        :x,
        Expr(
            :(::),
            :P,
            diff_name,
        ),
    )
    Expr(:(=), fn_def, args)
end

function define_phi(diff_name, name, args)
    fn_def = Expr(
        :call,
        :phi,
        Expr(
            :(::),
            Expr(
                :curly,
                :Val,
                QuoteNode(name),
            )
        ),
        :t,
        :x,
        Expr(
            :(::),
            :P,
            diff_name,
        ),
    )

    Expr(:(=), fn_def, args)
end

function define_phi_omission(diff_name, name)
    fn_def = Expr(
        :call,
        :ignore_for_cu,
        Expr(
            :(::),
            Expr(
                :curly,
                :Val,
                QuoteNode(name),
            )
        ),
        Expr(
            :(::),
            :P,
            diff_name,
        ),
    )
    Expr(
        :macrocall,
        Symbol("@inline"),
        :(#= none:1 =#),
        Expr(
            :(=),
            fn_def,
            Expr(
                :block,
                true
            )
        )
    )
end


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

struct PolicyDefinitionsApi <: SwaggerApi
    client::Swagger.Client
end

"""

Creates or updates a policy definition.
Param: policyDefinitionName::String (required)
Param: parameters::PolicyDefinition (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: PolicyDefinition
"""
function policyDefinitionsCreateOrUpdate(_api::PolicyDefinitionsApi, policyDefinitionName::String, parameters, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "PUT", PolicyDefinition, "/subscriptions/{subscriptionId}/providers/Microsoft.Authorization/policydefinitions/{policyDefinitionName}", ["azure_auth"], parameters)
    Swagger.set_param(_ctx.path, "policyDefinitionName", policyDefinitionName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json", "text/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json", "text/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

"""

Deletes a policy definition.
Param: policyDefinitionName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: Nothing
"""
function policyDefinitionsDelete(_api::PolicyDefinitionsApi, policyDefinitionName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "DELETE", Nothing, "/subscriptions/{subscriptionId}/providers/Microsoft.Authorization/policydefinitions/{policyDefinitionName}", ["azure_auth"])
    Swagger.set_param(_ctx.path, "policyDefinitionName", policyDefinitionName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json", "text/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json", "text/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

"""

Gets the policy definition.
Param: policyDefinitionName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: PolicyDefinition
"""
function policyDefinitionsGet(_api::PolicyDefinitionsApi, policyDefinitionName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", PolicyDefinition, "/subscriptions/{subscriptionId}/providers/Microsoft.Authorization/policydefinitions/{policyDefinitionName}", ["azure_auth"])
    Swagger.set_param(_ctx.path, "policyDefinitionName", policyDefinitionName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json", "text/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json", "text/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

"""

Gets all the policy definitions for a subscription.
Param: api_version::String (required)
Param: subscriptionId::String (required)
Param: filter::String
Return: PolicyDefinitionListResult
"""
function policyDefinitionsList(_api::PolicyDefinitionsApi, api_version::String, subscriptionId::String; filter=nothing, _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", PolicyDefinitionListResult, "/subscriptions/{subscriptionId}/providers/Microsoft.Authorization/policydefinitions", ["azure_auth"])
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "\$filter", filter)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json", "text/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json", "text/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

export policyDefinitionsCreateOrUpdate, policyDefinitionsDelete, policyDefinitionsGet, policyDefinitionsList


text:
function equals(obj::Object, arg0::Object)
    return jcall(obj, "equals", jboolean, (Object,), arg0)
end

function get_b(obj::BiDiagonalTransformer)
    return jcall(obj, "getB", RealMatrix, ())
end

function get_u(obj::BiDiagonalTransformer)
    return jcall(obj, "getU", RealMatrix, ())
end

function get_v(obj::BiDiagonalTransformer)
    return jcall(obj, "getV", RealMatrix, ())
end

function hash_code(obj::Object)
    return jcall(obj, "hashCode", jint, ())
end


text:
using MATLAB # To call spgl1

include("utils.jl")

##################################################
## Codebook update using SPGL1 and thresholding ##
##################################################

function update_codebooks_spgl1(
  X::Matrix{Float32}, # d-by-n matrix to update codebooks on
  B::Matrix{Int16},   # m-by-n matrix. X encoded
  h::Integer,         # number of entries per codebook
  tau::Float64,       # \tau argument in SPGL1. Sets sparsity.
  prevC::Vector{Matrix{Float32}}, # The previous set of codebooks, to use as init
  spgl1_path::AbstractString, # Path to spgl1
  V::Bool=false)      # whether to print progress

  if V print("Doing SPGL1 codebook update... "); st=time(); end

  d, n   = size(X);
  m, _   = size(B);
  C = sparsify_codes( B, h );
  K = SharedArray(Float32, d, size(C,2));

  prevK =  convert( Matrix{Float64}, hcat( prevC... ) );

  # Convert C and X for matlab-friendlines
  C = convert(SparseMatrixCSC{Float64,Int64}, C);
  X = convert( Matrix{Float64}, X );

  # Copy values to matlab
  @mput X
  @mput C
  @mput tau
  @mput prevK
  @mput spgl1_path

  mat"addpath( spgl1_path )"

  mat"options = spgSetParms()"
  mat"options.verbosity = 1"

  @time begin
  mat"
    Xt = X';
    prevKt = prevK';

    [d, n] = size( X );

    [n, mh] = size( C );
    A = @(x, mode) sparse_lsq_fun( C, x, n, d, mode );

    try
      [K, ~, ~, ~] = spgl1(A, Xt(:), tau, [], prevKt(:), options);
    catch
      warning('spgl1 failed for some reason. Redoing without init')
      [K, ~, ~, ~] = spgl1(C, Xt(:), tau, [], [], options);
    end

    K = reshape( K, size(C,2), d );
    K = K';
    "
  end

  @mget K

  K = convert(Matrix{Float32}, K);
  new_C = K2vec( K, m, h );

  if V @printf("done in %.2f seconds\n", time()-st); end

  return new_C;

end # function update_codebooks

function update_codebooks_spgl1_threshold(
  X::Matrix{Float32}, # d-by-n matrix to update codebooks on
  B::Matrix{Int16},   # m-by-n matrix. X encoded
  h::Integer,         # number of entries per codebook
  tau::Float64,       # \tau argument in SPGL1. Sets sparsity
  prevC::Vector{Matrix{Float32}}, # The previous set of codebooks, to use as init
  S::Integer, # The number of highest non-zero entries to keep
  spgl1_path::AbstractString, # Path to spgl1
  V::Bool=false)      # whether to print progress

  d, n   = size(X);
  m, _   = size(B);

  # Call the regular SPGL1 function
  C = update_codebooks_spgl1( X, B, h, tau, prevC, spgl1_path, V );
  K = hcat( C... );

  @printf("%d non-zero elements after update\n", sum( K .!= 0 ))

  # Keep the S max elements
  highest = sortperm( abs(K[:]), rev=true );

  nk = length( K );
  for i = S+1:nk
    K[ highest[i] ] = 0;
  end

  C = K2vec( K, m, h );

end

max_issues_repo_name: xiaosu-zhu/local-search-quantization

max_forks_repo_name: xiaosu-zhu/local-search-quantization


text:
const Maybe{T} = Union{T,Nothing}

const TupleN{T,N} = NTuple{N,T}

const AbsArr{T,N} = AbstractArray{T,N}
const AbsMat{T} = AbstractMatrix{T}
const AbsVec{T} = AbstractVector{T}

const AbsNestedArr{N} = AbstractArray{<:AbstractArray,N}
const AbsSimilarNestedArr{V,M,N} = AbstractArray{<:AbstractArray{V,M},N}

const StaticOrVal{X} = Union{Static{X}, Val{X}}


text:
#Kronecker Products
type KronProdMat{T} <: AbstractKronMat{T}
    terms::Vector{AbstractMatrix{T}}
end
type KronProdVec{T} <: AbstractKronVec{T}
    terms::Vector{AbstractVector{T}}
end

function KronProdMat(terms...)
    T = promote_type([eltype(term) for term in terms]...)
    return KronProdMat(AbstractMatrix{T}[terms...])
end
function KronProdVec(terms...)
    T = promote_type([eltype(term) for term in terms]...)
    return KronProdVec(AbstractVector{T}[terms...])
end

⊗(terms...) =  KronProdMat(terms...)

terms(C::KronProdMat) = C.terms
terms(C::KronProdVec) = C.terms

function getindex(C::KronProdMat, i::Integer,j::Integer)
    Ms = terms(C)
    hk,hl = kronindexes(C,i,j)
    return prod([M[k,l] for (M,k,l) in zip(Ms,hk,hl)])
end

function getindex(C::KronProdVec,i::Integer)
    Vs = terms(C)
    hk = kronindexes(C,i)
    return prod([V[k] for (V,k) in zip(Vs,hk)])
end

#unary operations
for f = (:ctranspose,:tranpose,:inv)
    @eval ($f)(C::KronProdMat) = KronProdMat(Any[($f)(term) for term in terms(C)]...)
end

for f = (:trace,:rank)
    @eval ($f)(C::KronProdMat) = prod([($f)(term) for term in terms(C)])
end

function det{T}(C::KronProdMat{T})
    Ms,Ns = [[S...] for S in zip(sizes(C)...)]
    bigM,bigN = prod(Ms),prod(Ns)
    bigM == bigN || throw(DimensionMismatch("matrix is not square"))
    if all(Ms .== Ns)
        return prod([det(term)^div(bigM,M) for (term,M) in zip(terms(C),Ms)])
    else
        return zero(T)
    end
end

^{T<:Integer}(C::KronProdMat,n::T) = KronProdMat(Any[term^n for term in terms(C)]...)

#nary ops

function *(Cs::KronProdMat...)
    return [*(factors...) for factors in zip(map(terms,Cs)...)]
end


text:
function parse_kw(ps::ParseState, ::Type{Val{Tokens.LET}})
    startbyte = ps.t.startbyte
    start_col = ps.t.startpos[2] + 4

    # Parsing
    ret = EXPR{Let}(EXPR[INSTANCE(ps)], -startbyte, Variable[], "")
    format_kw(ps)

    @default ps @closer ps comma @closer ps block while !closer(ps)
        @catcherror ps startbyte a = parse_expression(ps)
        push!(ret.args, a)
        if ps.nt.kind == Tokens.COMMA
            next(ps)
            push!(ret.args, INSTANCE(ps))
            format_comma(ps)
        end
    end
    block = EXPR{Block}(EXPR[], 0, Variable[], "")
    @catcherror ps startbyte @default ps parse_block(ps, block, start_col)

    # Construction
    push!(ret.args, block)
    next(ps)
    push!(ret.args, INSTANCE(ps))
    ret.span += ps.nt.startbyte

    return ret
end


text:
using ITensors
using ITensorInfiniteMPS
using Test
using Random

@testset "ITensorMap" begin
  i = Index(2)
  A1 = randomITensor(i'', dag(i))
  A2 = randomITensor(i', dag(i''))
  B1 = randomITensor(i'', dag(i))
  B2 = randomITensor(i', dag(i''))
  c1 = 2.3
  c2 = 3.4

  MA = ITensorMap([A1, A2])
  MB = ITensorMap([B1, B2])
  v = randomITensor(i)

  @test MA(v) ≈ noprime(A2 * A1 * v)
  @test MB(v) ≈ noprime(B2 * B1 * v)

  c1⨉MA = c1 * MA
  c2⨉MB = c2 * MB
  @test c1⨉MA isa ITensorMap
  @test c2⨉MB isa ITensorMap
  @test c1⨉MA(v) ≈ noprime(c1 * A2 * A1 * v)
  @test c2⨉MB(v) ≈ noprime(c2 * B2 * B1 * v)

  MA⁺MB = MA + MB
  @test MA⁺MB isa ITensorInfiniteMPS.ITensorMapSum
  @test MA⁺MB(v) ≈ noprime(A2 * A1 * v + B2 * B1 * v)

  c1MA⁺c2MB = c1 * MA + c2 * MB
  @test c1MA⁺c2MB isa ITensorInfiniteMPS.ITensorMapSum
  @test c1MA⁺c2MB(v) ≈ noprime(c1 * A2 * A1 * v + c2 * B2 * B1 * v)

  ⁻MA = -MA
  @test ⁻MA isa ITensorMap
  @test ⁻MA(v) ≈ noprime(-A2 * A1 * v)

  MA⁻MB = MA - MB
  @test MA⁻MB isa ITensorInfiniteMPS.ITensorMapSum
  @test MA⁻MB(v) ≈ noprime(A2 * A1 * v - B2 * B1 * v)
end


text:
module DynamicDiscreteModels

include("statisticalmodels-stopgap.jl")
include("statisticalmodels-sugar.jl")

#need to import StatsBase.StatisticalModel once I make the transition
import Optim
import Distributions: Dirichlet, wsample

export coef!, coef_jac!, rand, loglikelihood, loglikelihood_jac, dim, mle,
	DynamicDiscreteModel, estep, em, viterbi, hmm2ddm!,
	filtr

include("dynamicdiscretemodel.jl")
include("simulate.jl")
include("filtering.jl")
include("estep.jl")
include("emalgorithm.jl")
include("viterbi.jl")
include("loglikelihood.jl")

#non-core convenience function for front-end packages
include("hmm.jl")

end


text:
using Documenter, ImageFiltering

makedocs(
    modules  = [ImageFiltering, Kernel, KernelFactors, ImageFiltering.MapWindow],
    format   = Documenter.HTML(),
    sitename = "ImageFiltering",
    pages    = [
        "index.md",
        "Function reference" => "function_reference.md"
    ]
)

deploydocs(
    repo   = "github.com/JuliaImages/ImageFiltering.jl.git",
    push_preview = "true"
)


text:
########################################
## File Name: observation_test.jl
## Author: Haruki Nishimura (hnishimura@stanford.edu)
## Date Created: 2020/05/12
## Description: Test code for observation functions in
## src/state_transition_observation_models.jl
########################################

using LinearAlgebra
using Random

@testset "Observation Test" begin

# ObserveModel_Range
observeRange = ObserveModel_Range();
xPos = PhysPos(0.,zeros(2));
q_pos = [0.4,10.0];

o = observe(observeRange,xPos,q_pos)
@test isapprox(o, 10.0079968)

o2 = observe(observeRange,xPos,q_pos,Matrix(1.0I, 2, 2),MersenneTwister(123))

J = observe_jacobi_auto(observeRange,xPos,q_pos)
J2 = observe_jacobi(observeRange,xPos,q_pos)
@test all(isapprox.(J, J2))

# ObserveModel_Manipulate2D
observeMP = ObserveModel_Manipulate2D();
xMP = PhysManipulate2D(0.,[0.1,0.2],0.5,[0.3,0.4],0.6,3.,5.,[10.,11.],0.7);
uMP = MControl2D(0.,6.,7.,8.);

o3 = observe(observeMP,xMP,uMP)
@test all(isapprox.(o3, [3.60214
                        14.64766
                         0.5
                        -8.36860
                         2.50129
                         0.6
                       157.19801
                       -40.90407
                       -10.83419], atol=1e-5))

o4 = observe(observeMP, xMP, uMP, Matrix(1.0I, 9, 9), MersenneTwister(123))

J3 = observe_jacobi_auto(observeMP,xMP,uMP)
J4 = observe_jacobi(observeMP,xMP,uMP)
@test all(isapprox.(J3, J4))

end


text:
"""
gradient_forward(
  θ::AbstractVector{<:Real},
  vi::VarInfo,
  model::Function,
  spl::Union{Nothing, Sampler}=nothing,
  chunk_size::Int=CHUNKSIZE[],
)

Computes the gradient of the log joint of `θ` for the model specified by `(vi, spl, model)`
using forwards-mode AD from ForwardDiff.jl.
"""
function gradient_forward(
  θ::AbstractVector{<:Real},
  vi::VarInfo,
  model::Function,
  spl::Union{Nothing, Sampler}=nothing,
  chunk_size::Int=CHUNKSIZE[],
)
  # Record old parameters.
  θ_old = vi[spl]

  # Define function to compute log joint.
  function f(θ)
      vi[spl] = θ
      return -runmodel(model, vi, spl).logp
  end

  # Set chunk size and do ForwardMode.
  chunk = ForwardDiff.Chunk(min(length(θ), chunk_size))
  config = ForwardDiff.GradientConfig(f, θ, chunk)
  ∂l∂θ = ForwardDiff.gradient!(similar(θ), f, θ, config)

  # Replace old parameters to ensure this function doesn't mutate `vi`.
  vi[spl] = θ_old

  return getlogp(vi), ∂l∂θ
end

"""
gradient_reverse(
  θ::AbstractVector{<:Real},
  vi::VarInfo,
  model::Function,
  spl::Union{Nothing, Sampler}=nothing,
)

Computes the gradient of the log joint of `θ` for the model specified by `(vi, spl, model)`
using reverse-mode AD from Flux.jl.
"""
function gradient_reverse(
  θ::AbstractVector{<:Real},
  vi::Turing.VarInfo,
  model::Function,
  spl::Union{Nothing, Sampler}=nothing,
)
  # Specify objective function.
  function f(θ)
    vi[spl] = θ
    return -runmodel(model, vi, spl).logp
  end

  # Compute forward and reverse passes.
  l, ȳ = Tracker.forward(f, θ)
  ∂l∂θ = ȳ(1)[1]

  # Remove tracking info from variables in model (because mutable state).
  vi.logp = Tracker.data(vi.logp)
  vi[spl] .= Tracker.data.(vi[spl])

  # Return non-tracked gradient value
  return Tracker.data(l), Tracker.data(∂l∂θ)
end

function verifygrad(grad::AbstractVector{<:Real})
  if any(isnan.(grad)) || any(isinf.(grad))
    @warn("Numerical error has been found in gradients.")
    @warn("grad = $(grad)")
    return false
  else
    return true
  end
end


text:
module TestConstraintDet

using Test

using VecMathOptInterface
const MOI = VecMathOptInterface

function runtests()
    for name in names(@__MODULE__; all = true)
        if startswith("$(name)", "test_")
            @testset "$(name)" begin
                getfield(@__MODULE__, name)()
            end
        end
    end
    return
end

include("../utilities.jl")

function test_LogDet_basic()
    mock = MOI.Utilities.MockOptimizer(
        MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
    )
    bridged_mock = MOI.Bridges.Constraint.LogDet{Float64}(mock)
    MOI.Test.runtests(
        bridged_mock,
        MOI.Test.Config(),
        include = [
            "test_basic_VectorOfVariables_LogDetCone",
            "test_basic_VectorAffineFunction_LogDetCone",
            "test_basic_VectorQuadraticFunction_LogDetCone",
        ],
    )
    return
end

function _setup_LogDetConeTriangle(
    mock::MOI.Utilities.MockOptimizer,
    ::Type{F},
) where {F}
    var_primal = [0, 1, 0, 1, 1, 0, 1, 0, 0, 1]
    exp_duals = [[-1, -1, 1], [-1, -1, 1]]
    psd_dual = [1, 0, 1, -1, 0, 1, 0, -1, 0, 1]
    mock.optimize! =
        (mock::MOI.Utilities.MockOptimizer) -> MOI.Utilities.mock_optimize!(
            mock,
            var_primal,
            (MOI.VariableIndex, MOI.EqualTo{Float64}) => [2],
            (MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}) =>
                [-1],
            (MOI.VectorAffineFunction{Float64}, MOI.Nonnegatives) =>
                [[1, 1]],
            (F, MOI.ExponentialCone) => exp_duals,
            (
                MOI.VectorAffineFunction{Float64},
                MOI.PositiveSemidefiniteConeTriangle,
            ) => [psd_dual],
        )
    return
end

function test_conic_LogDetConeTriangle()
    mock = MOI.Utilities.MockOptimizer(
        MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
        eval_variable_constraint_dual = false,
    )
    bridged_mock = MOI.Bridges.Constraint.LogDet{Float64}(mock)
    config = MOI.Test.Config()
    _setup_LogDetConeTriangle(mock, MOI.VectorOfVariables)
    MOI.Test.test_conic_LogDetConeTriangle_VectorOfVariables(
        bridged_mock,
        config,
    )
    MOI.empty!(bridged_mock)
    _setup_LogDetConeTriangle(mock, MOI.VectorAffineFunction{Float64})
    MOI.Test.test_conic_LogDetConeTriangle_VectorAffineFunction(
        bridged_mock,
        config,
    )

    # set primal/dual start is not yet implemented for LogDet bridge
    ci = first(
        MOI.get(
            bridged_mock,
            MOI.ListOfConstraintIndices{
                MOI.VectorAffineFunction{Float64},
                MOI.LogDetConeTriangle,
            }(),
        ),
    )
    _test_delete_bridge(
        bridged_mock,
        ci,
        5,
        (
            (MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}, 0),
            (MOI.VectorAffineFunction{Float64}, MOI.ExponentialCone, 0),
            (
                MOI.VectorAffineFunction{Float64},
                MOI.PositiveSemidefiniteConeTriangle,
                0,
            ),
        ),
    )
    return
end

function test_conic_LogDetConeTriangle_2()
    mock = MOI.Utilities.MockOptimizer(
        MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
    )
    bridged_mock = MOI.Bridges.Constraint.LogDet{Float64}(mock)
    var_primal =
        [log(5), 3, 2, 2 / 3, 1, 1 / 3, 2.5, log(3), log(2 / 3), log(2.5)]
    exp_duals = [
        [-1, log(3) - 1, 1 / 3],
        [-1, log(2 / 3) - 1, 3 / 2],
        [-1, log(2.5) - 1, 0.4],
    ]
    psd_dual = [
        1,
        -1,
        1.6,
        0,
        -0.2,
        0.4,
        -1 / 3,
        0,
        0,
        1 / 3,
        1,
        -1.5,
        0,
        0,
        1.5,
        0,
        0.2,
        -0.4,
        0,
        0,
        0.4,
    ]
    mock.optimize! =
        (mock::MOI.Utilities.MockOptimizer) -> MOI.Utilities.mock_optimize!(
            mock,
            var_primal,
            (MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}) =>
                [-1.0],
            (MOI.VectorAffineFunction{Float64}, MOI.ExponentialCone) =>
                exp_duals,
            (
                MOI.VectorAffineFunction{Float64},
                MOI.PositiveSemidefiniteConeTriangle,
            ) => [psd_dual],
        )
    config = MOI.Test.Config()
    MOI.Test.test_conic_LogDetConeTriangle(bridged_mock, config)

    # set primal/dual start is not yet implemented for LogDet bridge
    ci = first(
        MOI.get(
            bridged_mock,
            MOI.ListOfConstraintIndices{
                MOI.VectorAffineFunction{Float64},
                MOI.LogDetConeTriangle,
            }(),
        ),
    )
    _test_delete_bridge(
        bridged_mock,
        ci,
        1,
        (
            (MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}, 0),
            (MOI.VectorAffineFunction{Float64}, MOI.ExponentialCone, 0),
            (
                MOI.VectorAffineFunction{Float64},
                MOI.PositiveSemidefiniteConeTriangle,
                0,
            ),
        ),
    )
    return
end

function test_conic_RootDetConeTriangle()
    mock = MOI.Utilities.MockOptimizer(
        MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
    )
    bridged_mock = MOI.Bridges.Constraint.RootDet{Float64}(mock)
    var_primal = [1, 1, 0, 1, 1, 0, 1]
    geomean_dual = [-1, 0.5, 0.5]
    psd_dual = [0.5, 0, 0.5, -0.5, 0, 0.5, 0, -0.5, 0, 0.5]
    mock.optimize! =
        (mock::MOI.Utilities.MockOptimizer) -> MOI.Utilities.mock_optimize!(
            mock,
            var_primal,
            (MOI.VectorAffineFunction{Float64}, MOI.Nonnegatives) =>
                [[0.5, 0.5]],
            (MOI.VectorAffineFunction{Float64}, MOI.GeometricMeanCone) =>
                [geomean_dual],
            (
                MOI.VectorAffineFunction{Float64},
                MOI.PositiveSemidefiniteConeTriangle,
            ) => [psd_dual],
        )
    config = MOI.Test.Config()
    MOI.Test.test_conic_RootDetConeTriangle_VectorOfVariables(
        bridged_mock,
        config,
    )
    MOI.empty!(bridged_mock)
    MOI.Test.test_conic_RootDetConeTriangle_VectorAffineFunction(
        bridged_mock,
        config,
    )

    # set primal/dual start is not yet implemented for RootDet bridge
    ci = first(
        MOI.get(
            bridged_mock,
            MOI.ListOfConstraintIndices{
                MOI.VectorAffineFunction{Float64},
                MOI.RootDetConeTriangle,
            }(),
        ),
    )
    _test_delete_bridge(
        bridged_mock,
        ci,
        4,
        (
            (MOI.VectorAffineFunction{Float64}, MOI.RotatedSecondOrderCone, 0),
            (MOI.VectorAffineFunction{Float64}, MOI.GeometricMeanCone, 0),
            (
                MOI.VectorAffineFunction{Float64},
                MOI.PositiveSemidefiniteConeTriangle,
                0,
            ),
        ),
    )
    return
end

function test_RootDet_basic()
    mock = MOI.Utilities.MockOptimizer(
        MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
    )
    bridged_mock = MOI.Bridges.Constraint.RootDet{Float64}(mock)
    MOI.Test.runtests(
        bridged_mock,
        MOI.Test.Config(),
        include = [
            "test_basic_VectorOfVariables_RootDetCone",
            "test_basic_VectorAffineFunction_RootDetCone",
            "test_basic_VectorQuadraticFunction_RootDetCone",
        ],
    )
    return
end

function test_conic_RootDetConeTriangle_2()
    mock = MOI.Utilities.MockOptimizer(
        MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
    )
    bridged_mock = MOI.Bridges.Constraint.RootDet{Float64}(mock)
    var_primal = [5^inv(3), 3, 2, 2 / 3, 1, 1 / 3, 2.5]
    multiplier = 5^inv(3) / 3
    geomean_dual = vcat(-1, multiplier * [inv(3), 1.5, 0.4])
    psd_dual =
        multiplier * [
            1,
            -1,
            1.6,
            0,
            -0.2,
            0.4,
            -1 / 3,
            0,
            0,
            1 / 3,
            1,
            -1.5,
            0,
            0,
            1.5,
            0,
            0.2,
            -0.4,
            0,
            0,
            0.4,
        ]
    mock.optimize! =
        (mock::MOI.Utilities.MockOptimizer) -> MOI.Utilities.mock_optimize!(
            mock,
            var_primal,
            (MOI.VectorAffineFunction{Float64}, MOI.GeometricMeanCone) =>
                [geomean_dual],
            (
                MOI.VectorAffineFunction{Float64},
                MOI.PositiveSemidefiniteConeTriangle,
            ) => [psd_dual],
        )
    MOI.Test.test_conic_RootDetConeTriangle(bridged_mock, MOI.Test.Config())

    # set primal/dual start is not yet implemented for RootDet bridge
    ci = first(
        MOI.get(
            bridged_mock,
            MOI.ListOfConstraintIndices{
                MOI.VectorAffineFunction{Float64},
                MOI.RootDetConeTriangle,
            }(),
        ),
    )
    _test_delete_bridge(
        bridged_mock,
        ci,
        1,
        (
            (MOI.VectorAffineFunction{Float64}, MOI.RotatedSecondOrderCone, 0),
            (MOI.VectorAffineFunction{Float64}, MOI.GeometricMeanCone, 0),
            (
                MOI.VectorAffineFunction{Float64},
                MOI.PositiveSemidefiniteConeTriangle,
                0,
            ),
        ),
    )
    return
end

end  # module

TestConstraintDet.runtests()