text:
# TODO: finish cpu implementation.

### Size chart (Julia sizes for CUDNN calls)
# Note: For Julia calls, x and y do not need the initial 1 dimension and B,T are optional.
#
# x: (1,X,B,T) where X = inputSize, B = miniBatch, T = seqLength
# xDesc: Array of T (1,X,B) descriptors
# y: (1,Y,B,T) where Y = hiddenSize * (bidirectional ? 2 : 1)
# yDesc: Array of T (1,Y,B) descriptors
# w: (1,1,W) where W = cudnnGetRNNParamsSize()
# hx,cx,hy,cy: (H,B,L) where H = hidden size, L = numLayers * (bidirectional ? 2 : 1)
#
# Note: cudnn docs say min tensor dims 4 but RNN_example.cu uses 3D tensors

"""
    rnn = RNN(inputSize, hiddenSize; opts...)
    rnn(x; batchSizes) => y
    rnn.h, rnn.c  # hidden and cell states

`RNN` returns a callable RNN object `rnn`. Given a minibatch of sequences `x`, `rnn(x)`
returns `y`, the hidden states of the final layer for each time step. `rnn.h` and `rnn.c`
fields can be used to set the initial hidden states and read the final hidden states of all
layers.  Note that the final time step of `y` always contains the final hidden state of the
last layer, equivalent to `rnn.h` for a single layer network.

**Dimensions:** The input `x` can be 1, 2, or 3 dimensional and `y` will have the same number
of dimensions as `x`. size(x)=(X,[B,T]) and size(y)=(H/2H,[B,T]) where X is inputSize, B is
batchSize, T is seqLength, H is hiddenSize, 2H is for bidirectional RNNs. By default a 1-D `x`
represents a single instance for a single time step, a 2-D `x` represents a single minibatch
for a single time step, and a 3-D `x` represents a sequence of identically sized minibatches
for multiple time steps. The output `y` gives the hidden state (of the final layer for
multi-layer RNNs) for each time step. The fields `rnn.h` and `rnn.c` represent the hidden
states of all layers in a single time step and have size (H,B,L/2L) where L is numLayers and
2L is for bidirectional RNNs.

**batchSizes:** If `batchSizes=nothing` (default), all sequences in a minibatch are assumed to
be the same length. If `batchSizes` is an array of (non-increasing) integers, it gives us the
batch size for each time step (allowing different sequences in the minibatch to have different
lengths). In this case `x` will typically be 2-D with the second dimension representing
variable size batches for time steps. If `batchSizes` is used, `sum(batchSizes)` should equal
`length(x) ÷ size(x,1)`. When the batch size is different in every time step, hidden states
will have size (H,B,L/2L) where B is always the size of the first (largest) minibatch.

**Hidden states:** The hidden and cell states are kept in `rnn.h` and `rnn.c` fields (the cell
state is only used by LSTM). They can be initialized during construction using the `h` and `c`
keyword arguments, or modified later by direct assignment. Valid values are `nothing`
(default), `0`, or an array of the right type and size possibly wrapped in a `Param`. If the
value is `nothing` the initial state is assumed to be zero and the final state is discarded
keeping the value `nothing`. If the value is `0` the initial state is assumed to be zero and
`0` is replaced by the final state on return. If the value is a valid state, it is used as the
initial state and is replaced by the final state on return.

In a differentiation context the returned final hidden states will be wrapped in `Result`
types. This is necessary if the same RNN object is to be called multiple times in a single
iteration. Between iterations (i.e. after diff/update) the hidden states need to be unboxed
with e.g. `rnn.h = value(rnn.h)` to prevent spurious dependencies. This happens automatically
during the backward pass for GPU RNNs but needs to be done manually for CPU RNNs. See the
[CharLM Tutorial](https://github.com/denizyuret/Knet.jl/blob/master/tutorial/80.charlm.ipynb)
for an example.

**Keyword arguments for RNN:**
- `h=nothing`: Initial hidden state.
- `c=nothing`: Initial cell state.
- `rnnType=:lstm` Type of RNN: One of :relu, :tanh, :lstm, :gru.
- `numLayers=1`: Number of RNN layers.
- `bidirectional=false`: Create a bidirectional RNN if `true`.
- `dropout=0`: Dropout probability. Applied to input and between layers.
- `skipInput=false`: Do not multiply the input with a matrix if `true`.
- `dataType=Float32`: Data type to use for weights.
- `algo=0`: Algorithm to use, see CUDNN docs for details.
- `seed=0`: Random number seed for dropout. Uses `time()` if 0.
- `winit=xavier`: Weight initialization method for matrices.
- `binit=zeros`: Weight initialization method for bias vectors.
- `usegpu=(gpu()>=0)`: GPU used by default if one exists.

**Formulas:** RNNs compute the output h[t] for a given iteration from the recurrent input
h[t-1] and the previous layer input x[t] given matrices W, R and biases bW, bR from the
following equations:

`:relu` and `:tanh`: Single gate RNN with activation function f:

    h[t] = f(W * x[t] .+ R * h[t-1] .+ bW .+ bR)

`:gru`: Gated recurrent unit:

    i[t] = sigm(Wi * x[t] .+ Ri * h[t-1] .+ bWi .+ bRi) # input gate
    r[t] = sigm(Wr * x[t] .+ Rr * h[t-1] .+ bWr .+ bRr) # reset gate
    n[t] = tanh(Wn * x[t] .+ r[t] .* (Rn * h[t-1] .+ bRn) .+ bWn) # new gate
    h[t] = (1 - i[t]) .* n[t] .+ i[t] .* h[t-1]

`:lstm`: Long short term memory unit with no peephole connections:

    i[t] = sigm(Wi * x[t] .+ Ri * h[t-1] .+ bWi .+ bRi) # input gate
    f[t] = sigm(Wf * x[t] .+ Rf * h[t-1] .+ bWf .+ bRf) # forget gate
    o[t] = sigm(Wo * x[t] .+ Ro * h[t-1] .+ bWo .+ bRo) # output gate
    n[t] = tanh(Wn * x[t] .+ Rn * h[t-1] .+ bWn .+ bRn) # new gate
    c[t] = f[t] .* c[t-1] .+ i[t] .* n[t]               # cell output
    h[t] = o[t] .* tanh(c[t])
"""
RNN

"RNN descriptor"
mutable struct RD; ptr; end

"Dropout descriptor"
mutable struct DD; ptr; states::KnetArray{UInt8,1}; end

mutable struct RNN{T}
    w::Param{T}
    h
    c
    inputSize::Cint
    hiddenSize::Cint
    numLayers::Cint
    dropout::Float64
    seed::Culonglong
    inputMode::Cint    # CUDNN_LINEAR_INPUT = 0, CUDNN_SKIP_INPUT = 1
    direction::Cint    # CUDNN_UNIDIRECTIONAL = 0, CUDNN_BIDIRECTIONAL = 1
    mode::Cint         # CUDNN_RNN_RELU = 0, CUDNN_RNN_TANH = 1, CUDNN_LSTM = 2, CUDNN_GRU = 3
    algo::Cint         # CUDNN_RNN_ALGO_STANDARD = 0, CUDNN_RNN_ALGO_PERSIST_STATIC = 1, CUDNN_RNN_ALGO_PERSIST_DYNAMIC = 2
    dataType::DataType # CUDNN_DATA_FLOAT  = 0, CUDNN_DATA_DOUBLE = 1, CUDNN_DATA_HALF   = 2
    rnnDesc::Union{RD,Nothing}
    dropoutDesc::Union{DD,Nothing}
    dx
    dhx
    dcx
end

function RNN(inputSize, hiddenSize;
             h=nothing, c=nothing,
             handle=gethandle(),
             numLayers=1,
             dropout=0.0,
             skipInput=false,     # CUDNN_LINEAR_INPUT = 0, CUDNN_SKIP_INPUT = 1
             bidirectional=false, # CUDNN_UNIDIRECTIONAL = 0, CUDNN_BIDIRECTIONAL = 1
             rnnType=:lstm,       # CUDNN_RNN_RELU = 0, CUDNN_RNN_TANH = 1, CUDNN_LSTM = 2, CUDNN_GRU = 3
             dataType=Float32,    # CUDNN_DATA_FLOAT  = 0, CUDNN_DATA_DOUBLE = 1, CUDNN_DATA_HALF   = 2
             algo=0,              # CUDNN_RNN_ALGO_STANDARD = 0, CUDNN_RNN_ALGO_PERSIST_STATIC = 1, CUDNN_RNN_ALGO_PERSIST_DYNAMIC = 2
             seed=0,              # seed=0 for random init, positive integer for replicability
             winit=xavier,
             binit=zeros,
             usegpu=(gpu()>=0),
             )
    inputSize = Cint(inputSize)
    hiddenSize = Cint(hiddenSize)
    numLayers = Cint(numLayers)
    dropout = Float64(dropout)
    seed = Culonglong(seed)
    inputMode = Cint(skipInput ? 1 : 0)
    direction = Cint(bidirectional ? 1 : 0)
    mode = findfirst(isequal(rnnType), (:relu,:tanh,:lstm,:gru))
    if mode == nothing; error("RNN: Valid modes are :relu,:tanh,:lstm,:gru"); end
    mode = Cint(mode - 1)
    algo = Cint(algo)
    @assert dataType isa DataType
    if usegpu
        dropoutDesc = DD(handle=handle,dropout=dropout,seed=seed) # Need to keep dropoutDesc in RNN so it does not get gc'ed.
        rnnDesc = RD(hiddenSize,numLayers,dropoutDesc,inputMode,direction,mode,algo,dataType)
        w = Param(KnetArray{dataType}(undef,1,1,1)) # to get the right type
        r = RNN(w,h,c,inputSize,hiddenSize,numLayers,dropout,seed,inputMode,direction,mode,algo,dataType,rnnDesc,dropoutDesc,nothing,nothing,nothing)
        r.w = Param(KnetArray{dataType}(undef,1,1,cudnnGetRNNParamsSize(r)))
    else
        w = Param(Array{dataType}(undef,1,1,1)) # to get the right type
        r = RNN(w,h,c,inputSize,hiddenSize,numLayers,dropout,seed,inputMode,direction,mode,algo,dataType,nothing,nothing,nothing,nothing,nothing)
        r.w = param(Array{dataType}(undef,1,1,getRNNParamsSize(r)))
    end
    for a in rnnparams(r; handle=handle, useview=true)
        if a == nothing
            continue
        elseif ndims(a) == 2
            copyto!(a, winit(dataType, size(a)))
        elseif ndims(a) == 1
            copyto!(a, binit(dataType, size(a)))
        else
            error("Invalid RNN param $(summary(a))")
        end
    end
    return r
end

function (r::RNN)(x; batchSizes=nothing)
    # Check type/dims of inputs
    WTYPE = typeof(vec(value(r.w)))
    @assert vec(value(x)) isa WTYPE
    @assert ndims(x) <= 3
    @assert size(x,1) == r.inputSize
    HSIZE = (r.hiddenSize, batchSizes == nothing ? size(x,2) : batchSizes[1], r.numLayers * (r.direction + 1))
    @assert r.h == nothing || r.h == 0 || (vec(value(r.h)) isa WTYPE && ndims(r.h) <= 3 && (size(r.h,1),size(r.h,2),size(r.h,3)) == HSIZE)
    @assert r.c == nothing || r.c == 0 || (vec(value(r.c)) isa WTYPE && ndims(r.c) <= 3 && (size(r.c,1),size(r.c,2),size(r.c,3)) == HSIZE)
    # apply dropout to input: rnnforw only applies it between layers.
    x = dropout(x,r.dropout)
    # use hidden inputs unless nothing or 0
    hx = (r.h == 0 ? nothing : r.h)
    cx = (r.c == 0 ? nothing : r.c)
    # produce hidden outputs if hidden inputs != nothing
    hy = (r.h != nothing)
    cy = (r.mode == 2 && r.c != nothing)
    (y, hyout, cyout, rs) = rnnforw(r, r.w, x, hx, cx; hy=hy, cy=cy, batchSizes=batchSizes)
    # In a @diff context hyout and cyout will be Results
    # This is necessary to keep dependencies if r is run several times during one forw pass
    # However unless Results are cleared up before next iteration we will run out of memory
    # I do the clearing in the back function when I know the iteration is done for sure.
    if hy; r.h = hyout; end
    if cy; r.c = cyout; end
    return y
end

function show(io::IO, r::RNN)
    print(io, ("RNNRELU","RNNTANH","LSTM","GRU")[r.mode+1], "(input=", r.inputSize, ",hidden=", r.hiddenSize)
    if r.direction == 1; print(io, ",bidirectional"); end
    if r.numLayers > 1; print(io, ",layers=", r.numLayers); end
    if r.dropout > 0; print(io, ",dropout=", r.dropout); end
    if r.inputMode == 1; print(io, ",skipinput"); end
    if r.dataType != Float32; print(io, ',', r.dataType); end
    print(io, ')')
end

Base.unsafe_convert(::Type{Cptr}, dd::DD)=dd.ptr

function DD(; handle=cudnnhandle(), dropout=0.0, seed=0, o...)
    if seed==0; seed=floor(Culonglong,time()); end
    d = Cptr[0]; s = Csize_t[0] # TODO: Can multiple RNNs share dropout descriptors? Can dropout probability be changed?
    @cudnn(cudnnCreateDropoutDescriptor,(Ptr{Cptr},),d)
    @cudnn(cudnnDropoutGetStatesSize,(Cptr,Ptr{Csize_t}),handle,s)
    states = KnetArray{UInt8}(undef,s[1]) # TODO: Can this be shared? 638976 bytes.
    @cudnn(cudnnSetDropoutDescriptor,(Cptr,Cptr,Cfloat,Cptr,Csize_t,Culonglong),
          d[1],handle,dropout,states,bytes(states),seed)
    dd = DD(d[1],states)
    finalizer(x->@cudnn(cudnnDestroyDropoutDescriptor,(Cptr,),x.ptr),dd)
    return dd
end

Base.unsafe_convert(::Type{Cptr}, rd::RD)=rd.ptr

function RD()
    d = Cptr[0]
    @cudnn(cudnnCreateRNNDescriptor,(Ptr{Cptr},),d)
    rd = RD(d[1])
    finalizer(x->@cudnn(cudnnDestroyRNNDescriptor,(Cptr,),x.ptr),rd)
    return rd
end

function RD(hiddenSize,numLayers,dropoutDesc,inputMode,direction,mode,algo,dataType; handle=gethandle())
    rnnDesc = RD()
    if cudnnVersion >= 7000
        @cudnn(cudnnSetRNNDescriptor,(Cptr,Cptr,Cint,Cint,Cptr,Cint,Cint,Cint,Cint,Cint),
              handle,rnnDesc,hiddenSize,numLayers,dropoutDesc,inputMode,direction,mode,algo,DT(dataType))
    elseif cudnnVersion >= 6000
        @cudnn(cudnnSetRNNDescriptor_v6,(Cptr,Cptr,Cint,Cint,Cptr,Cint,Cint,Cint,Cint,Cint),
              handle,rnnDesc,hiddenSize,numLayers,dropoutDesc,inputMode,direction,mode,algo,DT(dataType))
    elseif cudnnVersion >= 5000
        @cudnn(cudnnSetRNNDescriptor,(Cptr,Cint,Cint,Cptr,Cint,Cint,Cint,Cint),
              rnnDesc,hiddenSize,numLayers,dropoutDesc,inputMode,direction,mode,DT(dataType))
    else
        error("CUDNN $cudnnVersion does not support RNNs")
    end
    return rnnDesc
end

rnnids(r) = (r.mode == 2 ? 8 : r.mode == 3 ? 6 : 2)

function cudnnGetRNNParamsSize(r::RNN; handle=cudnnhandle())
    if r.rnnDesc != nothing
        res = Csize_t[0]
        xDesc = TD(r.dataType, 1, r.inputSize, 1)    # xDesc: (1,X,B) where X = inputSize, B is ignored, so assume 1
        @cudnn(cudnnGetRNNParamsSize,
              # handle, rnndesc, xdesc, result, dataType
              (Cptr,  Cptr, Cptr, Ptr{Csize_t}, UInt32),
              handle, r.rnnDesc, xDesc, res, DT(r.dataType))
        div(res[1], sizeof(r.dataType))
    else # on CPU, so we guess the size
        X,H,L,I = r.inputSize, r.hiddenSize, r.numLayers, rnnids(r)
        biases = L*I
        inputMatrices = (r.inputMode == 1 ? 0 : r.direction == 1 ? I : div(I,2))
        hiddenMatrices = (r.direction == 1 ? (L-1)*I : (L-1)*I + div(I,2))
        biases * H + inputMatrices * X * H + hiddenMatrices * H * H
    end
end

function getRNNParamsSize(r::RNN)
    whidden = r.hiddenSize * r.hiddenSize
    winput =  r.inputMode == 1 ? 0 : r.hiddenSize * r.inputSize
    bhidden = r.hiddenSize
    binput =  bhidden
    coef = (r.mode == 2 ? 4 : r.mode == 3 ? 3 : 1) * (1 + r.direction)
    nparams = 0
    for i = 1:r.numLayers
        nparams += coef * (whidden + winput + bhidden + binput)
        winput = (1 + r.direction) * whidden
        binput = bhidden
    end
    nparams
end

"Keeps an array of 3D tensor descriptors"
mutable struct TDs; pvec::Vector{Cptr}; xDesc::Vector{TD}; end     # Keep xDesc in TDs so it does not get gc'ed

Base.unsafe_convert(::Type{Ptr{Cptr}}, tds::TDs)=pointer(tds.pvec)
Base.length(tds::TDs)=length(tds.pvec)

function TDs(x::KnetArray{A},::Nothing) where {A} # Treat x: (X,B?,T?) as a 4D array: (1,X,B,T)
    xDesc = TD(A,1,size(x,1),size(x,2)) # we can use a single xDesc
    pvec = Vector{Cptr}(undef, size(x,3))
    pvec[:] .= xDesc.ptr
    return TDs(pvec, [xDesc])
end

function TDs(x::KnetArray{A},batchSizes) where {A} # x: (X,B*), batchSizes gives us Bt sizes
    @assert sum(batchSizes) == div(length(x),size(x,1))
    X = size(x,1)
    xs = [ TD(A,1,X,B) for B in batchSizes ]
    ps = [ xd.ptr for xd in xs ]
    return TDs(ps,xs)
end

function TD3(a::KnetArray) # Treat a as a 3D array, pad from right
    n = ndims(a)
    if n==3; TD(a)
    elseif n==2; TD(reshape(a, size(a,1), size(a,2), 1))
    elseif n==1; TD(reshape(a, size(a,1), 1, 1))
    else; throw(DimensionMismatch())
    end
end

function FD3(a::KnetArray) # Treat a as a 3D array, pad from left
    n = ndims(a)
    if n==3; FD(a)
    elseif n==2; FD(reshape(a, 1, size(a,1), size(a,2)))
    elseif n==1; FD(reshape(a, 1, 1, size(a,1)))
    else; throw(DimensionMismatch())
    end
end

function cudnnGetRNNWorkspaceSize(rd::RD, tds::TDs; handle=cudnnhandle())
    res = Csize_t[1]
    @cudnn(cudnnGetRNNWorkspaceSize,
          # handle, rnndesc, seqLength, xdesc, res        ,
          (Cptr,  Cptr, Cint, Ptr{Cptr}, Ptr{Csize_t}),
          handle, rd, length(tds), tds, res)
    return Int(res[1])
end

function cudnnGetRNNTrainingReserveSize(rd::RD, tds::TDs; handle=cudnnhandle())
    res = Csize_t[1]
    @cudnn(cudnnGetRNNTrainingReserveSize,
          # handle, rnndesc, seqLength, xdesc, res        ,
          (Cptr,  Cptr, Cint, Ptr{Cptr}, Ptr{Csize_t}),
          handle, rd, length(tds), tds, res)
    return Int(res[1])
end

# Return eltype,size
function cudnnGetFilterNdDescriptor(wDesc::FD; nbDimsRequested = 8)
    dataType = Cint[0]
    format = Cint[0]
    nbDims = Cint[0]
    filterDimA = Vector{Cint}(undef,nbDimsRequested)
    @cudnn(cudnnGetFilterNdDescriptor,
          (Cptr, Cint, Ptr{UInt32}, Ptr{UInt32}, Ptr{Cint}, Ptr{Cint}),
          wDesc, nbDimsRequested, dataType, format, nbDims, filterDimA)
    if nbDims[1] > nbDimsRequested
        cudnnGetFilterNdDescriptor(wDesc::FD; nbDimsRequested = nbDims[1])
    else
        (Float32,Float64,Float16)[1+dataType[1]],
        (filterDimA[nbDims[1]:-1:1]...,)
    end
end

# call gpu everytime to reflect device changes
gethandle() = gpu() >= 0 ? cudnnhandle() : nothing

"""
    rnnparam(r::RNN, layer, id, param)

Return a single weight matrix or bias vector as a slice of RNN weights.

Valid `layer` values:
* For unidirectional RNNs 1:numLayers
* For bidirectional RNNs 1:2*numLayers, forw and back layers alternate.

Valid `id` values:
* For RELU and TANH RNNs, input = 1, hidden = 2.
* For GRU reset = 1,4; update = 2,5; newmem = 3,6; 1:3 for input, 4:6 for hidden
* For LSTM inputgate = 1,5; forget = 2,6; newmem = 3,7; output = 4,8; 1:4 for input, 5:8 for hidden

Valid `param` values:
* Return the weight matrix (transposed!) if `param==1`.
* Return the bias vector if `param==2`.

The effect of skipInput: Let I=1 for RELU/TANH, 1:3 for GRU, 1:4 for LSTM
* For skipInput=false (default), rnnparam(r,1,I,1) is a (inputSize,hiddenSize) matrix.
* For skipInput=true, rnnparam(r,1,I,1) is `nothing`.
* For bidirectional, the same applies to rnnparam(r,2,I,1): the first back layer.
"""
function rnnparam(r::RNN, layer::Integer, id::Integer, par::Integer; handle=gethandle(), useview=false)
    ((1 <= par <= 2) &&
     ((r.direction == 0 && 1 <= layer <= r.numLayers) ||
      (r.direction == 1 && 1 <= layer <= 2*r.numLayers)) &&
     ((r.mode == 0 && 1 <= id <= 2) ||
      (r.mode == 1 && 1 <= id <= 2) ||
      (r.mode == 2 && 1 <= id <= 8) ||
      (r.mode == 3 && 1 <= id <= 6))) || error("Bad parameter index")
    i1 = i2 = len = 0
    w = value(r.w)
    if isa(w, KnetArray)
        T = eltype(w)
        xDesc = TD(T,1,r.inputSize,1)
        wDesc = FD(T,1,1,length(w))
        paramDesc = FD(T,1,1,1,1)
        param = Cptr[0]
        if par == 1 # matrix
            @cudnn(cudnnGetRNNLinLayerMatrixParams,
                  (Cptr, Cptr, Cint, #handle,rdesc, layer
                   Cptr, Cptr, Cptr, #xDesc, wDesc, w
                   Cint, Cptr, Ptr{Cptr}), #lid, lmatdesc, linlayermat
                  handle, r.rnnDesc, layer-1,
                  xDesc, wDesc, w,
                  id-1, paramDesc, param)
        else # bias
            @cudnn(cudnnGetRNNLinLayerBiasParams,
                  (Cptr, Cptr, Cint, #handle,rdesc, layer
                   Cptr, Cptr, Cptr, #xDesc, wDesc, w
                   Cint, Cptr, Ptr{Cptr}), #lid, lmatdesc, linlayermat
                  handle, r.rnnDesc, layer-1,
                  xDesc, wDesc, w,
                  id-1, paramDesc, param)
        end
        dt,sz = cudnnGetFilterNdDescriptor(paramDesc)
        len = prod(sz)
        if len > 1
            i1 = 1 + div(Int(param[1] - pointer(w)), sizeof(T))
            i2 = i1 + len - 1
        end
    else
        # guess i1,i2,len from layer,id,par and rnn specs
        ids = rnnids(r)
        if par == 1 # matrix
            # input matrices are (inputSize,hiddenSize) all others (hiddenSize,hiddenSize)
            # if inputMode==1 then input matrices are empty
            # I=1 for RELU/TANH, 1:3 for GRU, 1:4 for LSTM are input matrices with L=1 for uni and L=1,2 for bi
            inputLayer(r,l)=(l==1 || (l==2 && r.direction==1))
            X, H = r.inputSize, r.hiddenSize
            XH, HH = X*H, H*H
            inputIds = div(ids,2)
            for l = 1:layer, i = 1:ids
                if inputLayer(r,l) && i <= inputIds
                    len = (r.inputMode==0 ? XH : 0)
                elseif l>2 && r.direction==1 && i <= div(ids, 2)
                    # bidirectional weight
                    len = 2HH
                else
                    len = HH
                end
                i2 += len
                if l==layer && i==id; break; end
            end
            if len==0; len=1; end # cudnn uses size (1,1,1) for empty weights
            i1 = i2 - len + 1
        else # bias
            # all biases are length=hidden and there are always numLayers * numIds of them
            len = r.hiddenSize
            layers = r.numLayers * (r.direction == 1 ? 2 : 1)
            i1 = 1 + length(w) - layers * ids * len + ((id-1) + ids * (layer-1)) * len
            i2 = i1 + len - 1
        end
    end
    @inline access(a, rng) = (isa(a, KnetArray) || ~useview) ? a[rng] : view(a, rng)
    if len == 1 # empty weights when inputMode=1 show up as size (1,1,1)
        nothing
    elseif par == 1 # matrix
        h = Int(r.hiddenSize)
        reshape(access(r.w, i1:i2),
                (div(len,h),h)) # weight matrices are transposed
    else # bias
        access(r.w, i1:i2)
    end
end

"""
    rnnparams(r::RNN)

Return the RNN parameters as an Array{Any}.

The order of params returned (subject to change):
* All weight matrices come before all bias vectors.
* Matrices and biases are sorted lexically based on (layer,id).
* See @doc rnnparam for valid layer and id values.
* Input multiplying matrices are `nothing` if r.inputMode = 1.
"""
function rnnparams(r::RNN; handle=gethandle(), useview=false)
    layers = r.numLayers * (r.direction == 1 ? 2 : 1)
    ids = rnnids(r)
    ws = []
    for m in (1,2)
        for l in 1:layers
            for i in 1:ids
                push!(ws, rnnparam(r, l, i, m; handle=handle, useview=useview))
            end
        end
    end
    return ws
end

function rnnforw(r::RNN{KnetArray{T,3}}, w::KnetArray{T,3}, x::KnetArray{T},
                 hx::Union{KnetArray{T},Nothing}=nothing,
                 cx::Union{KnetArray{T},Nothing}=nothing;
                 handle=cudnnhandle(),
                 batchSizes=nothing,
                 hy = (hx != nothing),
                 cy = (cx != nothing && r.mode == 2),
                 ) where {T}
    @assert w === value(r.w)
    # Input descriptors
    if size(x,1) != r.inputSize
        throw(DimensionMismatch("size(x,1)=$(size(x,1)) does not match r.inputSize=$(r.inputSize)"))
    end
    seqLength = batchSizes==nothing ? size(x,3) : length(batchSizes) # (X,B,T) or (X,B+) with batchSizes
    wDesc = FD3(w)              # (1,1,W)
    xtds = TDs(x,batchSizes)    # (1,X,Bt) x T
    isnothing(a) = a == nothing || a == C_NULL
    if hx==nothing; hx=hxDesc=C_NULL; else; hxDesc=TD3(hx); end # (H,B,L/2L)
    if cx==nothing || r.mode != 2; cx=cxDesc=C_NULL; else; cxDesc=TD3(cx); end

    # Output arrays and descriptors
    ysize = collect(size(x))
    ysize[1] = r.hiddenSize * (r.direction == 1 ? 2 : 1)
    y = similar(x, ysize...)    # (H/2H,B,T) or (H/2H,B+) -- y mirrors x except for the first dimension
    ytds = TDs(y,batchSizes)    # (1,H/2H,Bt) x T

    # Optionally output hidden and cell of last step
    hyout = hyDesc = cyout = cyDesc = C_NULL
    if hy || cy
        firstBatchSize = batchSizes==nothing ? size(x,2) : batchSizes[1]
        hsize = (Int(r.hiddenSize), Int(firstBatchSize), Int(r.numLayers * (r.direction == 1 ? 2 : 1))) # (H,B,L/2L)
        if hy; hyout=similar(y,hsize); hyDesc=TD3(hyout); end
        if cy && r.mode==2; cyout=similar(y,hsize); cyDesc=TD3(cyout); end
        if !isnothing(hx) && any(size(hx,i)!=hsize[i] for i=1:3) # compare one by one in case hx is 1-D or 2-D
            throw(DimensionMismatch("size(hx)=$(size(hx)) does not match hsize=$(hsize)"))
        end
        if !isnothing(cx) && r.mode == 2 && any(size(cx,i)!=hsize[i] for i=1:3)
            throw(DimensionMismatch("size(cx)=$(size(cx)) does not match hsize=$(hsize)"))
        end
    end

    # workSpace and reserveSpace
    wss = cudnnGetRNNWorkspaceSize(r.rnnDesc, xtds; handle=handle)
    ws = cudnnWorkSpace(wss)

    if training()
        rss = cudnnGetRNNTrainingReserveSize(r.rnnDesc, xtds; handle=handle)
        rs = KnetArray{UInt8}(undef,rss)
        @cudnn(cudnnRNNForwardTraining,
              (Cptr, Cptr, Cint,  # handle,rnnDesc,seqLength
               Ptr{Cptr}, Ptr{T}, #x
               Cptr, Ptr{T}, #hx
               Cptr, Ptr{T}, #cx
               Cptr, Ptr{T}, #w
               Ptr{Cptr}, Ptr{T}, #y
               Cptr, Ptr{T}, #hy
               Cptr, Ptr{T}, #cy
               Cptr, Csize_t, #ws
               Cptr ,Csize_t#rs
               ),
              handle, r.rnnDesc, seqLength,
              xtds, x,
              hxDesc, hx,
              cxDesc, cx,
              wDesc, w,
              ytds, y,
              hyDesc, hyout,
              cyDesc, cyout,
              ws, wss,
              rs, rss)
    else
        rs = nothing
        @cudnn(cudnnRNNForwardInference,
              (Cptr, Cptr, Cint,  # handle,rnnDesc,seqLength
               Ptr{Cptr}, Ptr{T}, #x
               Cptr, Ptr{T}, #h
               Cptr, Ptr{T}, #c
               Cptr, Ptr{T}, #w
               Ptr{Cptr}, Ptr{T}, #y
               Cptr, Ptr{T}, #hy
               Cptr, Ptr{T}, #cy
               Cptr, Csize_t, #ws
               ),
              handle, r.rnnDesc, seqLength,
              xtds, x,
              hxDesc, hx,
              cxDesc, cx,
              wDesc, w,
              ytds, y,
              hyDesc, hyout,
              cyDesc, cyout,
              ws, wss)
    end
    if hyout == C_NULL; hyout = nothing; end
    if cyout == C_NULL; cyout = nothing; end
    return y, hyout, cyout, rs
end

@primitive rnnforw(r::RNN, w::KnetArray, x...; o...),dy,y nothing rnnback2(dy,y,r,w,x...;o...) value(r).dx value(r).dhx value(r).dcx

function rnnback2(dt, t, r, w, x, hx=nothing, cx=nothing; o...)
    @assert r.w === w
    y,hy,cy,rs = value(t)
    dy,dhy,dcy,drs = value(dt)
    r=value(r); w=value(w); x=value(x); hx=value(hx); cx=value(cx)
    # To prevent dependencies to next iteration we need to clear the Result type from r.h,r.c
    # We can't do this during forward, because another forward may be run within the same iteration.
    # Doing it here is safe, means the iteration is done and we are taking gradients.
    # Note that this does not work on the cpu and these have to be cleaned by hand.
    # The cpu version is not a primitive and has no back function. (TODO: find better solution)
    r.h = value(r.h); r.c = value(r.c)
    rnnback(r, w, x, y, dy, hx, cx, dhy, dcy, rs; o...)
end

function rnnback(r::RNN, w::KnetArray{T}, x::KnetArray{T}, y::KnetArray{T},
                 dy, hx, cx, dhy, dcy, rs; handle=cudnnhandle(), batchSizes=nothing, o...) where {T}
    @assert value(r.w) === w
    # Input descriptors:
    seqLength = batchSizes==nothing ? size(x,3) : length(batchSizes) # (X,B,T) or (X,B+) with batchSizes
    wDesc = FD3(w)              # (1,1,W)
    xtds = TDs(x,batchSizes)    # (X,B,T) -> (1,X,B) x T
    ytds = TDs(y,batchSizes)    # (H/2H,B,T) -> (1,H/2H,B) x T
    # dytds = TDs(dy,batchSizes)  # we use ytds for dytds
    if dy == nothing; dy=zero(y); end
    if hx == nothing; hx=hxDesc=C_NULL; else; hxDesc=TD3(hx); end
    if cx == nothing || r.mode != 2; cx=cxDesc=C_NULL; else; cxDesc=TD3(cx); end
    if dhy == nothing; dhy=dhyDesc=C_NULL; else; dhyDesc=TD3(dhy); end
    if dcy == nothing || r.mode != 2; dcy=dcyDesc=C_NULL; else; dcyDesc=TD3(dcy); end

    # Output arrays and descriptors:
    dx = similar(x)             # (X,B,T) or (X,B+) with batchSizes
    # dxtds = TDs(dx,batchSizes)  # we use xtds here
    dw = zero(w)               # dw is used additively, so we need zero
    dwDesc = FD3(dw)
    if hx == C_NULL; dhx=dhxDesc=C_NULL; else; dhx=similar(hx); dhxDesc=TD3(dhx); end
    if cx == C_NULL; dcx=dcxDesc=C_NULL; else; dcx=similar(cx); dcxDesc=TD3(dcx); end

    # workSpace and reserveSpace
    ws = cudnnWorkSpace()
    wss = bytes(ws)
    rss = bytes(rs)

    # data backward
    @cudnn(cudnnRNNBackwardData,
          (Cptr, Cptr, Cint,  # handle, rnnDesc, seqLength
           Ptr{Cptr}, Ptr{T}, #y
           Ptr{Cptr}, Ptr{T}, #dy
           Cptr, Ptr{T}, #dhy
           Cptr, Ptr{T}, #dcy
           Cptr, Ptr{T}, #w
           Cptr, Ptr{T}, #hx
           Cptr, Ptr{T}, #cx
           Ptr{Cptr}, Ptr{T}, #dx
           Cptr, Ptr{T}, #dhx
           Cptr, Ptr{T}, #dcx
           Cptr, Csize_t, #ws
           Cptr, Csize_t), #rs
          # Use rtd with nullables
          handle, r.rnnDesc, seqLength,
          ytds, y,
          ytds, dy,
          dhyDesc, dhy,
          dcyDesc, dcy,
          wDesc, w,
          hxDesc, hx,
          cxDesc, cx,
          xtds, dx,
          dhxDesc, dhx,
          dcxDesc, dcx,
          ws, wss,
          rs, rss)
    # weights backward
    @cudnn(cudnnRNNBackwardWeights,
          (Cptr, Cptr, Cint,  # handle, rnnDesc, seqLength
           Ptr{Cptr}, Ptr{T}, #x
           Cptr, Ptr{T}, #hx
           Ptr{Cptr}, Ptr{T}, #y
           Cptr, Csize_t, #ws
           Cptr, Ptr{T}, #dw
           Ptr{Cptr}, Csize_t), #rs
          handle, r.rnnDesc, seqLength,
          xtds, x,
          hxDesc, hx,
          ytds, y,
          ws, wss,
          dwDesc, dw,
          rs, rss)
    # Update the cache
    if dhx==C_NULL; dhx=nothing; end
    if dcx==C_NULL; dcx=nothing; end
    r.dx, r.dhx, r.dcx = dx, dhx, dcx
    return dw
end

# CPU version
function rnnforw(r::RNN, w::AbstractArray{T}, x::AbstractArray{T},
                 hx::Union{AbstractArray{T},Nothing}=nothing,
                 cx::Union{AbstractArray{T},Nothing}=nothing;
                 # handle=cudnnhandle(), training=false,
                 batchSizes=nothing,
                 hy = (hx != nothing),
                 cy = (cx != nothing && r.mode == 2),
                 o...) where {T}
    @assert w === value(r.w)
    rnntest(r,w,x,hx,cx;batchSizes=batchSizes,hy=hy,cy=cy)
end

# rnnforw is an AutoGrad primitive for KnetArray, but a regular function for AbstractArray:
rnnforw(r::RNN, w::AutoGrad.Value{<:AbstractArray}, x...; o...) = rnntest(r,w,x...;o...)

# non-CUDNN cpu/gpu version
function rnntest(r::RNN, ws, x, hx=nothing, cx=nothing;
                 batchSizes=nothing,
                 hy = (hx != nothing),
                 cy = (cx != nothing && r.mode == 2),
                 o...)
    if batchSizes != nothing
        return rnntest_bs(batchSizes,r, ws, x, hx, cx; hy=hy, cy=cy, o...)
    end
    @assert value(r.w) === value(ws)
    w = rnnparams(r)
    X,B,T = (size(x,i) for i=1:3) # ndims(x) may be 1,2 or 3
    @assert X == r.inputSize
    Y = Int(r.hiddenSize * (r.direction == 1 ? 2 : 1))
    ysize = ntuple(i->(i==1 ? Y : size(x,i)), ndims(x)) # to match ndims(y) to ndims(x)
    H = Int(r.hiddenSize)
    #@assert (r.inputMode == 0 || H == X)
    L = Int(r.numLayers) * (r.direction == 1 ? 2 : 1)
    hsize = (H, B, L)
    @assert hx == nothing || eqsize(size(hx), hsize)
    @assert cx == nothing || eqsize(size(cx), hsize)
    h = hx==nothing ? fill!(similar(value(x),hsize),0) : hx
    #  hs = Array{Any}[ h[:,:,l] for l=1:L ]
    hs = Array{Any}(undef,L)
    for l = 1:L
        hs[l] = h[:,:,l]
    end
    ys = []
    direction = r.direction
    pdrop = r.dropout
    #=
    All complexity of bidirectional execution
    is packed inside this inline function.
    This causes code repetition, but  works w/o
    touching the existing unidirectional test code
    =#
    @inline bidirect(update_h!) = begin
        xl = x
        for l = 1:(1+direction):L
            skip = l==1 && r.inputMode==1
            hts = []
            if l>1; xl = dropout(xl, pdrop); end
            for t = 1:T
                for (i,ti) in zip([l, l+1], [t, T-t+1])
                    # this function updates h[i]
                    update_h!(xl, i, ti, skip)
                    push!(hts, hs[i])
                end
            end
            # construct the next layer input
            yforw = Array{Any}(hts[1:2:end-1])
            yback = Array{Any}(reverse(hts[2:2:end]))
            ybs = []
            for (yf, yb) in zip(yforw, yback)
                push!(ybs, vcat(yf, yb))
            end
            # now ybs contans (2 * hiddenSize, batchSize) matrices
            # so cat them to add time dimension
            xl = reshape(hcat(ybs...), (2r.hiddenSize, size(x,2), length(ybs)))
        end
        ys = xl
    end

    if r.mode <= 1
        #@assert r.inputMode == 0 || all(w[1:1+r.direction] .== nothing)
        f = r.mode == 0 ? relu : tanh
        if direction == 0
            for t = 1:T
                for l = 1:L
                    xl = l>1 ? dropout(hs[l-1], pdrop) : nothing
                    wx,wh,bx,bh = w[2l-1],w[2l],w[2L+2l-1],w[2L+2l]
                    wxt = (l > 1 ? wx' * xl : r.inputMode==0 ? wx' * x[:,:,t] : x[:,:,t])
                    hs[l] = f.(wxt .+ wh' * hs[l] .+ bx .+ bh)
                end
                push!(ys, hs[L])
            end
        else
            bidirect() do xl, i, ti, skip
                wx,wh,bx,bh = w[2i-1],w[2i],w[2L+2i-1],w[2L+2i]
                wxt =  skip ? xl[:,:,ti] : wx' * xl[:,:,ti]
                hs[i] = f.(wxt .+ wh' * hs[i] .+ bx .+ bh)
            end
        end
    elseif r.mode == 2           # LSTM
        #@assert r.inputMode == 0 || all(w[1:4*(1+r.direction)] .== nothing)
        # it = σ(Wixt + Riht-1 + bWi + bRi)
        # ft = σ(Wfxt + Rfht-1 + bWf + bRf)
        # ot = σ(Woxt + Roht-1 + bWo + bRo)
        # c't = tanh(Wcxt + Rcht-1 + bWc + bRc)
        # ct = ft◦ct-1 + it◦c't
        # ht = ot◦tanh(ct)
        c = cx==nothing ? fill!(similar(value(x),hsize),0) : cx
        cs = Array{Any}(undef,L)
        for l = 1:L
            cs[l] = c[:,:,l]
        end
        if direction == 0
            for t = 1:T
                for l = 1:L
                    xl = l>1 ? dropout(hs[l-1], pdrop) : nothing
                    Wi,Wf,Wc,Wo,Ri,Rf,Rc,Ro = w[1+8*(l-1):8l]
                    bWi,bWf,bWc,bWo,bRi,bRf,bRc,bRo = w[8L+1+8*(l-1):8L+8l]
                    Wixt = (l > 1 ? Wi' * xl : r.inputMode==0 ? Wi' * x[:,:,t] : x[:,:,t])
                    Wfxt = (l > 1 ? Wf' * xl : r.inputMode==0 ? Wf' * x[:,:,t] : x[:,:,t])
                    Wcxt = (l > 1 ? Wc' * xl : r.inputMode==0 ? Wc' * x[:,:,t] : x[:,:,t])
                    Woxt = (l > 1 ? Wo' * xl : r.inputMode==0 ? Wo' * x[:,:,t] : x[:,:,t])
                    it = sigm.(Wixt .+ Ri' * hs[l] .+ bWi .+ bRi)
                    ft = sigm.(Wfxt .+ Rf' * hs[l] .+ bWf .+ bRf)
                    ot = sigm.(Woxt .+ Ro' * hs[l] .+ bWo .+ bRo)
                    cn = tanh.(Wcxt .+ Rc' * hs[l] .+ bWc .+ bRc)
                    cs[l] = ft .* cs[l] .+ it .* cn
                    hs[l] = ot .* tanh.(cs[l])
                end
                push!(ys, hs[L])
            end
        else
            bidirect() do xl, i, ti, skip
                Wi,Wf,Wc,Wo,Ri,Rf,Rc,Ro = w[1+8*(i-1):8i]
                bWi,bWf,bWc,bWo,bRi,bRf,bRc,bRo = w[8L+1+8*(i-1):8L+8i]
                Wixt = skip ? xl[:,:,ti] : Wi' * xl[:,:,ti]
                Wfxt = skip ? xl[:,:,ti] : Wf' * xl[:,:,ti]
                Wcxt = skip ? xl[:,:,ti] : Wc' * xl[:,:,ti]
                Woxt = skip ? xl[:,:,ti] : Wo' * xl[:,:,ti]
                it = sigm.(Wixt .+ Ri' * hs[i] .+ bWi .+ bRi)
                ft = sigm.(Wfxt .+ Rf' * hs[i] .+ bWf .+ bRf)
                ot = sigm.(Woxt .+ Ro' * hs[i] .+ bWo .+ bRo)
                cn = tanh.(Wcxt .+ Rc' * hs[i] .+ bWc .+ bRc)
                cs[i] = ft .* cs[i] .+ it .* cn
                hs[i] = ot .* tanh.(cs[i])
            end
        end
    elseif r.mode == 3           # GRU
        #@assert r.inputMode == 0 || all(w[1:3*(1+r.direction)] .== nothing)
        # rt = σ(Wrxt + Rrht-1 + bWr + bRr)
        # it = σ(Wixt + Riht-1 + bWi + bRu)
        # h't = tanh(Whxt + rt◦(Rhht-1 + bRh) + bWh)
        # ht = (1 - it)◦h't + it◦ht-1
        if direction == 0
            for t = 1:T
                for l = 1:L
                    xl = l>1 ? dropout(hs[l-1], pdrop) : nothing
                    Wr,Wi,Wh,Rr,Ri,Rh = w[1+6*(l-1):6l]
                    bWr,bWi,bWh,bRr,bRi,bRh = w[6L+1+6*(l-1):6L+6l]
                    Wrxt = (l > 1 ? Wr' * xl : r.inputMode==0 ? Wr' * x[:,:,t] : x[:,:,t])
                    Wixt = (l > 1 ? Wi' * xl : r.inputMode==0 ? Wi' * x[:,:,t] : x[:,:,t])
                    Whxt = (l > 1 ? Wh' * xl : r.inputMode==0 ? Wh' * x[:,:,t] : x[:,:,t])
                    rt = sigm.(Wrxt .+ Rr' * hs[l] .+ bWr .+ bRr)
                    it = sigm.(Wixt .+ Ri' * hs[l] .+ bWi .+ bRi)
                    ht = tanh.(Whxt .+ rt .* (Rh' * hs[l] .+ bRh) .+ bWh)
                    hs[l] = (1 .- it) .* ht .+ it .* hs[l]
                end
                push!(ys, hs[L])
            end
        else
            bidirect() do xl, i, ti, skip
                Wr,Wi,Wh,Rr,Ri,Rh = w[1+6*(i-1):6i]
                bWr,bWi,bWh,bRr,bRi,bRh = w[6L+1+6*(i-1):6L+6i]
                Wrxt = skip ? xl[:, :, ti] : Wr' * xl[:, :, ti]
                Wixt = skip ? xl[:, :, ti] : Wi' * xl[:, :, ti]
                Whxt = skip ? xl[:, :, ti] : Wh' * xl[:, :, ti]
                rt = sigm.(Wrxt .+ Rr' * hs[i] .+ bWr .+ bRr)
                it = sigm.(Wixt .+ Ri' * hs[i] .+ bWi .+ bRi)
                ht = tanh.(Whxt .+ rt .* (Rh' * hs[i] .+ bRh) .+ bWh)
                hs[i] = (1 .- it) .* ht .+ it .* hs[i]
            end
        end
    else
        error("RNN not supported")
    end
    y = r.direction == 0 ? reshape(hcat(ys...), ysize) : reshape(ys,ysize)
    hyout = hy ? reshape(hcat(hs...), hsize) : nothing
    cyout = cy && r.mode == 2 ? reshape(hcat(cs...), hsize) : nothing
    return (y,hyout,cyout,nothing)
end

# compare sizes ignoring trailing ones
function eqsize(a, b)
    na = length(a)
    nb = length(b)
    (na == nb ? a == b : na > nb ?
     a[1:nb] == b && all(a[nb+1:end] .== 1) :
     b[1:na] == a && all(b[na+1:end] .== 1))
end

# TODO: WIP
function rnntest_bs(batchSizes, r::RNN, w, x,
                    hx=nothing, cx=nothing;
                    # handle=cudnnhandle(), training=false,
                    hy = (hx != nothing),
                    cy = (cx != nothing && r.mode == 2),
                    o...)
    # TODO: fix this implementation
    error("Implementation of batchSizes is not completed in CPU")
    @assert value(r.w) === value(w)
    # Here needs reshaping hidden sizes
    if length(Set{Int}(batchSizes)) == 1
        x_ = reshape(x, (r.inputSize, div(size(x,2), batchSizes[1]), batchSizes[1]))
        y,hy,cy,rs = rnntest(r, w, x_, hx, cx; hy=hy, cy=cy, o...)
        y = reshape(y, (size(y, 1), size(y,2) * size(y,3)))
        return y, hy, cy, rs
    end
    hrem(h, bs1, bs2) = (bs2 < bs1) ? (h[:, 1:bs1-bs2+1, l] for l=1:size(h,3)) : nothing
    hnext(h, bs1, bs2) = (bs2 < bs1) ? h[:, 1:bs2, :] : h
    hrems = []
    ys = []
    crems = []
    ind = 1
    for i = 1:length(batchSizes)
        xt = x[:, ind:ind+batchSizes[i]-1]
        xt = reshape(xt, size(xt)..., 1)
        y, hx, cx = rnntest(r, w, xt, hx, cx;
                            hy=true, cy=true)
        if i > 1
            hr = hrem(hx,batchSizes[i],batchSizes[i+1])
            if hr !== nothing
                hy && push!(hrems, hr...)
                hx = hnext(hx, batchSizes[i],batchSizes[i+1])
            end
            if r.mode == 2
                cr = crem(h,batchSizes[i], batchSizes[i+1])
                if cr !== nothing
                    cy && push!(crems, cr...)
                    cx = hnext(cx, batchSizes[i],batchSizes[i+1])
                end
            end
        end
        ind += batchSizes[i]
        push!(ys, reshape(y, size(y,1,2))) #only last layer output
    end
    # reconstruct the output
    ## hx has size (h,bs[end],l)
    ## hrems has (hs,bsi) dimensional matrices
    hout(hx, hy, hrems) = begin
        nlayers = size(hx, 3)
        if hy
            hts = []
            for i = 1:nlayers
                push!(hts, hy[:,:,i])
            end
            hyouts = []
            for l = 1:nlayers
                push!(hyouts, hcat(hts[l], reverse(hrems[l:nlayers:end-nlayers+l])...))
            end
            hsize = (size(hyouts[end])..., nlayers)
            reshape(length(hyouts) > 1 ? hcat(hyouts...) : hyouts[1], hsize)
        else
            nothing
        end
    end
    return (hcat(ys...), hout(hx, hy, hrems), r.mode==2 ? hout(cx, cy, crems) : nothing, nothing)
end

## DEPRECATED:
function rnninit(x...; o...)
    @warn "rnninit is deprecated, use RNN instead" maxlog=1
    r=RNN(x...; o...)
    return (r,r.w)
end

function rnnparams(r,w;o...)
    @warn "rnnparams(r,w) is deprecated, use rnnparams(r) instead" maxlog=1
    @assert value(w)===value(r.w)
    rnnparams(r;o...)
end

function rnnparam(r,w,l,i,d;o...)
    @warn "rnnparam(r,w,l,i,d) is deprecated, use rnnparam(r,l,i,d) instead" maxlog=1
    @assert value(w)===value(r.w)
    rnnparam(r,l,i,d;o...)
end


text:
# generates a lcov.info file in the format generated by `geninfo`. This format
# can be parsed by a variety of useful utilities to display coverage info

export LCOV

"""
CoverageTools.LCOV Module

This module provides functionality to generate LCOV info format files from
Julia coverage data. It exports the `writefile` function.
"""
module LCOV

using CoverageTools
using CoverageTools: CovCount

export writefile, readfile

"""
    writefile(outfile::AbstractString, fcs)

Write the given coverage data to a file in LCOV info format. The data can either
be a `FileCoverage` instance or a vector of `FileCoverage` instances.
"""
function writefile(outfile::AbstractString, fcs)
    open(outfile, "w") do f
        write(f, fcs)
        nothing
    end
    nothing
end

"""
    write(io::IO, fcs)

Write the given coverage data to an `IO` stream in LCOV info format. The data
can either be a `FileCoverage` instance or a vector of `FileCoverage` instances.
"""
function write end

function write(io::IO, fcs::Vector{FileCoverage})
    for fc in fcs
        write(io, fc)
    end
end

function write(io::IO, fc::FileCoverage)
    instrumented = 0
    covered = 0
    println(io, "SF:$(fc.filename)")
    for (line, cov) in enumerate(fc.coverage)
        (lineinst, linecov) = writeline(io, line, cov)
        instrumented += lineinst
        covered += linecov > 0 ? 1 : 0
    end
    println(io, "LH:$covered")
    println(io, "LF:$instrumented")
    println(io, "end_of_record")
end

# document the writeline function instead of individual methods
"""
    writeline(io::IO, line::Int, count)

Write LCOV data for a single line to the given `IO` stream. Returns a 2-tuple
of the number of lines instrumented (0 or 1) and the count for how many times
the line was executed during testing.
"""
function writeline end

function writeline(io::IO, line::Int, count::Int)
    println(io, "DA:$line,$count")
    (1, count)
end
function writeline(io::IO, line::Int, count::Nothing)
    # skipped line, nothing to do here
    (0, 0)
end

"""
    readfile(infofile::AbstractString) -> Vector{FileCoverage}

Read coverage data from a file in LCOV info format.
"""
function readfile(infofile::AbstractString)
    source_files = FileCoverage[]
    coverage = nothing
    for line in eachline(infofile)
        if startswith(line, "end_of_record")
            coverage = nothing
        elseif (m = match(r"^SF:(.+)", line)) !== nothing
            sf = String(m[1])
            coverage = FileCoverage(sf, "", CovCount[])
            push!(source_files, coverage)
        elseif (m = match(r"^DA:(\d+),(-?\d+)(,[^,\s]+)?", line)) !== nothing
            if coverage !== nothing
                ln = parse(Int, m[1])
                da = parse(Int, m[2])
                cv = coverage.coverage
                lc = length(cv)
                if ln > 0
                    if lc < ln
                        resize!(cv, ln)
                        cv[(lc + 1):ln] .= nothing
                    end
                    cv[ln] = something(cv[ln], 0) + da
                end
            end
        end
        nothing
    end
    return source_files
end

"""
    readfolder(folder) -> Vector{FileCoverage}

Process the contents of a folder of LCOV files to collect coverage statistics.
Will recursively traverse child folders.
Post-process with `merge_coverage_counts(coverages)` to combine duplicates.
"""
function readfolder(folder)
    @info """CoverageTools.LCOV.readfolder: Searching $folder for .info files..."""
    source_files = FileCoverage[]
    files = readdir(folder)
    for file in files
        fullfile = joinpath(folder, file)
        if isfile(fullfile)
            # Is it a tracefile?
            if endswith(fullfile, ".info")
                append!(source_files, readfile(fullfile))
            else
                @debug "CoverageTools.LCOV.readfolder: Skipping $file, not a .info file"
            end
        elseif isdir(fullfile)
            # If it is a folder, recursively traverse
            append!(source_files, readfolder(fullfile))
        end
    end
    return source_files
end

end # module

ext: jl


text:
using ModelingToolkit

# Bounds
@variables u [bounds=(-1,1)]
@test getbounds(u) == (-1, 1)
@test hasbounds(u)

@variables y
@test !hasbounds(y)

# Disturbance
@variables u [disturbance=true]
@test isdisturbance(u)

@variables y
@test !isdisturbance(y)

# Tunable
@parameters u [tunable=true]
@test istunable(u)

@parameters y
@test !istunable(y)

# Distributions
struct FakeNormal end
d = FakeNormal()
@parameters u [dist=d]
@test hasdist(u)
@test getdist(u) == d

@parameters y
@test !hasdist(y)

## System interface
@parameters t
Dₜ = Differential(t)
@variables x(t)=0 u(t)=0 [input=true] y(t)=0 [output=true]
@parameters T [tunable = true, bounds = (0, Inf)]
@parameters k [tunable = true, bounds = (0, Inf)]
@parameters k2
eqs = [
    Dₜ(x) ~ (-k2*x + k*u) / T
    y ~ x
]
sys = ODESystem(eqs, t, name=:tunable_first_order)

p = tunable_parameters(sys)
sp = Set(p)
@test k ∈ sp
@test T ∈ sp
@test k2 ∉ sp
@test length(p) == 2

lb, ub = getbounds(p)
@test lb == [0,0]
@test ub == [Inf, Inf]

b = getbounds(sys)
@test b[T] == (0, Inf)

p = tunable_parameters(sys, default=true)
sp = Set(p)
@test k ∈ sp
@test T ∈ sp
@test k2 ∈ sp
@test length(p) == 3


text:
module HoJDiscordBot

using Discord
using Dates

const ENV["DISCORD_TOKEN"] = read("secrets/BotToken.txt", String)

function main()
    client = Client(ENV["DISCORD_TOKEN"]; prefix='!')

    include("handlers.jl")
    include("commands.jl")

    # --------------- Handlers ---------------
    # The bot listens to all events in the servers and DMs it is in
    # and handles them with a specified function in handlers.jl
    # example: add_hander!(client, Event, function)
    # List of events can be found @ https://posts.xh4h.com/Discord.jl/dev/events/
    # https://posts.xh4h.com/Discord.jl/dev/client/#Discord.add_handler!
    add_handler!(client, MessageCreate, handleping)
    add_handler!(client, MessageCreate, handlepong; tag=:messagepong) # function can be tagged to be deleted later
    delete_handler!(client, MessageCreate, :messagepong) # bot will not reply to "pong"

    # --------------- User Commands ---------------
    # The bot listens to user messages with the specified prefix (!)
    # and runs a function
    # example: add_command!(client, :commandname, function)
    # https://posts.xh4h.com/Discord.jl/dev/client/#Discord.add_command!
    add_command!(client, :ping, replypong; help="Replies pong")
    add_command!(client, :pong, replyping; help="Replies ping")
    add_command!(client, :top, replyworkleaderboard; help="Shows top workers")
    add_command!(client, :work, workcontroller; help="Start work")

    return client
end

end

client = HoJDiscordBot.main()
while true
    @info "Bot starting up"
    open(client)
    wait(client)
    isopen(c) && close(c)
    @info "Lost connection to Discord, restarting bot."
end

max_stars_repo_name: dzonimn/HoJDiscordBot

max_issues_repo_name: dzonimn/HoJDiscordBot

max_forks_repo_name: dzonimn/HoJDiscordBot


text:
export dmrg

function dmrg(H::MPO,
              psi0::MPS,
              sweeps::Sweeps;
              kwargs...)::Tuple{Float64,MPS}

  which_factorization::String = get(kwargs,:which_factorization,"automatic")

  psi = copy(psi0)
  N = length(psi)

  PH = ProjMPO(H)
  position!(PH,psi0,1)
  energy = 0.0

  for sw=1:nsweep(sweeps)
    sw_time = @elapsed begin

    @debug begin
      reset!(timer)
      println("\n\n~~~ This is sweep $sw ~~~\n\n")
    end

    for (b,ha) in sweepnext(N)
      position!(PH,psi,b)

      phi = psi[b]*psi[b+1]

      energy,phi = davidson(PH,phi;kwargs...)

      dir = ha==1 ? "fromleft" : "fromright"
      replaceBond!(psi,b,phi;
                   maxdim=maxdim(sweeps,sw),
                   mindim=mindim(sweeps,sw),
                   cutoff=cutoff(sweeps,sw),
                   dir=dir,
                   which_factorization=which_factorization)
    end
    end
    @printf "After sweep %d energy=%.12f maxDim=%d time=%.3f\n" sw energy maxDim(psi) sw_time
    @debug printTimes(timer)
  end
  return (energy,psi)
end

@doc """
Use the density matrix renormalization group (DMRG) algorithm to
optimize a matrix product state (MPS) to be the eigenvector
of the Hermitian matrix product operator (MPO) H with minimal
eigenvalue.

Inputs:
* `H::MPO` - a Hermitian MPO
* `psi0::MPS` - MPS used to initialize the optimization
* `sweeps::Sweeps` - `Sweeps` object used to control the algorithm

Returns:
* `energy::Float64` - eigenvalue of the optimized MPS
* `psi::MPS` - optimized MPS

""" dmrg


text:
"""
    module IrosBase

Essential tools used in the rest of Iros.
"""
module Common

include("Points.jl")
using .Points
export Point

include("Shapes.jl")
using .Shapes
export Interval
export Circle
export Ellipse
export DeformedDisk
export Annulus
export Square
export Rectangle

include("BoundaryLayers.jl")
using .BoundaryLayers
export PML
export cPML
export noBL

include("BoundaryConditions.jl")
using .BoundaryConditions
export noBC
export DirichletBC
export NeumannBC
export FloquetBC
export MatchedBC

include("Boundaries.jl")
using .Boundaries
export Boundary

include("DielectricFunctions.jl")
using .DielectricFunctions
export DielectricFunction
export PumpFunction

include("Dispersions.jl")
using .Dispersions
export TwoLevelSystem
export jacobian_lasing
export susceptability

include("Lattices.jl")
using .Lattices
export Lattices
export Lattice

include("Domains.jl")
using .Domains
export Domain
export Cavity
export Resonator
export Void
export Dielectric
export Waveguide
export Lead

include("Curlcurls.jl")
using .Curlcurls
export Curlcurl

include("SelfEnergies.jl")
using .SelfEnergies

include("Simulations.jl")
using .Simulations
export Simulation
export smooth!
export update_dielectric!
export update_pump!
export Unsymmetric
export Symmetric
export Hermitian

include("MaxwellOperators.jl")
using .MaxwellOperators
export maxwell
export maxwell_lep
export maxwell_nep

end # module


text:
using RobotDynamics
using StaticArrays
using LinearAlgebra
using BenchmarkTools
using Test

## Setup
include("random_linear.jl")

N = 11
n,m = 15,10
dt = 0.1
times = range(0,step=dt,length=N)
t = times[3]
Ad,Bd = gencontrollable(n,m)
@test isstabled(Ad)
x = @SVector rand(n)
u = @SVector rand(m)
z = KnotPoint(x,u,dt,t)

## Discrete Time Invariant
struct DLTI{TA,TB} <: DiscreteLTI
    A::TA
    B::TB
end
RobotDynamics.state_dim(::DLTI) = 15
RobotDynamics.control_dim(::DLTI) = 10

model = DLTI(SMatrix{n,n}(Ad),SMatrix{n,m}(Bd))
@test RobotDynamics.is_affine(model) == Val(false)
@test RobotDynamics.is_time_varying(model) == false

@test_throws ErrorException discrete_dynamics(DiscreteSystemQuadrature, model, x, u, t, dt)
RobotDynamics.get_A(model::DLTI) = model.A
RobotDynamics.get_B(model::DLTI) = model.B
x2 = discrete_dynamics(DiscreteSystemQuadrature, model, x, u, t, dt)
@test x2 ≈ Ad*x + Bd*u
@test (@ballocated discrete_dynamics(DiscreteSystemQuadrature, $model, $x, $u, $t, $dt) samples=2 evals=2) == 0

F = RobotDynamics.DynamicsJacobian(model)
discrete_jacobian!(DiscreteSystemQuadrature, F, model, z)
@test F ≈ [Ad Bd]

## DLTV
struct DLTV{TA,TB} <: DiscreteLTV
    A::Vector{TA}
    B::Vector{TB}
    times::Vector{Float64}
end
RobotDynamics.state_dim(::DLTV) = 15
RobotDynamics.control_dim(::DLTV) = 10
RobotDynamics.get_A(model::DLTV, k::Int) = model.A[k]
RobotDynamics.get_B(model::DLTV, k::Int) = model.B[k]

N = 11
systems = [gencontrollable(n,m) for k = 1:N]
A = [SMatrix{n,n}(s[1]) for s in systems]
B = [SMatrix{n,m}(s[2]) for s in systems]

model = DLTV(A,B,collect(times))

@test RobotDynamics.is_affine(model) == Val(false)
@test RobotDynamics.is_time_varying(model) == true

@test_throws ErrorException discrete_dynamics(DiscreteSystemQuadrature, model, x, u, t, dt)
RobotDynamics.get_times(model::DLTV) = model.times
x2 = discrete_dynamics(DiscreteSystemQuadrature, model, x, u, t, dt)
@test x2 ≈ A[3]*x + B[3]*u
@test !(x2 ≈ A[1]*x + B[1]*u)
@test (@ballocated discrete_dynamics(DiscreteSystemQuadrature, $model, $x, $u, $t, $dt)) == 0

F = RobotDynamics.DynamicsJacobian(model)
discrete_jacobian!(DiscreteSystemQuadrature, F, model, z)
@test F ≈ [A[3] B[3]]
@test !(F ≈ [A[1] B[1]])

## CLTI
struct CLTI{TA,TB} <: ContinuousLTI
    A::TA
    B::TB
end
RobotDynamics.state_dim(::CLTI) = 15
RobotDynamics.control_dim(::CLTI) = 10

model = CLTI(SMatrix{n,n}(Ad),SMatrix{n,m}(Bd))
@test_throws ErrorException dynamics(model, x, u, t)
RobotDynamics.get_A(model::CLTI) = model.A
RobotDynamics.get_B(model::CLTI) = model.B
xdot = dynamics(model, x, u, t)
@test xdot ≈ Ad*x + Bd*u
@test (@ballocated dynamics($model, $x, $u, $t)) == 0

jacobian!(F, model, z)
@test F ≈ [Ad Bd]

x2 = discrete_dynamics(RK3, model, z)
begin
    local k1_ = (Ad*x + Bd*u)*dt
    local k2_ = (Ad*(x + k1_/2) + Bd*u)*dt
    local k3_ = (Ad*(x - k1_ + 2k2_) + Bd*u)*dt
    @test x2 ≈ x + (k1_+4k2_+k3_)/6
end

discrete_jacobian!(RK3, F, model, z)
@test F.A ≈ I(n) + (Ad*dt + 4*Ad*(I + Ad*dt/2)*dt + Ad*(I - Ad*dt + 2*Ad*(I + Ad*dt/2)*dt)*dt)/6
@test F.B ≈ (Bd*dt + 4*(Ad'Bd*dt/2 + Bd)*dt + (Ad'*(-Bd*dt + 2*(Ad'Bd*dt/2 + Bd)*dt) + Bd)*dt)/6
@test (@ballocated discrete_jacobian!(RK3, $F, $model, $z)) == 0

## CLTV
struct CLTV{TA,TB} <: ContinuousLTV
    A::Vector{TA}
    B::Vector{TB}
    times::Vector{Float64}
end
RobotDynamics.state_dim(::CLTV) = 15
RobotDynamics.control_dim(::CLTV) = 10
RobotDynamics.get_A(model::CLTV, k::Int) = model.A[k]
RobotDynamics.get_B(model::CLTV, k::Int) = model.B[k]

model = CLTV(A,B,collect(times))
@test is_time_varying(model)
@test is_affine(model) == Val(false)

@test_throws ErrorException dynamics(model, x, u, t)
RobotDynamics.get_times(model::CLTV) = model.times
xdot = dynamics(model, x, u, t)
@test xdot ≈ dynamics(model, z)
@test xdot ≈ A[3]*x + B[3]*u

jacobian!(F, model, z)
@test F ≈ [A[3] B[3]]

x2 = discrete_dynamics(RK2, model, z)
A_,B_ = A[3], B[3]
begin
    local k1_ = (A_*x + B_*u)*dt
    local k2_ = (A_*(x + k1_/2) + B_*u)*dt
    @test x2 ≈ x + k2_
end

discrete_jacobian!(RK2, F, model, z)
@test F.A ≈ I(n) + (A_'*(I + A_*dt/2))*dt
@test F.B ≈ (A_'B_*dt/2 + B_)*dt

## Test an affine system
struct DLTI2{TA,TB,TD} <: DiscreteLTI
    A::TA
    B::TB
    d::TD
end
RobotDynamics.state_dim(::DLTI2) = 15
RobotDynamics.control_dim(::DLTI2) = 10
RobotDynamics.get_A(model::DLTI2) = model.A
RobotDynamics.get_B(model::DLTI2) = model.B
RobotDynamics.get_d(model::DLTI2) = model.d
RobotDynamics.is_affine(::DLTI2) = Val(true)

d = @SVector rand(n)
model = DLTI2(SMatrix{n,n}(Ad), SMatrix{n,m}(Bd), d)
x2 = discrete_dynamics(DiscreteSystemQuadrature, model, z)
@test x2 ≈ Ad*x + Bd*u + d
@test (@ballocated discrete_dynamics(DiscreteSystemQuadrature, $model, $z)) == 0

F = RobotDynamics.DynamicsJacobian(model)
discrete_jacobian!(DiscreteSystemQuadrature, F, model, z)
@test F ≈ [Ad Bd]


text:
using Pkg

Pkg.add("CUDA")
Pkg.add("Knet")
Pkg.add("AutoGrad")
Pkg.add("Images")
Pkg.add("ImageMagick")
Pkg.add("ImageDraw")
Pkg.add("MAT")
Pkg.add("JLD2")
Pkg.add("JLD")
Pkg.add("FileIO")
Pkg.add("PyCall")
Pkg.add("IJulia")
Pkg.add("ArgParse")
Pkg.add("Plots")

Pkg.precompile()

max_stars_repo_name: barisbatuhan/retinaface

max_issues_repo_name: barisbatuhan/retinaface

max_forks_repo_name: barisbatuhan/retinaface


text:
@pyimport wntr.network.controls as controls
@pyimport wntr.sim.epanet as sims
include("UpdateExtrema.jl")
include("ManipulateTanks.jl")
include("ManipulatePumps.jl")
function Q_extrema(wn_dict::Dict{Any,Any}, Q_lb::Float64, tight_coef::Float64)
    Qmax_dict = Dict{String, Float64}()
    Qmin_dict = Dict{String, Float64}()
    Hmin_dict = Dict{String,Float64}()
    Hmax_dict= Dict{String,Float64}()
    wn = wn_dict["wn"]
    controls_dict = wn_dict["controls"]
    control_names = keys(controls_dict)
    duration = wn["options"]["time"]["duration"]
    time_periods = wn["options"]["time"]["report_timestep"]
    time_periods = convert(Float64, time_periods)
    time_steps = 24 #duration/time_periods
    #time_steps = convert(Int64, time_steps)
    ############ Level Limits Only ##################################
    tank_level(wn, wn_dict, Qmax_dict, Qmin_dict, Hmax_dict, Hmin_dict, tight_coef, Q_lb)
    ############ Tank Controls and Level Limits #####################
    remove_controls(wn)

    tank_triggers(wn, wn_dict, controls_dict, control_names, Qmax_dict, Qmin_dict, Hmax_dict, Hmin_dict, tight_coef, Q_lb)
    ################ Tank Controls only #############################
    for (key, tank) in wn_dict["tanks"]
        wn[:get_node](key)[:init_level] = tank["init_level"]
    end
    UpdateExtrema(wn, wn_dict, Qmax_dict, Qmin_dict, Hmax_dict, Hmin_dict, 2, 2, tight_coef, Q_lb)
    ############## Restore wn ########################################
    remove_controls(wn)
    for (name, control) in controls_dict
        wn[:add_control](name, control)
    end
    ############## Random Level Limits Only ##########################

    tank_level_rand(wn, wn_dict, Qmax_dict, Qmin_dict, Hmax_dict, Hmin_dict, tight_coef, Q_lb)

    ############## Random Tank Controls and Random Level Limits ######
    remove_controls(wn)
    tank_triggers_rand(wn, wn_dict, controls_dict, control_names, Qmax_dict, Qmin_dict, Hmax_dict, Hmin_dict, tight_coef, Q_lb)

    ############# Random Tank Controls Only ##########################
    for (key, tank) in wn_dict["tanks"]
        wn[:get_node](key)[:init_level] = tank["init_level"]
    end
    UpdateExtrema(wn, wn_dict, Qmax_dict, Qmin_dict, Hmax_dict, Hmin_dict, 2, 2, tight_coef, Q_lb)
    ############## Restore wn ########################################
    remove_controls(wn)
    for (name, control) in controls_dict
        wn[:add_control](name, control)
    end
    ################### Manipulate Pumps ##############################
    ################### Pumps all on after TN only ####################
    remove_controls(wn_dict, wn) #only removes pump controls
    for Tn = 1:length(time_steps)
        OnAfterTn(wn_dict, wn, Tn, Qmax_dict, Qmin_dict, Hmax_dict, Hmin_dict, tight_coef, Q_lb)
    end
    ############## Restore wn ########################################
    remove_controls(wn)
    for (name, control) in controls_dict
        wn[:add_control](name, control)
    end
    ##################Pumps all Off after TN #########################
    remove_controls(wn_dict, wn) #only removes pump controls
    for Tn = 1:length(time_steps)
        OffAfterTn(wn_dict, wn, Tn, Qmax_dict, Qmin_dict, Hmax_dict, Hmin_dict, tight_coef, Q_lb)
    end
    ############## Restore wn ########################################
    remove_controls(wn)
    for (name, control) in controls_dict
        wn[:add_control](name, control)
    end
    ################### Pumps on and Off at Random Times #############
    remove_controls(wn_dict, wn) #only removes pump controls
    RandomTime(wn_dict, wn, time_steps, time_periods, Qmax_dict, Qmin_dict, Hmax_dict, Hmin_dict, tight_coef, Q_lb)
    ######### Random Initial Levels and Pump On/Off randomly #########
    tank_level_rand(wn, wn_dict, Qmax_dict, Qmin_dict, Hmax_dict, Hmin_dict, tight_coef, Q_lb)
    ############## Restore wn ########################################
    remove_controls(wn)
    for (name, control) in controls_dict
        wn[:add_control](name, control)
    end
    UpdateExtrema(wn, wn_dict, Qmax_dict, Qmin_dict, Hmax_dict, Hmin_dict, 2, 2, tight_coef, Q_lb) #test original system
    ############ Finalize Qmax/Qmin ##################################
    num_tanks = length(wn_dict["tanks"])
    num_sim = 2*num_tanks
    UpdateExtrema(wn, wn_dict, Qmax_dict, Qmin_dict, Hmax_dict, Hmin_dict, num_sim, num_sim +1, tight_coef, Q_lb)

    return Qmin_dict, Qmax_dict, Hmin_dict, Hmax_dict

end

function remove_controls(wn::PyCall.PyObject)
    for name in wn[:control_name_list]
        wn[:remove_control](name)
    end
end

function remove_controls(wn_dict::Dict{Any,Any}, wn::PyCall.PyObject)
    controls_dict = wn_dict["controls"]
    for name in wn[:control_name_list]
        target_obj = controls_dict[name][:_then_actions][1][:_target_obj]
        if target_obj[:link_type] == "Pump"
            wn[:remove_control](name)
        end
    end
end


text:
# Copyright (c) 2018: Matthew Wilhelm & Matthew Stuber.
# This code is licensed under MIT license (see LICENSE.md for full details)
#############################################################################
# EAGO
# A development environment for robust and global optimization
# See https://github.com/PSORLab/EAGO.jl
#############################################################################
# src/eago_optimizer/unsafe_utilities.jl
# Inbounds and non-allocating versions of simple Julia functions.
#############################################################################

"""
$(TYPEDSIGNATURES)

Performs `map!(f, y, x)` in an unsafe manner if y[i] is true, else no-op.
Assumes `n == length(x) == length(y)`. About 2x faster for small arrays (n < 1000).
"""
function unsafe_check_fill!(f, y::Vector{T}, x::T, n::Int) where T
    i = 1
    m = n + 1
    while i < m
        flag = f(@inbounds y[i])
        if flag
            @inbounds y[i] = x
        end
        i += 1
    end
    nothing
end

_rf_findmax((fm, m), (fx, x)) = Base.isless(fm, fx) ? (fx, x) : (fm, m)
map_findmax(f, itr) = mapfoldl(((k, v),) -> (f(v), k), _rf_findmax, pairs(itr))
map_argmax(f, itr) = map_findmax(f, itr)[2]

argmax(f, domain) = mapfoldl(x -> (f(x), x), _rf_findmax, domain)[2]

relative_gap(L::Float64, U::Float64) = ((L > -Inf) && (U < Inf)) ?  abs(U - L)/(max(abs(L), abs(U))) : Inf
relative_tolerance(L::Float64, U::Float64, tol::Float64) = relative_gap(L, U)  > tol || ~(L > -Inf)

#=
"""
    gershgorin_λmin

Computes a lower bound on the smallest eigenvalue of `x` by means of Gershgorin's
Circle Theorem. See the link provided for information
https://mathworld.wolfram.com/GershgorinCircleTheorem.html.
"""
function gershgorin_λmin(x::AbstractMatrix{T}) where T
    xsum = -sum(i -> abs(@inbounds x[i]), diagind(x))
    xmin = typemax(T)
    @inbounds for k in diagind(x)
        xk = x[k]
        xadd = abs(xk) + xk
        xsum += xadd
        if xsum < xmin
            xmin = xsum
        end
        xsum -= xadd
    end
    return xmin
end
=#

max_stars_repo_name: MatthewStuber/EAGO

max_issues_repo_name: MatthewStuber/EAGO

max_forks_repo_name: MatthewStuber/EAGO


text:
using DrWatson
@quickactivate
using Test
using GenerativeAD
using UCI
using ArgParse

s = ArgParseSettings()
@add_arg_table! s begin
    "--complete"
    	action = :store_true
        help = "run a more thorough set of tests, may be slow"
end
parsed_args = parse_args(ARGS, s)
@unpack complete = parsed_args

include("data.jl")
include("experiments.jl")
if complete
	include("models/runtests.jl")
end


text:
"""
    Econometrics

    Econometrics in Julia.
"""
module Econometrics
using Base.Iterators: flatten
using CategoricalArrays: CategoricalArrays, CategoricalValue, levels, levels!, isordered, ordered!
using Distributions:
    Distributions,
    cdf,
    ccdf,
    Categorical,
    Chisq,
    FDist,
    logpdf,
    Logistic,
    Normal,
    pdf,
    TDist,
# Statistics
    mean,
    quantile,
    var
using FillArrays: FillArrays, Ones
using ForwardDiff: ForwardDiff, Dual, value
using LinearAlgebra:
    LinearAlgebra,
    bunchkaufman,
    bunchkaufman!,
    cholesky!,
    diag,
    diagm,
    Diagonal,
    Hermitian,
    I,
    LowerTriangular,
    qr,
    UpperTriangular
using MLJModelInterface: MLJModelInterface, MLJModelInterface, Probabilistic
using Optim: Optim, hessian!, optimize, minimizer, TwiceDifferentiable
using Parameters: Parameters, @unpack, @pack!
using Printf: Printf, @sprintf
using StatsBase:
    StatsBase, aic, aicc, bic, harmmean, FrequencyWeights, CoefTable, ConvergenceException, Weights
using StatsFuns: StatsFuns, softmax
using StatsModels:
    StatsModels,
    AbstractContrasts,
    AbstractTerm,
    apply_schema,
    CategoricalTerm,
    ConstantTerm,
    ContrastsMatrix,
    DummyCoding,
    @formula,
    FormulaTerm,
    FunctionTerm,
    InteractionTerm,
    InterceptTerm,
    MatrixTerm,
    modelcols,
    schema,
    terms,
    termvars
using TableOperations:
    TableOperations,
    select,
# Tables
    Tables,
    Tables.columntable,
    Tables.columns,
    Tables.eachcolumn,
    Tables.materializer,
    Tables.rows
import Base: show
import StatsBase:
    coef,
    coefnames,
    coeftable,
    confint,
    deviance,
    islinear,
    nulldeviance,
    loglikelihood,
    nullloglikelihood,
    score,
    nobs,
    dof,
    dof_residual,
    mss,
    rss,
    informationmatrix,
    vcov,
    stderror,
    weights,
    isfitted,
    fit,
    fit!,
    r2,
    adjr2,
    fitted,
    response,
    meanresponse,
    modelmatrix,
    leverage,
    residuals,
    predict,
    predict!,
    dof_residual,
    RegressionModel,
    params
import StatsModels: hasintercept, implicit_intercept
import MLJModelInterface: clean!

"""
    const DATAPATH::String = realpath(joinpath(pkgdir(Econometrics), "data"))
Return the path to the data directory for the Econometrics.jl module.
"""
const DATAPATH = realpath(joinpath(dirname(@__FILE__), "..", "data"))
foreach(
    file -> include(joinpath(dirname(@__DIR__), "src", "$file.jl")),
    ["structs", "transformations", "formula", "main", "mlj-model-interface", "solvers", "statsbase", "wald"],
    )
export @formula,
    DummyCoding,
    aic,
    aicc,
    bic,
    coef,
    coefnames,
    coeftable,
    confint,
    deviance,
    islinear,
    nulldeviance,
    loglikelihood,
    nullloglikelihood,
    nobs,
    dof,
    mss,
    rss,
    informationmatrix,
    vcov,
    stderror,
    weights,
    isfitted,
    fit,
    fit!,
    r2,
    adjr2,
    fitted,
    response,
    meanresponse,
    modelmatrix,
    leverage,
    residuals,
    predict,
    dof_residual,
    params,
    hasintercept,
    EconometricModel,
    EconometricsMLJModel,
    absorb,
    BetweenEstimator,
    RandomEffectsEstimator,
    ContinuousResponse,
    OIM,
    HC0,
    HC1,
    HC2,
    HC3,
    HC4,
    levels!
end


text:
struct LinearCombination{T, As<:LinearMapTuple} <: LinearMap{T}
    maps::As
    function LinearCombination{T, As}(maps::As) where {T, As}
        N = length(maps)
        ax = axes(maps[1])
        for n in eachindex(maps)
            A = maps[n]
            axes(A) == ax || throw(DimensionMismatch("LinearCombination"))
            @assert promote_type(T, eltype(A)) == T  "eltype $(eltype(A)) cannot be promoted to $T in LinearCombination constructor"
        end
        new{T, As}(maps)
    end
end

LinearCombination{T}(maps::As) where {T, As} = LinearCombination{T, As}(maps)

MulStyle(A::LinearCombination) = MulStyle(A.maps...)

# basic methods
Base.size(A::LinearCombination) = size(A.maps[1])
Base.axes(A::LinearCombination) = axes(A.maps[1])
# following conditions are sufficient but not necessary
LinearAlgebra.issymmetric(A::LinearCombination) = all(issymmetric, A.maps)
LinearAlgebra.ishermitian(A::LinearCombination) = all(ishermitian, A.maps)
LinearAlgebra.isposdef(A::LinearCombination) = all(isposdef, A.maps)

# adding linear maps
"""
    +(A::LinearMap, B::LinearMap)::LinearCombination

Construct a (lazy) representation of the sum/linear combination of the two operators.
Sums of `LinearMap`/`LinearCombination` objects and
`LinearMap`/`LinearCombination` objects are reduced to a single `LinearCombination`.
In sums of `LinearMap`s and `AbstractMatrix`/`UniformScaling` objects, the latter
get promoted to `LinearMap`s automatically.

# Examples
```jldoctest; setup=(using LinearAlgebra, LinearMaps)
julia> CS = LinearMap{Int}(cumsum, 3)::LinearMaps.FunctionMap;

julia> LinearMap(ones(Int, 3, 3)) + CS + I + rand(3, 3);
```
"""
function Base.:(+)(A₁::LinearMap, A₂::LinearMap)
    size(A₁) == size(A₂) || throw(DimensionMismatch("+"))
    T = promote_type(eltype(A₁), eltype(A₂))
    return LinearCombination{T}(tuple(A₁, A₂))
end
function Base.:(+)(A₁::LinearMap, A₂::LinearCombination)
    size(A₁) == size(A₂) || throw(DimensionMismatch("+"))
    T = promote_type(eltype(A₁), eltype(A₂))
    return LinearCombination{T}(tuple(A₁, A₂.maps...))
end
Base.:(+)(A₁::LinearCombination, A₂::LinearMap) = +(A₂, A₁)
function Base.:(+)(A₁::LinearCombination, A₂::LinearCombination)
    size(A₁) == size(A₂) || throw(DimensionMismatch("+"))
    T = promote_type(eltype(A₁), eltype(A₂))
    return LinearCombination{T}(tuple(A₁.maps..., A₂.maps...))
end
Base.:(-)(A₁::LinearMap, A₂::LinearMap) = +(A₁, -A₂)

# comparison of LinearCombination objects, sufficient but not necessary
Base.:(==)(A::LinearCombination, B::LinearCombination) =
    (eltype(A) == eltype(B) && A.maps == B.maps)

# special transposition behavior
LinearAlgebra.transpose(A::LinearCombination) =
    LinearCombination{eltype(A)}(map(transpose, A.maps))
LinearAlgebra.adjoint(A::LinearCombination) =
    LinearCombination{eltype(A)}(map(adjoint, A.maps))

# multiplication with vectors & matrices
for (In, Out) in ((AbstractVector, AbstractVecOrMat), (AbstractMatrix, AbstractMatrix))
    @eval begin
        function _unsafe_mul!(y::$Out, A::LinearCombination, x::$In)
            _unsafe_mul!(y, first(A.maps), x)
            _mul!(MulStyle(A), y, A, x, true)
            return y
        end
        function _unsafe_mul!(y::$Out, A::LinearCombination, x::$In, α::Number, β::Number)
            if iszero(α) # trivial cases
                iszero(β) && return fill!(y, zero(eltype(y)))
                isone(β) && return y
                return rmul!(y, β)
            else
                A1 = first(A.maps)
                if MulStyle(A1) === ThreeArg() && !iszero(β)
                    # if we need an intermediate vector, allocate here and reuse in
                    # LinearCombination multiplication
                    !isone(β) && rmul!(y, β)
                    z = similar(y)
                    muladd!(ThreeArg(), y, A1, x, α, z)
                    __mul!(y, Base.tail(A.maps), x, α, z)
                else # MulStyle(A1) === FiveArg() || β == 0
                    # this is allocation-free
                    _unsafe_mul!(y, A1, x, α, β)
                    # let _mul! decide whether an intermediate vector needs to be allocated
                    _mul!(MulStyle(A), y, A, x, α)
                end
                return y
            end
        end
    end
end

function _mul!(::FiveArg, y, A::LinearCombination, x, α)
    return __mul!(y, Base.tail(A.maps), x, α, nothing)
end
function _mul!(::ThreeArg, y, A::LinearCombination, x, α)
    return __mul!(y, Base.tail(A.maps), x, α, similar(y))
end

__mul!(y, As::LinearMapTuple, x, α, z) =
    __mul!(__mul!(y, first(As), x, α, z), Base.tail(As), x, α, z)
__mul!(y, A::Tuple{LinearMap}, x, α, z) = __mul!(y, first(A), x, α, z)
__mul!(y, A::LinearMap, x, α, z) = muladd!(MulStyle(A), y, A, x, α, z)

muladd!(::FiveArg, y, A, x, α, _) = _unsafe_mul!(y, A, x, α, true)
function muladd!(::ThreeArg, y, A, x, α, z)
    _unsafe_mul!(z, A, x)
    if isone(α)
        y .+= z
    else
        y .+= z .* α
    end
    return y
end


text:
import Base.REPL: REPLDisplay
export Query, ask, form, answers, make_widget, on_done

immutable Query{T}
    w::Widget
end

function (::Type{Query{Int}})()
    t = TextInput(IOBuffer(copy(b"1234"),true,true))
    t.optheight = 1
    Query{Int}(t)
end

function on_done(f::Function,i::Query{Int})
    i.w.on_done = str->f(parse(Int,str))
end

function (::Type{Query{AbstractString}})()
    t = TextInput()
    t.optheight = 1
    Query{AbstractString}(t)
end
on_done(f::Function,i::Query{AbstractString}) = i.w.on_done = f

function display(d::REPLDisplay, w::Widget)
    wait(InlineDialog(w,d.repl.t))
end

immutable Question
    text::AbstractString
    input::Query
end

function ask(question, answerT)
    Question(question,Query{answerT}())
end

make_widget(i::Query) = i.w
function make_widget(q::Question)
    cd = TerminalUI.CellDisplay(
        Vector{VT100.Cell}[[
            VT100.Cell(TerminalUI.dummycell(:default,:default,0),content = '['),
            VT100.Cell(TerminalUI.dummycell(:green,:default,0),content = '?'),
            VT100.Cell(TerminalUI.dummycell(:default,:default,0),content = ']')]])

    cd2 = TerminalUI.CellDisplay(
        Vector{VT100.Cell}[[
        VT100.Cell(TerminalUI.dummycell(:default,:default,0),content = c)
        for c in q.text]])

    RowLayout([cd,cd2,q.input.w],[1,1,1],[4,strwidth(q.text)+1,0],3)
end

display(d::REPLDisplay, q::Question) = display(d,make_widget(q))

immutable Form
    questions
    chain::ScrollableChain
    keyidxs
    answers
end

function form(questions; style = :incremental)
    widgets = map(x->make_widget(x[2]),questions)
    if style == :incremental
        for i = 2:length(widgets)
            hide(widgets[i])
        end
    end
    f = Form(questions,ScrollableChain(widgets),
        Dict(k => i for (i,(k,_)) in enumerate(questions)),
        Dict{Symbol,Any}())
    for (k,q) in questions
        on_done(q.input) do v
            answer(f, k, v)
            focus_next(f.chain, show_invisible = true) || (invalidate(f.chain); return :done)
        end
    end
    f
end

display(d::REPLDisplay,f::Form) = display(d,f.chain)
FullScreenDialog(f::Form,tty) = FullScreenDialog(f.chain, tty)

function answer(form::Form, key, value)
    # First record the answer
    form.answers[key] = value
    # Create a widget to permanently display the answer
    str = sprint(print,value)
    cd = TerminalUI.CellDisplay(
        Vector{VT100.Cell}[[
        VT100.Cell(TerminalUI.dummycell(:cyan,:default,Bright),content = c)
        for c in str]])
    # Replace the question widget with an answer widget
    form.chain.children[form.keyidxs[key]].cols[3] = cd
end

answers(f::Form) = f.answers

make_widget(w::Widget) = w
make_widget(s::AbstractString) = Label(s)
function make_widget(a::Array)
    if ndims(a) == 1
        ScrollableChain(map(make_widget,a))
    elseif ndims(a) == 2
        if size(a)[1] > 1
            widgets = [make_widget(a[:,i]) for i = 1:size(a)[2]]
        else
            widgets = map(make_widget,vec(a))
        end
        RowLayout(widgets,[1 for _ in widgets],[0 for _ in widgets],length(widgets))
    else
        error("Higher dimensional widgets only supported for
        higher dimensional terminals")
    end
end


text:
"""
    distance(from::Position, to::Position, units::String="kilometers")

Calculate the distance between two Points or Positions in degrees, radians, miles, or kilometers.
The distance is calculated using the [Haversine formula](http://en.wikipedia.org/wiki/Haversine_formula) to account for global curvature.

# Examples
```jldoctest
julia> point1 = Point([35, 45])
Point([35.0, 45.0])

julia> point2 = Point([38, 43])
Point([38.0, 43.0])

julia> distance(point1, point2)
327.1241022693492

julia> distance(point1, point2, "miles")
203.26549343667307

julia> distance(point1, point2, "degrees")
2.9384603841845878
```
"""
function distance(from::Position, to::Position, units::String="kilometers")
    δlat = deg2rad((ycoord(to) - ycoord(from)))
    δlon = deg2rad((xcoord(to) - xcoord(from)))

    latFrom = deg2rad(ycoord(from))
    latTo = deg2rad(ycoord(to))

    res = sin(δlat / 2)^2 + sin(δlon / 2)^2 * cos(latFrom) * cos(latTo)

    return radians_to_length(2 * atan(sqrt(res), sqrt(1 - res)), units)
end

distance(from::Point, to::Point, units::String="kilometers") = distance(from.coordinates, to.coordinates, units)

"""
    rhumb_distance(from::Position, to::Position, units::String="kilometers")

Calculate the distance along a rhumb line between two Points or Positions in degrees, radians,
miles, or kilometers.

# Examples
```jldoctest
julia> point1 = Point([35, 45])
Point([35.0, 45.0])

julia> point2 = Point([38, 43])
Point([38.0, 43.0])

julia> rhumb_distance(point1, point2)
327.1421400706353

julia> rhumb_distance(point1, point2, "miles")
203.2767016067636

julia> rhumb_distance(point1, point2, "degrees")
2.9386224124916716
```
"""
function rhumb_distance(from::Position, to::Position, units::String="kilometers")
    toLon = xcoord(to)
    toLon += (xcoord(to) - xcoord(from) > 180) ? -360 : ((xcoord(from) - xcoord(to) > 180)) ? 360 : 0

    ϕ1 = ycoord(from) * pi / 180
    ϕ2 = ycoord(to) * pi / 180

    Δϕ = ϕ2 - ϕ1

    Δλ = abs((toLon - xcoord(from))) * pi / 180
    if Δλ > pi
        Δλ -= 2 * pi
    end

    Δψ = log(tan(ϕ2 / 2 + pi / 4) / tan(ϕ1 / 2 + pi / 4))
    q = abs(Δψ) > 10e-12 ?  Δϕ / Δψ : cos(ϕ1)

    Δ = sqrt(Δϕ * Δϕ + q * q * Δλ * Δλ)

    res =  Δ * earth_radius

    return convert_length(res, "metres", units)
end

rhumb_distance(from::Point, to::Point, units::String="kilometers") = rhumb_distance(from.coordinates, to.coordinates, units)

"""
    point_to_line_distance(point::Point, line::LineString, units::String="kilometers", method::String="geodesic")

Return the minimum distance between a Point and a LineString, being the distance from a line the
minimum distance between the point and any segment of the `LineString`.

# Examples
```jldoctest
julia> point = Point([35, 45])
Point([35.0, 45.0])

julia> line = LineString([[35, 42], [38, 42.5], [40, 43]])
LineString(Array{Float64,1}[[35.0, 42.0], [38.0, 42.5], [40.0, 43.0]])

julia> point_to_line_distance(point, line)
326.92768049303146

julia> point_to_line_distance(point, line, "degrees", "planar")
2.936700172184462

julia> point_to_line_distance(point, line, "degrees")
2.9366959846668
```
"""
function point_to_line_distance(point::Point, line::LineString, units::String="kilometers", method::String="geodesic")
    !(method in ["planar", "geodesic"]) && throw(error("$(method) is not a valid options. Please use either 'planar' or 'geodesic'."))

    dist = Inf

    for i in 1:length(line.coordinates) - 1
        segment = LineString([line.coordinates[i], line.coordinates[i+1]])

        a = segment.coordinates[1]
        b = segment.coordinates[2]
        d = distance_to_segment(point, Point(a), Point(b), units, method)
        d < dist && (dist = d)
    end
    return convert_length(dist, "degrees", units)
end

"""
    distance_to_segment(point::Point, first::Point, last::Point, units::String="degrees", method::String="planar")

Return the distance between a point P on a segment AB.
"""
function distance_to_segment(point::Point, first::Point, last::Point, units::String="degrees", method::String="geodesic")
    !(method in ["planar", "geodesic"]) && throw(error("$(method) is not a valid options. Please use either 'planar' or 'geodesic'."))

    v = [last.coordinates[1] - first.coordinates[1], last.coordinates[2] - first.coordinates[2]]
    w = [point.coordinates[1] - first.coordinates[1], point.coordinates[2] - first.coordinates[2]]

    c1 = sum(w.*v)
    if c1 <= 0
        return method === "planar" ? rhumb_distance(point.coordinates, first.coordinates, units) : distance(point.coordinates, first.coordinates, "degrees")
    end
    c2 = sum(v.*v)

    if c2 <= c1
        return method === "planar" ? rhumb_distance(point.coordinates, last.coordinates, units) : distance(point.coordinates, last.coordinates, "degrees")
    end

    b2 = c1 / c2

    p = [first.coordinates[1] + (b2 * v[1]), first.coordinates[2] + (b2 * v[2])]

    return method === "planar" ? rhumb_distance(point.coordinates, p, units) : distance(point.coordinates, p, "degrees")
end

"""
    nearestpoint(target::Point, points::Vector{Point})

Take a reference Point and an array of Points and return the point from the
array closest to the reference. This calculation is geodesic.

# Examples
```jldoctest
julia> point = Point([35, 45])
Point([35.0, 45.0])

julia> point1 = Point([33, 42])
Point([33.0, 42.0])

julia> point2 = Point([39, 41])
Point([39.0, 41.0])

julia> point3 = Point([35, 39])
Point([35.0, 39.0])

julia> nearestpoint(point, [point1, point2, point3])
Point([33.0, 42.0])
```
"""
function nearestpoint(target::Point, points::Vector{Point})
    minDistance = Inf
    index = 0

    for (i, point) in enumerate(points)
        dist = distance(target.coordinates, point.coordinates)
        if dist < minDistance
            index = i
            minDistance = dist
        end
    end
    nearest::Point = points[index]

    return nearest
end

"""
    distance_weight(; geojson::T, treshold::Real=10000, p::Integer=2, binary::Bool=false, alpha::Real=-1., standardization::Bool=false) where {T <: AbstractFeatureCollection}

Calculate the Minkowski p-norm distance between two Points.
"""
function distance_weight(; geojson::T, treshold::Real=10000, p::Integer=2, binary::Bool=false, alpha::Real=-1., standardization::Bool=false) where {T <: AbstractFeatureCollection}
    features = []

    for feat in geojson.features
        push!(features, centroid(feat.geometry))
    end

    weights = []
    for i in 1:length(geojson.features)
        weights[i] = []
    end

    for i in 1:length(geojson.features)
        for j in i:length(geojson.features)

            i === j && (weights[i][j] = 0.)

            dist = pnorm_distance(geojson.features[i], geojson.features[j], p)
            weights[i][j] = dist
            weights[j][i] = dist
        end
    end

    for i in 1:length(geojson.features)
        for j in 1:length(geojson.features)
            dist = weights[i][j]
            dist === 0 && continue

            if binary
                if dist <= treshold
                    weights[i][j] = 1.
                else
                    weights[i][j] = 0.
                end
            else
                if dist <= treshold
                    weights[i][j] = dist^alpha
                else
                    weights[i][j] = 0.
                end

            end
        end
    end

    if standardization
        for i in 1:length(geojson.features)
            rows = reduce(+, weights[i])

            for j in eachindex(geojson.features)
                weights[i][j] = weights[i][j] / rows
            end
        end
    end

    return weights
end

"""
    pnorm_distance(point1::Point, point2::Point, p::Integer=2)::Real

Calculate the Minkowski p-norm distance between two Points.

# Examples
```jldoctest
julia> point = Point([35, 45])
Point([35.0, 45.0])

julia> point1 = Point([33, 42])
Point([33.0, 42.0])

julia> pnorm_distance(point, point1)
3.605551275463989

julia> pnorm_distance(point, point1, 1)
5.0

julia> pnorm_distance(point, point1, 3)
3.2710663101885897
```
"""
function pnorm_distance(point1::Point, point2::Point, p::Integer=2)::Real
    coords1 = point1.coordinates
    coords2 = point2.coordinates

    Δx = coords1[1] - coords2[1]
    Δy = coords1[2] - coords2[2]

    p === 1 && return abs(Δx) + abs(Δy)

    return (Δx^p + Δy^p)^(1/p)

end

"""
    midpoint(first::Point, second::Point)

Take two Points and returns a point midway between them.
The midpoint is calculated geodesically, meaning the curvature of the earth is taken into account.

# Examples
```jldoctest
julia> point = Point([35, 45])
Point([35.0, 45.0])

julia> point1 = Point([33, 42])
Point([33.0, 42.0])

julia> midpoint(point, point1)
Point([33.9751, 43.5044])
```
"""
function midpoint(first::Point, second::Point)
    dist = distance(first, second)
    bear = bearing(first.coordinates, second.coordinates)

    return destination(first.coordinates, dist / 2, bear)
end


text:
module HydrogenicAtoms
include("DiracAnalytic.jl")
end # module


text:
# Test runner script
using Base.Test

# tests

@testset "Ogre tests" begin
    include("test_common.jl")
    include("test_config.jl")
        include("test_config_defaults.jl")
    include("test_constants.jl")
    include("test_eos.jl")
    include("test_heatcapacity.jl")
    include("test_integrator.jl")
    include("test_physicalvalues.jl")
    include("test_structure.jl")
    include("test_util.jl")
end


text:
compile(l::Leaf{K}, var_map::Dict{Vector{Int}, Symbol}) where K = l.term
compile(::Fail{K}, var_map::Dict{Vector{Int}, Symbol}) where K =
		:(throw(Exception("Pattern matching failed!")))
function compile(s::Switch{K, T}, var_map::Dict{Vector{Int}, Symbol}) where {K,T}
	cases = s.cases
	cstrs = constructors(T)
	cstr_keys = constructor_keys(T)
	occ = s.occ
	function make_conditional(case_idx)
		println(var_map[occ])
		# condition
		cond = :(headof($(var_map[occ])) == $(cstr_keys[first(cases[case_idx])]))
		# body
		body = :(begin end)
		arity = length(cstrs[first(cases[case_idx])])
		for i = 1:arity
			var=get!(()->gensym(), var_map, [occ; i])
			push!(body.args, :($var = argof($(var_map[occ]), $i)))
		end
		push!(body.args, compile(last(cases[case_idx]), var_map))
		return cond, body
	end

	expr = :(begin end)
	if length(cases) > 0
		out_expr = Expr(:if)
		cond,body = make_conditional(1)
		push!(out_expr.args, cond)
		push!(out_expr.args, body)
		curr_expr = out_expr
		caseno = 2
		while caseno <= length(cases)
			next_expr = Expr(:elseif)
			push!(curr_expr.args, next_expr)
			cond,body = make_conditional(caseno)
			push!(next_expr.args, cond)
			push!(next_expr.args, body)
			curr_expr = next_expr
			caseno += 1
		end
		expr = out_expr
	end
	if !isnothing(s.default)
		push!(expr.args, compile(s.default, var_map))
	end
	return expr
end
function toplevel_compile(t::DTree{K}, inputs::Vector{T}, var_map::Dict{Vector{Int}, Symbol}) where {T,K}
	out = :(begin end)
	for i=1:length(inputs)
		var = get!(()->gensym(), var_map, [i])
		push!(out.args, :($var = $(inputs[i])))
	end
	push!(out.args, compile(t, var_map))
	return out
end


text:
module MKL

using MKL_jll

# We cannot depend on LinearAlgebra since we need to be loaded before it
if Base.USE_BLAS64
    const BlasInt = Int64
else
    const BlasInt = Int32
end

@enum Threading begin
    THREADING_INTEL
    THREADING_SEQUENTIAL
    THREADING_PGI
    THREADING_GNU
    THREADING_TBB
end
@enum Interface begin
    INTERFACE_LP64
    INTERFACE_ILP64
    INTERFACE_GNU
end

function set_threading_layer(layer::Threading = THREADING_INTEL)
    err = ccall((:MKL_Set_Threading_Layer, libmkl_rt), Cint,
        (Cint,), layer)
    if err == -1
        throw(ErrorException("return value was -1"))
    end
    return nothing
end

function set_interface_layer(interface = BlasInt == Int64 ? INTERFACE_ILP64 : INTERFACE_LP64)
    err = ccall((:MKL_Set_Interface_Layer, libmkl_rt), Cint,
        (Cint,), interface)
    if err == -1
        throw(ErrorException("return value was -1"))
    end
    return nothing
end

function __init__()
    set_threading_layer()
    set_interface_layer()
end

mklnorm(x::Vector{Float64}) = ccall((:dnrm2_, libmkl_rt), Float64,
    (Ref{BlasInt}, Ptr{Float64}, Ref{BlasInt}),
    length(x), x, 1)

# Include installation routines as well
include("install.jl")

end


text:
module BenchmarkCIExample

export naive_sum

function naive_sum(xs)
    s = 0
    for x in xs
        s += x
    end
    return s
end

end # module


text:
## ----------------------------------------------------------------------------
function plot_inner_plots(METHODS, METHODS_LABELS = string.(METHODS))

    # This is just for Nanchen
    Data = Nd
    iJR = NiJR
    src = nameof(Data)

    DAT = ChE.load_DAT(src)

    # experiments with approx D
    Dintervals = Dict(
        0.05 => 1:3,
        0.1 => 4:9,
        0.15 => 4:11,
        0.2 => 10:11,
        0.25 => 10:15,
        0.3 => 12:15,
        0.35 => 12:17,
        0.4 => 16:17,
        :tot => 1:17
    )

    int_ps = Dict()
    for (avD, EXPS) in Dintervals

        Nd_rxns_map = iJR.load_rxns_map()
        iJR_ider_subs = iJR.load_inner_rxns_subs()

        internals = iJR.load_inners_idermap() |> keys
        color_pool = Plots.distinguishable_colors(length(internals))
        ider_colors = Dict(id => c for (id, c) in zip(internals, color_pool))

        ps = Plots.Plot[]
        for (method_label, method) in zip(METHODS_LABELS, METHODS)

            p = plot(;
                xlabel = _textbf("exp abs. flux"),
                ylabel = _textbf("model abs. flux")
            )

            vals = []

            for sub in NiJR.subs()
                iJR_iders = iJR_ider_subs[sub]

                (sub in ["others"]) && continue
                marker = (10, subs_markers[sub])

                for iJR_ider in iJR_iders

                    Nd_ider = Nd_rxns_map[iJR_ider]

                    @info "Doing" Nd_ider iJR_ider

                    exp_exglcs = abs.(DAT[:exp, :flx, "GLC", EXPS])
                    exp_vals = abs.(DAT[:exp, :flx, Nd_ider, EXPS]) ./ exp_exglcs
                    model_vals = abs.(DAT[method, :flx, iJR_ider, EXPS]) ./ exp_exglcs
                    model_errs = abs.(get.([DAT], 0.0, method, :err, iJR_ider, EXPS)) ./ exp_exglcs
                    model_errs = map(model_errs) do err
                        isnan(err) && return 0.0
                        err > 1.0 && return 0.0
                        err
                    end

                    push!(vals, exp_vals..., model_vals...)

                    scatter!(p, exp_vals , model_vals;
                        label = "",
                        yerr = model_errs,
                        color = ider_colors[iJR_ider],
                        alpha = 0.8,
                        marker,
                    )

                end # for iJR_ider

            end # for (sub, iJR_iders)

            sort!(vals)
            plot!(p, [vals], [vals]; label = "",
                ls = :dash, alpha = 0.8, lw = 3,
                color = :black,
            )

            # title
            m, M = minimum(vals), maximum(vals)
            x = m + (M - m) * 0.01
            y = M - (M - m) * 0.05
            Dstr = (avD == :tot) ? "all" : round(avD; sigdigits=2)
            text = _textbf(method_label)
            annotate!(p, [(x, y, (text, 26, :left, :top, :black))])
            text = _textbf("D=", Dstr)
            annotate!(p, [(x, y - y * 0.13, (text, 26, :left, :top, :black))])

            push!(ps, p)

        end # for (method_label

        # sfig(ChE, ps,
        #     @fileid, "inner_flx_corr", (;avD), ".png",
        #     layout = (1, length(ps))
        # )

        int_ps[avD] = ps
    end

    return int_ps
end


text:
using Test
using HierarchicalUtils
using Base.Iterators: product
using Combinatorics
using HierarchicalUtils: OrderedDict
import HierarchicalUtils: _children_pairs, _childsort, _iter
import HierarchicalUtils: _node_predicate, PairVec

using Random; Random.seed!(42)

# definitions of all trees needed
abstract type AbstractVertex end

mutable struct Leaf <: AbstractVertex
    n::Int64
end
# vector of pairs as children
mutable struct VectorVertex <: AbstractVertex
    n::Int64
    chs::Vector
end
# named tuple as children
mutable struct NTVertex{T, U} <: AbstractVertex
    n::Int64
    chs::NamedTuple{T, U}
end
# vector as children
mutable struct BinaryVertex{T, U} <: AbstractVertex
    n::Int64
    ch1::T
    ch2::U
end
# tuple as children
mutable struct SingletonVertex{T} <: AbstractVertex
    n::Int64
    ch::T
end

HierarchicalUtils.@primitives
# dictionary as children
HierarchicalUtils.@hierarchical_dict
# vector as children
HierarchicalUtils.@hierarchical_vector
# tuple as children
HierarchicalUtils.@hierarchical_tuple
# vector of pairs as children
HierarchicalUtils.@hierarchical_pairvector
# named tuple as children
HierarchicalUtils.@hierarchical_namedtuple

import HierarchicalUtils: NodeType, nodeshow, nodecommshow, children

NodeType(::Type{Leaf}) = HierarchicalUtils.LeafNode()

NodeType(::Type{<:VectorVertex}) = HierarchicalUtils.InnerNode()
# we shuffle each time for the sake of testing
children(t::VectorVertex) = shuffle([k => ch for (k,ch) in enumerate(t.chs)])

NodeType(::Type{<:NTVertex}) = HierarchicalUtils.InnerNode()
children(t::NTVertex) = t.chs

NodeType(::Type{<:BinaryVertex}) = HierarchicalUtils.InnerNode()
children(t::BinaryVertex) = [t.ch1, t.ch2]

NodeType(::Type{<:SingletonVertex}) = HierarchicalUtils.InnerNode()
children(t::SingletonVertex) = (t.ch,)

NodeType(::Type{<:SingletonVertex}) = HierarchicalUtils.InnerNode()
children(t::SingletonVertex) = (t.ch,)

Base.show(io::IO, t::T) where T <: AbstractVertex = print(io, Base.nameof(T), " (", t.n, ")")
nodecommshow(io::IO, t::T) where T <: AbstractVertex = print(io, "comm")

const SINGLE_NODE_1 = NTVertex(1, NamedTuple())
const SINGLE_NODE_2 = Dict()
const SINGLE_NODE_3 = Leaf(1)
const SINGLE_NODE_4 = VectorVertex(1, AbstractVertex[])
const SINGLE_NODE_5 = AbstractVertex[]
const SINGLE_NODE_6 = Pair[]
const SINGLE_NODE_7 = ()
const SINGLE_NODES = [SINGLE_NODE_1, SINGLE_NODE_2, SINGLE_NODE_3, SINGLE_NODE_4, SINGLE_NODE_5, SINGLE_NODE_6, SINGLE_NODE_7]

const LINEAR_TREE_1 = VectorVertex(1,[
                                      NTVertex(2, (;
                                                   a=SingletonVertex(3,
                                                                     Leaf(4)
                                                                    )
                                                  ))
                                     ])

const LINEAR_TREE_2 = SingletonVertex(1,
                                      SingletonVertex(2,
                                                      VectorVertex(3, [
                                                                       NTVertex(4, NamedTuple())
                                                                      ])
                                                     )
                                     )

const LINEAR_TREE_3 = [[[[[[[Leaf(1)]]]]]]]

const COMPLETE_BINARY_TREE_1 = BinaryVertex(1,
                                            BinaryVertex(2,
                                                         Leaf(4),
                                                         Leaf(5)
                                                        ),
                                            BinaryVertex(3,
                                                         Leaf(6),
                                                         Leaf(7)
                                                        )
                                           )

const COMPLETE_BINARY_TREE_2 = NTVertex(1, (
                                            b = VectorVertex(3, [Leaf(6), NTVertex(7, NamedTuple())]),
                                            a = NTVertex(2, (b = NTVertex(5, NamedTuple()), a = VectorVertex(4, [])))
                                           ))

const T1 = NTVertex(1, (
                        ch1 = NTVertex(2, (
                                           ch1 = Leaf(4),
                                           ch2 = NTVertex(5, NamedTuple())
                                          )),
                        ch2 = BinaryVertex(3,
                                           Leaf(6),
                                           NTVertex(7, NamedTuple())
                                          )
                       ))
const T2 = NTVertex(1, (
                        ch2 = BinaryVertex(3,
                                           Leaf(6),
                                           Leaf(7)
                                          )
                        ,
                        ch1 = NTVertex(2, (
                                           ch2 = NTVertex(5, NamedTuple()),
                                          ))
                       ))
const T3 = NTVertex(1, (
                        ch1 = NTVertex(2, (
                                           ch1 = VectorVertex(4, []),
                                           ch2 = NTVertex(5, NamedTuple())
                                          )),
                       ))
const T4 = NTVertex(1, (
                        ch2 = BinaryVertex(3,
                                           Leaf(6),
                                           NTVertex(7, NamedTuple())),
                       ))
const T5 = NTVertex(1, NamedTuple())

const TEST_TREES = [
                    SINGLE_NODES...,
                    LINEAR_TREE_1, LINEAR_TREE_2, LINEAR_TREE_3,
                    COMPLETE_BINARY_TREE_1, COMPLETE_BINARY_TREE_2,
                    T1, T2, T3, T4, T5
                   ]

const TYPES = [Leaf, VectorVertex, BinaryVertex, NTVertex, Vector, Dict, Tuple, PairVec]
const ORDERS = [PreOrder(), PostOrder(), LevelOrder()]

for test_f in readdir(".")
    (endswith(test_f, ".jl") && test_f != "runtests.jl") || continue
    @testset "$test_f" begin
        include(test_f)
    end
end


text:
import Base.@deprecate

import Base:
    joinpath,
    dirname,
    ispath,
    realpath,
    normpath,
    abspath,
    relpath,
    filemode,
    isabspath,
    mkpath,
    mv,
    rm

@deprecate joinpath(root::AbstractPath, pieces::Union{AbstractPath, AbstractString}...) join(root, pieces...)
@deprecate dirname(path::AbstractPath) parent(path)
@deprecate ispath(path::AbstractPath) exists(path)
@deprecate realpath(path::AbstractPath) real(path)
@deprecate normpath(path::AbstractPath) norm(path)
@deprecate abspath(path::AbstractPath) abs(path)
@deprecate relpath(path::AbstractPath) relative(path)
@deprecate filemode(path::AbstractPath) mode(path)
@deprecate isabspath(path::AbstractPath) isabs(path)
@deprecate mkpath(path::AbstractPath) mkdir(path; recursive=true)
@deprecate mv(src::AbstractPath, dest::AbstractPath; kwargs...) move(src, dest; kwargs...)
@deprecate rm(path::AbstractPath; kwargs...) remove(path; kwargs...)


text:
# eql_root.jl

using MIRT: eql_root

using Test: @test, @test_throws, @inferred

tests = [
1 1 -1
0 1 2
4 5 0
1 0 9 # 3 both times
# 1 0 -9; # -3 both times. omitted because out of prog scope
4 4 5
] # (2x+5)(2x-1), so 1/2.

predicted = [-1, 1, 0, 3, 1/2]
results = eql_root(tests[:,1], tests[:,2], tests[:,3])

@test results == predicted
@test eql_root(1, 1, -1) == -1