text:
# IntervalSets.jl currently has deprecation warnings on Julia 0.7

export Interval, .., lessthan, LessThan, greaterthan, GreaterThan

"""
    Interval(start, stop)

Construct a closed interval, which is a collection which contains (via `in`) elements
between `start` and `stop` (inclusive) according to `isless`. The collection is abstract
in nature and doesn't support iteration, indexing, etc.

Can be constructed via the `..` function, e.g. `1..3 === Interval(1, 3)`.

# Examples

```julia
julia> 2 in Interval(1, 3)
true

julia> 3 in Interval(1, 3)
true

julia> 4 in Interval(1, 3)
false
"""
struct Interval{T1, T2}
    start::T1
    stop::T2
end

"""
    ..(start, stop)
    start..stop

Constructs an `Interval(start, stop)`, which represents the closed interval between `start`
and `stop`. `Interval`s are abstract collections which support `in` but not iteration,
indexing, etc.

The interval includes both the `start` and `stop` poitns. To exclude the `start` or `stop`
from the `Interval`, use the `greaterthan` or `lessthan` function.

# Examples

```julia
julia> 2 in 1..3
true

julia> 3 in 1..3
true

julia> 4 in 1..3
false

julia> 1 in greaterthan(1)..3
false

julia> 3 in 1..lessthan(3)
false
"""
..(start::T1, stop::T2) where {T1, T2} = Interval{T1, T2}(start, stop)

function Base.in(x, interval::Interval)
    return !(isless(x, interval.start) || isless(interval.stop, x))
end

function Base.isempty(i::Interval)
    isless(i.stop, i.start)
end

function Base.:(==)(i1::Interval, i2::Interval) # isequal falls back to ==
    if isempty(i1)
        return isempty(i2)
    elseif isempty(i2)
        return false
    else
        return isequal(i1.start, i2.start) && isequal(i1.stop, i2.stop)
    end
end

function Base.isless(i1::Interval, i2::Interval) # < falls back to isless
    if isempty(i1)
        return false
    elseif isempty(i2)
        return true
    else
        return isless(i1.start, i2.start) || (isequal(i1.start, i2.start) && isless(i1.stop, i2.stop))
    end
end

function Base.hash(interval::Interval, h::UInt)
    h = hash(UInt === UInt64 ? 0x0c3a059de789f681 : 0x0c88d4c5, h)
    if isempty(interval)
        return h
    end
    return hash(interval.stop, hash(interval.start, h))
end

function Base.show(io::IO, interval::Interval)
    print(io, interval.start)
    print(io, "..")
    print(io, interval.stop)
end

function Base.intersect(i1::Interval, i2::Interval)
    max(i1.start, i2.start) .. min(i1.stop, i2.stop)
end

function Base.union(i1::Interval, i2::Interval)
    if isless(i1.start, i2.stop)
        error("Union of intervals is not contiguous")
    end
    return min(i1.start, i2.start) .. max(i1.stop, i2.stop)
end

function Base.issubset(i1::Interval, i2::Interval)
    !(isless(i1.start, i1.start) || isless(i2.stop, i1.stop))
end

function intersects(i1::Interval, i2::Interval)
    i1.start ∈ i2 | i1.stop ∈ i2
end

struct LessThan{T}
    value::T
end

function Base.show(io::IO, x::LessThan)
    print(io, "lessthan(")
    print(io, x.value)
    print(io, ")")
end

"""
    lessthan(x)

Return a value which is immediately smaller than `x`. Than value is almost, but not quite,
equal to `x` - there should be no other values (of any type) in between `x` and
`lessthan(x)` according to Julia's `isless` and `isequal` canonical total ordering.

Amongst other uses, this may be used to create `Interval`s that exclude the end point.

See also `greaterthan`.

# Examples

```julua
julia> isequal(lessthan(10), 10)
false

julia> isless(lessthan(10), 10)
true

julia> 10 ∈ 0..10
true

julia> 10 ∈ 0..lessthan(10)
false
```
"""
function lessthan(x)
    LessThan{typeof(x)}(x)
end

Base.isequal(x::LessThan, y::LessThan) = isequal(x.value, y.value)
Base.isequal(x::LessThan, y) = false
Base.isequal(x, y::LessThan) = false
Base.isequal(x::LessThan, y::Missing) = false
Base.isequal(x::Missing, y::LessThan) = false

Base.isless(x::LessThan, y::LessThan) = isless(x.value, y.value)
Base.isless(x::LessThan, y) = islessequal(x.value, y)
Base.isless(x, y::LessThan) = isless(x, y.value)
Base.isless(x::LessThan, y::Missing) = true
Base.isless(x::Missing, y::LessThan) = false

Base.hash(x::LessThan, h::UInt) = hash(x.value, hash(UInt === UInt64 ? 0x1f61aad02a1ec08b : 0xd6318b8a, h))

struct GreaterThan{T}
    value::T
end

function Base.show(io::IO, x::GreaterThan)
    print(io, "greaterthan(")
    print(io, x.value)
    print(io, ")")
end

"""
    greaterthan(x)

Return a value which is immediately larger than `x`. Than value is almost, but not quite,
equal to `x` - there should be no other values (of any type) in between `x` and
`greaterthan(x)` according to Julia's `isless` and `isequal` canonical total ordering.

Amongst other uses, this may be used to create `Interval`s that exclude the starting point.

See also `lessthan`.

# Examples

```julua
julia> isequal(10, greaterthan(10))
false

julia> isless(10, greaterthan(10))
true

julia> 0 ∈ 0..10
true

julia> 0 ∈ greaterthan(0)..10
false
```
"""
function greaterthan(x)
    GreaterThan{typeof(x)}(x)
end

Base.isequal(x::GreaterThan, y::GreaterThan) = isequal(x.value, y.value)
Base.isequal(x::GreaterThan, y) = false
Base.isequal(x, y::GreaterThan) = false
Base.isequal(x::GreaterThan, y::Missing) = false
Base.isequal(x::Missing, y::GreaterThan) = false
Base.isequal(x::GreaterThan, y::LessThan) = false
Base.isequal(x::LessThan, y::GreaterThan) = false

Base.isless(x::GreaterThan, y::GreaterThan) = isless(x.value, y.value)
Base.isless(x::GreaterThan, y) = isless(x.value, y)
Base.isless(x, y::GreaterThan) = islessequal(x, y.value)
Base.isless(x::GreaterThan, y::Missing) = true
Base.isless(x::Missing, y::GreaterThan) = false
Base.isless(x::GreaterThan, y::LessThan) = isless(x.value, y.value)
Base.isless(x::LessThan, y::GreaterThan) = islessequal(x.value, y.value)

Base.hash(x::GreaterThan, h::UInt) = hash(x.value, hash(UInt === UInt64 ? 0x4fda09326e00a582 : 0xe4028a11, h))


text:
"""
"""
forwardinterp

function forwardinterp(bc, n)
    opn = _forwardinterp.(bc, n)
    eye = Eye{Bool}.(n)
    _forwardinterp(opn, eye)
end

_forwardinterp(::Dirichlet, n::Int) =
    Bidiagonal(ones(n) / 2, ones(n-1) / 2, :U)

_forwardinterp(::Periodic, n::Int) =
    spdiagm(1-n => ones(1) / 2,
            0 => ones(n) / 2,
            1 => ones(n-1) / 2)

_forwardinterp(opn::NTuple{1}, eye::NTuple{1}) = opn

_forwardinterp(opn::NTuple{2,Any}, eye::NTuple{2,Any}) =
    ApplyArray(kron, eye[2], opn[1]),
    ApplyArray(kron, opn[2], eye[1])

_forwardinterp(opn::NTuple{3,Any}, eye::NTuple{3,Any}) =
    ApplyArray(kron, eye[3], eye[2], opn[1]),
    ApplyArray(kron, eye[3], opn[2], eye[1]),
    ApplyArray(kron, opn[3], eye[2], eye[1])

"""
"""
backwardinterp

function backwardinterp(bc, n)
    opn = _backwardinterp.(bc, n)
    eye = Eye{Bool}.(n)
    _backwardinterp(opn, eye)
end

_backwardinterp(::Dirichlet, n::Int) =
    Bidiagonal(ones(n) / 2, ones(n-1) / 2, :L)

_backwardinterp(::Periodic, n::Int) =
    spdiagm(-1 => ones(n-1) / 2,
            0 => ones(n) / 2,
            n - 1 => ones(1) / 2)

_backwardinterp(opn::NTuple{1}, eye::NTuple{1}) = opn

_backwardinterp(opn::NTuple{2,Any}, eye::NTuple{2,Any}) =
    ApplyArray(kron, eye[2], opn[1]),
    ApplyArray(kron, opn[2], eye[1])

_backwardinterp(opn::NTuple{3,Any}, eye::NTuple{3,Any}) =
    ApplyArray(kron, eye[3], eye[2], opn[1]),
    ApplyArray(kron, eye[3], opn[2], eye[1]),
    ApplyArray(kron, opn[3], eye[2], eye[1])

ext: jl


text:
export uniform_mutate, goldman_mutate, profiling_mutate

"create a child by randomly mutating genes"
function uniform_mutate(cfg::NamedTuple, ind::CGPInd; kwargs...)::CGPInd
    chromosome = copy(ind.chromosome)
    chance = rand(length(ind.chromosome))
    non_output = length(ind.chromosome) - length(ind.outputs)
    change = [chance[1:non_output] .<= cfg.m_rate;
              chance[(non_output+1):end] .<= cfg.out_m_rate]
    chromosome[change] = rand(sum(change))
    kwargs_dict = Dict(kwargs)
    if haskey(kwargs_dict, :init_function)
        return kwargs_dict[:init_function](cfg, chromosome)
    else
        return CGPInd(cfg, chromosome)
    end
end

"create a child that is structurally different from the parent"
function goldman_mutate(cfg::NamedTuple, ind::CGPInd; kwargs...)::CGPInd
    child = uniform_mutate(cfg, ind; kwargs...)
    while true
        if any(ind.outputs != child.outputs)
            return child
        else
            for i in eachindex(ind.nodes)
                if ind.nodes[i].active
                    if child.nodes[i].active
                        if (ind.nodes[i].f != child.nodes[i].f
                            || ind.nodes[i].x != child.nodes[i].x
                            || ind.nodes[i].y != child.nodes[i].y)
                            return child
                        end
                    else
                        return child
                    end
                end
            end
        end
        child = uniform_mutate(cfg, ind; kwargs...)
    end
    nothing
end

"create a child that gives different outputs based on the provided inputs (a 2D matrix of (n_in, n_samples))"
function profiling_mutate(cfg::NamedTuple, ind::CGPInd, inputs::AbstractArray; kwargs...)::CGPInd
    child = uniform_mutate(cfg, ind; kwargs...)
    while true
        for i in 1:size(inputs, 2)
            out_ind = process(ind, inputs[:, i])
            out_child = process(child, inputs[:, i])
            if any(out_ind .!= out_child)
                return child
            end
        end
        child = uniform_mutate(cfg, ind; kwargs...)
    end
    nothing
end


text:
"""
    AdjustQuasiGLM(model, ϕ; level)
Estimates dispersion parameter, adjusts original GLM to reflect the dispersion and returns results in a pretty DataFrame.
Usage:
```julia-repl
AdjustQuasiGLM(model, ϕ; level)
```
Arguments:
- `model` : The `GLM` model.
- `data` : The `DataFrame` containing data that was used as input to the model.
- `level` : The desired degree of confidence.
"""

function AdjustQuasiGLM(model::StatsModels.TableRegressionModel, data::DataFrame; level::Real=0.95)

    # Calculate Pearson residuals

    resids = PearsonResiduals(model, data)

    # Estimate dispersion parameter ϕ and take √ to convert to multiplier

    ϕ = √EstimateDispersionParameter(resids, model)

    # Correct standard errors and calculate updated test statistics, p-values, and confidence intervals

    CorrectedOutputs = coefarray(model, ϕ; level)
    levstr = isinteger(level * 100) ? string(Integer(level * 100)) : string(level * 100)
    header = (["Parameter", "Estimate", "Std. Error", "t value", "Pr(>|t|)", "Lower $levstr%", "Upper $levstr%"])

    #--------------------------------------------
    # Organise results in a neat coeftable format
    #--------------------------------------------

    # Table formatting

    ctf = TextFormat(
    up_right_corner = ' ',
    up_left_corner = ' ',
    bottom_left_corner = ' ',
    bottom_right_corner = ' ',
    up_intersection = '─',
    left_intersection = ' ',
    right_intersection = ' ',
    middle_intersection = '─',
    bottom_intersection = '─',
    column = ' ',
    hlines = [ :begin, :header, :end]
    )

    # Render table

    println("\nCoefficients:")
    CorrectedOutputsPretty = PrettyTables.pretty_table(CorrectedOutputs; header = header, tf = ctf)

    # Return results in a DataFrame for further use

    CorrectedOutputs = DataFrame(CorrectedOutputs, :auto)

    CorrectedOutputs = rename!(CorrectedOutputs, [:x1, :x2, :x3, :x4, :x5, :x6, :x7] .=>  [Symbol(header[1]), Symbol(header[2]), Symbol(header[3]), Symbol(header[4]), Symbol(header[5]), Symbol(header[6]), Symbol(header[7])])

    # Recode column types from `Any` to `String` for parameter names and `Float64` for values columns

    for i in 2:size(header, 1)
        CorrectedOutputs[!, i] = convert(Array{Float64, 1}, CorrectedOutputs[!, i])
    end

    CorrectedOutputs[!, 1] = convert(Array{String, 1}, CorrectedOutputs[!, 1])
    return CorrectedOutputs
end


text:
module Xu

#=
This is an implementation of the algorithm presented in

"Alternative mathematical technique to determine LS spectral terms"

by Xu Renjun and Dai Zhenwen, published in JPhysB, 2006.
doi:10.1088/0953-4075/39/16/007

In the following, \(S'=2S\in\mathbb{Z}\) and
\(M_S'=2M_S\in\mathbb{Z}\), as in the original article.

\[\begin{aligned}
X(N,\ell,S',L) = A(N,\ell,\ell,S',L) - A(N,\ell,\ell,S',L+1)\\
+A(N,\ell,\ell,S'+2,L+1) - A(N,\ell,\ell,S'+2,L)
\end{aligned}
\quad(1)\]

with \(A(N,\ell,\ell_b,M_S',M_L)\) obeying four different cases:

Case 1, \(M_S'=1\), \(|M_L|\leq\ell\), and \(N=1\):

\[A(1,\ell,\ell_b,1,M_L) = 1\qquad(2)\]

Case 2, \(\{M_S'\}={2-N,4-N,...,N-2}\), \(|M_L| \leq
f\left(\frac{N-M_S'}{2}-1\right)+f\left(\frac{N+M_S'}{2}-1\right)\),
and \(1<N\leq 2\ell+1\):

\[\begin{aligned}
A(N,\ell,\ell,M_S',M_L) =
\sum_{\displaystyle M_{L-}=\max\left\{-f\left(\frac{N-M_S'}{2}-1\right),M_L-f\left(\frac{N+M_S'}{2}-1\right)\right\}}
^{\displaystyle\min\left\{f\left(\frac{N-M_S'}{2}-1\right),M_L+f\left(\frac{N+M_S'}{2}-1\right)\right\}}\\
\Bigg\{A\left(\frac{N-M_S'}{2},\ell,\ell,\frac{N-M_S'}{2},M_{L-}\right)
\times
A\left(\frac{N+M_S'}{2},\ell,\ell,\frac{N+M_S'}{2},M_L-M_{L-}\right)\Bigg\}
\end{aligned}\quad (3)\]

Case 3, \(M_S'=N\), \(|M_L|\leq f(N-1)\), and \(1<N\leq 2\ell+1\):

\[A(N,\ell,\ell_b,N,M_L) =
\sum_{\displaystyle M_{L_I} = \left\lfloor{\frac{M_L-1}{N}+\frac{N+1}{2}}\right\rfloor}
^{\displaystyle\min\{\ell_b,M_L+f(N-2)\}}
A(N-1,\ell,M_{L_I}-1,N-1,M_L-M_{L_I})
\]

Case 4, else:
\[A(N,\ell,\ell_b,M_S',M_L) = 0 \qquad (5)\]

\[f(n)=\begin{cases}
\displaystyle\sum_{m=0}^n(\ell-m), & n\geq0\\
0, & n<0
\end{cases}\]

=#

f(n::I,ℓ::I) where {I<:Integer} = n >= 0 ? sum(ℓ-m for m in 0:n) : 0

function A(N::I, ℓ::I, ℓ_b::I, M_S′::I, M_L::I) where {I<:Integer}
    if M_S′ == 1 && abs(M_L) <= ℓ && N == 1
        1 # Case 1
    elseif 1 < N && N <= 2ℓ + 1
        # Cases 2 and 3
        # N ± M_S' is always an even number
        a = (N-M_S′) >> 1
        b = (N+M_S′) >> 1
        fa = f(a-1,ℓ)
        fb = f(b-1,ℓ)

        if M_S′ in 2-N:2:N-2 && abs(M_L) <= fa + fb
            mapreduce(+, max(-fa, M_L - fb):min(fa, M_L + fb)) do M_Lm
                A(a,ℓ,ℓ,a,M_Lm)*A(b,ℓ,ℓ,b,M_L-M_Lm)
            end
        elseif M_S′ == N && abs(M_L) <= f(N-1,ℓ)
            mapreduce(+, floor(Int, (M_L-1)//N + (N+1)//2):min(ℓ_b,M_L + f(N-2,ℓ))) do M_LI
                A(N - 1, ℓ, M_LI - 1, N - 1, M_L-M_LI)
            end
        else
            0 # Case 4
        end
    else
        0 # Case 4
    end
end

X(N::I, ℓ::I, S′::I, L::I) where {I<:Integer} = A(N,ℓ,ℓ,S′,L) - A(N,ℓ,ℓ,S′,L+1) + A(N,ℓ,ℓ,S′+2,L+1) - A(N,ℓ,ℓ,S′+2,L)

end


text:
using Secp256k1
using Base58: base58checkdecode

@testset "Transaction" begin
    raw_tx = hex2bytes("0100000001813f79011acb80925dfe69b3def355fe914bd1d96a3f5f71bf8303c6a989c7d1000000006b483045022100ed81ff192e75a3fd2304004dcadb746fa5e24c5031ccfcf21320b0277457c98f02207a986d955c6e0cb35d446a89d3f56100f4d7f67801c31967743a9c8e10615bed01210349fc4e631e3624a545de3f89f5d8684c7b8138bd94bdd531d2e213bf016b278afeffffff02a135ef01000000001976a914bc3b654dca7e56b04dca18f2566cdaf02e8d9ada88ac99c39800000000001976a9141c4bc762dd5423e332166702cb75f40df79fea1288ac19430600")
    @testset "Parsing" begin
        @testset "Version" begin
            stream = IOBuffer(raw_tx)
            tx = parse(stream)::Tx
            @test tx.version == 1
        end
        @testset "Inputs" begin
            stream = IOBuffer(raw_tx)
            tx = parse(stream)::Tx
            @test length(tx.tx_ins) == 1
            want = hex2bytes("d1c789a9c60383bf715f3f6ad9d14b91fe55f3deb369fe5d9280cb1a01793f81")
            @test tx.tx_ins[1].prev_tx == want
            @test tx.tx_ins[1].prev_index == 0
            want = hex2bytes("6b483045022100ed81ff192e75a3fd2304004dcadb746fa5e24c5031ccfcf21320b0277457c98f02207a986d955c6e0cb35d446a89d3f56100f4d7f67801c31967743a9c8e10615bed01210349fc4e631e3624a545de3f89f5d8684c7b8138bd94bdd531d2e213bf016b278a")
            @test Bitcoin.serialize(tx.tx_ins[1].script_sig) == want
            @test tx.tx_ins[1].sequence == 0xfffffffe
        end
        @testset "Outputs" begin
            stream = IOBuffer(raw_tx)
            tx = parse(stream)::Tx
            @test length(tx.tx_outs) == 2
            want = 32454049
            @test tx.tx_outs[1].amount == want
            want = hex2bytes("1976a914bc3b654dca7e56b04dca18f2566cdaf02e8d9ada88ac")
            @test Bitcoin.serialize(tx.tx_outs[1].script_pubkey) == want
            want = 10011545
            @test tx.tx_outs[2].amount == want
            want = hex2bytes("1976a9141c4bc762dd5423e332166702cb75f40df79fea1288ac")
            @test Bitcoin.serialize(tx.tx_outs[2].script_pubkey) == want
        end
        @testset "Locktime" begin
            stream = IOBuffer(raw_tx)
            tx = parse(stream)::Tx
            @test tx.locktime == 410393
        end
    end
    @testset "Serialize" begin
        stream = IOBuffer(raw_tx)
        tx = Bitcoin.parse_legacy(stream)
        @test serialize(tx) == raw_tx
    end
    @testset  "Input Value" begin
        tx_hash = "d1c789a9c60383bf715f3f6ad9d14b91fe55f3deb369fe5d9280cb1a01793f81"
        index = 0
        want = 42505594
        tx_in = Bitcoin.TxIn(hex2bytes(tx_hash), index)
        @test Bitcoin.value(tx_in) == want
    end
    @testset "Input PubKey" begin
        tx_hash = "d1c789a9c60383bf715f3f6ad9d14b91fe55f3deb369fe5d9280cb1a01793f81"
        index = 0
        want = hex2bytes("1976a914a802fc56c704ce87c42d7c92eb75e7896bdc41ae88ac")
        tx_in = Bitcoin.TxIn(hex2bytes(tx_hash), index)
        @test Bitcoin.serialize(Bitcoin.script_pubkey(tx_in)) == want
    end
    @testset "Fee" begin
        raw_tx = hex2bytes("0100000001813f79011acb80925dfe69b3def355fe914bd1d96a3f5f71bf8303c6a989c7d1000000006b483045022100ed81ff192e75a3fd2304004dcadb746fa5e24c5031ccfcf21320b0277457c98f02207a986d955c6e0cb35d446a89d3f56100f4d7f67801c31967743a9c8e10615bed01210349fc4e631e3624a545de3f89f5d8684c7b8138bd94bdd531d2e213bf016b278afeffffff02a135ef01000000001976a914bc3b654dca7e56b04dca18f2566cdaf02e8d9ada88ac99c39800000000001976a9141c4bc762dd5423e332166702cb75f40df79fea1288ac19430600")
        stream = IOBuffer(raw_tx)
        tx = parse(stream)::Tx
        @test Bitcoin.fee(tx) == 40000
        raw_tx = hex2bytes("010000000456919960ac691763688d3d3bcea9ad6ecaf875df5339e148a1fc61c6ed7a069e010000006a47304402204585bcdef85e6b1c6af5c2669d4830ff86e42dd205c0e089bc2a821657e951c002201024a10366077f87d6bce1f7100ad8cfa8a064b39d4e8fe4ea13a7b71aa8180f012102f0da57e85eec2934a82a585ea337ce2f4998b50ae699dd79f5880e253dafafb7feffffffeb8f51f4038dc17e6313cf831d4f02281c2a468bde0fafd37f1bf882729e7fd3000000006a47304402207899531a52d59a6de200179928ca900254a36b8dff8bb75f5f5d71b1cdc26125022008b422690b8461cb52c3cc30330b23d574351872b7c361e9aae3649071c1a7160121035d5c93d9ac96881f19ba1f686f15f009ded7c62efe85a872e6a19b43c15a2937feffffff567bf40595119d1bb8a3037c356efd56170b64cbcc160fb028fa10704b45d775000000006a47304402204c7c7818424c7f7911da6cddc59655a70af1cb5eaf17c69dadbfc74ffa0b662f02207599e08bc8023693ad4e9527dc42c34210f7a7d1d1ddfc8492b654a11e7620a0012102158b46fbdff65d0172b7989aec8850aa0dae49abfb84c81ae6e5b251a58ace5cfeffffffd63a5e6c16e620f86f375925b21cabaf736c779f88fd04dcad51d26690f7f345010000006a47304402200633ea0d3314bea0d95b3cd8dadb2ef79ea8331ffe1e61f762c0f6daea0fabde022029f23b3e9c30f080446150b23852028751635dcee2be669c2a1686a4b5edf304012103ffd6f4a67e94aba353a00882e563ff2722eb4cff0ad6006e86ee20dfe7520d55feffffff0251430f00000000001976a914ab0c0b2e98b1ab6dbf67d4750b0a56244948a87988ac005a6202000000001976a9143c82d7df364eb6c75be8c80df2b3eda8db57397088ac46430600")
        stream = IOBuffer(raw_tx)
        tx = parse(stream)::Tx
        @test Bitcoin.fee(tx) == 140500
    end
    @testset "Sig Hash" begin
        raw_tx = hex2bytes("0100000001813f79011acb80925dfe69b3def355fe914bd1d96a3f5f71bf8303c6a989c7d1000000006b483045022100ed81ff192e75a3fd2304004dcadb746fa5e24c5031ccfcf21320b0277457c98f02207a986d955c6e0cb35d446a89d3f56100f4d7f67801c31967743a9c8e10615bed01210349fc4e631e3624a545de3f89f5d8684c7b8138bd94bdd531d2e213bf016b278afeffffff02a135ef01000000001976a914bc3b654dca7e56b04dca18f2566cdaf02e8d9ada88ac99c39800000000001976a9141c4bc762dd5423e332166702cb75f40df79fea1288ac19430600")
        stream = IOBuffer(raw_tx)
        tx = parse(stream)::Tx
        want = hex2bytes("27e0c5994dec7824e56dec6b2fcb342eb7cdb0d0957c2fce9882f715e85d81a6")
        @test Bitcoin.sig_hash(tx, 0) == want
    end
    @testset "hash_prevouts" begin
        tx = get_tx("d869f854e1f8788bcff294cc83b280942a8c728de71eb709a2c29d10bfe21b7c", testnet=true)
        want = hex2bytes("7fe72cb4866759324e1f158d0482fcc49015b0987966d80624dcc2db24638999")
        @test Bitcoin.hash_prevouts(tx) == want
    end
    @testset "Verify" begin
        @testset "P2PKH" begin
            tx = get_tx("452c629d67e41baec3ac6f04fe744b4b9617f8f859c63b3002f8684e7a4fee03")
            @test verify(tx)
            tx = get_tx("5418099cc755cb9dd3ebc6cf1a7888ad53a1a3beb5a025bce89eb1bf7f1650a2", testnet=true)
            @test verify(tx)
        end
        @testset "P2SH" begin
            println("Testing P2SH...")
            tx = get_tx("46df1a9484d0a81d03ce0ee543ab6e1a23ed06175c104a178268fad381216c2b")
            @test verify(tx)
        end
        @testset "P2WPKH" begin
            println("Testing P2WPKH...")
            tx = get_tx("d869f854e1f8788bcff294cc83b280942a8c728de71eb709a2c29d10bfe21b7c", testnet=true)
            @test verify(tx)
        end
        @testset "P2SH-P2WPKH" begin
            println("Testing P2SH-P2WPKH...")
            tx = get_tx("c586389e5e4b3acb9d6c8be1c19ae8ab2795397633176f5a6442a261bbdefc3a")
            @test verify(tx)
        end
        @testset "P2WSH" begin
            println("Testing P2WSH...")
            tx = get_tx("78457666f82c28aa37b74b506745a7c7684dc7842a52a457b09f09446721e11c", testnet=true)
            @test verify(tx)
        end
        @testset "P2SH-P2WSH" begin
            println("Testing P2SH-P2WSH...")
            tx = get_tx("954f43dbb30ad8024981c07d1f5eb6c9fd461e2cf1760dd1283f052af746fc88", testnet=true)
            @test verify(tx)
        end
    end
    @testset "Sign Input" begin
        keypair = KeyPair{:ECDSA}(8675309)
        tx_ins = TxIn[]
        prev_tx = hex2bytes("0025bc3c0fa8b7eb55b9437fdbd016870d18e0df0ace7bc9864efc38414147c8")
        push!(tx_ins, TxIn(prev_tx, 0))
        tx_outs = TxOut[]
        h160 = base58checkdecode(b"mzx5YhAH9kNHtcN481u6WkjeHjYtVeKVh2")[2:end]
        push!(tx_outs, TxOut(Int(0.99 * 100000000), Bitcoin.p2pkh_script(h160)))
        h160 = base58checkdecode(b"mnrVtF8DWjMu839VW3rBfgYaAfKk8983Xf")[2:end]
        push!(tx_outs, TxOut(Int(0.1 * 100000000), Bitcoin.p2pkh_script(h160)))
        tx = Tx(1, tx_ins, tx_outs, 0, true)
        @test txsigninput(tx, 0, keypair)
    end
    @testset "Is CoinbaseTx" begin
        raw_tx = hex2bytes("01000000010000000000000000000000000000000000000000000000000000000000000000ffffffff5e03d71b07254d696e656420627920416e74506f6f6c20626a31312f4542312f4144362f43205914293101fabe6d6d678e2c8c34afc36896e7d9402824ed38e856676ee94bfdb0c6c4bcd8b2e5666a0400000000000000c7270000a5e00e00ffffffff01faf20b58000000001976a914338c84849423992471bffb1a54a8d9b1d69dc28a88ac00000000")
        stream = IOBuffer(raw_tx)
        tx = parse(stream)::Tx
        @test iscoinbase(tx)
    end
    @testset "coinbase_height" begin
        raw_tx = hex2bytes("01000000010000000000000000000000000000000000000000000000000000000000000000ffffffff5e03d71b07254d696e656420627920416e74506f6f6c20626a31312f4542312f4144362f43205914293101fabe6d6d678e2c8c34afc36896e7d9402824ed38e856676ee94bfdb0c6c4bcd8b2e5666a0400000000000000c7270000a5e00e00ffffffff01faf20b58000000001976a914338c84849423992471bffb1a54a8d9b1d69dc28a88ac00000000")
        stream = IOBuffer(raw_tx)
        tx = parse(stream)::Tx
        @test coinbase_height(tx) == 465879
        raw_tx = hex2bytes("0100000001813f79011acb80925dfe69b3def355fe914bd1d96a3f5f71bf8303c6a989c7d1000000006b483045022100ed81ff192e75a3fd2304004dcadb746fa5e24c5031ccfcf21320b0277457c98f02207a986d955c6e0cb35d446a89d3f56100f4d7f67801c31967743a9c8e10615bed01210349fc4e631e3624a545de3f89f5d8684c7b8138bd94bdd531d2e213bf016b278afeffffff02a135ef01000000001976a914bc3b654dca7e56b04dca18f2566cdaf02e8d9ada88ac99c39800000000001976a9141c4bc762dd5423e332166702cb75f40df79fea1288ac19430600")
        stream = IOBuffer(raw_tx)
        tx = parse(stream)::Tx
        @test coinbase_height(tx) == nothing
    end
end
@testset "SegWit Transactions" begin
    raw = hex2bytes("0100000000010115e180dc28a2327e687facc33f10f2a20da717e5548406f7ae8b4c811072f85601000000000fffffff0100b4f505000000001976a9141d7cd6c75c2e86f4cbf98eaed221b30bd9a0b92888ac02483045022100df7b7e5cda14ddf91290e02ea10786e03eb11ee36ec02dd862fe9a326bbcb7fd02203f5b4496b667e6e281cc654a2da9e4f08660c620a1051337fa8965f727eb19190121038262a6c6cec93c2d3ecd6c6072efea86d02ff8e3328bbd0242b20af3425990ac00000000")
    @testset "Parse Segwit" begin
        tx = Bitcoin.parse(IOBuffer(raw), false)::Tx
        @test tx.version == 1
        @test tx.flag == 0x01
        hex1 = "3045022100df7b7e5cda14ddf91290e02ea10786e03eb11ee36ec02dd862fe9a326bbcb7fd02203f5b4496b667e6e281cc654a2da9e4f08660c620a1051337fa8965f727eb191901"
        hex2 = "038262a6c6cec93c2d3ecd6c6072efea86d02ff8e3328bbd0242b20af3425990ac"
        @test tx.tx_ins[1].witness.instructions == Script([hex2bytes(hex1), hex2bytes(hex2)]).instructions
    end
    @testset "Serialize" begin
        tx = Bitcoin.parse(IOBuffer(raw), false)::Tx
        @test Bitcoin.serialize(tx) == raw
    end
end


text:
using ControlSystems
using PyPlot

"""
    plotPZ(H, color="b", markersize=5)

Plot the poles and zeros of a transfer function.
"""
function plotPZ(H, color="b", markersize=5)
    z, p, = zpkdata(H)

    # Plot x and o for poles and zeros, respectively
    plot(real(p[1]), imag(p[1]), color=color, linestyle="None", marker="x", markersize=markersize, markeredgewidth=0.5)
    plot(real(z[1]), imag(z[1]), color=color, linestyle="None", marker="o", markerfacecolor="none", markersize=markersize, markeredgewidth=0.5)

    # Draw unit circle, real axis and imag axis
    circle = exp.(im*2*pi*(0:0.01:1))
    plot(real(circle), imag(circle), linewidth=0.5, color=:darkorange)
    axis("equal")
    axis([-1, 1, -1, 1])
    plot([0, 0], [-2, 2], linestyle="dotted", linewidth=1, color=:black)
    plot([-2, 2], [0, 0], linestyle="dotted", linewidth=1, color=:black)
end

max_stars_repo_name: Mahmoud-Kharsa/DeltaSigma

max_issues_repo_name: Mahmoud-Kharsa/DeltaSigma

max_forks_repo_name: Mahmoud-Kharsa/DeltaSigma


text:
using SubHunt
using Base.Test

using POMDPs
using POMDPToolbox
using DiscreteValueIteration
using QMDP
using ParticleFilters

pomdp = SubHuntPOMDP()

# show(STDOUT, MIME("text/plain"), SubVis(pomdp))

rng = MersenneTwister(6)
policy = RandomPolicy(pomdp, rng=rng)

# solver = QMDPSolver()
# policy = solve(solver, pomdp, verbose=true)
#
# s = initial_state(pomdp, rng)
# @show value(policy, s)
# @show value(policy, SubState(s.own, s.target, s.goal, true))

filter = SIRParticleFilter(pomdp, 10000, rng=rng)

for (s, a, r, sp) in stepthrough(pomdp, policy, filter, "s,a,r,sp", max_steps=200, rng=rng)
    v = SubVis(pomdp, s=s, a=a, r=r)
    show(STDOUT, MIME("text/plain"), v)
end

dpomdp = DSubHuntPOMDP(SubHuntPOMDP(), 1.0)

policy = RandomPolicy(dpomdp, rng=rng)
filter = SIRParticleFilter(dpomdp, 10000, rng=rng)

for (s, a, r, sp) in stepthrough(dpomdp, policy, filter, "s,a,r,sp", max_steps=200, rng=rng)
    v = SubVis(dpomdp.cp, s=s, a=a, r=r)
    show(STDOUT, MIME("text/plain"), v)
end


text:
"""
    module Molecules

Holds information of molecules.

Main reason form this module is to implement identical atom information when fitting potential
"""
module Molecules

export AbstractMolecule,
       Molecule,
       MoleculeIdenticalInformation,
       makeidentical!

using ..IdenticalTools
using ..Atoms

abstract type AbstractMolecule end

"""
    Molecule{T<:AbstractAtom} <: AbstractMolecule

Molecule representation

# Field
- `atoms::Vector{T}` : Contains information of atoms
"""
struct Molecule{T<:AbstractAtom} <: AbstractMolecule
    atoms::Vector{T}
    function Molecule{T}(atoms::AbstractVector{T}) where T <:AbstractAtom
         new(atoms)
    end
    function Molecule{T}(atomnames) where T <:AbstractAtom
         new(convert.(T,atomnames))
    end
end

"""
    Molecule{T<:AbstractAtom} <: AbstractMolecule

Molecule representation with information of identical atoms

# Field
- `atoms::Vector{T}` : Contains information of atoms
- `identical::Identical` : identical atom information
"""
struct MoleculeIdenticalInformation{T<:AbstractAtom} <: AbstractMolecule
    atoms::Vector{T}
    identical::Identical
    function MoleculeIdenticalInformation{T}(atoms::Vector{T}) where T <:AbstractAtom
         new(atoms,Identical(length(atoms)))
    end
    function MoleculeIdenticalInformation{T}(atomnames) where T <:AbstractAtom
         new(convert.(T,atomnames),Identical(length(atomnames)))
    end
end

"""
    makeidentical!(Mol::MoleculeIdenticalInformation,x)

Adds information for identical atoms to molecule.

# Arguments
- `Mol::MoleculeIdenticalInformation` : moleculet to which identical information is added
- `x` : collection of indices of identical atoms

## Throws
Error if `x` is out of bounds of `Mol`
"""
function makeidentical!(Mol::MoleculeIdenticalInformation,x)
    if maximum(x) <= length(Mol.atoms) && minimum(x) >= 1
        push!(Mol.identical,x)
    else
        error("Molecule has only $(length(Mol.atoms)) atoms")
    end
end

"""
    areidentical(mol::MoleculeIdenticalInformation,x)

Checks if atoms with indices given in collection `x` are identical
"""
IdenticalTools.areidentical(mol::MoleculeIdenticalInformation,x) = IdenticalTools.areidentical(mol.identical,x)

Base.convert(t::Type{<:AbstractMolecule}, m::AbstractMolecule) = t(m.atoms)

Base.length(mol::AbstractMolecule) = length(mol.atoms)

Base.show(io::IO, mol::AbstractMolecule) = print(io,"Molecule of $(length(mol)) atoms")

end #module


text:
# using Pkg
# Pkg.test("Term",coverage=true)

using Coverage
coverage = process_folder() # defaults to src/; alternatively, supply the folder name as argument
covered_lines, total_lines = get_summary(coverage)

import Term: Panel

res = "[blue][green]Covered[/green][white bold]/[/white bold][white]total[/white]: [green bold]$covered_lines[/green bold][white bold]/[/white bold][white bold]$total_lines[/red bold] ([yellow underline]$(round(covered_lines/total_lines, digits=3)* 100)%[/yellow underline])"

print("\n\n")
print(
    Panel(
        res,
        width=44,
        title="Term.jl",
        style="bold",
        justify=:center,
    )
)


text:
using Peacock
using PyPlot
#using Peacock.FDFD

# Permittivity
mu1 = [14 12.4im 0; -12.4im 14 0; 0 0 1]
function epf(x,y)
    d = 1/2.8
    r0 = 0.12
    # equation of a circle with radius 0.2a
    if (x - d)^2 + (y - d)^2 < r0^2 || (x - d)^2 + (y + d)^2 < r0^2 || (x + d)^2 + (y - d)^2 < r0^2 || (x + d)^2 + (y + d)^2 < r0^2
        #    if x^2+y^2 < r0^2
        # dielectric inside the circle
        return 15
    else
        # air outside the circle
        return 1
    end
end

# Permeability is unity everywhere
function muf(x,y)
    d = 1/2.8
    r0 = 0.12
    # equation of a circle with radius 0.2a
    if (x - d)^2 + (y - d)^2 < r0^2 || (x - d)^2 + (y + d)^2 < r0^2 || (x + d)^2 + (y - d)^2 < r0^2 || (x + d)^2 + (y + d)^2 < r0^2
        # dielectric inside the circle
        return mu1
    else
        # air outside the circle
        return 1
    end
end

a1 = [1, 0]  # first lattice vector
a2 = [0, 1]  # second lattice vector
#a2 = 0.5*[sqrt(3), 1]

d1 = 0.005  # resolution along first lattice vector
d2 = 0.005  # resolution along second lattice vector

G = BrillouinZoneCoordinate(  0,   0, "Γ")
X = BrillouinZoneCoordinate(1/2,   0, "X")
M = BrillouinZoneCoordinate(1/2, 1/2, "M")
ks = [G,X,M,G]

figure(figsize=(4,3))
geometry_FDFD = Geometry(epf, muf, a1, a2, d1, d2, TE)
solver_FDFD = Peacock.FDFD.Solver(geometry_FDFD, [2*d1, 2*d2])
plot_band_diagram(solver_FDFD, ks, TE, color="red",
            bands=1:4, dk=0.4, frequency_scale=1/2pi)

geometry_FDFD = Geometry(epf, muf, a1, a2, d1, d2, TM)
solver_FDFD = Peacock.FDFD.Solver(geometry_FDFD, [2*d1, 2*d2])
plot_band_diagram(solver_FDFD, ks, TM, color="blue",
            bands=1:4, dk=0.4, frequency_scale=1/2pi)
ylim(0,0.8)

geometry_PWE = Geometry(epf, muf, a1, a2, d1, d2)
solver_PWE =  Solver(geometry_PWE, 7)
figure(figsize=(4,3))
plot_band_diagram(solver_PWE, ks, TE, color="red",
            bands=1:4, dk=0.1, frequency_scale=1/2pi)
plot_band_diagram(solver_PWE, ks, TM, color="blue";
            bands=1:4, dk=0.1, frequency_scale=1/2pi)
ylim(0,0.8)


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

using BAT
using Test
using Random

@testset "rng_init" begin
    @testset "BAT.AbstractRNGSeed" begin
        philox = @inferred BAT.AbstractRNGSeed()
        @test typeof(philox) <: BAT.Philox4xSeed
    end

    @testset "BAT.Philox4xSeed" begin
        philox = @inferred BAT.Philox4xSeed()
        @test typeof(philox) <: BAT.AbstractRNGSeed
        @test typeof(philox.seed) <: NTuple{2, UInt64}
        philox = @inferred BAT.Philox4xSeed{UInt32}(tuple(ones(UInt32,2)...))
        @test typeof(philox.seed) <: NTuple{2, UInt32}
        philox = @inferred BAT.Philox4xSeed{UInt64}(tuple(UInt64(1),UInt64(2)))
        res = @inferred philox()
        @test res.key1 == UInt64(1) && res.key2 == UInt64(2)
    end

    @testset "BAT.Threefry4xSeed" begin
        tf4x = @inferred BAT.Threefry4xSeed()
        @test typeof(tf4x) <: BAT.AbstractRNGSeed
        @test typeof(tf4x.seed) <: NTuple{4, UInt64}
        tf4x = @inferred BAT.Threefry4xSeed{UInt32}(tuple(ones(UInt32,4)...))
        @test typeof(tf4x.seed) <: NTuple{4, UInt32}
        tf4x = @inferred BAT.Threefry4xSeed{UInt64}(tuple(UInt64(2):UInt64(5)...))
        res = @inferred tf4x()
        @test res.key1 == UInt64(2) && res.key2 == UInt64(3) &&
              res.key3 == UInt64(4) && res.key4 == UInt64(5)
    end

    @testset "BAT.reset_rng_counters!" begin
        philoxS = @inferred BAT.Philox4xSeed()
        philox = @inferred philoxS()
        BAT.reset_rng_counters!(philox, (1, 2, 3))
        @test philox.ctr4 == 1 && philox.ctr3 == 2 && philox.ctr2 == 3
        BAT.reset_rng_counters!(philox, 4, 5, 6)
        @test philox.ctr4 == 4 && philox.ctr3 == 5 && philox.ctr2 == 6
    end

    @testset "BAT.MersenneTwisterSeed" begin
        mertwS = @inferred BAT.MersenneTwisterSeed()
        @test typeof(mertwS) <: BAT.AbstractRNGSeed
        @test typeof(@inferred mertwS()) <: MersenneTwister
    end
end


text:
using DataStructures: OrderedDict, PriorityQueue
import DataStructures: peek, dequeue!

import Base: length, isempty, push!

"""
    MemoryJobStore()

`MemoryJobStore` implements `AbstractJobStore`.

This is a data structure which is responsible of
storing into memory jobs that should be executed later.
"""
struct MemoryJobStore <: AbstractJobStore
    d::OrderedDict
    q::PriorityQueue

    MemoryJobStore() = new(
        OrderedDict{String,Job}(),
        PriorityQueue{String, Priority}(Base.Order.Reverse)
    )
end

"""
    length(jobstore)

Return the number of jobs stored into jobstore.

A given job can be executed several times
(periodical jobs for example).
"""
length(jobstore::MemoryJobStore) = length(jobstore.d)

"""
    isempty(jobstore)

Returns `true` when there isn't any job in `jobstore`.
"""
isempty(jobstore::MemoryJobStore) = length(jobstore) == 0
function push!(jobstore::MemoryJobStore, job::Job)
    if hasjob(jobstore, job.id)
        error("Can't add a job with same job_id")
    else
        schedule(jobstore, job, true)
    end
end

"""
    schedule(jobstore, job, first_time)

Schedule a job according its `Priority`.
"""
function schedule(jobstore::MemoryJobStore, job::Job, first_time::Bool)
    now_ = job.dt_updated
    # Calculate the next fire time if there is none defined
    if !hasnextfire(job) || !first_time
        job.dt_next_fire = get_next_dt_fire(job.trigger, DateTime(0), now_)
    end
    jobstore.d[job.id] = job
    priority = Priority(job.dt_next_fire, job.priority)
    jobstore.q[job.id] = priority
    println("[+] $now_ Schedule job $(job.id)")
    dump(job)
    println()
end

"""
    job = peek(jobstore)

Return the next `job` from a `jobstore` without removing it from this `jobstore`.
"""
function peek(jobstore::MemoryJobStore)
    job_id, priority = peek(jobstore.q)
    job = jobstore.d[job_id]
    job
end

"""
    job = dequeue!(jobstore)

Remove and return the next `job` from a `jobstore`.
"""
function dequeue!(jobstore::MemoryJobStore)
    job_id = dequeue!(jobstore.q)
    jobstore.d[job_id]
end

"""
    hasjob(jobstore, job_id)

Return `true` if a job with identifier `job_id`
is stored into `jobstore`.
"""
function hasjob(jobstore::MemoryJobStore, job_id::String)
    job_id in keys(jobstore.d)
end

"""
    get_job_id(jobstore)

Returns a job identifier (a `job_id`).

Parameter `jobstore` is passed to ensure that
no other job stored in `jobstore` have same `job_id`.
"""
function get_job_id(jobstore::MemoryJobStore)
    job_id = ""
    while job_id == "" || hasjob(jobstore, job_id)
        job_id = get_job_id()
    end
    job_id
end

"""
    update(jobstore, job, now_)

Update `jobstore` for a given `job`.
Schedule when a `job` should be fire or delete it from `jobstore`
if it shouldn't be fired again.
"""
function update(jobstore::MemoryJobStore, job::Job, now_)
    job.dt_updated = now_
    itr_size = IteratorSize(job.trigger)
    if itr_size == HasLength()
        remaining = length(job.trigger) - job.n_triggered
        if remaining > 0
            schedule(jobstore, job, false)
        elseif remaining == 0
            delete!(jobstore.d, job.id)
        else
            error("remaining should be positive or no")
        end
    elseif itr_size == IsInfinite()
        schedule(jobstore, job, false)
    else
        error("undefined iterator size $itr_size")
    end
end


text:
function merge_groups(g::Symbol...)
    return unique(vcat(getindex.(Ref(parameter_groups), g)...))
end

# initialize parameter groups with all basic sets and then add combined ones
const parameter_groups = Dict(
    :standard => ["molec_id", "local_iso_id", "nu", "sw", "a", "elower", "gamma_air",
                "delta_air","gamma_self","n_air","n_self","gp","gpp"],

    :labels => ["statep", "statepp"],

    :voigt_air => ["gamma_air","delta_air","deltap_air","n_air"],
    :voigt_self => ["gamma_self","delta_self","deltap_self","n_self"],
    :voigt_h2 => ["gamma_H2","delta_H2","deltap_H2","n_H2"],
    :voigt_co2 => ["gamma_CO2","delta_CO2","n_CO2"],
    :voigt_he => ["gamma_He","delta_He","n_He"],
    :voigt_h2o => ["gamma_H2O","n_H2O"],
    :voigt_linemixing => ["y_air","y_self"],

    :sdvoigt_air => ["gamma_air","delta_air","deltap_air","n_air", "SD_air"],
    :sdvoigt_self => ["gamma_self","delta_self","deltap_self","n_self", "SD_self"],
    :sdvoigt_linemixing => ["Y_SDV_air_296","Y_SDV_self_296"],

    :ht_self => ["gamma_HT_0_self_50","n_HT_self_50","gamma_HT_2_self_50",
                "delta_HT_0_self_50","deltap_HT_self_50","delta_HT_2_self_50",
                "gamma_HT_0_self_150","n_HT_self_150","gamma_HT_2_self_150",
                "delta_HT_0_self_150","deltap_HT_self_150","delta_HT_2_self_150",
                "gamma_HT_0_self_296","n_HT_self_296","gamma_HT_2_self_296",
                "delta_HT_0_self_296","deltap_HT_self_296","delta_HT_2_self_296",
                "gamma_HT_0_self_700","n_HT_self_700","gamma_HT_2_self_700",
                "delta_HT_0_self_700","deltap_HT_self_700","delta_HT_2_self_700",
                "nu_HT_self","kappa_HT_self","eta_HT_self","Y_HT_self_296"],
    :ht_air => ["gamma_HT_0_air_296","n_HT_air_296","gamma_HT_2_air_296",
                "delta_HT_0_air_296","deltap_HT_air_296","delta_HT_2_air_296",
                "nu_HT_air","kappa_HT_air","eta_HT_air","Y_HT_air_296"]

)
parameter_groups[:voigt_all] = merge_groups(:voigt_air, :voigt_self, :voigt_h2, :voigt_co2, :voigt_he, :voigt_h2o, :voigt_linemixing)
parameter_groups[:sdvoigt_all] = merge_groups(:sdvoigt_air, :sdvoigt_self, :sdvoigt_linemixing)
parameter_groups[:ht_all] = merge_groups(:ht_self, :ht_air)

parameter_groups[:all] = merge_groups(:standard, :labels, :voigt_all, :sdvoigt_all, :ht_all)


text:
using SearchLight, Chirps

### Your tests here
@test 1 == 1

@show "foo bar"

max_stars_repo_name: essenciary/chirper

max_issues_repo_name: essenciary/chirper

max_forks_repo_name: essenciary/chirper


text:
# epilogue.jl

function cxxrtl_get(handle, name)
    parts = Ref(Csize_t(0));
	object = cxxrtl_get_parts(handle, name, parts)
	@assert(object == C_NULL || parts[] == 1)
	if (object == C_NULL || parts[] == 1)
		return object
    end
	return C_NULL
end

function cxxrtl_design_create(lib)
    global libcxxrtl = lib
    ccall((:cxxrtl_design_create, libcxxrtl), cxxrtl_toplevel, ())
end

function synth(output, files...)
    cmd = Yosys_jll.yosys()
    push!(cmd.exec, "-q", "-p", "write_cxxrtl $(output)", files...)
    run(cmd)
end

function compile(output, mod)
    cmd = Clang_jll.clang()
    push!(cmd.exec, "-g", "-O3", "-fPIC", "-shared", "-std=c++14",
                    "-I$(inc)", capi, mod, "-o", output)
    run(cmd)
end

include("cosim.jl")
# end

max_stars_repo_name: JuliaComputing/CxxrtlCosim

max_issues_repo_name: JuliaComputing/CxxrtlCosim

max_forks_repo_name: JuliaComputing/CxxrtlCosim


text:
export BlockSparse,
       BlockSparseTensor,
       Block,
       nzblock,
       BlockOffset,
       BlockOffsets,
       blockoffsets,
       blockview,
       nnzblocks,
       nzblocks,
       nnz,
       findblock,
       isblocknz

#
# BlockOffsets
#

const Block{N} = NTuple{N,Int}
const Blocks{N} = Vector{Block{N}}
const BlockOffset{N} = Pair{Block{N},Int}
const BlockOffsets{N} = Vector{BlockOffset{N}}

BlockOffset(block::Block{N},
            offset::Int) where {N} = BlockOffset{N}(block, offset)

nzblock(bof::BlockOffset) = first(bof)

offset(bof::BlockOffset) = last(bof)

nzblock(block::Block) = block

# Get the offset if the nth block in the block-offsets
# list
offset(bofs::BlockOffsets,
       n::Int) = offset(bofs[n])

# TODO: rename nzblock?
nzblock(bofs::BlockOffsets,
        n::Int) = nzblock(bofs[n])

nnzblocks(bofs::BlockOffsets) = length(bofs)
nnzblocks(bs::Blocks) = length(bs)

# TODO: make an iterator eachnzblocks to avoid allocation
function nzblocks(bofs::BlockOffsets{N}) where {N}
  blocks = Blocks{N}(undef, nnzblocks(bofs))
  for i in 1:nnzblocks(bofs)
    blocks[i] = nzblock(bofs, i)
  end
  return blocks
end

# define block ordering with reverse lexographical order
function isblockless(b1::Block{N},
                     b2::Block{N}) where {N}
  return CartesianIndex(b1) < CartesianIndex(b2)
end

function isblockless(bof1::BlockOffset{N},
                     bof2::BlockOffset{N}) where {N}
  return isblockless(nzblock(bof1), nzblock(bof2))
end

function isblockless(bof1::BlockOffset{N},
                     b2::Block{N}) where {N}
  return isblockless(nzblock(bof1), b2)
end

function isblockless(b1::Block{N},
                     bof2::BlockOffset{N}) where {N}
  return isblockless(b1, nzblock(bof2))
end

function check_blocks_sorted(blockoffsets::BlockOffsets)
  for jj in 1:length(blockoffsets)-1
    block_jj = nzblock(blockoffsets[jj])
    block_jj1 = nzblock(blockoffsets[jj+1])
    if !isblockless(block_jj,block_jj1)
      error("Blocks in BlockOffsets not ordered")
    end
  end
  return
end

function offset(bofs::BlockOffsets{N},
                block::Block{N}) where {N}
  block_pos = findblock(bofs,block)
  isnothing(block_pos) && return nothing
  return offset(bofs,block_pos)
end

"""
findblock(::BlockOffsets,::Block)

Output the index of the specified block in the block-offsets
list.
If not found, return nothing.
Searches assuming the blocks are sorted.
If more than one block exists, throw an error.
"""
function findblock(bofs::BlockOffsets{N},
                   find_block::Block{N}; sorted=true) where {N}
  r = sorted ? searchsorted(bofs,find_block;lt=isblockless) : findall(i->block(i)==find_block,bofs)
  length(r)>1 && error("In findblock, more than one block found")
  length(r)==0 && return nothing
  return first(r)
end

function nnz(bofs::BlockOffsets,inds)
  nnzblocks(bofs) == 0 && return 0
  lastblock,lastoffset = bofs[end]
  return lastoffset + blockdim(inds,lastblock)
end

"""
new_block_pos(::BlockOffsets,::Block)

Output the index where the specified block should go in
the block-offsets list.
Searches assuming the blocks are sorted.
If the block already exists, throw an error.
"""
function new_block_pos(bofs::BlockOffsets{N},
                       block::Block{N}) where {N}
  r = searchsorted(bofs,block;lt=isblockless)
  length(r)>1 && error("In new_block_pos, more than one block found")
  length(r)==1 && error("In new_block_pos, block already found")
  return first(r)
end

# TODO: should this be a constructor?
function blockoffsets(blocks::Blocks{N},
                      inds; sorted = true) where {N}
  if sorted
    blocks = sort(blocks;lt=isblockless)
  end
  blockoffsets = BlockOffsets{N}(undef,length(blocks))
  nnz = 0
  for (i,block) in enumerate(blocks)
    blockoffsets[i] = block=>nnz
    current_block_dim = blockdim(inds,block)
    nnz += current_block_dim
  end
  return blockoffsets,nnz
end

"""
diagblockoffsets(blocks::Blocks,inds)

Get the blockoffsets only along the diagonal.
The offsets are along the diagonal.
"""
function diagblockoffsets(blocks::Blocks{N},
                          inds) where {N}
  blocks = sort(blocks;lt=isblockless)
  blockoffsets = BlockOffsets{N}(undef,length(blocks))
  nnzdiag = 0
  for (i,block) in enumerate(blocks)
    blockoffsets[i] = block=>nnzdiag
    current_block_diaglength = blockdiaglength(inds,block)
    nnzdiag += current_block_diaglength
  end
  return blockoffsets,nnzdiag
end

# Permute the blockoffsets and indices
function Base.permutedims(boffs::BlockOffsets{N},
                          inds,
                          perm::NTuple{N,Int}) where {N}
  blocksR = Blocks{N}(undef,nnzblocks(boffs))
  for (i,(block,offset)) in enumerate(boffs)
    blocksR[i] = permute(block,perm)
  end
  indsR = permute(inds,perm)
  blockoffsetsR,_ = blockoffsets(blocksR,indsR)
  return blockoffsetsR,indsR
end

function Base.permutedims(blocks::Blocks{N},
                          perm::NTuple{N,Int}) where {N}
  blocks_perm = Blocks{N}(undef,nnzblocks(blocks))
  for (i,block) in enumerate(blocks)
    blocks_perm[i] = permute(block,perm)
  end
  return blocks_perm
end

"""
blockdim(T::BlockOffsets,nnz::Int,pos::Int)

Get the block dimension of the block at position pos.
"""
function blockdim(boffs::BlockOffsets,
                  nnz::Int,
                  pos::Int)
  if nnzblocks(boffs)==0
    return 0
  elseif pos==nnzblocks(boffs)
    return nnz-offset(boffs,pos)
  end
  return offset(boffs,pos+1)-offset(boffs,pos)
end

function Base.union(boffs1::BlockOffsets{N},
                    nnz1::Int,
                    boffs2::BlockOffsets{N},
                    nnz2::Int) where {N}
  n1,n2 = 1,1
  boffsR = BlockOffset{N}[]
  current_offset = 0
  while n1 <= length(boffs1) && n2 <= length(boffs2)
    if isblockless(boffs1[n1],boffs2[n2])
      push!(boffsR, BlockOffset(nzblock(boffs1[n1]),
                                current_offset))
      current_offset += blockdim(boffs1, nnz1, n1)
      n1 += 1
    elseif isblockless(boffs2[n2], boffs1[n1])
      push!(boffsR, BlockOffset(nzblock(boffs2[n2]),
                                current_offset))
      current_offset += blockdim(boffs2, nnz2, n2)
      n2 += 1
    else
      push!(boffsR, BlockOffset(nzblock(boffs1[n1]),
                                current_offset))
      current_offset += blockdim(boffs1, nnz1, n1)
      n1 += 1
      n2 += 1
    end
  end
  if n1 <= length(boffs1)
    for n in n1:length(boffs1)
      push!(boffsR, BlockOffset(nzblock(boffs1[n]),
                                current_offset))
      current_offset += blockdim(boffs1, nnz1, n)
    end
  elseif n2 <= length(boffs2)
    for n in n2:length(bofss2)
      push!(boffsR, BlockOffset(nzblock(boffs2[n]),
                                current_offset))
      current_offset += blockdim(boffs2, nnz2, n)
    end
  end
  return boffsR, current_offset
end


text:
module ImGuiOpenGLBackend

using ModernGL
using LibCImGui

include("macros.jl")
include("context.jl")
include("font.jl")
include("device.jl")
include("platform.jl")
include("render.jl")
include("interface.jl")
include("image.jl")

end


text:
# Kendall got everything done before 6 PM (2022-02-01)
"""
    readWAVECAR(io::IO) -> ReciprocalWavefunction{3,Float64,Float32}

Reads a WAVECAR file output from a VASP 4.6 calcuation.

Information about VASP WAVECAR files and much of the code was pulled from the WaveTrans website
(originally written in FORTRAN): https://www.andrew.cmu.edu/user/feenstra/wavetrans/

This function is limited to WAVECAR files which have an RTAG value of 45200 (meaning the data is
given as a `Complex{Float64}`) and have only a collinear magnetic field applied, like WaveTrans.
"""
function readWAVECAR(io::IO)
    # Function to increment HKL values in place
    function incrementHKL!(hkl::AbstractVector{<:Integer}, bounds::AbstractVector{<:AbstractRange})
        # Loop through the vector indices, but in most cases we don't need them all
        for n in eachindex(hkl)
            # Increment the current vector component
            hkl[n] = (hkl[n] in bounds[n] ? hkl[n] + 1 : minimum(bounds[n]))
            # Only increment the next components if the current one is zero
            hkl[n] == 0 || break
        end
    end
    # Data entry counter (for the entries in the WAVECAR)
    count = 0
    # Number of bytes per record
    nrecl = Int(read(io, Float64))
    @debug "Record length: " * string(nrecl)
    # Number of spin components
    nspin = Int(read(io, Float64))
    # Check for the proper format
    rtag = Int(read(io, Float64))
    rtag == 45200 || error("Unsupported format: format value is " * string(rtag))
    # Jump to the next record
    count +=1; seek(io, count*nrecl)
    # Number of k-points
    nkpt = Int(read(io, Float64))
    # Number of bands
    nband = Int(read(io, Float64))
    # Energy cutoff
    ecut = read(io, Float64)
    # Real and reciprocal lattice vectors
    latt = BasisVectors{3}([read(io, Float64) for a = 1:3, b = 1:3])
    rlatt = dual(latt) * 2pi
    # Get HKL coefficient bounds (as done in WaveTrans)
    hklbounds = SVector{3,UnitRange{Int}}(-g:g for g in maxHKLindex(rlatt, ecut))
    # List of k-points
    klist = Vector{SVector{3,Float64}}(undef, nkpt)
    # Store band info in a vector (of size nkpt) containing vectors (of size nband)
    # containing tuples (containing band energy and occupancy)
    bands = [Vector{NTuple{2,Float64}}(undef, nband) for kp in 1:nkpt]
    # Plane wave coefficients
    # Vector (size nkpt) of vectors (size nband) of vectors (size npw) of ComplexF32
    # so planewave coeffs per band per k-point - that is a *lot* of data
    waves = [zeros(HKLData{3,Complex{Float32}}, hklbounds...) for kp in 1:nkpt, b in 1:nband]
    # Loop through the spins
    for s in 1:nspin
        # Loop through the k-points
        for kp in 1:nkpt
            # Seek to the next data
            count += 1; seek(io, count*nrecl)
            # Number of plane waves for this k-point
            pos = position(io)
            @debug string("File pointer at ", pos, " (", count, " * ", nrecl, ")")
            npw = Int(read(io, Float64))
            # Add the position of the k-point to the list
            klist[kp] = [read(io, Float64) for n = 1:3]
            # Get the bands associated with the k-point
            bands[kp] = [(read(io, Float64), (skip(io, 8); read(io, Float64))) for b in 1:nband]
            @info string(
                "Read in data for k-point ", kp, "/", nkpt, " (", npw, " planewaves/band)\n",
                "Reciprocal space coordinates: ", @sprintf("[%f %f %f]", klist[kp]...)
            )
            for b in 1:nband
                # Seek to the next entry
                count +=1; seek(io, count*nrecl)
                # Reset the HKL indices
                hkl = zeros(MVector{3,Int})
                for p in 1:npw
                    # Get the planewave component
                    pw = read(io, Complex{Float32})
                    # Increment the HKL indices
                    while true
                        # Get the energy of the vector
                        sumkg = [dot(klist[kp] + hkl, rlatt[:,n]) for n in 1:3]
                        etot = _selfdot(sumkg)/CVASP
                        # Break when the G-vector energy is below ecut
                        # This may occur immediately if the k-vector already meets the criteria
                        etot < ecut ? break : incrementHKL!(hkl, hklbounds)
                    end
                    # Store the data at the HKL index
                    # Note: data is stored first by k-points, then by bands
                    waves[kp, b][hkl...] = pw
                    # Increment it for the next iteration
                    incrementHKL!(hkl, hklbounds)
                end
            end
        end
    end
    # Now call the constructors
    return ReciprocalWavefunction{3,Float32}(rlatt, waves)
end

readWAVECAR(filename::AbstractString) = open(readWAVECAR, filename)
# Read a WAVECAR in the current directory
readWAVECAR() = readWAVECAR("WAVECAR")

"""
    readDOSCAR(io::IO) -> Tuple{DensityOfStates, Vector{ProjectedDensityOfStates}}

Reads a DOSCAR file from VASP and returns its data as a tuple containing the total and projected
density of states (if present).
"""
function readDOSCAR(io::IO)
    ln = readline(io)
    # Get number of ions and whether there is pdos data
    (nion, haspdos::Bool) = parse.(Int, split(ln))[[2,3]]
    for n in 2:6
        ln = readline(io)
    end
    # line6 = parse.(Float64, split(ln))
    # Get NEDOS and fermi energy
    nedos = parse(Int, split(ln)[3])
    fermi = parse(Float64, split(ln)[4])
    # Get the number of entries in the TDOS dataset
    ln = readline(io)
    entries = parse.(Float64, split(ln))
    tdos_raw = Matrix{Float64}(undef, length(entries), nedos)
    tdos_raw[:,1] = entries
    # Loop normally until end
    for n in 2:nedos
        ln = readline(io)
        tdos_raw[:,n] = parse.(Float64, split(ln))
    end
    # Create DensityOfStates struct for total DOS
    tdos = DensityOfStates(fermi, tdos_raw[1,:], tdos_raw[2,:], tdos_raw[3,:])
    pdos = Vector{ProjectedDensityOfStates}(undef, nion)
    if haspdos
        # Loop through all the ions
        for i in 1:nion
            readline(io)
            # Get the size of the matrix needed
            ln = readline(io)
            entries = parse.(Float64, split(ln))
            pdos_raw = Matrix{Float64}(undef, length(entries), nedos)
            pdos_raw[:,1] = entries
            # Add in the rest of the lines
            for n in 2:nedos
                ln = readline(io)
                pdos_raw[:,n] = parse.(Float64, split(ln))
            end
            # Add to the vector
            pdos[i] = ProjectedDensityOfStates(fermi, pdos_raw[1,:], pdos_raw[2:end,:])
        end
    else
        # Don't return a vector filled with a bunch of undefined objects...
        pdos = ProjectedDensityOfStates[]
    end
    return (tdos, pdos)
end

readDOSCAR(filename::AbstractString) = open(readDOSCAR, filename)

readDOSCAR() = open(readDOSCAR, "DOSCAR")

"""
    readPROCAR(io::IO) -> FatBands{3}

Reads an lm-decomposed PROCAR file from VASP and returns its data as a FatBands struct.
"""
function readPROCAR(io::IO)
    ln = readline(io)

    # Checks for lm-decomposed PROCAR
    if !(split(ln)[2] == "lm")
        error("Not a lm-decomposed PROCAR. Use a different PROCAR.")
    end
    has_phase = false
    if split(ln)[length(split(ln))] == "phase"
        has_phase = true
    end

    # Read header for number of kpoints, number of bands, number of ions
    (nkpt, nband, nion) = parse.(Int, split(readline(io))[[4, 8, 12]])

    # Initializes containers. Some are unused for now, with potentially useful information is commented out.
    # kptlist = KPointList{3}(Vector{SVector{3,Float64}}(undef, nkpt))
    # kptwt = Vector{Float64}(undef,nkpt)
    # occupancies = Matrix{Float64}(undef, nkpt, nband)
    energies = zeros(Float64, nkpt, nband)
    projband = zeros(Float64, 9, nion, nband, nkpt)
    phase_real = zeros(Float64, 9, nion, nband, nkpt)
    phase_imag = zeros(Float64, 9, nion, nband, nkpt)
    # Begins loop
    for i in 1:nkpt
        readuntil(io, "k-point")
        ln = readline(io)
        # kptlist[i] = parse.(Float64, split(ln)[3:5])
        # kptwt[i] = parse.(Float64, split(ln)[8])
        for j in 1:nband
            readuntil(io,"band")
            ln = readline(io)
            energies[i,j] = parse(Float64, split(ln)[4])
            #occupancies[i,j] = parse(Float64, split(ln)[7])
            readuntil(io,"tot\n")
            for k in 1:nion
                ln = readline(io)
                projband[:,k,j,i] = parse.(Float64, split(ln)[2:10])
            end
            readline(io)
            if has_phase
                readline(io)
                for k in 1:nion
                    ln = readline(io)
                    phase_real[:,k,j,i] = parse.(Float64, split(ln)[2:10])
                    ln = readline(io)
                    phase_imag[:,k,j,i] = parse.(Float64, split(ln)[2:10])
                end
            end
        end
    end

    # Returns a FatBands struct
    return FatBands{3}(
        energies,
        projband,
        phase_real,
        phase_imag,
    )
end

readPROCAR(filename::AbstractString) = open(readPROCAR, filename)

readPROCAR() = open(readPROCAR, "PROCAR")


text:
# This file is part of IntegerSequences.
# Copyright Peter Luschny. License is MIT.

(@__DIR__) ∉ LOAD_PATH && push!(LOAD_PATH, (@__DIR__))

module EulerTransforms

using Nemo, NumberTheory

export ModuleEulerTransforms
export EulerTransform
export V006171, L006171, V107895, L107895, V061256, L061256
export V190905, L190905, V275585, L275585, V290351, L290351
export V052847, L052847

"""
* V006171, L006171, V107895, L107895, V061256, L061256, V190905, L190905, V275585, L275585, V290351, L290351
"""
const ModuleEulerTransforms = ""

"""
Return the Euler transform of f.
"""
function EulerTransform(f::Function)
    function E(n)
        haskey(CacheET, (f, n)) && return CacheET[(f, n)]
        n == 0 && return ZZ(1)
        a(j) = ZZ(sum(d*f(Int(d)) for d in Divisors(j)))
        b = sum(a(j) * E(n-j) for j in 1:n)
        r = div(b, n)
        CacheET[(f, n)] = r
        return r
    end
end

const CacheET = Dict{Tuple{Function, Int}, fmpz}()

"""
Return the number of factorization patterns of polynomials of degree n over integers.
"""
V006171(n) = EulerTransform(τ)(n)

"""
Return a list of length len of the Euler transform of tau.
"""
L006171(len) = [V006171(n) for n in 0:len-1]

"""
Return the Euler transform of sigma.
"""
V061256(n) = EulerTransform(σ)(n)

"""
Return a list of length len of the Euler transform of sigma.
"""
L061256(len) = [V061256(n) for n in 0:len-1]

"""
Return the Euler transform of sigma_2.
"""
V275585(n) = EulerTransform(σ2)(n)

"""
Return a list of length len of the Euler transform of sigma_2.
"""
L275585(len) = [V275585(n) for n in 0:len-1]

"""
Return the Euler transform of the factorial.
"""
V107895(n) = EulerTransform(fac)(n)

"""
Return a list of length len of the Euler transform of the factorial.
"""
L107895(len) = [V107895(n) for n in 0:len-1]

"""
Return the Euler transform of the swinging factorial (A056040).
"""
V190905(n) = EulerTransform(n -> div(fac(n), fac(div(n, 2))^2))(n)

"""
Return a list of length len of the Euler transform of the swinging factorial.
"""
L190905(len) = [V190905(n) for n in 0:len-1]

"""
Return the Euler transform of the Bell numbers.
"""
V290351(n) = EulerTransform(Nemo.bell)(n)

"""
Return a list of length len of the Euler transform of the Bell numbers.
"""
L290351(len) = [V290351(n) for n in 0:len-1]

"""
Return the Euler transform of [0, 1, 2, 3, ...].
"""
V052847(n) = EulerTransform(k -> k-1)(n)

"""
Return a list of length len of the Euler transform of [0, 1, 2, 3, ...].
"""
L052847(len) = [V052847(n) for n in 0:len-1]

#START-TEST-########################################################

using Test, SeqTests

function test()

    @testset "EulerTransform" begin

        if is_oeis_installed()

            L = [L107895, L190905, L061256, L275585, L006171, L290351, L052847]
            SeqTest(L, 'L')

            V = [V107895, V190905, V061256, V275585, V006171]
            for v in V SeqTest(v, 'V', 0) end
        end
    end
end

function demo()
end

"""
for n in 0:150 V107895(n) end ::
    0.303605 seconds (1.01 M allocations: 32.276 MiB, 29.90% gc time)
"""
function perf()
    @time (for n in 0:150 V107895(n) end)
end

function main()
    test()
    demo()
    perf()
end

main()

end # module


text:
function gradient(f, xs...)
  xs = param.(xs)
  back!(f(xs...))
  grad.(xs)
end

function ngradient(f, xs::AbstractArray...)
  grads = zeros.(xs)
  for (x, Δ) in zip(xs, grads), i in 1:length(x)
    δ = sqrt(eps())
    tmp = x[i]
    x[i] = tmp - δ/2
    y1 = f(xs...)
    x[i] = tmp + δ/2
    y2 = f(xs...)
    x[i] = tmp
    Δ[i] = (y2-y1)/δ
  end
  return grads
end

gradcheck(f, xs...) =
  all(isapprox.(ngradient(f, xs...),
                gradient(f, xs...), rtol = 1e-5, atol = 1e-5))


text:
include("HydroEnergyCascade.jl")
include("HydroDispatchCascade.jl")

export get_upstream
export set_upstream!


text:
# Generate maps of species richness, predicted numbers of links (median and 68% PI) and stability score

## Map of species richness

# Function to read species richness from a dataset in .tif extension
# Usable with data from biodiversitymapping.org
function read_richness(data)
    # Read dataset in .tif extension using ArchGDAL
    raw_tif = ArchGDAL.read(ArchGDAL.read(data))

    # Convert in matrix format and transpose
    raw = Matrix(transpose(Float32.(raw_tif[:, :, 1])))
    raw = raw[end:-1:1, :]

    # We are working on the actual pixels, which results in a big map, but we will
    # coarsen it later
    sp = SimpleSDMPredictor(raw, -16671202.0, 16818798.0, -8418528.0, 8381472.0)

    # We want to remove the 255 (0xff) and 65535 (oxffff)(NA values)
    sp_255 = filter(x -> sp.grid[x] in [255, 65535], CartesianIndices(sp.grid))

    for i in LinearIndices(sp.grid)[sp_255]
        sp.grid[i] = 0
    end

    # Return raster
    return(sp)
end

# Data paths
# Species richness of terrestrial mammals, birds and amphibians
# Rasters obtained from QGIS after uniformizing their extents
mammals_path = "data/BiodiversityMapping/mammals.tif"
birds_path = "data/BiodiversityMapping/birds.tif"
amphibians_path = "data/BiodiversityMapping/amphibians.tif"

# Read species richness of specific taxons (terrestrial mammals, birds and amphibians)
sp_mammals = read_richness(mammals_path)
sp_birds = read_richness(birds_path)
sp_amphibians = read_richness(amphibians_path)

# Raster sizes
size(sp_mammals)
size(sp_birds)
size(sp_amphibians)

# Visualize species richness
heatmap(sp_mammals)
heatmap(sp_birds)
heatmap(sp_amphibians)

# New SimpleSDMResponse with same coordinates
total_sp_richness = similar(sp_birds)

# Sum of birds, mammals and amphibians
for i in LinearIndices(total_sp_richness.grid)
    total_sp_richness.grid[i] = sp_mammals.grid[i] + sp_birds.grid[i] + sp_amphibians.grid[i]
end

# Change 0 to NAs
sp_NA = filter(x -> total_sp_richness.grid[x] == 0, CartesianIndices(total_sp_richness.grid))

for i in LinearIndices(total_sp_richness.grid)[sp_NA]
    total_sp_richness.grid[i] = NaN
end

# Map
# TK: bounding box at the beginning of this script
xlim = (-10500000, -3800000)
ylim = (5000000, 8381472)

heatmap(total_sp_richness, dpi=300, clim=(0,400), xlim=xlim, ylim=ylim, frame=:box)
savefig("figures/species_richness.png")

## Maps of predicted numbers of links (median and 68% PI) and stability score

# Load measures
measures = CSV.read(joinpath("data", "measures.csv"), DataFrame)

# Where are there at least 5 species?
cidx = filter(x -> total_sp_richness.grid[x] >= 5, CartesianIndices(total_sp_richness.grid))
cidx0 = filter(x -> total_sp_richness.grid[x] < 5, CartesianIndices(total_sp_richness.grid))

# SimpleSDMResponse objects
L_median_SDM = similar(total_sp_richness)
L_PI_SDM = similar(total_sp_richness)
stab_SDM = similar(total_sp_richness)

# Change the values in the SimpleSDMResponse objects by their respective measures according to species richness
function map_measures(measure, SDMObject)
    for i in LinearIndices(total_sp_richness.grid)[cidx]
    SDMObject.grid[i] = measures[measures.S .== total_sp_richness.grid[i], measure][1]
    end
    for i in LinearIndices(total_sp_richness.grid)[cidx0]
        SDMObject.grid[i] = NaN
    end
end

measures[!, :L_median_log] .= log10.(measures.L_median)
measures[!, :L_PI_log] .= log10.(measures.L_PI)

map_measures("L_median_log", L_median_SDM)
map_measures("L_PI_log", L_PI_SDM)
map_measures("stab", stab_SDM)

# Maps

heatmap(L_median_SDM, color=:BrBG_11, clim=(0,4), xlim=xlim, ylim=ylim, dpi=300, framestyle=:box) # log scale
savefig(joinpath("figures", "maps_L_median.png"))

heatmap(L_PI_SDM, color=:BrBG_11, clim=(0,4), xlim=xlim, ylim=ylim, dpi=300, framestyle=:box) # log scale
savefig(joinpath("figures", "maps_L_PI.png"))

heatmap(stab_SDM, color=:RdBu_11, dpi=300,  xlim=xlim, ylim=ylim, framestyle=:box)
savefig(joinpath("figures", "maps_stab.png"))

max_stars_repo_name: PoisotLab/worldwide-links-prediction

max_issues_repo_name: PoisotLab/worldwide-links-prediction

max_forks_repo_name: PoisotLab/worldwide-links-prediction


text:
MuZero.action_space_size(::Type{TicTacToe}) = 9
MuZero.is_terminal(game::TicTacToe) = game |> MuZero.state |> is_terminal
MuZero.get_state(g::TicTacToe) = MuZero.state(g)
state(g::TicTacToe) = g.states[end]

function make_tictactoe_config()
    make_board_game_config(
        TicTacToe,
        network_type=Incremental,
        max_moves=9,
        root_dirichlet_alpha=0.3,
        action_space_size=9, # Number of all possible actions in the game
        support_size = 10, # Value and reward are scaled (with almost sqrt) and encoded on a vector with a range of -support_size to support_size. Choose it so that support_size <= sqrt(max(abs(discounted reward)))
        observation_shape = (3, 3, 3),
        stacked_observations = 0, # Number of previous observations and previous actions to add to the current observation
        hidden_state_size = 32, # The 1-dimensional size of the encoded game state
        checkpoint_interval = 100, #  Number of training steps before using the model for self-playing
        num_actors = 3,
        weight_decay = 1e-3,
        lr_init=0.002f0,
        save_interval = 10000,
        num_unroll_steps=1,
    )
end

MuZero.num_cols(::Type{TicTacToeState}) = 3
MuZero.num_rows(::Type{TicTacToeState}) = 3
MuZero.getstatetype(::Type{TicTacToe}) = TicTacToeState


text:
module ODBC

using Tables, CategoricalArrays, WeakRefStrings, DataFrames, Dates, DecFP

export DataFrame, odbcdf

include("API.jl")

"just a block of memory; T is the element type, `len` is total # of **bytes** pointed to, and `elsize` is size of each element"
mutable struct Block{T}
    ptr::Ptr{T}    # pointer to a block of memory
    len::Int       # total # of bytes in block
    elsize::Int    # size between elements in bytes
end

"""
Block allocator:
    -Takes an element type, and number of elements to allocate in a linear block
    -Optionally specify an extra dimension of elements that make up each element (i.e. container types)
"""
function Block(::Type{T}, elements::Int, extradim::Integer=1) where {T}
    len = sizeof(T) * elements * extradim
    block = Block{T}(convert(Ptr{T}, Libc.malloc(len)), len, sizeof(T) * extradim)
    finalizer(x->Libc.free(x.ptr), block)
    return block
end

# used for getting messages back from ODBC driver manager; SQLDrivers, SQLError, etc.
Base.string(block::Block, len::Integer) = String(transcode(UInt8, unsafe_wrap(Array, block.ptr, len, own=false)))

struct ODBCError <: Exception
    msg::String
end

const BUFLEN = 1024

function ODBCError(handle::Ptr{Cvoid}, handletype::Int16)
    i = Int16(1)
    state = ODBC.Block(ODBC.API.SQLWCHAR, 6)
    native = Ref{ODBC.API.SQLINTEGER}()
    error_msg = ODBC.Block(ODBC.API.SQLWCHAR, BUFLEN)
    msg_length = Ref{ODBC.API.SQLSMALLINT}()
    while ODBC.API.SQLGetDiagRec(handletype, handle, i, state.ptr, native, error_msg.ptr, BUFLEN, msg_length) == ODBC.API.SQL_SUCCESS
        st  = string(state, 5)
        msg = string(error_msg, msg_length[])
        println("[ODBC] $st: $msg")
        i += 1
    end
    return true
end

#Macros to to check if a function returned a success value or not
macro CHECK(handle, handletype, func)
    str = string(func)
    esc(quote
        ret = $func
        ret != ODBC.API.SQL_SUCCESS && ret != ODBC.API.SQL_SUCCESS_WITH_INFO && ODBCError($handle, $handletype) &&
            throw(ODBCError("$($str) failed; return code: $ret => $(ODBC.API.RETURN_VALUES[ret])"))
        nothing
    end)
end

Base.@deprecate listdrivers ODBC.drivers
Base.@deprecate listdsns ODBC.dsns

"List ODBC drivers that have been installed and registered"
function drivers()
    descriptions = String[]
    attributes   = String[]
    driver_desc = Block(ODBC.API.SQLWCHAR, BUFLEN)
    desc_length = Ref{ODBC.API.SQLSMALLINT}()
    driver_attr = Block(ODBC.API.SQLWCHAR, BUFLEN)
    attr_length = Ref{ODBC.API.SQLSMALLINT}()
    dir = ODBC.API.SQL_FETCH_FIRST
    while ODBC.API.SQLDrivers(ENV[], dir, driver_desc.ptr, BUFLEN, desc_length, driver_attr.ptr, BUFLEN, attr_length) == ODBC.API.SQL_SUCCESS
        push!(descriptions, string(driver_desc, desc_length[]))
        push!(attributes,   string(driver_attr, attr_length[]))
        dir = ODBC.API.SQL_FETCH_NEXT
    end
    return [descriptions attributes]
end

"List ODBC DSNs, both user and system, that have been previously defined"
function dsns()
    descriptions = String[]
    attributes   = String[]
    dsn_desc    = Block(ODBC.API.SQLWCHAR, BUFLEN)
    desc_length = Ref{ODBC.API.SQLSMALLINT}()
    dsn_attr    = Block(ODBC.API.SQLWCHAR, BUFLEN)
    attr_length = Ref{ODBC.API.SQLSMALLINT}()
    dir = ODBC.API.SQL_FETCH_FIRST
    while ODBC.API.SQLDataSources(ENV[], dir, dsn_desc.ptr, BUFLEN, desc_length, dsn_attr.ptr, BUFLEN, attr_length) == ODBC.API.SQL_SUCCESS
        push!(descriptions, string(dsn_desc, desc_length[]))
        push!(attributes,   string(dsn_attr, attr_length[]))
        dir = ODBC.API.SQL_FETCH_NEXT
    end
    return [descriptions attributes]
end

"""
A DSN represents an established ODBC connection.
It is passed to most other ODBC methods as a first argument
"""
mutable struct DSN
    dsn::String
    dbc_ptr::Ptr{Cvoid}
    stmt_ptr::Ptr{Cvoid}
    stmt_ptr2::Ptr{Cvoid}
end

Base.show(io::IO,conn::DSN) = print(io, "ODBC.DSN($(conn.dsn))")

const dsn = DSN("", C_NULL, C_NULL, C_NULL)

"""
Construct a `DSN` type by connecting to a valid ODBC DSN or by specifying a valid connection string.
Takes optional 2nd and 3rd arguments for `username` and `password`, respectively.
1st argument `dsn` can be either the name of a pre-defined ODBC DSN or a valid connection string.
A great resource for building valid connection strings is [http://www.connectionstrings.com/](http://www.connectionstrings.com/).
"""
function DSN(connectionstring::AbstractString, username::AbstractString=String(""), password::AbstractString=String(""); prompt::Bool=true)
    dbc = ODBC.ODBCAllocHandle(ODBC.API.SQL_HANDLE_DBC, ODBC.ENV[])
    dsns = ODBC.dsns()
    found = false
    for d in dsns[:,1]
        connectionstring == d && (found = true)
    end
    if found
        @CHECK dbc ODBC.API.SQL_HANDLE_DBC ODBC.API.SQLConnect(dbc, connectionstring, username, password)
    else
        connectionstring = ODBCDriverConnect!(dbc, connectionstring, prompt)
    end
    stmt = ODBCAllocHandle(ODBC.API.SQL_HANDLE_STMT, dbc)
    stmt2 = ODBCAllocHandle(ODBC.API.SQL_HANDLE_STMT, dbc)
    global dsn
    dsn.dsn = connectionstring
    dsn.dbc_ptr = dbc
    dsn.stmt_ptr = stmt
    dsn.stmt_ptr2 = stmt2
    return DSN(connectionstring, dbc, stmt, stmt2)
end

"disconnect a connected `DSN`"
function disconnect!(conn::DSN)
    ODBCFreeStmt!(conn.stmt_ptr)
    ODBCFreeStmt!(conn.stmt_ptr2)
    ODBC.API.SQLDisconnect(conn.dbc_ptr)
    return nothing
end

mutable struct Statement
    dsn::DSN
    stmt::Ptr{Cvoid}
    query::String
    task::Task
end

# used to 'clear' a statement of bound columns, resultsets,
# and other bound parameters in preparation for a subsequent query
function ODBCFreeStmt!(stmt)
    ODBC.API.SQLFreeStmt(stmt, ODBC.API.SQL_CLOSE)
    ODBC.API.SQLFreeStmt(stmt, ODBC.API.SQL_UNBIND)
    ODBC.API.SQLFreeStmt(stmt, ODBC.API.SQL_RESET_PARAMS)
end

# "Allocate ODBC handles for interacting with the ODBC Driver Manager"
function ODBCAllocHandle(handletype, parenthandle)
    handle = Ref{Ptr{Cvoid}}()
    API.SQLAllocHandle(handletype, parenthandle, handle)
    handle = handle[]
    if handletype == API.SQL_HANDLE_ENV
        API.SQLSetEnvAttr(handle, API.SQL_ATTR_ODBC_VERSION, API.SQL_OV_ODBC3)
    end
    return handle
end

# "Alternative connect function that allows user to create datasources on the fly through opening the ODBC admin"
function ODBCDriverConnect!(dbc::Ptr{Cvoid}, conn_string, prompt::Bool)
    @static if Sys.iswindows()
        driver_prompt = prompt ? API.SQL_DRIVER_PROMPT : API.SQL_DRIVER_NOPROMPT
        window_handle = prompt ? ccall((:GetForegroundWindow, :user32), Ptr{Cvoid}, () ) : C_NULL
    else
        driver_prompt = API.SQL_DRIVER_NOPROMPT
        window_handle = C_NULL
    end
    out_conn = Block(API.SQLWCHAR, BUFLEN)
    out_buff = Ref{Int16}()
    @CHECK dbc API.SQL_HANDLE_DBC API.SQLDriverConnect(dbc, window_handle, conn_string, out_conn.ptr, BUFLEN, out_buff, driver_prompt)
    connection_string = string(out_conn, out_buff[])
    return connection_string
end

"`prepare` prepares an SQL statement to be executed"
function prepare(dsn::DSN, query::AbstractString)
    stmt = ODBCAllocHandle(API.SQL_HANDLE_STMT, dsn.dbc_ptr)
    @CHECK stmt API.SQL_HANDLE_STMT API.SQLPrepare(stmt, query)
    return Statement(dsn, stmt, query, Task(1))
end

function execute!(statement::Statement, values)
    stmt = statement.stmt
    values2 = Any[cast(x) for x in values]
    pointers = Ptr[]
    types = map(typeof, values2)
    for (i, v) in enumerate(values2)
        if ismissing(v)
            @CHECK stmt API.SQL_HANDLE_STMT API.SQLBindParameter(stmt, i, API.SQL_PARAM_INPUT,
                API.SQL_C_CHAR, API.SQL_CHAR, 0, 0, C_NULL, 0, Ref(API.SQL_NULL_DATA))
        else
            ctype, sqltype = API.julia2C[types[i]], API.julia2SQL[types[i]]
            csize, len, dgts = sqllength(v), clength(v), digits(v)
            ptr = getpointer(types[i], values2, i)
            # println("ctype: $ctype, sqltype: $sqltype, digits: $dgts, len: $len, csize: $csize")
            push!(pointers, ptr)
            @CHECK stmt API.SQL_HANDLE_STMT API.SQLBindParameter(stmt, i, API.SQL_PARAM_INPUT,
                ctype, sqltype, csize, dgts, ptr, len, Ref(len))
        end
    end
    execute!(statement)
    return
end

function execute!(statement::Statement)
    stmt = statement.stmt
    @CHECK stmt API.SQL_HANDLE_STMT API.SQLExecute(stmt)
    return
end

"`execute!` is a minimal method for just executing an SQL `query` string. No results are checked for or returned."
function execute!(dsn::DSN, query::AbstractString, stmt=dsn.stmt_ptr)
    ODBCFreeStmt!(stmt)
    @CHECK stmt API.SQL_HANDLE_STMT API.SQLExecDirect(stmt, query)
    return
end

include("Sink.jl")
include("Query.jl")
# include("sqlreplmode.jl")

const ENV = Ref{Ptr{Cvoid}}()

function __init__()
    ENV[] = ODBC.ODBCAllocHandle(ODBC.API.SQL_HANDLE_ENV, ODBC.API.SQL_NULL_HANDLE)
end

end #ODBC module


text:
using Statistics, CUDA

"""
`bnmoments(;momentum=0.1, mean=nothing, var=nothing, meaninit=zeros, varinit=ones)` can be used
 directly load moments from data. `meaninit` and `varinit` are called if `mean` and `var`
are nothing. Type and size of the `mean` and `var` are determined automatically from the inputs
in the `batchnorm` calls. A `BNMoments` object is returned.

# BNMoments
A high-level data structure used to store running mean and running variance
of batch normalization with the following fields:

 `momentum::AbstractFloat`: A real number between 0 and 1 to be used as the scale of
  last mean and variance. The existing running mean or variance is multiplied by
  (1-momentum).

 `mean`: The running mean.

 `var`: The running variance.

 `meaninit`: The function used for initialize the running mean. Should either be `nothing` or
of the form `(eltype, dims...)->data`. `zeros` is a good option.

 `varinit`: The function used for initialize the running variance. Should either be `nothing` or
`(eltype, dims...)->data`. `ones` is a good option.

"""
bnmoments(;momentum=0.1, meaninit=zeros, varinit=ones, mean=nothing, var=nothing) =
    BNMoments(momentum, mean, var, meaninit, varinit)

mutable struct BNMoments
    momentum::AbstractFloat
    mean
    var
    meaninit
    varinit
end

#TODO: improve the documentation. Explain what batchnorm is.
"""

`batchnorm(x[, moments, params]; kwargs...)` performs batch normalization to `x`
with optional scaling factor and bias stored in `params`.

2d, 4d and 5d inputs are supported. Mean and variance are computed over
dimensions (2,), (1,2,4) and (1,2,3,5) for 2d, 4d and 5d arrays, respectively.

`moments` stores running mean and variance to be used in testing.
It is optional in the training mode, but mandatory in the test mode.
Training and test modes are controlled by the `training` keyword argument.

`params` stores the optional affine parameters gamma and beta.
`bnparams` function can be used to initialize `params`.

# Example

    # Inilization, C is an integer
    moments = bnmoments()
    params = bnparams(C)
    ...
    # size(x) -> (H, W, C, N)
    y = batchnorm(x, moments, params)
    # size(y) -> (H, W, C, N)

# Keywords

 `eps=1e-5`: The epsilon parameter added to the variance to avoid division by 0.

 `training`: When `training` is true, the mean and variance of `x` are used and `moments`
 argument is modified if it is provided. When `training` is false, mean and variance stored in
 the `moments` argument are used. Default value is `true` when at least one of `x` and `params`
 is `AutoGrad.Value`, `false` otherwise.

"""
function batchnorm(x, moments::Union{BNMoments, Nothing}=nothing, params=nothing;
                   training=training(), o...)
    xnd = ndims(x)
    a = (x,)
    if params !== nothing
        g = reshape(_bnscale(params), _wsize(x))
        b = reshape(_bnbias(params), _wsize(x))
        a = (g, b, x)
    end
    # if ~isa(training, Bool)
    #     training = isa(x, AutoGrad.Value) || isa(params, AutoGrad.Value)
    # end
    if xnd == 2
        return batchnorm2(a...; o...,
                          moments=moments,
                          training=training,
                          cache=BNCache())
    elseif xnd in [4, 5]
        return batchnorm4(a...; o...,
                          moments=moments,
                          training=training,
                          cache=BNCache())
    else
        error("Unsupported input dimension ", xnd)
    end
end

"""
`bnparams(etype, channels)` creates a single 1d array that contains both
scale and bias of batchnorm, where the first half is scale and the
second half is bias.

`bnparams(channels)` calls `bnparams` with `etype=Float64`, following Julia convention

"""
function bnparams(etype, channels::Integer)
    buf = Array{etype}(undef,2channels)
    buf[1:channels] .= 1
    buf[channels+1:end] .= 0
    return buf
end

bnparams(channels::Integer) = bnparams(Float64, channels)

#=
LOW-LEVEL API that won't be exported by default
=#
@inline _bnscale(param) = param[1:div(length(param), 2)]
@inline _bnbias(param) = param[div(length(param), 2)+1:end]

const BN_MODE_SPATIAL = 1
const BN_MODE_ACTIVATION = 0
const CUDNN_BN_MIN_EPS = 1e-5

# A black box data type for storing the bn state
mutable struct BNCache
    mean
    ivar
    dx
    dg
    db
end

BNCache() = BNCache(nothing, nothing, nothing, nothing, nothing)

# Dimension helpers
@inline _wsize(y) = ((1 for _=1:ndims(y)-2)..., size(y)[end-1], 1)
@inline _reddims(y) = ((i for i=1:ndims(y)-2)..., ndims(y))

# TODO: consider automatic type conversion
function _lazy_init!(m::BNMoments, x)
    x = value(x)
    buf_size = (ndims(x) > 2) ? _wsize(x) : (size(x,1), 1)
    tx = typeof(x)
    ex = eltype(x)
    m.momentum = ex(m.momentum)
    if m.mean == nothing
        m.mean = tx(m.meaninit(ex, buf_size...))
    end
    if m.var == nothing
        m.var = tx(m.varinit(ex, buf_size...))
    end
end

# Only spatial mode is supported
# TODO: support per-activation mode
function batchnorm4(g::KnetArray{T}, b::KnetArray{T}, x::KnetArray{T};
                    training=training(),
                    cache=nothing,
                    moments=nothing,
                    eps=1e-5,
                    alpha=1, beta=0,
                    handle = CUDNN.handle(),
                    cache_verbose=false, #reporting cache uses
                    o...) where {T}
    y = KnetArray{T}(undef,size(x))
    weight_size = _wsize(y)
    # TODO: implement other bn mode
    bnmode = BN_MODE_SPATIAL
    # epsilon fix
    if eps < CUDNN_BN_MIN_EPS
        eps = CUDNN_BN_MIN_EPS
        warn("eps ", eps,
             " is too small for cudnn, so it is set to ",
             CUDNN_BN_MIN_EPS)
    end
    # moments
    if moments !== nothing
        _lazy_init!(moments, x)
        running_mean = moments.mean
        running_var = moments.var
        momentum = moments.momentum
    else
        running_mean, running_var = CU_NULL, CU_NULL
        momentum = .1
    end
    # The training mode
    if training
        # Cache the mean and ivar for later
        if cache !== nothing
            mean = KnetArray{T}(undef,weight_size)
            ivar = KnetArray{T}(undef,weight_size)
        else
            mean = CU_NULL
            ivar = CU_NULL
        end
        # @cudnn(cudnnBatchNormalizationForwardTraining,
        #       # Types
        #       (Cptr, UInt32,
        #        Ptr{T}, Ptr{T}, #alpha and beta
        #        Cptr, Ptr{T}, #xdesc and x
        #        Cptr, Ptr{T}, #ydesc and y
        #        Cptr, Ptr{T}, Ptr{T}, #desc, weight and bias
        #        Cdouble, Ptr{T}, Ptr{T}, #Decay factor, Running mean and Running var
        #        Cdouble, # eps
        #        Ptr{T}, Ptr{T}), #Cached mean and ivar
        #       # Actual Arguments
        #       handle, bnmode,
        #       Ref(T(alpha)), Ref(T(beta)),
        #       TD(x), x, #x
        #       TD(y), y, #y
        #       TD(g), g, b, #params
        #       momentum, running_mean, running_var,
        #       eps, mean, ivar)
        bnmode = CUDNN.cudnnBatchNormMode_t(bnmode)
        CUDNN.cudnnBatchNormalizationForwardTraining(handle, bnmode, Ref(T(alpha)), Ref(T(beta)), TD(x), x, TD(y), y, TD(g), g, b, momentum, running_mean, running_var, eps, mean, ivar)

        # Cache the resulting mean and inverse variance
        if cache != nothing
            cache_verbose && info("mean and ivar data saved to cache")
            cache.mean = mean
            cache.ivar = ivar
        end
    else
        @assert (moments!==nothing) "You must provide moments for the test mode!"
        # @cudnn(cudnnBatchNormalizationForwardInference,
        #       # Types
        #       (Cptr, UInt32,
        #        Ptr{T}, Ptr{T},
        #        Cptr, Ptr{T}, #x
        #        Cptr, Ptr{T}, #y
        #        Cptr, Ptr{T}, Ptr{T}, #params
        #        Ptr{T}, Ptr{T}, #rm and rf
        #        Cdouble),
        #       handle, bnmode,
        #       Ref(T(alpha)), Ref(T(beta)), #alpha and bete
        #       TD(x), x, # xdesc and x
        #       TD(y), y, #ydesc and y
        #       TD(g), g, b, #desc, scale and bias
        #       running_mean, running_var, #estimated stuff
        #       eps) #epsilon
        bnmode = CUDNN.cudnnBatchNormMode_t(bnmode)
        CUDNN.cudnnBatchNormalizationForwardInference(handle, bnmode, Ref(T(alpha)), Ref(T(beta)), TD(x), x, TD(y), y, TD(g), g, b, running_mean, running_var, eps)
    end
    return y
end

function batchnorm4(x::KnetArray{T};o...) where {T}
    # Dummy buffers
    #  (cudnn doesn't support bn w/o affine
    #  although it is used in many applications)
    g = KnetArray{T}(ones(_wsize(x)...))
    b = KnetArray{T}(zeros(_wsize(x)...))
    return batchnorm4(g, b, x; o...)
end

function batchnorm4_back(g::Union{KnetArray{T}, Nothing},
                         x::KnetArray{T}, dy::KnetArray{T};
                         training=training(),
                         cache=nothing,
                         moments=nothing,
                         grad_cache_disabled=false,
                         eps=1e-5, alpha=1, beta=0,
                         dalpha=1, dbeta=0,
                         handle = CUDNN.handle(),
                         cache_verbose=false,
                         o...) where {T}
    if training
        dx = KnetArray{T}(undef,size(x))
        weight_size = _wsize(dy)
        if g==nothing; g=KnetArray{T}(ones(weight_size)); end
        dg = KnetArray{T}(undef,weight_size)
        db = KnetArray{T}(undef,weight_size)
        # TODO: support other modes
        bnmode = BN_MODE_SPATIAL
        if cache !== nothing # (Assume cache still exists)
            mean, ivar = cache.mean, cache.ivar
            cache_verbose && info("mean and ivar are fetched from the cache")
        else
            mean, ivar = CU_NULL, CU_NULL
        end
        # @cudnn(cudnnBatchNormalizationBackward,
        #       # C Types
        #       (Cptr, UInt32,
        #        Ptr{T}, Ptr{T}, #data difs
        #        Ptr{T}, Ptr{T}, #param difs
        #        Cptr, Ptr{T}, #x
        #        Cptr, Ptr{T}, #dy
        #        Cptr, Ptr{T}, #dx
        #        Cptr, Ptr{T}, Ptr{T}, Ptr{T}, #desc,g,dg,db
        #        Cdouble, Ptr{T}, Ptr{T}),
        #       # Actual arguments
        #       handle, bnmode,
        #       Ref(T(alpha)), Ref(T(beta)),
        #       Ref(T(dalpha)), Ref(T(dbeta)),
        #       TD(x), x, TD(dy), dy, TD(dx), dx,
        #       TD(g), g, dg, db,
        #       eps, mean, ivar)
        bnmode = CUDNN.cudnnBatchNormMode_t(bnmode)
        CUDNN.cudnnBatchNormalizationBackward(handle, bnmode, Ref(T(alpha)), Ref(T(beta)), Ref(T(dalpha)), Ref(T(dbeta)), TD(x), x, TD(dy), dy, TD(dx), dx, TD(g), g, dg, db, eps, mean, ivar)

    else
        # At test mode, g .*( x ./ sqrt(var) - mean ./ sqrt(var)) .+ beta
        # is performed;
        # so the derivative dx = dy .* g. / sqrt(var + eps) since mean and var
        # are constants
        # Note: moments must exist since otherwise forward pass fails
        ivar = 1 ./ sqrt.(moments.var .+ eps)
        dx = (g !== nothing) ? (dy .* g .* ivar) : (dy .* ivar)
        if g !== nothing
            dg = sum(dy .* (x .- moments.mean) .* ivar, dims=_reddims(dy))
            db = sum(dy, dims=_reddims(dy))
        else
            dg, db = nothing, nothing
        end
    end
    return dg, db, dx
end

function batchnorm4g(g::Union{KnetArray{T}, Array{T}},
                     x::Union{KnetArray{T}, Array{T}},
                     dy::Union{KnetArray{T}, Array{T}};
                     cache=nothing, o...) where {T}
    dg, db, dx = batchnorm4_back(g, x, dy; cache=cache, o...)
    if cache !== nothing
        cache.dx, cache.db = dx, db
    else
        warn("Calling batchnorm with affine without",
             " cache is not recommended for performance")
    end
    return dg
end

function batchnorm4b(dy::Union{KnetArray{T}, Array{T}};
                     cache=nothing, o...) where {T}
    (cache == nothing || cache.db == nothing) && return sum(dy, dims=_reddims(dy))
    return cache.db
end

function batchnorm4x(g::Union{KnetArray{T}, Array{T}},
                     x::Union{KnetArray{T}, Array{T}},
                     dy::Union{KnetArray{T}, Array{T}};
                     cache=nothing, o...) where {T}
    if cache !== nothing && cache.dx !== nothing
        return cache.dx
    end
    return batchnorm4_back(g, x, dy; cache=cache, o...)[3]
end

function batchnorm4x(x::Union{KnetArray{T}, Array{T}},
                     dy::Union{KnetArray{T}, Array{T}}
                     ;o...) where {T}
    return batchnorm4_back(nothing, x, dy; o...)[3]
end

# CPU Implementation
function batchnorm4(g::Array{T}, b::Array{T}, x::Array{T};
                    o...) where {T}
    return _batchnorm4_fused(g, b, x; o...)
end

function batchnorm4(x::Array{T};
                    o...) where {T}
    return _batchnorm4_fused(nothing,nothing,x; o...)
end

function _batchnorm4_fused(g, b, x::Array{T};
                           eps=1e-5, training=training(),
                           cache=nothing, moments=nothing,
                           o...) where {T}
    y = copy(x)
    eps = T(eps)
    dims = _reddims(y)
    if moments!==nothing
        _lazy_init!(moments, x)
    end
    affine = (g !== nothing)
    if training
        mu = mean(y, dims=dims)
        sigma2 = var(y, dims=dims, corrected=false, mean=mu)
        _update_moments!(moments, mu, sigma2)
        sigma2 .= 1 ./ sqrt.(sigma2 .+ eps)
        ivar = sigma2
        if affine
            y .= g .* (y .- mu) .* ivar .+ b
        else
            y .= (y .- mu) .* ivar
        end
        if cache !== nothing
            cache.ivar = ivar
            cache.mean = mu
        end
    else
        @assert (moments!==nothing)  "Moments must be provided in test mode"
        ivar = 1 ./ sqrt.(moments.var .+ eps)
        if affine
            y .= g .* (y .- moments.mean) .* ivar .+ b
        else
            y .= (y .- moments.mean) .* ivar
        end
    end
    return y
end

function _update_moments!(moments, mu, sigma2)
    moments === nothing && return
    moments.mean .= moments.momentum .* mu .+
        (1 .- moments.momentum) .* moments.mean
    moments.var .= moments.momentum .* sigma2 .+
        (1 .- moments.momentum) .* moments.var
end

# CPU backward
function batchnorm4_back(g::Union{Array{T}, Nothing}, x::Array{T}, dy::Array{T};
                         eps=1e-5, training=training(),
                         cache=nothing, moments=nothing,  o...) where {T}
    eps = T(eps)
    dims = _reddims(x)
    if training
        mu, ivar = _get_cache_data(cache, x, eps)
        x_mu = x .- mu
        # equations from the original paper
        dyivar = dy .* ivar
        if g !== nothing
            dg = sum(x_mu .* dyivar, dims=dims)
            db = sum(dy, dims=dims)
            dyivar .*= g
        else
            dg, db = nothing, nothing
        end
        m = prod(d->size(x,d), dims) # size(x, dims...))
        dsigma2 = -sum(dyivar .* x_mu .* ivar.^2, dims=dims) ./ 2
        dmu = -sum(dyivar, dims=dims) .- 2dsigma2 .* sum(x_mu, dims=dims) ./ m
        dx = dyivar .+ dsigma2 .* 2x_mu ./ m .+ dmu ./ m
    else #same reasoning with the gpu version
        ivar = 1 ./ sqrt.(moments.var .+ eps)
        dx = (g !== nothing) ? (dy .* g .* ivar) : (dy .* ivar)
        if g !== nothing
            dg = sum(dy .* (x .- moments.mean) .* ivar, dims=dims)
            db = sum(dy, dims=dims)
        else
            dg, db = nothing, nothing
        end
    end
    return dg, db, dx
end

function _get_cache_data(cache, x, eps)
    if cache !== nothing
        mu = cache.mean
        ivar = cache.ivar
    else
        mu = mean(x, dims=_reddims(x))
        ivar = 1 ./ sqrt.(var(x, dims=_reddims(x); mean=mu, correct=false) .+ eps)
    end
    return mu, ivar
end

@primitive batchnorm4(x;o...),dy batchnorm4x(x, dy; o...)
@zerograd batchnorm4x(x, dy; o...)
# 4d with affine
@primitive batchnorm4(g, b, x; o...),dy batchnorm4g(g, x, dy; o...) batchnorm4b(dy; o...) batchnorm4x(g, x, dy; o...)
@zerograd batchnorm4x(g, x, dy; o...)
@zerograd batchnorm4g(g, x, dy; o...)
@zerograd batchnorm4b(dy; o...)

# Implement batchnorm2 using batchnorm4, with autograd

function batchnorm2(g, b, x; moments=nothing, training=training(), o...)
    # TODO: This support should be added when needed
    if training == false && (isa(g, AutoGrad.Value) || isa(x, AutoGrad.Value) || isa(b, AutoGrad.Value))
        error("Test mode backward is not supported with 2d inputs")
    end
    @inline _pad4(x) = reshape(x, (1,1,size(x,1),size(x,2)))
     # process moments
    if moments !== nothing
        _lazy_init!(moments, x)
        moments.mean = _pad4(moments.mean)
        moments.var = _pad4(moments.var)
    end
    x = _pad4(x)
    args = (x,)
    if g !== nothing
        g = _pad4(g)
        b = _pad4(b)
        args = (g, b, x)
    end
    y = mat(batchnorm4(args...;
                       moments=moments,
                       training=training,
                       o...))
    if moments !== nothing
        moments.mean = mat(moments.mean)
        moments.var = mat(moments.var)
    end
    return y
end

batchnorm2(x;o...) = batchnorm2(nothing, nothing, x; o...)


text:
# https://stackoverflow.com/questions/58580436/how-to-make-use-of-threads-optional-in-a-julia-function
# https://slides.com/valentinchuravy/julia-parallelism#/

default_step(n) =  min(n, round(Int, 10*sqrt(n)))

function getrange(n)
    tid = Threads.threadid()
    nt = Threads.nthreads()
    d , r = divrem(n, nt)
    from = (tid - 1) * d + min(r, tid - 1) + 1
    to = from + d - 1 + (tid ≤ r ? 1 : 0)
    from:to
end

# --------------------
# Thread mappers
# --------------------

abstract type AbstractThreadMapper end

struct Mapper <: AbstractThreadMapper
    start::Threads.Atomic{Int}
    stop::Int
    step::Int
    function Mapper(n::Int, step::Int=default_step(n))
        0 < n < typemax(Int) || throw(DomainError())
        start = Threads.Atomic{Int}(0)
        stop = n
        return new(start, stop, step)
    end
end

struct Mapper2 <: AbstractThreadMapper
    iter::Threads.Atomic{Int}
    stop::Int
    step::Int
    function Mapper2(stop::Int, step::Int=default_step(stop))
        iter = Threads.Atomic{Int}(1)
        return new(iter, stop, step)
    end
end

stop(m::AbstractThreadMapper) = m.stop
step(m::AbstractThreadMapper) = m.step

start(m::Mapper) = m.start

iter(m::Mapper2) = m.iter
add_iter!(m::Mapper2) = Threads.atomic_add!(iter(m),1)

function nextrange!(m::Mapper)
    strt = Threads.atomic_add!(start(m), step(m))
    if strt >= stop(m)
        return nothing
    else
        a = strt+1
        b = min(strt + step(m), stop(m))
        return a:b
    end
end


text:
module DBI2DBAPI

using DBI
using DBAPI

type DBIDriver <: AbstractDatabaseInterface
end

type DBICursor <: AbstractDatabaseCursor{DBIDriver}
    handle::DatabaseHandle
end

type DBIConnection <: AbstractDatabaseConnection{DBIDriver}
    cursors::Tuple{DBICursor}
    dbargs::Tuple{Any}
end

@doc """
The first argument needs to be

Passes the trailing args list on directly

kwargs are completely ignored but included for completeness
""" ->
function connect(::Type{DBIDriver}, drv::Type{T}, args...; kwargs...)
    return DBIConnection(drv, args)
end

@doc """
Internally, closes all open cursors
""" ->
function close(conn::DBIConnection)
    map(close, conn.cursors)
end

@doc """
Does nothing since DBI is not transactional.
""" ->
function commit(conn::DBIConnection)
end

@doc """
Does nothing since DBI is not transactional
""" ->
function rollback(conn::DBIConnection)
end

@doc """
Emulates cursors by using seperate connections for each cursor
""" ->
function cursor(conn::DBIConnection)
    cur = DBICursor(connect(typeof(conn.handle), conn.dbargs...)
    push!(cursors, cur)
    return cur
end

# callproc is left unimplemented

@doc """
Internally closes the connection the cursor contains
""" ->
function close(cur::DBICursor)
    disconnect(cur.handle)
end

@doc """
Similar to MySQLdb from Python the string is used to format the data in the
tuple preventing SQL injection attacks
""" ->
function execute(cur::DBICursor, sql::String; parameters::Tuple{String})
    prepare(cur.handle, sql)
end

@doc """
Similar to MySQLdb from Python the string is used to format the data in the
tuple preventing SQL injection attacks

TODO, implement parameter
""" ->
function executemany(cur::DBICursor, sql::String, seq_of_parameters::Tuple{String})
    stmt = prepare(cur.handle, sql)
    map(pars -> execute(stmt), seq_of_parameters)
end

end # module


text:
"""
Flexible binning of univariate arguments. See `ordered_bins`.
"""
module OrderedBinning

export ordered_bins

using ArgCheck: @argcheck
using UnPack: @unpack
using DocStringExtensions: SIGNATURES

struct OrderedBins{V<:AbstractVector,T}
    boundaries::V
    strict::Bool
    tolerance::T
    function OrderedBins(boundaries::V, strict::Bool,
                         tolerance::T) where {V<:AbstractVector,T}
        @argcheck issorted(boundaries)
        @argcheck tolerance ≥ 0
        new{V,T}(boundaries, strict, tolerance)
    end
end

"""
$(SIGNATURES)

Return a callable that places its argument in to bins according to `boundaries`, which
should be sorted.

Specifically, let `mi = first(x)` and `ma = last(x)`, and `i` is the value returned by the
callable. Then

- if `mi ≤ x ≤ ma`, `boundaries[i] ≤ x ≤ boundaries[i + 1]` holds. The last inequality is
  strict when `x < ma`.

- otherwise, `tolerance` extends the bins on the edge, ie if `x ≥ mi - tolerance`,
  `i == firstindex(boundaries)`, while if `x ≤ ma + tolerance`,
  `i == lastindex(boundaries) - 1`.

- finally, outside these ranges an error is thrown if `strict == true`, otherwise
  `i == firstindex(boundaries) - 1` or `i == lastindex(boundaries)`.

`tolerance` is meant for containing floating point error.

# Examples

```jldoctest
julia> ob = ordered_bins(0:3; strict = true, tolerance = 0.5)
OrderedBinning.OrderedBins{UnitRange{Int64}, Float64}(0:3, true, 0.5)

julia> ob(0)
1

julia> ob(3.5)
3
"""
function ordered_bins(boundaries; strict::Bool = true, tolerance = 0)
    OrderedBins(boundaries, strict, tolerance)
end

function (ob::OrderedBins)(x)
    @unpack boundaries, strict, tolerance = ob
    lo, hi = firstindex(boundaries), lastindex(boundaries)
    mi, ma = first(boundaries), last(boundaries)
    if x < mi
        if x ≥ mi - tolerance
            lo
        else
            strict && throw(DomainError(x, "below lowest boundary - tolerance"))
            lo - oftype(lo, 1)
        end
    elseif x ≥ ma
        if x ≤ ma + tolerance
            hi - 1
        else
            strict && throw(DomainError(x, "above highest boundary + tolerance"))
            hi
        end
    else
        searchsortedlast(ob.boundaries, x)
    end
end

end # module