text:
import StructTypes: StructType, CustomStruct, DictType, ArrayType, StringType, NumberType, BoolType, NullType, NoStructType, Struct, OrderedStruct, Mutable, construct, AbstractType, subtypes, subtypekey

struct RawType <: StructType end

struct RawValue{S}
    bytes::S
    pos::Int
    len::Int
end

function rawbytes end

read(io::IO, ::Type{T}; kw...) where {T} = read(Base.read(io, String), T; kw...)
read(bytes::AbstractVector{UInt8}, ::Type{T}; kw...) where {T} = read(VectorString(bytes), T; kw...)

function _prepare_read(str::AbstractString, ::Type{T}) where {T}
    buf = codeunits(str)
    len = length(buf)
    if len == 0
        error = UnexpectedEOF
        pos = 0
        @goto invalid
    end
    pos = 1
    b = getbyte(buf, pos)
    @wh
    return buf, pos, len, b
@label invalid
    invalid(error, buf, pos, T)
end

function read(str::AbstractString, ::Type{T}; kw...) where {T}
    buf, pos, len, b = _prepare_read(str, T)
    pos, x = read(StructType(T), buf, pos, len, b, T; kw...)
    return x
end

"""
    JSON3.read!(json_str, x; kw...)

Incrementally update an instance of a mutable object `x` with the contents of `json_str`.  See [`JSON3.read`](@ref) for more details.
"""
function read!(str::AbstractString, x::T; kw...) where {T}
    buf, pos, len, b = _prepare_read(str, T)
    pos, x = read!(StructType(T), buf, pos, len, b, T, x; kw...)
    return x
end

read(::NoStructType, buf, pos, len, b, ::Type{T}; kw...) where {T} = throw(ArgumentError("$T doesn't have a defined `StructTypes.StructType`"))

# https://github.com/quinnj/JSON3.jl/issues/172
# treat Union{Bool, Real} as non-StructTypes.NumberType so parsing works as expected
read(::NumberType, buf, pos, len, b, S::Type{Union{Bool, T}}; kw...) where {T <: Real} =
    read(Struct(), buf, pos, len, b, S; kw...)

function read(::Struct, buf, pos, len, b, U::Union; kw...)
    # Julia implementation detail: Unions are sorted :)
    # This lets us avoid the below try-catch when U <: Union{Missing,T}
    if U.a === Nothing || U.a === Missing
        if buf[pos] == UInt8('n')
            return read(StructType(U.a), buf, pos, len, b, U.a)
        else
            return read(StructType(U.b), buf, pos, len, b, U.b; kw...)
        end
    end
    try
        return read(StructType(U.a), buf, pos, len, b, U.a; kw...)
    catch e
        return read(StructType(U.b), buf, pos, len, b, U.b; kw...)
    end
end

const GLOBAL_IGNORED_TAPE = zeros(UInt64, 1024)

@inline function read(::RawType, buf, pos, len, b, ::Type{T}; kw...) where {T}
    newpos, _ = read!(buf, pos, len, b, GLOBAL_IGNORED_TAPE, 1, Any; kw...)
    return newpos, construct(T, RawValue(buf, pos, newpos - pos))
end

@inline function read(::Struct, buf, pos, len, b, ::Type{Any}; allow_inf::Bool=false, kw...)
    if b == UInt8('{')
        return read(DictType(), buf, pos, len, b, Dict{String, Any}; allow_inf=allow_inf, kw...)
    elseif b == UInt8('[')
        return read(ArrayType(), buf, pos, len, b, Base.Array{Any}; allow_inf=allow_inf, kw...)
    elseif b == UInt8('"')
        return read(StringType(), buf, pos, len, b, String; kw...)
    elseif b == UInt8('n')
        return read(NullType(), buf, pos, len, b, Nothing; kw...)
    elseif b == UInt8('t')
        return read(BoolType(), buf, pos, len, b, Bool; kw...)
    elseif b == UInt8('f')
        return read(BoolType(), buf, pos, len, b, Bool; kw...)
    elseif (UInt8('0') <= b <= UInt8('9')) || b == UInt8('-') || b == UInt8('+') || (allow_inf && (b == UInt8('N') || b == UInt8('I')))
        float, code, pos = Parsers.typeparser(Float64, buf, pos, len, b, Int16(0), Parsers.OPTIONS)
        if code > 0
            int = unsafe_trunc(Int64, float)
            if int == float
                return pos, int
            else
                return pos, float
            end
        end
    end
@label invalid
    invalid(InvalidChar, buf, pos, Any)
end

function read(::StringType, buf, pos, len, b, ::Type{T}; kw...) where {T}
    if b != UInt8('"')
        error = ExpectedOpeningQuoteChar
        @goto invalid
    end
    pos += 1
    @eof
    strpos = pos
    strlen = 0
    escaped = false
    b = getbyte(buf, pos)
    while b != UInt8('"')
        if b == UInt8('\\')
            escaped = true
            # skip next character
            pos += 2
            strlen += 2
        else
            pos += 1
            strlen += 1
        end
        @eof
        b = getbyte(buf, pos)
    end
    ptr = pointer(buf, strpos)
    return pos + 1, escaped ? construct(T, unescape(PointerString(ptr, strlen)); kw...) : construct(T, ptr, strlen; kw...)

@label invalid
    invalid(error, buf, pos, T)
end

function read(::BoolType, buf, pos, len, b, ::Type{T}; kw...) where {T}
    if pos + 3 <= len &&
        b            == UInt8('t') &&
        buf[pos + 1] == UInt8('r') &&
        buf[pos + 2] == UInt8('u') &&
        buf[pos + 3] == UInt8('e')
        return pos + 4, construct(T, true; kw...)
    elseif pos + 4 <= len &&
        b            == UInt8('f') &&
        buf[pos + 1] == UInt8('a') &&
        buf[pos + 2] == UInt8('l') &&
        buf[pos + 3] == UInt8('s') &&
        buf[pos + 4] == UInt8('e')
        return pos + 5, construct(T, false; kw...)
    else
        invalid(InvalidChar, buf, pos, Bool)
    end
end

function read(::NullType, buf, pos, len, b, ::Type{T}; kw...) where {T}
    if pos + 3 <= len &&
        b            == UInt8('n') &&
        buf[pos + 1] == UInt8('u') &&
        buf[pos + 2] == UInt8('l') &&
        buf[pos + 3] == UInt8('l')
        return pos + 4, construct(T, nothing; kw...)
    else
        invalid(InvalidChar, buf, pos, T)
    end
end

function read(::NumberType, buf, pos, len, b, ::Type{T}; parsequoted::Bool=false, kw...) where {T}
    quoted = false
    if parsequoted && b == UInt8('"')
        pos += 1
        @eof
        b = getbyte(buf, pos)
        @wh
        quoted = true
    end
    x, code, pos = Parsers.typeparser(StructTypes.numbertype(T), buf, pos, len, b, Int16(0), Parsers.OPTIONS)
    if quoted
        b = getbyte(buf, pos)
        @assert b == UInt8('"')
        pos += 1
    end
    if code > 0
        return pos, construct(T, x; kw...)
    end
@label invalid
    invalid(InvalidChar, buf, pos, T)
end

@inline read(::ArrayType, buf, pos, len, b, ::Type{T}; kw...) where {T} = read(ArrayType(), buf, pos, len, b, T, Base.IteratorEltype(T) == Base.HasEltype() ? eltype(T) : Any; kw...)
@inline read(::ArrayType, buf, pos, len, b, ::Type{T}, ::Type{eT}; kw...) where {T, eT} = readarray(buf, pos, len, b, T, eT; kw...)
read(::ArrayType, buf, pos, len, b, ::Type{Tuple}, ::Type{eT}; kw...) where {eT} = readarray(buf, pos, len, b, Tuple, eT; kw...)

@inline function readarray(buf, pos, len, b, ::Type{T}, ::Type{eT}; kw...) where {T, eT}
    if b != UInt8('[')
        error = ExpectedOpeningArrayChar
        @goto invalid
    end
    pos += 1
    @eof
    b = getbyte(buf, pos)
    @wh
    vals = Vector{eT}(undef, 0)
    if b == UInt8(']')
        return pos + 1, construct(T, vals; kw...)
    end
    while true
        # positioned at start of value
        pos, y = read(StructType(eT), buf, pos, len, b, eT; kw...)
        push!(vals, y)
        @eof
        b = getbyte(buf, pos)
        @wh
        if b == UInt8(']')
            return pos + 1, construct(T, vals; kw...)
        elseif b != UInt8(',')
            error = ExpectedComma
            @goto invalid
        end
        pos += 1
        @eof
        b = getbyte(buf, pos)
        @wh
    end

@label invalid
    invalid(error, buf, pos, T)
end

mutable struct TupleClosure{T, KW}
    buf::T
    pos::Int
    len::Int
    b::UInt8
    kw::KW
end

@inline function (f::TupleClosure)(i, nm, TT)
    buf, pos, len, b = f.buf, f.pos, f.len, f.b
    pos, x = read(StructType(TT), buf, pos, len, b, TT; f.kw...)
    @eof
    b = getbyte(buf, pos)
    @wh
    if b == UInt8(']')
        f.pos = pos + 1
        return x
    elseif b == UInt8(',')
        pos += 1
        @eof
        b = getbyte(buf, pos)
        @wh
        f.pos = pos
        f.b = b
        return x
    else
        error = ExpectedComma
        @goto invalid
    end
@label invalid
    invalid(error, buf, pos, TT)
end

@inline function read(::ArrayType, buf, pos, len, b, ::Type{T}, ::Type{eT}; kw...) where {T <: Tuple, eT}
    if b != UInt8('[')
        error = ExpectedOpeningArrayChar
        @goto invalid
    end
    pos += 1
    @eof
    b = getbyte(buf, pos)
    @wh
    if b == UInt8(']')
        pos += 1
        return pos, T()
    end
    c = TupleClosure(buf, pos, len, b, values(kw))
    x = StructTypes.construct(c, T)

    return c.pos, x
@label invalid
    invalid(error, buf, pos, T)
end

keyvalue(::Type{Symbol}, escaped, ptr, len) = escaped ? Symbol(unescape(PointerString(ptr, len))) : _symbol(ptr, len)
keyvalue(::Type{T}, escaped, ptr, len) where {T} = escaped ? construct(T, unescape(PointerString(ptr, len))) : construct(T, unsafe_string(ptr, len))

@inline read(::DictType, buf, pos, len, b, ::Type{T}; kw...) where {T} = read(DictType(), buf, pos, len, b, T, Symbol, Any; kw...)
@inline read(::Struct, buf, pos, len, b, ::Type{NamedTuple}; kw...) = read(DictType(), buf, pos, len, b, NamedTuple, Symbol, Any; kw...)
@inline read(::Struct, buf, pos, len, b, ::Type{NamedTuple{names}}; kw...) where {names} = read(DictType(), buf, pos, len, b, NamedTuple{names}, Symbol, Any; kw...)
@inline read(::DictType, buf, pos, len, b, ::Type{Dict}; kw...) = read(DictType(), buf, pos, len, b, Dict, String, Any; kw...)
@inline read(::DictType, buf, pos, len, b, ::Type{T}; kw...) where {T <: AbstractDict} = read(DictType(), buf, pos, len, b, T, keytype(T), valtype(T); kw...)

@inline function read(::DictType, buf, pos, len, b, ::Type{T}, ::Type{K}, ::Type{V}; kw...) where {T, K, V}
    if b != UInt8('{')
        error = ExpectedOpeningObjectChar
        @goto invalid
    end
    pos += 1
    @eof
    b = getbyte(buf, pos)
    @wh
    x = Dict{K, V}()
    if b == UInt8('}')
        return pos + 1, construct(T, x; kw...)
    elseif b != UInt8('"')
        error = ExpectedOpeningQuoteChar
        @goto invalid
    end
    pos += 1
    @eof
    while true
        keypos = pos
        keylen = 0
        escaped = false
        # read first key character
        b = getbyte(buf, pos)
        # positioned at first character of object key
        while b != UInt8('"')
            if b == UInt8('\\')
                escaped = true
                # skip next character
                pos += 2
                keylen += 2
            else
                pos += 1
                keylen += 1
            end
            @eof
            b = getbyte(buf, pos)
        end
        key = keyvalue(K, escaped, pointer(buf, keypos), keylen)
        pos += 1
        @eof
        b = getbyte(buf, pos)
        @wh
        if b != UInt8(':')
            error = ExpectedSemiColon
            @goto invalid
        end
        pos += 1
        @eof
        b = getbyte(buf, pos)
        @wh
        # now positioned at start of value
        pos, y = read(StructType(V), buf, pos, len, b, V; kw...)
        x[key] = y
        @eof
        b = getbyte(buf, pos)
        @wh
        if b == UInt8('}')
            return pos + 1, construct(T, x; kw...)
        elseif b != UInt8(',')
            error = ExpectedComma
            @goto invalid
        end
        pos += 1
        @eof
        b = getbyte(buf, pos)
        @wh
        if b != UInt8('"')
            error = ExpectedOpeningQuoteChar
            @goto invalid
        end
        pos += 1
        @eof
    end

@label invalid
    invalid(error, buf, pos, Object)
end

@inline function read(::CustomStruct, buf, pos, len, b, ::Type{T}; kw...) where {T}
    S = StructTypes.lowertype(T)
    pos, x = read(StructType(S), buf, pos, len, b, S; kw...)
    return pos, StructTypes.construct(T, x)
end

mutable struct MutableClosure{T, KW}
    buf::T
    pos::Int
    len::Int
    b::UInt8
    kw::KW
end

@inline function (f::MutableClosure)(i, nm, TT; kw...)
    kw2 = merge(values(kw), f.kw)
    pos_i, y_i = read(StructType(TT), f.buf, f.pos, f.len, f.b, TT; kw2...)
    f.pos = pos_i
    return y_i
end

@inline function read(::Mutable, buf, pos, len, b, ::Type{T}; kw...) where {T}
    x = T()
    pos, x = read!(Mutable(), buf, pos, len, b, T, x; kw...)
    return pos, x
end

@inline function read!(::Any, buf, pos, len, b, ::Type{T}, x::T; kw...) where {T}
    throw(ArgumentError("read! is only defined when T is of the `Mutable` struct type"))
end

@inline function read!(::Mutable, buf, pos, len, b, ::Type{T}, x::T; kw...) where {T}
    if b != UInt8('{')
        error = ExpectedOpeningObjectChar
        @goto invalid
    end
    pos += 1
    @eof
    b = getbyte(buf, pos)
    @wh
    if b == UInt8('}')
        pos += 1
        return pos, x
    elseif b != UInt8('"')
        error = ExpectedOpeningQuoteChar
        @goto invalid
    end
    pos += 1
    @eof
    while true
        keypos = pos
        keylen = 0
        escaped = false
        b = getbyte(buf, pos)
        while b != UInt8('"')
            if b == UInt8('\\')
                escaped = true
                # skip next character
                pos += 2
                keylen += 2
            else
                pos += 1
                keylen += 1
            end
            @eof
            b = getbyte(buf, pos)
        end
        key = keyvalue(Symbol, escaped, pointer(buf, keypos), keylen)
        pos += 1
        @eof
        b = getbyte(buf, pos)
        @wh
        if b != UInt8(':')
            error = ExpectedSemiColon
            @goto invalid
        end
        pos += 1
        @eof
        b = getbyte(buf, pos)
        @wh
        c = MutableClosure(buf, pos, len, b, values(kw))
        if StructTypes.applyfield!(c, x, key)
            pos = c.pos
        else
            pos, _ = read(Struct(), buf, pos, len, b, Any)
        end
        @eof
        b = getbyte(buf, pos)
        @wh
        if b == UInt8('}')
            pos += 1
            return pos, x
        elseif b != UInt8(',')
            error = ExpectedComma
            @goto invalid
        end
        pos += 1
        @eof
        b = getbyte(buf, pos)
        @wh
        if b != UInt8('"')
            error = ExpectedOpeningQuoteChar
            @goto invalid
        end
        pos += 1
        @eof
    end

@label invalid
    invalid(error, buf, pos, T)
end

mutable struct StructClosure{T, KW}
    buf::T
    pos::Int
    len::Int
    b::UInt8
    values::Vector{Any}
    kw::KW
end

const DEFAULT_STRUCT_FIELD_COUNT = 32

@inline function (f::StructClosure)(i, nm, TT; kw...)
    kw2 = merge(values(kw), f.kw)
    pos_i, y_i = read(StructType(TT), f.buf, f.pos, f.len, f.b, TT; kw2...)
    f.pos = pos_i
    if i <= DEFAULT_STRUCT_FIELD_COUNT
        f.values[i] = y_i
    else
        push!(f.values, y_i)
    end
    return
end

@inline function read(::Struct, buf, pos, len, b, ::Type{T}; kw...) where {T}
    values = Vector{Any}(undef, DEFAULT_STRUCT_FIELD_COUNT)
    if b != UInt8('{')
        error = ExpectedOpeningObjectChar
        @goto invalid
    end
    pos += 1
    @eof
    b = getbyte(buf, pos)
    @wh
    if b == UInt8('}')
        pos += 1
        @goto done
    elseif b != UInt8('"')
        error = ExpectedOpeningQuoteChar
        @goto invalid
    end
    pos += 1
    @eof
    while true
        keypos = pos
        keylen = 0
        escaped = false
        b = getbyte(buf, pos)
        while b != UInt8('"')
            if b == UInt8('\\')
                escaped = true
                # skip next character
                pos += 2
                keylen += 2
            else
                pos += 1
                keylen += 1
            end
            @eof
            b = getbyte(buf, pos)
        end
        key = keyvalue(Symbol, escaped, pointer(buf, keypos), keylen)
        pos += 1
        @eof
        b = getbyte(buf, pos)
        @wh
        if b != UInt8(':')
            error = ExpectedSemiColon
            @goto invalid
        end
        pos += 1
        @eof
        b = getbyte(buf, pos)
        @wh
        c = StructClosure(buf, pos, len, b, values, Base.values(kw))
        if StructTypes.applyfield(c, T, key)
            pos = c.pos
        else
            pos, _ = read(Struct(), buf, pos, len, b, Any)
        end
        @eof
        b = getbyte(buf, pos)
        @wh
        if b == UInt8('}')
            pos += 1
            @goto done
        elseif b != UInt8(',')
            error = ExpectedComma
            @goto invalid
        end
        pos += 1
        @eof
        b = getbyte(buf, pos)
        @wh
        if b != UInt8('"')
            error = ExpectedOpeningQuoteChar
            @goto invalid
        end
        pos += 1
        @eof
    end

@label done
    return pos, StructTypes.construct(values, T)

@label invalid
    invalid(error, buf, pos, T)
end

mutable struct OrderedStructClosure{T, KW}
    buf::T
    pos::Int
    len::Int
    kw::KW
end

@inline function (f::OrderedStructClosure)(i, nm, TT)
    pos_i, x_i = readvalue(f.buf, f.pos, f.len, TT; f.kw...)
    f.pos = pos_i
    return x_i
end

@inline function read(::OrderedStruct, buf, pos, len, b, ::Type{T}; kw...) where {T}
    if b != UInt8('{')
        error = ExpectedOpeningObjectChar
        @goto invalid
    end
    pos += 1
    @eof
    b = getbyte(buf, pos)
    @wh
    if b == UInt8('}')
        pos += 1
        @goto done
    elseif b != UInt8('"')
        error = ExpectedOpeningQuoteChar
        @goto invalid
    end
    pos += 1
    @eof
    c = OrderedStructClosure(buf, pos, len, values(kw))
    x = StructTypes.construct(c, T)
    return c.pos, x

@label done
    return pos, StructTypes.construct(values, T)

@label invalid
    invalid(error, buf, pos, T)
end

@inline function readvalue(buf, pos, len, ::Type{T}; kw...) where {T}
    b = getbyte(buf, pos)
    while b != UInt8('"')
        pos += ifelse(b == UInt8('\\'), 2, 1)
        @eof
        b = getbyte(buf, pos)
    end
    pos += 1
    @eof
    b = getbyte(buf, pos)
    @wh
    if b != UInt8(':')
        error = ExpectedSemiColon
        @goto invalid
    end
    pos += 1
    @eof
    b = getbyte(buf, pos)
    @wh
    # read value
    pos, y = read(StructType(T), buf, pos, len, b, T; kw...)
    @eof
    b = getbyte(buf, pos)
    @wh
    if b == UInt8('}')
        pos += 1
        return pos, y
    elseif b != UInt8(',')
        error = ExpectedComma
        @goto invalid
    end
    pos += 1
    @eof
    b = getbyte(buf, pos)
    @wh
    if b != UInt8('"')
        error = ExpectedOpeningQuoteChar
        @goto invalid
    end
    pos += 1
    @eof
    return pos, y
@label invalid
    invalid(error, buf, pos, T)
end

@inline function read(::AbstractType, buf, pos, len, b, ::Type{T}; kw...) where {T}
    types = subtypes(T)
    if length(types) == 1
        only_subtype = types[1]
        return read(StructType(only_subtype), buf, pos, len, b, only_subtype; kw...)
    end
    return _read(AbstractType(), buf, pos, len, b, T; kw...)
end

@inline function _read(::AbstractType, buf, pos, len, b, ::Type{T}; kw...) where {T}
    startpos = pos
    startb = b
    if b != UInt8('{')
        error = ExpectedOpeningObjectChar
        @goto invalid
    end
    pos += 1
    @eof
    b = getbyte(buf, pos)
    @wh
    if b == UInt8('}')
        throw(ArgumentError("invalid json abstract type"))
    elseif b != UInt8('"')
        error = ExpectedOpeningQuoteChar
        @goto invalid
    end
    pos += 1
    @eof
    types = subtypes(T)
    skey = subtypekey(T)
    while true
        keypos = pos
        keylen = 0
        escaped = false
        b = getbyte(buf, pos)
        while b != UInt8('"')
            if b == UInt8('\\')
                escaped = true
                # skip next character
                pos += 2
                keylen += 2
            else
                pos += 1
                keylen += 1
            end
            @eof
            b = getbyte(buf, pos)
        end
        key = keyvalue(Symbol, escaped, pointer(buf, keypos), keylen)
        pos += 1
        @eof
        b = getbyte(buf, pos)
        @wh
        if b != UInt8(':')
            error = ExpectedSemiColon
            @goto invalid
        end
        pos += 1
        @eof
        b = getbyte(buf, pos)
        @wh
        # read value
        pos, val = read(Struct(), buf, pos, len, b, Any)
        if key == skey
            TT = types[Symbol(val)]
            return read(StructType(TT), buf, startpos, len, startb, TT; kw...)
        end
        @eof
        b = getbyte(buf, pos)
        @wh
        if b == UInt8('}')
            pos += 1
            throw(ArgumentError("invalid json abstract type: didn't find subtypekey"))
        elseif b != UInt8(',')
            error = ExpectedComma
            @goto invalid
        end
        pos += 1
        @eof
        b = getbyte(buf, pos)
        @wh
        if b != UInt8('"')
            error = ExpectedOpeningQuoteChar
            @goto invalid
        end
        pos += 1
        @eof
    end

@label invalid
    invalid(error, buf, pos, T)
end


text:
export normsum, norm01, normeuclid, normmean, normmeanstd
export normsum!, norm01!, normeuclid!, normmean!, normmeanstd!
export normquantile
export normunique, valuemap, unvaluemap
export normmean_
export pcawhitening, zcawhitening, WhiteningBase
export clamp, clamp!
export nanfunction, nanmean, nanstd, nanmedian
export distance
export randbase, randnormbase, randproj
export mean_, std_, var_, min_, max_, sum_, median_, dropdims_

#######################################
##  normsum, norm01, normeuclid, normmeanstd
##  normsum!, norm01!, normeuclid!, normmeanstd!

copyfloat(a::AbstractArray{T}) where {T<:Int} = float(a)
copyfloat(a::AbstractArray{T}) where {T} = copy(a)

normsum(a) = normsum!(copyfloat(a))
function normsum!(r)
    a = sum(r)
    a == 0 && return r
    for i = 1:length(r)
        r[i] /= a
    end
    r
end

norm01(a) = norm01!(copyfloat(a))
function norm01!(r)
    mi = minimum(r)
    ma = maximum(r)
    d = ma - mi
    for i = 1:length(r)
        r[i] -= mi
    end
    if d>0
        for i = 1:length(r)
            r[i] /= d
        end
    end
    r
end

normeuclid(a) = normeuclid!(copyfloat(a))
function normeuclid!(r)
    a = zeroel(r)
    for i = 1:length(r)
        a += r[i]*r[i]
    end
    a = sqrt(a)

    if a > 0
        for i = 1:length(r)
            r[i] /= a
        end
    end
    r
end

normmean(a) = normmean!(copyfloat(a))
function normmean!(r)
    a = mean(r)
    for i = 1:length(r)
        r[i] /= a
    end
    r
end

normmeanstd(a) = normmeanstd!(copyfloat(a))
function normmeanstd!(r)
    m = mean(r)
    s = std(r)
    for i = 1:length(r)
        r[i] -= m
    end
    if s > 0
        for i = 1:length(r)
            r[i] /= s
        end
    end
    r
end

normquantile(a, q1, q2) = normquantile(a, [q1,q2])
function normquantile(a, q = [0.1, 0.9])
    qs = quantile(vec(a), q)
    clamp((a .- qs[1]) ./ (qs[2]-qs[1]) .* (q[2]-q[1]) .+ q[1], 0, 1)
end

function normunique(a)  # FIXME
    u = unique(a)
    d = Dict()
    [d[u[i]] = i for i = 1:length(u)]
    d
    r = zero(a)
    for i = 1:length(a)
        r[i] = d[a[i]]
    end
    r
end

# #######################################
# ##  minimum and maximum with optional functions and number of values to find

import Base.maximum
maximum_ = maximum
maximum(f::Function, g::Function) = error("not defined")
maximum(x::Array,f::Function) = @p map x f | argmax | at x _

import Base.minimum
minimum_ = minimum
minimum(f::Function, g::Function) = error("not defined")
minimum(x::Array,f::Function) = @p map x f | argmin | at x _

function valuemap(a, m::Dict, default = 0)
    map(x->get(m,x,default),a)
end
function unvaluemap(a, m::Dict, default = 0)
    valuemap(a, flip(m), default)
end

function valuemap(data, mapping)
    r = Array{promote_type(eltype(mapping),(eltype(data))), ndims(data)}(undef,size(data))
    for i = 1:length(data)
        v = data[i]
        if !isa(v,Number) || (!isnan(v) && v>0)
            if isa(v,Number) && v>length(mapping)
                error("valuemap: v>length(mapping), with v==$v and length(mapping)==$(length(mapping))")
            end
            r[i] = mapping[v]
        else
            r[i] = data[i]
        end
    end
    r
end

import Base.clamp, Base.clamp!
clamp(a::AbstractArray, mi::Number, ma::Number) = clamp.(a, mi, ma)
clamp(a::AbstractArray{T}, mi::Union{AbstractArray,Tuple}, ma::Union{AbstractArray,Tuple}) where {T} = (r = Base.copy(a); clamp!(r, r, mi, ma); r)
clamp!(a::AbstractArray{T,2}, mi::Union{AbstractArray,Tuple}, ma::Union{AbstractArray,Tuple}) where {T} = clamp!(a, a, mi, ma)
function clamp!(r::AbstractArray{T,2}, a::AbstractArray{T,2}, mi::Union{AbstractArray,Tuple}, ma::Union{AbstractArray,Tuple}) where {T}
    if !(size(a,1)==length(mi)==length(ma))
        error("clamp!: size(a,1)==length(mi)==length(ma) was false: $(size(a,1))==$(length(mi))==$(length(ma))")
    end
    for j = 1:size(a,2), i = 1:size(a,1)
        r[i,j] = min(max(mi[i],a[i,j]), ma[i])
    end
    r
end
clamp(a::Tuple, mi::Union{AbstractArray,Tuple}, ma::Union{AbstractArray,Tuple}) = (r = Base.copy(a); clamp!(r,r,mi,ma); r)
clamp!(a, mi, ma) = clamp!(a, a, mi, ma)
function clamp!(r, a, mi, ma)
    if !(size(a,1)==length(mi)==length(ma))
        error("clamp!: size(a,1)==length(mi)==length(ma) was false: $(size(a,1))==$(length(mi))==$(length(ma))")
    end
    for i = 1:length(a)
        r[i] = min(max(mi[i],a[i]), ma[i])
    end
    r
end
clamp(a, v::AbstractArray) = clamp(a, ones(ndims(v),1), siz(v))

function nanfunction(f,a::Array{T,2},dim) where {T}
  if dim == 1
    r = NaN*zeros(1,size(a,2))
    for n = 1:size(a,2)
      v = a[:,n]
      v = v[.!isnan.(v)]
      if !isempty(v)
        r[1,n] = f(v)
      end
    end
    return r
  elseif dim == 2
    r = NaN*zeros(size(a,1),1)
    for m = 1:size(a,1)
      v = a[m,:]
      v = v[.!isnan.(v)]
      if !isempty(v)
        r[m,1] = f(v)
      end
    end
    return r
  else
    error("dim must be 1 or 2")
  end
end
nanfunction(f,a) = f(a[.!isnan.(a)])
nanmedian(a) = nanfunction(median,a)
nanmedian(a,dim) = nanfunction(median,a,dim)
nanmean(a) = nanfunction(mean,a)
nanmean(a,dim) = nanfunction(mean,a,dim)
nanstd(a) = nanfunction(std,a)
nanstd(a,dim) = nanfunction(std,a,dim)

# function nanfunction(f,a::Array{T,3},dim) where {T}
#   if dim == 1
#     r = NaN*zeros(1,size(a,2),size(a,3))
#     for n = 1:size(a,2), o = 1:size(a,3)
#       v = a[:,n,o]
#       v = v[!isnan(v)]
#       if !isempty(v)
#         r[1,n,o] = f(v)
#       end
#     end
#     return r
#   elseif dim == 2
#     r = NaN*zeros(size(a,1),1,size(a,3))
#     for m = 1:size(a,1), o = 1:size(a,3)
#       v = a[m,:,o]
#       v = v[!isnan(v)]
#       if !isempty(v)
#         r[m,1,o] = f(v)
#       end
#     end
#     return r
#   elseif dim == 3
#     r = NaN*zeros(size(a,1),size(a,2),1)
#     for m = 1:size(a,1), n = 1:size(a,2)
#       v = a[m,n,:]
#       v = v[!isnan(v)]
#       if !isempty(v)
#         r[m,n,1] = f(v)
#       end
#     end
#     return r
#   else
#     error("dim must be 1, 2 or 3")
#   end
# end
# nanfunction(f,a) = f(a[!isnan(a)])
# nanmedian(a) = nanfunction(median,a)
# nanmedian(a,dim) = nanfunction(median,a,dim)
# nanmean(a) = nanfunction(mean,a)
# nanmean(a,dim) = nanfunction(mean,a,dim)
# nanstd(a) = nanfunction(std,a)
# nanstd(a,dim) = nanfunction(std,a,dim)
# @test_equal nanmedian([1 2 3], 1)  [1 2 3]
# @test_equal nanmedian([1 2 NaN], 1)  [1 2 NaN]
# @test_equal nanmedian([1 2 3], 2)  [2]
# @test_equal nanmedian([1 2 NaN], 2)  [1.5]
# @test_equal nanmedian([1, 2, NaN])  1.5
# d = rand(2,3,4)
# @test_equal nanmean(d,1) mean(d,1)
# @test_equal nanmean(d,2) mean(d,2)
# @test_equal nanmean(d,3) mean(d,3)

mutable struct WhiteningBase
    mean
    base
    vars
end

function zcawhitening(a; kargs...)
    b = flatten(Base.map(col, row(a)))
    r = zcawhitening(b; kargs...)
    s = size(a[1])
    reshape([reshape(r[:,i], s) for i in 1:length(a)], size(a))
end

pcawhitening(a, base = []; kargs...) = zcawhitening(a, base; pcawhitening = true, kargs...)
function zcawhitening(a::Array{T}, base = [];  perpatchnormalization = false, pcawhitening = false, keepvar = 0.95) where {T<:Real}
    # a ... nDim x nSamples
    if perpatchnormalization
        a = a .- mean(a,1)
        s = std(a,1)
        ind = find(s .> 0.)
        a[:,ind] = a[:,ind] ./ s[1,ind]
    end
    if !isa(base, WhiteningBase)
        m = mean(a,2)
        a = a .- m
        sigma = a * a' ./ size(a, 2)
        (U,S,V) = svd(float(sigma))
        coeff = U' * a
        base = WhiteningBase(m, U, S)
    end
    coeff = base.base' * (a .- base.mean)
    S = base.vars
    r = diagm(1 ./sqrt(base.vars .+ 1e-5)) * coeff
    if !pcawhitening
        r = base.base * r
    end
    r, base
end

distance(a) = distance(a,a)
function distance(a::Vector{T},b::Vector{T}) where {T}
    r = zero(eltype(a))
    for i = 1:length(a)
        r += (a[i]-b[i])^2
    end
    sqrt(r)
end

distance(a::AbstractArray{T1,1},b::AbstractArray{T2,2}) where {T1,T2} = distance(col(a),b)
distance(a::AbstractArray{T1,2},b::AbstractArray{T2,1}) where {T1,T2} = distance(a,col(b))

function distance(a::AbstractArray{T1,2},b::AbstractArray{T2,2}) where {T1,T2}
    #% Author   : Roland Bunschoten
    #%            University of Amsterdam
    #%            Intelligent Autonomous Systems (IAS) group
    #%            Kruislaan 403  1098 SJ Amsterdam
    #%            tel.(+31)20-5257524
    #%            bunschot@wins.uva.nl
    #% Last Rev : Oct 29 16:35:48 MET DST 1999
    #% Tested   : PC Matlab v5.2 and Solaris Matlab v5.3
    #% Thanx    : Nikos Vlassis
    #
    #% Copyright notice: You are free to modify, extend and distribute
    #%    this code granted that the author of the original code is
    #%    mentioned as the original author of the code.
    #%% modified by eva.dittrich@student.tuwien.ac.at due to
    #% post by Oliver Woodford (Nov 2008)
    #% on http://www.mathworks.com/matlabcentral/fileexchange/71

    aa = sum(a.*a, dims = 1)
    bb = sum(b.*b, dims = 1)
    sqrt.(abs.((aa' .+ bb) - 2*a'*b))
end
distance(a::Number,b::Number) = abs(a-b)

randbase(basedim, origdim) = randn(MersenneTwister(0), basedim, origdim)
randnormbase(basedim, origdim) = @p randbase basedim origdim | transpose | map normeuclid | transpose

randproj(a, targetdim::Int) = randbase(targetdim, sizem(a))*a

dropdims_(a) = dropdims(a, dims = ndims(a))
mean_(a) = dropdims_(mean(a, ndims(a)))
median_(a) = dropdims_(median(a, ndims(a)))
std_(a) = dropdims_(std(a, ndims(a)))
var_(a) = dropdims_(var(a, ndims(a)))
min_(a) = dropdims_(minimum(a, ndims(a)))
max_(a) = dropdims_(maximum(a, ndims(a)))
sum_(a) = dropdims_(sum(a, ndims(a)))

normmean_(a) = a .- mean_(a)

ext: jl

max_stars_repo_name: rened/FDUtils

max_issues_repo_name: rened/FDUtils

max_forks_repo_name: rened/FDUtils


text:
abstract Shape{T} <: FieldVector{T}

immutable Cube{T} <: Shape{T}
    a::T
end

immutable Sphere{T} <: Shape{T}
    r::T
end

immutable SphericalCap{T} <: Shape{T}
    a::T
    c::T
    function SphericalCap(a, c)
        if a < 0 || c ≤ 0
            throw(ArgumentError("All lengths should be positive."))
        end
        return new(a, c)
    end
end

for shape in (:Ellipsoid, :EllipticCylinder, :Parallelepiped)
    @eval begin
        immutable $shape{T} <: Shape{T}
            a::T
            b::T
            c::T
            function $shape(a, b, c)
                if a ≤ 0 || b ≤ 0 || c ≤ 0
                    throw(ArgumentError("All lengths should be positive."))
                end
                return new(a, b, c)
            end
        end
    end
end

immutable Cylinder{T} <: Shape{T}
    r::T
    c::T
    function Cylinder(r, c)
        if r ≤ 0 || c ≤ 0
            throw(ArgumentError("All lengths should be positive."))
        end
        return new(r, c)
    end
end

immutable HollowCylinder{T} <: Shape{T}
    R::T
    r::T
    c::T
    function HollowCylinder(R, r, c)
        if r < 0 || c ≤ 0
            throw(ArgumentError("All lengths should be positive."))
        elseif R < r
            throw(ArgumentError("Outer radius `R` should be larger than `r`."))
        end
        return new(R, r, c)
    end
end

immutable SquarePyramid{T,S} <: Shape{T}
    a::T
    b::T
    m::S
    function SquarePyramid(a, b, m)
        if a ≤ 0 || b ≤ 0
            throw(ArgumentError("All lengths should be positive."))
        end
        @assert m == 2b / a
        return new(a, b, m)
    end
end

immutable RectangularPyramid{T,S} <: Shape{T}
    a::T
    b::T
    c::T
    ma::S
    mb::S
    function RectangularPyramid(a, b, c, ma, mb)
        if a ≤ 0 || b ≤ 0 || c ≤ 0
            throw(ArgumentError("All lengths should be positive."))
        end
        @assert ma == 2c / a
        @assert mb == 2c / b
        return new(a, b, c, ma, mb)
    end
end

immutable TriangularToroid{T} <: Shape{T}
    r::T
    b::T
    c::T
    function TriangularToroid(r, b, c)
        if r ≤ 0 || b ≤ 0 || c ≤ 0
            throw(ArgumentError("All lengths should be positive."))
        elseif r < b
            throw(ArgumentError("The radius must be equal or higher than " *
                                "the triangle base semi-length"))
        end
        return new(r, b, c)
    end
end

immutable TruncatedSquarePyramid{T,R,S} <: Shape{T}
    a::T
    b::T
    r::R
    m::S
    function TruncatedSquarePyramid(a, b, r, m)
        if a ≤ 0 || b ≤ 0
            throw(ArgumentError("All lengths should be positive."))
        end
        @assert 0 < r ≤ 1
        @assert m == 2b / a
        return new(a, b, r, m)
    end
end

immutable Ring{T} <: Shape{T}
    R::T
    a::T
    b::T
end

immutable Torus{T} <: Shape{T}
    R::T
    r::T
end

immutable Vec{T} <: FieldVector{T}
    x::T
    y::T
    z::T
end


text:
module AquaTesting

using Base: PkgId, UUID
using Pkg
using Test

# Taken from Test/test/runtests.jl
mutable struct NoThrowTestSet <: Test.AbstractTestSet
    results::Vector
    NoThrowTestSet(desc) = new([])
end
Test.record(ts::NoThrowTestSet, t::Test.Result) = (push!(ts.results, t); t)
Test.finish(ts::NoThrowTestSet) = ts.results

macro testtestset(args...)
    @gensym TestSet
    expr = quote
        $TestSet = $NoThrowTestSet
        $Test.@testset($TestSet, $(args...))
    end
    esc(expr)
end

const SAMPLE_PKGIDS = [
    PkgId(UUID("1649c42c-2196-4c52-9963-79822cd6227b"), "PkgWithIncompatibleTestProject"),
    PkgId(UUID("6e4a843a-fdff-4fa3-bb5a-e4ae67826963"), "PkgWithCompatibleTestProject"),
    PkgId(UUID("7231ce0e-e308-4079-b49f-19e33cc3ac6e"), "PkgWithPostJulia12Support"),
    PkgId(UUID("8981f3dd-97fd-4684-8ec7-7b0c42f64e2e"), "PkgWithoutTestProject"),
    PkgId(nothing, "PkgWithoutProject"),
]

const SAMPLE_PKG_BY_NAME = Dict(pkg.name => pkg for pkg in SAMPLE_PKGIDS)

function with_sample_pkgs(f)
    sampledir = joinpath(@__DIR__, "sample")

    original_load_path = copy(LOAD_PATH)
    try
        pushfirst!(LOAD_PATH, sampledir)
        f()
    finally
        append!(empty!(LOAD_PATH), original_load_path)
    end
end

end  # module


text:
using Test
using BinaryProvider, FileIO, Random, Pkg
include("../examples/library.jl")

record_reference_images = get(ENV, "RECORD_EXAMPLES", false) == "true"
version = v"0.0.7"

download_dir = joinpath(@__DIR__, "testimages")
tarfile = joinpath(download_dir, "images.zip")
url = "https://github.com/SimonDanisch/ReferenceImages/archive/v$(version).tar.gz"
refpath = joinpath(download_dir, "ReferenceImages-$(version)", "images")
recordpath = joinpath(homedir(), "ReferenceImages", "images")
if record_reference_images
    cd(homedir()) do
        isdir(dirname(recordpath)) || run(`git clone git@github.com:SimonDanisch/ReferenceImages.git`)
        isdir(recordpath) && rm(recordpath, recursive = true, force = true)
        mkdir(recordpath)
    end
end

# function url2hash(url::String)
#     path = download(url)
#     open(io-> bytes2hex(BinaryProvider.sha256(io)), path)
# end
# url2hash(url) |> println

if !record_reference_images
    if get(ENV, "USE_REFERENCE_IMAGES", "false") == "true"
        @info("Using Local reference image repository")
        refpath = recordpath
    elseif !isdir(refpath)
        download_images() = BinaryProvider.download_verify(
            url, "dabc4c35eb82a801708596964dc35c0acf1879b487a35b73322f094f611f19c3",
            tarfile
        )
        try
            download_images()
        catch e
            if isa(e, ErrorException) && occursin("Hash Mismatch", e.msg)
                rm(tarfile, force = true)
                download_images()
            else
                rethrow(e)
            end
        end
        BinaryProvider.unpack(tarfile, download_dir)
        # check again after download
        if !isdir(refpath)
            error("Something went wrong while downloading reference images. Plots can't be compared to references")
        end
    end
else
    refpath = recordpath
end

function toimages(f, example, x::Scene, record)
    image = Makie.scene2image(x)
    rpath = joinpath(refpath, "$(example.unique_name).jpg")
    if record
        FileIO.save(joinpath(recordpath, "$(example.unique_name).jpg"), image)
    else
        path = joinpath(refpath, "$(example.unique_name).jpg")
        if isfile(path)
            refimage = FileIO.load(path)
            f(image, refimage)
        else
            @warn("No reference image for $(example.unique_name) - skipping test")
            FileIO.save(joinpath(test_diff_path, "no_reference_$(basename(path))"), image)
        end
    end
end

is_image_file(path) = lowercase(splitext(path)[2]) in (".png", ".jpg", ".jpeg")

function toimages(f, example, s::Stepper, record)
    ispath(s.folder) || error("Not a path: $(s.folder)")
    if record
        # just copy the stepper files from s.folder into the recordpath
        rpath2 = joinpath(recordpath, basename(s.folder))
        cp(s.folder, rpath2, force = true)
    else
        for frame in readdir(s.folder)
            is_image_file(frame) || continue
            path = joinpath(s.folder, frame)
            pathref = joinpath(refpath, basename(s.folder), frame)
            if isfile(pathref)
                image = FileIO.load(path)
                refimage = FileIO.load(pathref)
                f(image, refimage)
            else
                warned || @warn("No reference frames for $(example.unique_name) found - skipping tests")
                cp(path, joinpath(test_diff_path, "no_reference_$(basename(path))"))
                break
            end
        end
    end
end
function toimages(f, example, path::String, record)
    isfile(path) || error("Not a file: $path")
    filepath, ext = splitext(path)
    rpath = joinpath(refpath, basename(filepath))

    if record
        rpath2 = joinpath(recordpath, basename(filepath))
        isdir(rpath2) || mkpath(rpath2)
        run(`ffmpeg -loglevel quiet -i $(abspath(path)) -y $rpath2\\frames%04d.jpg`)
    else
        filepath, ext = splitext(path)
        isdir(filepath) || mkdir(filepath)
        run(`ffmpeg -loglevel quiet -i $path -y $filepath\\frames%04d.jpg`)
        for frame in readdir(filepath)
            pathref = joinpath(refpath, basename(filepath), frame)
            if isfile(pathref)
                image = FileIO.load(joinpath(filepath, frame))
                refimage = FileIO.load(pathref)
                f(image, refimage)
            else
                warned || @warn("No reference frames for $(example.unique_name) found - skipping tests")
                cp(joinpath(filepath, frame), joinpath(test_diff_path, "no_reference_$(frame)"))
                break
            end
        end
    end
end

include("visualregression.jl")

test_diff_path = joinpath(dirname(pathof(Makie)), "..", "test", "testresults")

rm(test_diff_path, force = true, recursive = true)
mkpath(test_diff_path)

function test_examples(record, tags...; kw_args...)
    Random.seed!(42)
    @testset "Visual Regression" begin
        eval_examples(tags..., replace_nframes = true, outputfile = (entry, ending)-> "./media/" * string(entry.unique_name, ending); kw_args...) do example, value
            sigma = [1,1]; eps = 0.02
            maxdiff = 0.05
            toimages(example, value, record) do image, refimage
                @testset "$(example.title):" begin
                    diff = approx_difference(image, refimage, sigma, eps)
                    if diff >= maxdiff
                        save(
                            joinpath(test_diff_path, "$(example.unique_name)_differ.jpg"),
                            hcat(image, refimage)
                        )
                    end
                    @test diff < maxdiff
                end
            end
            # reset global states
            Random.seed!(42)
            AbstractPlotting.set_theme!(resolution = (500, 500))
        end
    end
end

cd(@__DIR__)
isdir("media") || mkdir("media")
isdir("testresults") || mkdir("testresults")
AbstractPlotting.set_theme!(resolution = (500, 500))

@info("Number of examples in database: $(length(database))")

exclude_tags = ["bigdata"]
@info("Excluding tags: $exclude_tags")

indices_excluded = []
for tag in exclude_tags
    global indices_excluded
    indices = find_indices(tag)
    indices_excluded = vcat(indices_excluded, indices)
end
num_excluded = length(unique(indices_excluded))
@info("Number of examples to be skipped: $(num_excluded)")

# run the tests
test_examples(record_reference_images; exclude_tags = exclude_tags)


text:
function shock_loading(m::HetDSGELag, TTT_jump::Matrix{Float64})
    exo = m.exogenous_shocks
    endo = m.endogenous_states

    _RRR = zeros(n_backward_looking_states(m), n_shocks_exogenous(m))
    RRR  = zeros(n_model_states(m), n_shocks_exogenous(m))

    _RRR[endo[:B′], exo[:b_sh]] .= 1.
    _RRR[endo[:G′], exo[:g_sh]] .= 1.
    _RRR[endo[:z′_t], exo[:z_sh]] .= 1.
    _RRR[endo[:MU′], exo[:μ_sh]] .= 1.
    _RRR[endo[:LAMW′], exo[:λ_w_sh]] .= 1.
    _RRR[endo[:LAMF′], exo[:λ_f_sh]] .= 1.
    _RRR[endo[:MON′], exo[:MON_sh]] .= 1.

    # Loading on states and jumps
    RRR[1:n_backward_looking_states(m), :] = _RRR
    RRR[n_backward_looking_states(m)+1:end, :] = TTT_jump*_RRR

    return RRR
end


text:
using DelimitedFiles
using Statistics

# define my colors -
negative_color = (1/255)*[255,51,51]
positive_color = (1/255)*[51,153,255]

function reduce(path_to_sensitivity_array::String; epsilon::Float64=1e-2)

    # load -
    sensitivity_array = readdlm(path_to_sensitivity_array)

    # parameters to kepp -
    index_keep_p = Int64[]

    # the parameters are on the rows - are any values abpove a threshold -
    (NR,NC) = size(sensitivity_array)
    for row_index = 1:NR

        # get the row data -
        row_data = sensitivity_array[row_index, :]

        # check -
        idx_check = findall(x->x>=epsilon, row_data)
        if (isempty(idx_check) == false)
            push!(index_keep_p, row_index)
        end
    end

    # pull out rows -
    reduced_array = sensitivity_array[index_keep_p, :]

    # last thing - lets put things in groups of 0,1,2,3 (or none,low,medium,high)
    (NR,NC) = size(reduced_array)
    category_array = zeros(NR,NC)
    for row_index = 1:NR
        for col_index = 1:NC

            old_value = reduced_array[row_index, col_index]
            order_of_magnitude = floor(Int,log10(old_value))
            oom_eps = floor(Int,log10(epsilon))

            # rules -
            if (order_of_magnitude<-4)                                  # below -4
                category_array[row_index, col_index] = 0
            elseif (order_of_magnitude == -4)                           # O(-3)
                category_array[row_index, col_index] = 1
            elseif (order_of_magnitude == -3)                           # O(-2)
                category_array[row_index, col_index] = 2
            elseif (order_of_magnitude == -2)                           # O(-1)
                category_array[row_index, col_index] = 3
            elseif (order_of_magnitude == -1 )                           # O(0)
                category_array[row_index, col_index] = 4
            elseif (order_of_magnitude == 0)                            # O(1)
                category_array[row_index, col_index] = 5
            elseif (order_of_magnitude == 1)                            # O(2)
                category_array[row_index, col_index] = 6
            elseif (order_of_magnitude == 2)                            # O(3)
                category_array[row_index, col_index] = 7
            else
                category_array[row_index, col_index] = 8
            end
        end
    end

    # return -
    return category_array
end

function visualize(sarray)

    epsilon = 0.2;
    (NR,NC) = size(sarray)

    for i ∈ 1:NC

        for j ∈ 1:NR

            # compute origin point -
            x = [
                (i - 1), (i - 1)*epsilon + 1, (i - 1)*epsilon + 1, (i - 1)*epsilon + 1, (i - 1), (i - 1), (i - 1), (i - 1)
            ];

            y = [
                (j - 1), (j - 1), (j - 1), (j - 1)*epsilon + 1, (j - 1)*epsilon + 1, (j - 1)*epsilon + 1, (j - 1), (j - 1)
            ];

            t = Shape(x,y)
            plot!(t, legend=false)
        end
    end
    current()

end

# setup path to sensitvity array -
path_to_sensitivity_array = "./Sens.dat"

# execute -
rsa = reduce(path_to_sensitivity_array; epsilon=1e-6)

max_stars_repo_name: varnerlab/UVM-TopDown-LegacyModel

max_issues_repo_name: varnerlab/UVM-TopDown-LegacyModel

max_forks_repo_name: varnerlab/UVM-TopDown-LegacyModel


text:
using CSV, DataFrames
#=
Load the Salary and Wages data and Components of Value Added data
=#
oews = CSV.read(joinpath("data", "oews_15-1256.csv"), DataFrame)
cva = CSV.read(joinpath("data", "comp_of_va.csv"), DataFrame)
bls = CSV.read(joinpath("data", "bls_salary_wages_to_total_compensation.csv"), DataFrame)
usetbl = CSV.read(joinpath("data", "usetbl5415.csv"), DataFrame)
klems = CSV.read(joinpath("data", "klems.csv"), DataFrame)

#=
This is the approach based on Carol.
=#
car_sw = transform(oews, [:tot_emp, :a_mean] => ByRow(*) => :s_w)
car_sw = combine(groupby(car_sw, :year),
                 [:tot_emp, :s_w] .=> sum,
                 renamecols = false)
select!(car_sw, :year => identity, [:s_w, :tot_emp] => ByRow(/) => :s_w, renamecols = false)
car = innerjoin(car_sw, bls, on = :year)
select!(car,
        :year => identity,
        [:s_w, :salary_to_total_compensation] => ByRow(*) => :coe,
        renamecols = false)
car = innerjoin(car, usetbl, on = :year)
transform!(car,
           [:coe, :ti_coe] => ByRow(*) => :ti,
           [:coe, :gos_coe] => ByRow(*) => :gos,
           [:coe, :top_coe] => ByRow(*) => :top)
select!(car,
        :year => identity,
        [:coe, :ti, :gos, :top] => ByRow(+) => :car,
        renamecols = false)
transform!(car, :car => ByRow(x -> x / 12), renamecols = false)
#=
Alternative Approach
=#
labor = innerjoin(oews, cva, on = [:year, :naics])
employment = combine(groupby(labor, :year), :tot_emp => sum, renamecols = false)
select!(labor,
        [:year, :naics] .=> identity,
        [:tot_emp, :a_mean, :l_w] => ByRow(*) => :coe, renamecols = false)
jbsc = innerjoin(labor, klems, on = [:year, :naics])
transform!(jbsc,
           [:coe, :top_coe] => ByRow((coe, top_coe) -> coe * max(top_coe, 0)) => :top,
           [:coe, :enepsi_coe] => ByRow(*) => :ii,
           [:coe, :gos_coe] => ByRow(*) => :gos)
transform!(jbsc, [:coe, :gos, :top, :ii] => ByRow(+) => :y)
jbsc = combine(groupby(jbsc, :year), :y => sum, renamecols = false)
jbsc = innerjoin(jbsc, employment, on = :year)
select!(jbsc,
        :year => identity,
        [:y, :tot_emp] => ByRow((y, tot_emp) -> y / tot_emp / 12) => :jbsc,
        renamecols = false)

resource_cost = innerjoin(car, jbsc, on = :year)
transform!(resource_cost,
           [:car, :jbsc] .=> ByRow(x -> round(Int, x)),
           renamecols = false)
#=
Save the final result 🎊
=#
CSV.write(joinpath("data", "person_monthly_resource_cost.csv"), resource_cost)

max_stars_repo_name: Nosferican/Estimating-Resource-Costs

max_issues_repo_name: Nosferican/Estimating-Resource-Costs

max_forks_repo_name: Nosferican/Estimating-Resource-Costs


text:
"""
    DiamondSquare

This type generates a neutral landscape using the diamond-squares
algorithm, which produces fractals with variable spatial autocorrelation.

https://en.wikipedia.org/wiki/Diamond-square_algorithm

The algorithm is named diamond-square because it is an iterative procedure of
"diamond" and "square" steps.

The degree of spatial autocorrelation is controlled by a parameter `H`,
which varies from 0.0 (low autocorrelation) to 1.0 (high autocorrelation) ---
note this is non-inclusive and H = 0 and H = 1 will not behave as expected.
The result of the diamond-square algorithm is a fractal with dimension D = 2 + H.

A similar algorithm, midpoint-displacement, functions almost
identically, except that in DiamondSquare, the square step interpolates
edge midpoints from the nearest two corners and the square's center, where as
midpoint-displacement only intepolates from the nearest corners (see `MidpointDisplacement`).

"""
struct DiamondSquare <: NeutralLandscapeMaker
    H::Float64
    function DiamondSquare(H::T) where {T <: Real}
        @assert 0 <=  H < 1
        new(H)
    end
end

"""
    DiamondSquare()

    Initializes DiamondSqaure with H = 0.5.
"""
DiamondSquare() = DiamondSquare(0.5)

"""
    MidpointDisplacement()

    Creates a midpoint-displacement algorithm object `MidpointDisplacement`. The degree of spatial autocorrelation is controlled by a parameter `H`,
    which varies from 0.0 (low autocorrelation) to 1.0 (high autocorrelation) --- note this is non-inclusive and H = 0 and H = 1 will not behavive as expected.

    A similar algorithm, diamond-square, functions almost
    identically, except that in diamond-square, the square step interpolates
    edge midpoints from the nearest two corners and the square's center, where as
    `MidpointDisplacement` only intepolates from the nearest corners (see `DiamondSquare`).
"""
struct MidpointDisplacement <: NeutralLandscapeMaker
    H::Float64
    function MidpointDisplacement(H::T) where {T <: Real}
        @assert 0 <= H < 1
        new(H)
    end
end
MidpointDisplacement() = MidpointDisplacement(0.5)

"""
    _landscape!(mat, alg::Union{DiamondSquare, MidpointDisplacement}; kw...)

    Check if `mat` is the right size.
    If mat is not the correct size (DiamondSquare can only run on a lattice of size NxN where N = (2^n)+1 for integer n),
    allocates the smallest lattice large enough to contain `mat` that can run DiamondSquare.
"""
function _landscape!(mat, alg::Union{DiamondSquare, MidpointDisplacement}; kw...) where {IT <: Integer}

    rightSize::Bool = _isPowerOfTwo(size(mat)[1]-1) && _isPowerOfTwo(size(mat)[2]-1)
    latticeSize::Int = size(mat)[1]

    dsMat = mat
    if !rightSize
        dim1, dim2 = size(mat)
        smallestContainingLattice::Int = 2^ceil(log2(max(dim1, dim2))) + 1
        @warn "$alg cannot be run on the input dimensions ($dim1 x $dim2),
                and will instead run on the next smallest valid size ($smallestContainingLattice x $smallestContainingLattice).
                This can slow performance as it involves additional memory allocation."
        dsMat = zeros(smallestContainingLattice, smallestContainingLattice)
    end
    _diamondsquare!(dsMat, alg)

    mat .= dsMat[1:size(mat)[1], 1:size(mat)[2]]
end

"""
    _diamondsquare!(mat, alg)

    Runs the diamond-square algorithm on a matrix `mat` of size
    `NxN`, where `N=(2^n)+1` for some integer `n`, i.e (N=5,9,17,33,65)

    Diamond-square is an iterative procedure, where the lattice is divided
    into subsquares in subsequent rounds. At each round, the subsquares shrink in size,
    as previously uninitialized values in the lattice are interpolated as a mean of nearby points plus random displacement.
    As the rounds increase, the magnitude of this displacement decreases. This creates spatioautocorrelation, which is controlled
    by a single parameter `H` which varies between `0` (no autocorrelation) and `1` (high autocorrelation)

"""
function _diamondsquare!(mat, alg)
    latticeSize = size(mat)[1]
    numberOfRounds::Int = log2(latticeSize-1)
    _initializeDiamondSquare!(mat, alg)

    for round in 0:(numberOfRounds-1)    # counting from 0 saves us a headache later
        subsquareSideLength::Int = 2^(numberOfRounds-(round))
        numberOfSubsquaresPerAxis::Int = ((latticeSize-1) / subsquareSideLength)-1

        for x in 0:numberOfSubsquaresPerAxis         # iterate over the subsquares within the lattice at this side length
            for y in 0:numberOfSubsquaresPerAxis
                subsquareCorners = _subsquareCornerCoordinates(x,y,subsquareSideLength)

                _diamond!(mat, alg, round, subsquareCorners)
                _square!(mat, alg, round, subsquareCorners)
            end
        end
    end
end

"""
    _initializeDiamondSquare!(mat, alg)

    Initialize's the `DiamondSquare` algorithm by displacing the four corners of the
    lattice using `displace`, scaled by the algorithm's autocorrelation `H`.
"""
function _initializeDiamondSquare!(mat, alg)
    latticeSize = size(mat)[1]
    corners = _subsquareCornerCoordinates(0,0, latticeSize-1)
    for mp in corners
        mat[mp...] = _displace(alg.H, 1)
    end
end

"""
    _subsquareCornerCoordinates(x::Int, y::Int, sideLength::Int)

    Returns the coordinates for the corners of the subsquare (x,y) given a side-length `sideLength`.
"""
function _subsquareCornerCoordinates(x::Int, y::Int, sideLength::Int)
    corners = [1 .+ sideLength.*i for i in [(x,y), (x+1, y), (x, y+1), (x+1, y+1)]]
end

"""
     _diamond!(mat, alg::DiamondSquare, round::Int, corners::AbstractVector{Tuple{Int, Int}})

     Runs the diamond step of the `DiamondSquare` algorithm on the square defined by
     `corners` on the matrix `mat`. The center of the square is interpolated from the
     four corners, and is displaced. The displacement is drawn according to `alg.H` and round using `displace`
"""
function _diamond!(mat, alg, round::Int, corners::AbstractVector{Tuple{Int, Int}})
    centerPt = _centerCoordinate(corners)
    mat[centerPt...] = _interpolate(mat, corners) + _displace(alg.H, round)
end

"""
    _square!(mat, alg::DiamondSquare, round::Int, corners::AbstractVector{Tuple{Int,Int}})

    Runs the square step of the `DiamondSquare` algorithm on the square defined
    by `corners` on the matrix `mat`. The midpoint of each edge of this square is interpolated
    by computing the mean value of the two corners on the edge and the center of the square, and the
    displacing it. The displacement is drawn according to `alg.H` and round using `displace`

"""
function _square!(mat, alg::DiamondSquare, round::Int, corners::AbstractVector{Tuple{Int, Int}})
    bottomLeft,bottomRight,topLeft,topRight = corners
    leftEdge, bottomEdge, topEdge, rightEdge = _edgeMidpointCoordinates(corners)
    centerPoint = _centerCoordinate(corners)

    mat[leftEdge...] = _interpolate(mat, [topLeft,bottomLeft,centerPoint]) + _displace(alg.H, round)
    mat[bottomEdge...] = _interpolate(mat, [bottomLeft,bottomRight,centerPoint]) + _displace(alg.H, round)
    mat[topEdge...] = _interpolate(mat, [topLeft,topRight,centerPoint]) + _displace(alg.H, round)
    mat[rightEdge...] = _interpolate(mat, [topRight,bottomRight,centerPoint]) + _displace(alg.H, round)
end

"""
    _square!(mat, alg::MidpointDisplacement, round::Int, corners::AbstractVector{Tuple{Int,Int}})

    Runs the square step of the `MidpointDisplacement` algorithm on the square defined
    by `corners` on the matrix `mat`. The midpoint of each edge of this square is interpolated
    by computing the mean value of the two corners on the edge and the center of the square, and the
    displacing it. The displacement is drawn according to `alg.H` and round using `displace`
"""
function _square!(mat, alg::MidpointDisplacement, round::Int, corners::AbstractVector{Tuple{Int, Int}})
    bottomLeft,bottomRight,topLeft,topRight = corners
    leftEdge, bottomEdge, topEdge, rightEdge = _edgeMidpointCoordinates(corners)
    mat[leftEdge...] = _interpolate(mat, [topLeft,bottomLeft]) + _displace(alg.H, round)
    mat[bottomEdge...] = _interpolate(mat, [bottomLeft,bottomRight]) + _displace(alg.H, round)
    mat[topEdge...] = _interpolate(mat, [topLeft,topRight]) + _displace(alg.H, round)
    mat[rightEdge...] = _interpolate(mat, [topRight,bottomRight]) + _displace(alg.H, round)
end

"""
    _interpolate(mat, points::AbstractVector{Tuple{Int,Int}})

    Computes the mean of a set of points, represented as a list of indecies to a matrix `mat`.
"""
function _interpolate(mat, points::AbstractVector{Tuple{Int,Int}})
    return mean(mat[pt...] for pt in points)
end

"""
    _displace(H::Float64, round::Int)

    `displace` produces a random value as a function of  `H`, which is the
    autocorrelation parameter used in `DiamondSquare` and must be between `0`
    and `1`, and `round` which describes the current tiling size for the
    DiamondSquare() algorithm.

    Random value are drawn from a Gaussian distribution using `Distribution.Normal`
    The standard deviation of this Gaussian, σ, is set to (1/2)^(round*H), which will
    move from 1.0 to 0 as `round` increases.

"""
function _displace(H::Float64, round::Int)
    σ = 0.5^(round*H)
    return(rand(Normal(0, σ)))
end

"""
    _centerCoordinate(corners::AbstractVector{Tuple{Int,Int}})

    Returns the center coordinate for a square defined by `corners` for the
    `DiamondSquare` algorithm.
"""
function _centerCoordinate(corners::AbstractVector{Tuple{Int,Int}})
    bottomLeft,bottomRight,topLeft,topRight = corners
    centerX::Int = (_xcoord(bottomLeft)+_xcoord(bottomRight)) ÷ 2
    centerY::Int = (_ycoord(topRight)+_ycoord(bottomRight)) ÷ 2
    return (centerX, centerY)
end

"""
    _xcoord(pt::Tuple{Int,Int})

    Returns the x-coordinate from a lattice coordinate `pt`.
"""
 _xcoord(pt::Tuple{Int,Int}) = pt[1]

 """
     _ycoord(pt::Tuple{Int,Int})

     Returns the y-coordinate from a lattice coordinate `pt`.
 """
 _ycoord(pt::Tuple{Int,Int}) = pt[2]

"""
    _edgeMidpointCoordinates(corners::AbstractVector{Tuple{Int,Int}})

    Returns an array of midpoints for a square defined by `corners` for the `DiamondSquare` algorithm.
"""
function _edgeMidpointCoordinates(corners::AbstractVector{Tuple{Int,Int}})
    # bottom left, bottom right, top left, top right
    bottomLeft,bottomRight,topLeft,topRight = corners

    leftEdgeMidpoint::Tuple{Int,Int} = (_xcoord(bottomLeft), (_ycoord(bottomLeft)+_ycoord(topLeft))÷2 )
    bottomEdgeMidpoint::Tuple{Int,Int} = ( (_xcoord(bottomLeft)+ _xcoord(bottomRight))÷2, _ycoord(bottomLeft) )
    topEdgeMidpoint::Tuple{Int,Int} = ( (_xcoord(topLeft)+_xcoord(topRight))÷2, _ycoord(topLeft))
    rightEdgeMidpoint::Tuple{Int,Int} = ( _xcoord(bottomRight), (_ycoord(bottomRight)+_ycoord(topRight))÷2)

    edgeMidpoints = [leftEdgeMidpoint, bottomEdgeMidpoint, topEdgeMidpoint, rightEdgeMidpoint]
    return edgeMidpoints
end

"""
     _isPowerOfTwo(x::IT) where {IT <: Integer}

     Determines if `x`, an integer, can be expressed as `2^n`, where `n` is also an integer.
"""
function _isPowerOfTwo(x::IT) where {IT <: Integer}
    return x & (x-1) == 0
end


text:
# test for hybridatnode
# still not automatic test function
# Claudia April 2016

include("../src/types.jl")
include("../src/functions.jl")

net = readTopologyLevel1("(Gma,(Gch,((Gmo_C7,((Gmo_N2,Gmo_I3):0.037,Gmo_I7):0.181):2.318,#H8:::0.03):2.892):2.347,(Bsa)#H8:::0.97);")
plot(net)
plot(net, showNodeNumber=true)
# want to change hybrid to -3

net2 = hybridatnode(net,-3)
plot(net2)
plot(net2, showEdgeNumber=true)
rootonedge!(net2,13) #Bsa outgroup
plot(net2)

hybridatnode!(net,-3)
plot(net)

net = readTopology("(A:1.0,(B,(((C:0.52,(E:0.5)#H2:0.02::0.7):0.6,(#H2:0.01::0.3,F:0.7):0.8):0.9,D:0.8):1.3):0.7):0.1;");
plot(net, showNodeNumber=true)
hybridatnode!(net, -7)
plot(net)


text:
#TODO: make ε a tensor
"""
    module DielectricFunctions

for defining dielectric and pump functions either as piecewise constant
or with user-supplied functions
"""
module DielectricFunctions

using ..Points

export DielectricFunction,
PumpFunction,
piecewise_constant_ε,
piecewise_constant_F

"""
    struct DielectricFunction

---

    DielectricFunction(n1=1,n2=0) -> dielectric
    DielectricFunction(n::Complex) -> dielectric

piecewise constant dielectric function with index given by `n1+1im*n2`

---

    DielectricFunction(ε::Function, parameters::Dict) -> dielectric

ε is any function with the signature ε(x,y,parameters::Dict{Symbol,Float64})
which evaluates to the dielectric at position `x,y`. For example:
    ε(x,y,params) = 2+sin(2πx/params[:period])

---

    (::DielectricFunction)(x,y) -> ε
"""
struct DielectricFunction{TF}
    ε::TF
    parameters::Dict{Symbol,Float64}

    function DielectricFunction(ε::TF, parameters::Dict{Symbol,T}) where {T<:Number,TF<:Function}
        if TF<:typeof(piecewise_constant_ε)
            if haskey(parameters,:n₁)
                get!(parameters,:n1,parameters[:n₁])
            elseif haskey(parameters,:n1)
                get!(parameters,:n₁,parameters[:n1])
            else
                throw(ErrorException("Dict parameters must contain at least one of :n₁, :n1"))
            end
            if haskey(parameters,:n₂)
                get!(parameters,:n2,parameters[:n₂])
            elseif haskey(parameters,:n2)
                get!(parameters,:n₂,parameters[:n2])
            else
                throw(ErrorException("Dict parameters must contain at least one of :n₂, :n2"))
            end
        end
        return new{TF}(ε,parameters)
    end
end
DielectricFunction() = DielectricFunction(piecewise_constant_ε,Dict(:n₁=>1,:n₂=>0))
DielectricFunction(n1::Real,n2::Real=0.0) = DielectricFunction(piecewise_constant_ε, Dict{Symbol,Float64}(:n₁=>n1,:n₂=>n2))
DielectricFunction(n::Complex) = DielectricFunction(real(n),imag(n))
DielectricFunction(ε::Function,n1::Real,n2::Real=0.0) = DielectricFunction(ε, Dict{Symbol,Float64}(:n₁=>n1,:n₂=>n2))
(df::DielectricFunction)(p::Point) = df.ε(p,df.parameters)

Base.getindex(de::DielectricFunction,sym::Symbol) = Base.getproperty(de,sym)
function Base.getproperty(de::DielectricFunction,sym::Symbol)
    if Base.sym_in(sym,(:n1,:n₁))
        return de.parameters[:n₁]
    elseif Base.sym_in(sym,(:n2,:n₂))
        return de.parameters[:n₂]
    else
        getfield(de,sym)
    end
end

function Base.show(io::IO,d::DielectricFunction)
    println(io, "DielectricFunction: ")
    if d.ε==piecewise_constant_ε
        print(io, "\tpiecewise_constant_ε: n = ",d.parameters[:n₁]+1im*d.parameters[:n₂])
    else
        for p ∈ d.parameters
            print(io, "\t", p.first, ": ")
            println(io, p.second)
        end
    end
end

"""
    struct PumpFunction

---

    PumpFunction(F=0) -> pump

piecewise constant pump function with pump parameter given by `F`

---

    PumpFunction(F::Function, parameters::Dict) -> pump

F is any real function with the signature F(x,y,parameters::Dict{Symbol,Float64})
which evaluates to the pump profile at position `x,y`. For example:
    F(x,y,params) = sin(2πx/params[:period])

---

    (::PumpFunction)(x,y) -> F
"""
struct PumpFunction{TF}
    F::TF
    parameters::Dict{Symbol,Float64}

    PumpFunction(F::TF, parameters::Dict) where TF<:Function = new{TF}(F,parameters)
end
PumpFunction() = PumpFunction(piecewise_constant_F,Dict{Symbol,Float64}(:F=>0))
PumpFunction(F::Real) = PumpFunction(piecewise_constant_F,Dict(:F=>F))
(pf::PumpFunction)(p::Point) = pf.F(p,pf.parameters)

function Base.show(io::IO,d::PumpFunction)
    println(io, "PumpFunction: ")
    if d.F==piecewise_constant_F
        print(io, "\tpiecewise_constant_F: ",d.parameters[:F])
    else
        for p ∈ d.parameters
            print(io, "\t", p.first, ": ")
            println(io, p.second)
        end
    end
end

"""
    piecewise_constant_ε(point,parameters) -> ε
"""
function piecewise_constant_ε(args...)
    parameters = args[end]
    n1 = get!(parameters,:n₁,nothing)
    n1 = get!(parameters,:n1,nothing)
    n2 = get!(parameters,:n₂,n1)
    n2 = get!(parameters,:n2,n2)
    return complex(n1,n2)^2
end

"""
    piecewise_constant_F(point,parameters) -> F
"""
piecewise_constant_F(args...) = get!(args[end],:F,nothing)

function Base.conj(de::DielectricFunction{typeof(piecewise_constant_ε)})
    return DielectricFunction(de.n1,-de.n2)
end

end # module


text:
build = `../build_cxx_dependencies.sh`
run(build);


text:
##################################################################################################
# Integration with UncertainData.jl (sampling from full supports of the furnishing distributions)
#################################################################################################

import HypothesisTests: OneSampleTTest

# Vectors of uncertain values
uvals_x = [UncertainValue(Normal, rand(Normal(0, 5)), abs(rand(Normal(0, 3)))) for i = 1:100]
uvals_y = [UncertainValue(Normal, rand(Normal(0, 5)), abs(rand(Normal(0, 3)))) for i = 1:100];

# UncertainValueDataset
UVX = UncertainValueDataset(uvals_x)
UVY = UncertainValueDataset(uvals_y)

# UncertainIndexDataset
UVX_idx = UncertainIndexDataset(uvals_x)
UVY_idx = UncertainIndexDataset(uvals_y)

# Real-valued vectors
x = resample.(uvals_x);
y = resample.(uvals_y);

jtest = JointDistanceDistributionTTest()

@test causality(x, y, jtest) isa OneSampleTTest
@test causality(uvals_x, uvals_y, jtest) isa OneSampleTTest
@test causality(x, uvals_y, jtest) isa OneSampleTTest
@test causality(uvals_x, y, jtest) isa OneSampleTTest
@test causality(UVX, UVY, jtest) isa OneSampleTTest
@test causality(x, UVY, jtest) isa OneSampleTTest
@test causality(UVX, y, jtest) isa OneSampleTTest

################################################################
# Integration with UncertainData.jl (with sampling constraints)
################################################################
onevar_constraints = ConstrainedResampling(TruncateStd(1))
twovar_constraints = ConstrainedResampling(TruncateStd(2), TruncateStd(1))

test_jddt = JointDistanceDistributionTTest()
@test causality(x, y, test_jddt, twovar_constraints) isa OneSampleTTest
@test causality(uvals_x, uvals_y, test_jddt, twovar_constraints) isa OneSampleTTest
@test causality(x, uvals_y, test_jddt, onevar_constraints) isa OneSampleTTest
@test causality(uvals_x, y, test_jddt, onevar_constraints) isa OneSampleTTest
@test causality(UVX, UVY, test_jddt, twovar_constraints) isa OneSampleTTest
@test causality(uvals_x, UVY, test_jddt, twovar_constraints) isa OneSampleTTest
@test causality(UVX, uvals_y, test_jddt, twovar_constraints) isa OneSampleTTest
@test causality(x, UVY, test_jddt, onevar_constraints) isa OneSampleTTest
@test causality(UVX, y, test_jddt, onevar_constraints) isa OneSampleTTest


text:
is_ch_type(::Val{:Date})  = true
result_type(::Val{:Date})  = Vector{Date}

function read_col_data(sock::ClickHouseSock, num_rows::VarUInt, ::Val{:Date})
    data = chread(sock, Vector{UInt16}, num_rows)
    return Date(1970) + Day.(data)

end

function write_col_data(sock::ClickHouseSock,
            data::AbstractVector{Date}, ::Val{:Date})
    d = Vector{UInt16}(undef, length(data))
    d .= Dates.value.(data .- Date(1970))
    chwrite(sock, d)
end


text:
"""
    icrs2g(dr, dd)

Transformation from ICRS to Galactic Coordinates.

### Given ###

- `dr`: ICRS right ascension (radians)
- `dd`: ICRS declination (radians)

### Returned ###

- `dl`: Galactic longitude (radians)
- `db`: Galactic latitude (radians)

### Notes ###

1. The IAU 1958 system of Galactic coordinates was defined with
   respect to the now obsolete reference system FK4 B1950.0.  When
   interpreting the system in a modern context, several factors have
   to be taken into account:

   - The inclusion in FK4 positions of the E-terms of aberration.

   - The distortion of the FK4 proper motion system by differential
     Galactic rotation.

   - The use of the B1950.0 equinox rather than the now-standard
     J2000.0.

   - The frame bias between ICRS and the J2000.0 mean place system.

   The Hipparcos Catalogue (Perryman & ESA 1997) provides a rotation
   matrix that transforms directly between ICRS and Galactic
   coordinates with the above factors taken into account.  The
   matrix is derived from three angles, namely the ICRS coordinates
   of the Galactic pole and the longitude of the ascending node of
   the galactic equator on the ICRS equator.  They are given in
   degrees to five decimal places and for canonical purposes are
   regarded as exact.  In the Hipparcos Catalogue the matrix
   elements are given to 10 decimal places (about 20 microarcsec).
   In the present ERFA function the matrix elements have been
   recomputed from the canonical three angles and are given to 30
   decimal places.

2. The inverse transformation is performed by the function [`g2icrs`](@ref).

### Called ###

- [`anp`](@ref): normalize angle into range 0 to 2pi
- [`anpm`](@ref): normalize angle into range +/- pi
- [`s2c`](@ref): spherical coordinates to unit vector
- [`rxp`](@ref): product of r-matrix and p-vector
- [`c2s`](@ref): p-vector to spherical

### Reference ###

- Perryman M.A.C. & ESA, 1997, ESA SP-1200, The Hipparcos and Tycho
    catalogues.  Astrometric and photometric star catalogues
    derived from the ESA Hipparcos Space Astrometry Mission.  ESA
    Publications Division, Noordwijk, Netherlands.

"""
function icrs2g(a, b)
    r1 = Ref{Cdouble}()
    r2 = Ref{Cdouble}()
    ccall((:eraIcrs2g, liberfa), Cvoid,
            (Cdouble, Cdouble, Ref{Cdouble}, Ref{Cdouble}),
            a, b, r1, r2)
    return r1[], r2[]
end

"""
    ir()

Initialize an r-matrix to the identity matrix.

!!! warning "Deprecated"
    Use `Array{Float64}(LinearAlgebra.I, 3, 3)` instead.

### Returned ###

- `r`: r-matrix
"""
ir

function _ir()
    r = Array{Cdouble}(undef, 3, 3)
    ccall((:eraIr, liberfa), Cvoid, (Ptr{Cdouble},), r)
    return r
end

@deprecate ir() Array{Float64}(I, 3, 3)


text:
#'---
#'author: Garrett Smith
#'title: Parameter recovery with a hierarchical model
#'---
#'
#' # Parameter recovery test using a hierarchical model
#'
#' The goal of this script is to test the recovery of parameters from a hierarchical model
#' with a `fpdistribution` likelihood.

#+ echo=false
using Distributions
using Turing
using Zygote
Turing.setadbackend(:zygote)
using Plots, StatsPlots
using Pkg
Pkg.activate("../../FirstPassageTools.jl/")
using FirstPassageTools

#' ## Generating the true data
#'
#' First, we set up the transition rate matrices for the first-passage time distribution we
#' want to fit.

#T = 4*[-1.0 0 0; 1 -1 1; 0 1 -2]
T = 4*[-1 1.0 0; 1 -2 1; 0 1 -2]
A = 4*[0 0 1.0]
p0 = [1.0, 0, 0]

#' Scaling the transition rate matrices by `τ = 2.5` should give mean first-passage times of
#' around 400ms. Generating and fitting the paramters (τ and the separate τᵢ) will be done
#' on the log scale and then exponentiated in order to keep the transition rates positive.

nparticipants = 50
true_tau = 2.5
true_sd = 0.2
true_tau_i = exp.(rand(Normal(0, true_sd), nparticipants));

#' The data will be saved in wide format: Each participant's data corresponds to a row, and
#' each column is a data point.

#+ echo=false
ndata = 10
data = zeros(nparticipants, ndata)
param = true_tau .* true_tau_i
for i = 1:nparticipants
    data[i,:] = rand(fpdistribution(param[i]*T, param[i]*A, p0), ndata)
end

#' ## Specifying the model
#'
#' Now, we can write the full model including the likelihood. Note that we're using a
#' non-centered parameterization. This is because pilot simulations suggested that
#' sampling was biased in the centered parameterization.

# Switch to param = exp(tau) + exp(tau_i). This will prevent really small params b/c multiplication.
@model function mod(y, Tmat=T, Amat=A, p0vec=p0)
    np, nd = size(y)
    # Priors
    # Using the non-centered parameterization for τ
    τ ~ Normal()
    τ̂ = 1 + 0.1*τ  # Corresponds to Normal(1, 0.1)
    sd ~ Exponential(0.25)
    τᵢ ~ filldist(Normal(), np)
    τ̂ᵢ = sd .* τᵢ  # Corresponds to MvNormal(0, sd)
    # Likelihood
    mult = exp.(τ̂ .+ τ̂ᵢ)
    y ~ filldist(arraydist([fpdistribution(mult[p]*Tmat, mult[p]*Amat, p0vec) for p in 1:np]), nd)
    return τ̂, τ̂ᵢ, mult
end

#+ echo=false
#@model function mod_centered(y)
#    np, nd = size(y)
#    # Priors
#    # Using the centered parameterization for τ
#    τ ~ Normal(1, 0.1)
#    sd ~ Exponential(0.25)
#    τᵢ ~ filldist(Normal(0, sd), np)
#    # Likelihood
#    mult = exp.(τ .+ τᵢ)
#    y ~ filldist(arraydist([fpdistribution(mult[p]*T, mult[p]*A, p0) for p in 1:np]), nd)
#end
#
#@model function mod_noncentered_tau(y)
#    np, nd = size(y)
#    # Priors
#    # Using the centered parameterization for τ
#    τ ~ Normal()
#    τ̂ = 1.0 + 0.1*τ
#    sd ~ Exponential(0.25)
#    τᵢ ~ filldist(Normal(0, sd), np)
#    # Likelihood
#    mult = exp.(τ̂ .+ τᵢ)
#    y ~ filldist(arraydist([fpdistribution(mult[p]*T, mult[p]*A, p0) for p in 1:np]), nd)
#end
#
#@model function mod_noncentered_tau_i(y)
#    np, nd = size(y)
#    # Priors
#    # Using the centered parameterization for τ
#    τ ~ Normal(1, 0.1)
#    sd ~ Exponential(0.25)
#    τᵢ ~ filldist(Normal(), np)
#    τ̂ᵢ = sd .* τᵢ
#    # Likelihood
#    mult = exp.(τ .+ τ̂ᵢ)
#    y ~ filldist(arraydist([fpdistribution(mult[p]*T, mult[p]*A, p0) for p in 1:np]), nd)
#end
#' ## Sampling
#'
#' Here, we'll use the NUTS sampler with a burnin of 100 samples and an acceptance rate of
#' 0.65 posterior. We'll use four chains of 1000 samples each. Make sure to execute this
#' script with `julia -t 4 HierarchicalParameterRecovery.jl`.

#posterior = sample(mod(data), NUTS(250, 0.65), MCMCThreads(), 1000, 4)
posterior = sample(mod(data), NUTS(250, 0.7), 1000);
#posterior_centered = sample(mod_centered(data), NUTS(100, 0.65), 1000);
#posterior_noncentered_tau = sample(mod_noncentered_tau(data), NUTS(100, 0.65), 1000);
#posterior_noncentered_tau_i = sample(mod_noncentered_tau_i(data), NUTS(100, 0.65), 1000);

#' ## Evaluating parameter recovery
#'
#' First, we summarize the chains:

posterior

#' Let's also look at the Gelman-Rubin statistic for the chains:
#gelmandiag(posterior)

# #' And the centered version:
#posterior_centered

# #' τ non-centered, τᵢ centered:
#posterior_noncentered_tau

# #' τ centered, τᵢ non-centered
#posterior_noncentered_tau_i

#' And plot histograms of the parameters on the millisecond scale:
#+ echo=false
gen = generated_quantities(mod(data), posterior);
gen = vcat(gen...);
gen_tau = [x[1] for x in gen];
gen_tau_i = transpose(hcat([x[2] for x in gen]...));
shrunken = transpose(hcat([x[3] for x in gen]...));

#+ echo=false; fig_cap="Posterior distribution of τ. Vertical line shows the true value"
histogram(exp.(gen_tau), xlabel="τ")
vline!([true_tau])

#+ echo=false; fig_cap="Posterior distribution of the SD. Vertical lines show the true value."
histogram(posterior[:sd][:], xlabel="Std. deviation of the τᵢ")
vline!([true_sd], label="True value")

#+ echo=false; fig_cap="Posterior distributions of the τᵢ. Vertical lines show the true value."
plot_vec = [];
for p = 1:nparticipants
    plt = histogram(exp.(gen_tau_i[:,p]), xticks=false, yticks=false)
    plt = vline!(plt, [true_tau_i[p]])
    push!(plot_vec, plt)
end
plot(plot_vec..., legend=false)

#+ echo=false; fig_cap="Posterior distributions of τ*τᵢ. Vertical lines show the true value."
plot_vec = [];
for p = 1:nparticipants
    plt = histogram(shrunken[:,p], xticks=false, yticks=false)
    plt = vline!(plt, [param[p]])
    push!(plot_vec, plt)
end
plot(plot_vec..., legend=false)

#' If the posterior contains the true values of the parameters, we can say the parameters
#' were recovered successfully.
#'


text:
module DataCollection

export AbstractParameters, AbstractData, Parameters, Data, addrow!, tokw

# stdlib
using Dates
using Serialization

# "external" dependencies
import DataFrames: DataFrame
import DataStructures: SortedDict

#####
##### Experiment Parameters
#####

# Add a hook for modifying values when we turn Parameters structs into dictionaries
lower(x) = x
lower(f::Function) = string(f)
lower(x::AbstractArray) = map(lower, x)

# Use experiment parameters to collect required information for each run.
abstract type AbstractParameters end
abstract type AbstractData end

const Bundle = Union{AbstractParameters,AbstractData}
kwexclude(::Bundle) = ()
paramexclude(x::Bundle) = kwexclude(x)

#####
##### GenericParameters and GenericData
#####

namesof(::NamedTuple{names}) where {names} = names
struct Parameters{names,T} <: AbstractParameters
    nt::NamedTuple{names,T}
end
Parameters(; kw...) = Parameters((; kw...))
const GenericParameters = Parameters

struct Data{names,T} <: AbstractData
    nt::NamedTuple{names,T}
end
Data(; kw...) = Data((; kw...))
const GenericData = Data

Base.propertynames(x::Union{Parameters,Data}) = namesof(x.nt)
function Base.getproperty(x::Union{Parameters,Data}, s::Symbol)
    nt = getfield(x, :nt)
    s == :nt && return nt
    return nt[s]
end

Base.merge(a::Parameters, b::Parameters) = Parameters(merge(a.nt, b.nt))
Base.merge(a::Data, b::Data) = Data(merge(a.nt, b.nt))

# DataFrame hooks
combine_error(a, b) = error("Found duplicate keys! Values are `$a` and `$b`.")
function dict(x::Bundle)
    return SortedDict(
        k => lower(getproperty(x, k)) for k in propertynames(x) if !in(k, paramexclude(x))
    )
end
lower(x::Bundle) = dict(x)

# Support for turning parameters into a collection of keyword arguments
function tokw(x::Bundle)
    return (k => getproperty(x, k) for k in propertynames(x) if !in(k, kwexclude(x)))
end

timestamp() = SortedDict(:timestamp => now())

_lower(x...) = merge(combine_error, map(lower, x)...)
function addrow!(
    df::DataFrame,
    y::AbstractData,
    x::Union{AbstractParameters,<:AbstractDict}...;
    cols = :setequal,
    sanitizer = Returns(nothing),
    force = false,
)
    # if the datacrame is empty, then we want to set `cols` to `:union` in order
    # to bootstrap DataFrame creation
    d = _lower(x...)
    row = findrow(df, d)

    # If a row for these parameters already exists, check if we're forcing update.
    # If so, delete the row and continue.
    #
    # Otherwise, we need to throw an error
    if !isnothing(row)
        if force
            delete!(df, row)
        else
            error("Row for parameters already exists!")
        end
    end

    cols = isempty(df) ? :union : cols
    newrow = _lower(timestamp(), y, d)
    sanitizer(newrow)
    return push!(df, newrow; promote = true, cols = cols)
end

function addrow!(
    savefile::AbstractString,
    y::AbstractData,
    x::Union{AbstractParameters,<:AbstractDict}...;
    kw...,
)
    df = load(savefile)
    addrow!(df, y, x...; kw...)
    save(df, savefile)
    return df
end

# Utilities for searching for existing entries.
struct RowSearch{T}
    dict::T
end

RowSearch(x::AbstractParameters...) = RowSearch(_lower(x...))

# Treat `missing` as falses.
unmissing(x::Bool) = x
unmissing(x::Missing) = false

function (R::RowSearch)(x)
    pred =
        kv -> hasproperty(x, first(kv)) && unmissing(getproperty(x, first(kv)) == last(kv))
    return all(pred, R.dict)
end

findrow(df::DataFrame, x::AbstractParameters...) = findrow(df, _lower(x...))
findrow(df::DataFrame, x::AbstractDict) = findrow(df, RowSearch(x))
function findrow(df::DataFrame, f::RowSearch)
    for (i, row) in enumerate(eachrow(df))
        f(row) && (return i)
    end
    return nothing
end

hasrow(df::DataFrame, x...) = findrow(df, x...) !== nothing

# Do a trick of creating a temp directory and only replacing the original when done.
# This is something I learned a while back.
#
# Sometimes, serialization can error due to interrupt or for some other reason.
#
# If that happens and we're serializing to the same file we deserialized from, that file
# become corrupt and we lose data - which is really sad.
function save(df::DataFrame, path)
    temp = tempname()
    serialize(temp, df)
    mv(temp, path; force = true)
    return nothing
end

load(path::AbstractString) = ispath(path) ? deserialize(path)::DataFrame : DataFrame()

end # module


text:
# define the api for users

"""
    $SIGNATURES

Approximate the pointwise the product of functionals on manifolds using KernelDensityEstimate.

Notes:
- Always pass full beliefs, for partials use for e.g. `partialDimsWorkaround=[1;3;6]`
- Can also multiply different partials together

Example
-------
```julia
# setup
M = TranslationGroup(3)
N = 75
p = manikde!(M, [randn(3) for _ in 1:N])
q = manikde!(M, [randn(3) .+ 1 for _ in 1:N])

# approximate the product between hybrid manifold densities
pq = manifoldProduct([p;q])

# direct histogram plot
using Gadfly
plot( x=getPoints(pq)[1,:], y=getPoints(pq)[2,:], Geom.histogram2d )

# TODO, convenient plotting (work in progress...)
```
"""
function manifoldProduct( ff::AbstractVector{<:ManifoldKernelDensity},
                          mani::M=ff[1].manifold;
                          makeCopy::Bool=false,
                          Niter::Int=1,
                          # partialDimsWorkaround=1:MB.manifold_dimension(mani),
                          ndims::Int=maximum(Ndim.(ff)),
                          N::Int = maximum(Npts.(ff)),
                          u0 = getPoints(ff[1], false)[1],
                          oldPoints::AbstractVector{P}= [identity_element(mani, u0) for i in 1:N],
                          addEntropy::Bool=true,
                          recordLabels::Bool=false,
                          selectedLabels::Vector{Vector{Int}}=Vector{Vector{Int}}(),
                          _randU=Vector{Float64}(),
                          _randN=Vector{Float64}(),
                          logger=ConsoleLogger()  ) where {M <: MB.AbstractManifold, P}
  #
  # check quick exit
  if 1 == length(ff)
    # @show Ndim(ff[1]), Npts(ff[1]), getPoints(ff[1],false)[1]
    return (makeCopy ? x->deepcopy(x) : x->x)(ff[1])
  end

  Ndens = length(ff)
  # Npartials = length(partials)
  Ndims = maximum(Ndim.(ff))
  with_logger(logger) do
    @debug "[x$(Ndens),d$(Ndims),N$(N)],"
  end

  glbs = KDE.makeEmptyGbGlb();
  glbs.recordChoosen = recordLabels

  # TODO DEPRECATE ::NTuple{Symbol} approach
  manif = convert(Tuple, mani)  #[partialDimsWorkaround]
  addopT, diffopT, getManiMu, _ = buildHybridManifoldCallbacks(manif)

  bws = ones(ndims)
  # MAKE SURE inplace ends up as matrix of coordinates from incoming ::Vector{P}
  oldpts = _pointsToMatrixCoords(mani, oldPoints)
  # FIXME currently assumes oldPoints are in coordinates...
  # @cast oldpts_[i,j] := oldPoints[j][i]
  # oldpts = collect(oldpts_)
  inplace = kde!(oldpts, bws, addopT, diffopT ); # rand(ndims,N)

  # TODO REMOVE
  _ff = (x->x.belief).(ff)
  partialDimMask = Vector{BitVector}(undef, length(ff))
  for (k,md) in enumerate(ff)
    partialDimMask[k] = ones(Int,ndims) .== 1
    if isPartial(md)
      for i in 1:ndims
        if !(i in md._partial)
          partialDimMask[k][i] = false
        end
      end
    end
  end

  ndims = maximum(Ndim.(_ff))
  Ndens = length(_ff)
  Np    = Npts(inplace)
  maxNp = maximum(Int[Np; Npts.(_ff)])
  Nlevels = floor(Int,(log(Float64(maxNp))/log(2.0))+1.0)
  if 0 == length(_randU)
    _len = Int(Np*Ndens*(Niter+2)*Nlevels)
    resize!(_randU, _len)
    _randU .= rand(_len)
  end
  if 0 == length(_randN)
    _len = Int(ndims*Np*(Nlevels+1))
    resize!(_randN, _len)
    _randN .= randn(_len)
  end

  ## TODO check both _ff and inplace use a matrix of coordinates (columns)
  # expects Matrix with columns as samples and rows are coordinate dimensions
  pGM, = prodAppxMSGibbsS(inplace, _ff,
                          nothing, nothing, Niter=Niter,
                          partialDimMask=partialDimMask,
                          addop=addopT,
                          diffop=diffopT,
                          getMu=getManiMu,
                          glbs=glbs,
                          addEntropy=addEntropy,
                          ndims=ndims,
                          Ndens=Ndens,
                          Np=Np,
                          maxNp=maxNp,
                          Nlevels=Nlevels,
                          randU=_randU,
                          randN=_randN  );
  #

  if recordLabels
    # how many levels in ball tree
    lc = glbs.labelsChoosen
    nLevels = maximum(keys(lc[1][1]) |> collect)

    # push final label selections onto selectedLabels
    resize!(selectedLabels, N)
    for i in 1:N
      selectedLabels[i] = Int[]
      for j in 1:length(ff)
        push!(selectedLabels[i], lc[i][j][nLevels])
      end
    end
  end

  # if only partials, then keep other dimension values from oldPoints
  otherDims = ones(ndims) .== 0
  for msk in partialDimMask
    otherDims .|= msk
  end

  # build new output ManifoldKernelDensity
  bws[:] = getKDEManifoldBandwidths(pGM, manif)
  bel = kde!(pGM, bws, addopT, diffopT)

  # FIXME u0 might not be representative of the partial information
  return ManifoldKernelDensity(mani, bel, otherDims, ff[1]._u0)
end

# NOTE, this product does not handle combinations of different partial beliefs properly yet
function *(PP::AbstractVector{<:MKD{M,B}}) where {M<:MB.AbstractManifold{MB.ℝ},B}
  manifoldProduct(PP, PP[1].manifold)
end

function *(P1::MKD{M,B}, P2::MKD{M,B}, P_...) where {M<:MB.AbstractManifold{MB.ℝ},B}
  manifoldProduct([P1;P2;P_...], P1.manifold)
end

#


text:
function vtk_grid(
        dtype::VTKRectilinearGrid, filename::AbstractString,
        xs::Tuple{Vararg{AbstractVector, 3}};
        extent = map(length, xs), whole_extent = extent, kwargs...,
    )
    Ns = map(length, xs)
    Npts = prod(Ns)
    Ncls = num_cells_structured(Ns)

    xvtk = XMLDocument()
    vtk = DatasetFile(dtype, xvtk, filename, Npts, Ncls; kwargs...)

    # VTKFile node
    xroot = vtk_xml_write_header(vtk)

    # RectilinearGrid node
    xGrid = new_child(xroot, vtk.grid_type)
    let ext = extent_attribute(whole_extent)
        set_attribute(xGrid, "WholeExtent", ext)
    end

    # Piece node
    xPiece = new_child(xGrid, "Piece")
    let ext = extent_attribute(extent)
        set_attribute(xPiece, "Extent", ext)
    end

    # Coordinates node
    xPoints = new_child(xPiece, "Coordinates")

    # DataArray node
    x, y, z = xs
    data_to_xml(vtk, xPoints, x, "x")
    data_to_xml(vtk, xPoints, y, "y")
    data_to_xml(vtk, xPoints, z, "z")

    vtk
end

"""
    vtk_grid(filename::AbstractString,
             x::AbstractVector{T}, y::AbstractVector{T}, [z::AbstractVector{T}];
             kwargs...)

Create 2D or 3D rectilinear grid (`.vtr`) file.

Coordinates are specified by separate vectors `x`, `y`, `z`.

# Examples

```jldoctest
julia> vtk = vtk_grid("abc", [0., 0.2, 0.5], collect(-2.:0.2:3), [1., 2.1, 2.3])
VTK file 'abc.vtr' (RectilinearGrid file, open)
```

"""
function vtk_grid(
        filename::AbstractString,
        x::AbstractVector{T}, y::AbstractVector{T},
        z::AbstractVector{T} = Zeros{T}(1);
        kwargs...,
    ) where {T}
    vtk_grid(VTKRectilinearGrid(), filename, (x, y, z); kwargs...)
end


text:
@info "Testing pcloud_func..."
@testset "Transforms pcloud_func" begin

    for (inplace, FUNC) in [(true, Flux3D.normalize!), (false, Flux3D.normalize)]
        @testset "$(FUNC) function" begin
            p = rand(Float32, 3, 8, 2)
            pc1 = PointCloud(p)
            pc2 = FUNC(pc1)
            if inplace
                @test pc1 == pc2
            else
                @test pc1 != pc2
            end
            @test all(isapprox.(mean(pc2.points; dims = 2), 0, rtol = 1e-5, atol = 1e-5))
            @test all(isapprox.(std(pc2.points; dims = 2), 1, rtol = 1e-5, atol = 1e-5))
        end
    end

    for (inplace, FUNC) in [(true, Flux3D.scale!), (false, Flux3D.scale)]
        @testset "$(FUNC) function" begin
            p = rand(Float32, 3, 8, 2)
            pc1 = PointCloud(p)
            pc2 = FUNC(FUNC(pc1, 2.0), 0.5)
            if inplace
                @test pc1 == pc2
            else
                @test pc1 != pc2
            end
            @test all(isapprox.(p, pc2.points, rtol = 1e-5, atol = 1e-5))
        end
    end

    for (inplace, FUNC) in [(true, Flux3D.rotate!), (false, Flux3D.rotate)]
        @testset "$(FUNC) function" begin
            p = rand(Float32, 3, 8, 2)
            rotmat = 2 .* one(rand(Float32, 3, 3))
            rotmat_inv = inv(rotmat)
            rotmat_b = cat(rotmat, rotmat, dims = 3)
            rotmat_inv_b = cat(rotmat_inv, rotmat_inv, dims = 3)
            pc1 = PointCloud(p)
            pc2 = FUNC(FUNC(pc1, rotmat), rotmat_inv)
            pc3 = FUNC(FUNC(pc1, rotmat_b), rotmat_inv_b)
            if inplace
                @test pc1 == pc2
                @test pc1 == pc3
            else
                @test pc1 != pc2
                @test pc1 != pc3
            end
            @test all(isapprox.(p, pc2.points, rtol = 1e-5, atol = 1e-5))
            @test all(isapprox.(p, pc3.points, rtol = 1e-5, atol = 1e-5))
        end
    end

    for (inplace, FUNC) in [(true, Flux3D.realign!), (false, Flux3D.realign)]
        @testset "$(FUNC) function" begin
            p1 = rand(Float32, 8, 3, 2)
            p2 = rand(Float32, 8, 3, 1)
            src = PointCloud(p1)
            tgt = PointCloud(p2)
            src_1 = FUNC(src, tgt)
            if inplace
                @test src == src_1
            else
                @test src != src_1
            end
            # Transformed PointCloud should be inside the bounding box defined by `tgt` PointCloud
            @test all(
                reshape(maximum(tgt.points, dims = 2), (:, 1)) .>=
                src_1.points .>=
                reshape(minimum(tgt.points, dims = 2), (:, 1)),
            )
        end
    end
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "alsa"
version = v"1.2.1-1"

# Collection of sources required to complete build
sources = [
    ArchiveSource("ftp://ftp.alsa-project.org/pub/lib/alsa-lib-1.2.1.1.tar.bz2",
                  "c95ac63c0aad43a6ac457d960569096b0b2ef72dc4e3737e77e3e2de87022cec"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/alsa-lib*/
./configure --prefix=${prefix} --build=${MACHTYPE} --host=${target}
make
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = filter(Sys.islinux, supported_platforms())

# The products that we will ensure are always built
products = [
    LibraryProduct("libasound", :libasound),
    LibraryProduct("libatopology", :libatopology),
]

# Dependencies that must be installed before this package can be built
dependencies = [
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: devmotion/Yggdrasil

max_issues_repo_name: devmotion/Yggdrasil

max_forks_repo_name: devmotion/Yggdrasil


text:
module Scenario
include("./DERs.jl")
import .DERs

export *

function renewable(sys, alpha)
    renewable_location = [26, 29, 13, 6]
    renewable_config = DERs.Config(false, true, 0.1, 0.0, alpha, 0.0)
    for location in renewable_location
        push!(sys.dgs, DERs.dg(renewable_config, location))
    end
    return sys
end
function ess(sys, alpha)
    ess_location = [4, 32]
    ess_config = DERs.Config(true, true, 0.05, 0.05, alpha, alpha)
    for location in ess_location
        push!(sys.dgs, DERs.ess(ess_config, location))
    end
    return sys
end
function ev(sys, alpha)
    ev_location = [17, 20]
    ev_config = DERs.Config(false, true, 0.05, 0.0, alpha, 0.0)
    for location in ev_location
        push!(sys.dgs, DERs.ev_charger(ev_config, location))
    end
    return sys
end
function dg_dispached(sys, alpha)
    dg_dispached_location = [18, 33, 22, 25]
    dg_dispached_config = DERs.Config(true, true, 0.1, 100.0, 1.0, 1.0)
    for location in dg_dispached_location
        push!(sys.dgs, DERs.dg(dg_dispached_config, location))
    end
    return sys
end
function ess_dispached(sys, alpha)
    ess_dispached_location = [9]
    ess_dispached_config = DERs.Config(true, true, 0.1, 0.1, 1.0, 1.0)
    for location in ess_dispached_location
        push!(sys.dgs, DERs.ess(ess_dispached_config, location))
    end
    return sys
end
function scenario(alpha::Float64, ders)
    null = (sys, alpha) -> sys
    return mapreduce((func) -> (sys) -> func(sys, alpha), ∘, [ders; null])
end
end


text:
#
# add a random field to hhld and person dataframes
#
using CSV,DataFrames,ScottishTaxBenefitModel.Randoms

people=CSV.File( "data/model_people_scotland.tab" )|>DataFrame
n = size(people)[1]
people.onerand=fill("",n)
for i in 1:n
    people[i,:onerand]=mybigrandstr()
end
CSV.write("data/model_people_scotland.tab", people, delim = "\t")

n = size(hhlds)[1]
hhld.onerand=fill("",n)
for i in 1:n
    hhlds[i,:onerand]=mybigrandstr()
end
CSV.write("data/model_households_scotland.tab", hhlds, delim = "\t")

# ... and so on for examples, UK data as needed.

max_stars_repo_name: grahamstark/ScottishTaxBenefitModel

max_issues_repo_name: grahamstark/ScottishTaxBenefitModel

max_forks_repo_name: grahamstark/ScottishTaxBenefitModel


text:
using Glissa
using Test

@testset "Glissa.jl" begin
    include("test_cubic_spline.jl")
end


text:
module MMWavenet

using BSON
using Shuffle

export main, build_wavenet, train!
export save_dataset, load_dataset
export convert_jsons_2_pickle_samples

include("network.jl")
include("dataset.jl")

@doc """
    main(sample_len::Int64, sample_features::Int64, path::String)

Creates neural network for sample shape ('sample_len', 'sample_features', batch).
Loads BSON dataset from 'path' and train network. Returns network and its parametres.
'itr' is number of training iterations (default 10).
""" ->
function main(sample_len::Int64, sample_features::Int64, path::String; itr::Int64=10)

    wavenet, θ = build_wavenet(sample_features, sample_len);

    trn_dataset, val_dataset = load_dataset(path)

    train!(θ, wavenet, trn_dataset, val_dataset, iter=itr)

    println("Confusion matrix for validation data.")
    print_conmat(Vector{Float64}(wavenet(val_dataset[1])'[:,1]), Vector{Int64}(val_dataset[2]), true)
    print_conmat(Vector{Float64}(wavenet(val_dataset[1])'[:,1]), Vector{Int64}(val_dataset[2]), false)

    return wavenet, θ
end

end

max_stars_repo_name: kozvojtex/MMWavenet

max_issues_repo_name: kozvojtex/MMWavenet

max_forks_repo_name: kozvojtex/MMWavenet


text:
#
# not syntax (to prime or tag the compliment
# of the specified indices/pattern)
#

struct Not{T}
  pattern::T
  Not(p::T) where {T} = new{T}(p)
end

"""
not(p)

Represents the compliment of the input
for pattern matching in priming, tagging
and other IndexSet related functions.
"""
function not end

"""
parent(n::Not)

Get the original pattern.
"""
Base.parent(n::Not) = n.pattern


text:
srand(1)
N = 1_000_000
v = randn(N)
dv = DataArray(v)
dvna = deepcopy(dv)
dvna[rand(1:N, 10_000)] = NA
idxv = shuffle([1:N])
idxdv = DataArray(idxv)

f1(v) = sum(v)
f2(v) = sum(removeNA(v))
f3(v) = sum(removeNA(v)) # Make this an iterator
f4(v) = mean(v)
f5(v) = mean(removeNA(v))
f6(v) = mean(removeNA(v)) # Make this an iterator
f7(v1, v2) = v1 + v2
f8(v1, v2) = v1 .> v2
f9(v, i) = v[i]

perf_test = Dict()

perf_test["sum(v): Vector with no NA's"] = () -> f1(v)
perf_test["sum(dv): DataVector with no NA's"] = () -> f1(dv)
perf_test["sum(removeNA(dv)): DataVector with no NA's"] = () -> f2(dv)
perf_test["sum(*removeNA(dv)): DataVector with no NA's"] = () -> f3(dv)

perf_test["sum(dvna): DataVector with NA's"] = () -> f4(dv)
perf_test["sum(removeNA(dvna)): DataVector with NA's"] = () -> f5(dv)
perf_test["sum(*removeNA(dvna)): DataVector with NA's"] = () -> f6(dv)

perf_test["mean(v): Vector with no NA's"] = () -> f4(v)
perf_test["mean(dv): DataVector with no NA's"] = () -> f4(dv)
perf_test["mean(removeNA(dv)): DataVector with no NA's"] = () -> f5(dv)
perf_test["mean(*removeNA(dv)): DataVector with no NA's"] = () -> f6(dv)

perf_test["mean(dvna): DataVector with NA's"] = () -> f4(dv)
perf_test["mean(removeNA(dvna)): DataVector with NA's"] = () -> f5(dv)
perf_test["mean(*removeNA(dvna)): DataVector with NA's"] = () -> f6(dv)

perf_test["v + 1.0 : Vector"] = () -> f7(v, 1.0)
perf_test["dv + 1.0 : DataVector with no NA's"] = () -> f7(dv, 1.0)
perf_test["dvna + 1.0 : DataVector with NA's"] = () -> f7(dvna, 1.0)

perf_test["v .> 1.0 : Vector"] = () -> f8(v, 1.0)
perf_test["dv .> 1.0 : DataVector with no NA's"] = () -> f8(dv, 1.0)
perf_test["dvna .> 1.0 : DataVector with NA's"] = () -> f8(dvna, 1.0)

perf_test["v[idxv] : Vector"] = () -> f9(v, idxv)
perf_test["dv[idxv] : DataVector and Vector indexing"] = () -> f9(dv, idxv)
perf_test["dv[idxdv] : DataVector and DataVector indexing"] = () -> f9(dv, idxdv)

for (name, f) in perf_test
    res = benchmark(f, "DataArray Operations", name, 10)
	# TODO: Keep permanent record
	printtable(res, header=false)
end


text:
# Kernels for Scalar,Array->Array

using Knet: broadcast_ops

function cuda01src(f, j=f, ex="$f(xi,yi)"; BLK=128, THR=128)
  sprint() do s
    for (T,F) in [("float","$(f)_32"),("double","$(f)_64")]
        print(s,
"""
__global__ void _$(F)_01(int n, $T xi, $T *y, $T *z) {
  int i = threadIdx.x + blockIdx.x * blockDim.x;
  while (i < n) {
    $T yi = y[i];
    z[i] = $ex;
    i += blockDim.x * gridDim.x;
  }
}
#ifdef __cplusplus
extern "C" {
#endif
  void $(F)_01(int n, $T xi, $T *y, $T *z) {
    _$(F)_01<<<$BLK,$THR>>>(n,xi,y,z);
  }
#ifdef __cplusplus
}
#endif
""")
    end
  end
end

for a in broadcast_ops
    if !isa(a,Tuple); a=(a,); end
    print(cuda01src(a...))
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule Unishox_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("Unishox")
JLLWrappers.@generate_main_file("Unishox", UUID("3d704607-676d-54ad-a58f-27926b53472b"))
end  # module Unishox_jll


text:
@testset "Simulated Annealing" begin
    Random.seed!(1)

    function f_s(x::Vector)
        (x[1] - 5.0)^4
    end
    options = Optim.Options(iterations=100_000)
    results = Optim.optimize(f_s, [0.0], SimulatedAnnealing(), options)
    @test norm(Optim.minimizer(results) - [5.0]) < 0.1

    function rosenbrock_s(x::Vector)
        (1.0 - x[1])^2 + 100.0 * (x[2] - x[1]^2)^2
    end
    options = Optim.Options(iterations=100_000)
    results = Optim.optimize(rosenbrock_s, [0.0, 0.0], SimulatedAnnealing(), options)
    @test norm(Optim.minimizer(results) - [1.0, 1.0]) < 0.1

    options = Optim.Options(iterations=10, show_trace=true, store_trace=true, extended_trace=true)
    results = Optim.optimize(rosenbrock_s, [0.0, 0.0], SimulatedAnnealing(), options)
end