text:
# Listen, Attend and Spell: arxiv.org/abs/1508.01211

module ListenAttendSpell

using CuArrays
using Flux
using Flux: reset!, onecold, @functor, Recur, LSTMCell
using Zygote
using Zygote: Buffer, @adjoint
using LinearAlgebra
using JLD2
using FileIO
using IterTools
using Base.Iterators: reverse
using OMEinsum
using Logging
using TensorBoardLogger
using StatsBase

CuArrays.allowscalar(false)

include("utils.jl")

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

Constructs a bidirectional LSTM layer.
"""
struct BLSTM{M <: DenseMatrix, V <: DenseVector}
   forward  :: Recur{LSTMCell{M,V}}
   backward :: Recur{LSTMCell{M,V}}
   dim_out  :: Int
end

@functor BLSTM (forward, backward)

function BLSTM(in::Integer, out::Integer)
   forward  = LSTM(in, out)
   backward = LSTM(in, out)
   return BLSTM(forward, backward, out)
end

function BLSTM(forward::Recur{LSTMCell{M,V}}, backward::Recur{LSTMCell{M,V}}) where {M <: DenseMatrix, V <: DenseVector}
    size(forward.cell.Wi, 2) == size(backward.cell.Wi, 2) || throw(DimensionMismatch("input dimension, $(size(forward.cell.Wi, 2)), of the forward-time LSTM layer does not match the input dimension, $(size(backward.cell.Wi, 2)), of the backward-time LSTM layer"))

    out_dim = length(forward.cell.h)
    out_dim == length(backward.cell.h) || throw(DimensionMismatch("output dimension, $out_dim, of the forward-time LSTM layer does not match the output dimension, $(length(backward.cell.h)), of the backward-time LSTM layer"))
    return BLSTM(forward, backward, out_dim)
end

Base.show(io::IO, l::BLSTM) = print(io,  "BLSTM(", size(l.forward.cell.Wi, 2), ", ", l.dim_out, ")")

function flip(f, xs::V) where V <: AbstractVector
   rev_time = reverse(eachindex(xs))
   return getindex.(Ref(
      f.(getindex.(Ref(xs), rev_time))::V
   ), rev_time)
   # the same as
   # flipped_xs = Buffer(xs)
   # @inbounds for t ∈ rev_time
   #    flipped_xs[t] = f(xs[t])
   # end
   # return copy(flipped_xs)
   # but implemented via broadcasting as Zygote differentiates loops much slower than broadcasting
end

"""
    (m::BLSTM)(xs::DenseVector{<:DenseVecOrMat}) -> DenseVector{<:DenseVecOrMat}

Forward pass of the bidirectional LSTM layer for a vector of matrices input.
Input must be a vector of length T (time duration), whose each element is a matrix of size D×B (input dimension × # of batches).
"""
function (m::BLSTM)(xs::VM) where VM <: AbstractVector{<:DenseVecOrMat}
   vcat.(m.forward.(xs), flip(m.backward, xs))::VM
end

"""
    (m::BLSTM)(Xs::DenseArray{<:Real,3}) -> DenseArray{<:Real,3}

Forward pass of the bidirectional LSTM layer for a 3D tensor input.
Input tensor must be arranged in D×T×B (input dimension × time duration × # batches) order.
"""
function (m::BLSTM)(Xs::DenseArray{<:Real,3})
   # preallocate output buffer
   Ys = Buffer(Xs, 2m.dim_out, size(Xs,2), size(Xs,3))
   axisYs₁ = axes(Ys, 1)
   time    = axes(Ys, 2)
   rev_time = reverse(time)
   @inbounds begin
      # get forward and backward slice indices
      slice_f = axisYs₁[1:m.dim_out]
      slice_b = axisYs₁[(m.dim_out+1):end]
      # bidirectional run step
      setindex!.(Ref(Ys),  m.forward.(view.(Ref(Xs), :, time, :)), Ref(slice_f), time, :)
      setindex!.(Ref(Ys), m.backward.(view.(Ref(Xs), :, rev_time, :)), Ref(slice_b), rev_time, :)
      # the same as
      # @views for (t_f, t_b) ∈ zip(time, rev_time)
      #    Ys[slice_f, t_f, :] =  m.forward(Xs[:, t_f, :])
      #    Ys[slice_b, t_b, :] = m.backward(Xs[:, t_b, :])
      # end
      # but implemented via broadcasting as Zygote differentiates loops much slower than broadcasting
   end
   return copy(Ys)
end

# Flux.reset!(m::BLSTM) = reset!((m.forward, m.backward)) # not needed as taken care of by @functor

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

Constructs pyramid BLSTM layer, which is the same as the BLSTM layer, with the addition that the input is first concatenated at every two consecutive time steps before feeding it to the usual BLSTM layer.
"""
struct PBLSTM{M <: DenseMatrix, V <: DenseVector}
   forward  :: Recur{LSTMCell{M,V}}
   backward :: Recur{LSTMCell{M,V}}
   dim_out  :: Int
end

@functor PBLSTM (forward, backward)

function PBLSTM(in::Integer, out::Integer)
   forward  = LSTM(2in, out)
   backward = LSTM(2in, out)
   return PBLSTM(forward, backward, out)
end

function PBLSTM(forward::Recur{LSTMCell{M,V}}, backward::Recur{LSTMCell{M,V}}) where {M <: DenseMatrix, V <: DenseVector}
   size(forward.cell.Wi, 2) == size(backward.cell.Wi, 2) || throw(DimensionMismatch("input dimension, $(size(forward.cell.Wi, 2)), of the forward-time LSTM layer does not match the input dimension, $(size(backward.cell.Wi, 2)), of the backward-time LSTM layer"))

   out_dim = length(forward.cell.h)
   out_dim == length(backward.cell.h) || throw(DimensionMismatch("output dimension, $out_dim, of the forward-time LSTM layer does not match the output dimension, $(length(backward.cell.h)), of the backward-time LSTM layer"))
   return PBLSTM(forward, backward, out_dim)
end

Base.show(io::IO, l::PBLSTM) = print(io, "PBLSTM(", size(l.forward.cell.Wi, 2)÷2, ", ", l.dim_out, ")")

"""
    (m::PBLSTM)(xs::DenseVector{<:DenseVecOrMat}) -> DenseVector{<:DenseVecOrMat}

Forward pass of the pyramid BLSTM layer for a vector of matrices input.
Input must be a vector of length T (time duration), whose each element is a matrix of size D×B (input dimension × # of batches).
"""
function (m::PBLSTM)(xs::VM) where VM <: AbstractVector{<:DenseVecOrMat}
   # reduce time duration by half by restacking consecutive pairs of input along the time dimension
   evenidxs = (firstindex(xs)+1):2:lastindex(xs)
   x̄s = (i -> [xs[i-1]; xs[i]]).(evenidxs)::VM
   # counterintuitively the gradient of the following version is not much faster (on par in fact),
   # even though it is implemented via broadcasting
   # x̄s = vcat.(getindex.(Ref(xs), 1:2:lastindex(xs)), getindex.(Ref(xs), 2:2:lastindex(xs)))
   # x̄s = @views @inbounds(vcat.(xs[1:2:end], xs[2:2:end]))
   # x̄s = vcat.(xs[1:2:end], xs[2:2:end])
   # bidirectional run step
   return vcat.(m.forward.(x̄s), flip(m.backward, x̄s))::VM
end

"""
    (m::PBLSTM)(Xs::DenseArray{<:Real,3}) -> DenseArray{<:Real,3}

Forward pass of the pyramid BLSTM layer for a 3D tensor input.
Input tensor must be arranged in D×T×B (input dimension × time duration × # batches) order.
"""
function (m::PBLSTM)(Xs::DenseArray{<:Real,3})
   D, T, B = size(Xs)
   T½ = T÷2
   # reduce time duration by half by restacking consecutive pairs of input along the time dimension
   X̄s = reshape(Xs, 2D, T½, B)
   # preallocate output buffer
   Ys = Buffer(Xs, 2m.dim_out, T½, B)
   axisYs₁ = axes(Ys, 1)
   time    = axes(Ys, 2)
   rev_time = reverse(time)
   # get forward and backward slice indices
   slice_f = axisYs₁[1:m.dim_out]
   slice_b = axisYs₁[(m.dim_out+1):end]
   # bidirectional run step
   setindex!.(Ref(Ys),  m.forward.(view.(Ref(X̄s), :, time, :)), Ref(slice_f), time, :)
   setindex!.(Ref(Ys), m.backward.(view.(Ref(X̄s), :, rev_time, :)), Ref(slice_b), rev_time, :)
   # the same as
   # @views for (t_f, t_b) ∈ zip(time, reverse(time))
   #    Ys[slice_f, t_f, :] =  m.forward(X̄s[:, t_f, :])
   #    Ys[slice_b, t_b, :] = m.backward(X̄s[:, t_b, :])
   # end
   # but implemented via broadcasting as Zygote differentiates loops much slower than broadcasting
   return copy(Ys)
end

# Flux.reset!(m::PBLSTM) = reset!((m.forward, m.backward)) # not needed as taken care of by @functor

"""
    Encoder(dims::NamedTuple{(:blstm, :pblstms_out),Tuple{NamedTuple{(:in, :out),Tuple{Int,Int}}, T}}) where T <: Union{Int, NTuple{N,Int} where N} -> Chain

Construct Encoder neural network a.k.a Listener from dimensions specified by named tuple `dims`.
`Encoder` consists BLSTM layer followed by a block of PBLSTM layers. It accepts filter bank spectra as input and acts as an acoustic model encoder.

# Examples
```jldoctest
julia> dims = (
          blstm       = (in = 4, out = 8),
          pblstms_out = (12, 16, 20)
       );

julia> Encoder(dims)
Chain(BLSTM(4, 8), PBLSTM(16, 12), PBLSTM(48, 16), PBLSTM(64, 20))
```
"""
function Encoder(dims::NamedTuple{(:blstm, :pblstms_out),Tuple{NamedTuple{(:in, :out),Tuple{Int,Int}}, T}}) where T <: Union{Int, NTuple{N,Int} where N}
   (length(dims.pblstms_out) >= 1) || throw("Encoder must have at least 1 pyramid BLSTM layer")

   pblstm_layers = ( PBLSTM(2in, out) for (in, out) ∈ partition(dims.pblstms_out, 2, 1) )
   model = Chain(
      BLSTM(dims.blstm.in, dims.blstm.out),
      PBLSTM(2dims.blstm.out, first(dims.pblstms_out)),
      pblstm_layers...
   )
   return model
end

function MLP(in::Integer, out::NTuple{N,Integer}, σs::NTuple{N,Function}) where N
   model = Dense(in, first(out), first(σs))
   if N > 1
      layers = ( Dense(in, out, σ) for ((in, out), σ) ∈ zip(partition(out, 2, 1), σs[2:end]) )
      model = Chain(model, layers...)
   end
   return model
end

function MLP(in::Integer, out::Union{Integer, NTuple{N,Integer} where N}, σ::Function=identity)
   model = Dense(in, first(out), σ)
   if length(out) > 1
      layers = ( Dense(in, out, σ) for (in, out) ∈ partition(out, 2, 1) )
      model = Chain(model, layers...)
   end
   return model
end

function Decoder(in::Integer, out::Union{Integer, NTuple{N,Integer} where N})
   model = LSTM(in, first(out))
   if length(out) > 1
      layers = ( LSTM(in, out) for (in, out) ∈ partition(out, 2, 1) )
      model = Chain(model, layers...)
   end
   return model
end

CharacterDistribution(in::Integer, out::Integer) = Chain(Dense(in, out), logsoftmax)

struct LAS{M <: DenseMatrix, E, Aϕ, Aψ, D, C}
   context₀  :: M  # trainable initial value of the context
   decoding₀ :: M  # trainable initial value of the decoding
   listen    :: E  # encoder function
   key_ψ     :: Aψ # keys attention context function
   query_ϕ   :: Aϕ # query attention context function
   spell     :: D  # LSTM decoder
   infer     :: C  # character distribution inference function
end

@functor LAS

function LAS(encoder_dims::NamedTuple,
             attention_dim::Integer,
             decoder_dims::Tuple{Integer,Integer},
             output_dim::Integer)

   encoding_dim = 2last(encoder_dims.pblstms_out)
   decoding_dim =  last(decoder_dims)
   context₀  = zeros(Float32, encoding_dim, 1)
   decoding₀ = zeros(Float32, decoding_dim, 1)
   listen  = Encoder(encoder_dims)
   key_ψ   = MLP(encoding_dim, attention_dim)
   query_ϕ = MLP(decoding_dim, attention_dim)
   spell   = Decoder(encoding_dim + decoding_dim + output_dim, decoder_dims)
   infer   = CharacterDistribution(encoding_dim + decoding_dim, output_dim)

   LAS(context₀, decoding₀, listen, key_ψ, query_ϕ, spell, infer) |> gpu
end

function LAS(encoder_dims::NamedTuple,
             attention_dim::Integer,
             decoder_out_dim::Integer,
             output_dim::Integer)
   decoder_out_dim₁ = last(encoder_dims.pblstms_out) + decoder_out_dim + out_dim÷2
   decoder_dims = (decoder_out_dim₁, decoder_out_dim)
   LAS(encoder_dims, attention_dim, decoder_dims, output_dim)
end

function Base.show(io::IO, m::LAS)
   input_dim     = size(first(m.listen).forward.cell.Wi, 2)
   encoding_dim  = size(m.context₀, 1)
   attention_dim = length(m.key_ψ.b)
   decoding_dim  = size(m.decoding₀, 1)
   output_dim    = length(first(m.infer).b)
   print(io,
      """
      LAS(
         # (input_dim=$input_dim, encoding_dim=$encoding_dim, attention_dim=$attention_dim, decoding_dim=$decoding_dim, output_dim=$output_dim)
         $(m.listen),
         $(m.key_ψ),
         $(m.query_ϕ),
         $(m.spell),
         $(m.infer)
      )
      """
   )
end

# Flux.reset!(m::LAS) = reset!((m.listen, m.spell)) # not needed as taken care of by @functor

time_squashing_factor(m::LAS) = 2^(length(m.listen) - 1)

# 0. rewrite `@ein context[d,b] := αᵢs[t,b] * Hs[d,t,b]` as `@ein context[d,b] := Hs[d,t,b] * αᵢs[t,b]` to elide the call to permutedims and ensure that `context` did not change.
# 1. include actual one-hot prediction instead of current just output of the network
# 2. add teacher forcing
# 3. add shufling batch order after every epoch and withing batches during construction

@inline function decode(m::LAS{M}, Hs::T, maxT::Integer) where {M <: DenseMatrix, T <: DenseArray{<:Real,3}}
   prediction′₀ = Zygote.bufferfrom(zeros(eltype(Hs), length(first(m.infer).b), 1))
   prediction′₀[1] = 1.0f0
   prediction₀ = M(copy(prediction′₀))
   M′ = addparent(M, T)
   return _decode(m, Hs, maxT, prediction₀, M′)
end

@inline function _decode(m::LAS{M}, Hs::DenseArray{<:Real,3}, maxT::Integer, prediction::M, M′::DataType) where M <: DenseMatrix
   batch_size = size(Hs, 3)
   # precompute keys ψ(H) by gluing the slices of Hs along the batch dimension into a single D×TB matrix, then
   # passing it through the ψ dense layer in a single pass and then reshaping the result back into D′×T×B tensor
   ψHs = reshape(m.key_ψ(reshape(Hs, size(Hs,1), :)), size(m.key_ψ.W, 1), :, batch_size)
   # ψhs = m.key_ψ.(getindex.(Ref(Hs), :, axes(Hs,2), :))
   # check: all(ψhs .≈ eachslice(ψHs; dims=2))
   # compute inital decoder state for the entire batch
   decoding = m.spell(repeat([m.decoding₀; prediction; m.context₀]::M, 1, batch_size))::M
   # allocate D×B×T output tensor
   Ŷs = Buffer(Hs, size(prediction, 1), batch_size, maxT)
   @inbounds for t ∈ axes(Ŷs, 3)
      # compute query ϕ(sᵢ)
      ϕsᵢ = m.query_ϕ(decoding)
      # compute energies via batch matrix multiplication
      # @ein Eᵢs[t,b] := ϕsᵢ[d,b] * ψHs[d,t,b]
      Eᵢs = einsum(EinCode{((1,2), (1,3,2)), (3,2)}(), (ϕsᵢ, ψHs))::M′ # dispatches to batched_contract
      # check: Eᵢs ≈ reduce(hcat, diag.((ϕsᵢ',) .* ψhs))'
      # compute attentions weights
      αᵢs = softmax(Eᵢs)
      # compute attended context using Einstein summation convention, i.e. contextᵢ = Σᵤαᵢᵤhᵤ
      # @ein context[d,b] := αᵢs[t,b] * Hs[d,t,b]
      context = einsum(EinCode{((1,2), (3,1,2)), (3,2)}(), (αᵢs, Hs))::M′ # dispatches to batched_contract
      # check: context ≈ reduce(hcat, [sum(αᵢs[t,b] *Hs[:,t,b] for t ∈ axes(αᵢs, 1)) for b ∈ axes(αᵢs,2)])
      # predict probability distribution over character alphabet
      Ŷs[:,:,t] = prediction = m.infer([decoding; context]::M)
      # compute decoder state
      decoding = m.spell([decoding; prediction; context]::M)::M
   end
   return copy(Ŷs)
end

function (m::LAS)(xs::AbstractVector{<:DenseMatrix}, maxT::Integer)
   # compute input encoding, which are also values for the attention layer
   hs = m.listen(xs)
   dim_out, batch_size = size(first(hs))
   # transform T-length sequence of D×B matrices into the D×T×B tensor by first conconcatenating matrices
   # along the 1st dimension and to get singe DT×B matrix and then reshaping it into D×T×B tensor
   Hs = reshape(reduce(vcat, hs), dim_out, :, batch_size)
   # perform attend and spell steps
   Ŷs = decode(m, Hs, maxT)
   return Ŷs
end

function (m::LAS)(Xs::DenseArray{<:Real,3}, maxT::Integer)
   # compute input encoding, which are also values for the attention layer
   Hs = m.listen(Xs)
   # perform attend and spell steps
   Ŷs = decode(m, Hs, maxT)
   return Ŷs
end

function (m::LAS)(x::AbstractVector{<:DenseVector})
   T = length(x)
   X = reshape(reduce(hcat, pad(x, time_squashing_factor(m))), Val(3))
   return m(X, T)
end

# encoding_dim  = (512, 512, 512, 512)
# attention_dim = 512
# decoding_dim  = 512
# dim_feed_forward = 128
# dim_LSTM_speller = 512
# initialize with uniform(-0.1, 0.1)

function loss(m::LAS, X::DenseArray{<:Real,3}, indices::DenseVector, maxT::Integer)
   reset!(m)
   Ŷs = m(X, maxT)
   l = -sum(Ŷs[indices])
   return l
end

function loss(m::LAS, dataset::AbstractVector{<:Tuple{DenseArray{<:Real,3}, DenseVector, Integer}})
   return sum(dataset) do (X, indices, maxT)
      loss(m, X, indices, maxT)
   end
end

# best path decoding
function predict(m::LAS, xs::DenseVector{<:DenseMatrix{<:Real}}, lengths::DenseVector{<:Integer}, labels)::DenseVector{<:DenseVector}
   maxT = maximum(lengths)
   Ŷs = m(xs, maxT) |> cpu
   predictions = [onecold(@view(Ŷs[:, 1:len, n]), labels) for (n, len) ∈ enumerate(lengths)]
   return predictions
end

function predict(m::LAS, xs::DenseVector{<:DenseVector{<:Real}}, labels)::DenseVector
   Ŷ = m(xs) |> cpu
   prediction = onecold(Ŷ, labels)
   return prediction
end

function main(; n_epochs::Integer=1, saved_results::Bool=false, log_weights::Bool=true)
# load data & construct the neural net
las, phonemes, data_trn,
data_evl, length_evl,
data_val, length_val =
let batch_size = 88, valsetsize = 352
   JLD2.@load "data/TIMIT/TIMIT_MFCC/data_train.jld2" Xs ys PHONEMES
   output_dim = length(PHONEMES)

   if saved_results
      JLD2.@load "ListenAttendSpell.jl/models/TIMIT/las.jld2" las loss_val_saved
   else
      # encoder_dims = (
      #    blstm       = (in = (length ∘ first ∘ first)(Xs), out = 17),
      #    pblstms_out = (7, 5, 3)
      # )
      # attention_dim = 11
      # decoder_out_dim = 13
      encoder_dims = (
         blstm       = (in = (length ∘ first ∘ first)(Xs), out = 240),
         pblstms_out = (230, 230, 230)
      )
      attention_dim   = 230
      decoder_out_dim = 230
      las = LAS(encoder_dims, attention_dim, decoder_out_dim, output_dim)
   end

   multiplicity = time_squashing_factor(las)
   data_trn = batch_dataset(Xs, ys, output_dim, batch_size, multiplicity)

   idxs_evl = sample(eachindex(ys), valsetsize; replace=false)
   ys_evl = ys[idxs_evl]
   data_evl = batch_dataset(Xs[idxs_evl], ys_evl, output_dim, batch_size, multiplicity)
   length_evl = sum(length, ys_evl)

   JLD2.@load "data/TIMIT/TIMIT_MFCC/data_test.jld2" Xs ys
   ys_val = ys[1:valsetsize]
   data_val = batch_dataset(Xs[1:valsetsize], ys_val, output_dim, batch_size, multiplicity)
   length_val = sum(length, ys_val)

   las, PHONEMES, data_trn,
   data_evl, length_evl,
   data_val, length_val
end

# initialize TensorBoard logger
tblogger = TBLogger("log", tb_overwrite)
loss_val_prev = loss_val_saved = loss(las, data_val)
optimiser = ADAM(0.0001)

if log_weights
   logweights = () -> begin
      params_dict = param_dict(las, "las")
      @info "model" params = params_dict log_step_increment=0
   end
else
   logweights = () -> nothing
end

function callback()
   loss_evl = loss(las, data_evl)
   loss_val = loss(las, data_val)
   truthprob_evl  = truthprob(loss_evl, length_evl)
   truthprob_val  = truthprob(loss_val, length_val)
   println()
   @show loss_evl loss_val truthprob_evl truthprob_val
   println()
   if loss_val < loss_val_saved
      loss_val_saved = loss_val
      save("ListenAttendSpell.jl/models/TIMIT/las.jld2", Dict("las" => cpu(las), "loss_val_saved" => loss_val_saved))
      @info "Saved results!"
      println()
   end

   # halve_η = false
   # if loss_val > loss_val_prev
   #    # if the validation error increased from previous epoch
   #    # then halve the learning rate
   #    halve_η = true
   # else # ask user if he still wants to halve the learning rate
   #    println("Do you want to halve the current learning rate of η = $(optimiser.eta)? [yes/\e[4mno\e[0m]")
   #    ans = "\n"
   #    @async(ans = readline(stdin))
   #    timedwait(() -> ans != "\n", 10.0; pollint=0.5)
   #    (lowercase(ans) ∈ ("yes", "y")) && (halve_η = true)
   # end
   # if halve_η
   #    # halve the learning rate
   #    η = max(optimiser.eta/2, 1e-5)
   #    optimiser.eta = η
   #    println()
   #    @info "Adjusted learning rate to:" η
   # end

   with_logger(tblogger) do
      logweights()
      @info "train" loss=loss_evl truthprob=truthprob_evl log_step_increment=0
      @info "valid" loss=loss_val truthprob=truthprob_val
   end
   loss_val_prev = loss_val
end

θ = Flux.params(las)

nds = ndigits(length(data_trn))
callback()
# main training loop
for epoch ∈ 1:n_epochs
   println()
   @info "Epoch $epoch:"
   duration = @elapsed for (n, (X, indices, maxT)) ∈ enumerate(data_trn)
      l, pb = Flux.pullback(θ) do
         loss(las, X, indices, maxT)
      end
      print("batch ", ' '^(nds - ndigits(n)), n, ": ")
      @show l
      θ̄ = pb(one(l))
      Flux.Optimise.update!(optimiser, θ, θ̄)
   end
   println("\nCompleted in ", round(duration/60; sigdigits=2), " minutes")
   callback()
end
end

"""
    levendist(seq₁::AbstractVector, seq₂::AbstractVector)::Int
    levendist(seq₁::AbstractString, seq₂::AbstractString)::Int

Levenshtein distance between sequences `seq₁` and `seq₂`.
"""
function levendist(seq₁::AbstractVector, seq₂::AbstractVector)::Int
   # ensure that length(seq₁) <= length(seq₂)
   if length(seq₁) > length(seq₂)
      seq₁, seq₂ = seq₂, seq₁
   end
   # ignore prefix common to both sequences
   start = length(seq₁) + 1
   for (i, (el₁, el₂)) ∈ enumerate(zip(seq₁, seq₂))
      if el₁ != el₂
         start = i
         break
      end
   end
   @views begin
      seq₁, seq₂ = seq₁[start:end], seq₂[start:end]
      # ignore suffix common to both sequences
      lenseq₁ = length(seq₁)
      offset = lenseq₁
      for (i, el₁, el₂) ∈ zip(0:lenseq₁, reverse(seq₁), reverse(seq₂))
         if el₁ != el₂
            offset = i
            break
         end
      end
      seq₁, seq₂ = seq₁[1:(end-offset)], seq₂[1:(end-offset)]
   end
   lenseq₁ = length(seq₁)
   dist = length(seq₂)
   # if all of shorter sequence matches prefix and/or suffix of longer sequence, then Levenshtein
   # distance is just the delete cost of the additional characters present in longer sequence
   lenseq₁ == 0 && return dist

   costs = collect(eachindex(seq₂))
   @inbounds for (i, el₁) ∈ zip(0:(lenseq₁-1), seq₁)
      left = dist = i
      for (j, el₂) ∈ enumerate(seq₂)
         # cost on diagonal (substitution)
         above, dist, left = dist, left, costs[j]
         if el₁ != el₂
            # minimum of substitution, insertion and deletion costs
            dist = 1 + min(dist, left, above)
         end
         costs[j] = dist
      end
      # @show costs, dist
   end
   return dist
end

levendist(seq₁::AbstractString, seq₂::AbstractString)::Int = levendist(collect(seq₁), collect(seq₂))

per(source_phoneme, target_phoneme)::Real = levendist(source_phoneme, target_phoneme)/length(target_phoneme)
cer(source_chars, target_chars)::Real = levendist(source_chars, target_chars)/length(target_chars)
wer(source_words, target_words)::Real = levendist(source_words, target_words)/length(target_words)

end # ListenAttendSpell


text:
using Test, PolynomialBases

ufunc(x) = sinpi(x)

for basis_type in subtypes(PolynomialBases.NodalBasis{PolynomialBases.Line})
    for p in 1:5
        basis = basis_type(p)
        u = ufunc.(basis.nodes)
        D, M, R, B, MinvRtB = utility_matrices(basis)

        @test integrate(u, basis) ≈ sum(M * u)
        Ru = R*u
        @test Ru[1] ≈ interpolate(-1, u, basis)
        @test Ru[2] ≈ interpolate(+1, u, basis)

        @inferred includes_boundaries(basis)
        if includes_boundaries(basis) == Val{true}()
            @test Ru[1] ≈ u[1]
            @test Ru[end] ≈ u[end]
        else
            @test !(Ru[1] ≈ u[1])
            @test !(Ru[end] ≈ u[end])
        end

        # SBP property
        if satisfies_sbp(basis) == Val{true}()
            @test M*D + D'*M ≈ R'*B*R
        end
    end
end

ext: jl


text:
using Documenter
makedocs(;
    sitename="数値計算"
)


text:
function xwalkforward(x)
  2x
end


text:
using LazySets, MathematicalPredicates

function linearSwitching_specification()
    # initial set in mode 1
    X0_m1 = Singleton([3.1, 4.0, 0.0, 0.0, 0.0])
    X0 = [(1, X0_m1)]

    # safety property: x1 > -1.2
    border = HalfSpace([1.0, 0.0, 0.0, 0.0, 0.0], -1.2)  # x1 ≤ -1.2
    property = is_disjoint_from(border)

    # time horizon: 1 time unit
    time_horizon = 1.0

    # specification
    O = Dict(:T => time_horizon, :property => property)

    return X0, O
end

max_stars_repo_name: JuliaReach/ReachabilityBenchmarks

max_issues_repo_name: JuliaReach/ReachabilityBenchmarks

max_forks_repo_name: JuliaReach/ReachabilityBenchmarks


text:
"""
TODO: Provide description of crossover module.

RULE: Number of inputs must be greater or equal to the number of outputs.
"""
module crossover
importall common, variation
using core, utility, representation, individual, species
export Crossover, CrossoverDefinition, num_required

"""
Base type used by all crossover operation definitions.
"""
abstract CrossoverDefinition <: VariationDefinition

"""
Base type used all by crossover operators.

All crossover operators MUST implement a +stage+ property, specifying the name
of the developmental stage that this operator acts upon.
"""
abstract Crossover <: Variation

operate!(c::Crossover, srcs::IndividualCollection, n::Int) =
  operate!(c, srcs.stages[c.stage], n)

function operate!{T}(c::Crossover, srcs::Vector{IndividualStage{T}}, n::Int)
  n_in = num_inputs(c)
  n_out = num_outputs(c)
  n_calls = ceil(n / n_out)
  n_in_from = n_in_to = n_out_from = n_out_to = 0
  for i = 1:n_calls
    # Calculate range of input chromosomes.
    n_in_from = n_in_to + 1
    n_in_to = n_in_to + n_in

    # Calculate range of output chromosomes.
    n_out_from = n_out_to + 1
    n_out_to = n_out_to + n_out

    operate!(c, srcs[n_out_from:n_out_to], srcs[n_in_from:n_in_to])
  end
end

function operate!{T}(
  c::Crossover,
  outputs::Vector{IndividualStage{T}},
  inputs::Vector{IndividualStage{T}}
)
  # Ensure that all the provided inputs are in a valid state.
  # If not, leave the output buffer unchanged.
  try; operate!(c, outputs, map(get, inputs)...); end
end

"""
By default, crossover operators accept two chromosomes as input.
"""
num_inputs(::Crossover) = 2

"""
By default, crossover operators produce two chromosomes as output.
"""
num_outputs(::Crossover) = 2

"""
Calculates the number of (internal) calls required to this operator to
produce a desired number of offspring.
"""
num_calls_required(c::Crossover, n::Int) = convert(Int, ceil(n / num_outputs(c)))

"""
Calculates the number of input chromosomes required to this operator in order
to produce a specified number of offspring.
"""
num_required(c::Crossover, n::Int) = num_calls_required(c, n) * num_inputs(c)

# Include each of the crossover operators.
include("crossover/null.jl")
include("crossover/one_point.jl")
end


text:
#Julia 1.4 Environment
using Pkg
pkg"add ForwardDiff; precompile;"
pkg"add Plots; precompile;"
pkg"add PyPlot; precompile;"
import ForwardDiff  # calcula derivadas usando automatic differentiation in forward mode
using LinearAlgebra #adiciona operações de algebra linear

pyplot()#usa PyPlot como backend engine
∇ = (h,x)->ForwardDiff.gradient(h,x) #gradiente de h(x)
Hessian = (h,x)->ForwardDiff.hessian(h,x) # hessian de h(x)
Df = (h,x)->ForwardDiff.derivative(h,x);  # primeira derivada de h'(x)
D_2f = (h,x)->ForwardDiff.derivative(z->ForwardDiff.derivative(h,z),x) # segunda derivada h"(x)

#função de busca linear
function linear_search(f::Function,x::Vector,d::Vector,method::Function)
    ϕ = t->f(x+t*d)
    α = method(f,x,d)
    return (ϕ(α),x+α*d)
end

function secao_aurea(f::Function,x::Vector,d::Vector;φ::Function = t->f(x+t*d),ρ = 1/2,ϵ = 1e-9)::Float64
    θ_1 = (3 - sqrt(5))/2.0
    θ_2 = (sqrt(5) - 1)/2.0
    #Obtenção intervalo
    a=0
    b=2*ρ
    s = b / 2
    while(φ(b)<φ(s))
        a = s
        s = b
        b *= 2
    end
    #Obtenção de t
    u = a + θ_1*(b - a)
    v = a + θ_2*(b - a)
    while(b - a > ϵ)
        if(φ(u)<φ(v))
            b = v
            v = u
            u = a + θ_1*(b - a)
        else
            a = u
            u = v
            v = a + θ_2*(b - a)
        end
    end
    return (u+v)/2
end

function secao_aurea(f::Function,x::Vector,d::Vector,(a,b)::Tuple{Float64,Float64};φ::Function = t->f(x+t*d),ϵ = 1e-9)::Float64
    θ_1 = (3 - sqrt(5))/2.0
    θ_2 = (sqrt(5) - 1)/2.0
    #Obtenção de t
    u = a + θ_1*(b - a)
    v = a + θ_2*(b - a)
    while(b - a > ϵ)
        if(φ(u)< φ(v))
            b = v
            v = u
            u = a + θ_1*(b - a)
        else
            a = u
            u = v
            v = a + θ_2*(b - a)
        end
    end
    return (u+v)/2
end

function newton(f::Function,x::Vector,d::Vector;φ::Function = y->f(x+y*d),t::Float64=0.25,ϵ = 1e-9,n_iter=10000)::Float64
   iter = 0
   while (abs(Df(φ,t)) > ϵ && iter < n_iter)
         t = t - (Df(φ,t)/D_2f(φ,t))
         iter = iter + 1
   end
   return t
end

function armijo(f::Function,x::Vector,d::Vector;φ::Function = y->f(x+y*d),t::Float64= 1.0,η::Float64 = 0.25)
    # f(x + td) > f(x) +t(∇f(x)⋅ d)
    while (φ(t)>φ(0)+ η*t*(∇(f,x)⋅d))
      t *=0.8
    end
    return t
end

max_stars_repo_name: Huebr/BuscaLinear

max_issues_repo_name: Huebr/BuscaLinear

max_forks_repo_name: Huebr/BuscaLinear


text:
export ElectronSite

const ElectronSite = TagType"Electron"

function siteinds(::ElectronSite,
                  N::Int; kwargs...)
  return [Index(4,"Site,Electron,n=$n") for n=1:N]
end

function state(::ElectronSite,
               st::AbstractString)
  if st == "0" || st == "Emp"
    return 1
  elseif st == "Up" || st == "↑"
    return 2
  elseif st == "Dn" || st == "↓"
    return 3
  elseif st == "UpDn" || st == "↑↓"
    return 4
  end
  throw(ArgumentError("State string \"$st\" not recognized for Electron site"))
  return 0
end

function op(::ElectronSite,
            s::Index,
            opname::AbstractString)::ITensor
  sP = prime(s)
  Emp   = s(1)
  EmpP  = sP(1)
  Up    = s(2)
  UpP   = sP(2)
  Dn    = s(3)
  DnP   = sP(3)
  UpDn  = s(4)
  UpDnP = sP(4)

  Op = ITensor(dag(s), s')

  if opname == "Nup"
    Op[Up, UpP] = 1.
    Op[UpDn, UpDnP] = 1.
  elseif opname == "Ndn"
    Op[Dn, DnP] = 1.
    Op[UpDn, UpDnP] = 1.
  elseif opname == "Ntot"
    Op[Up, UpP] = 1.
    Op[Dn, DnP] = 1.
    Op[UpDn, UpDnP] = 2.
  elseif opname == "Cup" || opname == "Aup"
    Op[Up, EmpP] = 1.
    Op[UpDn, DnP] = 1.
  elseif opname == "Cdagup" || opname == "Adagup"
    Op[Emp, UpP] = 1.
    Op[Dn, UpDnP] = 1.
  elseif opname == "Cdn" || opname == "Adn"
    Op[Dn, EmpP] = 1.
    Op[UpDn, UpP] = 1.
  elseif opname == "Cdagdn" || opname == "Adagdn"
    Op[Emp, DnP] = 1.
    Op[Up, UpDnP] = 1.
  elseif opname=="F" || opname=="FermiPhase" || opname=="FP"
    Op[Up, UpP] = -1.
    Op[Emp, EmpP] = 1.
    Op[Dn, DnP] = -1.
    Op[UpDn, UpDnP] = 1.
  elseif opname == "Fup"
    Op[Emp, EmpP] = 1.
    Op[Up, UpP] = -1.
    Op[Dn, DnP] = 1.
    Op[UpDn, UpDnP] = -1.
  elseif opname == "Fdn"
    Op[Emp, EmpP] = 1.
    Op[Up, UpP] = 1.
    Op[Dn, DnP] = -1.
    Op[UpDn, UpDnP] = -1.
  elseif opname == "Sᶻ" || opname == "Sz"
    Op[Up, UpP] = 0.5
    Op[Dn, DnP] = -0.5
  elseif opname == "Sˣ" || opname == "Sx"
    Op[Up, DnP] = 0.5
    Op[Dn, UpP] = 0.5
  elseif opname=="S+" || opname=="Sp" || opname == "S⁺" || opname == "Splus"
    Op[Dn, UpP] = 1.0
  elseif opname=="S-" || opname=="Sm" || opname == "S⁻" || opname == "Sminus"
    Op[Up, DnP] = 1.0
  elseif opname == "Emp" || opname == "0"
    pEmp = ITensor(s)
    pEmp[Emp] = 1.0
    return pEmp
  elseif opname == "Up" || opname == "↑"
    pU = ITensor(s)
    pU[Up] = 1.0
    return pU
  elseif opname == "Dn" || opname == "↓"
    pD = ITensor(s)
    pD[Dn] = 1.0
    return pD
  elseif opname == "UpDn" || opname == "↑↓"
    pUD = ITensor(s)
    pUD[UpDn] = 1.0
    return pUD
  else
    throw(ArgumentError("Operator name $opname not recognized for ElectronSite"))
  end
  return Op
end


text:
"""
In:

* `D`  [m],           Conductor diameter
* `p_f`  [kg/m^3],      Air density
* `V_w`  [m/s],         Wind speed
* `mu_f` [kg/m-s],      Dynamic viscosity of air

Out: `N_Re` [-], Reynolds number
"""
function eq2c_N_Re(D::Float64, p_f::Float64, V_w::Float64, mu_f::Float64)
    (D * p_f * V_w) / mu_f
end

"""
In:

* `K_angle` [deg],      Wind direction factor
* `N_Re` [-],           Reynolds number
* `k_f` [W/m-C],        Thermal conductivity of air
* `T_s` [C],            Conductor surface temperature
* `T_a` [C],            Ambient temperature

Out: `q_c1` [W/m], Forced convection at low wind speeds
"""
function eq3a_q_c1(
    K_angle::Float64,
    N_Re::Float64,
    k_f::Float64,
    T_s::Float64,
    T_a::Float64
    )
    K_angle * (1.01 + 1.35 * N_Re^0.52) * k_f * (T_s - T_a)
end

"""
In:

* `K_angle` [deg],      Wind direction factor
* `N_Re` [-],           Reynolds number
* `k_f` [W/m-C],        Thermal conductivity of air
* `T_s` [C],            Conductor surface temperature
* `T_a` [C],            Ambient temperature

Out: `q_c2` [W/m], Forced convection at high wind speeds
"""
function eq3b_q_c2(
    K_angle::Float64,
    N_Re::Float64,
    k_f::Float64,
    T_s::Float64,
    T_a::Float64
    )
    K_angle * 0.754 * (N_Re^0.6) * k_f * (T_s - T_a)
end

"""
In: `phi` [deg], Angle between wind and conductor axis

Out: K_angle [-], Wind direction factor
"""
function eq4a_K_angle(phi::Float64)
    1.194 - cosd(phi) + 0.194 * cosd(2 * phi) + 0.368 * sind(2 * phi)
end

"""
In:

* `T_s` [C],            Conductor surface temperature
* `T_a` [C],            Ambient temperature

Out: `T_film` [C], Film temperature
"""
function eq6_T_film(T_s::Float64, T_a::Float64)
    (T_s + T_a) / 2
end

"""
In:

* `alpha` [-],          Solar absorptivity
* `Q_se` [W/m^2],       Elevation-corrected heat flux density
* `theta` [deg],        Sun ray incidence angle
* `A_prime` [m^2],      Conductor area

Out: `q_s`/`eta_s` [W/m], Solar heat gain rate/coefficient
"""
function eq8_q_s(
    alpha::Float64,
    Q_se::Float64,
    theta::Float64,
    A_prime::Float64
    )
    alpha * Q_se * sind(theta) * A_prime
end

"""
In:

* `H_c` [deg], Solar altitude
* `Z_c` [deg], Solar azimuth
* `Z_l` [deg], Line azimuth; East-West => pi/2, North-South => 0

Out: `theta` [deg], Sun ray incidence angle
"""
function eq9_theta(
    H_c::Float64,
    Z_c::Float64,
    Z_l::Float64
    )
    acosd(cosd(H_c) * cosd(Z_c - Z_l))
end

"""
In: `T_film` [C], Film temperature

Out: `mu_f` [kg/m-s], Dynamic viscosity of air
"""
function eq13a_mu_f(T_film::Float64)
    (1.458e-6 * (T_film + 273)^1.5) / (T_film + 383.4)
end

"""
In:

* `H_e` [m],            Conductor elevation
* `T_film` [C],         Film temperature

Out: `p_f` [kg/m^3], Air density
"""
function eq14a_p_f(H_e::Float64, T_film::Float64)
    (1.293 - 1.525e-4 * H_e + 6.379e-9 * H_e^2) / (1 + 0.00367 * T_film)
end

"""
In: `T_film` [C], Film temperature

Out: `k_f` [W/m-C], Thermal conductivity of air
"""
function eq15a_k_f(T_film::Float64)
    @evalpoly(T_film, 2.424e-2, 7.477e-5, -4.407e-9)
end

"""
In: `hours_from_noon`

Out: `omega` [deg], Hour angle
"""
eq_omega(hours_from_noon::Float64) = hours_from_noon * 15.0

"""
In:

* `lat` [deg],          Latitude
* `delta` [deg],        Solar declination
* `omega` [deg],        Hour angle

Out: `H_c` [deg], Solar altitude
"""
function eq16a_H_c(lat::Float64, delta::Float64, omega::Float64)
    asind(cosd(lat) * cosd(delta) * cosd(omega) + sind(lat) * sind(delta))
end

"""
In: `N` [-], Day of year

Out: `delta` [deg], Solar declination
"""
function eq16b_delta(N::Number)
    23.46 * sind((284 + N) * 360 / 365)
end

"""
In:

* `C` [deg],            Solar azimuth constant
* `chi` [deg],          Solar azimuth variable

Out: `Z_c` [deg], Solar azimuth
"""
function eq17a_Z_c(C::Float64, chi::Float64)
    C + atand(chi)
end

"""
In:

* `omega` [deg],        Hour angle
* `lat` [deg],          Latitude
* `delta` [deg],        Solar declination

Out: `chi` [deg], Solar azimuth variable
"""
function eq17b_chi(omega::Float64, lat::Float64, delta::Float64)
    sind(omega) / (sind(lat) * cosd(omega) - cosd(lat) * tand(delta))
end

"""
In:

* `omega` [deg],        Hour angle
* `chi` [deg],          Solar azimuth variable

Out: `C` [deg], Solar azimuth constant
"""
function eqtable2_C(omega::Float64, chi::Float64)
    if -180 <= omega < 0
        C = (chi >= 0) ? 0.0 : 180.0
    elseif 0 <= omega < 180
        C = (chi >= 0) ? 180.0 : 360.0
    end
end

"""
In: `H_c` [deg], Solar altitude

Out: `Q_s` [W/m^2], Heat flux density

Note: Assumes clear atmosphere and SI units in Table 3
"""
function eq18_Q_s(H_c::Float64)
    @evalpoly(
        H_c, -42.2391, 63.8044, -1.9220, 3.46921e-2,
        -3.61118e-4, 1.94318e-6, -4.07608e-9
    )
end

"""
In:

* `K_solar` [-],        Solar altitude correction factor
* `Q_s` [W/m^2],        Heat flux density

Out: `Q_se` [W/m^2], Elevation-corrected heat flux density
"""
eq19_Q_se(K_solar::Float64, Q_s::Float64) = K_solar * Q_s

"""
In: `H_e` [m], Height above sea level

Out: `K_solar` [-], Solar altitude correction factor
"""
eq20_K_solar(H_e::Float64) = @evalpoly(H_e, 1.0, 1.148e-4, -1.108e-8)

"""
In:

* `k_f` [W/m-C],        Thermal conductivity of air
* `K_angle` [deg],      Wind direction factor
* `N_Re` [-],           Reynolds number

Out: `eta_c` [W/m-C], Conductive heat rate coefficient
"""
function eq_eta_c(k_f::Float64, K_angle::Float64, N_Re::Float64)
    K_angle * k_f * max(1.01 + 1.35 * N_Re^0.52, 0.754 * N_Re^0.6)
end

"""
In:

* `D` [m],      Conductor diameter
* `emm` [-],    Emissivity

Out: `eta_r` [W/m-C^4], Radiative heat loss rate coefficient
"""
eq_eta_r(D::Float64, emm::Float64) = 17.8 * D * emm / 1e8

"""
In:

* `Al_m` [kg/m],        Aluminum mass
* `St_m` [kg/m],        Steel mass

Out: `mCp` [J/m-C], Conductor heat capacity
"""
function eq_mCp(Al_m::Float64, St_m::Float64)
    Al_Cp = 955.0 # Aluminum heat capacity
    St_Cp = 476.0 # Steel heat capacity
    Al_m * Al_Cp + St_m * St_Cp
end

eq_mCp(acsr_spec::ACSRSpecsMetric) = eq_mCp(acsr_spec.Al_m, acsr_spec.St_m)


text:
using MixedModelsPlots
using Test

@testset "MixedModelsPlots.jl" begin
    # Write your tests here.
end


text:
# MIT License
#
# Copyright (c) 2020 Joshua E Gentry

# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

using Test
using JSON
using Base.Iterators

const DATA    = JSON.parsefile("$(@__DIR__)/test_vectors.json")
const KEY     = Vector{UInt8}(DATA["key"])
const CONTEXT = DATA["context_string"]
const CASES   = DATA["cases"]

struct Case
    input::Vector{UInt8}
    hash::Vector{UInt8}
    keyed::Vector{UInt8}
    derive::Vector{UInt8}
    Case(case::Dict{String, Any}) = new(
        [x for x = take(cycle(UInt8(0):UInt8(250)), case["input_len"]) ],
        hex2bytes(case["hash"]),
        hex2bytes(case["keyed_hash"]),
        hex2bytes(case["derive_key"])
    )
end

include("stress.jl")

@testset "SIMD128 Implmentation" begin
    include("test_simd_128.jl")
end

@testset "StaticArrays Implmentation" begin
    include("test_w_static_arrays.jl")
end

@testset "Reference Implmentation" begin
    include("test_ref.jl")
end

@testset "Default Usage" begin
    include("test_default.jl")
end


text:
# To add a new cell, type '# %%'
# To add a new markdown cell, type '# %% [markdown]'
# %%
using HDF5
using Plots
using GR
# using QuadGK
using Printf
using LaTeXStrings
using Statistics
ENV["GKSwstype"] = "100"
# %%

rootdirectory = "/home/oku/SimulationData/isogal"
tstart = 100
time = 100 # number of snapshots
H = 0.25
rlim = 10 # limit of radius (active only plane mode)
binwidthr = 0.25
binwidthrho = 0.2
rhomin = -30
rhomax = -18
binwidthsigma = 0.2
sigmamin = -3
sigmamax = 4
nSF = 10
TSF = 1e4
UnitLength_in_cm = 3.085678e21
UnitMass_in_g = 1.989e43
UnitDensity_in_cgs = UnitMass_in_g / UnitLength_in_cm^3
HydrogenMassFrac = 0.76
ProtonMass = 1.6726e-24

struct Param
  H::Float32
  rlim::Float32
  binwidthr::Float32
  binwidthrho::Float32
  binwidthsigma::Float32
  rhomin::Float32
  rhomax::Float32
  sigmamin::Float32
  sigmamax::Float32
  nSF::Float32
  TSF::Float32
  UnitDensity_in_cgs::Float32
  HydrogenMassFrac::Float32
  ProtonMass::Float32
end

param = Param(H, rlim, binwidthr, binwidthrho, binwidthsigma, rhomin, rhomax, sigmamin, sigmamax, nSF, TSF, UnitDensity_in_cgs, HydrogenMassFrac, ProtonMass)

models = [
    ["ss_model/ver10081/fiducial", "Fiducial", "green", "solid"],
    ["ss_model/ver10081/nonstochastic", "Non-stochastic", "blue", "dash"],
    ["ss_model/ver10081/sniikinetic", "SNII kinetic", "purple", "dot"],
    ["ss_model/ver10081/sniithermal", "SNII thermal", "red", "dashdot"],
    ["ss_model/ver10081/nofb", "No FB", "black", "dashdotdot"],
    # ["ss_model/ver10081/med", "High-reso", "red", "dash"],
    # ["ss_model/ver08261/med_T5", "08261 T5", "blue", "solid"],
    # ["ss_model/ver08261/med_r2", "08261 med_r2", "blue", "solid"],
    # ["ss_model/ver08261/med_p2", "08261 med_p2", "purple", "solid"],
    # ["original/Halov3", "Osaka", "black", "dashdotdot"],
        ]
# %%

function getcenter(sph)
    x = y = z = mtot = 0
    for i in 1:length(sph["mass"])
      x += sph["rho"][i] * sph["mass"][i] * sph["pos"][1, i]
      y += sph["rho"][i] * sph["mass"][i] * sph["pos"][2, i]
      z += sph["rho"][i] * sph["mass"][i] * sph["pos"][3, i]
      mtot += sph["rho"][i] * sph["mass"][i]
    end
    x /= mtot
    y /= mtot
    z /= mtot
    return [x, y, z]
end
function getvelocity(sph)
    x = y = z = mtot = 0
    for i in 1:length(sph["mass"])
      x += sph["mass"][i] * sph["vel"][1, i]
      y += sph["mass"][i] * sph["vel"][2, i]
      z += sph["mass"][i] * sph["vel"][3, i]
      mtot += sph["mass"][i]
    end
    x /= mtot
    y /= mtot
    z /= mtot
    return [x, y, z]
end
function padding(array)
    ret = []
    for i in 1:length(array)
        value = array[i] > 0 ? array[i] : 1e-10
        push!(ret, value)
    end
    return ret
end

function calcsigmaz(sph, galaxy)
    numbinr = param.rlim/param.binwidthr
    vz_r = zeros(Int(numbinr))
    m_r = zeros(Int(numbinr))
    for i in 1:length(sph["mass"])
        x = sph["pos"][1, i] - galaxy["pos"][1, 1]
        y = sph["pos"][2, i] - galaxy["pos"][2, 1]
        z = sph["pos"][3, i] - galaxy["pos"][3, 1]
        vz = sph["vel"][3, i] - galaxy["vel"][3, 1]
        m = sph["mass"][i]
        rho = sph["rho"][i]

        # nH = rho * param.UnitDensity_in_cgs * param.HydrogenMassFrac / param.ProtonMass
        # if nH > param.nSF && sph["temp"][i] < param.TSF

            r = sqrt(x^2 + y^2)
            if r < param.rlim && abs(z) < param.H
                ibinr = Int(ceil(r/param.binwidthr))
                vz_r[ibinr] += m * (vz)^2
                # vz_r[ibinr] += m * z^2
                m_r[ibinr] += m
            end
        # end
    end
    vz_r = vz_r./m_r
    vz_r = sqrt.(vz_r)

    return vz_r
end

function calcsigma(sph, galaxy)
    numbinr = param.rlim/param.binwidthr
    vr_r = zeros(Int(numbinr))
    m_r = zeros(Int(numbinr))
    for i in 1:length(sph["mass"])
        x = sph["pos"][1, i] - galaxy["pos"][1, 1]
        y = sph["pos"][2, i] - galaxy["pos"][2, 1]
        z = sph["pos"][3, i] - galaxy["pos"][3, 1]
        vx = sph["vel"][1, i] - galaxy["vel"][1, 1]
        vy = sph["vel"][2, i] - galaxy["vel"][2, 1]
        vz = sph["vel"][3, i] - galaxy["vel"][3, 1]
        m = sph["mass"][i]

        r = sqrt(x^2 + y^2)
        vr = (-y*vx + x*vy)/r

        if r < param.rlim && abs(z) < param.H
            ibinr = Int(ceil(r/param.binwidthr))
            vr_r[ibinr] += m*vr
            m_r[ibinr] += m
        end
    end
    vr_r = vr_r./m_r

    sph["sigma"] = zeros(length(sph["mass"]))
    for i in 1:length(sph["mass"])
        x = sph["pos"][1, i] - galaxy["pos"][1, 1]
        y = sph["pos"][2, i] - galaxy["pos"][2, 1]
        z = sph["pos"][3, i] - galaxy["pos"][3, 1]
        vx = sph["vel"][1, i] - galaxy["vel"][1, 1]
        vy = sph["vel"][2, i] - galaxy["vel"][2, 1]
        vz = sph["vel"][3, i] - galaxy["vel"][3, 1]
        m = sph["mass"][i]

        r = sqrt(x^2 + y^2)

        if r < param.rlim && abs(z) < param.H
            ibinr = Int(ceil(r/param.binwidthr))
            vr = vr_r[ibinr]
            vrx = vr * -y / r
            vry = vr * x / r

            sigma = (vx - vrx)^2 + (vy - vry)^2 + vz^2
            sph["sigma"][i] = sqrt(sigma)
        end
    end
    return vr_r
end

function main(sph, galaxy)
    ret = Dict()
    numbinrho = (param.rhomax - param.rhomin)/param.binwidthrho
    numbinsigma = (param.sigmamax - param.sigmamin)/param.binwidthsigma

    ret["dm/dlogrho"] = zeros(Int(numbinrho))
    ret["dm/dlogsigma"] = zeros(Int(numbinsigma))
    ret["dm/dlogsigmaSF"] = zeros(Int(numbinsigma))
    ret["ddm/dlogcsdlogsigma"] = zeros(Int(numbinsigma), Int(numbinsigma))
    ret["ddm/dlogcsdlogsigmaSF"] = zeros(Int(numbinsigma), Int(numbinsigma))
    Mtot = 0 # total gas mass
    MSF = 0 # total starforming gas mass
    for i in 1:length(sph["mass"])
        x = sph["pos"][1, i] - galaxy["pos"][1, 1]
        y = sph["pos"][2, i] - galaxy["pos"][2, 1]
        z = sph["pos"][3, i] - galaxy["pos"][3, 1]
        vx = sph["vel"][1, i] - galaxy["vel"][1, 1]
        vy = sph["vel"][2, i] - galaxy["vel"][2, 1]
        vz = sph["vel"][3, i] - galaxy["vel"][3, 1]
        m = sph["mass"][i]
        rho = sph["rho"][i]
        sigma = sph["sigma"][i]
        temp = sph["temp"][i]
        cs = sqrt(5/3 * sph["press"][i] / sph["dens"][i])
        Mach = sigma / cs

        H = abs(z)
        R = sqrt(x^2 + y^2)
        if H < param.H && R < param.rlim
            Mtot += m

            ibinrho = Int(ceil((log10(rho*param.UnitDensity_in_cgs) - param.rhomin)/ param.binwidthrho))
            ret["dm/dlogrho"][ibinrho] += m

            ibinMach = Int(ceil((log10(Mach) - param.sigmamin)/param.binwidthsigma))
            ret["dm/dlogsigma"][ibinMach] += m

            ibinsigma = Int(ceil((log10(sigma) - param.sigmamin)/param.binwidthsigma))
            ibincs = Int(ceil((log10(cs) - param.sigmamin)/param.binwidthsigma))

            ret["ddm/dlogcsdlogsigma"][ibinsigma, ibincs] += m

            nH = rho * param.UnitDensity_in_cgs * param.HydrogenMassFrac / param.ProtonMass
            if nH > param.nSF && temp < param.TSF
                MSF += m

                ret["dm/dlogsigmaSF"][ibinMach] += m

                ret["ddm/dlogcsdlogsigmaSF"][ibinsigma, ibincs] += m

            end
        end
    end

    # ret["dm/dlogrho"] = ret["dm/dlogrho"]./Mtot
    # ret["dm/dlogsigma"] = ret["dm/dlogsigma"]./Mtot
    # ret["dm/dlogsigmaSF"] = ret["dm/dlogsigmaSF"]./MSF
    ret["dm/dlogrho"] = ret["dm/dlogrho"]./param.binwidthrho
    ret["dm/dlogsigma"] = ret["dm/dlogsigma"]./param.binwidthsigma
    ret["dm/dlogsigmaSF"] = ret["dm/dlogsigmaSF"]./param.binwidthsigma
    ret["ddm/dlogcsdlogsigma"] = ret["ddm/dlogcsdlogsigma"]./param.binwidthsigma^2
    ret["ddm/dlogcsdlogsigmaSF"] = ret["ddm/dlogcsdlogsigmaSF"]./param.binwidthsigma^2
    return ret
end

## main part

# %%
using Plots
gr()

Plots.scalefontsizes(1.2)

for j in tstart:time
    Results = Dict()

    Results["dm/dlogrho"] = [[] for i in 1:length(models)]
    Results["dm/dlogsigma"] = [[] for i in 1:length(models)]
    Results["dm/dlogsigmaSF"] = [[] for i in 1:length(models)]
    Results["ddm/dlogcsdlogsigma"] = [zeros(1, 1) for i in 1:length(models)]
    Results["ddm/dlogcsdlogsigmaSF"] = [zeros(1, 1) for i in 1:length(models)]
    Results["rotationalvelocity"] = [[] for i in 1:length(models)]
    Results["sigmaz"] = [[] for i in 1:length(models)]

    @time for i in 1:length(models)
        ## load data
        path2snapshot = @sprintf("%s/%s/snapshot_%03d/snapshot_%03d.hdf5", rootdirectory, models[i][1], j, j)

        sph = Dict()
        galaxy = Dict()
        h5open(path2snapshot, "r") do fp
            sph["pos"] = read(fp,"PartType0/Coordinates")
            sph["vel"] = read(fp,"PartType0/Velocities")
            sph["hsml"] = read(fp, "PartType0/SmoothingLength")
            sph["mass"] = read(fp, "PartType0/Masses").*1e10
            sph["rho"] = read(fp, "PartType0/Density")
            sph["temp"] = read(fp, "PartType0/Temperature")
            sph["u"] = read(fp, "PartType0/InternalEnergy")
            sph["press"] = read(fp, "PartType0/Pressure")
            sph["dens"] = read(fp, "PartType0/EnergyWeightedDensity")
        end
        galaxy["pos"] = getcenter(sph)
        galaxy["vel"] = getvelocity(sph)

        Results["rotationalvelocity"][i] = calcsigma(sph, galaxy)

        Results["sigmaz"][i] = calcsigmaz(sph, galaxy)

        ret = main(sph, galaxy)

        ## save
        for l in collect(keys(ret))
            Results[l][i]= ret[l]
        end
    end

    RHO = 0:(param.rhomax - param.rhomin)/param.binwidthrho - 1
    RHO = RHO .* param.binwidthrho
    RHO = RHO .+ (param.rhomin + 0.5*param.binwidthrho)
    Plots.plot(ylim=(1e5, 1e10), xlim=(-26, -20), yscale=:log10, xlabel=latexstring("\$\\log\\,(\\rho\\,[\\mathrm{g cm}^{-3}])\$"), ylabel=latexstring("\$dm/d\\log \\rho\\,[M_\\odot]\$"))
    for i in 1:length(models)
        Plots.plot!(RHO, Results["dm/dlogrho"][i], label=models[i][2], color=Symbol(models[i][3]), linestyle=Symbol(models[i][4]))
    end
    figurename = @sprintf("movie/dmdlogrho%03d.png", j)
    Plots.savefig(figurename)

    SIGMA = 0:(param.sigmamax - param.sigmamin)/param.binwidthsigma - 1
    SIGMA = SIGMA .* param.binwidthsigma
    SIGMA = SIGMA .+ (param.sigmamin + 0.5*param.binwidthsigma)
    Plots.plot(ylim=(1e5, 1e10), yscale=:log10, xlabel=latexstring("\$\\log\\,(\\sigma / c_s)\$"), ylabel=latexstring("\$dm/d\\log \\sigma\\,[M_\\odot]\$"))
    # Plots.plot(ylim=(1e-10, 1e0), yscale=:log10, xlabel=latexstring("\$\\log\\,(\\sigma\\,[\\mathrm{km s}^{-1}])\$"), ylabel=latexstring("\$dm/d\\log \\sigma\$"))
    for i in 1:length(models)
        Plots.plot!(SIGMA, Results["dm/dlogsigma"][i], label=models[i][2], color=Symbol(models[i][3]), linestyle=Symbol(models[i][4]))
    end
    figurename = @sprintf("movie/dmdlogsigma%03d.png", j)
    Plots.savefig(figurename)

    Plots.plot(ylim=(1e5, 1e10), yscale=:log10, xlabel=latexstring("\$\\log\\,(\\sigma / c_s)\$"), ylabel=latexstring("\$dm/d\\log \\sigma\\,[M_\\odot]\$"))
    # Plots.plot(ylim=(1e-10, 1e0), yscale=:log10, xlabel=latexstring("\$\\log\\,(\\sigma\\,[\\mathrm{km s}^{-1}])\$"), ylabel=latexstring("\$dm/d\\log \\sigma\$"))
    for i in 1:length(models)
        Plots.plot!(SIGMA, Results["dm/dlogsigmaSF"][i], label=models[i][2], color=Symbol(models[i][3]), linestyle=Symbol(models[i][4]))
    end
    figurename = @sprintf("movie/dmdlogsigmaSF%03d.png", j)
    Plots.savefig(figurename)

    Plots.contour(xlim=(0, 3), ylim=(0, 3), xlabel=latexstring("\$\\log\\,c_s\$"), ylabel=latexstring("\$\\log\\,v\$"), levels=[6, 7, 8, 9])
    # Plots.plot(ylim=(1e-10, 1e0), yscale=:log10, xlabel=latexstring("\$\\log\\,(\\sigma\\,[\\mathrm{km s}^{-1}])\$"), ylabel=latexstring("\$dm/d\\log \\sigma\$"))
    for i in 1:length(models)
        Plots.contour!(SIGMA, SIGMA, log10.(Results["ddm/dlogcsdlogsigmaSF"][i]), label=models[i][2], color=Symbol(models[i][3]), linestyle=Symbol(models[i][4]), levels=[6, 7, 8, 9])
    end
    figurename = @sprintf("movie/ddmdlogcsdlogsigmaSF%03d.png", j)
    Plots.savefig(figurename)

    R = 0:param.rlim/param.binwidthr - 1
    R = R .* param.binwidthr
    R = R .+ 0.5*param.binwidthr
    Plots.plot(ylim=(0,300), xlim=(0, param.rlim), xlabel=latexstring("\$R\\,[\\mathrm{kpc}]\$"), ylabel=latexstring("\$ V_r\\,[\\mathrm{km s}^{-1}]\$"))
    for i in 1:length(models)
        Plots.plot!(R, Results["rotationalvelocity"][i], label=models[i][2], color=Symbol(models[i][3]), linestyle=Symbol(models[i][4]))
    end
    figurename = @sprintf("movie/rotationalvelocity%03d.png", j)
    Plots.savefig(figurename)

    Plots.plot(xlim=(0, param.rlim), xlabel=latexstring("\$R\\,[\\mathrm{kpc}]\$"), ylabel=latexstring("\$\\sigma_{z}\\,[\\mathrm{km s^{-1}}]\$"))
    for i in 1:length(models)
        Plots.plot!(R, Results["sigmaz"][i], label=models[i][2], color=Symbol(models[i][3]), linestyle=Symbol(models[i][4]))
    end
    figurename = @sprintf("movie/sigmaz%03d.png", j)
    Plots.savefig(figurename)
end


text:
export BrickletOLED128x64V2ReadPixelsLowLevel
struct BrickletOLED128x64V2ReadPixelsLowLevel
    pixels_length::Integer
    pixels_chunk_offset::Integer
    pixels_chunk_data::Vector{Bool}
end

export BrickletOLED128x64V2DisplayConfiguration
struct BrickletOLED128x64V2DisplayConfiguration
    contrast::Integer
    invert::Bool
    automatic_draw::Bool
end

export BrickletOLED128x64V2SPITFPErrorCount
struct BrickletOLED128x64V2SPITFPErrorCount
    error_count_ack_checksum::Integer
    error_count_message_checksum::Integer
    error_count_frame::Integer
    error_count_overflow::Integer
end

export BrickletOLED128x64V2Identity
struct BrickletOLED128x64V2Identity
    uid::String
    connected_uid::String
    position::Char
    hardware_version::Vector{Integer}
    firmware_version::Vector{Integer}
    device_identifier::Integer
end

export BrickletOLED128x64V2
"""
3.3cm (1.3") OLED display with 128x64 pixels
"""
mutable struct BrickletOLED128x64V2 <: TinkerforgeDevice
    ipcon::IPConnection
    deviceInternal::PyObject

    """
    Creates an object with the unique device ID *uid* and adds it to
    the IP Connection *ipcon*.
    """
    function BrickletOLED128x64V2(uid::String, ipcon::IPConnection)
        package = pyimport("tinkerforge.bricklet_oled_128x64_v2")
        deviceInternal = package.BrickletOLED128x64V2(uid, ipcon.ipconInternal)

        return new(ipcon, deviceInternal)
    end
end

export write_pixels_low_level
"""
    $(SIGNATURES)

Writes pixels to the specified window.

The pixels are written into the window line by line top to bottom
and each line is written from left to right.

If automatic draw is enabled (default) the pixels are directly written to
the screen. Only pixels that have actually changed are updated on the screen,
the rest stays the same.

If automatic draw is disabled the pixels are written to an internal buffer and
the buffer is transferred to the display only after :func:`Draw Buffered Frame`
is called. This can be used to avoid flicker when drawing a complex frame in
multiple steps.

Automatic draw can be configured with the :func:`Set Display Configuration`
function.
"""
function write_pixels_low_level(device::BrickletOLED128x64V2, x_start, y_start, x_end, y_end, pixels_length, pixels_chunk_offset, pixels_chunk_data)
    device.deviceInternal.write_pixels_low_level(x_start, y_start, x_end, y_end, pixels_length, pixels_chunk_offset, pixels_chunk_data)
end

export read_pixels_low_level
"""
    $(SIGNATURES)

Reads pixels from the specified window.

The pixels are read from the window line by line top to bottom
and each line is read from left to right.

If automatic draw is enabled (default) the pixels that are read are always the
same that are shown on the display.

If automatic draw is disabled the pixels are read from the internal buffer
(see :func:`Draw Buffered Frame`).

Automatic draw can be configured with the :func:`Set Display Configuration`
function.
"""
function read_pixels_low_level(device::BrickletOLED128x64V2, x_start, y_start, x_end, y_end)
    return device.deviceInternal.read_pixels_low_level(x_start, y_start, x_end, y_end)
end

export clear_display
"""
    $(SIGNATURES)

Clears the complete content of the display.

If automatic draw is enabled (default) the pixels are directly cleared.

If automatic draw is disabled the the internal buffer is cleared and
the buffer is transferred to the display only after :func:`Draw Buffered Frame`
is called. This can be used to avoid flicker when drawing a complex frame in
multiple steps.

Automatic draw can be configured with the :func:`Set Display Configuration`
function.
"""
function clear_display(device::BrickletOLED128x64V2)
    device.deviceInternal.clear_display()
end

export set_display_configuration
"""
    $(SIGNATURES)

Sets the configuration of the display.

You can set a contrast value from 0 to 255 and you can invert the color
(white/black) of the display.

If automatic draw is set to *true*, the display is automatically updated with every
call of :func:`Write Pixels` or :func:`Write Line`. If it is set to false, the
changes are written into an internal buffer and only shown on the display after
a call of :func:`Draw Buffered Frame`.
"""
function set_display_configuration(device::BrickletOLED128x64V2, contrast, invert, automatic_draw)
    device.deviceInternal.set_display_configuration(contrast, invert, automatic_draw)
end

export get_display_configuration
"""
    $(SIGNATURES)

Returns the configuration as set by :func:`Set Display Configuration`.
"""
function get_display_configuration(device::BrickletOLED128x64V2)
    return device.deviceInternal.get_display_configuration()
end

export write_line
"""
    $(SIGNATURES)

Writes text to a specific line with a specific position.
The text can have a maximum of 22 characters.

For example: (1, 10, "Hello") will write *Hello* in the middle of the
second line of the display.

The display uses a special 5x7 pixel charset. You can view the characters
of the charset in Brick Viewer.

If automatic draw is enabled (default) the text is directly written to
the screen. Only pixels that have actually changed are updated on the screen,
the rest stays the same.

If automatic draw is disabled the text is written to an internal buffer and
the buffer is transferred to the display only after :func:`Draw Buffered Frame`
is called. This can be used to avoid flicker when drawing a complex frame in
multiple steps.

Automatic draw can be configured with the :func:`Set Display Configuration`
function.
"""
function write_line(device::BrickletOLED128x64V2, line, position, text)
    device.deviceInternal.write_line(line, position, text)
end

export draw_buffered_frame
"""
    $(SIGNATURES)

Draws the currently buffered frame. Normally each call of :func:`Write Pixels` and
:func:`Write Line` draws directly onto the display. If you turn automatic draw off
(:func:`Set Display Configuration`), the data is written in an internal buffer and
only transferred to the display by calling this function. This can be used to
avoid flicker when drawing a complex frame in multiple steps.

Set the `force complete redraw` to *true* to redraw the whole display
instead of only the changed parts. Normally it should not be necessary to set this to
*true*. It may only become necessary in case of stuck pixels because of errors.
"""
function draw_buffered_frame(device::BrickletOLED128x64V2, force_complete_redraw)
    device.deviceInternal.draw_buffered_frame(force_complete_redraw)
end

export get_spitfp_error_count
"""
    $(SIGNATURES)

Returns the error count for the communication between Brick and Bricklet.

The errors are divided into

* ACK checksum errors,
* message checksum errors,
* framing errors and
* overflow errors.

The errors counts are for errors that occur on the Bricklet side. All
Bricks have a similar function that returns the errors on the Brick side.
"""
function get_spitfp_error_count(device::BrickletOLED128x64V2)
    return device.deviceInternal.get_spitfp_error_count()
end

export set_bootloader_mode
"""
    $(SIGNATURES)

Sets the bootloader mode and returns the status after the requested
mode change was instigated.

You can change from bootloader mode to firmware mode and vice versa. A change
from bootloader mode to firmware mode will only take place if the entry function,
device identifier and CRC are present and correct.

This function is used by Brick Viewer during flashing. It should not be
necessary to call it in a normal user program.
"""
function set_bootloader_mode(device::BrickletOLED128x64V2, mode)
    return device.deviceInternal.set_bootloader_mode(mode)
end

export get_bootloader_mode
"""
    $(SIGNATURES)

Returns the current bootloader mode, see :func:`Set Bootloader Mode`.
"""
function get_bootloader_mode(device::BrickletOLED128x64V2)
    return device.deviceInternal.get_bootloader_mode()
end

export set_write_firmware_pointer
"""
    $(SIGNATURES)

Sets the firmware pointer for :func:`Write Firmware`. The pointer has
to be increased by chunks of size 64. The data is written to flash
every 4 chunks (which equals to one page of size 256).

This function is used by Brick Viewer during flashing. It should not be
necessary to call it in a normal user program.
"""
function set_write_firmware_pointer(device::BrickletOLED128x64V2, pointer)
    device.deviceInternal.set_write_firmware_pointer(pointer)
end

export write_firmware
"""
    $(SIGNATURES)

Writes 64 Bytes of firmware at the position as written by
:func:`Set Write Firmware Pointer` before. The firmware is written
to flash every 4 chunks.

You can only write firmware in bootloader mode.

This function is used by Brick Viewer during flashing. It should not be
necessary to call it in a normal user program.
"""
function write_firmware(device::BrickletOLED128x64V2, data)
    return device.deviceInternal.write_firmware(data)
end

export set_status_led_config
"""
    $(SIGNATURES)

Sets the status LED configuration. By default the LED shows
communication traffic between Brick and Bricklet, it flickers once
for every 10 received data packets.

You can also turn the LED permanently on/off or show a heartbeat.

If the Bricklet is in bootloader mode, the LED is will show heartbeat by default.
"""
function set_status_led_config(device::BrickletOLED128x64V2, config)
    device.deviceInternal.set_status_led_config(config)
end

export get_status_led_config
"""
    $(SIGNATURES)

Returns the configuration as set by :func:`Set Status LED Config`
"""
function get_status_led_config(device::BrickletOLED128x64V2)
    return device.deviceInternal.get_status_led_config()
end

export get_chip_temperature
"""
    $(SIGNATURES)

Returns the temperature as measured inside the microcontroller. The
value returned is not the ambient temperature!

The temperature is only proportional to the real temperature and it has bad
accuracy. Practically it is only useful as an indicator for
temperature changes.
"""
function get_chip_temperature(device::BrickletOLED128x64V2)
    return device.deviceInternal.get_chip_temperature()
end

export reset
"""
    $(SIGNATURES)

Calling this function will reset the Bricklet. All configurations
will be lost.

After a reset you have to create new device objects,
calling functions on the existing ones will result in
undefined behavior!
"""
function reset(device::BrickletOLED128x64V2)
    device.deviceInternal.reset()
end

export write_uid
"""
    $(SIGNATURES)

Writes a new UID into flash. If you want to set a new UID
you have to decode the Base58 encoded UID string into an
integer first.

We recommend that you use Brick Viewer to change the UID.
"""
function write_uid(device::BrickletOLED128x64V2, uid)
    device.deviceInternal.write_uid(uid)
end

export read_uid
"""
    $(SIGNATURES)

Returns the current UID as an integer. Encode as
Base58 to get the usual string version.
"""
function read_uid(device::BrickletOLED128x64V2)
    return device.deviceInternal.read_uid()
end

export get_identity
"""
    $(SIGNATURES)

Returns the UID, the UID where the Bricklet is connected to,
the position, the hardware and firmware version as well as the
device identifier.

The position can be 'a', 'b', 'c', 'd', 'e', 'f', 'g' or 'h' (Bricklet Port).
A Bricklet connected to an :ref:`Isolator Bricklet <isolator_bricklet>` is always at
position 'z'.

The device identifier numbers can be found :ref:`here <device_identifier>`.
|device_identifier_constant|
"""
function get_identity(device::BrickletOLED128x64V2)
    return device.deviceInternal.get_identity()
end

export write_pixels
"""
Writes pixels to the specified window.

The pixels are written into the window line by line top to bottom
and each line is written from left to right.

If automatic draw is enabled (default) the pixels are directly written to
the screen. Only pixels that have actually changed are updated on the screen,
the rest stays the same.

If automatic draw is disabled the pixels are written to an internal buffer and
the buffer is transferred to the display only after :func:`Draw Buffered Frame`
is called. This can be used to avoid flicker when drawing a complex frame in
multiple steps.

Automatic draw can be configured with the :func:`Set Display Configuration`
function.
"""
function write_pixels(device::BrickletOLED128x64V2, x_start, y_start, x_end, y_end, pixels)

        x_start = UInt8(x_start)
    y_start = UInt8(y_start)
    x_end = UInt8(x_end)
    y_end = UInt8(y_end)
    pixels = convert(Vector{Bool}, pixels)

    if length(pixels) > 65535
        throw(TinkerforgeInvalidParameterError("Pixels can be at most 65535 items long"))
    end

    pixels_length = length(pixels)
    pixels_chunk_offset = 0

    if pixels_length == 0
        pixels_chunk_data = [False] * 448
        ret = write_pixels_low_level(device, x_start, y_start, x_end, y_end, pixels_length, pixels_chunk_offset, pixels_chunk_data)
    else
        lock(device.stream_lock) do
            while pixels_chunk_offset < pixels_length
                pixels_chunk_data = create_chunk_data(pixels, pixels_chunk_offset, 448, False)
                ret = write_pixels_low_level(device, x_start, y_start, x_end, y_end, pixels_length, pixels_chunk_offset, pixels_chunk_data)
                pixels_chunk_offset += 448
            end
        end
    end

        return ret
end

export read_pixels
"""
Reads pixels from the specified window.

The pixels are read from the window line by line top to bottom
and each line is read from left to right.

If automatic draw is enabled (default) the pixels that are read are always the
same that are shown on the display.

If automatic draw is disabled the pixels are read from the internal buffer
(see :func:`Draw Buffered Frame`).

Automatic draw can be configured with the :func:`Set Display Configuration`
function.
"""
function read_pixels(device::BrickletOLED128x64V2, x_start, y_start, x_end, y_end)

        x_start = UInt8(x_start)
    y_start = UInt8(y_start)
    x_end = UInt8(x_end)
    y_end = UInt8(y_end)

    lock(device.stream_lock) do
        ret = read_pixels_low_level(device, x_start, y_start, x_end, y_end)
            pixels_length = ret.pixels_length
        pixels_out_of_sync = ret.pixels_chunk_offset != 0

        pixels_data = ret.pixels_chunk_data

        while !pixels_out_of_sync && length(pixels_data) < pixels_length
            ret = read_pixels_low_level(device, x_start, y_start, x_end, y_end)
                pixels_length = ret.pixels_length
            pixels_out_of_sync = ret.pixels_chunk_offset != length(pixels_data)
            pixels_data += ret.pixels_chunk_data
        end

        if pixels_out_of_sync # discard remaining stream to bring it back in-sync
            while ret.pixels_chunk_offset + 480 < pixels_length
                ret = read_pixels_low_level(device, x_start, y_start, x_end, y_end)
                    pixels_length = ret.pixels_length
            end

            throw(TinkerforgeStreamOutOfSyncError("Pixels stream is out-of-sync"))
        end
    end

        return pixels_data[:pixels_length]
end


text:
include("../utils.jl")
include("group_utils.jl")

@testset "Product group" begin
    SOn = SpecialOrthogonal(3)
    Tn = TranslationGroup(2)
    Rn = Rotations(3)
    M = ProductManifold(SOn, Tn)
    G = ProductGroup(M)
    @test_throws ErrorException ProductGroup(ProductManifold(Rotations(3), Stiefel(3,2)))
    @test G isa ProductGroup
    @test submanifold(G, 1) === SOn
    @test submanifold(G, 2) === Tn
    @test base_manifold(G) === M
    @test sprint(show, G) == "ProductGroup($(SOn), $(Tn))"
    @test sprint(show, "text/plain", G) == "ProductGroup with 2 subgroups:\n $(SOn)\n $(Tn)"
    x = Matrix{Float64}(I, 3, 3)
    for f in [group_exp!, group_log!]
        @test Manifolds.decorator_transparent_dispatch(f,G,x,x) === Val{:transparent}()
    end
    t = Vector{Float64}.([1:2, 2:3, 3:4])
    ω = [[1.0, 2.0, 3.0], [3.0, 2.0, 1.0], [1.0, 3.0, 2.0]]
    tuple_pts = [(exp(Rn, x, hat(Rn, x, ωi)), ti) for (ωi, ti) in zip(ω, t)]
    tuple_v = (hat(Rn, x, [1.0, 0.5, -0.5]), [-1.0, 2.0])
    eA = [x, zeros(2)]
    shape_se = Manifolds.ShapeSpecification(Manifolds.ArrayReshaper(), M.manifolds...)
    e = Manifolds.prod_point(shape_se, eA...)

    @testset "identity specializations" begin
        @test inv(G, Identity(G,e)) === Identity(G,e)
        @test identity(G, Identity(G,e)) === Identity(G,e)
        @test submanifold_component(G, Identity(G,e), Val(1)) == Identity(SOn,x)
        @test submanifold_component(G, Identity(G,e), Val(2)) == Identity(Tn, zeros(2))
        @test submanifold_components(G, Identity(G,e)) == (Identity(SOn,x), Identity(Tn,zeros(2)))
        @test compose(G, Identity(G,e), Identity(G,e)) === Identity(G,e)
    end

    @testset "product point" begin
        reshapers = (Manifolds.ArrayReshaper(), Manifolds.StaticReshaper())
        for reshaper in reshapers
            shape_se = Manifolds.ShapeSpecification(reshaper, M.manifolds...)
            pts = [Manifolds.prod_point(shape_se, tp...) for tp in tuple_pts]
            v_pts = [Manifolds.prod_point(shape_se, tuple_v...)]
            @test compose(G, pts[1], Identity(G,e)) == pts[1]
            @test compose(G, Identity(G,e), pts[1]) == pts[1]
            test_group(G, pts, v_pts, v_pts; test_diff = true)
            @test isapprox(
                M,
                identity(M, pts[1]),
                group_exp(M, v_pts[1]),
                Manifolds.prod_point(
                    shape_se,
                    group_exp(SOn, v_pts[1].parts[1]),
                    group_exp(Tn, v_pts[1].parts[2]),
                ),
            )
            @test isapprox(
                M,
                identity(M, pts[1]),
                group_log(M, pts[1]),
                Manifolds.prod_point(
                    shape_se,
                    group_log(SOn, pts[1].parts[1]),
                    group_log(Tn, pts[1].parts[2]),
                ),
            )
        end
    end

    @testset "product repr" begin
        pts = [ProductRepr(tp...) for tp in tuple_pts]
        v_pts = [ProductRepr(tuple_v...)]
        @test compose(G, pts[1], Identity(G,e)) == pts[1]
        @test compose(G, Identity(G,e), pts[1]) == pts[1]
        test_group(G, pts, v_pts, v_pts; test_diff = true, test_mutating = false)
        @test isapprox(
            M,
            group_exp(M, v_pts[1]),
            ProductRepr(
                group_exp(SOn, v_pts[1].parts[1]),
                group_exp(Tn, v_pts[1].parts[2]),
            ),
        )
        @test isapprox(
            M,
            group_log(M, pts[1]),
            ProductRepr(group_log(SOn, pts[1].parts[1]), group_log(Tn, pts[1].parts[2])),
        )
    end
    @test sprint(show, "text/plain", G) === """
    ProductGroup with 2 subgroups:
     SpecialOrthogonal(3)
     TranslationGroup(2; field = ℝ)"""
end


text:
import RecipesBase
function recipetest(lib::Polyhedra.Library)
    v = convexhull([0.0, 0.0], [1.0, 0.0], [0.0, 1.0])
    p = polyhedron(v)
    @test RecipesBase.apply_recipe(Dict{Symbol, Any}(), p)[1].args == ([0.0, 0.0, 1.0, 0.0],
                                                                       [0.0, 1.0, 0.0, 0.0])
    @testset "Error for unbounded polyhedron" begin
        vr = convexhull(v, Ray([1.0, 1.0]))
        pr = polyhedron(vr)
        @test_throws ErrorException RecipesBase.apply_recipe(Dict{Symbol, Any}(), pr)
        vl = convexhull(v, Line([1.0, 1.0]))
        pl = polyhedron(vl)
        @test_throws ErrorException RecipesBase.apply_recipe(Dict{Symbol, Any}(), pl)
    end
end


text:
# This code needs to be reduced by several fold, storing it here for now

###########################  Blotter & Account & Portfolio

# single row
function getindex(b::Union(Blotter,Account,Portfolio), n::Int)
    isa(b,Blotter) ?
    Blotter(b.timestamp[n], b.values[n,:], b.colnames) :
    isa(b,Account) ?
    Account(b.timestamp[n], b.values[n,:], b.colnames) :
    Portfolio(b.timestamp[n], b.values[n,:], b.colnames)
end

# range of rows
function getindex(b::Union(Blotter,Account,Portfolio), r::Range1{Int})
    isa(b,Blotter) ?
    Blotter(b.timestamp[r], b.values[r,:], b.colnames) :
    isa(b,Account) ?
    Account(b.timestamp[r], b.values[r,:], b.colnames) :
    Portfolio(b.timestamp[r], b.values[r,:], b.colnames)
end

# array of rows
function getindex(b::Union(Blotter,Account,Portfolio), a::Array{Int})
    isa(b,Blotter) ?
    Blotter(b.timestamp[a], b.values[a,:], b.colnames) :
    isa(b,Account) ?
    Account(b.timestamp[a], b.values[a,:], b.colnames) :
    Portfolio(b.timestamp[a], b.values[a,:], b.colnames)
end

# single column by name
function getindex(b::Union(Blotter,Account,Portfolio), s::ASCIIString)
    n = findfirst(b.colnames, s)
    isa(b,Blotter) ?
    Blotter(b.timestamp, b.values[:, n], ASCIIString[s]) :
    isa(b,Account) ?
    Account(b.timestamp, b.values[:, n], ASCIIString[s]) :
    Portfolio(b.timestamp, b.values[:, n], ASCIIString[s])
end

# array of columns by name
function getindex(b::Union(Blotter,Account,Portfolio), args::ASCIIString...)
    ns = [findfirst(b.colnames, a) for a in args]
    isa(b,Blotter) ?
    Blotter(b.timestamp, b.values[:,ns], ASCIIString[a for a in args]) :
    isa(b,Account) ?
    Account(b.timestamp, b.values[:,ns], ASCIIString[a for a in args]) :
    Portfolio(b.timestamp, b.values[:,ns], ASCIIString[a for a in args])
end

# single date
function getindex(b::Union(Blotter,Account,Portfolio), d::Date)
   for i in 1:length(b)
     if [d] == b[i].timestamp
       return b[i]
     else
       nothing
     end
   end
 end

# range of dates
function getindex(b::Union(Blotter,Account,Portfolio), dates::Array{Date,1})
  counter = Int[]
  for i in 1:length(dates)
    if findfirst(b.timestamp, dates[i]) != 0
      push!(counter, findfirst(b.timestamp, dates[i]))
    end
  end
  b[counter]
end

function getindex(b::Union(Blotter,Account,Portfolio), r::StepRange)
    b[[r]]
end

# day of week
# getindex{T,N}(b::Union(Blotter,Account,Portfolio){T,N}, d::DAYOFWEEK) = b[dayofweek(b.timestamp) .== d]

# ###########################  Account
#
# # single row
# function getindex(b::Account, n::Int)
#     Account(b.timestamp[n], b.values[n,:], b.colnames)
# end
#
# # range of rows
# function getindex(b::Account, r::Range1{Int})
#     Account(b.timestamp[r], b.values[r,:], b.colnames)
# end
#
# # array of rows
# function getindex(b::Account, a::Array{Int})
#     Account(b.timestamp[a], b.values[a,:], b.colnames)
# end
#
# # single column by name
# function getindex(b::Account, s::ASCIIString)
#     n = findfirst(b.colnames, s)
#     Account(b.timestamp, b.values[:, n], ASCIIString[s])
# end
#
# # array of columns by name
# function getindex(b::Account, args::ASCIIString...)
#     ns = [findfirst(b.colnames, a) for a in args]
#     Account(b.timestamp, b.values[:,ns], ASCIIString[a for a in args])
# end
#
# # single date
# function getindex(b::Account, d::Date)
#    for i in 1:length(b)
#      if [d] == b[i].timestamp
#        return b[i]
#      else
#        nothing
#      end
#    end
#  end
#
# # range of dates
# function getindex(b::Account, dates::Array{Date,1})
#   counter = Int[]
#   for i in 1:length(dates)
#     if findfirst(b.timestamp, dates[i]) != 0
#       push!(counter, findfirst(b.timestamp, dates[i]))
#     end
#   end
#   b[counter]
# end
#
# function getindex(b::Account, r::StepRange)
#     b[[r]]
# end
#
# # day of week
# # getindex{T,N}(b::Account{T,N}, d::DAYOFWEEK) = b[dayofweek(b.timestamp) .== d]
#
# ###########################  Blotter
#
# # single row
# function getindex(b::Blotter, n::Int)
#     Blotter(b.timestamp[n], b.values[n,:], b.colnames, b.ticker)
# end
#
# # range of rows
# function getindex(b::Blotter, r::Range1{Int})
#     Blotter(b.timestamp[r], b.values[r,:], b.colnames, b.ticker)
# end
#
# # array of rows
# function getindex(b::Blotter, a::Array{Int})
#     Blotter(b.timestamp[a], b.values[a,:], b.colnames, b.ticker)
# end
#
# # single column by name
# function getindex(b::Blotter, s::ASCIIString)
#     n = findfirst(b.colnames, s)
#     Blotter(b.timestamp, b.values[:, n], ASCIIString[s], b.ticker)
# end
#
# # single date
# function getindex(b::Blotter, d::Date)
#    for i in 1:length(b)
#      if [d] == b[i].timestamp
#        return b[i]
#      else
#        nothing
#      end
#    end
#  end
#
# # range of dates
# function getindex(b::Blotter, dates::Array{Date,1})
#   counter = Int[]
#   for i in 1:length(dates)
#     if findfirst(b.timestamp, dates[i]) != 0
#       push!(counter, findfirst(b.timestamp, dates[i]))
#     end
#   end
#   b[counter]
# end
#
# ## DOESN'T WORK because of StepRange
# function getindex(b::Blotter, r::StepRange)
#     b[[r]]
# end
#
# # day of week
# # getindex{T,N}(b::Blotter{T,N}, d::DAYOFWEEK) = b[dayofweek(b.timestamp) .== d]
#
# ###########################  Portfolio
#
# # single row
# function getindex(b::Portfolio, n::Int)
#     Portfolio(b.timestamp[n], b.values[n,:], b.colnames)
# end
#
# # range of rows
# function getindex(b::Portfolio, r::Range1{Int})
#     Portfolio(b.timestamp[r], b.values[r,:], b.colnames)
# end
#
# # array of rows
# function getindex(b::Portfolio, a::Array{Int})
#     Portfolio(b.timestamp[a], b.values[a,:], b.colnames)
# end
#
# # single column by name
# function getindex(b::Portfolio, s::ASCIIString)
#     n = findfirst(b.colnames, s)
#     Portfolio(b.timestamp, b.values[:, n], ASCIIString[s])
# end
#
# # array of columns by name
# function getindex(b::Portfolio, args::ASCIIString...)
#     ns = [findfirst(b.colnames, a) for a in args]
#     Portfolio(b.timestamp, b.values[:,ns], ASCIIString[a for a in args])
# end
#
# # single date
# function getindex(b::Portfolio, d::Date)
#    for i in 1:length(b)
#      if [d] == b[i].timestamp
#        return b[i]
#      else
#        nothing
#      end
#    end
#  end
#
# # range of dates
# function getindex(b::Portfolio, dates::Array{Date,1})
#   counter = Int[]
# #  counter = int(zeros(length(dates)))
#   for i in 1:length(dates)
#     if findfirst(b.timestamp, dates[i]) != 0
#       #counter[i] = findfirst(b.timestamp, dates[i])
#       push!(counter, findfirst(b.timestamp, dates[i]))
#     end
#   end
#   b[counter]
# end
#
# # DOESN'T WORK #########
# function getindex(b::Portfolio, r::StepRange)
#     b[[r]]
# end
#
# # day of week
# # getindex{T,N}(b::Portfolio{T,N}, d::DAYOFWEEK) = b[dayofweek(b.timestamp) .== d]


text:
## arithmetic.jl operators for NamedArray

## (c) 2013--2014 David A. van Leeuwen

## This code is licensed under the MIT License
## See the file LICENSE.md in this distribution

import Base: +, -, *, /, .+, .-, .*, ./, \

-(n::NamedArray) = NamedArray(-n.array, n.dicts, n.dimnames)

## disambiguation magic
.*(n::NamedArray{Bool}, b::BitArray) = NamedArray(n.array .* b, n.dicts, n.dimnames)
.*{N}(n::NamedArray{Bool,N}, b::BitArray{N}) = NamedArray(n.array .* b, n.dicts, n.dimnames)
.*{N}(b::BitArray{N}, n::NamedArray{Bool,N}) = n .* b

## disambiguation (Argh...)
for op in (:+, :-)
    @eval ($op){T1<:Number,T2<:Number}(x::Range{T1}, y::NamedVector{T2}) = NamedArray(($op)(x,y.array), y.dicts, y.dimnames)
end

for op in (:+, :-, :.+, :.-, :.*, :./)
    ## named %op% named
    @eval begin
        function ($op){T1<:Number, T2<:Number}(x::NamedArray{T1}, y::NamedArray{T2})
            if allnames(x)==allnames(y) && x.dimnames==y.dimnames
                NamedArray(($op)(x.array, y.array), x.dicts, x.dimnames)
            else
                warn("Dropping mismatching names")
                ($op)(x.array, y.array)
            end
        end
        ($op){T1<:Number,T2<:Number,N}(x::NamedArray{T1,N}, y::AbstractArray{T2,N}) = NamedArray(($op)(x.array, y), x.dicts, x.dimnames)
        ($op){T1<:Number,T2<:Number,N}(x::AbstractArray{T1,N}, y::NamedArray{T2,N}) = NamedArray(($op)(x, y.array), y.dicts, y.dimnames)
    end
end

## scalar arithmetic
## disambiguate
for op in (:+, :-)
    @eval begin
        ($op)(x::NamedArray{Bool}, y::Bool) = NamedArray(($op)(x.array, y), x.dimnames, x.dicts)
        ($op)(x::Bool, y::NamedArray{Bool}) = NamedArray(($op)(x, y.array), y.dimnames, y.dicts)
    end
end

## NamedArray, Number
for op in (:+, :-, :*, :.+, :.-, :.*, :./)
    @eval begin
        ($op){T1<:Number,T2<:Number}(x::NamedArray{T1}, y::T2) = NamedArray(($op)(x.array, y), x.dicts, x.dimnames)
        ($op){T1<:Number,T2<:Number}(x::T1, y::NamedArray{T2}) = NamedArray(($op)(x, y.array), y.dicts, y.dimnames)
    end
end
/{T1<:Number,T2<:Number}(x::NamedArray{T1}, y::T2) = NamedArray(x.array / y, x.dicts, x.dimnames)
\{T1<:Number,T2<:Number}(x::T1, y::NamedArray{T2}) = NamedArray(x \ y.array, y.dicts, y.dimnames)

## disambiguation
-{T1<:Number,T2<:Number}(x::NamedVector{T1}, y::Range{T2}) = NamedArray(x.array - y, x.dicts, x.dimnames)
-{T1,T2}(x::NamedArray{T1}, y::Range{T2}) = NamedArray(x.array - y, x.dicts, x.dimnames)

## NamedArray, AbstractArray, same dimensions
for op in (:+, :-, :.+, :.-, :.*, :./)
    @eval begin
        ($op)(x::NamedArray, y::AbstractArray) = NamedArray(($op)(x,y), x.dicts, x.dimnames)
    end
end

import Base: A_mul_B!, A_mul_Bc!, A_mul_Bc, A_mul_Bt!, A_mul_Bt, Ac_mul_B, Ac_mul_B!, Ac_mul_Bc, Ac_mul_Bc!, At_mul_B, At_mul_B!, At_mul_Bt, At_mul_Bt!

if VERSION >= v"0.5.0-dev" ## v0.4 ambiguity-hell with AbstractTriangular c.s.
    ## Assume dimensions/names are correct
    for op in (:A_mul_B!, :A_mul_Bc!, :A_mul_Bt!, :Ac_mul_B!, :Ac_mul_Bc!, :At_mul_B!, :At_mul_Bt!)
        @eval ($op)(C::NamedMatrix, A::AbstractMatrix, B::AbstractMatrix) = ($op)(C.array, A, array, B.array)
    end
end

for op in (:A_mul_Bc, :A_mul_Bt)
    @eval ($op)(A::NamedMatrix, B::NamedMatrix) = NamedArray(($op)(A.array, B.array), (A.dicts[1], B.dicts[1]), (A.dimnames[1], B.dimnames[1]))
    for T in [Union{Base.LinAlg.QRCompactWYQ, Base.LinAlg.QRPackedQ}, StridedMatrix, AbstractMatrix] ## v0.4 ambiguity-hell
        @eval ($op)(A::NamedMatrix, B::$T) = NamedArray(($op)(A.array, B), (A.dicts[1], defaultnamesdict(size(B,1))), (A.dimnames[1], :B))
        @eval ($op)(A::$T, B::NamedMatrix) = NamedArray(($op)(A, B.array), (defaultnamesdict(size(A,1)), B.dicts[1]), (:A, B.dimnames[1]))
    end
end
for op in (:Ac_mul_B, :At_mul_B)
    @eval ($op)(A::NamedMatrix, B::NamedMatrix) = NamedArray(($op)(A.array, B.array), (A.dicts[2], B.dicts[2]), (A.dimnames[2], B.dimnames[2]))
    for T in [StridedMatrix, AbstractMatrix] ## v0.4 ambiguity-hell
        @eval ($op)(A::NamedMatrix, B::$T) = NamedArray(($op)(A.array, B), (A.dicts[2], defaultnamesdict(size(B,2))), (A.dimnames[2], :B))
        @eval ($op)(A::$T, B::NamedMatrix) = NamedArray(($op)(A, B.array), (defaultnamesdict(size(A,2)), B.dicts[2]), (:A, B.dimnames[2]))
    end
end
for op in (:Ac_mul_Bc, :At_mul_Bt)
    @eval ($op)(A::NamedMatrix, B::NamedMatrix) = NamedArray(($op)(A.array, B.array), (A.dicts[2], B.dicts[1]), (A.dimnames[2], B.dimnames[1]))
    for T in [StridedMatrix, AbstractMatrix] ## v0.4 ambiguity-hell
        @eval ($op)(A::NamedMatrix, B::$T) = NamedArray(($op)(A.array, B), (A.dicts[2], defaultnamesdict(size(B,1))), (A.dimnames[2], :B))
        @eval ($op)(A::$T, B::NamedMatrix) = NamedArray(($op)(A, B.array), (defaultnamesdict(size(A,2)), B.dicts[1]), (:A, B.dimnames[1]))
    end
end

import Base.LinAlg: Givens, BlasFloat, lufact!, LU, ipiv2perm, cholfact!, cholfact, qrfact!, qrfact, eigfact!, eigvals!,
    hessfact, hessfact!, schurfact!, schurfact, svdfact!, svdfact, svdvals!, svdvals, svd, diag, diagm, scale!, scale,
    cond, kron, linreg, lyap, sylvester, isposdef

## matmul
## ambiguity, this can somtimes be a pain to resolve...
*{Tx,TiA,Ty}(x::SparseMatrixCSC{Tx,TiA},y::NamedMatrix{Ty}) = x*y.array
*{Tx,S,Ty}(x::SparseMatrixCSC{Tx,S},y::NamedVector{Ty}) = x*y.array
for t in (:Tridiagonal, :AbstractTriangular, :Givens, :Bidiagonal)
    @eval *(x::$t, y::NamedMatrix) = NamedArray(x*y.array, ([string(i) for i in 1:size(x,1)],names(y,2)), y.dimnames)
    @eval *(x::$t, y::NamedVector) = x*y.array
end

## There is no such thing as a A_mul_B
## Named * Named
*(A::NamedMatrix, B::NamedMatrix) = NamedArray(A.array * B.array, (A.dicts[1], B.dicts[2]), (A.dimnames[1], B.dimnames[2]))
*(A::NamedMatrix, B::NamedVector) = NamedArray(A.array * B.array, (A.dicts[1],), (B.dimnames[1],))
## Named * Abstract
*(A::NamedMatrix, B::AbstractMatrix) = NamedArray(A.array * B, (A.dicts[1], defaultnamesdict(size(B,2))), A.dimnames)
*(A::AbstractMatrix, B::NamedMatrix) = NamedArray(A * B.array, (defaultnamesdict(size(A,1)), B.dicts[2]), B.dimnames)
*(A::NamedMatrix, B::AbstractVector) = NamedArray(A.array * B, (A.dicts[1],), (A.dimnames[1],))
*(A::AbstractMatrix, B::NamedVector) = A * B.array

if false

Base.Ac_mul_Bc!(A::Matrix, B::NamedMatrix, C::Matrix) = Ac_mul_Bc!(A, B.array, C)

end

## \ --- or should we overload A_div_B?
## Named \ Named
\(x::NamedVector, y::NamedVector) = x.array \ y.array
\(x::NamedMatrix, y::NamedVector) = NamedArray(x.array\y.array, (names(x,2),), (x.dimnames[2],))
\(x::NamedVector, y::NamedMatrix) = NamedArray(x.array\y.array, (["1"],names(y,2)), (:A, y.dimnames[2]))
\(x::NamedMatrix, y::NamedMatrix) = NamedArray(x.array\y.array, (names(x,2),names(y,2)), (x.dimnames[2], y.dimnames[2]))

## Named \ Abstract
\(x::NamedVector, y::AbstractVecOrMat) = x.array \ y
\(x::NamedMatrix, y::AbstractVector) = NamedArray(x.array \ y, (x.dicts[2],), (x.dimnames[2],))
\(x::NamedMatrix, y::AbstractMatrix) = NamedArray(x.array \ y, (names(x,2),[string(i) for i in 1:size(y,2)]), (x.dimnames[2],:B))
## Abstract \ Named
## ambiguity
\{Tx<:Number,Ty<:Number}(x::Diagonal{Tx}, y::NamedVector{Ty}) = x \ y.array
@compat \{Tx<:Number,Ty<:Number}(x::Union{Bidiagonal{Tx},AbstractTriangular{Tx}}, y::NamedVector{Ty}) = x \ y.array
@compat \{Tx<:Number,Ty<:Number}(x::Union{Bidiagonal{Tx},AbstractTriangular{Tx}}, y::NamedMatrix{Ty}) = NamedArray(x \ y.array, ([string(i) for i in 1:size(x,2)], names(y,2)), (:A, y.dimnames[2]))
if VERSION >= v"0.4.0-dev"
    \(x::Bidiagonal,y::NamedVector) = NamedArray(x \ y.array, ([string(i) for i in 1:size(x,2)], names(y,2)), (:A, y.dimnames[2]))
    \(x::Bidiagonal,y::NamedMatrix) = NamedArray(x \ y.array, ([string(i) for i in 1:size(x,2)], names(y,2)), (:A, y.dimnames[2]))
end
## AbstractVectorOrMat gives us more ambiguities than separate entries...
\(x::AbstractVector, y::NamedVector) = x \ y.array
\(x::AbstractMatrix, y::NamedVector) = x \ y.array
\(x::AbstractVector, y::NamedMatrix) = NamedArray(x \ y.array, (["1"],names(y,2)), (:A, y.dimnames[2]))
\(x::AbstractMatrix, y::NamedMatrix) = NamedArray(x \ y.array, ([string(i) for i in 1:size(x,2)], names(y,2)), (:A, y.dimnames[2]))

## keeping names for some matrix routines
for f in (:inv, :chol, :sqrtm, :pinv, :expm)
    eval(Expr(:import, :Base, f))
    @eval ($f)(n::NamedArray) = NamedArray(($f)(n.array), n.dicts, n.dimnames)
end

## tril, triu
Base.tril!(n::NamedMatrix, k::Integer) = (tril!(n.array, k); n)
Base.triu!(n::NamedMatrix, k::Integer) = (triu!(n.array, k); n)
Base.triu(n::NamedMatrix, k::Integer) = triu!(copy(n), k)
Base.tril(n::NamedMatrix, k::Integer) = tril!(copy(n), k)

## LU factorization
function lufact!{T}(n::NamedArray{T}; pivot=true)
    luf = lufact!(n.array; pivot=pivot)
    LU{T,typeof(n),}(n, luf.ipiv, luf.info)
end

## after lu.jl, this could be merged at Base.
function Base.getindex{T,DT,AT}(A::LU{T,NamedArray{T,2,AT,DT}}, d::Symbol)
    m, n = size(A)
    if d == :L
        L = tril!(A.factors[1:m, 1:min(m,n)])
        for i = 1:min(m,n); L[i,i] = one(T); end
        setnames!(L, defaultnames(L,2), 2)
        setdimnames!(L, :LU, 2)
        return L
    end
    if d == :U
        U = triu!(A.factors[1:min(m,n), 1:n])
        setnames!(U, defaultnames(U,1), 1)
        setdimnames!(U, :LU, 1)
        return U
    end
    d == :p && return ipiv2perm(A.ipiv, m)
    if d == :P
        p = A[:p]
        P = zeros(T, m, m)
        for i in 1:m
            P[i,p[i]] = one(T)
        end
        return P
    end
    throw(KeyError(d))
end

## TODO: wait until Cholesky and CholeskyPivoted contain abstractmatrix types
function cholfact!{T<:BlasFloat}(n::NamedArray{T}, uplo::Symbol=:U; pivot=false, tol=0.0)
    uplochar = string(uplo)[1]
    if pivot
        A, piv, rank, info = LAPACK.pstrf!(uplochar, n.array, tol)
        return CholeskyPivoted{T}(A, uplochar, piv, rank, tol, info)
    else
        C, info = LAPACK.potrf!(uplochar, n.array)
        return info==0 ?  Cholesky(C, uplochar) : throw(PosDefException(info))
    end
end

cholfact{T<:BlasFloat}(n::NamedArray{T}, uplo::Symbol=:U; pivot=false, tol=0.0) = cholfact!(copy(n), uplo, pivot=pivot, tol=tol)

## ldlt skipped

## TODO: Wait for factorization.jl to change type of QR: AbstractMatrix
## from factorization
qrfact!{T<:BlasFloat}(n::NamedMatrix{T}; pivot=false) = pivot ? QRPivoted{T}(LAPACK.geqp3!(n.array)...) : QRCompactWY(LAPACK.geqrt!(n.array, min(minimum(size(A)), 36))...)
qrfact{T<:Base.LinAlg.BlasFloat}(n::NamedMatrix{T}; pivot=false) = qrfact!(copy(n.array),pivot=pivot)

eigfact!(n::NamedMatrix; permute::Bool=true, scale::Bool=true) = eigfact!(n.array, permute=permute, scale=scale)

eigvals!(n::NamedMatrix; permute::Bool=true, scale::Bool=true) = eigvals!(n.array, permute=permute, scale=scale)

hessfact!(n::NamedMatrix) = hessfact!(n.array)

schurfact!(n::NamedMatrix) = schurfact!(n.array)
schurfact(n::NamedMatrix) = schurfact!(copy(n))
schurfact(A::NamedMatrix, B::AbstractMatrix) = schurfact!(A.array, B)

svdfact!(n::NamedMatrix; thin::Bool=true) = svdfact!(n.array; thin=thin)
svdfact{T<:BlasFloat}(A::NamedMatrix{T};thin=true) = svdfact!(copy(A),thin=thin)

svdvals!(n::NamedArray) = svdvals!(n.array)
svdvals(n::NamedArray) = svdvals(n.array)

diag(n::NamedMatrix) = NamedArray(diag(n.array), n.dicts[1:1], n.dimnames[1:1])

diagm(n::NamedVector) = NamedArray(diagm(n.array), n.dicts[[1,1]], n.dimnames[[1,1]])

scale!(C::NamedMatrix, b::AbstractVector, A::AbstractMatrix) = (scale!(C.array, b, A); C)
scale(A::NamedMatrix, b::AbstractVector) = NamedArray(scale(A.array, b), A.dicts, A.dimnames)
scale(b::AbstractVector, A::NamedMatrix) = NamedArray(scale(b, A.array), A.dicts, A.dimnames)

# rank, vecnorm, norm, condskeel, trace, det, logdet OK
cond(n::NamedArray) = cond(n.array)

# null(n::NamedArray) = null(n.array)

function kron(a::NamedArray, b::NamedArray)
    n = Array(typeof(AbstractString[]), 2)
    dn = AbstractString[]
    for dim in 1:2
        n[dim] = AbstractString[]
        for i in names(a, dim)
            for j in names(b, dim)
                push!(n[dim], string(i, "×", j))
            end
        end
        push!(dn, string(dimnames(a,dim), "×", dimnames(b,dim)))
    end
    NamedArray(kron(a.array, b.array), tuple(n...), tuple(dn...))
end

linreg(x::NamedVecOrMat, y::AbstractVector) = linreg(x.array, y)
linreg(x::NamedVecOrMat, y::AbstractVector, w::AbstractVector) = linreg(x.array, y, w)

lyap(A::NamedMatrix, C::AbstractMatrix) = NamedArray(lyap(A.array,C), A.dicts, A.dimnames)

sylvester(A::NamedMatrix, B::AbstractMatrix, C::AbstractMatrix) = NamedArray(sylvester(A.array, B, C), A.dicts, A.dimnames)

## issym, istriu, istril OK
isposdef(n::NamedArray) = isposdef(n.array)

## eigs OK


text:
module GMHRunnersTest

import GeneralizedMetropolisHastings:
    _prop2seg,AbstractMHRunner,MHRuntimePolicy,numaccepted,numproposed

import Base.Test: @test,@test_approx_eq

function gettwoprocessnumbers(p::Int)
    if p == 1
        return [1,1]
    elseif p == 2
        w = workers()[1]
        return [w,w]
    end
    return workers()[1:2] #if there are many workers
end

function testfrom(state,vals,ll)
    @test state.from.values == vals
    @test state.from.loglikelihood == ll
end

function testproposals(state,vals,ll)
    @test state.proposals.values == vals
    @test state.proposals.loglikelihood == ll
end

function teststore(chain,vals,ll,i,a,p)
    @test chain.values[:,i] == vals
    @test chain.loglikelihood[i] == ll
    @test numaccepted(chain) == a
    @test numproposed(chain) == p
end

function setfrom(state,vals,loglikelihood)
    copy!(state.from.values,vals)
    copy!(state.from.loglikelihood,loglikelihood)
end

function setproposals(state,vals,loglikelihood)
    copy!(state.proposals.values,vals)
    copy!(state.proposals.loglikelihood,loglikelihood)
end

function setindicator(indicator,probs,samples)
    copy!(indicator.stationary,probs)
    copy!(indicator.samples,samples)
end

function setsegment(rsegs,from,ll1,prop,ll2,j)
    l = length(prop)
    p,s = _prop2seg(rsegs,j)
    r = rsegs.remote[s]
    @sync begin
        @spawnat r.where copy!(fetch(r).samplerstate.from.values,from)
        @spawnat r.where copy!(fetch(r).samplerstate.from.loglikelihood,ll1)
        @spawnat r.where copy!(fetch(r).samplerstate.proposals.values,l*(p-1)+1,prop,1,l)
        @spawnat r.where copy!(fetch(r).samplerstate.proposals.loglikelihood,p,ll2,1,1)
    end
    rsegs.prop2collected[j] = (p,s)
end

function testtransitionprobability(states,indicator)
    acc = Float64[]
    for x in states
        acc = vcat(acc,acceptanceratio(x))
    end
    acc = vcat(acc,[0.0])
    acc = exp(acc - maximum(acc))
    acc /= sum(acc)
    @test_approx_eq indicator.stationary acc
end

function testupdatefrom(indicatorstate,states,indicator)
    indend = indicatorsamples(indicator)[end]
    if indend != numproposals(indicator) + 1
        s,t = ind2sub(GeneralizedMetropolisHastings._ind2subtuple(states),indend)
        @test from(indicatorstate).values == proposals(states[t]).values[:,s]
        @test from(indicatorstate).loglikelihood[1] == proposals(states[t]).loglikelihood[s]
        @test from(indicatorstate).logprior[1] == proposals(states[t]).logprior[s]
    end
end

function dumpiterate1(samplerstates,indicator)
    dump(samplerstates[1])
    println("===============================================")
    dump(indicator1)
    println("+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++")
end

function runiterateandstore1!(runner,model,samplerstates,indicator,chain)
    prevacc = numaccepted(chain)
    prevprop = numproposed(chain)
    iterateandstore!(runner,model,samplerstates,indicator,chain)
    teststore(samplerstates[1],samplerstates,indicator,chain,prevacc,prevprop)
    testupdatefrom(samplerstates1[1],samplerstates1,indicator1)
    dumpiterate1(samplerstates,indicator)
end

function dumpiterate2(indicatorstate,samplerstates,indicator)
    dump(indicatorstate)
    println("000000000000000000000000000000000000000000000000000000000000000000000000000")
    dump(samplerstates[1])
    println("111111111111111111111111111111111111111111111111111111111111111111111111111")
    dump(samplerstates[2])
    println("222222222222222222222222222222222222222222222222222222222222222222222222222")
    dump(samplerstates[3])
    println("333333333333333333333333333333333333333333333333333333333333333333333333333")
    dump(indicator)
    println("+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++")
end

function runiterateandstore2!(runner,model,indicatorstate,samplerstates,indicator,chain)
    prevacc = numaccepted(chain)
    prevprop = numproposed(chain)
    iterateandstore!(runner,model,indicatorstate,samplerstates,indicator,chain)
    teststore(indicatorstate,samplerstates,indicator,chain,prevacc,prevprop)
    testupdatefrom(indicatorstate,samplerstates,indicator)
    dumpiterate2(indicatorstate,samplerstates,indicator)
end

type TestMHRunner <: AbstractMHRunner
    numburnin::Int
    numiterations::Int
    policy::MHRuntimePolicy
end

end


text:
module FHIcicle

using ..Ahorn, Maple

@mapdef Entity "FrostHelper/Icicle" Icicle(x::Integer, y::Integer, directory::String="objects/FrostHelper/icicle/", speed::Number=15.0, breakSfx::String="event:/game/09_core/iceball_break")

const placements = Ahorn.PlacementDict(
    "Icicle (Frost Helper)" => Ahorn.EntityPlacement(
        Icicle
    )
)

function Ahorn.selection(entity::Icicle)
    x, y = Ahorn.position(entity)
    return Ahorn.getSpriteRectangle("$(get(entity.data, "directory", "objects/FrostHelper/icicle/"))idle00", x, y)
end

Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::Icicle, room::Maple.Room) = Ahorn.drawSprite(ctx, "$(get(entity.data, "directory", "objects/FrostHelper/icicle/"))idle00", 0, 0)

end

max_stars_repo_name: JaThePlayer/FrostHelper

max_issues_repo_name: JaThePlayer/FrostHelper

max_forks_repo_name: JaThePlayer/FrostHelper