text:
gate(h, n) = (1:h) .+ h*(n-1)
gate(x::AbstractVector, h, n) = @view x[gate(h,n)]
gate(x::AbstractMatrix, h, n) = x[gate(h,n),:]

# Stateful recurrence

"""
    Recur(cell)

`Recur` takes a recurrent cell and makes it stateful, managing the hidden state
in the background. `cell` should be a model of the form:

    h, y = cell(h, x...)

For example, here's a recurrent network that keeps a running total of its inputs:

```julia
accum(h, x) = (h + x, x)
rnn = Flux.Recur(accum, 0)
rnn(2)      # 2
rnn(3)      # 3
rnn.state   # 5
rnn.(1:10)  # apply to a sequence
rnn.state   # 60
```

Folding over a 3d Array of dimensions `(features, batch, time)` is also supported:

```julia
accum(h, x) = (h .+ x, x)
rnn = Flux.Recur(accum, zeros(Int, 1, 1))
rnn([2])                    # 2
rnn([3])                    # 3
rnn.state                   # 5
rnn(reshape(1:10, 1, 1, :)) # apply to a sequence of (features, batch, time)
rnn.state                   # 60
```

"""
mutable struct Recur{T,S}
  cell::T
  state::S
end

function (m::Recur)(x)
  m.state, y = m.cell(m.state, x)
  return y
end

@functor Recur
trainable(a::Recur) = (a.cell,)

Base.show(io::IO, m::Recur) = print(io, "Recur(", m.cell, ")")

"""
    reset!(rnn)

Reset the hidden state of a recurrent layer back to its original value.

Assuming you have a `Recur` layer `rnn`, this is roughly equivalent to:
```julia
rnn.state = hidden(rnn.cell)
```
"""
reset!(m::Recur) = (m.state = m.cell.state0)
reset!(m) = foreach(reset!, functor(m)[1])

# TODO remove in v0.13
function Base.getproperty(m::Recur, sym::Symbol)
  if sym === :init
    Zygote.ignore() do
      @warn "Recur field :init has been deprecated. To access initial state weights, use m::Recur.cell.state0 instead."
    end
    return getfield(m.cell, :state0)
  else
    return getfield(m, sym)
  end
end

flip(f, xs) = reverse(f.(reverse(xs)))

function (m::Recur)(x::AbstractArray{T, 3}) where T
  h = [m(view(x, :, :, i)) for i in 1:size(x, 3)]
  sze = size(h[1])
  reshape(reduce(hcat, h), sze[1], sze[2], length(h))
end

# Vanilla RNN

struct RNNCell{F,A,V,S}
  σ::F
  Wi::A
  Wh::A
  b::V
  state0::S
end

RNNCell(in::Integer, out::Integer, σ=tanh; init=Flux.glorot_uniform, initb=zeros32, init_state=zeros32) =
  RNNCell(σ, init(out, in), init(out, out), initb(out), init_state(out,1))

function (m::RNNCell{F,A,V,<:AbstractMatrix{T}})(h, x::Union{AbstractVecOrMat{T},OneHotArray}) where {F,A,V,T}
  σ, Wi, Wh, b = m.σ, m.Wi, m.Wh, m.b
  h = σ.(Wi*x .+ Wh*h .+ b)
  sz = size(x)
  return h, reshape(h, :, sz[2:end]...)
end

@functor RNNCell

function Base.show(io::IO, l::RNNCell)
  print(io, "RNNCell(", size(l.Wi, 2), ", ", size(l.Wi, 1))
  l.σ == identity || print(io, ", ", l.σ)
  print(io, ")")
end

"""
    RNN(in::Integer, out::Integer, σ = tanh)

The most basic recurrent layer; essentially acts as a `Dense` layer, but with the
output fed back into the input each time step.

The parameters `in` and `out` describe the size of the feature vectors passed as input and as output. That is, it accepts a vector of length `in` or a batch of vectors represented as a `in x B` matrix and outputs a vector of length `out` or a batch of vectors of size `out x B`.

This constructor is syntactic sugar for `Recur(RNNCell(a...))`, and so RNNs are stateful. Note that the state shape can change depending on the inputs, and so it is good to `reset!` the model between inference calls if the batch size changes. See the examples below.

# Examples
```jldoctest
julia> r = RNN(3, 5)
Recur(
  RNNCell(3, 5, tanh),                  # 50 parameters
)         # Total: 4 trainable arrays, 50 parameters,
          # plus 1 non-trainable, 5 parameters, summarysize 432 bytes.

julia> r(rand(Float32, 3)) |> size
(5,)

julia> Flux.reset!(r);

julia> r(rand(Float32, 3, 10)) |> size # batch size of 10
(5, 10)
```

!!! warning "Batch size changes"

    Failing to call `reset!` when the input batch size changes can lead to unexpected behavior. See the following example:

    ```julia
    julia> r = RNN(3, 5)
    Recur(
      RNNCell(3, 5, tanh),                  # 50 parameters
    )         # Total: 4 trainable arrays, 50 parameters,
              # plus 1 non-trainable, 5 parameters, summarysize 432 bytes.

    julia> r.state |> size
    (5, 1)

    julia> r(rand(Float32, 3)) |> size
    (5,)

    julia> r.state |> size
    (5, 1)

    julia> r(rand(Float32, 3, 10)) |> size # batch size of 10
    (5, 10)

    julia> r.state |> size # state shape has changed
    (5, 10)

    julia> r(rand(Float32, 3)) |> size # erroneously outputs a length 5*10 = 50 vector.
    (50,)
    ```
"""
RNN(a...; ka...) = Recur(RNNCell(a...; ka...))
Recur(m::RNNCell) = Recur(m, m.state0)

# TODO remove in v0.13
function Base.getproperty(m::RNNCell, sym::Symbol)
  if sym === :h
    Zygote.ignore() do
      @warn "RNNCell field :h has been deprecated. Use m::RNNCell.state0 instead."
    end
    return getfield(m, :state0)
  else
    return getfield(m, sym)
  end
end

# LSTM

struct LSTMCell{A,V,S}
  Wi::A
  Wh::A
  b::V
  state0::S
end

function LSTMCell(in::Integer, out::Integer;
                  init = glorot_uniform,
                  initb = zeros32,
                  init_state = zeros32)
  cell = LSTMCell(init(out * 4, in), init(out * 4, out), initb(out * 4), (init_state(out,1), init_state(out,1)))
  cell.b[gate(out, 2)] .= 1
  return cell
end

function (m::LSTMCell{A,V,<:NTuple{2,AbstractMatrix{T}}})((h, c), x::Union{AbstractVecOrMat{T},OneHotArray}) where {A,V,T}
  b, o = m.b, size(h, 1)
  g = m.Wi*x .+ m.Wh*h .+ b
  input = σ.(gate(g, o, 1))
  forget = σ.(gate(g, o, 2))
  cell = tanh.(gate(g, o, 3))
  output = σ.(gate(g, o, 4))
  c = forget .* c .+ input .* cell
  h′ = output .* tanh.(c)
  sz = size(x)
  return (h′, c), reshape(h′, :, sz[2:end]...)
end

@functor LSTMCell

Base.show(io::IO, l::LSTMCell) =
  print(io, "LSTMCell(", size(l.Wi, 2), ", ", size(l.Wi, 1)÷4, ")")

"""
    LSTM(in::Integer, out::Integer)

[Long Short Term Memory](https://www.researchgate.net/publication/13853244_Long_Short-term_Memory)
recurrent layer. Behaves like an RNN but generally exhibits a longer memory span over sequences.

The parameters `in` and `out` describe the size of the feature vectors passed as input and as output. That is, it accepts a vector of length `in` or a batch of vectors represented as a `in x B` matrix and outputs a vector of length `out` or a batch of vectors of size `out x B`.

This constructor is syntactic sugar for `Recur(LSTMCell(a...))`, and so LSTMs are stateful. Note that the state shape can change depending on the inputs, and so it is good to `reset!` the model between inference calls if the batch size changes. See the examples below.

See [this article](https://colah.github.io/posts/2015-08-Understanding-LSTMs/)
for a good overview of the internals.

# Examples
```jldoctest
julia> l = LSTM(3, 5)
Recur(
  LSTMCell(3, 5),                       # 190 parameters
)         # Total: 5 trainable arrays, 190 parameters,
          # plus 2 non-trainable, 10 parameters, summarysize 1.062 KiB.

julia> l(rand(Float32, 3)) |> size
(5,)

julia> Flux.reset!(l);

julia> l(rand(Float32, 3, 10)) |> size # batch size of 10
(5, 10)
```

!!! warning "Batch size changes"
    Failing to call `reset!` when the input batch size changes can lead to unexpected behavior. See the example in [`RNN`](@ref).
"""
LSTM(a...; ka...) = Recur(LSTMCell(a...; ka...))
Recur(m::LSTMCell) = Recur(m, m.state0)

# TODO remove in v0.13
function Base.getproperty(m::LSTMCell, sym::Symbol)
  if sym === :h
    Zygote.ignore() do
      @warn "LSTMCell field :h has been deprecated. Use m::LSTMCell.state0[1] instead."
    end
    return getfield(m, :state0)[1]
  elseif sym === :c
    Zygote.ignore() do
      @warn "LSTMCell field :c has been deprecated. Use m::LSTMCell.state0[2] instead."
    end
    return getfield(m, :state0)[2]
  else
    return getfield(m, sym)
  end
end

# GRU

function _gru_output(Wi, Wh, b, x, h)
  o = size(h, 1)
  gx, gh = Wi*x, Wh*h
  r = σ.(gate(gx, o, 1) .+ gate(gh, o, 1) .+ gate(b, o, 1))
  z = σ.(gate(gx, o, 2) .+ gate(gh, o, 2) .+ gate(b, o, 2))

  return gx, gh, r, z
end

struct GRUCell{A,V,S}
  Wi::A
  Wh::A
  b::V
  state0::S
end

GRUCell(in, out; init = glorot_uniform, initb = zeros32, init_state = zeros32) =
  GRUCell(init(out * 3, in), init(out * 3, out), initb(out * 3), init_state(out,1))

function (m::GRUCell{A,V,<:AbstractMatrix{T}})(h, x::Union{AbstractVecOrMat{T},OneHotArray}) where {A,V,T}
  b, o = m.b, size(h, 1)
  gx, gh, r, z = _gru_output(m.Wi, m.Wh, b, x, h)
  h̃ = tanh.(gate(gx, o, 3) .+ r .* gate(gh, o, 3) .+ gate(b, o, 3))
  h′ = (1 .- z) .* h̃ .+ z .* h
  sz = size(x)
  return h′, reshape(h′, :, sz[2:end]...)
end

@functor GRUCell

Base.show(io::IO, l::GRUCell) =
  print(io, "GRUCell(", size(l.Wi, 2), ", ", size(l.Wi, 1)÷3, ")")

"""
    GRU(in::Integer, out::Integer)

[Gated Recurrent Unit](https://arxiv.org/abs/1406.1078v1) layer. Behaves like an
RNN but generally exhibits a longer memory span over sequences. This implements
the variant proposed in v1 of the referenced paper.

The parameters `in` and `out` describe the size of the feature vectors passed as input and as output. That is, it accepts a vector of length `in` or a batch of vectors represented as a `in x B` matrix and outputs a vector of length `out` or a batch of vectors of size `out x B`.

This constructor is syntactic sugar for `Recur(GRUCell(a...))`, and so GRUs are stateful. Note that the state shape can change depending on the inputs, and so it is good to `reset!` the model between inference calls if the batch size changes. See the examples below.

See [this article](https://colah.github.io/posts/2015-08-Understanding-LSTMs/)
for a good overview of the internals.

# Examples
```jldoctest
julia> g = GRU(3, 5)
Recur(
  GRUCell(3, 5),                        # 140 parameters
)         # Total: 4 trainable arrays, 140 parameters,
          # plus 1 non-trainable, 5 parameters, summarysize 792 bytes.

julia> g(rand(Float32, 3)) |> size
(5,)

julia> Flux.reset!(g);

julia> g(rand(Float32, 3, 10)) |> size # batch size of 10
(5, 10)
```

!!! warning "Batch size changes"
    Failing to call `reset!` when the input batch size changes can lead to unexpected behavior. See the example in [`RNN`](@ref).
"""
GRU(a...; ka...) = Recur(GRUCell(a...; ka...))
Recur(m::GRUCell) = Recur(m, m.state0)

# TODO remove in v0.13
function Base.getproperty(m::GRUCell, sym::Symbol)
  if sym === :h
    Zygote.ignore() do
      @warn "GRUCell field :h has been deprecated. Use m::GRUCell.state0 instead."
    end
    return getfield(m, :state0)
  else
    return getfield(m, sym)
  end
end

# GRU v3

struct GRUv3Cell{A,V,S}
  Wi::A
  Wh::A
  b::V
  Wh_h̃::A
  state0::S
end

GRUv3Cell(in, out; init = glorot_uniform, initb = zeros32, init_state = zeros32) =
  GRUv3Cell(init(out * 3, in), init(out * 2, out), initb(out * 3),
            init(out, out), init_state(out,1))

function (m::GRUv3Cell{A,V,<:AbstractMatrix{T}})(h, x::Union{AbstractVecOrMat{T},OneHotArray}) where {A,V,T}
  b, o = m.b, size(h, 1)
  gx, gh, r, z = _gru_output(m.Wi, m.Wh, b, x, h)
  h̃ = tanh.(gate(gx, o, 3) .+ (m.Wh_h̃ * (r .* h)) .+ gate(b, o, 3))
  h′ = (1 .- z) .* h̃ .+ z .* h
  sz = size(x)
  return h′, reshape(h′, :, sz[2:end]...)
end

@functor GRUv3Cell

Base.show(io::IO, l::GRUv3Cell) =
  print(io, "GRUv3Cell(", size(l.Wi, 2), ", ", size(l.Wi, 1)÷3, ")")

"""
    GRUv3(in::Integer, out::Integer)

[Gated Recurrent Unit](https://arxiv.org/abs/1406.1078v3) layer. Behaves like an
RNN but generally exhibits a longer memory span over sequences. This implements
the variant proposed in v3 of the referenced paper.

The parameters `in` and `out` describe the size of the feature vectors passed as input and as output. That is, it accepts a vector of length `in` or a batch of vectors represented as a `in x B` matrix and outputs a vector of length `out` or a batch of vectors of size `out x B`.

This constructor is syntactic sugar for `Recur(GRUv3Cell(a...))`, and so GRUv3s are stateful. Note that the state shape can change depending on the inputs, and so it is good to `reset!` the model between inference calls if the batch size changes. See the examples below.

See [this article](https://colah.github.io/posts/2015-08-Understanding-LSTMs/)
for a good overview of the internals.

# Examples
```jldoctest
julia> g = GRUv3(3, 5)
Recur(
  GRUv3Cell(3, 5),                      # 140 parameters
)         # Total: 5 trainable arrays, 140 parameters,
          # plus 1 non-trainable, 5 parameters, summarysize 848 bytes.

julia> g(rand(Float32, 3)) |> size
(5,)

julia> Flux.reset!(g);

julia> g(rand(Float32, 3, 10)) |> size # batch size of 10
(5, 10)
```

!!! warning "Batch size changes"
    Failing to call `reset!` when the input batch size changes can lead to unexpected behavior. See the example in [`RNN`](@ref).
"""
GRUv3(a...; ka...) = Recur(GRUv3Cell(a...; ka...))
Recur(m::GRUv3Cell) = Recur(m, m.state0)

@adjoint function Broadcast.broadcasted(f::Recur, args...)
  Zygote.∇map(__context__, f, args...)
end

"""
    Bidirectional{A,B}

A wrapper layer that allows the use of [bidirectional](https://ieeexplore.ieee.org/document/650093) layers. It contains two parts that are Flux layers: `forward` and `backward` where
the forward layer weights are concatenated with the reversed order of the backward layer weights.

It is intended to be used with recurrent layers such as `LSTM`, `GRU` or `RNN` to benefit from the sequential information that recurrent
layers have, but it will not raise an error if used with a different layer such as `Dense`, as long as the layer is compatible with the concatenation function `vcat`.

For flexibility, it is possible to use the contructor `Bidirectional(rnn, in::Int, out::Int, a...; ka...)` by passing the input and output dimensions
together with the desired recurrent layers (one of `LSTM`, `GRU`, or `RNN`). Check the examples below for more details.

# Examples

1. Using the flexible constructor to create a bidirectional LSTM layer (BiLSTM):
```jldoctest
julia> BLSTM = Bidirectional(LSTM, 3, 5)
Bidirectional(
  Recur(
    LSTMCell(3, 5),                     # 190 parameters
  ),
  Recur(
    LSTMCell(3, 5),                     # 190 parameters
  ),
)         # Total: 10 trainable arrays, 380 parameters,
          # plus 4 non-trainable, 20 parameters, summarysize 2.141 KiB.
```
2. Checking the dimension after running the bidirectional layer on an input vector. Shows that the
dimension of the output is twice the dimension of the input:
```jldoctest
julia> Bidirectional(LSTM, 3, 5)(rand(Float32, 3)) |> size
(10,)
```
3. It is possible to use the bidirectional layer inside the `Chain` container:
```
julia> model = Chain(Embedding(10000, 200), Bidirectional(LSTM, 200, 128), Dense(256, 5), softmax)
Chain(
  Embedding(10000, 200),                # 2_000_000 parameters
  Bidirectional(
    Recur(
      LSTMCell(200, 128),               # 168_704 parameters
    ),
    Recur(
      LSTMCell(200, 128),               # 168_704 parameters
    ),
  ),
  Dense(256, 5),                        # 1_285 parameters
  NNlib.softmax,
)         # Total: 13 trainable arrays, 2_338_693 parameters,
          # plus 4 non-trainable, 512 parameters, summarysize 8.922 MiB.
```
4. It is also possible to use the default constructor
```jldoctest
julia> BiLSTM = Bidirectional(GRU(3, 5), GRU(3, 5))
Bidirectional(
  Recur(
    GRUCell(3, 5),                      # 140 parameters
  ),
  Recur(
    GRUCell(3, 5),                      # 140 parameters
  ),
)         # Total: 8 trainable arrays, 280 parameters,
          # plus 2 non-trainable, 10 parameters, summarysize 1.562 KiB.
```
5. And use other parameters available on the recurrent layers
```jldoctest
julia> BiLSTM = Bidirectional(RNN(3, 5, tanh; init=glorot_normal), LSTM(3, 5; initb=zeros32, init_state=zeros32))
Bidirectional(
  Recur(
    RNNCell(3, 5, tanh),                # 50 parameters
  ),
  Recur(
    LSTMCell(3, 5),                     # 190 parameters
  ),
)         # Total: 9 trainable arrays, 240 parameters,
          # plus 3 non-trainable, 15 parameters, summarysize 1.500 KiB.
```
"""
struct Bidirectional{A,B}
  forward::A
  backward::B
end

# Constructor that creates a bidirectional with the same layer for forward and backward
# Needs to have `in` explicitly declared to avoid conflicts with the default construtor
Bidirectional(rnn, in::Integer, args...; ka...) = Bidirectional(rnn(in, args...; ka...), rnn(in, args...; ka...))

# Concatenate the forward and reversed backward weights
function (m::Bidirectional)(x::Union{AbstractVecOrMat{T},OneHotArray}) where {T}
  return vcat(m.forward(x), reverse(m.backward(reverse(x; dims=1)); dims=1))
end

@functor Bidirectional

function Base.show(io::IO, b::Bidirectional)
  print(io, "Bidirectional(")
  show(io, b.forward.cell)
  print(io, ", ")
  show(io, b.backward.cell)
  print(io, ")")
end;


text:
module Sim3
include("./lib/probs.jl")
include("./lib/pccf.jl")
include("./detector.jl")
include("./lib/roc.jl")
include("./lib/utils.jl")
#using Probs, Detector, Pccf

function calcTmax(mu, n)
    round(Int,(mu*(n+1)))
end

function calcWidth(mu)
    1 + round(Int, 1.33 * mu)
end

# 1
function changesLocations(n, mu, sigma)
    changes1::Array{Int,1}
    tMax = calcTmax(mu,n)
    width = calcWidth(mu)
    changes1 = Probs.generateChanges(n,mu,sigma,width)
end
# 2
function errorsLocations(n,mu,errRate)
    tMax = calcTmax(mu,n)
    avgNumOfErrors = errRate * n
    #round(Int, rand(1:tMax, avgNumOfErrors))
    rand(1:tMax, avgNumOfErrors)
end

## ;;;;;;;;;;;;;;;;;Pccf/v1.0/src/recurrency_modelling2/sim3.clj;;;;;;;
## (defn gen-vals-in-interval-uniform [n from-v to-v]
##   "Function to generate values within interval inclusive"
##   (let [rng (inc (- to-v from-v))]
##     (map #(+ % from-v) (repeatedly n #(rand-int rng)))))
##
## (defn gen-vals-in-interval-gauss [n from-v to-v]
##   (let [mu (/ (+ to-v from-v) 2.0) sigma (inc (- to-v from-v))]
##     (utils/rnorm n mu sigma)))
## ;; (gen-vals-in-interval 5 7 8)
## ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

function genValsInIntervalGauss(n,st,en)
    # used to simulate changes
    mu = (en+st)/2 #old: (n+st)/2
    sigma = en-st+1
    randn(n) * sigma + mu
end

function updVec(vec, ids::Array{Int,1}, newVals)
    v=copy(vec)
    for i in 1:length(ids)
        v[ids[i]]=newVals[i]
    end
    return v
end

function genSig(n,mu,sigma,errRate,tvf,tvt,evf,evt)
    # Generate signal. For changes values are from tvf-tvt; errors in rvf-evt
    tMax        = calcTmax(mu,n)
    changesLoc  = changesLocations(n,mu,sigma)
    errorsLoc   = errorsLocations(n,mu,errRate)
    nChps       = length(changesLoc)
    nErrors     = length(errorsLoc)
    valsChanges = genValsInIntervalGauss(nChps, tvf, tvt)
    valsErrors  = genValsInIntervalGauss(nErrors, evf, evt)
    vals        = zeros(tMax)
    changes1    = updVec(vals, errorsLoc, valsErrors)
    changes2    = updVec(changes1, changesLoc, valsChanges)
    Dict("changes"=>changesLoc, "sig"=>changes2)
end

function avgPairs(listOfTuples)
    f=[x[1] for x in listOfTuples]
    s=[x[2] for x in listOfTuples]
    (mean(f),mean(s))
end
# julia> avgPairs([(1,2),(5,7)])-> (3.0,4.5)

function simForH(nSim,nChps, mu,sigma,errRate,tvFrom, tvTo, evFrom,evTo,h)
    # Perform n simulations for a given threshold h
    performance=[]
    for i in 1:nSim
        dat=genSig(nChps,mu,sigma,errRate,tvFrom,tvTo,evFrom,evTo)
        newsignal=dat["sig"]
        tChps=dat["changes"]
        maxDelay=1
        detections=Detector.naiveDetector(newsignal,h)
        tprate=Roc.tpRate(detections,tChps,maxDelay)
        farate=Roc.faRate(detections,tChps,maxDelay,length(newsignal))
        push!(performance,(tprate,farate))
    end
    avgPairs(performance)
end

function simDynamicH(probsFile,nSim,nChps,mu,sigma,errRate,tvFrom,tvTo,
                     evFrom,evTo,delayInMu,hmin,hmax)
    # Calculate TP/FA rates for a given pair of hmin hmax
    maxDelay=1
    tMax=calcTmax(mu,nChps)
    hForProbsFunction=0.1
    T=round(Int, mu*nChps*2)
    probsNoConfirm=Pccf.pccf(T,mu,sigma)
    performance=[]
    for i in 1:nSim
        dat = genSig(nChps,mu,sigma,errRate,tvFrom,tvTo,evFrom,evTo)
        newsignal = dat["sig"]
        tChps = dat["changes"]
        subProbsNoConfirm = probsNoConfirm[1:length(newsignal)]
        probs = Pccf.pccfWithConfirmation(subProbsNoConfirm,delayInMu,tChps)
        dynamicHvec = Probs.dynamicH(probs,hForProbsFunction,hmin,hmax, tMax)

        detections = Detector.naiveDetectorDynamic(newsignal,dynamicHvec)

        tprate = Roc.tpRate(detections, tChps, maxDelay)
        farate = Roc.faRate(detections, tChps, maxDelay,length(newsignal))
        push!(performance, (tprate,farate))
    end
    avgPairs(performance)
end

function simNotAdaptive(nSim,nChps,mu,sigma,
                        errRate,tvFrom,tvTo,evFrom,
                        evTo,h1,h2)
    # For every threshold h from h1 to h2 perform n-sim simulations
    performance=[]
    h=h1
    while h<=h2
        avgPerf=simForH(nSim,nChps,mu,sigma,errRate,tvFrom,tvTo,evFrom,evTo, float(h))
        push!(performance, (h, avgPerf[1],avgPerf[2]))
        h=h+1
    end
    performance
end

function runSimNotAdaptive(nSim, mainpath, settings)
    s=settings
    notAdaptivePerf=simNotAdaptive(nSim,s["nchps"],s["mu"],s["sigma"],
                                   s["errRate"],s["tvFrom"],s["tvTo"],s["evFrom"],
                                   s["evTo"],20,50)

    write_items(mainpath*"NotadAptivePerf.txt", notAdaptivePerf, "Delay,TP,FA")
end

function runSimAdaptive(nSim,delay1,delay2,hmin,hmax,mainpath,settings)
    s=settings
    performance=[]
    for confirmDelay in delay1:delay2
        adaptivPerf=simDynamicH(mainpath*"recprobs.txt",nSim,s["nchps"],s["mu"],
                                s["sigma"],s["errRate"],s["tvFrom"],s["tvTo"],
                                s["evFrom"],s["evTo"],confirmDelay,hmin,hmax)
        push!(performance, (confirmDelay,adaptivPerf[1],adaptivPerf[2]))
    end
    outpath=mainpath*"DelayVsPerformance.txt"
    write_items(outpath, performance, "Delay,TP,FA")
    println("... runSimAdaptive(): results are written into $outpath")
end

function run(settings)
    mainpath="D:/_tmp/out/"
    stCommon=Dict("nsim"=>200)
    stAdaptive=Dict("delay1"=>1, "delay2"=>3,"hmin"=>27.0,"hmax"=>30.0)
    runSimNotAdaptive(stCommon["nsim"], mainpath, settings)
    runSimAdaptive(stCommon["nsim"],stAdaptive["delay1"],stAdaptive["delay2"],
                   stAdaptive["hmin"],stAdaptive["hmax"], mainpath, settings)
end

end

ext: jl

max_stars_repo_name: av-maslov/Recurrency

max_issues_repo_name: av-maslov/Recurrency

max_forks_repo_name: av-maslov/Recurrency


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct Map <: Edit
        file::String
        from::Array{Symbol,1}
        to::Array{Symbol,1}
        input::Array{Symbol,1}
        output::Array{Symbol,1}
        kind::Symbol
    end

Define an `output` column containing values based on those in an `input` column. The mapping columns `from` -> `to` are contained in a .csv `file` in the coremaps directory. The columns `input` and `from` should contain the same values, as should `output` and `to`.

# Arguments
- `file::String`: mapping .csv file name in the coremaps directory
- `from::Array{Symbol,1}`: name of the mapping column containing input values
- `to::Array{Symbol,1}`: name of the mapping column containing output values
- `input::Array{Symbol,1}`: name of the input column to map
- `output::Array{Symbol,1}`: name of the output column created
- `kind::Symbol`: type of join to perform.
"""
mutable struct Map <: Edit
    "mapping .csv file name in the coremaps directory"
    file::String
    "name of the mapping column containing input values"
    from::Array{Symbol,1}
    "name of the mapping column containing output values"
    to::Array{Symbol,1}
    "name of the input column to map"
    input::Array{Symbol,1}
    "name of the output column created"
    output::Array{Symbol,1}
    "type of join to perform."
    kind::Symbol
end

function Map(; file, from, to, input, output, kind, )
    Map(file, from, to, input, output, kind, )
end

"""Get [`Map`](@ref) `file`."""
get_file(value::Map) = value.file
"""Get [`Map`](@ref) `from`."""
get_from(value::Map) = value.from
"""Get [`Map`](@ref) `to`."""
get_to(value::Map) = value.to
"""Get [`Map`](@ref) `input`."""
get_input(value::Map) = value.input
"""Get [`Map`](@ref) `output`."""
get_output(value::Map) = value.output
"""Get [`Map`](@ref) `kind`."""
get_kind(value::Map) = value.kind

"""Set [`Map`](@ref) `file`."""
set_file!(value::Map, val) = value.file = val
"""Set [`Map`](@ref) `from`."""
set_from!(value::Map, val) = value.from = val
"""Set [`Map`](@ref) `to`."""
set_to!(value::Map, val) = value.to = val
"""Set [`Map`](@ref) `input`."""
set_input!(value::Map, val) = value.input = val
"""Set [`Map`](@ref) `output`."""
set_output!(value::Map, val) = value.output = val
"""Set [`Map`](@ref) `kind`."""
set_kind!(value::Map, val) = value.kind = val

max_stars_repo_name: matbesancon/SLiDE

max_issues_repo_name: matbesancon/SLiDE

max_forks_repo_name: matbesancon/SLiDE


text:
function bar_PV(Data_df, date)
    # bars
    @df Data_df[(Data_df[:,:day].==Dates.day(date)) .&(Data_df[:,:month].==Dates.month(date)),:] groupedbar(
        [:PV_HP :PV_GR :PV_B :PV_DE], color =[:firebrick :grey :purple :orange],
        ylim=(0,15), yticks=0:5:15, xticks=0:6:24,
        annotations=(2,12, text("$(Dates.month(date))/$(Dates.day(date))", 10)), legendfontsize=6,
        label=["PV_HP" "PV_GR" "PV_B" "PV_DE"], legend=false, bar_position = :stack, alpha=0.8);
    # battery cap
    @df Data_df[(Data_df[:,:day].==Dates.day(date)) .&(Data_df[:,:month].==Dates.month(date)),:] plot!(
        [:Soc_B], color =[:purple], ylim=(0,15), yticks=0:2:15, xticks=0:6:24,
        label=["SOC_b"], legend=false, linewidth= 2.0);
    # PV generation
    @df Data_df[(Data_df[:,:day].==Dates.day(date)) .&(Data_df[:,:month].==Dates.month(date)),:] plot!(
        [:PV_generation], ylim=(0,15), yticks=0:5:15, xticks=0:6:24, color =[:gold],
        label=["g_e"], legend=false, linewidth= 2.0);

    plot!(0:23, ones(24,1)*13.5, linestyle=:dot, linewidth= 2, color=:purple)
end

function bar_demand(Data_df, date)
    # bars
    @df Data_df[(Data_df[:,:day].==Dates.day(date)) .&(Data_df[:,:month].==Dates.month(date)),:] groupedbar(
        [:B_DE :GR_DE :PV_DE ], color =[:purple :grey :gold],
        ylim=(0,2), yticks=0:0.5:2, xticks=0:6:24,
        annotations=(2,1.8, text("$(Dates.month(date))/$(Dates.day(date))", 10)), legendfontsize=6,
        label=["B_DE" "GR_DE" "PV_DE"], legend=false, bar_position = :stack, alpha=0.8);
    # electricity demand
    @df Data_df[(Data_df[:,:day].==Dates.day(date)) .&(Data_df[:,:month].==Dates.month(date)),:] plot!(
        [:electkwh], ylim=(0,2), yticks=0:0.5:2, xticks=0:6:24, color =[:orange],
        label=["d_e"], legend=false, linewidth= 2.0);

    plot!(0:23, ones(24,1)*13.5, linestyle=:dash, color=:purple)
end

function bar_heat(Data_df, date)
    # hot water demand
    @df Data_df[(Data_df[:,:day].==Dates.day(date)) .&(Data_df[:,:month].==Dates.month(date)),:] plot(
        [:hotwaterkwh]+[:heatingkwh], fillrange=[:heatingkwh], ylim=(0,15), yticks=0:5:15, xticks=0:6:24,
        color =[:steelblue], label=["d_fh"], legend=false, alpha=0.2);
    # heat demand
    @df Data_df[(Data_df[:,:day].==Dates.day(date)) .&(Data_df[:,:month].==Dates.month(date)),:] plot!(
        [:heatingkwh], fillrange=zeros(24), ylim=(0,15), yticks=0:5:15, xticks=0:6:24,
        color =[:firebrick], label=["d_fh"], legend=false, alpha=0.2);
    # bars
    @df Data_df[(Data_df[:,:day].==Dates.day(date)) .&(Data_df[:,:month].==Dates.month(date)),:] groupedbar!(
        [:B_HP :GR_HP :PV_HP], color =[:purple :grey :gold],
        ylim=(0,15), yticks=0:5:15, xticks=0:6:24,
        annotations=(2, 14, text("$(Dates.month(date))/$(Dates.day(date))", 10)), legendfontsize=6,
        label=["B_DE" "GR_DE" "PV_DE"], legend=false, bar_position = :stack, alpha=0.8);

    plot!(0:23, ones(24,1)*3, linestyle=:dot, linewidth= 2, color=:grey)
end

function bar_comfort(Data_df, date)
    @df Data_df[(Data_df[:,:day].==Dates.day(date)) .&(Data_df[:,:month].==Dates.month(date)),:] groupedbar(
        ([:HP_FH :HP_HW].>0.005)*200, color =[:firebrick :steelblue],
        ylim=(19,23), yticks=20:1:22, xticks=0:6:24,
        annotations=(2, 22.5, text("$(Dates.month(date))/$(Dates.day(date))", 10)), legendfontsize=6,
        label=["Mod_fh" "Mod_hw"], legend=false, bar_position = :stack, alpha=0.15);

    # state-of-charge fh
    @df Data_df[(Data_df[:,:day].==Dates.day(date)) .&(Data_df[:,:month].==Dates.month(date)),:] plot!(
        [:Temp_FH], ylim=(19,23), yticks=20:1:22, xticks=0:6:24, color =[:firebrick],
        label=["T_fh"], legend=false, linewidth= 2.0);

    plot!(0:23, ones(24,1)*20, linestyle=:dot, linewidth= 2, color=:firebrick)
    plot!(0:23, ones(24,1)*22, linestyle=:dot, linewidth= 2, color=:firebrick)

    plt = twinx()
    # state-of-charge hw
    @df Data_df[(Data_df[:,:day].==Dates.day(date)) .&(Data_df[:,:month].==Dates.month(date)),:] plot!(
        plt, [:Vol_HW], ylim=(10,190), yticks=20:40:180, xticks=0:6:24, color =[:steelblue],
        label=["V_fh"], legend=false, linewidth= 2.0, linestyle=:dash, grid=false);

    plot!(plt, 0:23, ones(24,1)*20, linestyle=:dot, linewidth= 2, color=:steelblue)
    plot!(plt, 0:23, ones(24,1)*180, linestyle=:dot, linewidth= 2, color=:steelblue)
end

function bar_row(plotfunction, Data_df, date, length)
    plot([plotfunction(Data_df, (date+Day(i-1))) for i in 1:length]..., layout=(1,length), size=(300*length,200));
end


text:
# Tests to simulate multivariate traits

## Get a network and make it ultrametric
net = readTopology("(((Ag:5,(#H1:1::0.056,((Ak:2,(E:1,#H2:1::0.004):1):1,(M:2)#H2:1::0.996):1):1):1,(((((Az:1,Ag2:1):1,As:2):1)#H1:1::0.944,Ap:4):1,Ar:5):1):1,(P:4,20:4):3,165:7);");
trait_dim = 3

@testset "Simulate data and check means and dimensions" begin

## Simulate a MBD
Random.seed!(17920921); # fix the seed
μ = randn(trait_dim)
Σ = randn(trait_dim, trait_dim)
Σ = Σ * Σ' # needs to be positive definite
pars = ParamsMultiBM(μ, Σ); # params of a MBD
@test_logs show(devnull, pars)

sim = simulate(net, pars); # simulate according to a BM
@test_logs show(devnull, sim)

# Extract simulated values
traitsTips = sim[:Tips];
traitsNodes = sim[:InternalNodes];

# Check dimensions
@test size(traitsTips) == (trait_dim, net.numTaxa)
@test size(traitsNodes) == (trait_dim, net.numNodes - net.numTaxa)

expectations = sim[:All, :Exp]

# Check means (no shifts)
max_diff = maximum(expectations - μ * ones(net.numNodes)')
@test max_diff ≈ 0.0 atol=1e-10

end

###############################################################################
## Test of distibution
###############################################################################
@testset "Simulate test distribution" begin

## Generate some values
Random.seed!(18480224); # fix the seed
μ = randn(trait_dim)
Σ = randn(trait_dim, trait_dim)
Σ = Σ * Σ' # needs to be positive definite
pars = ParamsMultiBM(μ, Σ); # params of a MBD

N = 50000
S = length(tipLabels(net));
μ_sim = zeros(trait_dim, S)
Σ_sim = zeros(trait_dim * S, trait_dim * S)
for i = 1:N
    tips = simulate(net, pars)[:Tips]
    μ_sim .+= tips
    v_sim = vec(tips)
    Σ_sim += v_sim * v_sim'
end

μ_sim ./= N
Σ_sim ./= N
Σ_sim = Σ_sim - vec(μ_sim) * vec(μ_sim)'

## Check means

μ_true = μ * ones(S)'

μ_max = maximum(abs.(μ_true - μ_sim))
@test μ_max < 1e-1

## Check covariance

Ψ = Matrix(vcv(net))
Σ_true = kron(Ψ, Σ)
Σ_max = maximum(abs.(Σ_true - Σ_sim))
@test Σ_max < 2e-1

end


text:
struct CyclicContainer <: AbstractVector{String}
    c::Vector{String}
    n::Int
end
CyclicContainer(c) = CyclicContainer(c, 0)

Base.length(c::CyclicContainer) = length(c.c)
Base.size(c::CyclicContainer) = size(c.c)
Base.getindex(c::CyclicContainer, i) = c.c[mod1(i, length(c.c))]
function Base.getindex(c::CyclicContainer)
    c.n += 1
    c[c.n]
end
Base.iterate(c::CyclicContainer, i = 1) = iterate(c.c, i)

COLORS = [
    "#6D44D0",
    "#2CB3BF",
    "#1B1B1B",
    "#DA5210",
    "#866373",
    "#03502A"
]

COLORSCHEME = COLORS

CCOLORS = CyclicContainer(COLORS)
NONBLACK = CyclicContainer(COLORS[[1,2,4,5]])
BLACK = COLORS[3]

LINESTYLES = CyclicContainer(["-", "--", "-.",  ":",])

figx = 16 # default width correspoding to full text width, in inches
figy = 5  # default height corresponding to 1 row of plots, in inches

export COLORSCHEME, COLORS, CCOLORS, LINESTYLES, figx, figy

max_stars_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_issues_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_forks_repo_name: JuliaDynamics/NonlinearDynamicsTextbook


text:
module Oracle

import Tables
using Dates

const DEPS_FILE = joinpath(@__DIR__, "..", "deps", "deps.jl")
if !isfile(DEPS_FILE)
    error("Oracle.jl is not installed properly, run Pkg.build(\"Oracle\") and restart Julia.")
end
include(DEPS_FILE)

# see issue #21
#include("oranumbers/oranumbers.jl")
#import .OraNumbers.OraNumber

include("macros.jl")
include("constants.jl")
include("enums.jl")
include("exec_mode.jl")
include("types.jl")
include("odpi.jl")
include("timestamps.jl")
include("ora_timestamp.jl")
include("context.jl")
include("connection.jl")
include("object.jl")
include("oracle_value.jl")
include("stmt.jl")
include("bind.jl")
include("lob.jl")
include("cursor.jl")
include("pool.jl")
include("variable.jl")
include("tables_interface.jl")
include("messages.jl")
include("queue.jl")

function __init__()
    # this function is defined in DEPS_FILE
    check_deps()

    # check size of structs affected by C unions
    @assert sizeof(OraDataBuffer) == sizeof_dpiDataBuffer() "OraDataBuffer should have sizeof $(sizeof_dpiDataBuffer()) bytes"
    @assert sizeof(OraData) == sizeof_dpiData() "OraData should have sizeof $(sizeof_dpiData()) bytes"
    @assert sizeof(OraTimestamp) == sizeof_dpiTimestamp() "OraTimestamp should have sizeof $(sizeof_dpiTimestamp()) bytes"
    @assert sizeof(OraErrorInfo) == sizeof_dpiErrorInfo() "OraErrorInfo should have sizeof $(sizeof_dpiErrorInfo()) bytes"
    @assert sizeof(OraCommonCreateParams) == sizeof_dpiCommonCreateParams() "OraCommonCreateParams should have size of $(sizeof_dpiCommonCreateParams()) bytes"
    @assert sizeof(OraAppContext) == sizeof_dpiAppContext() "OraAppContext should have size of $(sizeof_dpiAppContext()) bytes"
    @assert sizeof(OraPoolCreateParams) == sizeof_dpiPoolCreateParams() "OraPoolCreateParams should have sizeof $(sizeof_dpiPoolCreateParams()) bytes"
    @assert sizeof(OraConnCreateParams) == sizeof_dpiConnCreateParams() "OraConnCreateParams should have sizeof $(sizeof_dpiConnCreateParams()) bytes"
    @assert sizeof(OraDataTypeInfo) == sizeof_dpiDataTypeInfo() "OraDataTypeInfo should have size of $(sizeof_dpiDataTypeInfo()) bytes"
    @assert sizeof(OraQueryInfo) == sizeof_dpiQueryInfo() "OraQueryInfo should have sizeof $(sizeof_dpiQueryInfo()) bytes"
    @assert sizeof(OraStmtInfo) == sizeof_dpiStmtInfo() "OraStmtInfo should have sizeof $(sizeof_dpiStmtInfo()) bytes"
    @assert sizeof(OraVersionInfo) == sizeof_dpiVersionInfo() "OraVersionInfo should have sizeof $(sizeof_dpiVersionInfo()) bytes"
    @assert sizeof(OraBytes) == sizeof_dpiBytes() "OraBytes should have sizeof $(sizeof_dpiBytes()) bytes"
    @assert sizeof(OraEncodingInfo) == sizeof_dpiEncodingInfo() "OraEncodingInfo should have sizeof $(sizeof_dpiEncodingInfo()) bytes"
    @assert sizeof(OraObjectTypeInfo) == sizeof_dpiObjectTypeInfo() "OraObjectTypeInfo should have sizeof $(sizeof_dpiObjectTypeInfo()) bytes"
    #@assert sizeof(OraNumber) == sizeof_dpiNumber() # see issue #21
end

@inline function error_check(ctx::Context, result::OraResult)
    if result == ORA_FAILURE
        error_info_ref = Ref{OraErrorInfo}()
        dpiContext_getError(ctx.handle, error_info_ref)
        error_info = error_info_ref[]
        throw(error_info)
    end
    @assert result == ORA_SUCCESS

    nothing
end

function odpi_version(vnum::Integer=odpi_version_number()) :: VersionNumber
    t1 = 100 # threshold 1
    t2 = 10000 # threshold 2
    major = div(vnum, t2)
    minor = div(mod(vnum, t2), t1)
    patch = mod(vnum, t1)
    return VersionNumber("$major.$minor.$patch")
end

end # module Oracle


text:
function inference_svi_gp(m::GDTM, num_inducing; rand_inducing::Bool=false, normalize_timestamps::Bool=false)
	tic()

	e_step_likelihoods = Float64[]
	m_step_likelihoods = Float64[]
	full_elbo = Float64[]

	test_schedule = 1
	cur_count = 0

	# size of the minibatch used to compute updates and samples
	mini_batch_size = m.batch_size

	# set the number of epochs (needs some kind of manipulation, monitor elbo, map ?)
	# println("handling $(m.D) docs")
	epochs = round(Int64, m.D/mini_batch_size)
	info("doing $epochs epochs with minibatch size $mini_batch_size")
	info("parameter:\nKernel: $(m.krn)\nalpha: $(m.alpha)\nK: $(m.K)")

	# timestamp normalization is disabled by default
	if normalize_timestamps
		m.times .-= minimum(m.times) - 1
		m.times = m.times ./ maximum(m.times)
	end

	# random inducing point selection is disabled by default
	# with a roughly uniform distribution of documents over timestamps, even spacing works better
	if num_inducing >= m.T
		# if number of inducing points is larger than total timestamps, use all (degenerating into a full rank model)
		m.inducing_points = m.times
		num_inducing = m.T
	elseif rand_inducing
		# randomly select induing points from time points
		m.inducing_points = m.times[randperm(m.T)][1:num_inducing]
	else
		m.inducing_points = collect(linspace(minimum(m.times),maximum(m.times), num_inducing))
	end
	# compute number of all tokens in the training set, used for likelihood smoothing
	all_tokens = sum(map(d->sum(collect(values(d.value))), m.corpus.documents[m.training_doc_ids]))

	# compute inducing points covariance matrix
	m.Kmm = Hermitian(Kernels.computeCovarianceMatrix(m.krn, m.inducing_points))

	m.jitter = Diagonal(ones(num_inducing) .* 1e-7)
	m.KmmInv = inv(Hermitian(m.Kmm + m.jitter))
	jitter_full = Diagonal(ones(m.T) .* 1e-7)

	# compute training covariance matrix
	m.Knn = Hermitian(Kernels.computeCovarianceMatrix(m.krn, m.times) + eye(m.T).*m.s_x + jitter_full)

	# compute cross covariance and some intermediate results
	m.Knm = Kernels.computeCrossCovarianceMatrix(m.krn, m.times, m.inducing_points)
	m.KnmKmmInv = m.Knm*m.KmmInv

	m.K_tilde = Hermitian(m.Knn - m.KnmKmmInv*m.Knm')
	m.K_tilde_diag = diag(m.K_tilde)

	# optionally plot heatmaps of covariance matrices for debugging
	if m.visualize
		display(Plots.heatmap(m.Knn, title="original covariance"))
		display(Plots.heatmap(m.Knm*m.KnmKmmInv', title="reconstructed covariance"))
		display(Plots.heatmap(m.Kmm, title="inducing point covariance"))
		display(Plots.heatmap(m.K_tilde, title="reconstruction error"))
	end

	info("initializing the model...")

	# initialize means, covariances and natural parameters
	m.μ = Vector{Matrix{Float64}}(m.K)
	m.s = Matrix{Matrix{Float64}}(m.K, m.V)
	m.eta_1 = Vector{Matrix{Float64}}(m.K)
	m.eta_2 = Matrix{Matrix{Float64}}(m.K, m.V)

	L = chol(m.Knn)'

	# for mean initialization, use document seeding (one doc per topic and time) or draw randomly from a GP with covariance as given by the timestamps
	if m.use_seeding
		Threads.@threads for t in 1:m.T
			all_docs_t = DynamicData.get_datapoints_for_timestamp(m.corpus.documents, m.timestamps[t])
			p = zeros(Int64, m.K)
			if length(all_docs_t) < m.K
				for i in 1:m.K
					p[i] = rand(1:length(all_docs_t))
				end
			else
				p = randperm(length(all_docs_t))[1:m.K]
			end
			mean = zeros(m.V)
			for k in 1:m.K
				doc = m.corpus.documents[p[k]].value
				(words, freqs) =  CorpusUtils.get_words_and_counts(doc)
				mean = randn(m.V).*m.s_0 + m.m_0
				mean[words] = log(freqs .+ rand(length(words)))
				m.μ[k] = m.KnmKmmInv[t,:]' \ mean'
			end
		end
	end

	# initialize variational covariances and natural parameters to topic distributions
	@inbounds for k in 1:m.K
		if !m.use_seeding
			m.μ[k] = zeros(num_inducing, m.V)
		end
		m.eta_1[k] = zeros(num_inducing, m.V)

		Threads.@threads for w in 1:m.V
			if !m.use_seeding
				m.μ[k][:,w] = m.KnmKmmInv \ (L*randn(m.T))
			end
			m.s[k,w] = eye(num_inducing)
			m.eta_2[k,w] = -.5*eye(num_inducing)
			m.eta_1[k][:,w] = -2*m.eta_2[k,w] * m.μ[k][:,w]
		end
	end

	# assume s is identity matrix for init
	# compute the Λ diagonals and update 𝜁 parameters
	init_val = diag(m.KnmKmmInv*m.KnmKmmInv')
	Λ_diags = zeros(m.T, m.V)
	for k in 1:m.K
		Threads.@threads for w in 1:m.V
			Λ_diags[:,w] = copy(init_val)
		end
		means = m.KnmKmmInv*m.μ[k]
		svi_update_zeta(m, k, means, Λ_diags)
	end

	# parameters to steer the step size in stochastic gradient updates
	a = .1
	b = 10
	gamma = .7

	# helpers for performing inference loop
 	iter = 0
  	last_e_step_ll = -1e100
  	svi_counter = 1
	mb_idx_rand = m.training_doc_ids[randperm(m.D)]
	e_step_time_agg = 0.
	m_step_time_agg = 0.
	info("done, starting inference")
	@fastmath @inbounds for e in 1:epochs
		tic()

	    iter += 1
		# determine the minibatch to operate on
		end_pos = svi_counter + mini_batch_size - 1
		if end_pos > m.D
			end_pos = m.D
		end
		mb_idx = mb_idx_rand[svi_counter:end_pos]

		# keep track of documents seen
		cur_count += end_pos - svi_counter + 1

		# when reading the end of the training set, reshuffle the documents and start at beginning again
		svi_counter += mini_batch_size
		if svi_counter > m.D
			svi_counter = 1
			mb_idx_rand = m.training_doc_ids[randperm(m.D)]
			iter = 1
		end
		# token count in the minibatch, gives a more realistic estimate of the fraction N/|S| for the gradient (instead of fractions of document counts)
		mb_wordcount = sum([sum(collect(values(doc.value))) for doc in m.corpus.documents[mb_idx]])
		mult = all_tokens/mb_wordcount

		# compute the RM learning rate (using a common form)
	    lr = a*(b+iter)^-gamma
	    # keep track of learning rates
		push!(m.learning_rates, lr)

		#do local udpate step (termed "e-step")
		(online_bound, ss_tk, ss_tkx, words_seen, t_mb, _) = e_step(m, mb_idx)
		# multiply bound estimate by multiplier
		online_bound *= mult
		# estimate of document specific bound for whole corpus
		push!(e_step_likelihoods, online_bound)

		e_step_time = toq()
		e_step_time_agg += e_step_time

		# m-step, update global variables
		tic()
		# compute natural gradients according to sufficient statistics (gradient w.r.t. to canonical params is natural gradient w.r.t. natural params)
		m_step_bound = 0.
		# only look at words actually observed in the minibatch
		KnmKmmInvEff = m.KnmKmmInv[t_mb,:]

		Λ_diags = zeros(m.T, m.V)
		means = zeros(m.T, m.V)
		word_bounds = zeros(m.V)
		for k in 1:m.K
			ɸ_k = mult.*ss_tk[k][t_mb] #(T x 1)
			Ξ_k = mult.*ss_tkx[k][t_mb, :] #(T x V)

			Threads.@threads for w in 1:m.V
				Λ_diags[t_mb,w] = vec(sum(KnmKmmInvEff*m.s[k,w] .* KnmKmmInvEff, 2))
			end

			means[t_mb, :] = KnmKmmInvEff*m.μ[k]

			B_tilde_k = ɸ_k .* exp(means[t_mb,:] .+ .5(Λ_diags[t_mb,:] .+ m.K_tilde_diag[t_mb]) .- m.zeta[k,t_mb])
			# euclidean gradient for mean
			dL_dm = KnmKmmInvEff' * (Ξ_k + B_tilde_k .* (means[t_mb,:] .- 1))

			m.eta_1[k] = (1-lr)*m.eta_1[k] + lr.*dL_dm
			Threads.@threads for w in 1:m.V
				# euclidean gradient for variance
				dL_dS = -.5*(m.KmmInv + (KnmKmmInvEff .* B_tilde_k[:,w])'*KnmKmmInvEff)

				# update step for second natural parameter
				m.eta_2[k,w] = (1-lr)*m.eta_2[k,w] + lr.*dL_dS

				# compute inverse and determinant using cholesky decomposition if possible
				eta_inv, det_eta_inv = chol_inv(m.eta_2[k,w] + m.jitter)
				# compute new covariance matrix and determinant (inducing points)
				m.s[k,w] = -.5*eta_inv
				det_s = (-.5)^num_inducing * det_eta_inv

				# compute new mean (inducing points)
				m.μ[k][:,w] = m.s[k,w]*m.eta_1[k][:,w]

				Λ_diags[:,w] = vec(sum(m.KnmKmmInv*m.s[k,w] .* m.KnmKmmInv,2))
				p_u = -.5*sum(m.KmmInv.*m.s[k,w])
				q_u = -.5*log(abs(det_s))
				word_bounds[w] = p_u-q_u
			end
			# compute means for all timestamps (needed to recompute zeta)
			means = m.KnmKmmInv*m.μ[k]
			m_step_bound += -.5*sum(m.KmmInv*m.μ[k] .* m.μ[k])
			# update zetas again
			svi_update_zeta(m, k, means, Λ_diags)

			m_step_bound += sum(word_bounds)
		end

		push!(m_step_likelihoods, m_step_bound)
		online_bound += m_step_bound
		m_step_time = toq()
		m_step_time_agg += m_step_time
		info("epoch $e, current elbo: $online_bound, e-step: $e_step_time, m-step: $m_step_time")

		# keep track of likelihoods
		push!(m.likelihoods, online_bound)
		push!(m.likelihood_counts, cur_count)

		# test on unseen data according to schedule
		if length(m.test_doc_ids) > 0 && e % test_schedule == 0
			test_ppx = test(m)
			push!(m.test_counts, cur_count)
			push!(m.test_perplexities, test_ppx)
			if test_schedule < 256 test_schedule *= 2	end
		end
	end
	if m.visualize
		display(Plots.plot(m.learning_rates, label="learning rates"))
		display(Plots.plot(m.likelihood_counts, e_step_likelihoods, label="e-steps"))
		display(Plots.plot(m.likelihood_counts, m_step_likelihoods, label="m-steps"))
		display(Plots.plot(m.test_counts, m.test_perplexities, title="Out-of-Sample Perplexities"))
	end
	time = toq()

	info("time for GP SVI: $time seconds, average e-step time: $(e_step_time_agg/epochs), average m-step time: $(m_step_time_agg/epochs)")
	info("e-step full: $e_step_time_agg, m-step full: $m_step_time_agg")
	time
end

function svi_update_zeta(m::GDTM, k::Int64, means::Matrix{Float64}, Λ_diags::Matrix{Float64})
	Threads.@threads for t in 1:m.T
		m.zeta[k,t] = log_sum(means[t,:] + .5*(Λ_diags[t,:] .+ m.K_tilde[t,t]))
	end
end

function svi_update_zeta(m::GDTM)
	for k in 1:m.K
		svi_update_zeta(m, k)
	end
end


text:
# ---
# title: 766. Toeplitz Matrix
# id: problem766
# author: Qling
# date: 2021-02-22
# difficulty: Easy
# categories: Array
# link: <https://leetcode.com/problems/toeplitz-matrix/description/>
# hidden: true
# ---
#
# Given an `m x n` `matrix`, return  _`true` if the matrix is Toeplitz.
# Otherwise, return `false`._
#
# A matrix is **Toeplitz** if every diagonal from top-left to bottom-right has
# the same elements.
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2020/11/04/ex1.jpg)
#
#
#
#     Input: matrix = [[1,2,3,4],[5,1,2,3],[9,5,1,2]]
#     Output: true
#     Explanation:
#     In the above grid, the diagonals are:
#     "[9]", "[5, 5]", "[1, 1, 1]", "[2, 2, 2]", "[3, 3]", "[4]".
#     In each diagonal all elements are the same, so the answer is True.
#
#
# **Example 2:**
#
# ![](https://assets.leetcode.com/uploads/2020/11/04/ex2.jpg)
#
#
#
#     Input: matrix = [[1,2],[2,2]]
#     Output: false
#     Explanation:
#     The diagonal "[1, 2]" has different elements.
#
#
#
#
# **Constraints:**
#
#   * `m == matrix.length`
#   * `n == matrix[i].length`
#   * `1 <= m, n <= 20`
#   * `0 <= matrix[i][j] <= 99`
#
#
#
# **Follow up:**
#
#   * What if the `matrix` is stored on disk, and the memory is limited such that you can only load at most one row of the matrix into the memory at once?
#   * What if the `matrix` is so large that you can only load up a partial row into the memory at once?
#
#
## @lc code=start
using LeetCode

function is_toeplitz_matrix(matrix::Vector{Vector{Int}})
    m, n = length(matrix), length(matrix[1])

    for i in 1:(m - 1)
        (matrix[i][1:(end - 1)] != matrix[i + 1][2:end]) && return false
    end

    return true
end

## @lc code=end


text:
module JuliaFEM

type Beam
    L::Float64
    rhoA::Float64
    EI::Float64
    Beam(l, ra, ei) = new(l, ra, ei)
end

function mass_matrix(beam::Beam)
    a = beam.L / 2
    beam.rhoA * a / 105 * [78. 22  27 -13;
                           22   8  13  -6;
                           27  13  78 -22;
                          -13  -6 -22   8]
end

function stiff_matrix(beam::Beam)
    a = beam.L / 2

    beam.EI/(2*a^3) * [3.  3a  -3  3a;
                       3a  4a^2 -3a 2a^2;
                       -3  -3a  3 -3a;
                       3a 2a^2 -3a 4a^2]
end

function rhs_function(beam, p)
    a = beam.L / 2
    [p*a, p*a*a/3, p*a, -p*a*a/3]
end

type Beam2d
    L::Float64
    rhoA::Float64
    EI::Float64
    EA::Float64
    Beam(l, ra, ei, ea) = new(l, ra, ei, ea)
end

function mass_matrix(b::Beam2d)
    a = b.L/2
    b.rho*b.A*a/105.0 * [70.0  0   0  35  0  0;
                         0    78 22a   0 27 -13a;
                         0  22a 8a^2 0 13a -6a^2;
                         35 0     0   70  0  0;
                         0  27  13a  0  78 -22a;
                         0 -13a -6a^2 0 -22a 8a^2]
end

function stiff_matrix(b::Beam2d)
    a = b.L/2
    r2 = b.EI / b.EA

    f = a^2 / r2

    b.E*b.I/(2a^3) * [ f  0  0 -f  0  0;
                       0  3 3a  0 -3 3a;
                       0 3a 4a^2 0 -3a 2a^2;
                      -f  0  0  f  0   0;
                       0 -3  -3a  0  3 -3a;
                       0 3a 2a^2 0 -3a 4a^2]
end

end


text:
include("../src/sir_small.jl")
using PyPlot

function no_of_particles()
	K = 50
	σ_o = 0.1
	σ_d = 0.0
	Δ = 5
	Nth = 500
	obsfun = obs1

	Npts = 20
	Np = 100

	trj_len = K*Δ
	mean_rmse = zeros(Npts)
	Np_arr = zeros(Int64, Npts)

	for i = 1:Npts
	   @show Np
	   x_trj, w_trj, y, x_true = assimilate(K, Np,
							σ_o, σ_d,
							Δ, Nth,
							obsfun)
	    x_trj = permutedims(x_trj,[3,2,1])

	    x_true = permutedims(x_true,[2,1])

        filter_mean = reshape(sum(x_trj,dims=2)/Np, trj_len, d)
        norm2(x,dim) = sum(x -> x^2, x, dims=dim)
        rmse = norm2(x_true .- filter_mean, 2)
	    rmse .= sqrt.(rmse)/d
	    mean_rmse[i] = sum(rmse[1:Δ:end])/K
		Np_arr[i] = Np
		Np = Np + 500
	end
	save("../data/mean_rmse_vs_Np.jld", "Np", Np_arr, "mean_rmse", mean_rmse)
end


text:
#generate community
p = random_micrm_params(5,10,0.5)

#check dimensions
@testset "random_micrm_params dimensions" begin
    @test p[:N] == 5
    @test p[:M] == 10
    @test size(p[:u]) == (5,10)
    @test size(p[:m]) == (5,)
    @test size(p[:ρ]) == (10,)
    @test size(p[:ω]) == (10,)
    @test size(p[:l]) == (10,10)
end

@testset "random_micrm_params values" begin
    @test p[:N] == 5
    @test p[:M] == 10
    @test all(sum(p[:l], dims = 2) .≈ 0.5)
    @test all(sum(p[:u], dims = 2) .≈ 1.0)
end


text:
"""
    rotation matrix: z-axis
"""
function rotation_matrix(x)
    [cos(x) -sin(x); sin(x) cos(x)]
end

"""
    rotation matrix: vector to vector
"""
function vector_rotation_matrix(ns; nw=[0.0; 1.0])
	a = atan(nw[2], nw[1]) - atan(ns[2], ns[1])
    rotation_matrix(a)
end

"""
    scientific notation
"""
function scn(a::Number; digits::Int=1, exp_digits::Int=1)
	typeof(a) <: Float64 ? nothing : return nothing
end

function scn(a::Float64; digits::Int=1, exp_digits::Int=1)
	@assert digits >= 0
    # a = m x 10^e
    if a == 0
        e = 0
        m = 0.0
    else
        e = Int(floor(log(abs(a))/log(10)))
        m = a*exp(-e*log(10))
    end

    m = round(m, digits=digits)
	if m == 10.0
		m = 1.0
		e += 1
	end
    if digits == 0
        m = Int(floor(m))
		strm = string(m)
	else
		strm = string(m)
		is_neg = m < 0.
		strm = strm*"0"^max(0, 2+digits+is_neg-length(strm))
    end
    sgn = a >= 0 ? " " : ""
    sgne = e >= 0 ? "+" : "-"

	stre = string(abs(e))
	stre = "0"^max(0, exp_digits - length(stre)) * stre
    return "$sgn$(strm)e$sgne$(stre)"
end


text:
# convert_a3d

# Method that allows federation by setting the `output_format`  in the Stanmodel().

"""

# convert_a3d

Convert the output file created by cmdstan to the shape of choice.

### Method
```julia
convert_a3d(a3d_array, cnames, ::Val{Symbol}; start=1)
```
### Required arguments
```julia
* `a3d_array::Array{Float64, 3},`      : Read in from output files created by cmdstan
* `cnames::Vector{AbstractString}`     : Monitored variable names
```

### Optional arguments
```julia
* `::Val{Symbol}`                      : Output format
* `::start=1`                          : First draw for MCMCChains.Chains
```
Method called is based on the output_format defined in the stanmodel, e.g.:

   stanmodel = Stanmodel(`num_samples`=1200, thin=2, name="bernoulli",
     model=bernoullimodel, `output_format`=:mcmcchains);

Current formats supported are:

1. :array (a3d_array format, the default for CmdStan)
2. :namedarray (NamedArrays object)
3. :dataframe (DataFrames object)
4. :mambachains (Mamba.Chains object)
5. :mcmcchains (TuringLang/MCMCChains.Chains object)

Options 3 through 5 are respectively provided by the packages StanDataFrames,
StanMamba, StanMCMCChains and StanMCMCChains.
```

### Return values
```julia
* `res`                       : Draws converted to the specified format.
```
"""
convert_a3d(a3d_array, cnames, ::Val{:array}; start=1) = a3d_array

convert_a3d(a3d_array, cnames, ::Val{:namedarray}; start=1) =
  [NamedArray(a3d_array[:,:,i], (collect(1:size(a3d_array, 1)), Symbol.(cnames)))
    for i in 1:size(a3d_array, 3)]

function convert_a3d(a3d_array, cnames, ::Val{:mcmcchains}; start=1)
  pi = filter(p -> length(p) > 2 && p[end-1:end] == "__", cnames)
  p = filter(p -> !(p in  pi), cnames)

  MCMCChains.Chains(a3d_array,
    cnames,
    Dict(
      :parameters => p,
      :internals => pi
    );
    start=start
  )
end


text:
using FTCTests
using JLD2
using UnPack
using DataFrames
using Statistics
using Transducers

"""
# Notes
- `simulation_height_success`: once height becomes lower than a prescribed threshold, `false`. Otherwise, `true`.
- `simulation_actual_time_success`: If a simulation time (real time) takes too long, `false`. Otherwise, `true`.
"""
function _file_path_to_nt(file_path::String, t1, threshold)
    jld2 = JLD2.load(file_path)
    @unpack faults, fdi, simulation_height_success, simulation_actual_time_success = jld2
    if simulation_actual_time_success == true
        if simulation_height_success == false
            is_success = false
        else
            is_success = evaluate(jld2, t1, threshold)
        end
    else
        is_success = missing  # not considered as either true or false
    end
    _nt = (;
           fault=extract_fault_property(faults),
           fault_kind=get_fault_kind(faults),
           fdi_delay=fdi.τ,
           is_success=is_success,
           sim_time_success=simulation_actual_time_success,
          )
end

function compute_recovery_rate(; _dir_log="data", t1=15.0, threshold=0.1)
    df = DataFrame()
    for manoeuvre in [:hovering, :forward]
        for method in [:adaptive, :adaptive2optim]
            _df = DataFrame()
            dir_log = joinpath(joinpath(_dir_log, String(manoeuvre)), String(method))
            file_paths = readdir(dir_log; join=true)
            _ = file_paths |> Map(file_path ->
                                  push!(_df, _file_path_to_nt(file_path, t1, threshold); promote=true)
                                 ) |> collect
            _df_simtime = filter(:sim_time_success => t -> t == true, _df)
            # group by FDI delay
            groups_delay = groupby(_df_simtime, :fdi_delay)
            for df_delay in groups_delay
                # group by fault kind
                groups_fault_kind = groupby(df_delay, :fault_kind)
                for df_fault_kind in groups_fault_kind
                    recovery_rate = mean(df_fault_kind.is_success)
                    push!(df, (;
                               manoeuvre=manoeuvre,
                               method=method,
                               delay=df_delay.fdi_delay[1],
                               fault_kind=df_fault_kind.fault_kind[1],  # the same fault kind
                               recovery_rate=recovery_rate,
                               num=length(df_fault_kind.is_success)
                              ); promote=true)
                end
            end
        end
    end
    df
end


text:
using JSON

scheme_prob = JSON.Parser.parsefile(".\\src\\scheme.json")

dielectrics = JSON.Parser.parsefile(".\\src\\dielectrics.json")

scheme_prob[1]["simul_params"]["freq"]

scheme_prob["waveguide"]["dimension"]["a"]


text:
using Documenter
using GroupTheoreticML

makedocs(
    modules = [GroupTheoreticML], doctest = false
)

# Documenter can also automatically deploy documentation to gh-pages.
# See "Hosting Documentation" and deploydocs() in the Documenter manual
# for more information.
deploydocs(
    deps = Deps.pip("pygments", "mkdocs", "python-markdown-math"),
    repo = "github.com/nchaudhr/GroupTheoreticML.jl.git",
    julia  = "0.5",
    osname = "linux"
    # ...
)

max_stars_repo_name: nchaudhr/GroupTheoreticML

max_issues_repo_name: nchaudhr/GroupTheoreticML

max_forks_repo_name: nchaudhr/GroupTheoreticML


text:
#include("imports.jl")
using FastAI
import FastAI: Block, Encoding, encode, decode, checkblock, encodedblock, decodedblock
using FastAI: Label, LabelMulti, Mask, Image, ImageTensor, testencoding
using FastAI: OneHot
using Test
using StaticArrays
using Images
using FastAI: grabbounds
using Images

##
struct ABlock <: Block
    end
checkblock(::ABlock, ::Int) = true

struct BBlock <: Block
end
checkblock(::BBlock, ::String) = true

struct AtoB <: Encoding end

encode(::AtoB, _, ::ABlock, data) = string(data)
decode(::AtoB, _, ::BBlock, data) = parse(Int, data)
encodedblock(::AtoB, ::ABlock) = BBlock()
decodedblock(::AtoB, ::BBlock) = ABlock()

Test.@testset "Encoding API" begin

    enc = AtoB()
    testencoding(enc, ABlock(), 100)
    testencoding(enc, (ABlock(), ABlock()), (100, 100))
    testencoding((enc,), ABlock(), 100)

    Test.@testset "tuple of encodings" begin
        @test encodedblock((enc, enc,), ABlock()) isa BBlock
        @test encodedblock((enc, enc,), BBlock()) isa Nothing
        @test decodedblock((enc, enc,), BBlock()) isa ABlock
        @test decodedblock((enc, enc,), ABlock()) isa Nothing
    end
end

@testset "OneHot" begin
    enc = OneHot()
    testencoding(enc, Label(1:10), 1)
    testencoding(enc, LabelMulti(1:10), [1])
    testencoding(enc, Mask{2}(1:10), rand(1:10, 50, 50))
end

@testset "Composition" begin
    encodings = (ProjectiveTransforms((5, 5)), ImagePreprocessing(), OneHot())
    blocks = (Image{2}(), Label(1:10))
    data = (rand(RGB{N0f8}, 10, 10), 7)
    testencoding(encodings, blocks, data)
end

@testset "ProjectiveTransforms" begin
    enc = ProjectiveTransforms((32, 32), buffered = false)
    block = Image{2}()
    image = rand(RGB, 100, 50)

    testencoding(enc, block, image)
    @testset "randstate is shared" begin
        im1, im2 = encode(enc, Training(), (block, block), (image, image))
        im3 = encode(enc, Training(), block, image)
        @test im1 ≈ im2
        @test !(im1 == im3)
    end

    @testset "don't transform data that doesn't need to be resized" begin
        imagesmall = rand(RGB, 32, 32)
        @test imagesmall ≈ encode(enc, Validation(), block, imagesmall)
    end

    @testset "3D" begin

        testencoding(ProjectiveTransforms((16, 16, 16)), Image{3}(), rand(RGB{N0f8}, 32, 24, 24))
    end
end


text:
"""
    checkpoint(f, args::Tuple)

Treat `f` as a primitive, whose gradient is computed by running and forwards- and reverse-
pass of ∇ inside the outer pass. This can have the effect of significantly lowering the
peak memory usage of a gradient computation at the expense of additional computation.

`f` musn't have any fields. In particular, this precludes `f` from being a closure.
"""
@generated function checkpoint(f, args::Tuple)
    fieldcount(f) > 0 && error("f mustn't have fields.")
    return :(f(args...))
end
@generated function checkpoint(f, args::Tuple, kwargs::NamedTuple)
    fieldcount(f) > 0 && error("f mustn't have fields. Can not checkpoint a closure or functor. Please reformulate as a plain function that does not close over any parameters.")
    return :(f(args...; kwargs...))
end

@explicit_intercepts checkpoint Tuple{Any, Tuple} [false, true]
@explicit_intercepts checkpoint Tuple{Any, Tuple, NamedTuple} [false, true, false]

function ∇(::typeof(checkpoint), ::Type{Arg{2}}, p, y, ȳ, f, args::Tuple)
    args_ = Leaf.(Tape(), args)
    y = f(args_...)
    if y isa Node
        ∇f = ∇(y, ȳ)
        ∇args = map(args_, args) do arg_, arg
            isassigned(∇f, arg_) ? ∇f[arg_] : zero(arg)
        end
    else
        ∇args = zero.(args)
    end
    return ∇args
end

function ∇(
    ::typeof(checkpoint),
    ::Type{Arg{2}},
    p, y, ȳ, f,
    args::Tuple,
    kwargs::NamedTuple,
)
    args_ = Leaf.(Tape(), args)
    y = f(args_...; kwargs...)
    if y isa Node
        ∇f = ∇(y, ȳ)
        ∇args = map(args_, args) do arg_, arg
            isassigned(∇f, arg_) ? ∇f[arg_] : zero(arg)
        end
    else
        ∇args = zero.(args)
    end
    return ∇args
end


text:
module Server

using Sockets
using JSON

using ..Simulation

mutable struct SockServer
    socket::Sockets.TCPSocket
    server::Sockets.TCPServer

    # The channel is a bridge between the socket and the Simulation.
    # For example, a message from the client to the Simulation would flow as:
    #
    # |----- client app --------|    |---- server app----------|
    # client -> channel -> socket -> socket -> channel -> server
    chan::Channel{String}

    running::Bool

    function SockServer()
        o = new()
        o.chan = Channel{String}(10)
        o.running = false
        o
    end
end

include("handlers.jl")

const PORT = 2001

function start(soc::SockServer)
    soc.server = listen(PORT)

    # Start task that scans the socket for incomming JSON messages
    # **Socket Task**
    @async while true
        soc.socket = accept(soc.server)
        open = true

        while open
            open = isopen(soc.socket)
            if open
                data = readline(soc.socket)
                msg = String(data)
                # Immediately send it to the channel
                # println("Server putting msg on chan: ", msg)
                put!(soc.chan, msg)
            else
                println("Socket closed!")
            end
        end
    end

    soc.running = true
    println("Listening on Channel...")

    # This is the server's main task that listens to the channel
    while soc.running
        # println("Taking...")
        if isready(soc.chan)
            json = take!(soc.chan)
            # println("Msg from channel: ", json)

            if json ≠ ""
                data = JSON.parse(json)

                handle_msg(soc, data)
            end
        else
            # The **Socket Task** above needs time to run so we yield.
            yield()
            # ~_-~_-~_-~_-~_-~_-~_-~_-~_-~_-~_-~_-~_-~_-~_-~_-~_-~_-
            # Just for simulating long process times.
            # Remove when done debugging and testing.
            sleep(0.1)
            # ~_-~_-~_-~_-~_-~_-~_-~_-~_-~_-~_-~_-~_-~_-~_-~_-~_-~_-
        end
    end

    println("Closing socket")
    Sockets.close(soc.server)

    println("Server is shutdown")
    println("Goodbye.")
end

end # Module -------------------------------------------


text:
import Oracle
using Dates

@assert VERSION >= v"0.7-"

@assert isfile(joinpath(@__DIR__, "..", "test","credentials.jl")) """
Before running tests, create a file `test/credentials.jl` with the content:

username = "your-username"
password = "your-password"
connect_string = "your-connect-string"
auth_mode = Oracle.ORA_MODE_AUTH_DEFAULT # or Oracle.ORA_MODE_AUTH_SYSDBA if user is SYSDBA
"""
include(joinpath(@__DIR__, "..", "test","credentials.jl"))

const NUM_CONNECTIONS = 10

let
    # JIT
    ctx = Oracle.Context()
    conn = Oracle.Connection(ctx, username, password, connect_string)
    pool = Oracle.Pool(ctx, username, password, connect_string)

    Oracle.close(conn)
    Oracle.close(pool)
end

println("Create $NUM_CONNECTIONS Connections")
@time let
    ctx = Oracle.Context()
    v = Vector{Oracle.Connection}()

    for i in 1:NUM_CONNECTIONS
        push!(v, Oracle.Connection(ctx, username, password, connect_string))
    end

    for c in v
        Oracle.close(c)
    end
end

let
    ctx = Oracle.Context()
    println("Time to create a Pool")
    @time pool = Oracle.Pool(ctx, username, password, connect_string, max_sessions=NUM_CONNECTIONS, session_increment=1)

    println("Time to acquire $NUM_CONNECTIONS connections from the pool 1st time")

    @time let
        v = Vector{Oracle.Connection}()

        for i in 1:NUM_CONNECTIONS
            push!(v, Oracle.Connection(pool))
        end

        for c in v
            Oracle.close(c)
        end
    end

    println("Time to acquire $NUM_CONNECTIONS connections from the pool 2nd time")

    @time let
        v = Vector{Oracle.Connection}()

        for i in 1:NUM_CONNECTIONS
            push!(v, Oracle.Connection(pool))
        end

        for c in v
            Oracle.close(c)
        end
    end

    Oracle.close(pool)
end


text:
#
# http://en.wikipedia.org/wiki/Allen%27s_interval_algebra
#

function before{T1, T2, N}(b1::HyperRectangle{N, T1}, b2::HyperRectangle{N, T2})
    for i = 1:N
        maximum(b1)[i] < minimum(b2)[i] || return false
    end
    true
end

@inline function meets{T1, T2, N}(b1::HyperRectangle{N, T1}, b2::HyperRectangle{N, T2})
    maximum(b1) == minimum(b2)
end

function overlaps{T1, T2, N}(b1::HyperRectangle{N, T1}, b2::HyperRectangle{N, T2})
    for i = 1:N
        maximum(b2)[i] > maximum(b1)[i] > minimum(b2)[i] &&
        minimum(b1)[i] < minimum(b2)[i] || return false
    end
    true
end

function starts{T1, T2, N}(b1::HyperRectangle{N, T1}, b2::HyperRectangle{N, T2})
    if minimum(b1) == minimum(b2)
        for i = 1:N
            maximum(b1)[i] < maximum(b2)[i] || return false
        end
        return true
    else
        return false
    end
end

function during{T1, T2, N}(b1::HyperRectangle{N, T1}, b2::HyperRectangle{N, T2})
    for i = 1:N
        maximum(b1)[i] < maximum(b2)[i] && minimum(b1)[i] > minimum(b2)[i] ||
        return false
    end
    true
end

function finishes{T1, T2, N}(b1::HyperRectangle{N, T1}, b2::HyperRectangle{N, T2})
    if maximum(b1) == maximum(b2)
        for i = 1:N
            minimum(b1)[i] > minimum(b2)[i] || return false
        end
        return true
    else
        return false
    end
end

#
# Containment
#

function isinside(rect::SimpleRectangle, x::Real, y::Real)
    rect.x <= x && rect.y <= y && rect.x + rect.w >= x && rect.y + rect.h >= y
end

function isinside(c::Circle, x::Real, y::Real)
    @inbounds ox,oy = origin(c)
    xD = abs(ox - x)
    yD = abs(oy - y)
    xD <= c.r && yD <= c.r
end

"""
contains(s::Simplex, pt; atol=0., rtol=defaulrtol(eltype(pt)))

Check if a point is contained inside a Simplex. If the intrinsic dimension
of the simplex is smaller then the dimension of the surrounding space
(for example a triangle in 3d) one needs the atol, rtol parameters.
"""
function contains(s::Simplex, pt; atol=0., rtol=Base.rtoldefault(eltype(pt)))
    w = weights(pt, s)
    all(w .>= -atol) || return false  # projection lies outside of s
    pt_proj = vertexmat(s) * w
    return isapprox(pt_proj, pt, atol=atol, rtol=rtol)
end

"""
Check if HyperRectangles are contained in each other. This does not use
strict inequality, so HyperRectangles may share faces and this will still
return true.
"""
function contains{T1, T2, N}(b1::HyperRectangle{N, T1}, b2::HyperRectangle{N, T2})
    for i = 1:N
        maximum(b2)[i] <= maximum(b1)[i] && minimum(b2)[i] >= minimum(b1)[i] ||
        return false
    end
    true
end

"""
Check if a point is contained in a HyperRectangle. This will return true if
the point is on a face of the HyperRectangle.
"""
function contains{T, N}(b1::HyperRectangle{N, T}, pt::Union{FixedVector, AbstractVector})
    for i = 1:N
        pt[i] <= maximum(b1)[i] && pt[i] >= minimum(b1)[i] || return false
    end
    true
end

"""
Check if HyperRectangles are contained in each other. This does not use
strict inequality, so HyperRectangles may share faces and this will still
return true.
"""
in(b1::HyperRectangle, b2::HyperRectangle) = contains(b2, b1)

"""
Check if a point is contained in a HyperRectangle. This will return true if
the point is on a face of the HyperRectangle.
"""
in(pt::Union{FixedVector, AbstractVector}, b1::HyperRectangle) = contains(b1, pt)

#
# Equality
#

==(a::AbstractMesh, b::AbstractMesh) = false
function =={M <: AbstractMesh}(a::M, b::M)
    for ((ka, va), (kb, vb)) in zip(all_attributes(a), all_attributes(b))
        va != vb && return false
    end
    true
end

(==){T1, T2, N}(b1::HyperRectangle{N, T1}, b2::HyperRectangle{N, T2}) =
    minimum(b1) == minimum(b2) && widths(b1) == widths(b2)

isequal(b1::HyperRectangle, b2::HyperRectangle) = b1 == b2

isless(a::SimpleRectangle, b::SimpleRectangle) = isless(area(a), area(b))


text:
Base.@kwdef struct DensityAnalysis{D, K, B}
    datalimits::D=automatic
    npoints::Int=200
    kernel::K=automatic
    bandwidth::B=automatic
end

# Work around lack of length 1 tuple method
# TODO: also add weights support here
_kde(data::NTuple{1, Any}; kwargs...) = kde(data...; kwargs...)
_kde(data::Tuple; kwargs...) = kde(data; kwargs...)

defaultdatalimits(positional) = map(v -> mapreduce(extrema, extend_extrema, v), Tuple(positional))

applydatalimits(f::Function, d) = map(f, d)
applydatalimits(limits::Tuple, _) = limits

function _density(vs::Tuple; datalimits, npoints, kwargs...)
    k = _kde(vs; kwargs...)
    intervals = applydatalimits(datalimits, vs)
    rgs = map(intervals) do (min, max)
        return range(min, max; length=npoints)
    end
    res = pdf(k, rgs...)
    return (rgs..., res)
end

function (d::DensityAnalysis)(input::ProcessedLayer)
    datalimits = d.datalimits === automatic ? defaultdatalimits(input.positional) : d.datalimits
    options = valid_options(; datalimits, d.npoints, d.kernel, d.bandwidth)
    output = map(input) do p, n
        return _density(Tuple(p); pairs(n)..., pairs(options)...), (;)
    end
    N = length(input.positional)
    labels = set(input.labels, N+1 => "pdf")
    plottypes = [LinesFill, Heatmap, Volume]
    default_plottype = plottypes[N]
    plottype = Makie.plottype(input.plottype, default_plottype)
    return ProcessedLayer(output; plottype, labels)
end

"""
    density(; datalimits=automatic, kernel=automatic, bandwidth=automatic, npoints=200)

Fit a kernel density estimation of `data`.
Here, `datalimits` specifies the range for which the density should be calculated,
and `kernel` and `bandwidth` are forwarded to `KernelDensity.kde`.
`npoints` is the number of points used by Makie to draw the line
"""
density(; options...) = transformation(DensityAnalysis(; options...))


text:
using Test
using Suppressor
using LogicCircuits

include("helper/plain_logic_circuits.jl")

@testset "Queries test" begin

    r1 = fully_factorized_circuit(PlainLogicCircuit,10)

    @test isdecomposable(r1)
    @test isdecomposable(compile(PlainLogicCircuit, Lit(1)))

    @test isdeterministic(r1)
    @test isdeterministic(compile(PlainLogicCircuit, Lit(1)))

    @test isstruct_decomposable(r1)
    @test isstruct_decomposable(compile(PlainLogicCircuit, Lit(1)))

    @test variables(r1) == BitSet(1:10)
    @test variables_by_node(r1)[r1] == BitSet(1:10)
    @test variables_by_node(r1)[children(children(r1)[1])[5]] == BitSet(5)

    @test num_variables(r1) == 10
    @test issmooth(r1)

    @test isone(sat_prob(r1))
    @test model_count(r1) == BigInt(2)^10

    @test r1(true, false, true, false, true, false, true, false, true, false)
    @test r1([true, false, true, false, true, false, true, false, true, false])

    #####################
    ors = map(1:10) do v
        pos = compile(PlainLogicCircuit, var2lit(Var(v)))
        neg = compile(PlainLogicCircuit, -var2lit(Var(v)))
        pos | neg
    end
    and1 = conjoin(ors[1], ors[2], ors[1], ors[3])
    and2 = conjoin(ors[5:10])
    or1 = and1 | and2
    @test !isdecomposable(or1)
    @test !isdeterministic(or1)
    @test !isstruct_decomposable(or1)

    #######################
    ors = map(1:10) do v
        pos = compile(PlainLogicCircuit, var2lit(Var(v)))
        neg = compile(PlainLogicCircuit, -var2lit(Var(v)))
        pos | neg
    end
    and1 = conjoin(ors[1:3])
    and2 = conjoin(ors[3:6])
    and3 = conjoin(ors[6:8])
    and4 = conjoin(ors[8:10])
    or1 = and1 | and2
    or2 = and3 | and4
    and5 = or1 & or2
    @test !isdecomposable(and5)
    @test !isdeterministic(and5)
    @test !isstruct_decomposable(and5)

    #######################
    leaf1 = compile(PlainLogicCircuit, Lit(1))
    leaf2 = compile(PlainLogicCircuit, Lit(-1))
    and = leaf1 & leaf2
    @test !isdecomposable(and)
    @test isdeterministic(and)
    @test !isstruct_decomposable(and)

    #######################
    n0c = little_3var_constants()
    @test isdecomposable(n0c)
    @test isdeterministic(n0c)
    @test isstruct_decomposable(n0c)
    @test !isdecomposable(n0c & little_2var())
    @test !isstruct_decomposable(n0c & little_2var())
    @test !issmooth(n0c | little_2var())
    @test !isdeterministic(n0c | little_2var())
    @test issmooth(n0c)
    @test @suppress_out !iscanonical(n0c & little_3var_constants(), 5, verbose=true)

    #######################
    lits = map(1:3) do var
        compile(PlainLogicCircuit, Lit(var))
    end

    or1 = lits[1] & lits[2] | lits[1] & lits[2]
    or2 = lits[1] & lits[3] | lits[1] & lits[3]
    and1 = or1 & lits[3]
    and2 = or2 & lits[2]
    circuit = and1 | and2
    @test isstruct_decomposable(or1)
    @test isstruct_decomposable(and2)
    @test isdecomposable(circuit)
    @test !isstruct_decomposable(circuit)
    @test !isstruct_decomposable(or1 & lits[3] | lits[3] & or1)
    @test isstruct_decomposable(or1 & lits[3] | or1 & lits[3])

end

@testset "Infer Vtrees Test" begin

    little = load_logic_circuit(zoo_psdd_file("little_4var.psdd"))
    little_vtree = zoo_vtree("little_4var.vtree")

    inferred = infer_vtree(little)::Vtree

    @test respects_vtree(little, little_vtree)
    @test respects_vtree(little, inferred)

    random = zoo_sdd("random.sdd")
    random_vtree = zoo_vtree("random.vtree")
    @test_throws String random_inferred = infer_vtree(random) # remove test_throws after "issues/#47" is fixed

    @test respects_vtree(random, random_vtree)
    # @test respects_vtree(random, random_inferred) # "issues/#47"

end


text:
const DEFINE = Dict{Symbol,Symbol}()

DEFINE[:nvmlInit] = :nvmlInit_v2
DEFINE[:nvmlDeviceGetPciInfo] = :nvmlDeviceGetPciInfo_v3
DEFINE[:nvmlDeviceGetCount] = :nvmlDeviceGetCount_v2
DEFINE[:nvmlDeviceGetHandleByIndex] = :nvmlDeviceGetHandleByIndex_v2
DEFINE[:nvmlDeviceGetHandleByPciBusId] = :nvmlDeviceGetHandleByPciBusId_v2
DEFINE[:nvmlDeviceGetNvLinkRemotePciInfo] = :nvmlDeviceGetNvLinkRemotePciInfo_v2


text:
# %%
include("CellularBase.jl")
using .CellularBase
println("using CB")
@time using Plots
println("using Plots")
@time using Images
println("using Images")
# %%
mutable struct LifeGrid <: AbstractGrid{Int}
    state::Array{Int}
    otherstate::Array{Int}
    bc::BoundaryCondition
    neighborhood::AbstractNeighborhood
    τi::Int64
    τr::Int64
    max::Float64
    min::Float64
    function LifeGrid(
            size::Tuple{Int,Int}, τi::Int, τr::Int;
            max::Float64, min::Float64, init=nothing)
        z = zeros(size)
        if isnothing(init) z[(size .÷ 2)...] = 1
        else z[rand(CartesianIndices(z), init)] .= 1
        end
        new(z, copy(z), FixedMax, VonNeumannNeighborhood(3, 2), τi, τr, max, min)
    end
end

CellularBase.possible_states(l::LifeGrid) = 0:(l.τi + l.τr)

CellularBase.newstate(grid::LifeGrid) = grid.otherstate;

function CellularBase.state!(grid::LifeGrid, newstate)
    grid.otherstate = grid.state
    grid.state = newstate
end
# %%
infect(fraction::Float64, g::LifeGrid) = rand() < max(((g.max - g.min) * fraction + g.min), 0)
susceptible(i::Int, g::LifeGrid) = i == 0
infectious(i::Int, g::LifeGrid) = i ∈ 1:g.τi
refactory(i::Int, g::LifeGrid) = i ∈ g.τi:(g.τr + g.τi)
recovered(i::Int, g::LifeGrid) = i == (g.τr + g.τi)
function (g::LifeGrid)(neighbors::Array{Int}; kwargs...)
	current = neighbors[1]
    if susceptible(current, g)
        fraction = count(i -> infectious(i, g), neighbors) / length(neighbors)
        if infect(fraction, g)
            return 1
        else
            return 0
        end
    else
        return (current + 1) % (g.τi + g.τr)
    end
end
function color(i, grid)
    if susceptible(i, grid)
        RGB(0.9, 0.9, 0.9)
    elseif infectious(i, grid)
        RGB(1., 0., 0.)
    else
        RGB(0., 0., 0.)
    end
end;
# %%
function lockdown(grid, step; step_dict, radii)
    push!(radii, radius(grid.neighborhood))
    if step in keys(step_dict)
        grid.neighborhood = VonNeumannNeighborhood(step_dict[step], 2);
    end
end
step_dict = Dict(50 => 2, 100 => 1, 200 => 2, 300 => 3)
radii = []
grid = LifeGrid((200, 200), 4, 9, min=0., max=0.5, init=5)
k = simulate!(grid, 400; postrunhook=lockdown, step_dict=step_dict, radii=radii);
# %%
begin
    gif = @gif for frame ∈ k
        plot(frame .|> p -> color(p, grid))
    end every 2
    display(gif)
end
# %%
@time begin
    p = plot(1:length(k), k .|> frame -> count(frame .|> p -> susceptible(p, grid)), color=:green, label="Sus")
    plot!(1:length(k), k .|> frame -> count(frame .|> p -> infectious(p, grid)), color=:red, label="Inf")
    plot!(1:length(k), k .|> frame -> count(frame .|> p -> refactory(p, grid)), color=:blue, label="Ref")
    q = plot(radii, label="Neighbourhood Radius", legend=:bottomright)
    plot(p, q, size=(600, 500), layout=(2, 1))
end
# %%


text:
@diffusion_process FavettoSamson{T,K} begin
    :dimensions
    process --> 2
    wiener --> 2

    :parameters
    (α, β, λ, μ, σ1, σ2) --> T
    dose --> K

    :additional
    diagdiff --> true
    constdiff --> true
end

DiffusionDefinition.b(t, x, P::FavettoSamson) = @SVector [
    P.α * P.dose(t) - (P.λ +P.β)*x[1] + P.μ*x[2],
    P.λ*x[1] - P.μ*x[2]
]

DiffusionDefinition.σ(t, x, P::FavettoSamson) = SDiagonal(P.σ1, P.σ2)


text:
using DSP, BenchmarkTools

function model(t, β::Array{Float64,1}, ρ::Float64)

    Rt1 = Array{Float64}(undef, length(t))
    Rt2 = Array{Float64}(undef, length(t))
	A::Float64 = - 1/β[2]
    B::Float64 = 2/β[2] + 1/β[2]

	Threads.@threads for n in eachindex(t)

   		Rt1[n] = -exp(-(ρ^2/(t[n])) - β[1]*t[n])
   		Rt1[n] = Rt1[n]/(t[n]^(5/2))

   		Rt2[n] = A*exp(-(A^2/(t[n]))) - B*exp(-(B^2/(t[n])))

	end

    Rt = replace!(Rt1.*Rt2, NaN => 0)
    return Rt./maximum(Rt)

 end

 function conv_Rt(t, β::Array{Float64,1},  ρ::Float64, It)

	Rt = Array{Float64}(undef, length(t))
	convRt = Array{Float64}(undef, 2*length(t) - 1)

	Rt = model(t, β, ρ)

	convRT = conv(It, Rt)
    convRt = convRt./maximum(convRt)

end

function _conv_kern_fft!(out, u, v, su, sv, outsize, nffts)
    upad = _zeropad(u, nffts)
    vpad = _zeropad(v, nffts)
    p! = plan_fft!(upad)
    p! * upad # Operates in place on upad
    p! * vpad
    upad .*= vpad
    ifft!(upad)
    copyto!(out,
            CartesianIndices(out),
            upad,
            CartesianIndices(UnitRange.(1, outsize)))
end

@inline function _zeropad!(
    padded::AbstractVector,
    u::AbstractVector,
    padded_axes = axes(padded),
    data_dest::Tuple = (first(padded_axes[1]),),
    data_region = CartesianIndices(u),
)
    datasize = length(data_region)
    # Use axes to accommodate arrays that do not start at index 1
    data_first_i = first(data_region)[1]
    dest_first_i = data_dest[1]
    copyto!(padded, dest_first_i, u, data_first_i, datasize)
    padded[first(padded_axes[1]):dest_first_i - 1] .= 0
    padded[dest_first_i + datasize : end] .= 0

    padded
end

@inline function _zeropad!(
    padded::AbstractArray,
    u::AbstractArray,
    padded_axes = axes(padded),
    data_dest::Tuple = first.(padded_axes),
    data_region = CartesianIndices(u),
)
    fill!(padded, zero(eltype(padded)))
    dest_axes = UnitRange.(data_dest, data_dest .+ size(data_region) .- 1)
    dest_region = CartesianIndices(dest_axes)
    copyto!(padded, dest_region, u, data_region)

    padded
end

IRFpad = zeros(2*length(IRF) -1)
copyto!(IRFpad, IRF)

p! = plan_fft!(IRFpad)
p! * IRFpad

su = size(u)
sv = size(v)
outsize = su .+ sv .-1
axesu = axes(u)
axesv = axes(v)
out_offsets = first.(axesu) .+ first.(axesv)
out_axes = UnitRange.(out_offsets, out_offsets .+outsize .- 1)
out = similar(u, out_axes)

nffts = nextfastfft(outsize)

padded = similar(u, nffts)
padded_axes = axes(padded)
data_region = CartesianIndices(u)

fill!(padded, zero(eltype(padded)))
data_dest = first.(padded_axes)
dest_axes = UnitRange.(data_dest, data_dest .+ size(data_region) .- 1)
dest_region = CartesianIndices(dest_axes)
copyto!(padded, dest_region, u, data_region)

p! = plan_fft!(convert(Array{Complex{Float64}}, upad))
p! * upad
p! * vpad
upad .*= vpad
ifft!(upad)

function pad_fft(u, v)

    su = size(u)
    sv = size(v)
    outsize = su .+ sv .-1
    nffts = nextfastfft(outsize)

    upad = similar(u, nffts)
    vpad = similar(v, nffts)
    fill!(upad, zero(eltype(upad)))
    fill!(vpad, zero(eltype(vpad)))

    copyto!(upad, u)
    copyto!(vpad, v)

    p! = plan_fft!(convert(Array{Complex{Float64}}, upad))
    upad = p! * upad

function testing(upad, vpad, p!)

    vpad = p! * vpad
    vpad .*= upad
    real(ifft!(vpad))
end

outsize = 2 * length(IRF) .- 1
nffts = nextfastfft(outsize)

upad = zeros(Float64, nffts)
vpad = zeros(Float64, nffts)
convpad = zeros(Float64, nffts)

copyto!(upad, IRF)
copyto!(vpad, RtDT)

pl = plan_rfft(upad)

function testing1(pl, nffts, upad, vpad)
upadRfft = zeros(Complex{Float64}, 1 .+ div(nffts, 2, RoundDown))
vpadRfft = zeros(Complex{Float64}, 1 .+ div(nffts, 2, RoundDown))
u_vpadRfft = zeros(Complex{Float64}, 1 .+ div(nffts, 2, RoundDown))

upadRfft = pl * upad
vpadRfft = pl * vpad

u_vpadRfft = upadRfft .* vpadRfft

irfft(u_vpadRfft, length(upad))
end

outsize = 2 * length(IRF) .- 1
nffts = nextfastfft(outsize)

upad = zeros(Float64, nffts)
vpad = zeros(Float64, nffts)
convpad = zeros(Float64, nffts)

copyto!(upad, IRF)
pl = plan_rfft(upad)

upadRfft = zeros(Complex{Float64}, 1 .+ div(length(upad), 2, RoundDown))
upadRfft = pl * upad

RtDT = Array{Float64}(undef, length(IRF))

vpadRfft = zeros(Complex{Float64}, 1 .+ div(length(upad), 2, RoundDown))
u_vpadRfft = zeros(Complex{Float64}, 1 .+ div(length(upad), 2, RoundDown))

#copyto!(vpad, RtDT)

function testing1(t, pl, upadRfft, vpad, RtDT, u_vpadRfft, vpadRfft)

    RtDT = DT_model(t, [0.1,10.], 1., 1., 1.)

    copyto!(vpad, RtDT)

    vpadRfft = pl * vpad

    u_vpadRfft = upadRfft .* vpadRfft

    irfft(u_vpadRfft, length(upad))
end

function conv_DT(t, IRF)

	RtDT = Array{Float64}(undef, length(t))
	convDT = Array{Float64}(undef, 2*length(t) - 1)

    RtDT = DT_model(t, [0.1,10.], 1., 1., 1.)

	convDT = conv(IRF, RtDT)
	convDT = convDT./maximum(convDT)
end

max_stars_repo_name: sammzer/LightPropagation

max_issues_repo_name: sammzer/LightPropagation

max_forks_repo_name: sammzer/LightPropagation


text:
using BinaryBuilder, Compat

# Collection of sources required to build HDF5.  Use the CMake download because
# it includes appropriate zlib and szip sources, letting us build them into libhdf5.
sources = [
    "https://support.hdfgroup.org/ftp/HDF5/current18/src/CMake-hdf5-1.8.20.tar.gz" =>
    "1edc45f5065597a72ecf0663786c07aa89ecb4cfbaf025496c399b5ab1b6085e"
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/CMake-hdf5-*
mkdir build
cd build
cmake -DCMAKE_INSTALL_PREFIX=$prefix \
      -DCMAKE_TOOLCHAIN_FILE=/opt/$target/$target.toolchain -DCMAKE_SYSTEM_PROCESSOR=$(uname -p) \
      -DBUILD_SHARED_LIBS=ON -DH5_ENABLE_STATIC_LIB=OFF \
      -DHDF5_BUILD_CPP_LIB=OFF -DHDF5_BUILD_FORTRAN=OFF -DHDF5_ENABLE_PARALLEL=OFF \
      -DBUILD_TESTING=OFF -DHDF5_BUILD_TOOLS=OFF -DHDF5_BUILD_EXAMPLES=OFF -DHDF5_BUILD_HL_LIB=OFF \
      -DHDF5_PACKAGE_EXTLIBS=ON -DHDF5_ALLOW_EXTERNAL_SUPPORT=TGZ -DTGZPATH=`pwd`/.. \
      -DHDF5_ENABLE_SZIP_SUPPORT=ON -DSZIP_TGZ_NAME="SZip.tar.gz" -DSZIP_PACKAGE_NAME=szip \
      -DHDF5_ENABLE_Z_LIB_SUPPORT=ON -DZLIB_TGZ_NAME="ZLib.tar.gz" -DZLIB_PACKAGE_NAME=zlib \
      ../hdf5-*
make || cp bin/libszip-static.a bin/libszip.a # work around apparent bug in their makefile
make && make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms() # build on all supported platforms

# The products that we will ensure are always built
products(prefix) = [
    LibraryProduct(prefix, "libhdf5-shared", :libhdf5),
]

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

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


text:
using Test
using MLStyle

include("match.jl")
include("pattern.jl")
include("adt.jl")
include("fn.jl")
include("typelevel.jl")

# WARNING: typelevel.jl must be at the end of the test for once `type level` feature
#          is activated in a module, it won't be able to disable in this module.

# TODO:
# include("data.jl")


text:
module TestIteration
    using Base.Test, DataFrames, Compat

    dv = NullableArray(Nullable{Int}[1, 2, Nullable()])
    dm = NullableArray([1 2; 3 4])
    dt = NullableArray(zeros(2, 2, 2))

    df = DataFrame(A = 1:2, B = 2:3)

    for row in eachrow(df)
        @test isa(row, DataFrameRow)
        @test isequal(row[:B]-row[:A], Nullable(1))

        # issue #683 (https://github.com/JuliaStats/DataFrames.jl/pull/683)
        @test typeof(collect(row)) == @compat Array{Tuple{Symbol, Any}, 1}
    end

    for col in eachcol(df)
        @test isa(col, @compat Tuple{Symbol, NullableVector})
    end

    @test isequal(map(x -> minimum(convert(Array, x)), eachrow(df)), Any[1,2])
    @test isequal(map(minimum, eachcol(df)), DataFrame(A = [1], B = [2]))

    row = DataFrameRow(df, 1)

    row[:A] = 100
    @test isequal(df[1, :A], Nullable(100))

    row[1] = 101
    @test isequal(df[1, :A], Nullable(101))

    df = DataFrame(A = 1:4, B = ["M", "F", "F", "M"])

    s1 = sub(df, 1:3)
    s1[2,:A] = 4
    @test isequal(df[2, :A], Nullable(4))
    @test isequal(sub(s1, 1:2), sub(df, 1:2))

    s2 = sub(df, 1:2:3)
    s2[2, :B] = "M"
    @test isequal(df[3, :B], Nullable("M"))
    @test isequal(sub(s2, 1:1:2), sub(df, [1,3]))

    # @test_fail for x in df; end # Raises an error
end


text:
# ------------ Operation state alias and APIs ------------ #

const OperationState = MLIR.API.MlirOperationState

create_operation_state(name::StringRef, l::Location) = MLIR.API.mlirOperationStateGet(name, l)
create_operation_state(name::String, l::Location) = create_operation_state(StringRef(name), l)

# Constructor.
OperationState(name::String, l::Location) = create_operation_state(name, l)

function Base.display(state::OperationState)
    println("  __________________________________\n")
    println("             OperationState\n")
    for f in fieldnames(OperationState)
        println(" $f : $(getfield(state, f))")
    end
    println("  __________________________________\n")
end

@doc(
"""
const OperationState = MLIR.API.MlirOperationState
""", OperationState)

# ------------ Operation and APIs ------------ #

struct Operation
    ref::MLIR.API.MlirOperation
end

unwrap(o::Operation) = o.ref

create_operation(state::OperationState) = Operation(MLIR.API.mlirOperationCreate(Ref(state)))
destroy!(op::Operation) = MLIR.API.mlirOperationDestroy(unwrap(op))
is_null(op::Operation) = MLIR.API.mlirOperationIsNull(unwrap(op))
Base.:(==)(op1::Operation, op2::Operation) = MLIR.API.mlirOperationEqual(unwrap(op1), unwrap(op2))
get_block(op::Operation) = MLIR.API.mlirOperationGetBlock(unwrap(op))
get_parent_operation(op::Operation) = MLIR.API.mlirOperationGetParentOperation(unwrap(op))
get_num_regions(op::Operation) = MLIR.API.mlirOperationGetNumRegions(unwrap(op))
get_region(op::Operation, pos::Int) = MLIR.API.mlirOperationGetRegion(unwrap(op), pos)
get_next_in_block(op::Operation) = MLIR.API.mlirOperationGetNextInBlock(unwrap(op))
get_num_operands(op::Operation) = MLIR.API.mlirOperationGetNumOperands(unwrap(op))
get_operand(op::Operation, pos::Int) = MLIR.API.mlirOperationGetOperand(unwrap(op), pos)
get_num_results(op::Operation) = MLIR.API.mlirOperationGetNumResults(unwrap(op))
get_result(op::Operation, pos::Int) = MLIR.API.mlirOperationGetResult(unwrap(op), pos)
verify(op::Operation) = MLIR.API.mlirOperationVerify(unwrap(op))
dump(op::Operation) = MLIR.API.mlirOperationDump(unwrap(op))

# Constructor.
@inline Operation(state::OperationState) = create_operation(state)
function Operation(name::String, l::Location,
                   num_results::Int, result_types::Vector{Type},
                   num_operands::Int, operands::Vector{Value},
                   num_regions::Int, regions::Vector,
                   num_successors::Int, successors::Vector,
                   num_attributes::Int, attributes::Vector{NamedAttribute})
    state = OperationState(name, l, num_results, result_types, num_operands, operands, num_regions, regions, num_successors, successors, num_attributes, attributes)
    Operation(state)
end

@doc(
"""
const Operation = MLIR.API.MlirOperation
""", Operation)


text:
using Vulkan
using BenchmarkTools

# Benchmarks:
# - Creation of minimalist create info structs:
#   - directly
#   - from the high-level struct passed as a literal
#   - via a conversion from an existing high-level struct
# - Creation of high-level structs

println("InstanceCreateInfo")
@btime _InstanceCreateInfo([], []; application_info = _ApplicationInfo(v"1.0", v"1.5", v"1.2"; application_name="MyApp", engine_name="MyEngine"))
@btime _InstanceCreateInfo(InstanceCreateInfo([], []; application_info = ApplicationInfo(v"1.0", v"1.5", v"1.2"; application_name="MyApp", engine_name="MyEngine")))
ci = InstanceCreateInfo([], []; application_info = ApplicationInfo(v"1.0", v"1.5", v"1.2"; application_name="MyApp", engine_name="MyEngine"))
@btime _InstanceCreateInfo($ci)
@btime InstanceCreateInfo([], []; application_info = ApplicationInfo(v"1.0", v"1.5", v"1.2"; application_name="MyApp", engine_name="MyEngine"))

println("DescriptorPoolCreateInfo")
@btime _DescriptorPoolCreateInfo(3, _DescriptorPoolSize.([DESCRIPTOR_TYPE_SAMPLER, DESCRIPTOR_TYPE_STORAGE_BUFFER], [4, 6]))
@btime _DescriptorPoolCreateInfo(DescriptorPoolCreateInfo(3, DescriptorPoolSize.([DESCRIPTOR_TYPE_SAMPLER, DESCRIPTOR_TYPE_STORAGE_BUFFER], [4, 6])))
ci = DescriptorPoolCreateInfo(3, DescriptorPoolSize.([DESCRIPTOR_TYPE_SAMPLER, DESCRIPTOR_TYPE_STORAGE_BUFFER], [4, 6]))
@btime _DescriptorPoolCreateInfo($ci)
@btime DescriptorPoolCreateInfo(3, DescriptorPoolSize.([DESCRIPTOR_TYPE_SAMPLER, DESCRIPTOR_TYPE_STORAGE_BUFFER], [4, 6]))

# println("DeviceCreateInfo")
# @btime _DeviceCreateInfo([_DeviceQueueCreateInfo(find_queue_family(pdevice, QUEUE_GRAPHICS_BIT & QUEUE_COMPUTE_BIT), [1.0])], [], [])
# @btime _DeviceCreateInfo(DeviceCreateInfo([DeviceQueueCreateInfo(find_queue_family(pdevice, QUEUE_GRAPHICS_BIT & QUEUE_COMPUTE_BIT), [1.0])], [], []))
# ci = DeviceCreateInfo([DeviceQueueCreateInfo(find_queue_family(pdevice, QUEUE_GRAPHICS_BIT & QUEUE_COMPUTE_BIT), [1.0])], [], [])
# @btime _DeviceCreateInfo($ci)
# @btime DeviceCreateInfo([DeviceQueueCreateInfo(find_queue_family(pdevice, QUEUE_GRAPHICS_BIT & QUEUE_COMPUTE_BIT), [1.0])], [], [])

# Initialization (heavy on compile times)
println("Struct initialization")
@time using Vulkan
t = time_ns()
f1 = Vk.initialize(_PhysicalDeviceFeatures2, _PhysicalDeviceVulkan12Features, _PhysicalDeviceVulkanMemoryModelFeatures);
# f1 = Vk.initialize(_PhysicalDeviceFeatures2, _PhysicalDeviceVulkan12Features, _PhysicalDeviceVulkanMemoryModelFeatures, _PhysicalDeviceFragmentDensityMapFeaturesEXT);
println((time_ns() - t) / 1e9)
@btime Vk.initialize(_PhysicalDeviceFeatures2, _PhysicalDeviceVulkan12Features, _PhysicalDeviceVulkanMemoryModelFeatures);

max_stars_repo_name: JuliaGPU/VulkanAbstraction

max_issues_repo_name: JuliaGPU/VulkanAbstraction

max_forks_repo_name: JuliaGPU/VulkanAbstraction


text:
function fact(acc::BigInt,n)
    (n == 0) ? acc : fact(acc*n,n-1)
end

println(fact(BigInt(1),60000))