text:
#
# Upsampling
#

# GPU based bilinear upsampling including its gradient
#
# Based on the Caffe2 implementation at:
# The code is a translation from the following files:
# - https://github.com/pytorch/pytorch/blob/v1.8.0-rc1/caffe2/operators/upsample_op.cu
# - https://github.com/pytorch/pytorch/blob/v1.8.0-rc1/caffe2/core/common_gpu.h
#
# Copyright (c) 2016-2021 Facebook Inc.
# Copyright (c) 2015 Google Inc.
# Copyright (c) 2015 Yangqing Jia
# Copyright 2019-2020 Kakao Brain
#
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without modification, are
# permitted provided that the following conditions are met:
#
# 1. Redistributions of source code must retain the above copyright notice, this list of
#    conditions and the following disclaimer.
#
# 2. Redistributions in binary form must reproduce the above copyright notice, this list of
#    conditions and the following disclaimer in the documentation and/or other materials
#    provided with the distribution.
#
# 3. Neither the names of Facebook, Deepmind Technologies, NYU, NEC Laboratories America and
#    IDIAP Research Institute nor the names of its contributors may be used to endorse or
#    promote products derived from this software without specific prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY
# EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
# MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
# COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
# EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
# SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
# HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR
# TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
# SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

# Forward and backward pass have been tested to produce the same output
# as pytorch with align_corners=True - it works modulo bit noise.

@inline function compute_source_index(ratio::T, dst_index, align_corners) where T
    if align_corners
        return ratio*dst_index
    else
        src_idx = ratio * (dst_index + T(0.5)) - T(0.5)
        return max(zero(T), src_idx)
    end
end

###########
# linear
###########
function upsample_linear_wcn_kernel!(n_elem, rwidth, x, y, align_corners)
    index = (threadIdx().x - 1) + (blockIdx().x - 1) * blockDim().x

    if index < n_elem
        in_w, channels, batchsize = size(x)
        out_w, _, _ = size(y)

        ow = index % out_w

        # real_index = rwidth*ow
        real_index = compute_source_index(rwidth, ow, align_corners)
        iw0 = Base.floor(Int, real_index)
        offset = (iw0 < in_w-1) ? 1 : 0
        iw1 = iw0 + offset + 1
        w1lambda = real_index - iw0
        w0lambda = 1 - w1lambda
        iw0 += 1

        @inbounds for n in 1:batchsize
            for c in 1:channels
                val = (w0lambda * x[iw0, c, n]  + # w0 * i00
                       w1lambda * x[iw1, c, n])   # w1 * i01
                y[ow+1, c, n] = val
            end
        end
    end
    return nothing
end

# Δ is the gradient backpropagated from downstream layers
function ∇upsample_linear_wcn_kernel!(n_elem, rwidth, Δ, dx, align_corners)
    index = (threadIdx().x - 1) + (blockIdx().x - 1) * blockDim().x

    if index < n_elem
        in_width, channels, batchsize = size(Δ)
        out_width, _, _ = size(dx)

        iw = index % in_width

        # real_index_w = rwidth * iw
        real_index_w = compute_source_index(rwidth, iw, align_corners)
        ow0 = Base.floor(Int, real_index_w)
        offset = (ow0 < out_width - 1) ? 1 : 0
        ow1 = ow0 + offset + 1
        w1lambda = real_index_w - ow0
        w0lambda = 1 - w1lambda
        ow0 += 1

        @inbounds for n in 1:batchsize
            for c in 1:channels
                val = Δ[iw+1, c, n]
                CUDA.@atomic dx[ow0, c, n] += w0lambda * val
                CUDA.@atomic dx[ow1, c, n] += w1lambda * val
            end
        end
    end # if
    return nothing
end

function NNlib.upsample_linear_wcn!(y::CuArray{T,3}, x::CuArray{T,3}; align_corners=true) where T
    out_size = size(y)[1] # w

    if align_corners
        ratios = ntuple(i -> T((size(x,i)-1) / (size(y,i)-1)), 1)
    else
        ratios = ntuple(i -> T(size(x,i) / size(y,i)), 1)
    end

    kernel = @cuda launch=false upsample_linear_wcn_kernel!(out_size, ratios..., x, y, align_corners)
    config = launch_configuration(kernel.fun; max_threads=256)
    threads = Base.min(out_size, config.threads)
    blocks = cld(out_size, threads)
    kernel(out_size, ratios..., x, y, align_corners; threads=threads, blocks=blocks)
    return y
end

function NNlib.∇upsample_linear_wcn!(dx::CuArray{T,3}, Δ::CuArray{T,3}; align_corners=true) where T
    in_size = size(Δ)[1]
    if align_corners
        ratios = ntuple(i -> T((size(dx,i)-1) / (size(Δ,i)-1)), 1)  # reversed compared to forward pass
    else
        ratios = ntuple(i -> T(size(dx,i) / size(Δ,i)), 1)
    end

    kernel = @cuda launch=false ∇upsample_linear_wcn_kernel!(in_size, ratios..., Δ, dx, align_corners)
    config = launch_configuration(kernel.fun; max_threads=256)
    threads = Base.min(in_size, config.threads)
    blocks = cld(in_size, threads)
    kernel(in_size, ratios..., Δ, dx, align_corners; threads=threads, blocks=blocks)
    return dx
end

###########
# bilinear
###########
function upsample_bilinear_whcn_kernel!(n_elem, rwidth, rheight, x, y, align_corners)
    index = (threadIdx().x - 1) + (blockIdx().x - 1) * blockDim().x

    if index < n_elem
        in_w, in_h, channels, batchsize = size(x)
        out_w, out_h, _, _ = size(y)

        ow = index % out_w
        oh = index ÷ out_w

        # real_index = rheight*oh
        real_index = compute_source_index(rheight, oh, align_corners)
        ih0 = Base.floor(Int, real_index)
        offset = (ih0 < in_h-1) ? 1 : 0
        ih1 = ih0 + offset + 1
        h1lambda = real_index - ih0
        h0lambda = 1 - h1lambda
        ih0 += 1

        # real_index = rwidth*ow
        real_index = compute_source_index(rwidth, ow, align_corners)
        iw0 = Base.floor(Int, real_index)
        offset = (iw0 < in_w-1) ? 1 : 0
        iw1 = iw0 + offset + 1
        w1lambda = real_index - iw0
        w0lambda = 1 - w1lambda
        iw0 += 1

        @inbounds for n in 1:batchsize
            for c in 1:channels
                val = h0lambda * (w0lambda * x[iw0, ih0, c, n]  + # h0 * w0 * i00
                                  w1lambda * x[iw1, ih0, c, n]) + # h0 * w1 * i01
                      h1lambda * (w0lambda * x[iw0, ih1, c, n]  + # h1 * w0 * i10
                                  w1lambda * x[iw1, ih1, c, n])   # h1 * w1 * i11
                y[ow+1, oh+1, c, n] = val
            end
        end
    end
    return nothing
end

# Δ is the gradient backpropagated from downstream layers
function ∇upsample_bilinear_whcn_kernel!(n_elem, rwidth, rheight, Δ, dx, align_corners)
    index = (threadIdx().x - 1) + (blockIdx().x - 1) * blockDim().x

    if index < n_elem
        in_width, in_height, channels, batchsize = size(Δ)
        out_width, out_height, _, _ = size(dx)

        iw = index % in_width
        ih = index ÷ in_width

        # Compute Y axis lambdas
        # real_index_h = rheight*ih
        real_index_h = compute_source_index(rheight, ih, align_corners)
        oh0 = Base.floor(Int, real_index_h)
        offset = (oh0 < out_height-1) ? 1 : 0
        oh1 = oh0 + offset + 1
        h1lambda = real_index_h - oh0
        h0lambda = 1 - h1lambda
        oh0 += 1

        # # Compute X axis lambdas
        # real_index_w = rwidth * iw
        real_index_w = compute_source_index(rwidth, iw, align_corners)
        ow0 = Base.floor(Int, real_index_w)
        offset = (ow0 < out_width - 1) ? 1 : 0
        ow1 = ow0 + offset + 1
        w1lambda = real_index_w - ow0
        w0lambda = 1 - w1lambda
        ow0 += 1

        @inbounds for n in 1:batchsize
            for c in 1:channels
                val = Δ[iw+1, ih+1, c, n]
                CUDA.@atomic dx[ow0, oh0, c, n] += h0lambda * w0lambda * val
                CUDA.@atomic dx[ow1, oh0, c, n] += h0lambda * w1lambda * val
                CUDA.@atomic dx[ow0, oh1, c, n] += h1lambda * w0lambda * val
                CUDA.@atomic dx[ow1, oh1, c, n] += h1lambda * w1lambda * val
            end
        end
    end # if
    return nothing
end

function NNlib.upsample_bilinear_whcn!(y::CuArray{T,4}, x::CuArray{T,4}; align_corners=true) where T
    out_size = prod(size(y)[1:2]) # w*h

    if align_corners
        ratios = ntuple(i -> T((size(x,i)-1) / (size(y,i)-1)), 2)
    else
        ratios = ntuple(i -> T(size(x,i) / size(y,i)), 2)
    end

    kernel = @cuda launch=false upsample_bilinear_whcn_kernel!(out_size, ratios..., x, y, align_corners)
    config = launch_configuration(kernel.fun; max_threads=256)
    threads = Base.min(out_size, config.threads)
    blocks = cld(out_size, threads)
    kernel(out_size, ratios..., x, y, align_corners; threads=threads, blocks=blocks)
    return y
end

function NNlib.∇upsample_bilinear_whcn!(dx::CuArray{T,4}, Δ::CuArray{T,4}; align_corners=true) where T
    in_size = prod(size(Δ)[1:2])
    if align_corners
        ratios = ntuple(i -> T((size(dx,i)-1) / (size(Δ,i)-1)), 2)  # reversed compared to forward pass
    else
        ratios = ntuple(i -> T(size(dx,i) / size(Δ,i)), 2)
    end

    kernel = @cuda launch=false ∇upsample_bilinear_whcn_kernel!(in_size, ratios..., Δ, dx, align_corners)
    config = launch_configuration(kernel.fun; max_threads=256)
    threads = Base.min(in_size, config.threads)
    blocks = cld(in_size, threads)
    kernel(in_size, ratios..., Δ, dx, align_corners; threads=threads, blocks=blocks)
    return dx
end

###########
# trilinear
###########
function upsample_trilinear_whdcn_kernel!(n_elem, rwidth, rheight, rdepth, x, y, align_corners)
    index = (threadIdx().x - 1) + (blockIdx().x - 1) * blockDim().x

    if index < n_elem
        in_w, in_h, in_d, channels, batchsize = size(x)
        out_w, out_h, out_d, _, _ = size(y)

        ow = (index % (out_w * out_h)) % out_w
        oh = (index % (out_w * out_h)) ÷ out_w
        od = index ÷ (out_w * out_h)

        # real_index = rwidth*ow
        real_index = compute_source_index(rwidth, ow, align_corners)
        iw0 = Base.floor(Int, real_index)
        offset = (iw0 < in_w-1) ? 1 : 0
        iw1 = iw0 + offset + 1
        w1lambda = real_index - iw0
        w0lambda = 1 - w1lambda
        iw0 += 1

        # real_index = rheight*oh
        real_index = compute_source_index(rheight, oh, align_corners)
        ih0 = Base.floor(Int, real_index)
        offset = (ih0 < in_h-1) ? 1 : 0
        ih1 = ih0 + offset + 1
        h1lambda = real_index - ih0
        h0lambda = 1 - h1lambda
        ih0 += 1

        # real_index = rdepth*od
        real_index = compute_source_index(rdepth, od, align_corners)
        id0 = Base.floor(Int, real_index)
        offset = (id0 < in_d-1) ? 1 : 0
        id1 = id0 + offset + 1
        d1lambda = real_index - id0
        d0lambda = 1 - d1lambda
        id0 += 1

        @inbounds for n in 1:batchsize
            for c in 1:channels
                val = d0lambda *
                         (h0lambda *
                           (w0lambda * x[iw0, ih0, id0, c, n]  +
                            w1lambda * x[iw1, ih0, id0, c, n]) +
                          h1lambda *
                           (w0lambda * x[iw0, ih1, id0, c, n]   +
                            w1lambda * x[iw1, ih1, id0, c, n])) +
                      d1lambda *
                         (h0lambda *
                           (w0lambda * x[iw0, ih0, id1, c, n]  +
                            w1lambda * x[iw1, ih0, id1, c, n]) +
                          h1lambda *
                           (w0lambda * x[iw0, ih1, id1, c, n] +
                            w1lambda * x[iw1, ih1, id1, c, n]))

                y[ow+1, oh+1, od+1, c, n] = val
            end
        end
    end
    return nothing
end

# Δ is the gradient backpropagated from downstream layers
function ∇upsample_trilinear_whdcn_kernel!(n_elem, rwidth, rheight, rdepth, Δ, dx, align_corners)
    index = (threadIdx().x - 1) + (blockIdx().x - 1) * blockDim().x

    if index < n_elem
        in_width, in_height, in_depth, channels, batchsize = size(Δ)
        out_width, out_height, out_depth, _, _ = size(dx)

        iw = (index % (in_height * in_width)) % in_width
        ih = (index % (in_height * in_width)) ÷ in_width
        id = index ÷  (in_height * in_width)

        real_index_w = compute_source_index(rwidth, iw, align_corners)
        ow0 = Base.floor(Int, real_index_w)
        offset = (ow0 < out_width - 1) ? 1 : 0
        ow1 = ow0 + offset + 1
        w1lambda = real_index_w - ow0
        w0lambda = 1 - w1lambda
        ow0 += 1

        real_index_h = compute_source_index(rheight, ih, align_corners)
        oh0 = Base.floor(Int, real_index_h)
        offset = (oh0 < out_height-1) ? 1 : 0
        oh1 = oh0 + offset + 1
        h1lambda = real_index_h - oh0
        h0lambda = 1 - h1lambda
        oh0 += 1

        real_index_d = compute_source_index(rdepth, id, align_corners)
        od0 = Base.floor(Int, real_index_d)
        offset = (od0 < out_depth-1) ? 1 : 0
        od1 = od0 + offset + 1
        d1lambda = real_index_d - od0
        d0lambda = 1 - d1lambda
        od0 += 1

        @inbounds for n in 1:batchsize
            for c in 1:channels
                val = Δ[iw+1, ih+1, id+1, c, n]
                CUDA.@atomic dx[ow0, oh0, od0, c, n] += w0lambda * h0lambda * d0lambda * val
                CUDA.@atomic dx[ow1, oh0, od0, c, n] += w1lambda * h0lambda * d0lambda * val
                CUDA.@atomic dx[ow0, oh1, od0, c, n] += w0lambda * h1lambda * d0lambda * val
                CUDA.@atomic dx[ow1, oh1, od0, c, n] += w1lambda * h1lambda * d0lambda * val

                CUDA.@atomic dx[ow0, oh0, od1, c, n] += w0lambda * h0lambda * d1lambda * val
                CUDA.@atomic dx[ow1, oh0, od1, c, n] += w1lambda * h0lambda * d1lambda * val
                CUDA.@atomic dx[ow0, oh1, od1, c, n] += w0lambda * h1lambda * d1lambda * val
                CUDA.@atomic dx[ow1, oh1, od1, c, n] += w1lambda * h1lambda * d1lambda * val
            end
        end
    end # if
    return nothing
end

function NNlib.upsample_trilinear_whdcn!(y::CuArray{T,5}, x::CuArray{T,5}; align_corners=true) where T
    out_size = prod(size(y)[1:3]) # w*h*d

    if align_corners
        ratios = ntuple(i -> T((size(x,i)-1) / (size(y,i)-1)), 3)
    else
        ratios = ntuple(i -> T(size(x,i) / size(y,i)), 3)
    end

    kernel = @cuda launch=false upsample_trilinear_whdcn_kernel!(out_size, ratios..., x, y, align_corners)
    config = launch_configuration(kernel.fun; max_threads=256)
    threads = Base.min(out_size, config.threads)
    blocks = cld(out_size, threads)
    kernel(out_size, ratios..., x, y, align_corners; threads=threads, blocks=blocks)
    return y
end

function NNlib.∇upsample_trilinear_whdcn!(dx::CuArray{T,5}, Δ::CuArray{T,5}; align_corners=true) where T
    in_size = prod(size(Δ)[1:3])

    if align_corners
        ratios = ntuple(i -> T((size(dx,i)-1) / (size(Δ,i)-1)), 3)  # reversed compared to forward pass
    else
        ratios = ntuple(i -> T(size(dx,i) / size(Δ,i)), 3)
    end

    kernel = @cuda launch=false ∇upsample_trilinear_whdcn_kernel!(in_size, ratios..., Δ, dx, align_corners)
    config = launch_configuration(kernel.fun; max_threads=256)
    threads = Base.min(in_size, config.threads)
    blocks = cld(in_size, threads)
    kernel(in_size, ratios..., Δ, dx, align_corners; threads=threads, blocks=blocks)
    return dx
end


text:
const _rng = create_generator()
atexit(() -> destroy_generator(_rng))

# uniform
"""Generate n uniformly distributed numbers"""
curand(rng::RNG, ::Type{Float32}, n::Int) = generate_uniform(rng, UInt(n))
curand(::Type{Float32}, n::Int) = curand(_rng, Float64, n)

curand(rng::RNG, ::Type{Float64}, n::Int) = generate_uniform_double(rng, UInt(n))
curand(::Type{Float64}, n::Int) = curand(_rng, Float64, n)
curand(n::Int) = curand(Float64, n)

# normal
"""
Generate n normally distributed numbers with specified mean and
standard deviation
"""
curandn(rng::RNG, ::Type{Float32}, n::Int, mean::Float32, stddev::Float32) =
    generate_normal(rng, UInt(n), mean, stddev)
curandn(::Type{Float32}, n::Int, mean::Float32, stddev::Float32) =
    curandn(_rng, Float32, n, mean, stddev)

curandn(rng::RNG, ::Type{Float64}, n::Int, mean::Float64, stddev::Float64) =
    generate_normal_double(rng, UInt(n), mean, stddev)
curandn(::Type{Float64}, n::Int, mean::Float64, stddev::Float64) =
    curandn(_rng, Float64, n, mean, stddev)
curandn(n::Int, mean::Float64, stddev::Float64) =
    curandn(Float64, n, mean, stddev)

# log-normal
"""
Generate n log-normally distributed numbers with specified mean and
standard deviation
"""
curand_logn(rng::RNG, ::Type{Float32}, n::Int, mean::Float32, stddev::Float32) =
    generate_log_normal(rng, UInt(n), mean, stddev)
curand_logn(::Type{Float32}, n::Int, mean::Float32, stddev::Float32) =
    curand_logn(_rng, Float32, n, mean, stddev)

curand_logn(rng::RNG, ::Type{Float64}, n::Int, mean::Float64, stddev::Float64) =
    generate_log_normal_double(rng, UInt(n), mean, stddev)
curand_logn(::Type{Float64}, n::Int, mean::Float64, stddev::Float64) =
    curand_logn(_rng, Float64, n, mean, stddev)
curand_logn(n::Int, mean::Float64, stddev::Float64) =
    curand_logn(Float64, n, mean, stddev)

# poisson
"""Generate n numbers according to Poisson distribution"""
curand_poisson(rng::RNG, n::Int, lambda::Float64) =
    generate_poisson(rng, UInt(n), lambda)
curand_poisson(n::Int, lambda::Float64) = curand_poisson(_rng, n, lambda)

ext: jl


text:
using RDatasets, DataArrays, Gadfly

plot(dataset("car", "Womenlf"), x=:HIncome, y=:Region, Geom.histogram2d)


text:
export Adam

immutable Adam <: SolverMethod
end

function make_solver_parameters(solver::Adam; kwargs...)
  p = make_solver_parameters(lr_policy=LRPolicy.Inv(0.003, 0.0001, 0.5),
                             beta1=0.9,
                             beta2=0.999,
                             epsilon=1e-8)
  merge(p, SolverParameters(kwargs))
end

validate_parameters(solver::Adam, params::SolverParameters) = begin
    validate_parameters(params, :lr_policy, :beta1, :beta2, :epsilon)
end

type AdamSolverState <: InternalSolverState
  param_states        :: Vector{LayerState}
  grad_1st_moment_est :: Vector{Vector{Blob}} # Exponentially weighted moving average - biased estimate of 1st moment of gradient
  grad_2nd_moment_est :: Vector{Vector{Blob}} # Exponentially weighted moving average - biased estimate of raw 2nd moment of gradient
  t                   :: Float64  # timestep since estimates initialized
  learning_rate       :: Float64
end

type AdamSolverStateSnapshot <: SolverStateSnapshot
    iter                :: Int
    obj_val             :: Float64
    grad_1st_moment_est :: Vector{Vector{Array}}
    grad_2nd_moment_est :: Vector{Vector{Array}}
    t                   :: Float64  # timestep since estimates initialized
    learning_rate       :: Float64
end

function blobs_clone(blobs::Vector{Vector{Blob}})
    out = Array{Vector{Array}}(length(blobs))
    for (i, vecblobs) in enumerate(blobs)
        out[i] = [Array{eltype(b)}(size(b)) for b in vecblobs]
        for (dst, b) in zip(out[i], vecblobs)
            copy!(dst, b)
        end
    end
    return out
end

snapshot(state::SolverState{AdamSolverState}) = begin
    AdamSolverStateSnapshot(state.iter, state.obj_val,
                            blobs_clone(state.internal.grad_1st_moment_est),
                            blobs_clone(state.internal.grad_2nd_moment_est),
                            state.internal.t,
                            state.internal.learning_rate)
end

function solver_state(net::Net, snapshot::AdamSolverStateSnapshot)
    i_state = AdamSolverState(snapshot.learning_rate, net)
    for i in 1:length(i_state.param_states)
        for (dst, src) in zip(i_state.grad_1st_moment_est[i],
                              snapshot.grad_1st_moment_est[i])
            copy!(dst, src)
        end
        for (dst, src) in zip(i_state.grad_2nd_moment_est[i],
                              snapshot.grad_2nd_moment_est[i])
            copy!(dst, src)
        end
    end
    i_state.t = snapshot.t
    SolverState(snapshot.iter, snapshot.obj_val, Dict(), i_state)
end

function solver_state(method::Adam, net::Net, params::SolverParameters)
    learning_rate = get_learning_rate(params[:lr_policy])
    SolverState(AdamSolverState(learning_rate, net))
end

AdamSolverState(learning_rate::Float64, net::Net) = begin
    param_states = updatable_layer_states(net)

    grad_1st_moment_est = Array{Vector{Blob}}(length(param_states))
    grad_2nd_moment_est = Array{Vector{Blob}}(length(param_states))

    for i = 1:length(param_states)
        layerstate = param_states[i]
        grad_1st_moment_est[i] = [make_zero_blob(net.backend, eltype(x.blob),size(x.blob)...) for x in layerstate.parameters]
        grad_2nd_moment_est[i] = [make_zero_blob(net.backend, eltype(x.blob),size(x.blob)...) for x in layerstate.parameters]
    end

    return AdamSolverState(param_states, grad_1st_moment_est, grad_2nd_moment_est, 0, learning_rate)
end

function shutdown(state::SolverState{AdamSolverState})
  map(x -> map(destroy, x), state.internal.grad_1st_moment_est)
  map(x -> map(destroy, x), state.internal.grad_2nd_moment_est)
end

function update(solver::Solver{Adam}, net::Net, state::SolverState{AdamSolverState})
  state.internal.t += 1
  state.internal.learning_rate = get_learning_rate(solver.params[:lr_policy], state)
  for i = 1:length(state.internal.param_states)
    layer_state   = state.internal.param_states[i]
    m       = state.internal.grad_1st_moment_est[i]
    v       = state.internal.grad_2nd_moment_est[i]
    for j = 1:length(layer_state.parameters)
      gradient  = layer_state.parameters[j].gradient
      data_type = eltype(m[j])
      # N.B. we are ignoring the parameter-specific learning rate multipliers
      # since they ought to adapt automatically.
      update_parameters!(net, solver.method,
                         state.internal.learning_rate,
                         solver.params[:epsilon],
                         solver.params[:beta1],
                         solver.params[:beta2],
                         m[j], v[j],
                         gradient,
                         layer_state.parameters[j].blob,
                         state.internal.t, data_type)
    end
  end
end

function update_parameters!(net::Net{CPUBackend}, method::Adam,
                            alpha::Float64, epsilon::Float64, beta1::Float64, beta2::Float64,
                            m, v, gradient, param_blob, t, data_type)

  # update biased gradient moment estimates, m and v
  # m_t <- beta1*m_{t-1} + (1-beta1)*g_t

  BLAS.scal!(length(m), convert(data_type, beta1), pointer(m.data), 1)
  BLAS.axpy!(length(m), convert(data_type, 1-beta1), pointer(gradient.data), 1, pointer(m.data), 1)

  # now we need g2 = g.^2
  # Are we better to preallocate (heavy memory usage) or allocate each time (slower?)
  g2 = gradient.data .^ 2 # TODO more efficient way of making a new buffer, copying and squaring?

  BLAS.scal!(length(v), convert(data_type, beta2), pointer(g2), 1)
  BLAS.axpy!(length(v), convert(data_type, 1-beta2), pointer(g2), 1, pointer(v.data), 1)

  # Correct bias and calculate effective stepsize for timestep t
  alpha_t = alpha * sqrt(1 - beta2^t) / (1 - beta1^t)
  # Update param <- param - alpha_t * m/(sqrt(v) + epsilon)
  delta = m.data ./ (sqrt.(v.data) .+ epsilon)
  BLAS.axpy!(length(param_blob), convert(data_type, -alpha_t), pointer(delta), 1, pointer(param_blob.data), 1)

end


text:
using ImageReconstruction
using Test

@testset "radon: pixel-driven" begin
    pixels = 129
    views = 200
    I = zeros(pixels, pixels)
    I[pixels÷2, pixels÷2] = 1
    θ = range(0, step = 2π / views, length = views)
    t = -100:100

    P = radon(I, θ, t)

    @test all(P[101, :] .== 1)
end


text:
module MLStyle

export @match, Many, Do, @data, @use, use, @used
export def_pattern_def_app_pattern_def_gapp_pattern_mk_pattern_mk_app_pattern_mk_gapp_pattern
export PatternUnsolvedException, InternalException, SyntaxError, UnknownExtension, @syntax_err
export @active

include("Err.jl")
using MLStyle.Err

include("Extension.jl")
using MLStyle.Extension

include("toolz.jl")

include("Render.jl")

include("MatchCore.jl")
using MLStyle.MatchCore

include("Infras.jl")
using MLStyle.Infras

include("Pervasives.jl")
using MLStyle.Pervasives

include("StandardPatterns.jl")
using MLStyle.StandardPatterns

include("DataType.jl")
using MLStyle.DataType

export @λ
"""
Lambda cases.
e.g.
map((@λ (1, x) -> x), xs)

(2, 3) |> @λ begin
    1 -> 2
    2 -> 7
    (a, b) -> a + b
end

# 5
"""
macro λ(cases)
    TARGET = mangle(__module__)
    @match cases begin
        :($a -> $(b...)) =>
                esc(quote
                    function ($TARGET, )
                        @match $TARGET begin
                            $a => begin $(b...) end
                        end
                    end
                end)

        Do(stmts=[]) &&
        :(begin $(Many(:($a -> $(b...)) && Do(push!(stmts, :($a => begin $(b...) end))) ||
                       (a :: LineNumberNode) && Do(push!(stmts , a))
                       )...)
          end) =>
            esc(quote
                function ($TARGET, )
                    @match $TARGET begin
                        $(stmts...)
                    end
                end
            end)
        _ => @syntax_err "Syntax error in lambda case definition. Check if your arrow is `->` but not `=>`!"
    end
end

export @when
function when(let_expr)
    @match let_expr begin
       Expr(:let, Expr(:block, bindings...) ||  a && Do(bindings = [a]), ret) =>
            foldr(bindings, init=ret) do each, last
                @match each begin
                    :($a = $b) =>
                        :(
                            $MLStyle.@match $b begin
                                $a =>  $last
                                _  =>  nothing
                            end
                        )
                    a => :(let $a; $last end)
                end
            end
       Expr(a, _...) => @syntax_err "Expect a let-binding, but found a `$a` expression."
       _ => @syntax_err "Expect a let-binding."
    end
end

"""

1. Allow destructuring in binding sequences of let syntax.

In binding sequences,
- For the bindings with the form `a = b`, you can use destructuring here.
- For others like `@inline f(x) = 1`, it's the same as the original let binding.

@when let (a, 1) = x,
          [b, c, 5] = y
        (a, b, c)
end

2. For a regular assignment, like
```
@when (a, 2) = x begin
    # dosomething
end
```
It's nothing different with
```
@match x begin
    (a, 2) => # dosomething
    _ => nothing
end
```
"""
macro when(let_expr)
    when(let_expr) |> esc
end

macro when(assignment, ret)
    @match assignment begin
        :($_ = $_) =>
            when(Expr(:let, Expr(:block, assignment), ret)) |> esc
        _ => @syntax_err "Not match the form of `@when a = b expr`"
    end
end

macro stagedexpr(exp)
    __module__.eval(exp)
end

end # module


text:
using LinearAlgebraicRepresentation
using LARVIEW

VOID = [[]], [[1]]

# example 1

V = [[0,0] [1,0] [2,0] [0,1] [1,1] [2,1] [0,2] [1,2] [2,2]];
FV = [[1,2,4],[2,3,5],[3,5,6],[4,5,7],[5,7,8],[6,8,9]];
pattern = repeat([1,2,-3],outer=4);
model = (V,FV);
W,FW = LinearAlgebraicRepresentation.extrudeSimplicial(model, pattern);
LARVIEW.view(W,FW)

# example 2

model = LinearAlgebraicRepresentation.extrudeSimplicial( VOID, ones(10) )
LARVIEW.view(model)
model = LinearAlgebraicRepresentation.extrudeSimplicial( model, ones(10) )
LARVIEW.view(model)
model = LinearAlgebraicRepresentation.extrudeSimplicial( model, ones(10) )
LARVIEW.view(model)

# example 3

model = LinearAlgebraicRepresentation.extrudeSimplicial( VOID, repeat([1,-1],outer=10) )
LARVIEW.view(model)
model = LinearAlgebraicRepresentation.extrudeSimplicial( model, repeat([1,-1],outer=10) )
LARVIEW.view(model)

# example 4

grid_2d = LinearAlgebraicRepresentation.simplexGrid([3,3])
LARVIEW.view(grid_2d)
grid_3d = LinearAlgebraicRepresentation.simplexGrid([2,3,4])
LARVIEW.view(grid_3d)
V,CV = LinearAlgebraicRepresentation.simplexGrid([1,1,1])
LARVIEW.view(V,CV)

# example 5

SK2 = LinearAlgebraicRepresentation.simplexFacets(CV)
LARVIEW.view(V, SK2)
SK1 = LinearAlgebraicRepresentation.simplexFacets(SK2)
LARVIEW.view(V, SK1)

max_stars_repo_name: EmaLoprevite/Simplexn

max_issues_repo_name: EmaLoprevite/Simplexn

max_forks_repo_name: EmaLoprevite/Simplexn


text:
function test_measures()
    m1 = Lebesgue()
    @test !isdiscrete(m1)
    @test iscontinuous(m1)
    @test !isnormalized(m1)
    @test domaintype(m1) == Float64
    @test DomainIntegrals.unsafe_weightfun(m1, 0.4) == 1

    m2 = DomainIntegrals.LebesgueUnit()
    @test !isdiscrete(m2)
    @test iscontinuous(m2)
    @test isnormalized(m2)
    @test domaintype(m2) == Float64
    @test support(m2) == UnitInterval()
    @test weightfun(m2, 0.4) == 1
    @test weightfun(m2, 1.4) == 0
    @test m2 == dx(UnitInterval())

    m3 = LegendreWeight()
    @test !isdiscrete(m3)
    @test iscontinuous(m3)
    @test !isnormalized(m3)
    @test domaintype(m3) == Float64
    @test support(m3) == ChebyshevInterval()
    @test weightfun(m3, 0.4) == 1
    @test weightfun(m3, -0.4) == 1
    @test weightfun(m3, 1.4) == 0
    @test weightfun(m3, -1.4) == 0
    @test jacobi_α(m3) == 0
    @test jacobi_β(m3) == 0

    x = 0.5
    m4 = DiracWeight(x)
    @test !isdiscrete(m4)
    @test iscontinuous(m4)
    @test isnormalized(m4)
    @test domaintype(m4) == Float64
    @test support(m4) == Point(x)
    @test point(m4) == x
    @test weightfun(m4, x) == Inf
    @test weightfun(m4, x+1) == 0
    @test weightfun(m4, big(x)) == Inf
    @test weightfun(m4, big(x+1)) == 0

    m5 = GaussianWeight{SVector{2,Float64}}()
    @test !isdiscrete(m5)
    @test iscontinuous(m5)
    @test isnormalized(m5)
    @test domaintype(m5) == SVector{2,Float64}
    @test support(m5) == FullSpace{SVector{2,Float64}}()
    @test weightfun(m5, SVector(0.0,0.0)) ≈ 1/(2pi)
    @test weightfun(m5, SVector(0,0)) ≈ 1/(2pi)
    @test weightfun(m5, SVector(big(0),big(0))) ≈ 1/(2pi)
    @test weightfun(m5, SVector(big(0.0),big(0.0))) ≈ 1/(2pi)

    m6 = LaguerreWeight(0.0)
    @test !isdiscrete(m6)
    @test iscontinuous(m6)
    @test isnormalized(m6)
    @test !isnormalized(LaguerreWeight(0.1))
    @test domaintype(m6) == Float64
    @test support(m6) == HalfLine()
    @test weightfun(m6, 0.4) == exp(-0.4)
    @test weightfun(m6, -0.4) == 0
    @test weightfun(m6, big(0.4)) == exp(-big(0.4))

    m7 = ChebyshevTWeight()
    @test support(m7) isa ChebyshevInterval
    @test !isdiscrete(m7)
    @test iscontinuous(m7)
    @test !isnormalized(m7)
    @test domaintype(m7) == Float64
    @test jacobi_α(m7) == -1/2
    @test jacobi_β(m7) == -1/2
    @test m7 == JacobiWeight(-1/2, -1/2)

    m8 = ChebyshevUWeight()
    @test support(m8) isa ChebyshevInterval
    @test !isdiscrete(m8)
    @test iscontinuous(m8)
    @test !isnormalized(m8)
    @test domaintype(m8) == Float64
    @test jacobi_α(m8) == 1/2
    @test jacobi_β(m8) == 1/2
    @test m8 == JacobiWeight(1/2, 1/2)
end


text:
using Documenter
import InterestRates

makedocs(
    sitename = "InterestRates.jl",
    modules = [ InterestRates ],
    pages = [ "Home" => "index.md",
              "API Reference" => "api.md" ]
)

deploydocs(
    repo = "github.com/felipenoris/InterestRates.jl.git",
    target = "build",
)


text:
function wavematrix3D_allocate(ω::T, medium::Medium{T}, species::Species{T};
        dim = 3, tol::T = 1e-4,
        basis_order::Int = maximum_basis_order(ω, medium, species; tol=tol),
        radius_multiplier::T = T(1.005),
        t_vecs::Vector{Vector{Complex{T}}} = t_vectors(ω, medium, species; basis_order = basis_order, dim = dim),
        kws...) where T<:Number

    k = real(ω/medium.c)
    S = length(species)
    ho = basis_order
    len = (ho+1)^4 * S
    MM_mat = Matrix{Complex{T}}(undef,len,len)

    as = radius_multiplier*[(s1.r + s2.r) for s1 in species, s2 in species]
    cs = reshape(
        [
            gaunt_coefficient(l1,m1,l2,m2,l3,m3)
        for l1 = 0:ho for m1 = -l1:l1 for l2 = 0:ho for m2 = -l2:l2 for l3 = 0:ho for m3 = -l3:l3]
    , ((ho+1)^2,(ho+1)^2,(ho+1)^2));

    # NOTE that cs[n3,n2,n1] == gaunt_coefficient(l1,m1,l2,m2,l3,m3)
    ls, ms = spherical_harmonics_indices(ho);
    n_max = length(ls)

    function M_component2(keff,Ns,n,dn,n1,n2,s1,s2)::Complex{T}
        (n == dn && n1 == n2 && s1 == s2 ? 1.0 : 0.0) +
        as[s1,s2] * species[s2].num_density * t_vecs[s1][ls[n]+ho+1] *
        sum(cs[:,n,dn] .* cs[:,n1,n2] .* Ns[ls[:].+1,s1,s2]) / (keff^2.0 - k^2.0)
    end
    function MM2(keff::Complex{T})
        Ns = [kernelN2D(l3,k*as[s1,s2],keff*as[s1,s2]; dim = dim) for l3 = 0:ho, s1 = 1:S, s2 = 1:S]
        ind2 = 1
        for s2 = 1:S for dn = 1:n_max for n1 = 1:n_max
            ind1 = 1
            for s1 = 1:S for n = 1:n_max for n2 = 1:n_max
                MM_mat[ind1, ind2] = M_component2(keff,Ns,n,dn,n1,n2,s1,s2)
                ind1 += 1
            end end end
            ind2 += 1
        end end end
        return MM_mat
    end
    # keff = 1.0 + rand()*im
    # s1 = s2 = 1

    # norm([
    #     M_component2(keff,Ns,n,dn,n1,n2,s1,s2) - M_component(keff,Ns,ls[n],ms[n],ls[n2],ms[n2],s1,ls[dn],ms[dn],ls[n1],ms[n1],s2)
    # for n = 1:n_max for n2 = 1:n_max for dn = 1:n_max for n1 = 1:n_max ])
    # #
    # @time [
    #     M_component2(keff,Ns,n,dn,n1,n2,s1,s2)
    # for n = 1:n_max for dn = 1:n_max for n1 = 1:n_max for n2 = 1:n_max ];
    # #
    # @time [
    #      M_component(keff,Ns,ls[n],ms[n],ls[n2],ms[n2],s1,ls[dn],ms[dn],ls[n1],ms[n1],s2)
    # for n = 1:n_max for n2 = 1:n_max for dn = 1:n_max for n1 = 1:n_max ];
    return MM2
end

max_stars_repo_name: arturgower/EffectiveWaves


text:
using SpecialFunctions
lgamma_(x) = logabsgamma(x)[1]

logsumexp(a,b) = (m = max(a,b); m == -Inf ? -Inf : log(exp(a-m) + exp(b-m)) + m)

function randp(p,k)
    s = 0.; for j = 1:k; s += p[j]; end
    u = rand()*s
    j = 1
    C = p[1]
    while u > C
        j += 1
        C += p[j]
    end
    @assert(j <= k)
    return j
end

function randlogp!(log_p,k)
    log_s = -Inf; for j = 1:k; log_s = logsumexp(log_s,log_p[j]); end
    p = log_p
    for j = 1:k; p[j] = exp(log_p[j]-log_s); end
    return randp(p,k)
end

function ordered_insert!(index,list,t)
    j = t
    while (j>0) && (list[j]>index)
        list[j+1] = list[j]
        j -= 1
    end
    list[j+1] = index
end

function ordered_remove!(index,list,t)
    for j = 1:t
        if list[j]>=index; list[j] = list[j+1]; end
    end
end

function ordered_next(list)
    j = 1
    while list[j]==j; j += 1; end
    return j
end

function restricted_Gibbs!(zsa,zsb,tia,tib,tja,tjb,cia,cib,cja,cjb,ni,nj,i,j,S,ns,x,b,H,active)
# NOTE: The sufficient statistics of tia and tja must be in sync with zsa.
# Also, note that the sufficient statistics of tib and tjb are not updated in this procedure.
    if !active; tia,tja = deepcopy(tia),deepcopy(tja); end
    log_p = 0.
    for ks = 1:ns
        k = S[ks]
        if k!=i && k!=j
            if zsa[k]==cia
                ni -= 1; Theta_remove!(tia,x[k])
            else
                nj -= 1; Theta_remove!(tja,x[k])
            end
            Li = log_marginal(x[k],tia,H) - log_marginal(tia,H)
            Lj = log_marginal(x[k],tja,H) - log_marginal(tja,H)
            Pi = exp(log(ni+b)+Li - logsumexp(log(ni+b)+Li,log(nj+b)+Lj))

            if active
                if rand()<Pi
                    zsb[k] = cib
                else
                    zsb[k] = cjb
                end
            end
            if zsb[k]==cib
                ni += 1
                Theta_adjoin!(tia,x[k])
                log_p += log(Pi)
            else
                nj += 1
                Theta_adjoin!(tja,x[k])
                log_p += log(1-Pi)
            end
        end
    end
    return log_p,ni,nj
end

function split_merge!(x,z,zs,S,theta,list,N,t,H,a,b,log_v,n_split,n_merge)
    n = length(x)

    # randomly choose a pair of indices
    i = round(Int,ceil(rand()*n))
    j = round(Int,ceil(rand()*(n-1))); if j>=i; j += 1; end
    ci0,cj0 = z[i],z[j]
    ti0,tj0 = theta[ci0],theta[cj0]

    # set S[1],...,S[ns] to the indices of the points in clusters ci0 and cj0
    ns = 0
    for k = 1:n
        if z[k]==ci0 || z[k]==cj0; ns += 1; S[ns] = k; end
    end

    # find available cluster IDs for merge and split parameters
    k = 1
    while list[k]==k; k += 1; end; cm = k
    while list[k]==k+1; k += 1; end; ci = k+1
    while list[k]==k+2; k += 1; end; cj = k+2
    tm,ti,tj = theta[cm],theta[ci],theta[cj]

    # merge state
    for ks = 1:ns; Theta_adjoin!(tm,x[S[ks]]); end # get the sufficient statistics

    # randomly choose the split launch state
    zs[i] = ci; Theta_adjoin!(ti,x[i]); ni = 1
    zs[j] = cj; Theta_adjoin!(tj,x[j]); nj = 1
    for ks = 1:ns  # start with a uniformly chosen split
        k = S[ks]
        if k!=i && k!=j
            if rand()<0.5; zs[k] = ci; Theta_adjoin!(ti,x[k]); ni += 1
            else;          zs[k] = cj; Theta_adjoin!(tj,x[k]); nj += 1
            end
        end
    end
    for rep = 1:n_split  # make several moves
        log_p,ni,nj = restricted_Gibbs!(zs,zs,ti,ti,tj,tj,ci,ci,cj,cj,ni,nj,i,j,S,ns,x,b,H,true)
    end

    # make proposal
    if ci0==cj0  # propose a split
        # make one final sweep and compute it's probability density
        log_prop_ab,ni,nj = restricted_Gibbs!(zs,zs,ti,ti,tj,tj,ci,ci,cj,cj,ni,nj,i,j,S,ns,x,b,H,true)

        # probability of going from merge state to original state
        log_prop_ba = 0.0  # log(1)

        # compute acceptance probability
        log_prior_b = log_v[t+1] + lgamma_(ni+b)+lgamma_(nj+b)-2*lgamma_(a)
        log_prior_a = log_v[t] + lgamma_(ns+b)-lgamma_(a)
        log_lik_ratio = log_marginal(ti,H) + log_marginal(tj,H) - log_marginal(ti0,H)
        p_accept = min(1.0, exp(log_prop_ba-log_prop_ab + log_prior_b-log_prior_a + log_lik_ratio))
        #println("split proposal: ",p_accept)

        # accept or reject
        if rand()<p_accept # accept split
            #  z, list, N, theta, t
            for ks = 1:ns; z[S[ks]] = zs[S[ks]]; end
            ordered_remove!(ci0,list,t)
            ordered_insert!(ci,list,t-1)
            ordered_insert!(cj,list,t)
            N[ci0],N[ci],N[cj] = 0,ni,nj
            t += 1
            Theta_clear!(ti0)
        else # reject split
            Theta_clear!(ti)
            Theta_clear!(tj)
        end
        Theta_clear!(tm)

    else  # propose a merge
        # probability of going to merge state
        log_prop_ab = 0.0  # log(1)

        # compute probability density of going from split launch state to original state
        log_prop_ba,ni,nj = restricted_Gibbs!(zs,z,ti,ti0,tj,tj0,ci,ci0,cj,cj0,ni,nj,i,j,S,ns,x,b,H,false)

        # compute acceptance probability
        log_prior_b = log_v[t-1] + lgamma_(ns+b)-lgamma_(a)
        log_prior_a = log_v[t] + lgamma_(ni+b)+lgamma_(nj+b)-2*lgamma_(a)
        log_lik_ratio = log_marginal(tm,H) - log_marginal(ti0,H) - log_marginal(tj0,H)
        p_accept = min(1.0, exp(log_prop_ba-log_prop_ab + log_prior_b-log_prior_a + log_lik_ratio))
        #println("merge proposal: ",p_accept)

        # accept or reject
        if rand()<p_accept # accept merge
            #  z, list, N, theta, t
            for ks = 1:ns; z[S[ks]] = cm; end
            ordered_remove!(ci0,list,t)
            ordered_remove!(cj0,list,t-1)
            ordered_insert!(cm,list,t-2)
            N[cm],N[ci0],N[cj0] = ns,0,0
            t -= 1
            Theta_clear!(ti0)
            Theta_clear!(tj0)
        else # reject merge
            Theta_clear!(tm)
        end
        Theta_clear!(ti)
        Theta_clear!(tj)
    end
    return t
end

function sampler(options,n_total,n_keep)
    x,n = options.x,options.n
    t_max = options.t_max
    a,b,log_v = options.a,options.b,options.log_v
    use_splitmerge,n_split,n_merge = options.use_splitmerge,options.n_split,options.n_merge
    use_hyperprior = options.use_hyperprior

    model_type = options.model_type
    alpha_random,alpha = options.alpha_random,options.alpha
    sigma_alpha = 0.1 # scale for MH proposals in alpha move

    @assert(n==length(x))
    keepers = zeros(Int,n_keep)
    keepers[:] = round.(Int,range(round(Int,n_total/n_keep),stop=n_total,length=n_keep))
    keep_index = 0

    t = 1  # number of clusters
    z = ones(Int,n)  # z[i] = the cluster ID associated with data point i
    list = zeros(Int,t_max+3); list[1] = 1  # list[1:t] = the list of active cluster IDs
                                            # list is maintained in increasing order for 1:t, and is 0 after that.
                                            # Further, maximum(list[1:t]) is always <= t_max+3
    c_next = 2  # an available cluster ID to be used next
    N = zeros(Int,t_max+3); N[1] = n  # N[c] = size of cluster c

    H = construct_hyperparameters(options)
    theta = [new_theta(H)::Theta for c = 1:t_max+3]  # theta[c] = parameters for cluster c

    for i = 1:n; Theta_adjoin!(theta[1],x[i]); end

    log_p = zeros(n+1)
    zs = ones(Int,n)  # temporary variable used for split-merge assignments
    S = zeros(Int,n)  # temporary variable used for split-merge indices

    log_Nb = log.((1:n) .+ b)

    # Record-keeping variables
    t_r = zeros(Int8,n_total); @assert(t_max < 2^7)
    N_r = zeros(Int16,t_max+3,n_total); @assert(n < 2^15)
    z_r = zeros(Int8,n,n_keep); @assert(t_max < 2^7)
    theta_r = Array{Theta}(undef,0,0)

    for iteration = 1:n_total
        #  -------------- Resample H --------------
        if use_hyperprior
            update_hyperparameters!(H,theta,list,t,x,z)
        end
        if model_type=="DPM" && alpha_random
            # Metropolis-Hastings move for DP concentration parameter (using p_alpha(a) = exp(-a) = Exp(a|1))
            aprop = alpha*exp(randn()*sigma_alpha)
            top = t*log(aprop) - lgamma_(aprop+n) + lgamma_(aprop) - aprop + log(aprop)
            bot = t*log(alpha) - lgamma_(alpha+n) + lgamma_(alpha) - alpha + log(alpha)
            if rand() < min(1.0,exp(top-bot))
                alpha = aprop
            end
            #for i=1:t_max+1; log_v[i] = i*log(alpha) - lgamma_(alpha+n) + lgamma_(alpha); end
            log_v = float(1:t_max+1)*log(alpha) .- lgamma_(alpha+n) .+ lgamma_(alpha)
        end

        # -------------- Resample z's --------------
        for i = 1:n
            # remove point i from it's cluster
            c = z[i]
            N[c] -= 1
            Theta_remove!(theta[c],x[i])
            if N[c]>0
                c_prop = c_next
            else
                c_prop = c
                # remove cluster {i}, keeping list in proper order
                ordered_remove!(c,list,t)
                t -= 1
            end

            # compute probabilities for resampling
            for j = 1:t; cc = list[j]
                log_p[j] = log_Nb[N[cc]] + log_marginal(x[i],theta[cc],H) - log_marginal(theta[cc],H)
            end
            log_p[t+1] = log_v[t+1]-log_v[t] + log(a) + log_marginal(x[i],theta[c_prop],H)

            # sample a new cluster for it
            j = randlogp!(log_p,t+1)

            # add point i to it's new cluster
            if j<=t
                c = list[j]
            else
                c = c_prop
                ordered_insert!(c,list,t)
                t += 1
                c_next = ordered_next(list)
                @assert(t<=t_max, "Sampled t has exceeded t_max. Increase t_max and retry.")
            end
            # update sufficient statistics
            Theta_adjoin!(theta[c],x[i])
            z[i] = c
            N[c] += 1
        end

        # -------------- Split/merge move --------------
        if use_splitmerge
            t = split_merge!(x,z,zs,S,theta,list,N,t,H,a,b,log_v,n_split,n_merge)
            c_next = ordered_next(list)
            @assert(t<=t_max, "Sampled t has exceeded t_max. Increase t_max and retry.")
        end

        # -------------- Record results --------------
        t_r[iteration] = t
        for j = 1:t
            N_r[list[j],iteration] = N[list[j]]
        end
        if iteration==keepers[keep_index+1]
            keep_index += 1
            for i = 1:n; z_r[i,keep_index] = z[i]; end
        end
    end

    return t_r,N_r,z_r,theta_r,keepers
end


text:
using Dates
using Distributions

@testset "Black Scholes" begin
    #start_time = now()
    #end_time = now()
    run_with_job("scholes") do job
        size = 512000000  # 256000000
	price = BanyanArrays.fill(4.0, size)
	strike = BanyanArrays.fill(4.0, size)
	t = BanyanArrays.fill(4.0, size)
	rate = BanyanArrays.fill(4.0, size)
	vol = BanyanArrays.fill(4.0, size)

	d1 = map(
		(p, s, t, r, v)->(
			(log(p / s) + (r + v ^ 2 * 0.5) * t) / (v * sqrt(t))
		),
		price,
		strike,
		t,
		rate,
		vol
	)

	d2 = map(
		(d1, v, t)->(
			d1 - (v * sqrt(t))
		),
		d1,
		vol,
		t
	)

	call = map(
		  (d1, d2, p, s, t, r)->(
			(cdf(Normal(), d1) * p) - (cdf(Normal(), d2) * s * exp(-r * t))
		  ),
		  d1,
		  d2,
		  price,
		  strike,
		  t,
		  rate
	)

	call_sum = sum(call)
	res = collect(call_sum)
        #res = collect(res)
	#end = now()
        #@test typeof(res) == Base.Vector{Float64}
        #@test all(v->v==3.999999985812889, res)
    end
    #end_time = now()
    #println(end_time - start_time)
end

#function run_bs(size::Integer)

    # c05 = Float64(3.0)
    # c10 = Float64(1.5)
    # invsqrt2 = 1.0 / sqrt(2.0)

    # # rsig = rate + (vol.^2) * c05
    # rsig = vol .^ 2 # TODO: Fix issue in FutureArray where operands are same
    # rsig = rsig .* c05
    # rsig = rsig .+ rate

    # # rsig

    # # vol_sqrt = vol .* sqrt.(t)
    # vol_sqrt = sqrt.(t)
    # vol_sqrt = vol_sqrt .* vol

    # # d1 = (log.(price ./ strike) + rsig .* t) ./ vol_sqrt
    # d1 = price ./ strike
    # d1 = log.(d1)
    # tmp = rsig .* t
    # d1 = d1 .+ tmp
    # d1 = d1 ./ vol_sqrt

    # # d1

    # d2 = d1 .- vol_sqrt

    # # d2

    # # d1 = c05 .+ c05 .* exp.(d1 .* invsqrt2)
    # # d2 = c05 .+ c05 .* exp.(d2 .* invsqrt2)
    # d1 = d1 * invsqrt2
    # d1 = exp.(d1)
    # d1 = d1 .* c05
    # d1 = d1 .+ c05
    # d2 = d2 .* invsqrt2
    # d2 = exp.(d2)
    # d2 = d2 .* c05
    # d2 = d2 .+ c05

    # # e_rt = exp.((-rate) .* t)
    # e_rt = -rate
    # e_rt = e_rt .* t
    # e_rt = exp.(e_rt)

    # # call = price .* d1 - e_rt .* strike .* d2
    # call = price .* d1
    # tmp = e_rt .* strike
    # tmp = tmp .* d2
    # call = call .- tmp

    # put = e_rt .* strike .* (c10 .- d2) - price .* (c10 .- d1)

    # evaluate(call)
#end

max_stars_repo_name: banyan-team/banyan-julia

max_issues_repo_name: banyan-team/banyan-julia

max_forks_repo_name: banyan-team/banyan-julia


text:
push!(LOAD_PATH, ".")

module CITypes # based on the CI5 interpreter

using Error
using Lexer
export parse, type_of_expr, NumType, BoolType, FunType, NListType

# ===================================================

abstract AE
abstract TypeVal
abstract TypeEnvironment

# ===================================================

# AST nodes for expressions

type Num <: AE
  n::Real
end

type Boolean <: AE
  v::Bool
end

type Plus <: AE
  lhs::AE
  rhs::AE
end

type Minus <: AE
  lhs::AE
  rhs::AE
end

type IsZero <: AE
  arg::AE
end

type IfB <: AE
  condition::AE
  true_branch::AE
  false_branch::AE
end

type With <: AE
  name::Symbol
  binding_expr::AE
  body::AE
end

type Id <: AE
  name::Symbol
end

type FunDef <: AE
  formal_parameter::Symbol
  formal_type::TypeVal
  fun_body::AE
end

type FunApp <: AE
  fun_expr::AE
  arg_expr::AE
end

type NEmpty <: AE
end

type NIsEmpty <: AE
  list::AE
end

type NCons <: AE
  f::AE
  r::AE
end

type NFirst <: AE
  list::AE
end

type NRest <: AE
  list::AE
end

# ===================================================

# AST nodes for types (used for type information as well)

type NumType <: TypeVal
end

type BoolType <: TypeVal
end

type FunType <: TypeVal
  arg_type::TypeVal
  result_type::TypeVal
end

type NListType <: TypeVal
end

# ===================================================

# Type Environments

type mtTypeEnvironment <: TypeEnvironment
end

type CTypeEnvironment <: TypeEnvironment
  name::Symbol
  value::TypeVal
  parent::TypeEnvironment
end

# ===================================================

# Parser for expressions
# functional for valid input, doesn't fully reject bad input)

function parse( expr::Real )
  return Num( expr )
end

function parse( expr::Bool )
  return Boolean( expr )
end

function parse( expr::Symbol )
  if expr == :nempty
    return NEmpty()
  else
    return Id( expr )
  end
end

function parse( expr::Array{Any} )

  op_symbol = expr[1]

  if op_symbol == :+
    lhs = parse( expr[2] )
    rhs = parse( expr[3] )
    return Plus( lhs, rhs )

  elseif op_symbol == :-
    lhs = parse( expr[2] )
    rhs = parse( expr[3] )
    return Minus( lhs, rhs )

  elseif op_symbol == :iszero
    arg = parse( expr[2] )
    return IsZero( arg )

  elseif op_symbol == :ifb
    condition = parse( expr[2] )
    true_branch = parse( expr[3] )
    false_branch = parse( expr[4] )
    return IfB( condition, true_branch, false_branch )

  elseif op_symbol == :with    # (with x (+ 5 1) (+ x x) )
    sym = expr[2]
    binding_expr = parse( expr[3] )
    body = parse( expr[4] )
    return With( sym, binding_expr, body )

  elseif op_symbol == :lambda
    formal = expr[2]
    formal_type = parse_type(expr[4])
    body = parse(expr[5])
    return FunDef( formal, formal_type, body )

  elseif op_symbol == :ncons
    f = parse(expr[2])
    r = parse(expr[3])
    return NCons( f, r )

  elseif op_symbol == :nisempty
    list = parse(expr[2])
    return NIsEmpty( list )

  elseif op_symbol == :nfirst
    list = parse(expr[2])
    return NFirst( list )

  elseif op_symbol == :nrest
    list = parse(expr[2])
    return NRest( list )

  else
    return FunApp( parse(expr[1]), parse(expr[2]) )

  end
end

function parse( expr::Any )
  throw( LispError("Invalid expression $expr") )
end

# ===================================================

# Parser for type expressions

function parse_type( t::Symbol )
  if (t == :number)
    return NumType()
  elseif (t == :boolean)
    return BoolType()
  elseif (t == :nlist)
    return NListType()
  end
end

function parse_type( t :: Array{Any} )
  return FunType( parse_type(t[1]),
                  parse_type(t[3]))
end

function parse_type( expr::Any )
  throw( LispError("Invalid type $expr") )
end

# ===================================================

# Type checking functions (modeled after the earlier calc)

function type_of_expr( ast::AE )
  return type_of_expr( ast, mtTypeEnvironment() )
end

function type_of_expr( ae::Num, env::TypeEnvironment )
  return NumType()
end

function type_of_expr( ae::Boolean, env::TypeEnvironment )
  return BoolType()
end

function type_of_expr( ae::Plus, env::TypeEnvironment )
  left = type_of_expr( ae.lhs, env )
  right = type_of_expr( ae.rhs, env )
  return type_of_math_expr( left, right )
end

function type_of_expr( ae::Minus, env::TypeEnvironment )
  left = type_of_expr( ae.lhs, env )
  right = type_of_expr( ae.rhs, env )
  return type_of_math_expr( left, right )
end

# Start of student code

function type_of_expr(ae::IsZero, env::TypeEnvironment)
  if same_type(type_of_expr(ae.arg, env), NumType())
    return BoolType()
  else
    throw(LispError("Conditional in $(ae) must evaluate to a number"))
  end
end

function type_of_expr(ae::IfB, env::TypeEnvironment)
  if same_type(type_of_expr(ae.condition, env), BoolType())
    true_type = type_of_expr(ae.true_branch, env)
    false_type = type_of_expr(ae.false_branch, env)
    if same_type(true_type, false_type)
      return true_type
    else
      throw(LispError("Branches in $(ae) must be of same type"))
    end
  else
    throw(LispError("Conditional in $(ae) must evaluate to boolean"))
  end
end

function type_of_expr(ae::With, env::TypeEnvironment)
  binding_type = type_of_expr(ae.binding_expr, env)
  return type_of_expr(ae.body, CTypeEnvironment(ae.name, binding_type, env))
end

function type_of_expr(ae::Id, env::TypeEnvironment)
  if typeof(env) != mtTypeEnvironment
    if ae.name == env.name
      return env.value
    else
      return type_of_expr(ae, env.parent)
    end
  else
    throw(LispError("$(ae) of type ID could not be resolved in environment"))
  end
end

function type_of_expr(ae::FunDef, env::TypeEnvironment)
  new_env = CTypeEnvironment(ae.formal_parameter, ae.formal_type, env)
  return FunType(ae.formal_type, type_of_expr(ae.fun_body, new_env))
end

function type_of_expr(ae::FunApp, env::TypeEnvironment)
  func_type = type_of_expr(ae.fun_expr, env)
  arg_type = type_of_expr(ae.arg_expr, env)
  if typeof(func_type) == FunType
    if same_type(func_type.arg_type, arg_type)
      return func_type.result_type
    else
      throw(LispError("Argument expression in $(ae) must match argument type"))
    end
  else
    throw(LispError("First parameter in $(ae) must be function definition"))
  end
end

function type_of_expr(ae::NEmpty, env::TypeEnvironment)
  return NListType()
end

function type_of_expr(ae::NCons, env::TypeEnvironment)
  f_type = type_of_expr(ae.f, env)
  r_type = type_of_expr(ae.r, env)
  if same_type(f_type, NumType()) && same_type(r_type, NListType())
    return NListType()
  else
    throw(LispError("Parameters of $(ae) must be of type number and nlist"))
  end
end

function type_of_expr(ae::NIsEmpty, env::TypeEnvironment)
  list_type = type_of_expr(ae.list, env)
  if same_type(list_type, NListType())
    return BoolType()
  else
    throw(LispError("Parameter of $(ae) must be of nlist"))
  end
end

function type_of_expr(ae::NFirst, env::TypeEnvironment)
  list_type = type_of_expr(ae.list, env)
  if same_type(list_type, NListType())
    return NumType()
  else
    throw(LispError("Parameter of $(ae) must be of type nlist"))
  end
end

function type_of_expr(ae::NRest, env::TypeEnvironment)
  list_type = type_of_expr(ae.list, env)
  if same_type(list_type, NListType())
    return list_type
  else
    throw(LispError("Parameter of $(ae) must be of type nlist"))
  end
end

# ===================================================

# Helper function for comparing two type values recursively if necessary

same_type( t1::FunType, t2::FunType ) =
    (same_type( t1.arg_type, t2.arg_type )
  && same_type( t1.result_type, t2.result_type ))

same_type{ T <: TypeVal }( t1::T, t2::T ) = true

same_type( t1::TypeVal, t2::TypeVal ) = false

# ===================================================

# Type judgments (could be folded into type checking functions)

function type_of_math_expr( left::NumType, right::NumType )
  return NumType()
end

function type_of_math_expr( left::Any, right::Any )
  throw( LispError("Operands for + or - must be numbers") )
end

# the rest of your type-judgement functions (if you choose to separate them) go here...

# ===================================================

# convenience function to make everything easier
function type_of_expr( expr::AbstractString )
  return type_of_expr( parse( Lexer.lex(expr) ) )
end

function test(line::AbstractString)
  new_line = chomp(line)
  if length(new_line) != 0
    println("Attempting to compute line \""*new_line*"\"")
    try
      type_out=type_of_expr(new_line)
      println("Process succeeded with type:")
      println(type_out)
    catch error_obj_a
      println("Found error in line")
      try
        println("Printing lexer output")
        lexer_out = Lexer.lex(new_line)
        println(lexer_out)
      catch error_obj_b
        println("Printing lexer error")
        #println(error_obj_b.msg)
      end
      try
        println("Printing parser output")
        parser_out = parse(lexer_out)
        println(parser_out)
      catch error_obj_c
        println("Printing parser error")
        #println(error_obj_c.msg)
      end
      println("Printing error object")
      println(error_obj_a.msg)
    end
  else
    println("Line is empty.")
  end
  println()
end

# ===================================================

end # module

# ===================================================

# Begin test cases
#=
println("Beginning test cases\n")
# These should evaluate
CITypes.test("1")
CITypes.test("false")
CITypes.test("(+ 1 5)")
CITypes.test("(- 1 5)")
CITypes.test("(iszero 5)")
CITypes.test("(ifb false 1 2)")
CITypes.test("(with x 5 x)")
CITypes.test("(lambda x : number x)")
CITypes.test("(lambda f : (number : number) (f 3))")
CITypes.test("(lambda f : (nlist : number) f nempty)")
CITypes.test("nempty")
CITypes.test("(ncons 3 nempty)")
CITypes.test("(nisempty nempty)")
CITypes.test("(nfirst nempty)")
CITypes.test("(nrest nempty)")
# These shouldn't evaluate
CITypes.test("x")
CITypes.test("(+ nempty 5)")
CITypes.test("(- nempty 5)")
CITypes.test("(iszero nempty)")
CITypes.test("(ifb 1 1 2)")
CITypes.test("(ifb true 1 nempty)")
CITypes.test("(ncons 3 3)")
CITypes.test("(nisempty 1)")
CITypes.test("(nfirst 1)")
CITypes.test("(nrest 1)")
println("Finished with analysis\n")
=#

max_stars_repo_name: scottcrossen/JuliaDevelopment

max_issues_repo_name: scottcrossen/JuliaDevelopment

max_forks_repo_name: scottcrossen/JuliaDevelopment


text:
using DelimitedFiles

# Run the CMD from the command function which will create Temp.wav file and extract the raw string information as a Tuple of vectors.

function feature(command)

    print("""Starting CMD "$command"\n""")
    # run command on CMD
    run(`cmd /k $command "&" exit`)

    # remove Temp.wav file
    rm(raw".\extraFiles\Temp.wav", force=true)

    # temporarily solution: use readdlm, remove commas, convert the strings to floats

    temp = DelimitedFiles.readdlm(raw".\extraFiles\Temp.features") # this contains strings, commas, and numbers

    for i in eachindex(temp)
           if(typeof(temp[i]) == SubString{String})
               temp[i] = parse(Float64,replace(temp[i], ","=>""))
           end
    end

    track = (t=temp[:,1],voicing=temp[:,2],f0=temp[:,3],f1=temp[:,4],f2=temp[:,5],f3=temp[:,6],f4=temp[:,7],f5=temp[:,8],f6=temp[:,9],f7=temp[:,10],f8=temp[:,11],f9=temp[:,12],f10=temp[:,13],f11=temp[:,14],f12=temp[:,15],f13=temp[:,16],f14=temp[:,17],f15=temp[:,18],Formants=temp[:,4:end])

end


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

using JSON, Base.Threads, MultiThreadingTools

export DLEnv
type DLEnv # immutable
  dir::AbstractString
  name::AbstractString
  config::Dict
  _gpus::Vector{Int}
  _ext_event_libs::Dict{String,DLData}
  _verbosity::Integer # 0: nothing, 1: only errors, 2: default, 3: all

  DLEnv() = DLEnv("config")
  DLEnv(name::AbstractString) = DLEnv(abspath(""), name)

  function DLEnv(dir::AbstractString, name::AbstractString)
    env = new(dir, name, Dict(), Int[], Dict(), 0)
    reload(env; replace=true)
    setup(env)
    return env
  end
end

type ConfigurationException <: Exception end

function Base.reload(env::DLEnv; replace=false)
  f = open(joinpath(env.dir,"$(env.name).json"), "r")
  dicttxt = readstring(f)  # file information to string
  env.config = JSON.parse(dicttxt)  # parse and transform data
  if replace
    env._verbosity = get(env.config, "verbosity", 2)
  end
end

export get_properties
function get_properties(env::DLEnv, name::AbstractString)
  return get(env.config, name, Dict())
end

export set_properties!
function set_properties!(env::DLEnv, name::AbstractString, props::Dict)
  env.config[name] = props
end

export new_properties!
function new_properties!(modifier, env::DLEnv, template_name::AbstractString, new_name::AbstractString)
  d = copy(get_properties(env, template_name))
  modifier(d)
  set_properties!(env, new_name, d)
  return d
end

export resolvepath
function resolvepath(env::DLEnv, path::AbstractString...)
  if isabspath(path[1])
    return joinpath(path...)
  else
    return joinpath(env.dir, path...)
  end
end

function Base.joinpath(env::DLEnv, elements::String...)
  return joinpath(env.dir, elements...)
end

function Base.info(env::DLEnv, level::Integer, msg::AbstractString)
  if env._verbosity >= level
    threadsafe_info(msg)
    flush(STDOUT)
  end
end

export set_verbosity!
function set_verbosity!(env::DLEnv, verbosity::Integer)
  env._verbosity = verbosity
end

export get_verbosity
function get_verbosity(env::DLEnv)
  return env._verbosity
end

export with_verbosity
function with_verbosity(action, env::DLEnv, verb)
    default_verbosity = env._verbosity
    env._verbosity = verb
    result = action()
    env._verbosity = default_verbosity
    return result
end

export verbose
verbose(action, env::DLEnv) = with_verbosity(action, env, 1000)

export silent
silent(action, env::DLEnv) = with_verbosity(action, env, 0)

function Base.getindex(env::DLEnv, key::String)
  return env.config[key]
end

function Base.setindex!(env::DLEnv, key::String, value)
  env.config[key] = value
end

function Base.haskey(env::DLEnv, key::String)
  return haskey(env.config, key)
end

export detectors
function detectors(env::DLEnv, dettype::AbstractString)
  if dettype == "BEGe"
    return BEGes_GERDA_II()
  elseif dettype == "coax" || dettype == "semi-coaxial"
    return Coax_GERDA_II()
  elseif dettype == "natural"
    return Natural_GERDA_II()
  elseif dettype == "used"
    return Used_GERDA_II()
  else
    throw(ArgumentError("Unknown detector type: $dettype"))
  end
end

function detectors(env::DLEnv, keywords::AbstractString...)
    sets = [detectors(env, keyword) for keyword in keywords]
    return intersect(sets...)
end

function detectors(env::DLEnv)
  return phase2_detectors
end

export _create_h5data
function _create_h5data(env::DLEnv, output_dir)
    info(env, 2, "Reading original data from $(env.config["path"])")
    isdir(output_dir) || mkdir(output_dir)
    if haskey(env.config, "rawformat")
        _seg_to_hdf5(env, output_dir)
    else
        _mgdo_to_hdf5(env, output_dir)
    end
    info(env, 3, "Converted raw data, HDF5 stored in $output_dir.")
end

function _mgdo_to_hdf5(env::DLEnv, output_dir)
  keylists = KeyList[]
  for (i,keylist_path) in enumerate(env.config["keylists"])
    if !endswith(keylist_path, ".txt")
      keylist_path = keylist_path*".txt"
    end
    push!(keylists, parse_keylist(resolvepath(env, keylist_path), keylist_path))
  end
    load_root(; verbosity=get_verbosity(env))
    mgdo_to_hdf5(env.config["path"], output_dir, keylists; verbosity=get_verbosity(env))
end

function _seg_to_hdf5(env::DLEnv, output_dir)
    src_dirs = env.config["path"]
    keylists = Vector{AbstractString}[]
    if env.config["rawformat"] == "SegHDF5"
      for src_dir in src_dirs
          content = readdir(src_dir)
          content = joinpath.(src_dir, content[find(f->endswith(f, ".hdf5"), content)])
          push!(keylists, content)
      end
      segh5_to_hdf5(env.config["rawformat"], keylists, output_dir, get_verbosity(env))
    else
      for src_dir in src_dirs
          content = readdir(src_dir)
          content = src_dir .* "/" .* content[find(f->endswith(f, ".root"), content)]
          push!(keylists, content)
      end
      load_root(; verbosity=get_verbosity(env))
      seg_to_hdf5(env.config["rawformat"], keylists, output_dir, get_verbosity(env))
    end
end

export getdata
function getdata(env::DLEnv; preprocessing::Union{AbstractString,Void}=nothing, targets::Array{String}=String[])
  if preprocessing==nothing
    return get(env, "raw"; targets=targets) do
      _get_raw_data(env; targets=targets)
    end
  else
    data = _get_raw_data(env; targets=[preprocessing])
    preprocessed = get(env, preprocessing; targets=targets) do
      preprocess(env, data, preprocessing)
    end
    if preprocessed == nothing
      return nothing
    end
    # Else check whether cache is up to date
    steps = get_properties(env, preprocessing)["preprocessing"]
    cache_up_to_date = true
    for lib in preprocessed
      cached_steps = lib[:preprocessing]
      if cached_steps != steps
        cache_up_to_date = false
      end
    end
    if cache_up_to_date
      return preprocessed
    else
      info(env, 2, "Refreshing cache of 'preprocessed'.")
      delete!(env, preprocessing)
      return getdata(env; preprocessing=preprocessing, targets=targets)
    end
  end
end

function _get_raw_data(env::DLEnv; targets::Array{String}=String[])
  raw_dir = resolvepath(env, "data", "raw")
  if !isdir(raw_dir)
    _create_h5data(env, raw_dir)
  end
  return lazy_read_all(raw_dir)
end

export preprocess
function preprocess(env::DLEnv, data::DLData, config_name)
  config = env.config[config_name]
  @assert isa(config, Dict)

  select_channels=parse_detectors(config["detectors"])
  info(env, 3, "Selected channels: $select_channels")
  if isa(select_channels,Vector) && length(select_channels) > 0
    filter!(data, :detector_name, select_channels)
  end

  _builtin_filter(env, config, "test-pulses", data, :isTP, isTP -> isTP == 0)
  _builtin_filter(env, config, "baseline-events", data, :isBL, isBL -> isBL == 0)
  _builtin_filter(env, config, "unphysical-events", data, :E, E -> (E > 0) && (E < 9999))
  _builtin_filter(env, config, "low-energy-events", data, :E, E -> (E > 0) && (E < 1000))

  N_dset = length(parse_datasets(env, config["sets"]))
  result = DLData(fill(EventLibrary(zeros(0,0)), length(data)*N_dset))
  result.dir = _cachedir(env, "tmp-preprocessed-$config_name")

  steps = convert(Array{String}, config["preprocessing"])

  for i in 1:length(data)
    lib = data.entries[i]
    info(env,2, "Preprocessing $(lib[:name])")
    lib_t = preprocess_transform(env, lib, steps; copyf=identity)
    _builtin_filter(env, config, "failed-preprocessing", lib_t, :FailedPreprocessing, fail -> fail == 0)
    part_data = DLData(collect(values(split(env, lib_t, config["sets"]))))
    write_all_sequentially(part_data, result.dir, true)
    info(env,3, "Wrote datasets of $(lib[:name]) and released allocated memory.")
    dispose(lib)
    dispose(lib_t)
    for j in 1:length(part_data)
      result.entries[N_dset*(i-1) + j] = part_data.entries[j]
    end
    @assert length(lib.waveforms) == 0 && length(lib_t.waveforms) == 0
  end

  return result
end

function _builtin_filter(env::DLEnv, config::Dict, ftype_key::String, data::EventCollection, label, exclude_prededicate)
  # TODO per detector to save memory
  ftype = config[ftype_key]
  info(env, 3, "Filter $ftype_key is set to $ftype.")
  if ftype == "exclude"
    result = filter!(data, label, exclude_prededicate)
    return result
  elseif ftype == "include"
    return data
  elseif ftype == "only"
    return filter!(data, label, x -> !exclude_prededicate(x))
  else
    throw(ArgumentError("Unknown filter keyword in configuration $ftype_key: $ftype"))
  end
end

function parse_datasets(env::DLEnv, strdict::Dict)
  result = Dict{AbstractString,Vector{AbstractFloat}}()
  if haskey(env.config, "keylists")
    requiredlength = length(env["keylists"])
  else
    requiredlength = 1
  end
  for (key,value) in strdict
    if isa(value, Vector)
      @assert length(value) == requiredlength
      result[key] = value
    elseif isa(value, Real)
      result[key] = fill(value, requiredlength)
    else
      throw(ConfigurationException())
    end
  end
  return result
end

Base.split(env::DLEnv, data::EventCollection, strdict::Dict) = split(data, parse_datasets(env, strdict))

function setup(env::DLEnv)
  e_mkdir(env.dir)
  e_mkdir(joinpath(env.dir, "data"))
  e_mkdir(joinpath(env.dir, "models"))
  e_mkdir(joinpath(env.dir, "plots"))
end
e_mkdir(dir) = !isdir(dir) && mkdir(dir)

function Base.get(compute, env::DLEnv, lib_name::String; targets::Array{String}=String[], uninitialize=true)
  if !isempty(targets) && containsall(env, targets)
    info(env, 2, "Skipping retrieval of '$lib_name'.")
    return nothing
  end
  if contains(env, lib_name)
    info(env, 2, "Retrieving '$lib_name' from cache.")
    return get(env, lib_name)
  else
    info(env, 2, "Computing '$lib_name'...")
    data = compute()
    info(env, 3, "Computation of '$lib_name' finished.")

    # check type
    if !isa(data, DLData)
      throw(TypeError(Symbol(compute), "get_or_compute", DLData, data))
    end

    if data.dir == nothing
      info(env, 3, "Writing computed data to env as '$lib_name' (uninitialize=$uninitialize).")
      push!(env, lib_name, data; uninitialize=uninitialize)
      return data
    else
      # Rename directory
      targetdir = _cachedir(env, lib_name)
          if targetdir != data.dir
              info(env, 3, "Moving data from $(data.dir) to $targetdir and reloading.")
              mv(data.dir, targetdir)
              return get(env, lib_name)
            else
                info(env, 3, "Data is already in the right location ($(data.dir)).")
                return data
            end
    end
  end
end

function Base.get(env::DLEnv, lib_name::String)
  _ensure_ext_loaded(env, lib_name)
  return env._ext_event_libs[lib_name]
end

Base.push!(env::DLEnv, lib_name::String, lib::EventLibrary; uninitialize::Bool=false) = push!(env, lib_name, DLData([lib]); uninitialize=uninitialize)

function Base.push!(env::DLEnv, lib_name::String, data::DLData; uninitialize::Bool=false)
  _ensure_ext_loaded(env, lib_name)
  env._ext_event_libs[lib_name] = data
  if env.config["cache"]
    write_all_multithreaded(data, joinpath(env.dir, "data", lib_name), uninitialize)
  end
end

function Base.contains(env::DLEnv, lib_name::String)
  if haskey(env._ext_event_libs, lib_name)
    return true
  end
  return env.config["cache"] && isdir(_cachedir(env, lib_name))
end

export containsall
function containsall(env::DLEnv, lib_names::Array{String})
  for lib_name in lib_names
    if !contains(env, lib_name) return false end
  end
  return true
end

function _ensure_ext_loaded(env::DLEnv, lib_name::String)
  if !haskey(env._ext_event_libs, lib_name)
    c_dir = _cachedir(env, lib_name)
    if isdir(c_dir)
      env._ext_event_libs[lib_name] = lazy_read_all(c_dir)
    end
  end
end

function Base.delete!(env::DLEnv, lib_name::String)
  if haskey(env._ext_event_libs, lib_name)
    delete!(env._ext_event_libs, lib_name)
  end
  # Delete cached files (even if cache is set to false)
  c_dir = _cachedir(env, lib_name)
  if isdir(c_dir)
    rm(c_dir; recursive=true)
    info(env, 2, "Deleted cached $lib_name")
  end
end

function _cachedir(env::DLEnv, lib_name::String)
  return joinpath(env.dir, "data", lib_name)
end

function network(env::DLEnv, name::String)
    dir = joinpath(env.dir, "models", name)
    isdir(dir) || mkdir(dir)
    return NetworkInfo(name, dir, get_properties(env, name), to_context(env._gpus))
end

export use_gpus
function use_gpus(env, gpus::Int...)
  env._gpus = [gpus...]
end


text:
import Statistics

"""
"""
function simple_moving_average(
        x::AbstractVector,
        window::Integer,
        ) where T
    if window < 0
        error("window must be >=0")
    end
    n = length(x)
    result = similar(x)
    for i = 1:n
        lower_bound = max(1, i - window)
        upper_bound = min(n, i + window)
        result[i] = Statistics.mean(x[lower_bound:upper_bound])
    end
    return result
end


text:
immutable PoolingLayer <: Mocha.Layer
    bottoms::Vector{Symbol}
    name::AbstractString
    pathkeys::Vector{PathKey}
    pooling::Function
    tops::Vector{Symbol}
end

function PoolingLayer( ;
        bottoms::Vector{Symbol} = Symbol[],
        name::AbstractString = "pooling",
        pathkeys::Vector{PathKey} = PathKey[],
        pooling::Function = mean,
        tops::Vector{Symbol} = Symbol[])
    @assert length(bottoms) > 0
    @assert length(tops) == length(bottoms)
    PoolingLayer(bottoms, name, pathkeys, pooling, tops)
end

Mocha.can_do_bp(::PoolingLayer) = false
Mocha.has_neuron(::PoolingLayer) = false
Mocha.has_param(::PoolingLayer) = false


text:
import DoubleFloats: Double64

issquare(A::AbstractMatrix) = isequal(size(A)...)

"""
    row_scaling!(A, D::AbstractVector, tol=maximum(size(A))^2 * eps(real(eltype(A))))

Divide each row of `A` by its ∞-norm. This minimizes the ∞ condition number
for all possible scalings. In order to avoid scaling near zero rows back to
1 we only scale a row if its ∞-norm is larger than `tol`.
"""
function row_scaling!(A::AbstractMatrix{T}, D::AbstractVector, tol=max(size(A)...)^2 * eps(real(T))) where {T}
    D .= 0.0
    d_max = convert(real(T), -Inf)
    @inbounds for j=1:size(A, 2), i=1:size(A, 1)
        D[i] = @fastmath max(D[i], abs(A[i,j]))
        d_max = @fastmath max(d_max, D[i])
    end

    # In the case that the row norms are all less than 1 we have a lower bound
    # the maximal row norm
    # To avoid that we scale a zero matrix, we need a tol for wich we don't scale any more
    bound = clamp(d_max, tol, 1.0)
    D .= max.(D, bound)

    @inbounds for j=1:size(A, 2), i=1:size(A, 1)
        A[i, j] /= D[i]
    end
    A
end

@enum ActiveFactorization begin
    QR_FACTORIZATION
    LU_FACTORIZATION
end

"""
    Jacobian{T, <:Factorization}

Data structure holding a Jacobian matrix `J` with it's factorization.
Additional contains a vector `D` of scalings which row equilibarate `J`.

    Jacobian(J)

Setup a Jacobian with Jacobian matrices of the form `J`.
"""
mutable struct Jacobian{T}
    J::Matrix{T} # Jacobian matrix
    D::Vector{Float64} # Inverse of scaling factors D
    # lu is only there for square systems
    lu::Union{Nothing, LinearAlgebra.LU{T,Matrix{T}}} # LU Factorization of D * J
    qr::LinearAlgebra.QRPivoted{T,Matrix{T}} # Pivoted QR Factorization of D * J
    active_factorization::ActiveFactorization
    b::Vector{T} # Vector for overdetermined systems
    r::Vector{T} # Vector for iterative refinement
    perm::Vector{Int} # work copy of the permutation in the qr
    qr_work::Vector{ComplexF64} # preallocated memory for pivoted qr
    qr_rwork::Vector{Float64} # preallocated memory for pivoted qr
    ormqr_work::Vector{ComplexF64}
    # Numerical Informations
    corank::Int # The numerical corank
    corank_proposal::Int
    cond::Float64 # Estimate of the Jacobian
    # The relative number of digits lost during the solution of the linear systems
    # in Newton's method. See `solve_with_digits_lost!` in utilities/linear_algebra.jl
    # for how this is computed.
    digits_lost::Union{Nothing,Float64}
end

function Jacobian(A::AbstractMatrix, corank::Int=0)
    m, n = size(A)
    lu = m ≠ n ? nothing : LinearAlgebra.lu!(Random.randn!(complex.(float.(A))))

    work = Vector{ComplexF64}(undef, 1)
    ormqr_work = Vector{ComplexF64}(undef, 1)
    rwork = Vector{Float64}(undef, 2n)
    jpvt = zeros(BlasInt, n)
    B = ComplexF64.(A)
    tau = similar(B, min(m, n))
    geqp3!(B, jpvt, tau, work, rwork)
    qr = LinearAlgebra.QRPivoted(B, tau, jpvt)
    active_factorization = lu === nothing ? QR_FACTORIZATION : LU_FACTORIZATION
    J = copy(B)
    D = ones(m)
    r = similar(J, size(J, 1))
    b = copy(r)
    perm = zeros(Int, length(qr.p))
    cond = 1.0
    digits_lost = nothing
    Jacobian(J, D, lu, qr, active_factorization, b, r, perm, work, rwork, ormqr_work,
        corank, 0, cond, digits_lost)
end

function update_rank!(Jac::Jacobian)
    Jac.corank = Jac.corank_proposal
    Jac
end

function reset!(jacobian::Jacobian)
    jacobian.cond = 1.0
    jacobian.corank = 0
    jacobian.corank_proposal = 0
    jacobian.digits_lost = nothing
    jacobian
end

"""
    update!(jacobian, J, corank::Int=jacobian.corank)

Update the `jacobian` struct with the new Matrix `J`.
"""
function update!(Jac::Jacobian, J::AbstractMatrix; update_infos::Bool=false)
    copyto!(Jac.J, J)
    updated_jacobian!(Jac; update_infos=update_infos)
    Jac
end

"""
    update_jacobian!(jacobian)

This computes a factorization of the stored Jacobian matrix.
Call this instead of `update!` if `jacobian.J` got updated.
"""
function updated_jacobian!(Jac::Jacobian{ComplexF64}; update_infos::Bool=false) where {T}
    if !update_infos && issquare(Jac.J) && Jac.corank == 0
        Jac.lu !== nothing || return Jac
        copyto!(Jac.lu.factors, Jac.J)
        lu!(Jac.lu.factors, nothing, Jac.lu.ipiv)
        Jac.active_factorization = LU_FACTORIZATION
    else
        copyto!(Jac.qr.factors, Jac.J)

        m, n = size(Jac.J)
        if m == n
            # only apply row scaling for square matrices since
            # otherwise we change the problem
            row_scaling!(Jac.qr.factors, Jac.D, 1)
        end
        # this computes a pivoted qr factorization, i.e.,
        # qr!(Jac.qr.factors, Val(true)) but without allocating new memory
        geqp3!(Jac.qr.factors, Jac.qr.jpvt, Jac.qr.τ, Jac.qr_work, Jac.qr_rwork)

        ε = max(n,m) * eps() * 10
        # check rank 0
        rnm = min(n,m)
        r₁ = abs(real(Jac.qr.factors[1,1]))
        corank = 0
        if r₁ < ε
            corank = rnm
            Jac.cond = inv(ε)
        else
            for i in 2:rnm
                rᵢ = abs(real(Jac.qr.factors[i,i]))
                if ε * r₁ > rᵢ
                    corank = rnm - i + 1
                    break
                end
            end
            Jac.cond = max(1.0, r₁) / abs(real(Jac.qr.factors[rnm, rnm]))
        end
        Jac.corank_proposal = corank

        Jac.active_factorization = QR_FACTORIZATION
    end
    Jac
end

"""
    solve!([x, ], Jac::Jacobian, b)

Solve `Jac.J * x = b` inplace. Assumes `Jac` contains already the correct factorization.
This stores the result in `x`.
"""
function solve!(x, Jac::Jacobian, b::AbstractVector; update_digits_lost::Bool=false, refinement_step::Bool=false)
    if Jac.active_factorization == LU_FACTORIZATION
        lu = Jac.lu
        if lu !== nothing
            lu_ldiv!(x, lu, b)
        end
    elseif Jac.active_factorization == QR_FACTORIZATION
        if issquare(Jac.J)
            # we applied row scaling for square matrices
            @inbounds for i in eachindex(Jac.b)
                Jac.b[i] = b[i] / Jac.D[i]
            end
            custom_ldiv!(x, Jac.qr, Jac.b, Jac.corank, Jac.perm, Jac.ormqr_work )
        else
            Jac.b .= b
            custom_ldiv!(x, Jac.qr, Jac.b, Jac.corank, Jac.perm, Jac.ormqr_work )
        end
    end
    if issquare(Jac.J) && Jac.corank == 0 && update_digits_lost && !refinement_step
        norm_x = euclidean_norm(x)
        norm_δx = iterative_refinement_step!(x, Jac, b)
        Jac.digits_lost = log₁₀(norm_δx / eps(norm_x))
    elseif !issquare(Jac.J)
        Jac.digits_lost = nothing
    end
    x
end
solve!(A::Jacobian, b) = solve!(b, A::Jacobian, b)

const LA = LinearAlgebra

# This is an adoption of the Julia implementation of ldiv! for QRPivoted
function custom_ldiv!(x, A::LinearAlgebra.QRPivoted{ComplexF64},
                      b::Vector{ComplexF64}, corank::Int,
                      perm::Vector{Int}, ormqr_work::Vector{ComplexF64})
    mA, nA = size(A.factors)
    nr = min(mA,nA)
    nrhs = length(b)
    if nr == 0
        return x
    end
    if corank == nr
        fill!(x, zero(eltype(x)))
        return x
    end
    rank = nr - corank

    if rank == nA == nr
        # The following is equivalent to LA.lmul!(LA.adjoint(A.Q), b) but
        # we can also pass the preallocated work vector
        ormqr!('L', 'C', A.factors, A.τ, b, ormqr_work)
        ldiv_upper!(A.factors, b; singular_exception=false)
        @inbounds for i in 1:nr
            x[i] = b[i]
            perm[i] = A.p[i]
        end
        Base.invpermute!!(x, perm)
    else
        C, τ = LA.LAPACK.tzrzf!(A.factors[1:rank,:])
        LA.ldiv!(LA.UpperTriangular(C[1:rank,1:rank]),view(LA.lmul!(LA.adjoint(A.Q), b), 1:rank))
        for i in 1:rank
            x[i] = b[i]
        end
        for i in rank+1:nA
            x[i] = zero(eltype(x))
        end
        LA.LAPACK.ormrz!('L', eltype(x)<:Complex ? 'C' : 'T', C, τ, reshape(x, length(x), 1))
        invpermute!(x, A.p)
    end
    return x
end

"""
    iterative_refinement_step!(x, A, b, fac, T, r)

Apply one step of iterative refinement where the residual is computed with precision `T`.
Returns the euclidean norm of the update.
"""
function iterative_refinement_step!(x, Jac::Jacobian, b, ::Type{T}=eltype(x)) where T
    residual!(Jac.r, Jac.J, x, b, T)
    δx = solve!(Jac.r, Jac, Jac.r; refinement_step=true)
    norm_δx = euclidean_norm(δx)
    @inbounds for i in eachindex(x)
        x[i] = convert(T, x[i]) - convert(T, δx[i])
    end

    return norm_δx
end


text:
using Test, GeneticAlgorithm

STRING_TO_FIND = "hellohellohello"

function createGeneExample()
    return rand(STRING_TO_FIND)
end

function fitnessExample(ind)
    allDifferent = filter(((a,b),) -> a != b, collect(zip(STRING_TO_FIND, ind)))
    return length(allDifferent)
end

@test join(runGA(fitnessExample, createGeneExample, length(STRING_TO_FIND), maxNumberOfIterations=10, logging=false)[1]) == "hello"


text:
module Shared_types

  typealias Params Dict{Any,Any}
  export Params

end

max_stars_repo_name: jadnohra/mlcp-lib

max_issues_repo_name: jadnohra/mlcp-lib

max_forks_repo_name: jadnohra/mlcp-lib