text:
tc{N}(r1::NTuple{N}, r2::NTuple{N}) = all(map(x->tc(x...), [zip(r1,r2)...]))
tc{N}(r1::BitArray{N}, r2::Union(BitArray{N},Array{Bool,N})) = true
tc{T}(r1::T, r2::T) = true
tc(r1,r2) = false

function check_bitop(ret_type, func, args...)
    r1 = func(args...)
    r2 = func(map(x->(isa(x, BitArray) ? bitunpack(x) : x), args)...)
    @test isa(r1, ret_type)
    @test tc(r1, r2)
    @test isequal(r1, convert(ret_type, r2))
end

macro check_bit_operation(ex, ret_type)
    @assert Meta.isexpr(ex, :call)
    Expr(:call, :check_bitop, esc(ret_type), map(esc,ex.args)...)
end

let t0 = time()
    global timesofar
    function timesofar(str)
        return # no-op, comment to see timings
        t1 = time()
        println(str, ": ", t1-t0, " seconds")
        t0 = t1
    end
end

# vectors size
v1 = 260
# matrices size
n1, n2 = 17, 20
# arrays size
s1, s2, s3, s4 = 5, 8, 3, 4

allsizes = [((), BitArray{0}), ((v1,), BitVector),
            ((n1,n2), BitMatrix), ((s1,s2,s3,s4), BitArray{4})]

## Conversions ##

for (sz,T) in allsizes
    b1 = randbool!(falses(sz...))
    @test isequal(bitpack(bitunpack(b1)), b1)
    @test isequal(convert(Array{Float64,ndims(b1)}, b1),
                  convert(Array{Float64,ndims(b1)}, bitunpack(b1)))
    @test isequal(convert(AbstractArray{Float64,ndims(b1)}, b1),
                  convert(AbstractArray{Float64,ndims(b1)}, bitunpack(b1)))

    i1 = rand!(false:true, zeros(Bool, sz...))
    @test isequal(bitunpack(bitpack(i1)), i1)
end

timesofar("conversions")

## utility functions ##

b1 = randbool(v1)
@test isequal(fill!(b1, true), trues(size(b1)))
@test isequal(fill!(b1, false), falses(size(b1)))

for (sz,T) in allsizes
    @test isequal(bitunpack(trues(sz...)), ones(Bool, sz...))
    @test isequal(bitunpack(falses(sz...)), zeros(Bool, sz...))

    b1 = randbool!(falses(sz...))
    @test isa(b1, T)

    @check_bit_operation length(b1) Int
    @check_bit_operation ndims(b1)  Int
    @check_bit_operation size(b1)   (Int...)

    b2 = similar(b1)
    @check_bit_operation copy!(b2, b1) T
end

timesofar("utils")

## Indexing ##

# 0d
for (sz,T) in allsizes
    b1 = randbool!(falses(sz...))
    @check_bit_operation getindex(b1)         Bool
    @check_bit_operation setindex!(b1, true)  T
    @check_bit_operation setindex!(b1, false) T
    @check_bit_operation setindex!(b1, 1.0)  T
end

# linear
for (sz,T) in allsizes[2:end]
    l = *(sz...)
    b1 = randbool!(falses(sz...))
    for j = 1:l
        @check_bit_operation getindex(b1, j) Bool
    end
    @check_bit_operation getindex(b1, 100.0) Bool

    for j in [0, 1, 63, 64, 65, 127, 128, 129, 191, 192, 193, l-1, l]
        @check_bit_operation getindex(b1, 1:j)   BitVector
        @check_bit_operation getindex(b1, j+1:l) BitVector
    end
    for j in [1, 63, 64, 65, 127, 128, 129, div(l,2)]
        m1 = j:(l-j)
        @check_bit_operation getindex(b1, m1) BitVector
    end
    @check_bit_operation getindex(b1, 1.0:100.0) BitVector

    t1 = find(randbool(l))
    @check_bit_operation getindex(b1, t1)        BitVector
    @check_bit_operation getindex(b1, float(t1)) BitVector

    for j = 1:l
        x = randbool()
        @check_bit_operation setindex!(b1, x, j) T
    end

    x = randbool()
    @check_bit_operation setindex!(b1, x, 100.0) T
    y = rand(0.0:1.0)
    @check_bit_operation setindex!(b1, y, 100) T
    y = rand(0.0:1.0)
    @check_bit_operation setindex!(b1, y, 100.0) T

    for j in [1, 63, 64, 65, 127, 128, 129, 191, 192, 193, l-1]
        x = randbool()
        @check_bit_operation setindex!(b1, x, 1:j) T
        b2 = randbool(j)
        @check_bit_operation setindex!(b1, b2, 1:j) T
        x = randbool()
        @check_bit_operation setindex!(b1, x, j+1:l) T
        b2 = randbool(l-j)
        @check_bit_operation setindex!(b1, b2, j+1:l) T
    end
    for j in [1, 63, 64, 65, 127, 128, 129, div(l,2)]
        m1 = j:(l-j)
        x = randbool()
        @check_bit_operation setindex!(b1, x, m1) T
        b2 = randbool(length(m1))
        @check_bit_operation setindex!(b1, b2, m1) T
    end
    x = randbool()
    @check_bit_operation setindex!(b1, x, 1.0:100.0) T
    b2 = randbool(100)
    @check_bit_operation setindex!(b1, b2, 1.0:100.0) T

    y = rand(0.0:1.0)
    @check_bit_operation setindex!(b1, y, 1:100) T
    f2 = float(randbool(100))
    @check_bit_operation setindex!(b1, f2, 1:100) T
    f2 = float(randbool(100))
    @check_bit_operation setindex!(b1, f2, 1.0:100.0) T

    t1 = find(randbool(l))
    x = randbool()
    @check_bit_operation setindex!(b1, x, t1) T
    b2 = randbool(length(t1))
    @check_bit_operation setindex!(b1, b2, t1) T

    y = rand(0.0:1.0)
    @check_bit_operation setindex!(b1, y, t1) T
    f2 = float(randbool(length(t1)))
    @check_bit_operation setindex!(b1, f2, t1) T

    ft1 = float(t1)
    x = randbool()
    @check_bit_operation setindex!(b1, x, ft1) T
    b2 = randbool(length(t1))
    @check_bit_operation setindex!(b1, b2, ft1) T

    y = rand(0.0:1.0)
    @check_bit_operation setindex!(b1, y, ft1) T
    f2 = float(randbool(length(t1)))
    @check_bit_operation setindex!(b1, f2, ft1) T
end

# multidimensional

rand_m1m2() = rand(1:n1), rand(1:n2)

b1 = randbool(n1, n2)

m1, m2 = rand_m1m2()
b2 = randbool(m1, m2)
@check_bit_operation copy!(b1, b2) BitMatrix

function gen_getindex_data()
    m1, m2 = rand_m1m2()
    produce((m1, m2, Bool))
    m1, m2 = rand_m1m2()
    produce((m1, 1:m2, BitMatrix))
    m1, m2 = rand_m1m2()
    produce((m1, randperm(m2), BitMatrix))
    m1, m2 = rand_m1m2()
    produce((1:m1, m2, BitVector))
    m1, m2 = rand_m1m2()
    produce((1:m1, 1:m2, BitMatrix))
    m1, m2 = rand_m1m2()
    produce((1:m1, randperm(m2), BitMatrix))
    m1, m2 = rand_m1m2()
    produce((randperm(m1), m2, BitVector))
    m1, m2 = rand_m1m2()
    produce((randperm(m1), 1:m2, BitMatrix))
    m1, m2 = rand_m1m2()
    produce((randperm(m1), randperm(m2), BitMatrix))
end

for (k1, k2, T) in Task(gen_getindex_data)
    # println(typeof(k1), " ", typeof(k2), " ", T) # uncomment to debug
    @check_bit_operation getindex(b1, k1, k2) T
    @check_bit_operation getindex(b1, k1, k2, 1) T

    #@check_bit_operation getindex(b1, float(k1), k2) T
    #@check_bit_operation getindex(b1, k1, float(k2)) T

    @check_bit_operation getindex(b1, float(k1), float(k2)) T

    @check_bit_operation getindex(b1, k1, k2, 1.0) T
    #@check_bit_operation getindex(b1, float(k1), float(k2), 1.0) T
end

function gen_setindex_data()
    m1, m2 = rand_m1m2()
    produce((randbool(), m1, m2))
    m1, m2 = rand_m1m2()
    produce((randbool(), m1, 1:m2))
    produce((randbool(m2), m1, 1:m2))
    m1, m2 = rand_m1m2()
    produce((randbool(), m1, randperm(m2)))
    produce((randbool(m2), m1, randperm(m2)))
    m1, m2 = rand_m1m2()
    produce((randbool(), 1:m1, m2))
    produce((randbool(m1), 1:m1, m2))
    m1, m2 = rand_m1m2()
    produce((randbool(), 1:m1, 1:m2))
    produce((randbool(m1, m2), 1:m1, 1:m2))
    m1, m2 = rand_m1m2()
    produce((randbool(), 1:m1, randperm(m2)))
    produce((randbool(m1, m2), 1:m1, randperm(m2)))
    m1, m2 = rand_m1m2()
    produce((randbool(), randperm(m1), m2))
    produce((randbool(m1), randperm(m1), m2))
    m1, m2 = rand_m1m2()
    produce((randbool(), randperm(m1), 1:m2))
    produce((randbool(m1,m2), randperm(m1), 1:m2))
    m1, m2 = rand_m1m2()
    produce((randbool(), randperm(m1), randperm(m2)))
    produce((randbool(m1,m2), randperm(m1), randperm(m2)))
end

for (b2, k1, k2) in Task(gen_setindex_data)
    # println(typeof(b2), " ", typeof(k1), " ", typeof(k2)) # uncomment to debug
    @check_bit_operation setindex!(b1, b2, k1, k2) BitMatrix

    @check_bit_operation setindex!(b1, float(b2), k1, k2) BitMatrix
    @check_bit_operation setindex!(b1, b2, float(k1), k2) BitMatrix
    #@check_bit_operation setindex!(b1, b2, k1, float(k2)) BitMatrix

    #@check_bit_operation setindex!(b1, float(b2), float(k1), k2) BitMatrix
    #@check_bit_operation setindex!(b1, float(b2), k1, float(k2)) BitMatrix
    #@check_bit_operation setindex!(b1, b2, float(k1), float(k2)) BitMatrix

    @check_bit_operation setindex!(b1, float(b2), float(k1), float(k2)) BitMatrix
end

m1, m2 = rand_m1m2()
b2 = randbool(1, 1, m2)
@check_bit_operation setindex!(b1, b2, m1, 1:m2) BitMatrix
x = randbool()
b2 = randbool(1, m2, 1)
@check_bit_operation setindex!(b1, x, m1, 1:m2, 1)  BitMatrix
@check_bit_operation setindex!(b1, b2, m1, 1:m2, 1) BitMatrix

for p1 = [rand(1:v1) 1 63 64 65 191 192 193]
    for p2 = [rand(1:v1) 1 63 64 65 191 192 193]
        for n = 0 : min(v1 - p1 + 1, v1 - p2 + 1)
            b1 = randbool(v1)
            b2 = randbool(v1)
            @check_bit_operation copy!(b1, p1, b2, p2, n) BitVector
        end
    end
end

# logical indexing
b1 = randbool(n1, n2)
t1 = randbool(n1, n2)
@test isequal(bitunpack(b1[t1]), bitunpack(b1)[t1])
@test isequal(bitunpack(b1[t1]), bitunpack(b1)[bitunpack(t1)])

t1 = randbool(n1)
t2 = randbool(n2)
@test isequal(bitunpack(b1[t1, t2]), bitunpack(b1)[t1, t2])
@test isequal(bitunpack(b1[t1, t2]), bitunpack(b1)[bitunpack(t1), bitunpack(t2)])

b1 = randbool(n1, n2)
t1 = randbool(n1, n2)
@check_bit_operation setindex!(b1, true, t1) BitMatrix

t1 = randbool(n1, n2)
b2 = randbool(countnz(t1))
@check_bit_operation setindex!(b1, b2, t1) BitMatrix

m1 = rand(1:n1)
m2 = rand(1:n2)

t1 = randbool(n1)
b2 = randbool(countnz(t1), m2)
k2 = randperm(m2)
@check_bit_operation setindex!(b1, b2, t1, 1:m2)       BitMatrix
@check_bit_operation setindex!(b1, b2, t1, n2-m2+1:n2) BitMatrix
@check_bit_operation setindex!(b1, b2, t1, k2)         BitMatrix

t2 = randbool(n2)
b2 = randbool(m1, countnz(t2))
k1 = randperm(m1)
@check_bit_operation setindex!(b1, b2, 1:m1, t2)       BitMatrix
@check_bit_operation setindex!(b1, b2, n1-m1+1:n1, t2) BitMatrix
@check_bit_operation setindex!(b1, b2, k1, t2)         BitMatrix

timesofar("indexing")

## Dequeue functionality ##

b1 = BitArray(0)
i1 = Bool[]
for m = 1 : v1
    x = randbool()
    push!(b1, x)
    push!(i1, x)
    @test isequal(bitunpack(b1), i1)
end

for m1 = 0 : v1
    for m2 = [0, 1, 63, 64, 65, 127, 128, 129]
        b1 = randbool(m1)
        b2 = randbool(m2)
        i1 = bitunpack(b1)
        i2 = bitunpack(b2)
        @test isequal(bitunpack(append!(b1, b2)), append!(i1, i2))
    end
end

for m1 = 0 : v1
    for m2 = [0, 1, 63, 64, 65, 127, 128, 129]
        b1 = randbool(m1)
        b2 = randbool(m2)
        i1 = bitunpack(b1)
        i2 = bitunpack(b2)
        @test isequal(bitunpack(prepend!(b1, b2)), prepend!(i1, i2))
    end
end

b1 = randbool(v1)
i1 = bitunpack(b1)
for m = 1 : v1
    jb = pop!(b1)
    ji = pop!(i1)
    @test jb == ji
    @test isequal(bitunpack(b1), i1)
end
@test length(b1) == 0

b1 = BitArray(0)
i1 = Bool[]
for m = 1 : v1
    x = randbool()
    unshift!(b1, x)
    unshift!(i1, x)
    @test isequal(bitunpack(b1), i1)
end

b1 = randbool(v1)
i1 = bitunpack(b1)
for m = 1 : v1
    jb = shift!(b1)
    ji = shift!(i1)
    @test jb == ji
    @test isequal(bitunpack(b1), i1)
end
@test length(b1) == 0

b1 = BitArray(0)
i1 = bitunpack(b1)
for m = 1 : v1
    j = rand(1:m)
    x = randbool()
    insert!(b1, j, x)
    insert!(i1, j, x)
    @test isequal(bitunpack(b1), i1)
end

b1 = randbool(v1)
i1 = bitunpack(b1)
for j in [63, 64, 65, 127, 128, 129, 191, 192, 193]
    x = rand(0:1)
    insert!(b1, j, x)
    insert!(i1, j, x)
    @test isequal(bitunpack(b1), i1)
end

b1 = randbool(v1)
i1 = bitunpack(b1)
for m = v1 : -1 : 1
    j = rand(1:m)
    b = splice!(b1, j)
    i = splice!(i1, j)
    @test isequal(bitunpack(b1), i1)
    @test b == i
end
@test length(b1) == 0

b1 = randbool(v1)
i1 = bitunpack(b1)
for m = v1 : -1 : 1
    j = rand(1:m)
    deleteat!(b1, j)
    deleteat!(i1, j)
    @test isequal(bitunpack(b1), i1)
end
@test length(b1) == 0

b1 = randbool(v1)
i1 = bitunpack(b1)
for j in [63, 64, 65, 127, 128, 129, 191, 192, 193]
    b = splice!(b1, j)
    i = splice!(i1, j)
    @test isequal(bitunpack(b1), i1)
    @test b == i
end

b1 = randbool(v1)
i1 = bitunpack(b1)
for j in [63, 64, 65, 127, 128, 129, 191, 192, 193]
    deleteat!(b1, j)
    deleteat!(i1, j)
    @test isequal(bitunpack(b1), i1)
end

b1 = randbool(v1)
i1 = bitunpack(b1)
for m1 = 1 : v1
    for m2 = m1 : v1
        b2 = copy(b1)
        i2 = copy(i1)
        b = splice!(b2, m1:m2)
        i = splice!(i2, m1:m2)
        @test isequal(bitunpack(b2), i2)
        @test b == i
    end
end

b1 = randbool(v1)
i1 = bitunpack(b1)
for m1 = 1 : v1
    for m2 = m1 : v1
        b2 = copy(b1)
        i2 = copy(i1)
        deleteat!(b2, m1:m2)
        deleteat!(i2, m1:m2)
        @test isequal(bitunpack(b2), i2)
    end
end

b1 = randbool(v1)
i1 = bitunpack(b1)
for m1 = 1 : v1 + 1
    for m2 = m1 - 1 : v1
        for v2 = [0, 1, 63, 64, 65, 127, 128, 129, 191, 192, 193, rand(1:v1)]
            b2 = copy(b1)
            i2 = copy(i1)
            b3 = randbool(v2)
            i3 = bitunpack(b3)
            b = splice!(b2, m1:m2, b3)
            i = splice!(i2, m1:m2, i3)
            @test isequal(bitunpack(b2), i2)
            @test b == i
        end
    end
end

b1 = randbool(v1)
i1 = bitunpack(b1)
for m1 = 1 : v1
    for v2 = [0, 1, 63, 64, 65, 127, 128, 129, 191, 192, 193, rand(1:v1)]
        b2 = copy(b1)
        i2 = copy(i1)
        b3 = randbool(v2)
        i3 = bitunpack(b3)
        b = splice!(b2, m1, b3)
        i = splice!(i2, m1, i3)
        @test isequal(bitunpack(b2), i2)
        @test b == i
    end
end

b1 = randbool(v1)
i1 = bitunpack(b1)
for m1 = 1 : v1 - 1
    for m2 = m1 + 1 : v1
        locs = randbool(m2-m1+1)
        m = [m1:m2...][locs]
        b2 = copy(b1)
        i2 = copy(i1)
        deleteat!(b2, m)
        deleteat!(i2, m)
        @test isequal(bitunpack(b2), i2)
    end
end

b1 = randbool(v1)
i1 = bitunpack(b1)
empty!(b1)
empty!(i1)
@test isequal(bitunpack(b1), i1)

timesofar("dequeue")

## Unary operators ##

b1 = randbool(n1, n2)
@check_bit_operation (~)(b1)  BitMatrix
@check_bit_operation (!)(b1)  BitMatrix
@check_bit_operation (-)(b1)  Matrix{Int}
@check_bit_operation sign(b1) BitMatrix
@check_bit_operation real(b1) BitMatrix
@check_bit_operation imag(b1) BitMatrix
@check_bit_operation conj(b1) BitMatrix

b0 = falses(0)
@check_bit_operation (~)(b0)  BitVector
@check_bit_operation (!)(b0)  BitVector
@check_bit_operation (-)(b0)  Vector{Int}
@check_bit_operation sign(b0) BitVector

timesofar("unary arithmetic")

## Binary arithmetic operators ##

# Matrix{Bool}/Matrix{Bool}

b1 = randbool(n1, n2)
b2 = randbool(n1, n2)
@check_bit_operation (&)(b1, b2)  BitMatrix
@check_bit_operation (|)(b1, b2)  BitMatrix
@check_bit_operation ($)(b1, b2)  BitMatrix
@check_bit_operation (+)(b1, b2)  Matrix{Int}
@check_bit_operation (-)(b1, b2)  Matrix{Int}
@check_bit_operation (.*)(b1, b2) BitMatrix
@check_bit_operation (./)(b1, b2) Matrix{Float32}
@check_bit_operation (.^)(b1, b2) BitMatrix

b2 = trues(n1, n2)
@check_bit_operation div(b1, b2) BitMatrix
@check_bit_operation mod(b1, b2) BitMatrix

while true
    global b1
    b1 = randbool(n1, n1)
    if abs(det(float64(b1))) > 1e-6
        break
    end
end
b2 = randbool(n1, n1)

@check_bit_operation (*)(b1, b2) Matrix{Int}
@check_bit_operation (/)(b1, b1) Matrix{Float32}
@check_bit_operation (\)(b1, b1) Matrix{Float32}

b0 = falses(0)
@check_bit_operation (&)(b0, b0)  BitVector
@check_bit_operation (|)(b0, b0)  BitVector
@check_bit_operation ($)(b0, b0)  BitVector
@check_bit_operation (.*)(b0, b0) BitVector
@check_bit_operation (*)(b0, b0') Matrix{Int}

# Matrix{Bool}/Matrix{Int}
b1 = randbool(n1, n2)
i2 = rand(1:10, n1, n2)
@check_bit_operation (&)(b1, i2)  Matrix{Int}
@check_bit_operation (|)(b1, i2)  Matrix{Int}
@check_bit_operation ($)(b1, i2)  Matrix{Int}
@check_bit_operation (+)(b1, i2)  Matrix{Int}
@check_bit_operation (-)(b1, i2)  Matrix{Int}
@check_bit_operation (.*)(b1, i2) Matrix{Int}
@check_bit_operation (./)(b1, i2) Matrix{Float64}
@check_bit_operation (.^)(b1, i2) BitMatrix
@check_bit_operation div(b1, i2)  Matrix{Int}
@check_bit_operation mod(b1, i2)  Matrix{Int}

# Matrix{Bool}/Matrix{Float64}
b1 = randbool(n1, n2)
f2 = 1.0 .+ rand(n1, n2)
@check_bit_operation (.*)(b1, f2) Matrix{Float64}
@check_bit_operation (./)(b1, f2) Matrix{Float64}
@check_bit_operation (.^)(b1, f2) Matrix{Float64}
@check_bit_operation div(b1, f2)  Matrix{Float64}
@check_bit_operation mod(b1, f2)  Matrix{Float64}

# Number/Matrix
b2 = randbool(n1, n2)
i1 = rand(1:10)
u1 = uint8(i1)
f1 = float64(i1)
ci1 = complex(i1)
cu1 = complex(u1)
cf1 = complex(f1)

@check_bit_operation (&)(i1, b2)  Matrix{Int}
@check_bit_operation (|)(i1, b2)  Matrix{Int}
@check_bit_operation ($)(i1, b2)  Matrix{Int}
@check_bit_operation (.+)(i1, b2)  Matrix{Int}
@check_bit_operation (.-)(i1, b2)  Matrix{Int}
@check_bit_operation (.*)(i1, b2) Matrix{Int}

@check_bit_operation (&)(u1, b2)  Matrix{Uint8}
@check_bit_operation (|)(u1, b2)  Matrix{Uint8}
@check_bit_operation ($)(u1, b2)  Matrix{Uint8}
@check_bit_operation (.+)(u1, b2)  Matrix{Uint8}
@check_bit_operation (.-)(u1, b2)  Matrix{Uint8}
@check_bit_operation (.*)(u1, b2) Matrix{Uint8}

for (x1,t1) = {(f1, Float64),
               (ci1, Complex{Int}),
               (cu1, Complex{Uint8}),
               (cf1, Complex128)}
    @check_bit_operation (.+)(x1, b2)  Matrix{t1}
    @check_bit_operation (.-)(x1, b2)  Matrix{t1}
    @check_bit_operation (.*)(x1, b2) Matrix{t1}
end

b2 = trues(n1, n2)
@check_bit_operation (./)(true, b2)  Matrix{Float32}
@check_bit_operation div(true, b2)   BitMatrix
@check_bit_operation mod(true, b2)   BitMatrix
@check_bit_operation (./)(false, b2) Matrix{Float32}
@check_bit_operation div(false, b2)  BitMatrix
@check_bit_operation mod(false, b2)  BitMatrix

@check_bit_operation (./)(i1, b2) Matrix{Float64}
@check_bit_operation div(i1, b2)  Matrix{Int}
@check_bit_operation mod(i1, b2)  Matrix{Int}

@check_bit_operation (./)(u1, b2) Matrix{Float32}
@check_bit_operation div(u1, b2)  Matrix{Uint8}
@check_bit_operation mod(u1, b2)  Matrix{Uint8}

@check_bit_operation (./)(f1, b2) Matrix{Float64}
@check_bit_operation div(f1, b2)  Matrix{Float64}
@check_bit_operation mod(f1, b2)  Matrix{Float64}

@check_bit_operation (./)(ci1, b2) Matrix{Complex128}
@check_bit_operation (./)(cu1, b2) Matrix{Complex64}
@check_bit_operation (./)(cf1, b2) Matrix{Complex128}

b2 = randbool(n1, n2)
@check_bit_operation (.^)(false, b2) BitMatrix
@check_bit_operation (.^)(true, b2)  BitMatrix
@check_bit_operation (.^)(0x0, b2)   Matrix{Uint8}
@check_bit_operation (.^)(0x1, b2)   Matrix{Uint8}
@check_bit_operation (.^)(-1, b2)    Matrix{Int}
@check_bit_operation (.^)(0, b2)     Matrix{Int}
@check_bit_operation (.^)(1, b2)     Matrix{Int}
@check_bit_operation (.^)(0.0, b2)   Matrix{Float64}
@check_bit_operation (.^)(1.0, b2)   Matrix{Float64}
@check_bit_operation (.^)(0.0im, b2) Matrix{Complex128}
@check_bit_operation (.^)(1.0im, b2) Matrix{Complex128}
@check_bit_operation (.^)(0im, b2)   Matrix{Complex{Int}}
@check_bit_operation (.^)(1im, b2)   Matrix{Complex{Int}}
@check_bit_operation (.^)(0x0im, b2) Matrix{Complex{Uint8}}
@check_bit_operation (.^)(0x1im, b2) Matrix{Complex{Uint8}}

# Matrix/Number
b1 = randbool(n1, n2)
i2 = rand(1:10)
u2 = uint8(i2)
f2 = float64(i2)
ci2 = complex(i2)
cu2 = complex(u2)
cf2 = complex(f2)

@check_bit_operation (&)(b1, true)   BitMatrix
@check_bit_operation (&)(b1, false)  BitMatrix
@check_bit_operation (|)(b1, true)   BitMatrix
@check_bit_operation (|)(b1, false)  BitMatrix
@check_bit_operation ($)(b1, true)   BitMatrix
@check_bit_operation ($)(b1, false)  BitMatrix
@check_bit_operation (.+)(b1, true)   Matrix{Int}
@check_bit_operation (.+)(b1, false)  Matrix{Int}
@check_bit_operation (.-)(b1, true)   Matrix{Int}
@check_bit_operation (.-)(b1, false)  Matrix{Int}
@check_bit_operation (.*)(b1, true)  BitMatrix
@check_bit_operation (.*)(b1, false) BitMatrix
@check_bit_operation (./)(b1, true)  Matrix{Float32}
@check_bit_operation (./)(b1, false) Matrix{Float32}
@check_bit_operation div(b1, true)   BitMatrix
@check_bit_operation mod(b1, true)   BitMatrix

@check_bit_operation (&)(b1, i2)  Matrix{Int}
@check_bit_operation (|)(b1, i2)  Matrix{Int}
@check_bit_operation ($)(b1, i2)  Matrix{Int}
@check_bit_operation (.+)(b1, i2)  Matrix{Int}
@check_bit_operation (.-)(b1, i2)  Matrix{Int}
@check_bit_operation (.*)(b1, i2) Matrix{Int}
@check_bit_operation (./)(b1, i2) Matrix{Float64}
@check_bit_operation div(b1, i2)  Matrix{Int}
@check_bit_operation mod(b1, i2)  Matrix{Int}

@check_bit_operation (&)(b1, u2)  Matrix{Uint8}
@check_bit_operation (|)(b1, u2)  Matrix{Uint8}
@check_bit_operation ($)(b1, u2)  Matrix{Uint8}
@check_bit_operation (.+)(b1, u2)  Matrix{Uint8}
@check_bit_operation (.-)(b1, u2)  Matrix{Uint8}
@check_bit_operation (.*)(b1, u2) Matrix{Uint8}
@check_bit_operation (./)(b1, u2) Matrix{Float32}
@check_bit_operation div(b1, u2)  Matrix{Uint8}
@check_bit_operation mod(b1, u2)  Matrix{Uint8}

@check_bit_operation (.+)(b1, f2)  Matrix{Float64}
@check_bit_operation (.-)(b1, f2)  Matrix{Float64}
@check_bit_operation (.*)(b1, f2) Matrix{Float64}
@check_bit_operation (./)(b1, f2) Matrix{Float64}
@check_bit_operation div(b1, f2)  Matrix{Float64}
@check_bit_operation mod(b1, f2)  Matrix{Float64}

@check_bit_operation (.+)(b1, ci2)  Matrix{Complex{Int}}
@check_bit_operation (.-)(b1, ci2)  Matrix{Complex{Int}}
@check_bit_operation (.*)(b1, ci2) Matrix{Complex{Int}}
@check_bit_operation (./)(b1, ci2) Matrix{Complex128}

@check_bit_operation (.+)(b1, cu2)  Matrix{Complex{Uint8}}
@check_bit_operation (.-)(b1, cu2)  Matrix{Complex{Uint8}}
@check_bit_operation (.*)(b1, cu2) Matrix{Complex{Uint8}}
@check_bit_operation (./)(b1, cu2) Matrix{Complex64}

@check_bit_operation (.+)(b1, cf2)  Matrix{Complex128}
@check_bit_operation (.-)(b1, cf2)  Matrix{Complex128}
@check_bit_operation (.*)(b1, cf2) Matrix{Complex128}
@check_bit_operation (./)(b1, cf2) Matrix{Complex128}

@check_bit_operation (.^)(b1, false) BitMatrix
@check_bit_operation (.^)(b1, true)  BitMatrix
@check_bit_operation (.^)(b1, 0x0)   BitMatrix
@check_bit_operation (.^)(b1, 0x1)   BitMatrix
@check_bit_operation (.^)(b1, 0)     BitMatrix
@check_bit_operation (.^)(b1, 1)     BitMatrix
@check_bit_operation (.^)(b1, -1.0)  Matrix{Float64}
@check_bit_operation (.^)(b1, 0.0)   Matrix{Float64}
@check_bit_operation (.^)(b1, 1.0)   Matrix{Float64}
@check_bit_operation (.^)(b1, 0.0im) Matrix{Complex128}
@check_bit_operation (.^)(b1, 0x0im) Matrix{Complex64}
@check_bit_operation (.^)(b1, 0im)   Matrix{Complex128}

b1 = trues(n1, n2)
@check_bit_operation (.^)(b1, -1.0im) Matrix{Complex128}
@check_bit_operation (.^)(b1, 1.0im)  Matrix{Complex128}
@check_bit_operation (.^)(b1, -1im)   Matrix{Complex128}
@check_bit_operation (.^)(b1, 1im)    Matrix{Complex128}
@check_bit_operation (.^)(b1, 0x1im)  Matrix{Complex64}

timesofar("binary arithmetic")

## Binary comparison operators ##

b1 = randbool(n1, n2)
b2 = randbool(n1, n2)
@check_bit_operation (.==)(b1, b2) BitMatrix
@check_bit_operation (.!=)(b1, b2) BitMatrix
@check_bit_operation (.<)(b1, b2) BitMatrix
@check_bit_operation (.<=)(b1, b2) BitMatrix

timesofar("binary comparison")

## Data movement ##

b1 = randbool(s1, s2, s3, s4)
for d = 1 : 4
    j = rand(1:size(b1, d))
    #for j = 1 : size(b1, d)
        @check_bit_operation slicedim(b1, d, j) BitArray{4}
    #end
    @check_bit_operation flipdim(b1, d) BitArray{4}
end

b1 = randbool(n1, n2)
for k = 1 : 4
    @check_bit_operation rotl90(b1, k) BitMatrix
end

for m = 0 : v1
    b1 = randbool(m)
    @check_bit_operation reverse(b1) BitVector
end

b1 = randbool(v1)
for m = [rand(1:v1)-1 0 1 63 64 65 191 192 193 v1-1]
    @test isequal(b1 << m, [ b1[m+1:end]; falses(m) ])
    @test isequal(b1 >>> m, [ falses(m); b1[1:end-m] ])
    @test isequal(rol(b1, m), [ b1[m+1:end]; b1[1:m] ])
    @test isequal(ror(b1, m), [ b1[end-m+1:end]; b1[1:end-m] ])
    @test isequal(ror(b1, m), rol(b1, -m))
    @test isequal(rol(b1, m), ror(b1, -m))
end

timesofar("datamove")

## countnz & find ##

b1 = randbool(v1)
@check_bit_operation countnz(b1) Int

@check_bit_operation findfirst(b1)         Int
@check_bit_operation findfirst(trues(v1))  Int
@check_bit_operation findfirst(falses(v1)) Int

@check_bit_operation findfirst(b1, true)  Int
@check_bit_operation findfirst(b1, false) Int
@check_bit_operation findfirst(b1, 3)     Int

@check_bit_operation findfirst(x->x, b1)     Int
@check_bit_operation findfirst(x->!x, b1)    Int
@check_bit_operation findfirst(x->true, b1)  Int
@check_bit_operation findfirst(x->false, b1) Int

@check_bit_operation find(b1) Vector{Int}

b1 = trues(v1)
for i = 0:v1-1
    @test findfirst(b1 >> i) == i+1
    @test Base.findfirstnot(~(b1 >> i)) == i+1
end

for i = 3:v1-1
    for j = 2:i
        submask = b1 << (v1-j+1)
        @test findnext((b1 >> i) | submask,j) == i+1
        @test Base.findnextnot((~(b1 >> i)) $ submask,j) == i+1
    end
end

b1 = randbool(n1, n2)
@check_bit_operation findnz(b1) (Vector{Int}, Vector{Int}, BitArray)

timesofar("nnz&find")

## Reductions ##

b1 = randbool(s1, s2, s3, s4)
m1 = 1
m2 = 3
@check_bit_operation maximum(b1, (m1, m2)) BitArray{4}
@check_bit_operation minimum(b1, (m1, m2)) BitArray{4}
@check_bit_operation sum(b1, (m1, m2)) Array{Int,4}

@check_bit_operation maximum(b1) Bool
@check_bit_operation minimum(b1) Bool
@check_bit_operation any(b1) Bool
@check_bit_operation all(b1) Bool
@check_bit_operation sum(b1) Int

b0 = falses(0)
@check_bit_operation any(b0) Bool
@check_bit_operation all(b0) Bool
@check_bit_operation sum(b0) Int

timesofar("reductions")

## map over bitarrays ##

# TODO (not implemented)

## Filter ##

# TODO

## Transpose ##

b1 = randbool(v1)
@check_bit_operation transpose(b1) BitMatrix

for m1 = 0 : n1
    for m2 = 0 : n2
        b1 = randbool(m1, m2)
        @check_bit_operation transpose(b1) BitMatrix
    end
end

timesofar("transpose")

## Permutedims ##

b1 = randbool(s1, s2, s3, s4)
p = randperm(4)
@check_bit_operation permutedims(b1, p) BitArray{4}
@check_bit_operation permutedims(b1, tuple(p...)) BitArray{4}

timesofar("permutedims")

## Concatenation ##

b1 = randbool(v1)
b2 = randbool(v1)
@check_bit_operation hcat(b1, b2) BitMatrix
for m = 1 : v1 - 1
    @check_bit_operation vcat(b1[1:m], b1[m+1:end]) BitVector
end

b1 = randbool(n1, n2)
b2 = randbool(n1)
b3 = randbool(n1, n2)
b4 = randbool(1, n2)
@check_bit_operation hcat(b1, b2, b3) BitMatrix
@check_bit_operation vcat(b1, b4, b3) BitMatrix

b1 = randbool(s1, s2, s3, s4)
b2 = randbool(s1, s3, s3, s4)
b3 = randbool(s1, s2, s3, s1)
@check_bit_operation cat(2, b1, b2) BitArray{4}
@check_bit_operation cat(4, b1, b3) BitArray{4}
@check_bit_operation cat(6, b1, b1) BitArray{6}

b1 = randbool(1, v1, 1)
@check_bit_operation cat(2, 0, b1, 1, 1, b1) Array{Int,3}
@check_bit_operation cat(2, 3, b1, 4, 5, b1) Array{Int,3}
@check_bit_operation cat(2, false, b1, true, true, b1) BitArray{3}

b1 = randbool(n1, n2)
for m1 = 1 : n1 - 1
    for m2 = 1 : n2 - 1
        @test isequal([b1[1:m1,1:m2] b1[1:m1,m2+1:end]; b1[m1+1:end,1:m2] b1[m1+1:end,m2+1:end]], b1)
    end
end

timesofar("cat")

# Linear algebra

b1 = randbool(v1)
b2 = randbool(v1)
@check_bit_operation dot(b1, b2) Int

b1 = randbool(n1, n2)
for k = -max(n1,n2) : max(n1,n2)
    @check_bit_operation tril(b1, k) BitMatrix
    @check_bit_operation triu(b1, k) BitMatrix
end

b1 = randbool(n1, n1)
@check_bit_operation istril(b1) Bool
b1 = randbool(n1, n2)
@check_bit_operation istril(b1) Bool
b1 = randbool(n2, n1)
@check_bit_operation istril(b1) Bool

b1 = tril(randbool(n1, n1))
@check_bit_operation istril(b1) Bool
b1 = tril(randbool(n1, n2))
@check_bit_operation istril(b1) Bool
b1 = tril(randbool(n2, n1))
@check_bit_operation istril(b1) Bool

b1 = randbool(n1, n1)
@check_bit_operation istriu(b1) Bool
b1 = randbool(n1, n2)
@check_bit_operation istriu(b1) Bool
b1 = randbool(n2, n1)
@check_bit_operation istriu(b1) Bool

b1 = triu(randbool(n1, n1))
@check_bit_operation istriu(b1) Bool
b1 = triu(randbool(n1, n2))
@check_bit_operation istriu(b1) Bool
b1 = triu(randbool(n2, n1))
@check_bit_operation istriu(b1) Bool

b1 = randbool(n1)
b2 = randbool(n2)
@check_bit_operation kron(b1, b2) BitVector

b1 = randbool(s1, s2)
b2 = randbool(s3, s4)
@check_bit_operation kron(b1, b2) BitMatrix

#b1 = randbool(v1)
#@check_bit_operation diff(b1) Vector{Int}
#b1 = randbool(n1, n2)
#@check_bit_operation diff(b1) Vector{Int}

timesofar("linalg")

max_stars_repo_name: jakebolewski/julia

max_issues_repo_name: jakebolewski/julia

max_forks_repo_name: jakebolewski/julia


text:
module FluorescentSeries
using AxisArrays, ImageAxes

include("core.jl")
include("algorithms.jl")

export
FluorescentSerie,
deltaFF

end # module

ext: jl


text:
include("custom_graph.jl")

"""
    colorGraph(g)

Find a coloring of the graph g such that no two vertices connected
by an edge have the same color.
"""
function colorGraph(G::CGraph)
    colornumber = 0
    V = num_vertices(G)
    colors = zeros(Int64,V)

    while (V > 0)
        x = max_degree_vertex(G)
        colornumber = colornumber + 1
        colors[x] = colornumber
        nn = non_neighbors(G,x)
        while (length(nn) > 0)
            maxcn = -1
            ydegree = -1
            for z in nn
                cn = length_common_neighbor(G,z,x)
                if (cn > maxcn) || (cn == maxcn && vertex_degree(G,z) < ydegree)
                    y = z
                    ydegree = vertex_degree(G,y)
                    maxcn = cn
                end
            end
            if (maxcn == 0)
                y = max_degree_vertex(G,nn)
            end
            colors[y] = colornumber
            contract!(G,y,x)
            nn = non_neighbors(G,x)
        end
        rem_vertex!(G,x)
        V = num_vertices(G)
    end
    return colors
end

"""
    num_vertices(G)

Find the total number of vertices present in graph G.
"""
function num_vertices(G::CGraph)
    return length(vertices(G))
end

"""
    max_degree_vertex(G, nn)

Find the vertex in the group nn of vertices belonging to the
graph G which has the highest degree.
"""
function max_degree_vertex(G::CGraph,nn::Array{Int64,1})
    max_degree = -1
    max_degree_vertex = -1
    for v in nn
        v_degree = length(neighbors(G,v))
        if v_degree > max_degree
            max_degree = v_degree
            max_degree_vertex = v
        end
    end
    return max_degree_vertex
end

"""
    max_degree_vertex(G)

Find the vertex in graph with highest degree.
"""
function max_degree_vertex(G::CGraph)
    max_degree = -1
    max_degree_vertex = -1
    for v in vertices(G)
        v_degree = length(neighbors(G,v))
        if v_degree > max_degree
            max_degree = v_degree
            max_degree_vertex = v
        end
    end
    return max_degree_vertex
end

"""
    non_neighbors(G,x)

Find the set of vertices belonging to the graph G which do
not share an edge with the vertex x.
"""
function non_neighbors(G::CGraph, x::Int64)
    num_non_neighbors = num_vertices(G) - vertex_degree(G,x) - 1
    nn = zeros(Int64,num_non_neighbors)
    index = 1

    for i in vertices(G)
        if ((indexin([i], neighbors(G,x)))[1] == nothing) && (i != x)
            nn[index] = i
            index += 1
        end
    end
    return nn
end

"""
    length_common_neighbor(G,z,x)

Find the number of vertices that share an edge with both the
vertices z and x belonging to the graph G.
"""
function length_common_neighbor(G::CGraph,z::Int64, x::Int64)
    z_neighbors = neighbors(G,z)
    x_neighbors = neighbors(G,x)
    common_vertices = indexin(z_neighbors, x_neighbors)
    num_common_vertices = 0
    for i in common_vertices
        if i != nothing
            num_common_vertices += 1
        end
    end
    return num_common_vertices
end

"""
    vertex_degree(G,z)

Find the degree of the vertex z which belongs to the graph G.
"""
function vertex_degree(G::CGraph,z::Int64)
    return length(neighbors(G,z))
end

"""
    contract!(G,y,x)

Contract the vertex y to x, both of which belong to graph G, that is
delete vertex y and join x with the neighbors of y if they are not
already connected with an edge.
"""
function contract!(G::CGraph,y::Int64, x::Int64)
    for v in neighbors(G,y)
        if has_edge(G,v,x) == false
            add_edge!(G,v,x)
        end
    end
    rem_vertex!(G,y)

end


text:
module Econ890

using DocStringExtensions, Interpolations, Roots

export AbstractUtility, UtilityCRRA, UtilityLog
export euler_dev, utility, marg_utility, inv_utility, inv_marg_utility, c_growth

# Policy function iteration
export Model, init_test_model, euler_dev, solve

include("utility.jl");
include("vfi.jl");

end # module


text:
module GitLab

import ..GitForge: endpoint, into, postprocessor

using ..GitForge
using ..GitForge:
    @json,
    AStr,
    DoNothing,
    DoSomething,
    Endpoint,
    Forge,
    JSON,
    OnRateLimit,
    RateLimiter,
    HEADERS,
    ORL_THROW,
    @not_implemented

using Dates
using HTTP
using JSON2

export GitLabAPI, NoToken, OAuth2Token, PersonalAccessToken

const DEFAULT_URL = "https://gitlab.com/api/v4"
const JSON_OPTS = (
    dateformat=dateformat"y-m-d",
    read_datetimeformats=[dateformat"y-m-dTH:M:S.sZ", dateformat"y-m-dTH:M:S.s+ss:ss", dateformat"y-m-dTH:M:S.s-ss:ss"],
    write_datetimeformat=dateformat"y-m-dTH:M:S.sZ",
)

abstract type AbstractToken end

"""
    NoToken()

Represents no authentication.
Only public data will be available.
"""
struct NoToken <: AbstractToken end

"""
    OAuth2Token(token::$AStr)

An [OAuth2 bearer token](https://docs.gitlab.com/ce/api/#oauth2-tokens).
"""
struct OAuth2Token <: AbstractToken
    token::String
end

"""
    PersonalAccessToken(token::$AStr)

A [private access token](https://docs.gitlab.com/ce/api/#personal-access-tokens).
"""
struct PersonalAccessToken <: AbstractToken
    token::String
end

auth_headers(::NoToken) = []
auth_headers(t::OAuth2Token) = ["Authorization" => "Bearer $(t.token)"]
auth_headers(t::PersonalAccessToken) = ["Private-Token" => t.token]

"""
    GitLabAPI(;
        token::AbstractToken=NoToken(),
        url::$AStr="$DEFAULT_URL",
        has_rate_limits::Bool=true,
        on_rate_limit::OnRateLimit=ORL_THROW,
    )

Create a GitLab API client.

## Keywords
- `token::AbstractToken=NoToken()`: Authorization token (or lack thereof).
- `url::$AStr"$DEFAULT_URL"`: Base URL of the target GitLab instance.
- `has_rate_limits::Bool=true`: Whether or not the GitLab server has rate limits.
- `on_rate_limit::OnRateLimit=ORL_THROW`: Behaviour on exceeded rate limits.
"""
struct GitLabAPI <: Forge
    token::AbstractToken
    url::String
    hasrl::Bool
    orl::OnRateLimit
    rl::RateLimiter

    function GitLabAPI(;
        token::AbstractToken=NoToken(),
        url::AStr=DEFAULT_URL,
        has_rate_limits::Bool=true,
        on_rate_limit::OnRateLimit=ORL_THROW,
    )
        return new(token, url, has_rate_limits, on_rate_limit, RateLimiter())
    end
end

GitForge.base_url(g::GitLabAPI) = g.url
GitForge.request_headers(g::GitLabAPI, ::Function) = [HEADERS; auth_headers(g.token)]
GitForge.postprocessor(::GitLabAPI, ::Function) = JSON()
GitForge.has_rate_limits(g::GitLabAPI, ::Function) = g.hasrl
GitForge.rate_limit_check(g::GitLabAPI, ::Function) = GitForge.rate_limit_check(g.rl)
GitForge.on_rate_limit(g::GitLabAPI, ::Function) = g.orl
GitForge.rate_limit_wait(g::GitLabAPI, ::Function) = GitForge.rate_limit_wait(g.rl)
GitForge.rate_limit_period(g::GitLabAPI, ::Function) = GitForge.rate_limit_period(g.rl)
GitForge.rate_limit_update!(g::GitLabAPI, ::Function, r::HTTP.Response) =
    GitForge.rate_limit_update!(g.rl, r)

include("users.jl")
include("projects.jl")
include("merge_requests.jl")
include("groups.jl")
include("commits.jl")
include("branches.jl")
include("tags.jl")
include("notes.jl")
include("pipeline_schedules.jl")

encode(owner::AStr, repo::AStr) = HTTP.escapeuri("$owner/$repo")
encode(org::AStr) = HTTP.escapeuri(org)
encode(owner::AStr, subgroup::AStr, repo::AStr) = HTTP.escapeuri("$owner/$subgroup/$repo")

function ismember(r::HTTP.Response)
    r.status == 404 && return false
    m = JSON2.read(IOBuffer(r.body), Member)
    return m.access_level !== nothing && m.access_level >= 30
end

end


text:
# This file contains the inversion procedures for LAI (and ALIA)
# calculation. The default is the ellipsoidal optimization method.

# Constants

"Ring width for zenith57."
const RING_WIDTH = 5 / 180 * π
"""Number of grouped consecutive ρ² rings for Miller's approach
 (each typically with 4 or 8 pixels)."""
const MILLER_GROUPS = 10
"Minimum number of rings for PolarImages with slope."
const NRINGS_MIN_SLOPE = 25#10
"""Start view angle for use in Lang's regression method from Weiss
et al 2004 paragraph 2.2.2.2."""
const LANG_START = 25/180*pi
"""Stop view angle for use in Lang's regression method from Weiss
et al 2004 paragraph 2.2.2.2."""
const LANG_END = 65/180*pi
"Maximum viewing angle used."
const θMAX = π/2
"Number of points parameters for Lookup Table method as in Weiss et al 2004."
const LUT_POINTS = 10_000 # number of points in LUT
"Sample number for LAI result median for Lookup Table method as in Weiss et al 2004."
const LUT_NMEDIAN = 25 # number of points to sample (median) LAI from Weiss 2004
"""Default LAI starting point for `ellips_opt` method.
Better to use a simple method such as `zenith57` or `lang` for a realistic
LAI starting point for keyword argument `LAI_init`."""
const DEFAULT_LAI_INIT = 3.0
# min LAI value for warning and box constrained optimization. Was 0.2 for upward LAI, but reduced to 0.05 for crops.
const LAI_MIN = 0.05
const LAI_MAX = 10.0

abstract type InversionMethod end
struct EllipsOpt   <: InversionMethod end
struct Zenith57    <: InversionMethod end
struct Lang        <: InversionMethod end
struct Miller      <: InversionMethod end
struct MillerRings <: InversionMethod end
struct EllipsLUT   <: InversionMethod end

## Generic function
## ----------------

"Generic inversion function. Default method is EllipsOpt."
function inverse(polim::PolarImage, thresh = threshold(polim), im::InversionMethod=EllipsOpt(); kwargs...)
    inverse(polim, thresh, im; kwargs...)
end

# "Default inversion function for set of gapfractions using EllipsOpt."
# function inverse(θedges::AbstractArray, θmid::Vector{Float64}, K::Vector{Float64};
#                   kwargs...)
#     inverse(θedges, θmid, K, EllipsOpt(); kwargs...)
# end

"Calculate the default number of rings for integration of a Polar Image.
Because the logarithm of the gap fraction is calculated, there's a delicate trade off
between more rings for more algorithmic accuracy and more pixels per ring
(i.e. less rings) for more accurate log gap fraction per ring."
function Nrings_def(polim::PolarImage)::Int
    N =  polim.cl.fθρ(pi/2) / MILLER_GROUPS
    # hasslope(polim) && (N = max(NRINGS_MIN_SLOPE, N / INCIDENCE_GROUPS))
    return ceil(Int, N)
end
function maxviewangle(polim::PolarImage, θmax=θMAX)
    if hasslope(polim)
        α, ε = params(polim.slope)
        return min(θmax, pi/2 - α)
    else
        return θmax
    end
end

## Constant angle
## --------------

""" At constant zenith view angle of 1 rad the gapfraction is almost independent
of the leaf angle distribution."""
function inverse(polim::PolarImage, thresh, ::Zenith57; ringwidth = RING_WIDTH, args...)
    θedges, θmid, K = contactfreqs(polim, 1 - ringwidth, 1 + ringwidth, 1, thresh)
    G = 0.5
    K[1] / G
end

function inverse(θedges::AbstractArray, θmid::Vector{Float64}, K::Vector{Float64}, ::Zenith57)
    2 * K[1] # i.e. K[1] / 0.5
end

## Lang's method
## -------------

"Approximate the contact frequency with a first order regression."
function inverse(polim::PolarImage, thresh, ::Lang; θ1::Real = LANG_START, θ2::Real = LANG_END,
              Nrings = round(Int, Nrings_def(polim)*(θ2-θ1)/(pi/2)), kwargs...)
    checkθ1θ2(θ1, θ2)
    θedges, θmid, K = contactfreqs(polim, θ1, θ2, Nrings, thresh)
    inverse(θedges, θmid, K, Lang())
end

function inverse(θedges::AbstractArray, θmid::Vector{Float64}, K::Vector{Float64}, ::Lang)
    2 * sum(linreg(θmid, K))
end

## Miller's formula
## ----------------

"""
Miller's formula assumes a constant leaf angle. We divide the zenith range
in N rings each with an equal number of pixels.
"""
function inverse(polim::PolarImage, thresh, ::Miller; kwargs...)
    inverse(polim, thresh, MillerRings(); kwargs...)
end

"Create a N rings and integrate. This is the most robust way for Miller's method."
function inverse(polim::PolarImage, thresh, ::MillerRings;
                     N::Integer = Nrings_def(polim), θmax = θMAX, kwargs...)
    θedges, θmid, K = contactfreqs(polim, 0, θmax, N, thresh)
    inverse(θedges, θmid, K, MillerRings())
end

function inverse(θedges::AbstractArray, θmid::Vector{Float64}, K::Vector{Float64}, ::MillerRings)
    dθ = diff(θedges)
    2 * sum(K .* sin(θmid) .* dθ)
end

## Ellipsoidal
## -----------

"Ellipsoidal projection function, formula (A.6) Thimonier et al. 2010"
G_ellips(θᵥ, χ) = cos(θᵥ) * sqrt(χ^2 + tan(θᵥ)^2) / (χ + 1.702 * (χ + 1.12)^-0.708)

"function to convert ALIA to ellipsoidal parameter χ, formula (30) in Wang et al 2007"
ALIA_to_x(ALIA) = (ALIA / 9.65).^-0.6061 - 3

"The model that links ALIA and LAI with contact frequency K"
function model_ellips(θmid::Vector, params::Vector)
    alia, LAI = params
    K = [LAI * G_ellips(θᵥ, ALIA_to_x(alia)) for θᵥ in θmid]
end

"""
Assuming an ellipsoidal leaf angle distribution (with one parameter), we use the inverse model
to estimate the average leaf inclination angle (ALIA) and Leaf Area Index (LAI)
from the observed gap fraction per view zenith angle.

This method uses a curve fitting technique to find the optimal values for the leaf angle
distribution parameter ALIA and the LAI.
"""
function inverse(θedges::AbstractArray, θmid::Vector{Float64}, K::Vector{Float64}, ::EllipsOpt;
                 LAI_init::Float64 = DEFAULT_LAI_INIT)

    LAI_init == DEFAULT_LAI_INIT && @warn("Default value detected for `LAI_init`, it's better to use an estimate from `zenith57`.")
    if !(LAI_MIN < LAI_init < LAI_MAX)
        @warn("LAI starting point ($LAI_init) for optimization outside limits [$LAI_MIN, $LAI_MAX] and will be reset.")
        LAI_init = DEFAULT_LAI_INIT # max(min(LAI_init, 9), 0.2)
    end

    # Find an initial value for ALIA
    fitfunalia(alia) = sum((model_ellips(θmid, [alia, LAI_init]) .- K).^2)
    aliares = optimize(fitfunalia, 0.1, pi/2 - 0.1)
    ALIA_init = minimizer(aliares)

    # Optimize both ALIA and LAI at same time
    fitfun(x) = sum((K .- model_ellips(θmid, x)).^2)
    initial = [ALIA_init, LAI_init]
    fitdf = Optim.OnceDifferentiable(fitfun, initial; autodiff = :forward)

    try
        res = optimize(fitdf, initial, Optim.LBFGS())
        ALIA, LAI = minimizer(res)
        return LAI, ALIA
    catch y
        if isa(y, DomainError)
            # in case LsqFit.curve_fit does not converge and wanders out of the
            # parameter space, use box constrained optimization.
            lower = [0.05, LAI_MIN]
            upper = [pi/2 - 0.05, LAI_MAX]
            res = optimize(fitfun, lower, upper, initial, Optim.Fminbox(Optim.LBFGS()))
            ALIA, LAI = minimizer(res)
            return LAI, ALIA
        else
            throw(y)
        end
    end
    error("inverse EllipsOpt did not terminate normally")
end

function inverse(polim::PolarImage, thresh, ::EllipsOpt;
                    Nrings = Nrings_def(polim),
                    θmax = maxviewangle(polim),
                    Nτ = INCIDENCE_GROUPS)

    θedges, θmid, K = contactfreqs(polim, 0.0, θmax, Nrings, thresh;Nτ=Nτ)
    # Find inital value for LAI for optimization
    LAI_init = inverse(polim, thresh, Zenith57())
    inverse(θedges, θmid, K, EllipsOpt(); LAI_init=LAI_init)
end

## Ellipsoidal LUT
## ---------------

"Holds each element of the ellipsoidal Lookup Table method."
struct LUTel
    alia::Float64
    LAI::Float64
    modelled::Array{Float64}
end

"""Populates a Lookup Table (Weiss 2004) with random contact frequencies for a
range of ALIA and LAI values"""
function populateLUT(θmid::Vector{Float64}; Nlut = LUT_POINTS)
    LAI_max = LAI_MAX
    LUT = Array{LUTel}(Nlut)
    alia_max = pi/2 - .001 #against possible instability at π/2
    # As in paper [Weiss2004], populate LUT randomly.
    # TODO consider using Sobol pseudorandom numbers for consistency.
    for i = 1:Nlut
        alia = rand() * alia_max
        LAI = rand() * LAI_max
        modelled = model_ellips(θmid, [alia, LAI])
        LUT[i] = LUTel(alia, LAI, modelled)
    end
    return LUT
end

"""
Assuming an ellipsoidal leaf angle distribution (with one parameter), we use the
inverse model to estimate the average leaf inclination angle (ALIA) and Leaf Area Index (LAI)
from the observed gap fraction per view zenith angle.

This method uses a Lookup Table (LUT) for the leaf angle distribution parameter
ALIA and the LAI, following [Weiss2004](http://www.researchgate.net/profile/Inge_Jonckheere/publication/222931516_Review_of_methods_for_in_situ_leaf_area_index_(LAI)_determination_Part_II._Estimation_of_LAI_errors_and_sampling/links/09e4150cefe5a4fea5000000.pdf).
"""
function inverse(θedges::AbstractArray, θmid::Vector{Float64}, K::Vector{Float64}, ::EllipsLUT;
                    Nlut::Int = LUT_POINTS, Nmed::Integer=LUT_NMEDIAN, kwargs...)

    LUT = populateLUT(θmid; Nlut = Nlut)

    # create fitness function against observed contact frequencies
    LUTdiff = zeros(Float64, Nlut)
    for i in eachindex(LUT)
        # errfun as keyword argument does not accept FastAnonymous, so use SSE.
        #LUTdiff[i] = sum(map(errfun, LUT[i].modelled .- K))
        LUTdiff[i] = sum((LUT[i].modelled .- K).^2)
    end
    # get median of closest parameters
    LUTsortind = sortperm(LUTdiff)
    #LUT_alia = median([el.alia for el in LUT[LUTsortind[1:Nmed]]])
    LUT_LAI = median([el.LAI for el in LUT[LUTsortind[1:Nmed]]])
end

function inverse(polim::PolarImage, thresh, ::EllipsLUT;
                    Nrings = Nrings_def(polim), θmax = θMAX, Nlut::Int = LUT_POINTS,
                    Nmed::Integer = LUT_NMEDIAN, kwargs...)

    θedges, θmid, K = contactfreqs(polim, 0.0, θmax, Nrings, thresh)
    inverse(θedges, θmid, K, EllipsLUT(); Nlut = Nlut, Nmed = Nmed)
end


text:
#############################################################################
# multiply_divide.jl
# Handles scalar multiplication, matrix multiplication, and scalar division
# of variables, constants and expressions.
# All expressions and atoms are subtpyes of AbstractExpr.
# Please read expressions.jl first.
#############################################################################

import Base.Broadcast.broadcasted
export sign, monotonicity, curvature, evaluate, conic_form!

### Scalar and matrix multiplication

struct MultiplyAtom <: AbstractExpr
    head::Symbol
    id_hash::UInt64
    children::Tuple{AbstractExpr, AbstractExpr}
    size::Tuple{Int, Int}

    function MultiplyAtom(x::AbstractExpr, y::AbstractExpr)
        if x.size == (1, 1)
            sz = y.size
        elseif y.size == (1, 1)
            sz = x.size
        elseif x.size[2] ==  y.size[1]
            sz = (x.size[1], y.size[2])
        else
            error("Cannot multiply two expressions of sizes $(x.size) and $(y.size)")
        end
        children = (x, y)
        return new(:*, hash(children), children, sz)
    end
end

function sign(x::MultiplyAtom)
    return sign(x.children[1]) * sign(x.children[2])
end

function monotonicity(x::MultiplyAtom)
    return (sign(x.children[2]) * Nondecreasing(), sign(x.children[1]) * Nondecreasing())
end

# Multiplication has an indefinite hessian, so if neither children are constants,
# the curvature of the atom will violate DCP.
function curvature(x::MultiplyAtom)
    if vexity(x.children[1]) != ConstVexity() && vexity(x.children[2]) != ConstVexity()
        return NotDcp()
    else
        return ConstVexity()
    end
end

function evaluate(x::MultiplyAtom)
    return evaluate(x.children[1]) * evaluate(x.children[2])
end

function conic_form!(x::MultiplyAtom, unique_conic_forms::UniqueConicForms)
    if !has_conic_form(unique_conic_forms, x)
        # scalar multiplication
        if x.children[1].size == (1, 1) || x.children[2].size == (1, 1)
            if vexity(x.children[1]) == ConstVexity()
                const_child = x.children[1]
                expr_child = x.children[2]
            elseif vexity(x.children[2]) == ConstVexity()
                const_child = x.children[2]
                expr_child = x.children[1]
            else
                error("multiplication of two non-constant expressions is not DCP compliant")
            end
            objective = conic_form!(expr_child, unique_conic_forms)

            # make sure all 1x1 sized objects are interpreted as scalars, since
            # [1] * [1, 2, 3] is illegal in julia, but 1 * [1, 2, 3] is ok
            if const_child.size == (1, 1)
                const_multiplier = evaluate(const_child)[1]
            else
                const_multiplier = reshape(evaluate(const_child), length(const_child), 1)
            end

            objective = const_multiplier * objective

        # left matrix multiplication
        elseif vexity(x.children[1]) == ConstVexity()
            objective = conic_form!(x.children[2], unique_conic_forms)
            objective = kron(sparse(1.0I, x.size[2], x.size[2]), evaluate(x.children[1])) * objective
        # right matrix multiplication
        else
            objective = conic_form!(x.children[1], unique_conic_forms)
            objective = kron(transpose(evaluate(x.children[2])), sparse(1.0I, x.size[1], x.size[1])) * objective
        end
        cache_conic_form!(unique_conic_forms, x, objective)
    end
    return get_conic_form(unique_conic_forms, x)
end

function *(x::AbstractExpr, y::AbstractExpr)
    if hash(x) == hash(y) && x.size == (1, 1)
        return square(x)
    end
    return MultiplyAtom(x, y)
end

*(x::Value, y::AbstractExpr) = MultiplyAtom(Constant(x), y)
*(x::AbstractExpr, y::Value) = MultiplyAtom(x, Constant(y))
/(x::AbstractExpr, y::Value) = MultiplyAtom(x, Constant(1 ./ y))

### .*
# All constructors of this check (and so this function requires)
# that the first child be constant to have the expression be DCP
struct DotMultiplyAtom <: AbstractExpr
    head::Symbol
    id_hash::UInt64
    children::Tuple{AbstractExpr, AbstractExpr}
    size::Tuple{Int, Int}

    function DotMultiplyAtom(x::AbstractExpr, y::AbstractExpr)
        # check that the sizes of x and y are compatible
        try
            ones(size(x)) .* ones(size(y))
        catch
            error("cannot compute $x .* $y: sizes are not compatible")
        end
        children = (x, y)
        return new(:.*, hash(children), children, y.size)
    end
end

function sign(x::DotMultiplyAtom)
    return sign(x.children[1]) * sign(x.children[2])
end

function monotonicity(x::DotMultiplyAtom)
    return (sign(x.children[2]) * Nondecreasing(), sign(x.children[1]) * Nondecreasing())
end

function curvature(x::DotMultiplyAtom)
    if vexity(x.children[1]) == ConstVexity()
        return ConstVexity()
    else
        return NotDcp()
    end
end

function evaluate(x::DotMultiplyAtom)
    return evaluate(x.children[1]) .* evaluate(x.children[2])
end

function conic_form!(x::DotMultiplyAtom, unique_conic_forms::UniqueConicForms)
    if !has_conic_form(unique_conic_forms, x)
        if vexity(x.children[1]) != ConstVexity()
            if vexity(x.children[2]) != ConstVexity()
                error("multiplication of two non-constant expressions is not DCP compliant")
            else
                # make sure first child is the one that's constant
                x.children[1], x.children[2] = x.children[2], x.children[1]
            end
        end
        # promote the size of the coefficient matrix, so eg
        # 3 .* x
        # works regardless of the size of x
        coeff = evaluate(x.children[1]) .* ones(size(x.children[2]))
        # promote the size of the variable
        # we've previously ensured neither x nor y is 1x1
        # and that the sizes are compatible,
        # so if the sizes aren't equal the smaller one is size 1
        var = x.children[2]
        if size(var, 1) < size(coeff, 1)
            var = ones(size(coeff, 1)) * var
        elseif size(var, 2) < size(coeff, 2)
            var = var * ones(1, size(coeff, 1))
        end

        const_multiplier = spdiagm(0 => vec(coeff))
        objective = const_multiplier * conic_form!(var, unique_conic_forms)
        cache_conic_form!(unique_conic_forms, x, objective)
    end
    return get_conic_form(unique_conic_forms, x)
end

function broadcasted(::typeof(*), x::Constant, y::AbstractExpr)
    if x.size == (1, 1) || y.size == (1, 1)
        return x * y
    elseif size(y, 1) < size(x, 1) && size(y, 1) == 1
        return DotMultiplyAtom(x, ones(size(x, 1)) * y)
    elseif size(y, 2) < size(x, 2) && size(y, 2) == 1
        return DotMultiplyAtom(x, y * ones(1, size(x, 1)))
    else
        return DotMultiplyAtom(x, y)
    end
end
broadcasted(::typeof(*), y::AbstractExpr, x::Constant) = DotMultiplyAtom(x, y)

# if neither is a constant it's not DCP, but might be nice to support anyway for eg MultiConvex
function broadcasted(::typeof(*), x::AbstractExpr, y::AbstractExpr)
    if x.size == (1, 1) || y.size == (1, 1)
        return x * y
    elseif vexity(x) == ConstVexity()
        return DotMultiplyAtom(x, y)
    elseif hash(x) == hash(y)
        return square(x)
    else
        return DotMultiplyAtom(y, x)
    end
end
broadcasted(::typeof(*), x::Value, y::AbstractExpr) = DotMultiplyAtom(Constant(x), y)
broadcasted(::typeof(*), x::AbstractExpr, y::Value) = DotMultiplyAtom(Constant(y), x)
broadcasted(::typeof(/), x::AbstractExpr, y::Value) = DotMultiplyAtom(Constant(1 ./ y), x)
# x ./ y and x / y for x constant, y variable is defined in second_order_cone.qol_elemwise.jl


text:
module  WindID
using Statistics
using Printf
#External Alg Libray :
# AlfFileFolder="../AlarmAnalyssys/RAW3Alg-RTWind20190417/" #
# Base.MainInclude.include( "$(AlfFileFolder)cumSumAlg.jl")
# Base.MainInclude.include( "$(AlfFileFolder)AlarmGen.jl")
# Base.MainInclude.include( "$(AlfFileFolder)C01-RAW3AlgEntryRT.jl" )

include("../../RawFileofflineProc/BaseLib.jl")
##internal Use Libaray
# Base.include(Main, "ReadFiberData.jl")
include("../../RawFileofflineProc/dataStruct.jl")
# Base.MainInclude.
include( "../../RawFileofflineProc/ReadFiberData.jl")
# Base.MainInclude.
include("../../RawFileofflineProc/ConfigFileLib.jl")
# Base.MainInclude.
include("cumSumAlg.jl")
# Base.MainInclude.
include("AlarmGen.jl")

# @warn "OK？？？？？"
lastWriteDT= UInt(0)
# cfgObj = cfgStruct()
CalChn1Obj = spaceCumSumData( )
CalChn2Obj = spaceCumSumData( )
function Init( cfgObj )
    curpath = pwd()
    jlPath =joinpath(curpath ,"../Jl-Aux/glbCfg/WindIDCfg.jl")
    getCfg( Base.@__MODULE__,jlPath)
    # println("\n",windAlarmSwitchOFFMinTIME,"\n")

    InitObj(CalChn1Obj, Int(1),Int(cfgObj.chn2SegBegIdx-1), Int(1),Int(cfgObj.chn2SegBegIdx-1),-1 );
    InitObj(CalChn2Obj, Int(cfgObj.chn2SegBegIdx),Int( 8192) ,Int(cfgObj.chn2SegBegIdx), Int(8192),-1 );
    CalChn1Obj.CfgInfoObj = cfgObj
    CalChn2Obj.CfgInfoObj = cfgObj
end
function WindProc(row,startDTStr,machineName)
    global lastWriteDT

    # a = String(startDTStr)
    # @warn  startDTStr,String(machineName) ,"-------->Wind------"
    CalChn1Obj.dt = startDTStr
    CalChn2Obj.dt = startDTStr
    # Name = getStringFromArray( machineName )
    # evs ="machineName =\"$Name\""
    # reps = Base.Meta.parse( evs )
    # Core.eval(Main,reps )
    postProcFeatureData(CalChn2Obj ,row ) #include("AlarmGen.jl")
    postProcFeatureData(CalChn1Obj ,row ) #include("AlarmGen.jl")
    curDt = DT2UInt( startDTStr )
    a1 = showCurDT2File(lastWriteDT, curDt)
    a1 != 0 ?  lastWriteDT = UInt( a1 ) : 0 ;
end
end

max_stars_repo_name: wangzhongtian/fiberAlarmAnalyseTool

max_issues_repo_name: wangzhongtian/fiberAlarmAnalyseTool

max_forks_repo_name: wangzhongtian/fiberAlarmAnalyseTool


text:
function value_estimate(π::DiscreteNetwork, px, s, 𝒫)
        pdfs = logits(px, s)
        sum(value(π, s) .* pdfs, dims=1)
end

# function value_estimate(π::ActorCritic, px, s, 𝒫)
#         mean([value(π, s, action(px, s)) for _ in 1:𝒫[:N_samples]])
# end

#NOTE: Currently all of these assume that we get a reward (cost) ONCE at the end of an episode

# function E_target(π, 𝒫, 𝒟, γ::Float32; kwargs...)
#         px = 𝒫[:px]
#         if 𝒫[:use_likelihood_weights]
#                 return 𝒟[:likelihoodweight] .* (𝒟[:done] .* 𝒟[:r] .+ (1.f0 .- 𝒟[:done]) .* value_estimate(π, px, 𝒟[:sp], 𝒫))
#         else
#                 return 𝒟[:done] .* 𝒟[:r] .+ (1.f0 .- 𝒟[:done]) .* value_estimate(π, px, 𝒟[:sp], 𝒫)
#         end
# end

function CDF_target(π, 𝒫, 𝒟, γ::Float32; kwargs...)
        rα = 𝒫[:rα][1]
        # rs = 𝒫[:rs]
        # stdrα = 𝒫[:std_rα][1]
        # vard = Normal(rα, stdrα)
        px = 𝒫[:px]
        # B = length(𝒟[:r])

        # y = hcat([𝒟[:done] .* (𝒟[:r] .> rα) .+ (1.f0 .- 𝒟[:done]) .* value_estimate(π.policy, px, vcat(repeat([rα], 1, B), 𝒟[:sp]), 𝒫) for rα in rs]...)
        # return y
        # return (𝒟[:var_prob] .+ y) ./ 2f0
        # return 𝒟[:var_prob]

        # if 𝒫[:use_likelihood_weights]
        #         return 𝒟[:likelihoodweight] .* (𝒟[:done] .* (𝒟[:r] .> rα) .+ (1.f0 .- 𝒟[:done]) .* value_estimate(π, px, 𝒟[:sp], 𝒫))
        # else
         return 𝒟[:done] .* (𝒟[:r] .> rα) .+ (1.f0 .- 𝒟[:done]) .* value_estimate(π, px, 𝒟[:sp], 𝒫)
        #         # return 𝒟[:done] .* cdf.(vard, 𝒟[:r]) .+ (1.f0 .- 𝒟[:done]) .* value_estimate(π, px, 𝒟[:sp], 𝒫)
        # end
end

function CVaR_target(π, 𝒫, 𝒟, γ::Float32; kwargs...)
        rα = 𝒫[:rα][1]
        # stdrα = 𝒫[:std_rα][1]
        # vard = Normal(rα, stdrα)
        px = 𝒫[:px]

        # y = 𝒟[:done] .* 𝒟[:r] .* (𝒟[:r] .> rα) .+ (1.f0 .- 𝒟[:done]) .* value_estimate(π, px, 𝒟[:sp], 𝒫)
        # return y
        # return (𝒟[:cvar_prob] .+ y) ./ 2f0
        # return 𝒟[:cvar_prob]

        # if 𝒫[:use_likelihood_weights]
        #         return 𝒟[:likelihoodweight] .* (𝒟[:done] .* 𝒟[:r] .* (𝒟[:r] .> rα) .+ (1.f0 .- 𝒟[:done]) .* value_estimate(π, px, 𝒟[:sp], 𝒫))
        # else
        return 𝒟[:done] .* 𝒟[:r] .* (𝒟[:r] .> rα) .+ (1.f0 .- 𝒟[:done]) .* value_estimate(π, px, 𝒟[:sp], 𝒫)
        #         # return 𝒟[:done] .* 𝒟[:r] .* cdf.(vard, 𝒟[:r]) .+ (1.f0 .- 𝒟[:done]) .* value_estimate(π, px, 𝒟[:sp], 𝒫)
        # end
end

# function E_VaR_CVaR_target(π, 𝒫, 𝒟, γ::Float32; kwargs...)
        # [CDF_target(π.networks[1], 𝒫, 𝒟, γ; kwargs...), CVaR_target(π.networks[2], 𝒫, 𝒟, γ; kwargs...), E_target(π.networks[3], 𝒫, 𝒟, γ; kwargs...)]
# end

function VaR_CVaR_target(π, 𝒫, 𝒟, γ::Float32; kwargs...)
        [CDF_target(π.networks[1], 𝒫, 𝒟, γ; kwargs...), CVaR_target(π.networks[2], 𝒫, 𝒟, γ; kwargs...)]
end

# function IS_L_KL_log(π, 𝒫, 𝒟; info=Dict(), kwargs...)
#         x, logqx = exploration(π, 𝒟[:s])
#         Q_s = Zygote.ignore() do
#                 value_estimate(π, 𝒫[:px], 𝒟[:s], 𝒫)
#         end
#         px = exp.(logpdf(𝒫[:px], 𝒟[:s], x))
#         Q_sx = value(π, 𝒟[:s], x)
#
#         qstar = px .* Q_sx #./ Q_s
#
#
#         -mean(qstar .* logqx ./ exp.(logqx))
# end

# function fill_probs(D; 𝒮, info=Dict())
#         rα = 𝒮.𝒫[:rα][1]
#         epies = !(D isa ExperienceBuffer) ? episodes(ExperienceBuffer(D)) : episodes(D)
#         for ep in epies
#                 episode_range = ep[1]:ep[2]
#                 r = sum(D[:r][1, episode_range]) # total return
#                 varprob = r >= rα
#                 cvarprob = r * varprob
#
#                 for i in reverse(episode_range)
#                         D[:var_prob][:, i] .= varprob
#                         D[:cvar_prob][:, i] .= cvarprob
#
#                         varprob = D[:likelihoodweight][1,i] * varprob
#                         cvarprob = D[:likelihoodweight][1,i] * cvarprob
#                 end
#         end
# end

function ISDRL_Discrete(;π,
                        S,
                        N,
                        px,
                        prioritized=true,
                        use_likelihood_weights=true,
                        α,
                        target_fn=VaR_CVaR_target,
                        𝒫=(;),
                        buffer_size=N,
                        ΔN=4,
                        pre_train_callback,
                        π_explore=π,
                        c_opt::NamedTuple=(;),
                        log::NamedTuple=(;),
                        c_loss,
                        kwargs...)

                    𝒫 = (;px, rα=Float32[NaN], α, use_likelihood_weights, 𝒫...)
                    required_columns=[:logprob, :likelihoodweight, :var_prob, :cvar_prob]
                    agent = PolicyParams(π=π, π_explore=π_explore, π⁻=deepcopy(π), pa=px)
                    OffPolicySolver(;agent=agent,
                                     S=S,
                                     log=LoggerParams(;dir="log/isdrl_dqn", period=100, fns=[log_episode_averages([:r], 100)], log...),
                                     𝒫=𝒫,
                                     N=N,
                                     ΔN=ΔN,
                                     # post_sample_callback=fill_probs,
                                     pre_train_callback=pre_train_callback,
                                     buffer=ExperienceBuffer(S, agent.space, buffer_size, required_columns, prioritized=prioritized),
                                     c_opt = TrainingParams(;loss=c_loss, name="critic_", epochs=ΔN, c_opt...),
                                     target_fn=target_fn,
                                     kwargs...)
end

# function ISDRL_Continuous(;π::MixtureNetwork,
#                         S,
#                         N,
#                         px,
#                         N_samples=10,
#                         prioritized=true,
#                         use_likelihood_weights=true,
#                         α,
#                         𝒫=(;),
#                         buffer_size=N,
#                         ΔN=4,
#                         pre_train_callback,
#                         π_explore=π,
#                         c_opt::NamedTuple=(;),
#                         a_opt::NamedTuple=(;),
#                         log::NamedTuple=(;),
#                         c_loss,
#                         a_loss = IS_L_KL_log,
#                         kwargs...)
#
#                     𝒫 = (;px, rα=[NaN], α, use_likelihood_weights, N_samples, 𝒫...)
#                     required_columns=[:logprob, :likelihoodweight]
#                     agent = PolicyParams(π=π, π_explore=π_explore, π⁻=deepcopy(π), pa=px)
#                     OffPolicySolver(;agent=agent,
#                                      S=S,
#                                      log=LoggerParams(;dir="log/cerl_dqn", period=100, fns=[log_episode_averages([:r], 100)], log...),
#                                      𝒫=𝒫,
#                                      N=N,
#                                      ΔN=ΔN,
#                                      pre_train_callback=pre_train_callback,
#                                      buffer=ExperienceBuffer(S, agent.space, buffer_size, required_columns, prioritized=prioritized),
#                                      c_opt = TrainingParams(;loss=c_loss, name="critic_", epochs=ΔN, c_opt...),
#                                      a_opt = TrainingParams(;loss=a_loss, name="actor_", a_opt...),
#                                      target_fn=VaR_CVaR_target,
#                                      kwargs...)
# end
#


text:
"""
# Decomposition of graphs

- Why graphs?
    - How is a graph represented?
- Depth-first search in undirected graphs
    - Exploring mazes
    - Depth-first search
    - Connectivity in undirected graphs
    - Previsit and postvisit orderings
- Depth-first search in directed graphs
    - Types of edges
    - Directed acyclic graphs
- Strongly connected components
    - Defining connectivity for directed graphs
    - An efficient algorithm
"""
module Chapter3

export
    explore,
    dfs

import
    ..Graph

mutable struct Clock
    value::Int
end

function previsit!(v, pre::Dict, clock::Clock)
    pre[v] = clock.value
    clock.value += 1
end

function postvisit!(v, post::Dict, clock::Clock)
    post[v] = clock.value
    clock.value += 1
end

"""
    explore(G, v)

Find all nodes reachable from a particular node.
"""
function explore(G::Graph{T}, v::T,
                 visited::Dict=Dict(collect(keys(G)) .=> false),
                 pre::Dict=Dict(collect(keys(G)) .=> 0),
                 post::Dict=Dict(collect(keys(G)) .=> 0),
                 clock::Clock=Clock(0)) where {T}
    visited[v] = true
    previsit!(v, pre, clock)
    for edge ∈ G[v]
        (u, w) = edge
        if !visited[u]
            explore(G, u, visited, pre, post, clock)
        end
    end
    postvisit!(v, post, clock)
    return (pre=pre, post=post)
end

"""
    dfs(G, v)

Depth-first search algorithm.
"""
function dfs(G::Graph{T}, v::T, visited::Dict=Dict(collect(keys(G)) .=> false)) where {T}
    explored = []
    for (v, e) ∈ G
        if !visited[v]
            push!(explored, (v, explore(G, v, visited)))
        end
    end
    return explored
end

end # module


text:
include("default_params.jl")
include("model_def.jl")

ks_errors = DataFrame(DG=[],DG_limit=[],Order_1=[],QBDRAP=[])
l1_errors = DataFrame(DG=[],DG_limit=[],Order_1=[],QBDRAP=[])
l1_cell_probs_errors = DataFrame(DG=[],DG_limit=[],Order_1=[],QBDRAP=[])
ks_error_row = zeros(4)
l1_error_row = zeros(4)
l1_cell_probs_row = zeros(4)

sim_cdf = CSV.read(
    (@__DIR__)*"/data/sims_cdf_evaluated.csv",
    DataFrame,
)

for o in orders
    cols = [1;3;4]
    for (c,m) in enumerate(string.(approx_types[[1;3;4]]).*string.([1;3;4]))
        approx_data = CSV.read(
            (@__DIR__)*"/data/first_return_cdf_approximations/order_"*string(o)*"_model_"*m*".csv",
            DataFrame,
        )
        ks_error = max(
            maximum(abs.(approx_data.phase_2-sim_cdf.phase_2)),
            maximum(abs.(approx_data.phase_4-sim_cdf.phase_4)),
        )
        ks_error_row[cols[c]] = ks_error
        l1_cdf = sum(abs.(approx_data.phase_2[1:end-1]+approx_data.phase_2[2:end]-
                sim_cdf.phase_2[1:end-1]-sim_cdf.phase_2[2:end]).*diff(approx_data.x))/2 +
            sum(abs.(approx_data.phase_4[1:end-1]+approx_data.phase_4[2:end]-
                sim_cdf.phase_4[1:end-1]-sim_cdf.phase_4[2:end]).*diff(approx_data.x))/2
        l1_error_row[cols[c]] = l1_cdf
        cell_edges_idx = findall(mod.(approx_data.x,0.4).==0.0)
        sim_cell_probs = DataFrame(
            x=[approx_data.x[1]; approx_data.x[cell_edges_idx[2:end]] - approx_data.x[cell_edges_idx[1:end-1]]],
            phase_2=[sim_cdf.phase_2[1]; sim_cdf.phase_2[cell_edges_idx[2:end]] - sim_cdf.phase_2[cell_edges_idx[1:end-1]]],
            phase_4=[sim_cdf.phase_4[1]; sim_cdf.phase_4[cell_edges_idx[2:end]] - sim_cdf.phase_4[cell_edges_idx[1:end-1]]],
        )
        approx_cell_probs = DataFrame(
            x=[approx_data.x[1]; approx_data.x[cell_edges_idx[2:end]] - approx_data.x[cell_edges_idx[1:end-1]]],
            phase_2=[approx_data.phase_2[1]; approx_data.phase_2[cell_edges_idx[2:end]] - approx_data.phase_2[cell_edges_idx[1:end-1]]],
            phase_4=[approx_data.phase_4[1]; approx_data.phase_4[cell_edges_idx[2:end]] - approx_data.phase_4[cell_edges_idx[1:end-1]]],
        )
        l1_cell_probs_row[cols[c]] = sum(abs.(approx_cell_probs.phase_2-sim_cell_probs.phase_2)) +
            sum(abs.(approx_cell_probs.phase_4-sim_cell_probs.phase_4))
    end
    push!(ks_errors,ks_error_row)
    push!(l1_errors,l1_error_row)
    push!(l1_cell_probs_errors,l1_cell_probs_row)
end

CSV.write((@__DIR__)*"/data/errors/ks_error.csv",ks_errors)
CSV.write((@__DIR__)*"/data/errors/l1_error.csv",l1_errors)
CSV.write((@__DIR__)*"/data/errors/l1_cell_probs_error.csv",l1_cell_probs_errors)

ks_errors_lwr = DataFrame(DG=[],DG_limit=[],Order_1=[],QBDRAP=[])
l1_errors_lwr = DataFrame(DG=[],DG_limit=[],Order_1=[],QBDRAP=[])
l1_cell_probs_errors_lwr = DataFrame(DG=[],DG_limit=[],Order_1=[],QBDRAP=[])
ks_error_row_lwr = zeros(4)
l1_error_row_lwr = zeros(4)
l1_cell_probs_error_row_lwr = zeros(4)

ks_errors_upr = DataFrame(DG=[],DG_limit=[],Order_1=[],QBDRAP=[])
l1_errors_upr = DataFrame(DG=[],DG_limit=[],Order_1=[],QBDRAP=[])
l1_cell_probs_errors_upr = DataFrame(DG=[],DG_limit=[],Order_1=[],QBDRAP=[])
ks_error_row_upr = zeros(4)
l1_error_row_upr = zeros(4)
l1_cell_probs_error_row_upr = zeros(4)

ks_error_boot = zeros(n_boot)
l1_error_boot = zeros(n_boot)
l1_cell_probs_error_boot = zeros(n_boot)
for o in orders
    cols = [1;3;4]
    for (c,m) in enumerate(string.(approx_types[[1;3;4]]).*string.([1;3;4]))
        approx_data = CSV.read(
            (@__DIR__)*"/data/first_return_cdf_approximations/order_"*string(o)*"_model_"*m*".csv",
            DataFrame,
        )
        cell_edges_idx = findall(mod.(approx_data.x,0.4).==0.0)
        approx_cell_probs = DataFrame(
                x=[approx_data.x[1]; approx_data.x[cell_edges_idx[2:end]] - approx_data.x[cell_edges_idx[1:end-1]]],
                phase_2=[approx_data.phase_2[1]; approx_data.phase_2[cell_edges_idx[2:end]] - approx_data.phase_2[cell_edges_idx[1:end-1]]],
                phase_4=[approx_data.phase_4[1]; approx_data.phase_4[cell_edges_idx[2:end]] - approx_data.phase_4[cell_edges_idx[1:end-1]]],
            )
        for n in 1:n_boot
            sim_cdf = CSV.read(
                (@__DIR__)*"/data/bootstrap_first_return/sample_"*string(n)*".csv",
                DataFrame,
            )
            ks_error = max(
                maximum(abs.(approx_data.phase_2-sim_cdf.phase_2)),
                maximum(abs.(approx_data.phase_4-sim_cdf.phase_4)),
            )
            ks_error_boot[n] = ks_error
            l1_cdf = sum(abs.(approx_data.phase_2[1:end-1]+approx_data.phase_2[2:end]-
                    sim_cdf.phase_2[1:end-1]-sim_cdf.phase_2[2:end]).*diff(approx_data.x))/2 +
                sum(abs.(approx_data.phase_4[1:end-1]+approx_data.phase_4[2:end]-
                    sim_cdf.phase_4[1:end-1]-sim_cdf.phase_4[2:end]).*diff(approx_data.x))/2
            l1_error_boot[n] = l1_cdf
            sim_cell_probs = DataFrame(
                x=[approx_data.x[1]; approx_data.x[cell_edges_idx[2:end]] - approx_data.x[cell_edges_idx[1:end-1]]],
                phase_2=[sim_cdf.phase_2[1]; sim_cdf.phase_2[cell_edges_idx[2:end]] - sim_cdf.phase_2[cell_edges_idx[1:end-1]]],
                phase_4=[sim_cdf.phase_4[1]; sim_cdf.phase_4[cell_edges_idx[2:end]] - sim_cdf.phase_4[cell_edges_idx[1:end-1]]],
            )
            l1_cell_probs_error_boot[n] = sum(abs.(approx_cell_probs.phase_2-sim_cell_probs.phase_2)) +
                sum(abs.(approx_cell_probs.phase_4-sim_cell_probs.phase_4))
        end
        ks_error_row_lwr[cols[c]] = quantile(ks_error_boot,qtiles[1])
        l1_error_row_lwr[cols[c]] = quantile(l1_error_boot,qtiles[1])
        l1_cell_probs_error_row_lwr[cols[c]] = quantile(l1_cell_probs_error_boot,qtiles[1])
        ks_error_row_upr[cols[c]] = quantile(ks_error_boot,qtiles[2])
        l1_error_row_upr[cols[c]] = quantile(l1_error_boot,qtiles[2])
        l1_cell_probs_error_row_upr[cols[c]] = quantile(l1_cell_probs_error_boot,qtiles[2])
    end
    push!(ks_errors_lwr,ks_error_row_lwr)
    push!(l1_errors_lwr,l1_error_row_lwr)
    push!(l1_cell_probs_errors_lwr,l1_cell_probs_error_row_lwr)
    push!(ks_errors_upr,ks_error_row_upr)
    push!(l1_errors_upr,l1_error_row_upr)
    push!(l1_cell_probs_errors_upr,l1_cell_probs_error_row_upr)
end

CSV.write((@__DIR__)*"/data/errors/ks_error_lwr.csv",ks_errors_lwr)
CSV.write((@__DIR__)*"/data/errors/l1_error_lwr.csv",l1_errors_lwr)
CSV.write((@__DIR__)*"/data/errors/l1_cell_probs_error_lwr.csv",l1_cell_probs_errors_lwr)
CSV.write((@__DIR__)*"/data/errors/ks_error_upr.csv",ks_errors_upr)
CSV.write((@__DIR__)*"/data/errors/l1_error_upr.csv",l1_errors_upr)
CSV.write((@__DIR__)*"/data/errors/l1_cell_probs_error_upr.csv",l1_cell_probs_errors_upr)


text:
using FeynmanKleinert
using Base.Test

# write your own tests here
include("Feynman-Kleinert-1986-PRA-Table2.jl")


text:
module test_26_TooManyLegends

using  ModiaPlot
include("test_base/test_26_TooManyLegends.jl")

end


text:
# Written in Julia by Andreas Noack Jensen
# January 2013
#
# Translation of Fortran file tvpack.f authored by
#
# Alan Genz
# Department of Mathematics
# Washington State University
# Pullman, WA 99164-3113
# Email : alangenz@wsu.edu
#
# Original source available from
# http://www.math.wsu.edu/faculty/genz/software/fort77/tvpack.f

function tvtcdf(nu::Int, h::Vector{Float64}, r::Vector{Float64})
	pt = 0.5 * pi
	h1 = h[1]
	h2 = h[2]
	h3 = h[3]
	r12 = r[1]
	r13 = r[2]
	r23 = r[3]

	# Sort R's and check for special cases

	if abs(r12) > abs(r13)
	   h2 = h3
	   h3 = h[2]
	   r12 = r13
	   r13 = r[1]
	end
	if abs(r13) > abs(r23)
	   h1 = h2
	   h2 = h[1]
	   r23 = r13
	   r13 = r[3]
	end
	tvt = 0.0
	if abs(h1) + abs(h2) + abs(h3) < eps()
	   	tvt = (1.0 + (asin(r12) + asin(r13) + asin(r23)) / pt) * 0.125
	elseif nu < 1 && abs(r12) + abs(r13) < eps()
	   	tvt = cdf(Normal(), h1) * bvtcdf(nu, h2, h3, r23)
	elseif nu < 1 && abs(r13) + abs(r23) < eps()
	   	tvt = cdf(Normal(), h3) * bvtcdf(nu, h1, h2, r12)
	elseif nu < 1 && abs(r12) + abs(r23) < eps()
	   	tvt = cdf(Normal(), h2) * bvtcdf(nu, h1, h3, r13)
	elseif 1.0 - r23 < eps()
	   	tvt = bvtnorm(nu, h1, min(h2, h3), r12)
	elseif r23 + 1.0 < eps()
	   	if h2 > -h3 tvt = bvtcdf(nu, h1, h2, r12) - bvtcdf(nu, h1, -h3, r12) end
	else

	   	# Compute singular TVT value

	   	if nu < 1
	      	tvt = bvtcdf(nu, h2, h3, r23) * cdf(Normal(), h1)
	   	elseif r23 > 0
	      	tvt = bvtcdf( nu, h1, min( h2, h3 ), 0.0)
	   	elseif h2 > -h3
	      	tvt = bvtcdf(nu, h1, h2, 0.0) - bvtcdf(nu, h1, -h3, 0.0)
	   	end
	   	# Use numerical integration to compute probability

	   	rua = asin(r12)
	   	rub = asin(r13)
	   	ar = asin(r23)
	   	ruc = pt * sign(ar) - ar
	   	tvt += adonet(x->tvtmfn(x, rua, rub, nu, h1, h2, h3, r23, ar, ruc), 0.0, 1.0, sqrt(eps())) / (4.0*pt)
	end
	return max(0.0, min(tvt, 1.0))
end

#Computes Plackett formula integrands
function tvtmfn(x::Float64, rua::Float64, rub::Float64, nu::Int, h1::Float64, h2::Float64, h3::Float64, r23::Float64, ar::Float64, ruc::Float64)
	retval = 0.0
	r12, rr2 = sincs(rua * x)
	r13, rr3 = sincs(rub * x)
	if abs(rua) > 0
     	retval += rua * pntgnd(nu, h1, h2, h3, r13, r23, r12, rr2)
	end
	if abs(rub) > 0
     	retval += rub * pntgnd(nu, h1, h3, h2, r12, r23, r13, rr3)
	end
	if nu > 0
   		r, rr = sincs(ar + ruc * x)
   		retval -= ruc * pntgnd(nu, h2, h3, h1, 0.0, 0.0, r, rr)
	end
	return retval
end

# Computes SIN(X), COS(X)^2, with series approx. for |X| near PI/2
function sincs(x::Float64)
	ee = (0.5 * pi - abs(x))^2
	if ee < 5e-5
   		sx = abs(1.0 - ee * (1.0 - ee / 12.0) * 0.5) * sign(x)
   		cs = ee * (1.0 - ee * (1.0 - 2.0 * ee / 15.0) / 3.0)
	else
   		sx = sin(x)
   		cs = 1.0 - sx*sx
	end
	return sx, cs
end

# Computes Plackett formula integrand
function pntgnd(nu::Int, ba::Float64, bb::Float64, bc::Float64, ra::Float64, rb::Float64, r::Float64, rr::Float64)
	retval = 0.0
	dt = rr * (rr - (ra - rb)^2 - 2.0 * ra * rb * (1.0 - r))
	if dt > 0
   		bt = (bc * rr + ba * (r * rb - ra) + bb * (r * ra - rb)) / sqrt(dt)
   		ft = (ba - r * bb)^2 / rr + bb * bb
   		if nu < 1
      		if bt > -10 && ft < 100
         		retval = exp(-ft * 0.5)
         		if bt < 10; retval *= cdf(Normal(), bt) end
      		end
   		else
      		ft = sqrt(1.0 + ft / nu)
      		retval = tcdf(nu, bt / ft) / ft^nu
   		end
	end
	return retval
end

# # One Dimensional Globally Adaptive Integration Function

function adonet(f::Function, a::Float64, b::Float64, tol::Float64)
	ai = Array(Float64, 100)
	bi = Array(Float64, 100)
	ei = Array(Float64, 100)
	fi = Array(Float64, 100)
	ai[1] = a
	bi[1] = b
	err = 1.0
	ip = 1
	im = 1
	fin = 0.0
	while 4.0err > tol && im < 100
	   	im += 1
	   	bi[im] = bi[ip]
	   	ai[im] = (ai[ip] + bi[ip]) * 0.5
	   	bi[ip] = ai[im]
	   	fi[ip], ei[ip] = krnrdt(ai[ip], bi[ip], f)
		fi[im], ei[im] = krnrdt(ai[im], bi[im], f)
	   	err = 0.0
	   	fin = 0.0
	   	for i = 1:im
	      	if ei[i] > ei[ip]; ip = i; end
	      	fin += fi[i]
	      	err += ei[i]^2
	   	end
	   	err = sqrt(err)
	end
	return fin
end

# # Kronrod Rule

const kronrod_wg = [0.2729250867779007e+00,
					0.5566856711617449e-01,
					0.1255803694649048e+00,
					0.1862902109277352e+00,
					0.2331937645919914e+00,
					0.2628045445102478e+00]

const kronrod_xgk = [0.0000000000000000e+00,
					 0.9963696138895427e+00,
					 0.9782286581460570e+00,
					 0.9416771085780681e+00,
					 0.8870625997680953e+00,
					 0.8160574566562211e+00,
					 0.7301520055740492e+00,
					 0.6305995201619651e+00,
					 0.5190961292068118e+00,
					 0.3979441409523776e+00,
					 0.2695431559523450e+00,
					 0.1361130007993617e+00]

const kronrod_wgk = [0.1365777947111183e+00,
					 0.9765441045961290e-02,
					 0.2715655468210443e-01,
					 0.4582937856442671e-01,
					 0.6309742475037484e-01,
					 0.7866457193222764e-01,
					 0.9295309859690074e-01,
					 0.1058720744813894e+00,
					 0.1167395024610472e+00,
					 0.1251587991003195e+00,
					 0.1312806842298057e+00,
					 0.1351935727998845e+00]

function krnrdt(a::Float64, b::Float64, f::Function)
	wid = (b - a) * 0.5
	cen = (b + a) * 0.5
	fc = f(cen)
	resg = fc * kronrod_wg[1]
	resk = fc * kronrod_wgk[1]
	for j = 1:11
	   	t = wid * kronrod_xgk[j + 1]
	   	fc = f(cen - t) + f(cen + t)
	   	resk += kronrod_wgk[j + 1] * fc
	   	if mod(j, 2) == 0; resg += kronrod_wg[div(j,2) + 1] * fc; end
	end
	retval = wid * resk
	err = abs(wid * (resk - resg))
	return retval, err
end

# # Student t Distribution Function

function tcdf(nu::Int, t::Float64)
	if nu < 1
   		studnt = cdf(Normal(), t)
	elseif nu == 1
   		studnt = (1.0 + 2.0atan(t) / pi) * 0.5
	elseif nu == 2
   		studnt = (1.0 + t / sqrt(2.0 + t*t)) * 0.5
	else
   		tt = t*t
   		cssthe = 1.0 / (1.0 + tt / nu)
   		polyn = 1.0
   		for j = nu-2:-2:2
    	  	polyn = 1.0 + (j - 1) * cssthe * polyn / j
   		end
   		if mod(nu, 2) == 1
    	  	rn = nu
    	  	ts = t / sqrt(rn)
    	  	studnt = (1.0 + 2.0*(atan(ts) + ts * cssthe * polyn) / pi) * 0.5
   		else
    	  	snthe = t / sqrt(nu + tt)
    	  	studnt = (1.0 + snthe * polyn) * 0.5
   		end
   		studnt = max(0.0, min(studnt, 1.0))
	end
	return studnt
end

function bvtcdf(nu::Int, dh::Float64, dk::Float64, r::Float64)
	if nu < 1
		return bvnuppercdf(-dh, -dk, r)
	elseif 1.0 - r < eps()
		return tcdf(nu, min(dh, dk))
	elseif r + 1.0 < eps()
		if dh > -dk
			return tcdf(nu, dh) - tcdf(nu, -dk)
		else
			return 0.0
		end
	else
		snu = sqrt(nu)
		ors = 1.0 - r * r
		hrk = dh - r * dk
		krh = dk - r * dh
		if abs(hrk) + ors > 0
			xnhk = hrk^2 / (hrk^2 + ors * (nu + dk^2))
			xnkh = krh^2 / (krh^2 + ors * (nu + dh^2))
		else
			xnhk = 0.0
			xnkh = 0.0
		end
		hs = sign(hrk)
		ks = sign(krh)
		if mod(nu, 2) == 0
			bvt = atan2(sqrt(ors), -r) / (2.0pi)
			gmph = dh / sqrt(16.0 * (nu + dh^2))
			gmpk = dk / sqrt(16.0 * (nu + dk^2))
			btnckh = 2.0atan2(sqrt(xnkh), sqrt(1.0 - xnkh)) / pi
            btpdkh = 2.0sqrt(xnkh * (1.0 - xnkh)) / pi
            btnchk = 2.0atan2(sqrt(xnhk), sqrt(1.0 - xnhk)) / pi
            btpdhk = 2.0sqrt(xnhk * (1.0 - xnhk)) / pi
            for j = 1:div(nu, 2)
               	bvt += gmph * (1.0 + ks * btnckh)
               	bvt += gmpk * (1.0 + hs * btnchk)
               	btnckh += btpdkh
               	btpdkh *= 2.0j * (1.0 - xnkh) / (2j + 1)
               	btnchk += btpdhk
               	btpdhk *= 2.0j * (1.0 - xnhk) / (2j + 1)
               	gmph *= (2j - 1) / (2j * (1.0 + dh * dh / nu))
               	gmpk *= (2j - 1) / (2j * (1.0 + dk * dk / nu))
            end
        else
            qhrk = sqrt(dh * dh + dk * dk - 2.0r*dh*dk + nu*ors)
            hkrn = dh*dk + r*nu
            hkn = dh*dk - nu
            hpk = dh + dk
            bvt = atan2(-snu * (hkn * qhrk + hpk * hkrn), hkn * hkrn - nu * hpk * qhrk) / (2.0pi)
            if bvt < -eps()
            	bvt += 1.0
            end
            gmph = dh / (2.0pi * snu * (1.0 + dh * dh / nu))
            gmpk = dk / (2.0pi * snu * (1.0 + dk * dk / nu))
            btnckh = sqrt(xnkh)
            btpdkh = btnckh
            btnchk = sqrt(xnhk)
            btpdhk = btnchk
            for j = 1:div(nu - 1, 2)
               	bvt += gmph * (1.0 + ks * btnckh)
               	bvt += gmpk * (1.0 + hs * btnchk)
               	btpdkh *= (2j - 1) * (1.0 - xnkh) / (2*j)
               	btnckh += btpdkh
               	btpdhk *= (2j - 1) * (1.0 - xnhk) / (2*j)
               	btnchk += btpdhk
               	gmph *= 2.0j / ((2j + 1) * (1.0 + dh * dh / nu ))
               	gmpk *= 2.0j / ((2j + 1) * (1.0 + dk * dk / nu ))
            end
        end
        return bvt
    end
end

# This function is based on the method described by
#     Drezner, Z and G.O. Wesolowsky, (1989),
#     On the computation of the bivariate normal integral,
#     Journal of Statist. Comput. Simul. 35, pp. 101-107,
# with major modifications for double precision, and for |R| close to 1.

const bvncdf_w_array = [0.1713244923791705e+00 0.4717533638651177e-01 0.1761400713915212e-01;
     					0.3607615730481384e+00 0.1069393259953183e+00 0.4060142980038694e-01;
     					0.4679139345726904e+00 0.1600783285433464e+00 0.6267204833410906e-01;
     					0.0 				   0.2031674267230659e+00 0.8327674157670475e-01;
     					0.0					   0.2334925365383547e+00 0.1019301198172404e+00;
     					0.0					   0.2491470458134029e+00 0.1181945319615184e+00;
     					0.0					   0.0					  0.1316886384491766e+00;
     					0.0					   0.0					  0.1420961093183821e+00;
     					0.0					   0.0					  0.1491729864726037e+00;
     					0.0					   0.0					  0.1527533871307259e+00]

const bvncdf_x_array = [-0.9324695142031522e+00 -0.9815606342467191e+00 -0.9931285991850949e+00;
						-0.6612093864662647e+00 -0.9041172563704750e+00 -0.9639719272779138e+00;
						-0.2386191860831970e+00 -0.7699026741943050e+00 -0.9122344282513259e+00;
						 0.0 				    -0.5873179542866171e+00 -0.8391169718222188e+00;
						 0.0 				    -0.3678314989981802e+00 -0.7463319064601508e+00;
						 0.0 				    -0.1252334085114692e+00 -0.6360536807265150e+00;
						 0.0 				    0.0 				    -0.5108670019508271e+00;
						 0.0 				    0.0 				    -0.3737060887154196e+00;
						 0.0 				    0.0 				    -0.2277858511416451e+00;
						 0.0 				    0.0 				    -0.7652652113349733e-01]

function bvnuppercdf(dh::Float64, dk::Float64, r::Float64)
	if abs(r) < 0.3
	   ng = 1
	   lg = 3
	elseif abs(r) < 0.75
	   ng = 2
	   lg = 6
	else
	   ng = 3
	   lg = 10
	end
	h = dh
	k = dk
	hk = h*k
	bvn = 0.0
	if abs(r) < 0.925
	   	if abs(r) > 0
	      	hs = (h * h + k * k) * 0.5
	      	asr = asin(r)
	      	for i = 1:lg
	         	for j = -1:2:1
	            	sn = sin(asr * (j * bvncdf_x_array[i, ng] + 1.0) * 0.5)
	            	bvn += bvncdf_w_array[i, ng] * exp((sn * hk - hs) / (1.0 - sn*sn))
	        	end
	      	end
	      	bvn *= asr / (4.0pi)
	   	end
	   	bvn += cdf(Normal(), -h) * cdf(Normal(), -k)
	else
	   	if r < 0
	      	k = -k
	      	hk = -hk
	   	end
	   	if abs(r) < 1
	      	as = (1.0 - r) * (1.0 + r)
	      	a = sqrt(as)
	      	bs = (h - k)^2
	      	c = (4.0 - hk) * 0.125
	      	d = (12.0 - hk) * 0.0625
	      	asr = -(bs / as + hk) * 0.5
	      	if ( asr .gt. -100 )
	      		bvn = a * exp(asr) * (1.0 - c * (bs - as) * (1.0 - d * bs / 5.0) / 3.0 + c * d * as * as / 5.0)
	      	end
	      	if -hk < 100
	         	b = sqrt(bs)
	         	bvn -= exp(-hk * 0.5) * sqrt(2.0pi) * cdf(Normal(), -b / a) * b * (1.0 - c * bs * (1.0 - d * bs / 5.0) / 3.0)
	      	end
	     	a /= 2.0
		    for i = 1:lg
	         	for j = -1:2:1
	            	xs = (a * (j*bvncdf_x_array[i, ng] + 1.0))^2
	            	rs = sqrt(1.0 - xs)
	            	asr = -(bs / xs + hk) * 0.5
	            	if asr > -100
	               		bvn += a * bvncdf_w_array[i, ng] * exp(asr) * (exp(-hk * (1.0 - rs) / (2.0 * (1.0 + rs))) / rs - (1.0 + c * xs * (1.0 + d * xs)))
	            	end
	         	end
	        end
	      	bvn /= -2.0pi
	   	end
	   	if r > 0
	      	bvn += cdf(Normal(), -max(h, k))
	   	else
	      	bvn = -bvn
	      	if k > h
	      		bvn += cdf(Normal(), k) - cdf(Normal(), h)
	      	end
		end
	end
	return bvn
end


text:
using DynamicalSystemsBase
using Test, StaticArrays
using DynamicalSystemsBase: CDS, DDS
using DynamicalSystemsBase.Systems: hoop, hoop_jac, hiip, hiip_jac
using DynamicalSystemsBase.Systems: loop, loop_jac, liip, liip_jac

u0 = [0, 10.0, 0]
p = [10, 28, 8/3]
u0h = zeros(2)
ph = [1.4, 0.3]

FUNCTIONS = [liip, liip_jac, loop, loop_jac, hiip, hiip_jac, hoop, hoop_jac]
INITCOD = [u0, u0h]
PARAMS = [p, ph]

#######################################################################################
#                                      TIMINGS                                        #
#######################################################################################
println("\nStarting timings (using @time)..")

println("\nAutodiff CDS (create/jacobian):")
println("IIP")
@time lip = CDS(liip, u0, p)
jacobian(lip); @time jacobian(lip)
println("OOP")
@time lop = CDS(loop, u0, p)
jacobian(lop); @time jacobian(lop)

println("\n CDS (create/jacobian):")
println("IIP")
@time lipjac = CDS(liip, u0, p, liip_jac)
jacobian(lipjac); @time jacobian(lipjac)
println("OOP")
@time lopjac = CDS(loop, u0, p, loop_jac)
jacobian(lopjac); @time jacobian(lopjac)

println("\nAutodiff DDS (create/jacobian):")
println("IIP")
@time hip = DDS(hiip, u0h, ph)
jacobian(hip); @time jacobian(hip)
println("OOP")
@time hop = DDS(hoop, u0h, ph)
jacobian(hop); @time jacobian(hop)

println("\nDDS (create/jacobian):")
println("IIP")
@time hipjac = DDS(hiip, u0h, ph, hiip_jac)
jacobian(hipjac); @time jacobian(hipjac)
println("OOP")
@time hopjac = DDS(hoop, u0h, ph, hoop_jac)
jacobian(hopjac); @time jacobian(hopjac)


text:
using Random
using Functors
using Zygote

include("grid.jl")
include("conv.jl")
include("radfuncs.jl")
Random.seed!(1)

mutable struct Op
    l
    kernel::AbstractArray
    grid::Grid
    convfunc
    kwargs
    radfunc
    rmin
    rmax
end
@functor Op
Flux.trainable(m::Op) = [m.radfunc]

function makekernel(radfunc, rmin, rmax, l, grid)
    @unpack cell, rhat,R, dv = grid
    n = size(cell, 1)
    f = r-> rmin <= r <= rmax ? radfunc(r) : 0.0
    rscalars = f.(R)
    if l == 0
        kernel = rscalars
    elseif l == 1
        kernel = rscalars.*rhat
    end
    dv * kernel
end
function Op(
    radfunc,
    rmin,
    rmax,
    cell::Matrix;
    l = 0,
    convfunc=conv,
    boundary =:same
)
    grid = Grid(cell, rmax)
    kernel = makekernel(radfunc, rmin, rmax, l, grid)
    kwargs=(;boundary)

    Op(l, kernel, grid,convfunc,kwargs, radfunc,rmin, rmax)
end

"""
    Op(
        name::Union{Symbol,String},
        cell;
        boundary =:same,
        rmin = 0,
        rmax = Inf,
        l = 0,
        σ = 1.0,
    )
    Op(
        radfunc,
        rmin::AbstractFloat,
        rmax::AbstractFloat,
        cell;
        l = 0,
        boundary =:same,
    )

`Op` constructs finite difference operators. Prebuilt operators like differential operators (`▽`) & common Green's functions can be specified by name. Custom equivariant operators can be made by specifying radial function.
"""
function Op(
    name::Union{Symbol,String},
    cell;
    rmin = 0.,
    rmax = Inf,
    convfunc=conv,
    boundary =:same,
    l = 0,
    σ = 1.0,
)
    name = Symbol(name)
    if name == :▽
        name = :grad
    end

    dims = size(cell, 1)
    radfunc=nothing
if name == :Gaussian
    radfunc = r -> exp(-r^2 / (2 * σ^2)) / sqrt(2π * σ^(2n))
    rmax=2σ
    return Op(radfunc,rmin,rmax,cell)
end

    if name == :neural
    elseif name == :grad
        boundary = :smooth
        l = 1
        grid = Grid(cell, fill(3, dims))

            kernel = [
                sum(abs.(v)) > 1 ? zeros(dims) :
                    -cell' \ collect(v) / 2 for v in Iterators.product(fill(-1:1,dims)...)
            ]
    end
    kwargs=(;boundary)
    Op(l, kernel, grid, convfunc,kwargs,radfunc,rmin, rmax)
end

"""
    function (m::Op)(x::AbstractArray, )

"""
function (m::Op)(x::AbstractArray,product=*)
    @unpack grid, kernel, convfunc,kwargs = m
    convfunc(x,kernel;product,kwargs...)
end

function LinearAlgebra.:⋅(m::Op, x)
    m(x,⋅)
end
function ⨉(m::Op, x)
    m(x,⨉)
end
function remake!(m)
    @unpack radfunc,rmin, rmax,  grid,l = m
    m.kernel = makekernel(radfunc,rmin, rmax, l, grid)
end
function Base.abs(x::AbstractArray  )
    sum(abs.(x))
end
function nae(yhat, y; sumy = sum(abs.(y)))
    if sumy == 0
        error()
    end
    sum(abs.(yhat .- y)) / sumy
end


text:
module TSAD
    using ..OutlierDetectionData: ARFF
    using OutlierDetection: CLASS_NORMAL, CLASS_OUTLIER, to_categorical
    using DataFrames: DataFrame, Not, select, convert, disallowmissing!, rename!
    using DataDeps
    using DataFrames

    export list, load

    include("datasets.jl")

    const TSAD_DATASETS = [
        (
            name = "univariate",
            md5 = "0eb0d3d5fd107ea0e92e1a76d25a662e73de93e6c6bcdfc0ef693246f10401f4",
            url = "http://www.timeseriesclassification.com/Downloads/Archives/Univariate2018_arff.zip",
        ),
        (
            name = "multivariate",
            md5 =  "37508dd78a3683ec7aaf26e787f4f5890cd20872539a0c4cb28bcfacc7a02536",
            url =  "http://www.timeseriesclassification.com/Downloads/Archives/Multivariate2018_arff.zip",
        )
    ]

    list() = vcat(TSAD_UNIVARIATE, TSAD_MULTIVARIATE)
    to_name(dataset::AbstractString) = "tsad-$dataset"
    dataset_not_found(dataset) = ArgumentError("$dataset was not found in the available datasets.")
    check_dataset(dataset) = dataset in list() || throw(dataset_not_found(dataset))

    # some datasets use the '"' character in the comments, which leads to CSV parsing errors
    univariate_read(dataset) = last(splitdir(dataset)) in ["Coffee_TRAIN.arff",
                                                           "Coffee_TEST.arff",
                                                           "OliveOil_TRAIN.arff",
                                                           "OliveOil_TEST.arff"] ?
                               ARFF.read(dataset, quotechar=UInt8(''')) |> DataFrame :
                               ARFF.read(dataset) |> DataFrame

    function load(dataset::AbstractString, anomaly_class)
        dep_path = to_dep(dataset)
        path_train = joinpath(dep_path, dataset, "$(dataset)_TRAIN.arff")
        path_test = joinpath(dep_path, dataset, "$(dataset)_TEST.arff")

        # load training and test data and combine
        data_train, data_test = univariate_read(path_train), univariate_read(path_test)
        data_train[!, :is_train] .= true
        data_test[!, :is_train] .= false
        data = vcat(data_train, data_test)

        # check if anomaly class is in the targets
        unique_targets = sort(unique(data.target))
        @assert anomaly_class in unique_targets "The anomaly class $anomaly_class was not found in "*
                                                "the targets $unique_targets."

        # prepare data ready to return
        labels = to_categorical(ifelse.(data.target .== anomaly_class, CLASS_OUTLIER, CLASS_NORMAL))
        X, y = select(data, Not([:target])), labels
        return X, y
    end

    function load(dataset::AbstractString)
        check_dataset(dataset)
        load(dataset, TSAD_UNIVARIATE_DEFAULTS[dataset])
    end

    function to_dep(dataset::AbstractString)
        if dataset in TSAD_UNIVARIATE
            dep = @datadep_str to_name("univariate")
            return "$dep/Univariate_arff"
        elseif dataset in TSAD_MULTIVARIATE
            dep = @datadep_str to_name("multivariate")
            return "$dep/Multivariate_arff"
        else
            throw(dataset_not_found(dataset))
        end
    end

    function __init__()
        for (name, md5, url) in TSAD_DATASETS
            register(DataDep(
                to_name(name),
                """Dataset: $name
                Collection: The UCR Time Series Archive
                Authors: Dau et al.
                Website: https://www.timeseriesclassification.com""",
                url,
                md5,
                post_fetch_method = unpack
            ))
        end
    end
end


text:
using BesselK, SpecialFunctions, ArbNumerics
include("shared.jl")

function wbesselk(v,x)
  try
    return BesselK._besselk(v,x)
  catch
    return NaN
  end
end

rbesselk(v,x) =  Float64(ArbNumerics.besselk(ArbFloat(v), ArbFloat(x)))
abesselk(v,x) =  SpecialFunctions.besselk(v, x)

const BASELINE = [rbesselk(z[1], z[2]) for z in Iterators.product(VGRID, XGRID)]
const AMOS     = [abesselk(z[1], z[2]) for z in Iterators.product(VGRID, XGRID)]
const OURSOL   = [wbesselk(z[1], z[2]) for z in Iterators.product(VGRID, XGRID)]
const TOLS_A   = atolfun.(zip(BASELINE, AMOS))
const TOLS_U   = atolfun.(zip(BASELINE, OURSOL))
const TOLS_AU  = rtolfun.(zip(AMOS,     OURSOL))
const RTOLS_A  = rtolfun.(zip(BASELINE, AMOS))
const RTOLS_U  = rtolfun.(zip(BASELINE, OURSOL))

@assert iszero(length(findall(isnan, OURSOL))) "There were NaNs in our attempts!"

# quick simple test to find the worst atol when B(x) < 1:
let ix = findall(x->x<=one(x), BASELINE)
  res = findmax(abs,  TOLS_U[ix])
  _ix = res[2]
  println("Worst atol when true besselk(v,x) <= 1: $(res[1])")
  println("True value:                             $(BASELINE[ix][_ix])")
end

gnuplot_save_matrix!("../plotdata/atols_amos.csv",      TOLS_A,  VGRID, XGRID)
gnuplot_save_matrix!("../plotdata/atols_ours.csv",      TOLS_U,  VGRID, XGRID)
gnuplot_save_matrix!("../plotdata/rtols_amos.csv",      RTOLS_A, VGRID, XGRID)
gnuplot_save_matrix!("../plotdata/rtols_ours.csv",      RTOLS_U, VGRID, XGRID)
gnuplot_save_matrix!("../plotdata/rtols_amos_ours.csv", TOLS_AU, VGRID, XGRID)


text:
function get_name(obj::TimeScale)
    return jcall(obj, "getName", JString, ())
end

function offset_from_tai(obj::TimeScale, arg0::AbsoluteDate)
    return jcall(obj, "offsetFromTAI", jdouble, (AbsoluteDate,), arg0)
end

function offset_from_tai(obj::TimeScale, arg0::FieldAbsoluteDate)
    return jcall(obj, "offsetFromTAI", RealFieldElement, (FieldAbsoluteDate,), arg0)
end


text:
import JuMP: GenericAffExpr

init_neurons(model::Model, layers::Vector{Layer})     = init_variables(model, layers, :Cont, include_input = true)
init_deltas(model::Model, layers::Vector{Layer})      = init_variables(model, layers, :Bin)
init_multipliers(model::Model, layers::Vector{Layer}) = init_variables(model, layers, :Cont)
# for reluplex:
init_forward_facing_vars(model::Model, layers::Vector{Layer}) = init_variables(model, layers, :Cont)[1:end-1] #discard the last one
# Allow ::Network input also (NOTE for legacy purposes mostly...)
init_neurons(m,     network::Network) = init_neurons(m, network.layers)
init_deltas(m,      network::Network) = init_deltas(m,  network.layers)
init_multipliers(m, network::Network) = init_multipliers(m, network.layers)

function init_variables(model::Model, layers::Vector{Layer}, vartype::Symbol; include_input::Bool = false)
    # TODO: only neurons get offset array
    vars = Vector{Vector{Variable}}(undef, length(layers))
    all_layers_n = n_nodes.(layers)

    if include_input
        # input to the first layer also gets variables
        # essentially an input constraint
        input_layer_n = size(first(layers).weights, 2)
        prepend!(all_layers_n, input_layer_n)
        push!(vars, Vector{Variable}())        # expand vars by one to account
    end

    for (i, n) in enumerate(all_layers_n)
        vars[i] = @variable(model, [1:n], category = vartype)
    end
    return vars
end

function symbolic_max(m::Model, a, b)
    aux = @variable(m)
    @constraint(m, aux >= a)
    @constraint(m, aux >= b)
    return aux
end
symbolic_max(a::Variable, b::Variable)                               = symbolic_max(a.m, a, b)
symbolic_max(a::GenericAffExpr, b::GenericAffExpr)                   = symbolic_max(first(a.vars).m, a, b)
symbolic_max(a::Array{<:GenericAffExpr}, b::Array{<:GenericAffExpr}) = symbolic_max.(a, b)

function symbolic_abs(m::Model, v)
    aux = @variable(m) #get an anonymous variable
    @constraint(m, aux >= 0)
    @constraint(m, aux >= v)
    @constraint(m, aux >= -v)
    return aux
end
symbolic_abs(v::Variable)                = symbolic_abs(v.m, v)
symbolic_abs(v::GenericAffExpr)          = symbolic_abs(first(v.vars).m, v)
symbolic_abs(v::Array{<:GenericAffExpr}) = symbolic_abs.(v)

function symbolic_infty_norm(m::Model, v::Array{<:GenericAffExpr})
    aux = @variable(m)
    @constraint(m, aux >= 0)
    @constraint(m, aux .>= v)
    @constraint(m, aux .>= -v)
    return aux
end
# # in general, default to symbolic_abs behavior:
# symbolic_infty_norm(v) = symbolic_abs(v)
# only Array{<:GenericAffExpr} is needed
symbolic_infty_norm(v::Array{<:GenericAffExpr}) = symbolic_infty_norm(first(first(v).vars).m, v)


text:
import java.io.IOException;
import java.io.InputStream;
import java.net.URL;

import javax.swing.JEditorPane;

public class CachedJEditorPane extends JEditorPane {
    private static final long serialVersionUID = -2704024459339695654L;
    private final URLCache cache;

    public CachedJEditorPane() {
        super();
        cache = new URLCache();
    }

    /**
     * Fetches a stream for the given URL, which is about to
     * be loaded by the setPage method.  By
     * default, this simply opens the URL and returns the
     * stream.  This can be reimplemented to do useful things
     * like fetch the stream from a cache, monitor the progress
     * of the stream, etc.
     *

     * This method is expected to have the the side effect of
     * establishing the content type, and therefore setting the
     * appropriate EditorKit to use for loading the stream.
     *

     * If this the stream was an http connection, redirects
     * will be followed and the resulting URL will be set as
     * the Document.StreamDescriptionProperty so that relative
     * URL's can be properly resolved.
     *
     * @param page  the URL of the page
     */
    protected InputStream getStream(URL page) throws IOException {
	//	InputStream is = cache.get(page);
        /*if (is != null) {
            System.out.println("Cache hit for " + page);
        }
        else {*/
	//            System.out.println("Cache miss for " + page);
	URL{L} pageCheck = page;
	InputStream{L} data = super.getStream(page);
	cache.put(page, data);
	    //    is = cache.get(page);
	    //}
	return data;
	    //        return is;
    }

    public InputStream resolveStream(URL page) throws IOException {
	return getStream(page);
    }
}

max_stars_repo_name: HarvardPL/cryptoerase

max_issues_repo_name: HarvardPL/cryptoerase

max_forks_repo_name: HarvardPL/cryptoerase