text:
function GAIL_D_loss(gan_loss)
    (D, 𝒫, 𝒟_ex, 𝒟_π; info = Dict()) ->begin
        Lᴰ(gan_loss, D, 𝒟_ex[:a], 𝒟_π[:a], yD = (𝒟_ex[:s],), yG = (𝒟_π[:s],))
    end
end

function OnPolicyGAIL(;π,
                       S,
                       𝒟_demo,
                       normalize_demo::Bool=true,
                       D::ContinuousNetwork,
                       solver=PPO,
                       gan_loss::GANLoss=GAN_BCELoss(),
                       d_opt::NamedTuple=(;),
                       log::NamedTuple=(;),
                       kwargs...)

    d_opt = TrainingParams(;loss = GAIL_D_loss(gan_loss), name="discriminator_", d_opt...)
    normalize_demo && (𝒟_demo = normalize!(deepcopy(𝒟_demo), S, action_space(π)))
    𝒟_demo = 𝒟_demo |> device(π)

    function GAIL_callback(𝒟; info=Dict())
        batch_train!(D, d_opt, (;), 𝒟_demo, 𝒟, info=info)

        discriminator_signal = haskey(𝒟, :advantage) ? :advantage : :return
        D_out = value(D, 𝒟[:a], 𝒟[:s]) # This is swapped because a->x and s->y and the convention for GANs is D(x,y)
        r = Base.log.(sigmoid.(D_out) .+ 1f-5) .- Base.log.(1f0 .- sigmoid.(D_out) .+ 1f-5)
        𝒟[discriminator_signal] .= r
    end
    𝒮 = solver(;π=π, S=S, post_batch_callback=GAIL_callback, log=(dir="log/onpolicygail", period=500, log...), kwargs...)
    𝒮.c_opt = nothing # disable the critic
    𝒮
end


text:
export simple_hessian, nhessian, jacobian, local_nhessian

function (_::Type{Inv{GVar}})(x::GVar{<:GVar,<:GVar})
    Partial{:x}(x)
end

@i function ⊕(*)(out!::Partial, x::Partial, y::Partial)
    ⊕(*)(value(out!), value(x), value(y))
end

@i function ⊕(identity)(out!::Partial, x::Partial)
    ⊕(identity)(value(out!), value(x))
end

@i function hessian1(f, args; kwargs, index::Int)
    @safe @assert count(x -> x isa Loss, args) == 1
    @anc iloss = 0
    @routine getiloss begin
        for i=1:length(args)
            if (tget(args,i) isa Loss, iloss==i)
                iloss += identity(i)
            end
        end
    end

    # forward
    f'(args...; kwargs...)

    (~Loss)(tget(args,iloss))
    for i = 1:length(args)
        GVar(grad(tget(args,i)))
        GVar(value(tget(args,i)))
    end
    grad(grad(tget(args,index))) ⊕ 1.0
    # backward#2
    (Loss)(tget(args,iloss))
    (~f')(args...; kwargs...)

    ~@routine getiloss
end

function simple_hessian(f, args::Tuple; kwargs=())
    N = length(args)
    hmat = zeros(N, N)
    for i=1:N
        res = hessian1(f, args; kwargs=kwargs, index=i)
        hmat[:,i] .= map(x->x isa Loss ? grad(value(value(x))) :  grad(value(x)), res[2])
    end
    hmat
end

function nhessian(f, args; kwargs=(), η=1e-5)
    largs = Any[args...]
    narg = length(largs)
    res = zeros(narg, narg)
    for i = 1:narg
        @instr value(largs[i]) ⊕ η/2
        gpos = gradient(f, (largs...,); kwargs=kwargs)
        @instr value(largs[i]) ⊖ η
        gneg = gradient(f, (largs...,); kwargs=kwargs)
        @instr value(largs[i]) ⊕ η/2
        res[:,i] .= (gpos .- gneg)./η
    end
    return res
end

function local_nhessian(f, args; kwargs=())
    nargs = length(args)
    hes = zeros(nargs,nargs,nargs)
    for j=1:nargs
        @instr Loss(tget(args, j))
        hes[:,:,j] .= nhessian(f, args; kwargs=kwargs)
        @instr (~Loss)(tget(args, j))
    end
    hes
end

"""J(in, out)"""
function jacobian(f, args; kwargs=())
    narg = length(args)
    res = zeros(narg, narg)
    for i = 1:narg
        @instr Loss(tget(args, i))
        res[:,i] .= gradient(f, args; kwargs=kwargs)
        @instr (~Loss)(tget(args, i))
    end
    return res
end

ext: jl


text:
function rot_elem(N::Int, v::Vector{Int})
    elem = 0.0

    for i in 1:N
        l = v[2i-1]
        elem += l*(l+1)
    end

    elem
end

function sparse_rot(block::Block{N}) where {N}
    rows = Int[]
    elems = Float64[]

    for row in 1:block.size
        v = block.vectors[row]
        elem = rot_elem(N, v)

        push!(rows, row)
        push!(elems, elem)
    end

    sparse(rows, rows, elems, block.size, block.size)
end

"""
    sparse_rot(basis::SingleBlockBasis)

Compute a sparse representation of the dimensionless rotational energy matrix.
"""
sparse_rot(basis::SingleBlockBasis) = sparse_rot(basis.block)

"""
    diag_rot(basis::SingleBlockBasis)

Compute a diagonal representation of the dimensionless rotational energy
matrix.
"""
diag_rot(basis::SingleBlockBasis) = diag(sparse_rot(basis.block))

"""
    diag_rot(basis::MultiBlockBasis)

Compute a diagonal representation of the dimensionless rotational energy
matrix.
"""
function diag_rot(basis::MultiBlockBasis)
    result = Dict{BlockLabel,Vector{Float64}}()

    for (label, block) in basis.blocks
        result[label] = diag(sparse_rot(block))
    end

    result
end

function C(l::Int, lp::Int, m::Int)
    if lp == l-1
        sqrt((l-m)*(l+m)/((2l-1)*(2l+1)))
    elseif lp == l+1
        sqrt((lp-m)*(lp+m)/((2lp-1)*(2lp+1)))
    end
end

function Sp(l::Int, lp::Int, m::Int)
    if lp == l-1
        -sqrt((l-m-1)*(l-m)/((2l-1)*(2l+1)))
    elseif lp == l+1
        sqrt((lp+m)*(lp+m+1)/((2lp-1)*(2lp+1)))
    end
end

function Sm(l::Int, lp::Int, m::Int)
    if lp == l-1
        sqrt((l+m-1)*(l+m)/((2l-1)*(2l+1)))
    elseif lp == l+1
        -sqrt((lp-m)*(lp-m+1)/((2lp-1)*(2lp+1)))
    end
end

function pot_elem(l1::Int, m1::Int, l2::Int, m2::Int, l1p::Int, m1p::Int, l2p::Int, m2p::Int)
    if m1p == m1-1 && m2p == m2+1
        Sm(l1, l1p, m1) * Sp(l2, l2p, m2)
    elseif m1p == m1+1 && m2p == m2-1
        Sp(l1, l1p, m1) * Sm(l2, l2p, m2)
    elseif m1p == m1 && m2p == m2
        -4.0 * C(l1, l1p, m1) * C(l2, l2p, m2)
    end
end

function pot_rows(f::Function, N::Int, l_max::Int, l_total_max::Int, v::Vector{Int}, w::Vector{Int})
    l_total = 0
    for i in 1:N
        l_total += v[2i-1]
    end

    # All rotor pairs.
    for i in 1:(N-1)
        for j in (i+1):N
            # All potentially allowed transitions.
            for delta_li in [-1, 1]
                for delta_lj in [-1, 1]
                    for delta_m in -1:1
                        w .= v
                        w[2i-1] += delta_li
                        w[2i] += delta_m
                        w[2j-1] += delta_lj
                        w[2j] -= delta_m

                        # Bounds on l.
                        0 <= w[2i-1] <= l_max || continue
                        0 <= w[2j-1] <= l_max || continue
                        l_total + delta_li + delta_lj <= l_total_max || continue
                        # Bounds on m.
                        -w[2i-1] <= w[2i] <= w[2i-1] || continue
                        -w[2j-1] <= w[2j] <= w[2j-1] || continue

                        unscaled = pot_elem(v[2i-1], v[2i], v[2j-1], v[2j], w[2i-1], w[2i], w[2j-1], w[2j])
                        elem = 0.5 * unscaled / (j-i)^3

                        f(elem)
                    end
                end
            end
        end
    end

    nothing
end

function sparse_pot(basis::AbstractBasis{N}, block::Block{N}) where {N}
    rows = Int[]
    cols = Int[]
    elems = Float64[]

    w = Array{Int}(undef, 2N)

    for col in 1:block.size
        v = block.vectors[col]

        pot_rows(N, basis.l_max, basis.l_total_max, v, w) do elem
            row = block.lookup[w]

            push!(rows, row)
            push!(cols, col)
            push!(elems, elem)
        end
    end

    sparse(rows, cols, elems, block.size, block.size)
end

"""
    sparse_pot(basis::SingleBlockBasis)

Compute a sparse representation of the dimensionless potential energy matrix
for dipole-dipole interactions between all rotor pairs.
"""
sparse_pot(basis::SingleBlockBasis) = sparse_pot(basis, basis.block)

"""
    dense_pot(basis::SingleBlockBasis)

Compute a dense representation of the dimensionless potential energy matrix for
dipole-dipole interactions between all rotor pairs.
"""
dense_pot(basis::SingleBlockBasis) = Matrix(sparse_pot(basis, basis.block))

"""
    dense_pot(basis::MultiBlockBasis)

Compute a dense representation of the dimensionless potential energy matrix for
dipole-dipole interactions between all rotor pairs.
"""
function dense_pot(basis::MultiBlockBasis)
    result = Dict{BlockLabel,Matrix{Float64}}()

    for (label, block) in basis.blocks
        result[label] = Matrix(sparse_pot(basis, block))
    end

    result
end

"""
    scale_dense_pot(pot::Matrix{Float64}, k::Float64)

Scale the matrix elements of `pot` in place by a factor of `k`.
"""
function scale_dense_pot(pot::Matrix{Float64}, k::Float64)
    pot .*= k
end

"""
    scale_dense_pot(pot::Dict{BlockLabel,Matrix{Float64}}, k::Float64)

Scale the matrix elements of `pot` in place by a factor of `k`.
"""
function scale_dense_pot(pot::Dict{BlockLabel,Matrix{Float64}}, k::Float64)
    for block in values(pot)
        block .*= k
    end
end

"""
    dense_x(basis::SingleBlockBasis)

Compute a dense representation of the single-rotor x operator.
"""
function dense_x(basis::SingleBlockBasis)
    l_max = basis.l_max

    result = zeros(Float64, (l_max+1)^2, (l_max+1)^2)

    for l in 0:l_max
        for m in (-l):l
            col = l*(l+1) + m + 1

            for lp in 0:l_max
                for mp in (-lp):lp
                    row = lp*(lp+1) + mp + 1

                    if lp == l+1 && mp == m+1
                        result[row, col] -= 0.5 * sqrt((l+m+1)*(l+m+2)/((2l+1)*(2l+3)))
                    elseif lp == l-1 && mp == m+1
                        result[row, col] += 0.5 * sqrt((l-m-1)*(l-m)/((2l-1)*(2l+1)))
                    elseif lp == l+1 && mp == m-1
                        result[row, col] += 0.5 * sqrt((l-m+1)*(l-m+2)/((2l+1)*(2l+3)))
                    elseif lp == l-1 && mp == m-1
                        result[row, col] -= 0.5 * sqrt((l+m-1)*(l+m)/((2l-1)*(2l+1)))
                    end
                end
            end
        end
    end

    result
end

"""
    dense_iy(basis::SingleBlockBasis)

Compute a dense representation of the single-rotor iy operator.
"""
function dense_iy(basis::SingleBlockBasis)
    l_max = basis.l_max

    result = zeros(Float64, (l_max+1)^2, (l_max+1)^2)

    for l in 0:l_max
        for m in (-l):l
            col = l*(l+1) + m + 1

            for lp in 0:l_max
                for mp in (-lp):lp
                    row = lp*(lp+1) + mp + 1

                    if lp == l+1 && mp == m+1
                        result[row, col] -= 0.5 * sqrt((l+m+1)*(l+m+2)/((2l+1)*(2l+3)))
                    elseif lp == l-1 && mp == m+1
                        result[row, col] += 0.5 * sqrt((l-m-1)*(l-m)/((2l-1)*(2l+1)))
                    elseif lp == l+1 && mp == m-1
                        result[row, col] -= 0.5 * sqrt((l-m+1)*(l-m+2)/((2l+1)*(2l+3)))
                    elseif lp == l-1 && mp == m-1
                        result[row, col] += 0.5 * sqrt((l+m-1)*(l+m)/((2l-1)*(2l+1)))
                    end
                end
            end
        end
    end

    result
end

"""
    dense_z(basis::SingleBlockBasis)

Compute a dense representation of the single-rotor z operator.
"""
function dense_z(basis::SingleBlockBasis)
    l_max = basis.l_max

    result = zeros(Float64, (l_max+1)^2, (l_max+1)^2)

    for l in 0:l_max
        for m in (-l):l
            col = l*(l+1) + m + 1

            for lp in 0:l_max
                for mp in (-lp):lp
                    row = lp*(lp+1) + mp + 1

                    if lp == l+1 && mp == m
                        result[row, col] += sqrt((l-m+1)*(l+m+1)/((2l+1)*(2l+3)))
                    elseif lp == l-1 && mp == m
                        result[row, col] += sqrt((l-m)*(l+m)/((2l-1)*(2l+1)))
                    end
                end
            end
        end
    end

    result
end


text:
"""
```
imgb = binarize(Yen(), img)
binarize!(Yen(), img)
```

Uses Yen's maximum correlation criterion to threshold the histogram of a given
image.

# Arguments

The function arguments are described in more detail below.

##  `img`

An `AbstractGray` image which is to be binarized into background (`Gray(0)`)
and foreground (`Gray(1)`) pixels.

## Reference

Yen JC, Chang FJ, Chang S (1995), "A New Criterion for Automatic Multilevel Thresholding", IEEE Trans. on Image Processing 4 (3): 370-378, [doi:10.1109/83.366472](https://doi.org/10.1109/83.366472)
"""
function binarize(algorithm::Yen,  img::AbstractArray{T,2}) where T <: AbstractGray
  binarize!(algorithm, copy(img))
end

"""
```
t = binarize!(Yen(), img)
```
Same as [`binarize`](@ref binarize(::Yen, ::AbstractArray{<:AbstractGray,2})) except that it modifies the image that was passed as an argument.
"""
function binarize!(algorithm::Yen,  img::AbstractArray{T,2}) where T <: AbstractGray
  edges, counts = build_histogram(img,  256)
  t = find_threshold(HistogramThresholding.Yen(), counts[1:end], edges)
  map!(img,img) do val
    val < t ? Gray(0) : Gray(1)
  end
end


text:
# This is the pixel-by-pixel MNIST problem from: Le, Q. V., Jaitly,
# N., & Hinton, G. E. (2015). A Simple Way to Initialize Recurrent
# Networks of Rectified Linear Units. arXiv preprint arXiv:1504.00941.

module MNISTPixels
using Main, Knet, ArgParse
using Knet: nextidx, stack_isempty, mat2d

isdefined(:MNIST) || include("mnist.jl")

function main(args=ARGS)
    isa(args, AbstractString) && (args=split(args))
    opts = parse_commandline(args)
    info("Pixel-by-pixel MNIST problem from Le et al. 2015.")
    println(opts)
    for (k,v) in opts; @eval ($(symbol(k))=$v); end
    seed > 0 && setseed(seed)

    global trn = Pixels(mat2d(MNIST.xtrn), MNIST.ytrn; batch=batchsize, epoch=epochsize, bootstrap=true)
    global tst = Pixels(mat2d(MNIST.xtst), MNIST.ytst; batch=batchsize)
    global net = compile(:mpixels; rnn=symbol(nettype), hidden=hidden, nclass=10, winit=Gaussian(0,winit), fbias=fbias)
    setp(net; lr=lrate)

    l = [0f0,0f0]; m = [0f0,0f0]; acc = 0
    for epoch=1:epochs
        train(net, trn, softloss; gclip=gclip, losscnt=fill!(l,0), maxnorm=fill!(m,0))
        println(tuple(:trn,epoch*trn.epochsize,l[1]/l[2],m...))
        gcheck > 0 && gradcheck(net, f->gradloss(f,trn,softloss;grad=true), f->gradloss(f,trn,softloss); gcheck=gcheck)
        if epoch % testfreq == 0
            acc = 1-test(net, tst, zeroone)
            println(tuple(:tst,epoch*trn.epochsize,acc))
        end
        flush(STDOUT)
    end
    return (acc, l[1]/l[2], m...)
end

@knet function mpixels(x; rnn=nothing, hidden=0, nclass=0, o...)
    h = rnn(x; o..., out=hidden)
    if predict
        return wbf(h; o..., out=nclass, f=:soft)
    end
end

# @knet function p1(x; rnn=nothing, hidden=0, o...)
#     y = rnn(x; o..., out=hidden)
# end

# @knet function p2(x; nclass=0, o...)
#     y = wbf(x; o..., out=nclass, f=soft)
# end

function train(f, data, loss; gclip=0, losscnt=nothing, maxnorm=nothing)
    reset!(f)
    for (x,ygold) in data
        if ygold == nothing
            sforw(f, x; predict=false)
        else
            ypred = sforw(f, x; predict=true)
            losscnt[1] += loss(ypred, ygold); losscnt[2] += 1
            sback(f, ygold, loss)
            while !stack_isempty(f); sback(f); end
            g = gnorm(f); g > maxnorm[2] && (maxnorm[2]=g)
            gscale = (g > gclip > 0 ? gclip/g : 1)
            update!(f; gscale=gscale)
            w = wnorm(f); w > maxnorm[1] && (maxnorm[1]=w) # TODO: take this maxnorm calculation elsewhere, have a statistics optional function? just be smart about scaling?  batch-normalization, msra etc.
            reset!(f)
        end
    end
end

function test(f, data, loss; gclip=0, losscnt=nothing, maxnorm=nothing)
    sumloss = numloss = 0
    reset!(f)
    for (x,ygold) in data
        if ygold == nothing
            forw(f, x; predict=false)
        else
            ypred = forw(f, x; predict=true)
            sumloss += loss(ypred, ygold); numloss += 1
            reset!(f)
        end
    end
    sumloss / numloss
end

function gradloss(f, data, loss; grad=false, seed=42)
    data_rng = data.rng
    data.rng = MersenneTwister(seed)
    reset!(f)
    myforw = grad ? sforw : forw
    loss1 = 0
    for (x,ygold) in data
        if ygold == nothing
            myforw(f, x; predict=false)
        else
            ypred = myforw(f, x; predict=true)
            loss1 = loss(ypred, ygold)
            if grad
                sback(f, ygold, loss)
                while !stack_isempty(f); sback(f); end
            end
            break
        end
    end
    data.rng = data_rng
    return loss1
end

# Data generator:

import Base: start, next, done

# input comes in as xtrn(784,60000), ytrn(10,60000)
# the batch arg determines how many images we present in parallel
# each call to the next should output batch pixels, i.e. x(1,batch), y=nothing
# the last pixel should be served as x(1,batch), y(10,batch)

type Pixels; x; y; rng; datasize; epochsize; batchsize; bootstrap; shuffle; xbatch; ybatch; images;
    function Pixels(x, y; rng=MersenneTwister(0), epoch=ccount(x), batch=16, bootstrap=false, shuffle=false)
        nx = ccount(x)
        nx == ccount(y) || error("Item count mismatch")
        shuf = (shuffle ? shuffle!(rng,[1:nx;]) : nothing)
        xbatch = similar(x, (1,batch))
        ybatch = similar(y, (clength(y),batch))
        new(x, y, rng, nx, epoch, batch, bootstrap, shuf, xbatch, ybatch, nothing)
    end
end

# state is an image/pixel pair
start(d::Pixels)=(d.shuffle != nothing && shuffle!(d.rng, d.shuffle); (0,0))

# epochsize is given as image count, not pixel count
done(d::Pixels, s)=(s[1] >= d.epochsize)

function next(d::Pixels, s)
    (n,t) = s                   # image and pixel count
    t==0 && (d.images = nextidx(d,n)) # image indices for current batch
    nb = length(d.images)       # batch size, xbatch[1:nb] needs to be filled
    nt = clength(d.x)           # pixels per image
    t += 1                      # next pixel to serve
    for b=1:nb                  # batch index
        i=d.images[b]           # image index
        d.xbatch[b] = d.x[t,i]
        t == nt && (d.ybatch[:,b] = d.y[:,i])
    end
    return t < nt ?
    ((d.xbatch, nothing), (n, t)) :
    ((d.xbatch, d.ybatch), (n+nb, 0))
end

function parse_commandline(args)
    s = ArgParseSettings()
    @add_arg_table s begin
        "--epochs"
        help = "Number of epochs to train"
        arg_type = Int
        default = 1
        "--epochsize"
        help = "number of training examples per epoch"
        arg_type = Int
        default = 100 # 10000
        "--testfreq"
        help = "Compute test accuracy every acc epochs"
        arg_type = Int
        default = 1 # 10
        "--batchsize"
        help = "minibatch size"
        arg_type = Int
        default = 20 # 16
        "--hidden"
        help = "number of hidden units"
        arg_type = Int
        default = 100
        "--lrate"
        help = "learning rate"
        arg_type = Float64
        default = 0.005          # paper says 1e-8?
        "--gclip"
        help = "gradient clip"
        arg_type = Float64
        default = 1.0
        "--gcheck"              # TODO: gives terrible results, check.
        help = "gradient check"
        arg_type = Int
        default = 0
        "--nettype"
        help = "type of network"
        default = "irnn" # "lstm"
        "--fbias"
        help = "forget gate bias for lstm"
        arg_type = Float64
        default =  1.0
        "--winit"
        help = "stdev for weight initialization"
        arg_type = Float64
        default =  0.01
        "--rtol"
        help = "rtol parameter for gradient checks"
        arg_type = Float64
        default =  0.01
        "--atol"
        help = "atol parameter for gradient checks"
        arg_type = Float64
        default =  0.01
        "--seed"
        help = "Random seed"
        arg_type = Int
        default = 1003
    end
    parse_args(args,s)
end

!isinteractive() && !isdefined(Core.Main,:load_only) && main(ARGS)

end # module

# NOTES:

# batchsize=16 is too slow.  I got good progress with the following:
# Dict{AbstractString,Any}("hidden"=>100,"rtol"=>0.01,"batchsize"=>200,"lrate"=>0.001,"train"=>10000,"gclip"=>1.0,"acc"=>10,"gcheck"=>0,"fbias"=>1.0,"epochs"=>100,"atol"=>0.01,"seed"=>1003,"type"=>"irnn")

# DONE: reconsider the data interface.
# we should make the model interface not insist on Data, any iterable should do, in fact remove the Data type.
# x field should probably be named something else like data.
# if we are ever going to need zero or more than one input, then we need to decide how to represent:
# (x1,x2,...,y) or (y,x1,x2,...): probably the first
# (y,) or y for no input tuples: probably the first

# Sample output for debugging:
# julia> include("mnistpixels.jl")
# Dict{AbstractString,Any}("hidden"=>100,"rtol"=>0.01,"batchsize"=>200,"lrate"=>0.005,"train"=>10000,"gclip"=>1.0,"acc"=>10,"gcheck"=>0,"fbias"=>1.0,"epochs"=>100,"atol"=>0.01,"seed"=>1003,"type"=>"irnn","std"=>0.01)
# (:trn,10000,2.261733271443444,10.417845f0,537.297f0)
# (:trn,20000,2.3187815377575203,10.421104f0,5327.3f0)
# (:trn,30000,2.509523073669509,10.423339f0,10088.415f0)
# (:trn,40000,2.412966892292024,10.425221f0,9402.696f0)
# (:trn,50000,2.3695122707481397,10.4281435f0,13412.184f0)
# (:trn,60000,2.268378955795277,10.4292965f0,5621.7524f0)
# (:trn,70000,2.1479991531190734,10.431348f0,2633.2693f0)
# (:trn,80000,2.1210844235681883,10.433363f0,1978.7836f0)
# (:trn,90000,2.153266387169849,10.435466f0,3289.0789f0)
# (:trn,100000,2.1011612425016724,10.437336f0,3303.134f0)
# (:tst,100000,0.1735)

# why isn't this giving the same results?
# Could be softloss vs xentloss or Float32 vs Float64.

# Here is the new sample output:
# julia> include("mnistpixels.jl")
# INFO: Loading MNIST...
#   5.863069 seconds (279.54 k allocations: 498.444 MB, 1.61% gc time)
# Dict{AbstractString,Any}("hidden"=>100,"rtol"=>0.01,"batchsize"=>200,"lrate"=>0.005,"train"=>10000,"gclip"=>1.0,"acc"=>10,"std"=>0.01,"gcheck"=>0,"fbias"=>1.0,"epochs"=>100,"atol"=>0.01,"seed"=>1003,"type"=>"irnn")
# (:trn,10000,2.262306f0,10.417742f0,537.452f0)
# (:trn,20000,2.277638f0,10.4206705f0,2240.0566f0)
# (:trn,30000,2.45132f0,10.423019f0,11971.223f0)
# (:trn,40000,2.2885273f0,10.424828f0,7045.522f0)
# (:trn,50000,2.3827329f0,10.427834f0,12723.293f0)
# (:trn,60000,2.2246392f0,10.429581f0,3678.2678f0)
# (:trn,70000,2.157727f0,10.431644f0,2867.3608f0)
# (:trn,80000,2.1695452f0,10.434076f0,4772.0884f0)
# (:trn,90000,2.2076025f0,10.435675f0,13315.624f0)
# (:trn,100000,2.174602f0,10.437442f0,7863.392f0)
# (:tst,100000,0.2939)

# sample output for lstm:
#
# [dy_052@hpc3013 examples]$ julia mnistpixels.jl --type lstm
# WARNING: requiring "Options" in module "Main" did not define a corresponding module.
# INFO: Loading MNIST...
#   6.725308 seconds (363.19 k allocations: 358.143 MB, 0.62% gc time)
# Dict{AbstractString,Any}("hidden"=>100,"rtol"=>0.01,"batchsize"=>200,"lrate"=>0.005,"train"=>10000,"gclip"=>1.0,"acc"=>10,"std"=>0.01,"gcheck"=>0,"fbias"=>1.0,"epochs"=>100,"atol"=>0.01,"seed"=>1003,"type"=>"lstm")
# (:trn,10000,2.3025522f0,14.711951f0,0.1208294f0)
# (:trn,20000,2.3024626f0,14.717022f0,0.12132876f0)
# (:trn,30000,2.3024027f0,14.721309f0,0.10962964f0)
# (:trn,40000,2.3023393f0,14.725827f0,0.14956589f0)

# S2C no longer accepts Net, it expects kfun:
    # p1 = (nettype == "irnn" ? Net(irnn; out=hidden, winit=Gaussian(0,winit)) :
    #       nettype == "lstm" ? Net(lstm; out=hidden, fbias=fbias) :
    #       error("Unknown network type "*nettype))
    # p2 = Net(wbf; out=10, winit=Gaussian(0,winit), f=soft)


text:
using JuMP, EAGO, Gurobi

m = Model(with_optimizer(EAGO.Optimizer, relaxed_optimizer = Gurobi.Optimizer(OutputFlag=0)))

# ----- Variables ----- #
x_Idx = Any[1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22]
@variable(m, x[x_Idx])
JuMP.set_lower_bound(x[21], 0.0)
JuMP.set_lower_bound(x[1], -0.5)
JuMP.set_upper_bound(x[1], 0.5)
JuMP.set_lower_bound(x[2], 5.4)
JuMP.set_upper_bound(x[2], 6.4)
JuMP.set_lower_bound(x[3], 0.4)
JuMP.set_upper_bound(x[3], 1.4)
JuMP.set_lower_bound(x[4], 4.9)
JuMP.set_upper_bound(x[4], 5.9)
JuMP.set_lower_bound(x[5], 1.3)
JuMP.set_upper_bound(x[5], 2.3)
JuMP.set_lower_bound(x[6], 3.9)
JuMP.set_upper_bound(x[6], 4.9)
JuMP.set_lower_bound(x[7], 2.1)
JuMP.set_upper_bound(x[7], 3.1)
JuMP.set_lower_bound(x[8], 4.1)
JuMP.set_upper_bound(x[8], 5.1)
JuMP.set_lower_bound(x[9], 2.8)
JuMP.set_upper_bound(x[9], 3.8)
JuMP.set_lower_bound(x[10], 3.0)
JuMP.set_upper_bound(x[10], 4.0)
JuMP.set_lower_bound(x[11], 3.9)
JuMP.set_upper_bound(x[11], 4.9)
JuMP.set_lower_bound(x[12], 3.2)
JuMP.set_upper_bound(x[12], 4.2)
JuMP.set_lower_bound(x[13], 4.7)
JuMP.set_upper_bound(x[13], 5.7)
JuMP.set_lower_bound(x[14], 2.3)
JuMP.set_upper_bound(x[14], 3.3)
JuMP.set_lower_bound(x[15], 5.6)
JuMP.set_upper_bound(x[15], 6.6)
JuMP.set_lower_bound(x[16], 2.3)
JuMP.set_upper_bound(x[16], 3.3)
JuMP.set_lower_bound(x[17], 6.0)
JuMP.set_upper_bound(x[17], 7.0)
JuMP.set_lower_bound(x[18], 1.9)
JuMP.set_upper_bound(x[18], 2.9)
JuMP.set_lower_bound(x[19], 6.9)
JuMP.set_upper_bound(x[19], 7.9)
JuMP.set_lower_bound(x[20], 1.0)
JuMP.set_upper_bound(x[20], 2.0)
JuMP.set_upper_bound(x[21], 10.0)
JuMP.set_lower_bound(x[22], -2.0)
JuMP.set_upper_bound(x[22], 2.0)

# ----- Constraints ----- #
@NLconstraint(m, e2, x[22]*x[1]-x[2]+x[21] == 0.0)
@NLconstraint(m, e3, x[22]*x[3]-x[4]+x[21] == 0.0)
@NLconstraint(m, e4, x[22]*x[5]-x[6]+x[21] == 0.0)
@NLconstraint(m, e5, x[22]*x[7]-x[8]+x[21] == 0.0)
@NLconstraint(m, e6, x[22]*x[9]-x[10]+x[21] == 0.0)
@NLconstraint(m, e7, x[22]*x[11]-x[12]+x[21] == 0.0)
@NLconstraint(m, e8, x[22]*x[13]-x[14]+x[21] == 0.0)
@NLconstraint(m, e9, x[22]*x[15]-x[16]+x[21] == 0.0)
@NLconstraint(m, e10, x[22]*x[17]-x[18]+x[21] == 0.0)
@NLconstraint(m, e11, x[22]*x[19]-x[20]+x[21] == 0.0)

# ----- Objective ----- #
@NLobjective(m, Min, ( (x[1])^2+ (x[2]-5.9)^2+ (x[3]-0.9)^2+ (x[4]-5.4)^2+ (x[5]-1.8)^2+ (x[6]-4.4)^2+ (x[7]-2.6)^2+ (x[8]-4.6)^2+ (x[9]-3.3)^2+ (x[10]-3.5)^2+ (x[11]-4.4)^2+ (x[12]-3.7)^2+ (x[13]-5.2)^2+ (x[14]-2.8)^2+ (x[15]-6.1)^2+ (x[16]-2.8)^2+ (x[17]-6.5)^2+ (x[18]-2.4)^2+ (x[19]-7.4)^2+ (x[20]-1.5)^2))

JuMP.optimize!(m)

backend_opt = backend(m).optimizer.model.optimizer
last_relaxed_opt = backend_opt.relaxed_optimizer
run_time = backend_opt._run_time
println("run time: $run_time")


text:
using Base.Test
using SparseVectors

import SparseVectors: exact_equal

### BLAS Level-1

let x = sprand(16, 0.5), x2 = sprand(16, 0.4)
    xf = full(x)
    xf2 = full(x2)

    # axpy!
    for c in [1.0, -1.0, 2.0, -2.0]
        y = full(x)
        @test is(axpy!(c, x2, y), y)
        @test y == full(x2 * c + x)
    end

    # scale
    let sx = SparseVector(x.n, x.nzind, x.nzval * 2.5)
        @test exact_equal(scale(x, 2.5), sx)
        @test exact_equal(scale(2.5, x), sx)
        @test exact_equal(x * 2.5, sx)
        @test exact_equal(2.5 * x, sx)
        @test exact_equal(x .* 2.5, sx)
        @test exact_equal(2.5 .* x, sx)

        xc = copy(x)
        @test is(scale!(xc, 2.5), xc)
        @test exact_equal(xc, sx)
    end

    # dot
    let dv = dot(xf, xf2)
        @test dot(x, x) == sumabs2(x)
        @test dot(x2, x2) == sumabs2(x2)
        @test_approx_eq dot(x, x2) dv
        @test_approx_eq dot(x2, x) dv
        @test_approx_eq dot(full(x), x2) dv
        @test_approx_eq dot(x, full(x2)) dv
    end
end

let x = complex(sprand(32, 0.6), sprand(32, 0.6)),
    y = complex(sprand(32, 0.6), sprand(32, 0.6))
    xf = full(x)::Vector{Complex128}
    yf = full(y)::Vector{Complex128}
    @test_approx_eq dot(x, x) dot(xf, xf)
    @test_approx_eq dot(x, y) dot(xf, yf)
end

### BLAS Level-2:

## dense A * sparse x -> dense y

let A = randn(9, 16), x = sprand(16, 0.7)
    xf = full(x)
    for α in [0.0, 1.0, 2.0], β in [0.0, 0.5, 1.0]
        y = rand(9)
        rr = α * A * xf + β * y
        @test is(A_mul_B!(α, A, x, β, y), y)
        @test_approx_eq y rr
    end
    y = A * x
    @test isa(y, Vector{Float64})
    @test_approx_eq A * x A * xf
end

let A = randn(16, 9), x = sprand(16, 0.7)
    xf = full(x)
    for α in [0.0, 1.0, 2.0], β in [0.0, 0.5, 1.0]
        y = rand(9)
        rr = α * A'xf + β * y
        @test is(At_mul_B!(α, A, x, β, y), y)
        @test_approx_eq y rr
    end
    y = At_mul_B(A, x)
    @test isa(y, Vector{Float64})
    @test_approx_eq y At_mul_B(A, xf)
end

## sparse A * sparse x -> dense y

let A = sprandn(9, 16, 0.5), x = sprand(16, 0.7)
    Af = full(A)
    xf = full(x)
    for α in [0.0, 1.0, 2.0], β in [0.0, 0.5, 1.0]
        y = rand(9)
        rr = α * Af * xf + β * y
        @test is(A_mul_B!(α, A, x, β, y), y)
        @test_approx_eq y rr
    end
    y = densemv(A, x)
    @test isa(y, Vector{Float64})
    @test_approx_eq y Af * xf
end

let A = sprandn(16, 9, 0.5), x = sprand(16, 0.7)
    Af = full(A)
    xf = full(x)
    for α in [0.0, 1.0, 2.0], β in [0.0, 0.5, 1.0]
        y = rand(9)
        rr = α * Af'xf + β * y
        @test is(At_mul_B!(α, A, x, β, y), y)
        @test_approx_eq y rr
    end
    y = densemv(A, x; trans='T')
    @test isa(y, Vector{Float64})
    @test_approx_eq y At_mul_B(Af, xf)
end

let A = complex(sprandn(7, 8, 0.5), sprandn(7, 8, 0.5)),
    x = complex(sprandn(8, 0.6), sprandn(8, 0.6)),
    x2 = complex(sprandn(7, 0.75), sprandn(7, 0.75))
    Af = full(A)
    xf = full(x)
    x2f = full(x2)
    @test_approx_eq densemv(A, x; trans='N') Af * xf
    @test_approx_eq densemv(A, x2; trans='T') Af.' * x2f
    @test_approx_eq densemv(A, x2; trans='C') Af'x2f
end

## sparse A * sparse x -> sparse y

let A = sprandn(9, 16, 0.5), x = sprand(16, 0.7), x2 = sprand(9, 0.7)
    Af = full(A)
    xf = full(x)
    x2f = full(x2)

    y = A * x
    @test isa(y, SparseVector{Float64,Int})
    @test all(nonzeros(y) .!= 0.0)
    @test_approx_eq full(y) Af * xf

    y = At_mul_B(A, x2)
    @test isa(y, SparseVector{Float64,Int})
    @test all(nonzeros(y) .!= 0.0)
    @test_approx_eq full(y) Af'x2f
end

let A = complex(sprandn(7, 8, 0.5), sprandn(7, 8, 0.5)),
    x = complex(sprandn(8, 0.6), sprandn(8, 0.6)),
    x2 = complex(sprandn(7, 0.75), sprandn(7, 0.75))
    Af = full(A)
    xf = full(x)
    x2f = full(x2)

    y = A * x
    @test isa(y, SparseVector{Complex128,Int})
    @test_approx_eq full(y) Af * xf

    y = At_mul_B(A, x2)
    @test isa(y, SparseVector{Complex128,Int})
    @test_approx_eq full(y) Af.' * x2f

    y = Ac_mul_B(A, x2)
    @test isa(y, SparseVector{Complex128,Int})
    @test_approx_eq full(y) Af'x2f
end


text:
const hdfs_fsstore = Dict{@compat(Tuple{AbstractString,Integer,AbstractString}), Vector{Any}}()
typealias IPv4v6 Union(IPv4,IPv6)

function finalize_hdfs_fs(fs::HdfsFS)
    (C_NULL == fs.ptr) && return
    key, arr = _get_ptr_ref(fs.host, fs.port, fs.user, false)
    if((arr[1] -= 1) == 0)
        ccall((:hdfsDisconnect, _libhdfs), Int32, (Ptr{Void},), fs.ptr)
        arr[2] = fs.ptr = C_NULL
        # TODO: remove from hdfs_fsstore to save space
    end
end

# hdfs_connect returns pointers to the same handle across multiple calls for the same file system
# so we must have an abstraction with reference count
function _get_ptr_ref(host::AbstractString, port::Integer, user::AbstractString="", incr::Bool=true)
    key = (host, port, user)
    if(haskey(hdfs_fsstore, key))
        arr = hdfs_fsstore[key]
        if(arr[1] > 0)
            incr && (arr[1] += 1)
            return key, arr
        end
    end
    key, [0, C_NULL]
end

function hdfs_connect(host::AbstractString="default", port::Integer=0, user::AbstractString="")
    key, arr = _get_ptr_ref(host, port, user)
    (0 != arr[1]) && return HdfsFS(host, @compat(Int(port)), user, arr[2])
    porti32 = @compat Int32(port)
    ptr = (user == "") ?
            ccall((:hdfsConnect, _libhdfs), Ptr{Void}, (Ptr{UInt8}, Int32), bytestring(host), porti32) :
            ccall((:hdfsConnectAsUser, _libhdfs), Ptr{Void}, (Ptr{UInt8}, Int32, Ptr{UInt8}), bytestring(host), porti32, bytestring(user))
    (C_NULL == ptr) && error("hdfs connect failed")
    hdfs_fsstore[key] = [1, ptr]
    HdfsFS(host, @compat(Int(port)), "", ptr)
end

function hdfs_connect(url::HdfsURL)
    uri = URI(url.url)
    if ':' in uri.userinfo
        uname,_passwd = userinfo(uri.userinfo)
    else
        uname = uri.userinfo
        _passwd = ""
    end
    hname = uri.host
    portnum = uri.port
    hdfs_connect(hname, portnum, (nothing == uname)?"":uname)
end

# file control flags need to be done manually
function hdfs_open(fs::HdfsFS, path::AbstractString, mode::AbstractString, buffer_sz::Integer=0, replication::Integer=0, bsz::Integer=0)
    flags = -1
    (mode == "r") && (flags = Base.JL_O_RDONLY)
    (mode == "w") && (flags = Base.JL_O_WRONLY)
    (mode == "a") && (flags = (Base.JL_O_WRONLY | Base.JL_O_APPEND))
    (flags == -1) && error("unknown open mode $(mode)")

    file = ccall((:hdfsOpenFile, _libhdfs), Ptr{Void}, (Ptr{Void}, Ptr{UInt8}, Int32, Int32, Int16, Int32), fs.ptr, bytestring(path),
        convert(Int32, flags),
        convert(Int32, buffer_sz),
        convert(Int16, replication),
        convert(Int32, bsz))
    (C_NULL == file) && error("error opening file $(path)")
    return HdfsFile(fs, path, file)
end

function hdfs_open(f::HdfsFile, mode::AbstractString, buffer_sz::Integer=0, replication::Integer=0, bsz::Integer=0)
    (C_NULL != f.ptr) && (0 != hdfs_close(f)) && error("error closing file")
    fnew = hdfs_open(f.fs, f.path, mode, buffer_sz, replication, bsz)
    f.ptr = fnew.ptr
    fnew.ptr = C_NULL
    f
end

function hdfs_open(url::AbstractString, mode::AbstractString, buffer_sz::Integer=0, replication::Integer=0, bsz::Integer=0)
    uri = defrag(URI(url))
    (uri.schema != "hdfs") && error("not a HDFS URL")
    uname,_passwd = userinfo(uri.userinfo)
    hname = uri.host
    portnum = uri.port
    fname = uri.path

    fs = hdfs_connect(hname, portnum, (nothing == uname)?"":uname)
    hdfs_open(fs, fname, mode, buffer_sz, replication, bsz)
end
hdfs_open(url::HdfsURL, mode::AbstractString, buffer_sz::Integer=0, replication::Integer=0, bsz::Integer=0) = hdfs_open(url.url, mode, buffer_sz, replication, bsz)

function hdfs_close(file::HdfsFile)
    ret = (file.fs.ptr == C_NULL) ? 0 : ccall((:hdfsCloseFile, _libhdfs), Int32, (Ptr{Void}, Ptr{Void}), file.fs.ptr, file.ptr)
    (0 == ret) && (file.ptr = C_NULL)
    ret
end

function hdfs_exists(fs::HdfsFS, path::AbstractString)
    ret = ccall((:hdfsExists, _libhdfs), Int32, (Ptr{Void}, Ptr{UInt8}), fs.ptr, bytestring(path))
    # TODO: examine why hdfsExists returns -1 when file does not exist. as per docs it should return 1
    #(ret < 0) && println("error checking file $(path) exists. error: $(ret)")
    (ret == 0)
end

hdfs_seek(file::HdfsFile, desired_pos) = ccall((:hdfsSeek, _libhdfs), Int32, (Ptr{Void}, Ptr{Void}, Int64), file.fs.ptr, file.ptr, convert(Int64,desired_pos))

hdfs_tell(file::HdfsFile) = ccall((:hdfsTell, _libhdfs), Int64, (Ptr{Void}, Ptr{Void}), file.fs.ptr, file.ptr)

hdfs_read(file::HdfsFile, buff::Ptr, len::Integer) = ccall((:hdfsRead, _libhdfs), Int32, (Ptr{Void}, Ptr{Void}, Ptr{Void}, Int32), file.fs.ptr, file.ptr, buff, convert(Int32,len))
function hdfs_read(file::HdfsFile, len::Integer)
    buff = Array(UInt8, len)
    (-1 == (r = hdfs_read(file, convert(Ptr{Void}, pointer(buff)), len))) && error("error reading file: -1")
    (buff, r)
end

hdfs_pread(file::HdfsFile, position::Int64, buff::Ptr, len::Integer) = ccall((:hdfsPread, _libhdfs), Int32, (Ptr{Void}, Ptr{Void}, Int64, Ptr{Void}, Int32), file.fs.ptr, file.ptr, position, buff, len)
function hdfs_pread(file::HdfsFile, position::Int64, len::Integer)
    buff = Array(UInt8, len)
    (-1 == (r = hdfs_pread(file, position, buff, len))) && error("error reading file: -1")
    (buff, r)
end

#can be passed an ASCIIString (length not necessary in that case)
hdfs_write(file::HdfsFile, buff::Ptr, len::Integer) = ccall((:hdfsWrite, _libhdfs), Int32, (Ptr{Void}, Ptr{Void}, Ptr{Void}, Int32), file.fs.ptr, file.ptr, buff, convert(Int32,len))
hdfs_write(file::HdfsFile, buff::ASCIIString, len::Integer) = hdfs_write(file, convert(Ptr{Void}, pointer(buff.data)), len)
hdfs_write(file::HdfsFile, buff::ASCIIString) = hdfs_write(file, buff, length(buff))

hdfs_flush(file::HdfsFile) = ccall((:hdfsFlush, _libhdfs), Int32, (Ptr{Void}, Ptr{Void}), file.fs.ptr, file.ptr)

hdfs_available(file::HdfsFile) = ccall((:hdfsAvailable, _libhdfs), Int32, (Ptr{Void},Ptr{Void}), file.fs.ptr, file.ptr)

hdfs_copy(srcFS::HdfsFS, src::AbstractString, dstFS::HdfsFS, dst::AbstractString) = ccall((:hdfsCopy, _libhdfs), Int32, (Ptr{Void}, Ptr{UInt8}, Ptr{Void}, Ptr{UInt8}), srcFS.ptr, bytestring(src), dstFS.ptr, bytestring(dst))

hdfs_move(srcFS::HdfsFS, src::AbstractString, dstFS::HdfsFS, dst::AbstractString) = ccall((:hdfsMove, _libhdfs), Int32, (Ptr{Void}, Ptr{UInt8}, Ptr{Void}, Ptr{UInt8}), srcFS.ptr, bytestring(src), dstFS.ptr, bytestring(dst))

hdfs_delete(fs::HdfsFS, path::AbstractString) = ccall((:hdfsDelete, _libhdfs), Int32, (Ptr{Void}, Ptr{UInt8}), fs.ptr, bytestring(path))

hdfs_rename(fs::HdfsFS, old_path::AbstractString, new_path::AbstractString) = ccall((:hdfsRename, _libhdfs), Int32, (Ptr{Void}, Ptr{UInt8}, Ptr{UInt8}), fs.ptr, bytestring(old_path), bytestring(new_path))

function hdfs_pwd(fs::HdfsFS, buff::Vector{UInt8})
    ptr = ccall((:hdfsGetWorkingDirectory, _libhdfs), Ptr{UInt8}, (Ptr{Void}, Ptr{UInt8}, Int32), fs.ptr, buff, convert(Int32,length(buff)))
    (C_NULL == ptr) && error("error getting working directory")
    len = ccall((:strlen, "libc"), UInt, (Ptr{UInt8},), buff)
    ASCIIString(buff[1:len])
end
hdfs_pwd(fs::HdfsFS, buff_sz::Integer) = hdfs_pwd(fs, Array(UInt8, buff_sz))
hdfs_pwd(fs::HdfsFS) = hdfs_pwd(fs, 1024)

hdfs_cd(fs::HdfsFS, path::AbstractString) = ccall((:hdfsSetWorkingDirectory, _libhdfs), Int32, (Ptr{Void}, Ptr{UInt8}), fs.ptr, bytestring(path))

hdfs_mkdir(fs::HdfsFS, path::AbstractString) = ccall((:hdfsCreateDirectory, _libhdfs), Int32, (Ptr{Void}, Ptr{UInt8}), fs.ptr, bytestring(path))

hdfs_set_replication(fs::HdfsFS, path::AbstractString, replication::Integer)=ccall((:hdfsSetReplication, _libhdfs), Int32, (Ptr{Void}, Ptr{UInt8}, Int16), fs.ptr, bytestring(path), @compat(Int16(replication)))

function hdfs_list_directory(fs::HdfsFS, path::AbstractString)
    num_entries = Int32[-1]
    info_ptr = ccall((:hdfsListDirectory, _libhdfs), Ptr{c_hdfsfileinfo}, (Ptr{Void}, Ptr{UInt8}, Ptr{Int32}), fs.ptr, bytestring(path), num_entries)
    (0 == num_entries[1]) && (return [])
    (C_NULL == info_ptr) && error(string("Error listing path ", path))

    n = @compat Int(num_entries[1])
    ret = [HdfsFileInfo(x) for x in pointer_to_array(info_ptr, (n,))]
    ccall((:hdfsFreeFileInfo, _libhdfs), Void, (Ptr{Void}, Int32), info_ptr, num_entries[1])
    ret
end

function hdfs_get_path_info(fs::HdfsFS, path::AbstractString)
    info_ptr = ccall((:hdfsGetPathInfo, _libhdfs), Ptr{c_hdfsfileinfo}, (Ptr{Void}, Ptr{UInt8}), fs.ptr, bytestring(path))
    (C_NULL == info_ptr) && error(string("Error getting path ", path))

    ret = HdfsFileInfo(info_ptr)
    ccall((:hdfsFreeFileInfo, _libhdfs), Void, (Ptr{Void}, Int32), info_ptr, 1)
    ret
end

hdfs_is_directory(fs::HdfsFS, path::AbstractString) = (hdfs_get_path_info(fs,path).kind == HDFS_OBJ_DIR)

function hdfs_blocks(fs::HdfsFS, path::AbstractString, start::Integer=1, len::Integer=0, as_ip::Bool=false)
    (len == 0) && (len = hdfs_get_path_info(fs, path).size)
    c_ptr = ccall((:hdfsGetHosts, _libhdfs), Ptr{Ptr{Ptr{UInt8}}}, (Ptr{Void}, Ptr{UInt8}, Int64, Int64), fs.ptr, bytestring(path),
                    convert(Int64,start),
                    convert(Int64,len))
    (C_NULL == c_ptr) && error("Error getting hosts for file $(path)")

    i = 1
    ret_vals = Array(as_ip ? Array{IPv4v6,1} : Array{ASCIIString,1}, 0)
    while true
        h_list = unsafe_load(c_ptr, i)
        (h_list == C_NULL) && break

        j = 1
        arr = Array(as_ip ? IPv4v6 : ASCIIString, 0)
        while true
            hname = unsafe_load(h_list, j)
            (hname == C_NULL) && break
            shname = bytestring(hname)
            push!(arr, as_ip ? getaddrinfo(shname) : shname)
            #print(bytestring(hname), ", ")
            j += 1
        end
        push!(ret_vals, arr)
        #println("")
        i += 1
    end

    ccall((:hdfsFreeHosts, _libhdfs), Void, (Ptr{Ptr{Ptr{UInt8}}},), c_ptr)
    ret_vals
end
hdfs_blocks(f::HdfsFile, start::Integer=1, len::Integer=0, as_ip::Bool=false) = hdfs_blocks(f.fs, f.path, start, len, as_ip)
hdfs_blocks(url::HdfsURL, start::Integer=1, len::Integer=0, as_ip::Bool=false) = hdfs_blocks(hdfs_connect(url), URI(url.url).path, start, len, as_ip)

hdfs_get_default_block_size(fs::HdfsFS) = ccall((:hdfsGetDefaultBlockSize, _libhdfs), Int64, (Ptr{Void},), fs.ptr)

hdfs_get_capacity(fs::HdfsFS) = ccall((:hdfsGetCapacity, _libhdfs), Int64, (Ptr{Void},), fs.ptr)

hdfs_get_used(fs::HdfsFS) = ccall((:hdfsGetUsed, _libhdfs), Int64, (Ptr{Void},), fs.ptr)

hdfs_chown(fs::HdfsFS, path::AbstractString, owner::AbstractString, group::AbstractString) = ccall((:hdfsChown, _libhdfs), Int32, (Ptr{Void}, Ptr{UInt8}, Ptr{UInt8}, Ptr{UInt8}), fs.ptr, bytestring(path), bytestring(owner), bytestring(group))

hdfs_chmod(fs::HdfsFS, path::AbstractString, mode::Int16) = ccall((:hdfsChmod, _libhdfs), Int32, (Ptr{Void}, Ptr{UInt8}, Int16), fs.ptr, bytestring(path), mode)

hdfs_utime(fs::HdfsFS, path::AbstractString, mtime::Integer, atime::Integer) = ccall((:hdfsUtime, _libhdfs), Int32, (Ptr{Void}, Ptr{UInt8}, TimeT, TimeT), fs.ptr, path, convert(TimeT, mtime), convert(TimeT, atime))

isequal(u1::HdfsURL, u2::HdfsURL) = isequal(u1.url, u2.url)
==(u1::HdfsURL, u2::HdfsURL) = (u1.url == u2.url)
hash(u::HdfsURL) = hash(u.url)

##
# File system implementation for HdfsFS
pwd(fs::HdfsFS) = URI(hdfs_pwd(fs)).path
function cd(fs::HdfsFS, path::AbstractString)
    hdfs_exists(fs, path) && (0 == hdfs_cd(fs, path)) && (return nothing)
    error("no such file or directory $path")
end
mkdir(fs::HdfsFS, path::AbstractString) = (0 == hdfs_mkdir(fs, path)) ? nothing : error("error creating directory $path")
mv(srcFS::HdfsFS, src::AbstractString, dstFS::HdfsFS, dst::AbstractString) = (srcFS == dstFS) ? mv(srcFS, src, dst) : (0 == hdfs_move(srcFS, src, dstFS, dst)) ? nothing : error("error moving file")
mv(fs::HdfsFS, src::AbstractString, dst::AbstractString) = (0 == hdfs_rename(fs, src, dst)) ? nothing : error("error renaming file")
rm(fs::HdfsFS, path::AbstractString) = (0 == hdfs_delete(fs, path)) ? nothing : error("error deleting $path")
cp(srcFS::HdfsFS, src::AbstractString, dstFS::HdfsFS, dst::AbstractString) = (0 == hdfs_copy(srcFS, src, dstFS, dst)) ? nothing : error("error copying file")
cp(fs::HdfsFS, src::AbstractString, dst::AbstractString) = cp(fs, src, fs, dst)
rmdir(fs::HdfsFS, path::AbstractString) = rm(fs, path)
function readdir(fs::HdfsFS)
    offset = length(pwd(fs))+2
    AbstractString[URI(fi.name).path[offset:end] for fi in hdfs_list_directory(fs, ".")]
end
function readdir(fs::HdfsFS, path::AbstractString)
    (path == ".") && (return readdir(fs))
    p = pwd(fs)
    cd(fs, path)
    ret = readdir(fs)
    cd(fs, p)
    ret
end

isdir(fs::HdfsFS, path::AbstractString) = hdfs_is_directory(fs, path)
function isdir(f::HdfsURL)
    uri = URI(f.url)
    isdir(hdfs_connect(f), uri.path)
end

##
# IO implementation for HdfsFile
# this is crap. should implement the array interface also
open(fs::HdfsFS, path::AbstractString, mode::AbstractString="r", buffer_sz::Integer=0, replication::Integer=0, bsz::Integer=0) = hdfs_open(fs, path, mode, buffer_sz, replication, bsz)
open(f::HdfsFile, mode::AbstractString="r", buffer_sz::Integer=0, replication::Integer=0, bsz::Integer=0) = hdfs_open(f, mode, buffer_sz, replication, bsz)
open(url::HdfsURL, mode::AbstractString="r", buffer_sz::Integer=0, replication::Integer=0, bsz::Integer=0) = hdfs_open(url.url, mode, buffer_sz, replication, bsz)
close(f::HdfsFile) = hdfs_close(f)
eof(f::HdfsFile) = (position(f) == filesize(f))

const _sb = Array(UInt8, 1)
read(f::HdfsFile, x::Type{UInt8}) = read(f, _sb)[1]
function read{T}(f::HdfsFile, a::Array{T})
    remaining = length(a)*sizeof(T)
    avlb = nb_available(f)
    while((remaining > 0) && (avlb > 0))
        ret = hdfs_read(f, pointer(a, length(a)-remaining+1), min(remaining, avlb))
        remaining -= ret
        avlb -= ret
    end
    (remaining > 0) && throw(EOFError())
    a
end
readbytes(f::HdfsFile, nb::Integer) = bytestring(read(f, Array(UInt8, nb)))
readall(f::HdfsFile) = readbytes(f, nb_available(f))

function peek(f::HdfsFile)
    eof(f) && (return 0xff)
    ret = read(f, UInt8)
    skip(f, -1)
    ret
end

write(f::HdfsFile, p::Ptr, nb::Integer) = write(f, p, convert(Int,nb))
write(f::HdfsFile, p::Ptr, nb::Int) = hdfs_write(f, p, nb)
write(f::HdfsFile, x::UInt8) = write(f, UInt8[x])
write{T}(f::HdfsFile, a::Array{T}, len) = write_sub(f, a, 1, length(a))
write{T}(f::HdfsFile, a::Array{T}) = write(f, a, length(a))
write_sub{T}(f::HdfsFile, a::Array{T}, offs, len) = isbits(T) ? write(f, pointer(a,offs), len*sizeof(T)) : error("$T is not bits type")

flush(f::HdfsFile) = (0 == hdfs_flush(f)) ? nothing : error("error flushing")
nb_available(f::HdfsFile) = hdfs_available(f)
function position(f::HdfsFile)
    p = hdfs_tell(f)
    (p >= 0) ? p : error("error getting current position")
end

stat(fs::HdfsFS, path::AbstractString) = hdfs_get_path_info(fs, path)
stat(f::HdfsFile) = stat(f.fs, f.path)
stat(url::HdfsURL) = stat(hdfs_connect(url), URI(url.url).path)

filesize(fs::HdfsFS, path::AbstractString) = (stat(fs,path)).size
filesize(f::HdfsFile) = filesize(f.fs, f.path)
filesize(url::HdfsURL) = filesize(hdfs_connect(url), URI(url.url).path)

seek(f::HdfsFile, n::Integer) = ((n >= 0) ? (0 == hdfs_seek(f, n)) : error("invalid position"))
seekend(f::HdfsFile) = seek(f, filesize(f))
seekstart(f::HdfsFile) = seek(f, 0)
skip(f::HdfsFile, n::Integer) = seek(f, n+position(f))


text:
# Libraries
include("./get_h10_dataset.jl");
using MessyTimeSeriesOptim;
using FileIO, JLD;
using Random, Statistics;

#=
Load arguments passed through the command line
=#

# VAR or VMA
is_var = parse(Bool, ARGS[1]);
log_folder_path = ARGS[2];

#=
Load benchmark grid of candidates from any run
=#

candidates = load("./VAR_output/err_type_2.jld")["candidates"];

#=
Presample: 53 weeks (1999) -> compute weights
Out-of-sample sample: from the 54th week onwards (2000 to 2020) -> run pseudo out-of-sample exercise
Test sample: from the 158th week onwards (2002 to 2020) -> compute realised error
=#

# Download estimation and test samples
df = get_h10_dataset(fred_tickers, mnemonics, transform, to_include, "1999-01-01", "2020-12-25");

# Extract relevant data
presample = df[1:53, 2:end] |> JMatrix{Float64};
oos_sample = df[54:end, 2:end] |> JMatrix{Float64};
test_sample = df[158:end, 2:end] |> JMatrix{Float64};
test_dates = df[158:end, :date];

# Periods not indicated as NBER recession in the test sample
nber_recessions = vcat(collect(Date(2007, 12, 07):Week(1):Date(2009, 05, 29)), Date(2020, 02, 07):Week(1):Date(2020, 03, 27));
test_periods_excl_recessions = [findfirst(test_dates .== x) for x in setdiff(test_dates, nber_recessions)];

# Remove NaNs
nan_to_missing!(presample);
nan_to_missing!(oos_sample);
nan_to_missing!(test_sample);

# Transpose data
presample = permutedims(presample);
oos_sample = permutedims(oos_sample);
test_sample = permutedims(test_sample);

# Compute benchmark weights
weights = 1 ./ (std_skipmissing(presample).^2);

#=
Out-of-sample exercise
=#

# Validation settings
model_kwargs = (tol=1e-3, check_quantile=true, verb=false);

# Dimensions
t0 = 158-54;
n, T = size(oos_sample);
grid_length = size(candidates, 2);

# One-step ahead forecast output
weighted_errors = zeros(size(candidates,2));
weighted_errors_excl_recessions = zeros(size(candidates,2));
forecast_per_series = Array{FloatMatrix,1}(undef, size(candidates,2));

# Setup data for estimation and forecast
estimation_sample = @view oos_sample[:, 1:t0];
estimation_sample_mean = mean_skipmissing(estimation_sample);
estimation_sample_std = std_skipmissing(estimation_sample);
zscored_oos_sample = (oos_sample .- estimation_sample_mean) ./ estimation_sample_std;

# Loop over each candidate vector of hyperparameters
for i in axes(candidates, 2)

    @info("Iteration $(i) out of $(grid_length)")

    # Current candidate vector of hyperparameters
    p_float, λ, α, β = candidates[:,i];
    p = Int64(p_float);

    # Select appropriate EstimSettings
    if is_var
        estim = VARSettings(zscored_oos_sample[:, 1:t0], p, λ, α, β; model_kwargs...);
    else
        estim = VMASettings(zscored_oos_sample[:, 1:t0], p, λ, α, β; model_kwargs...);
    end

    # Estimate model
    sspace = ecm(estim, output_sspace_data=zscored_oos_sample);

    # Run Kalman filter
    status = kfilter_full_sample(sspace);

    # Forecast
    X_prior = mapreduce(Xt -> sspace.B*Xt, hcat, status.history_X_prior);
    forecast_per_series[i] = (X_prior[:, t0+1:end] .* estimation_sample_std) .+ estimation_sample_mean;
    weighted_se = weights .* (test_sample .- forecast_per_series[i]).^2;
    weighted_errors[i] = MessyTimeSeriesOptim.compute_loss(weighted_se)[1];
    weighted_errors_excl_recessions[i] = MessyTimeSeriesOptim.compute_loss(weighted_se[:, test_periods_excl_recessions])[1];
end

# Save output to JLD
save("$(log_folder_path)/realised_error.jld", Dict("candidates" => candidates, "weighted_errors" => weighted_errors, "weighted_errors_excl_recessions" => weighted_errors_excl_recessions, "test_periods_excl_recessions" => test_periods_excl_recessions, "forecast_per_series" => forecast_per_series));


text:
include("common.jl")
@testset "Return code" begin
    @testset for composite in (true, false)
        alg = MethodOfSteps(composite ? AutoTsit5(Rosenbrock23()) : Tsit5();
                            constrained=false)

        sol1 = solve(prob_dde_1delay, alg)
        @test sol1.retcode == :Success

        sol2 = solve(prob_dde_1delay, alg; maxiters = 1, verbose = false)
        @test sol2.retcode == :MaxIters

        sol3 = solve(prob_dde_1delay, alg; dtmin = 5, verbose = false)
        @test sol3.retcode == :DtLessThanMin
    end
end


text:
import Base.diag
#=
diag : Output the diagonal of input Assoc A.
Outputs the Assoc with only the diagonal elements of A.
=#
function diag(A::Assoc)
    return deepCondense(Assoc(A.row,A.col,A.val,sparse(diagm(diag(A.A)))))
end

########################################################
# D4M: Dynamic Distributed Dimensional Data Model
# Architect: Dr. Jeremy Kepner (kepner@ll.mit.edu)
# Software Engineer: Alexander Chen (alexc89@mit.edu)
########################################################


text:
#Copyright Nathan Parrish 2018
#University of Colorado, Boulder
#Colorado Center for Astrodynamics Research
#Celestial and Spaceflight Mechanics Laboratory
#
#Derivatives for states & costates -- for CRTBP

#For use with NEW version of DifferentialEquations.jl
function CRTBP_stateCostate_deriv!(dstate_costate, statesCostates, params, t)
  #Constant mass.

  #pull out parameters
  (MU, DU, TU, thrustLimit, mass, time_direction, p, rho) = params

  λv = statesCostates[10:12]

  # States
  X1 = statesCostates[1]
  X2 = statesCostates[2]
  X3 = statesCostates[3]
  X4 = statesCostates[4]
  X5 = statesCostates[5]
  X6 = statesCostates[6]

  #Costates:
  L1 = statesCostates[7]
  L2 = statesCostates[8]
  L3 = statesCostates[9]
  L4 = statesCostates[10]
  L5 = statesCostates[11]
  L6 = statesCostates[12]

  accelLimit = thrustLimit / mass / 1e3 * TU^2 / DU #(kg*m/s^2) -> (km/s^2) -> (DU/TU^2) acceleration

  #Thrust magnitude from control law:
  if p == 0
    #Force thrust to always be on, in the optimal direction. This limits how big
    #of a search we need to make for the neural net stuff.
    umag = accelLimit #N

  elseif p == 1 #minimum mass problem
      g = norm(λv) - 1
      umag = 1/2 * (1 + tanh( g / (2*rho))) * accelLimit

  elseif p > 1 #p is in the range (1, 2]
      umag = (1/p * norm(λv))^(1 / (p-1));

      if umag > accelLimit
        umag = accelLimit
      end
  else
    error("Invalid value of p!")
  end

  #Control vector, scaled properly:
  control_accel = -umag * λv ./ norm(λv) #DU/TU^2

  if isnan(control_accel[1])
      control_accel[:] = 0.
  end

  #calculate some stuff that shows up a lot:
  #calculate r1^3 and r2^3 once each
  r1_3 = ((X1+MU)^2 + X2^2 + X3^2)^(3/2);
  r2_3 = ((X1+MU-1)^2 + X2^2 + X3^2)^(3/2);

  temp1 = ((MU + X1 - 1)^2 + X2^2 + X3^2);
  temp2 = ((MU + X1)^2 + X2^2 + X3^2);
  temp3 = (2*MU + 2*X1 - 2);

  #################################################### Derivatives:
  #Output in place:
  dstate_costate[1:3] = statesCostates[4:6]
  dstate_costate[4] = -(1-MU)*(X1+MU)/r1_3 - MU*(X1-1+MU)/r2_3 + 2*time_direction*X5 + X1 + control_accel[1];
  dstate_costate[5] = -(1-MU)*X2/r1_3 - MU*X2/r2_3 - 2*time_direction*X4 + X2 + control_accel[2];
  dstate_costate[6] = -(1-MU)*X3/r1_3 - MU*X3/r2_3 + control_accel[3];

  dstate_costate[7] = - L5*((3*MU*X2*temp3)/(2*temp1^(5/2)) - (3*X2*(MU - 1)*(2*MU + 2*X1))/(2*temp2^(5/2))) - L6*((3*MU*X3*temp3)/(2*temp1^(5/2)) - (3*X3*(MU - 1)*(2*MU + 2*X1))/(2*temp2^(5/2))) - L4*((MU - 1)/temp2^(3/2) - MU/temp1^(3/2) + (3*MU*(MU + X1 - 1)*temp3)/(2*temp1^(5/2)) - (3*(MU + X1)*(MU - 1)*(2*MU + 2*X1))/(2*temp2^(5/2)) + 1);
  dstate_costate[8] = L6*((3*X2*X3*(MU - 1))/temp2^(5/2) - (3*MU*X2*X3)/temp1^(5/2)) - L5*((MU - 1)/temp2^(3/2) - MU/temp1^(3/2) - (3*X2^2*(MU - 1))/temp2^(5/2) + (3*MU*X2^2)/temp1^(5/2) + 1) - L4*((3*MU*X2*(MU + X1 - 1))/temp1^(5/2) - (3*X2*(MU + X1)*(MU - 1))/temp2^(5/2));
  dstate_costate[9] = L6*(MU/temp1^(3/2) - (MU - 1)/temp2^(3/2) + (3*X3^2*(MU - 1))/temp2^(5/2) - (3*MU*X3^2)/temp1^(5/2)) + L5*((3*X2*X3*(MU - 1))/temp2^(5/2) - (3*MU*X2*X3)/temp1^(5/2)) - L4*((3*MU*X3*(MU + X1 - 1))/temp1^(5/2) - (3*X3*(MU + X1)*(MU - 1))/temp2^(5/2));
  dstate_costate[10] = 2*L5*time_direction - L1;
  dstate_costate[11] = - L2 - 2*L4*time_direction;
  dstate_costate[12] = -L3

end

max_stars_repo_name: travelingspaceman/LowThrustOpt

max_issues_repo_name: travelingspaceman/LowThrustOpt

max_forks_repo_name: travelingspaceman/LowThrustOpt


text:
struct TypedSignal{T}
    s::Signal
end

add_child!(arg::TypedSignal,s::Signal) = add_child!(arg.s,s)

(s::TypedSignal)() = pull!(s)

(ts::TypedSignal{T})((@nospecialize val)) where T = push!(ts.s, isa(val,T) ?
    val : convert(T,val))

pull!(ts::TypedSignal{T}) where T = begin
    val = pull!(ts.s)
    (typeof(val) == T) ? val : convert(T,val)
end::T

TypedSignal(val) = TypedSignal{typeof(val)}(Signal(val))

TypedSignal(s::Signal) = TypedSignal{typeof(s())}(s)

TypedSignal(f::Function,args...;kwargs...) = begin
    s = Signal(f,args...;kwargs...)
    ts = TypedSignal(s)
end

TypedSignal{T}(f::Function,args...;kwargs...) where T = TypedSignal{T}(Signal(f,args...;kwargs...))

Base.show(io::IO, ::MIME"text/plain",@nospecialize ts::TypedSignal{T}) where T = begin
    state_str = try
        convert(T,ts[])
        ""
    catch err
        "(type mismatch: $(typeof(ts[])))"
    end
    valid_str = valid(ts.s) ? "" : "(invalidated)"
    printstyled(io, "$(typeof(ts))"; bold = true,color = 200)
    print(io, "$valid_str $state_str \nvalue: ", ts[])
end

export TypedSignal

function getindex(ts::TypedSignal)
    value(ts.s)
end

import Base.setindex!
function setindex!(ts::TypedSignal,val)
    set_value!(ts.s, val)
end


text:
module entropy
export velocity_histogram, kB, entropy_2D_MB, plot_velocity_distr

using SmoothedParticles
using Plots: plot, savefig
using Ipopt
using JuMP

"""
	Histogram(xs::Vector{Float64}, ys::Vector{Float64},N,dx)

Histogram structure storing ``N`` x values ``xs`` with uniform bin width ``dx`` and ``N`` y values ``ys``.
"""
struct Histogram
	xs::Vector{Float64}
	ys::Vector{Float64}
	N::Int64
	dx::Float64
end

"""
	velocity_histogram(sys::ParticleSystem; v_max = 0, N = 10)

Building the histrogram of 2D velocities (norms) with ``v_max`` the maximum velocity in the histogram and ``N`` bins.
"""
function velocity_histogram(sys::ParticleSystem; v_max = 0.0, N = 100)::Histogram
	if v_max == 0.0 # if v_max = 0, find the maximum velocity of the particles
		for k in 1:length(sys.particles)
			v = norm(sys.particles[k].v)
			if v > v_max
				v_max = v
			end
		end
	end

	# Find the heights of the histogram bins
	dv = v_max/N # velocity increment between the bins
	vs = 0.:dv:v_max
	ns = zeros(length(vs))
	for k in 1:length(sys.particles)
		v = norm(sys.particles[k].v)
		n = Int64(round(v/dv))
		if 1 <= n <= length(ns)
			ns[n] += 1.0/(dv*length(sys.particles))
		end
	end

	return Histogram(vs,ns,100,dv)
end

"""
	kB

Boltzmann constant (in the SI units)
"""
const kB = 1.380649e-23

"""
	entropy(fMB::Histogram)::Float64

Calculate Boltzmann entropy of a 2D Maxwell-Boltzmann distribution approximated by an ``fMB`` histogram.
"""
function entropy_2D_MB(fMB::Histogram)::Float64
	@assert(fMB.xs[1] == 0) # Assuming that the histogram starts at zero velocity

	S = 0.0

	# Approximating the reduced entropy near v=0, where a numerical singularity could appear
	fMBder = (fMB.ys[2]-fMB.ys[1])/fMB.dx
	if fMBder > 0
		S = - fMB.ys[1] * (log(fMBder)*fMB.dx - fMBder*(fMB.dx^3)/6)
	end

	# Approximating the rest of entropy
	for k in 2:length(fMB.xs)
		if fMB.xs[k] != 0
			if fMB.ys[k] > 0
				S += -fMB.ys[k] * log(fMB.ys[k]/fMB.xs[k]) * fMB.dx
			end
		end
	end

	return S
end

"""
	plot_velocity_distr(sys::ParticleSystem, name::String)

Plots the distribution of velocity magnitudes among the particles of ``sys`` and saves the resulting pdf,
together with a fit of the Maxwell-Boltzmann distribution, to file ``name``. ``m`` is the mass of each particle.
Returns the temperature.
"""
function plot_velocity_distr(sys::ParticleSystem, m::Float64, name::String; v_max = 0.0)::Float64
	distr = velocity_histogram(sys, v_max=v_max, N = 100)

	# fitting the histogram to a 2D Maxwell-Boltzmann distribution
	model = Model(Ipopt.Optimizer)
	@variable(model, beta)
	@NLobjective(
			model,
			Min,
			sum((distr.ys[i] - m*beta*distr.xs[i]*exp(-0.5*m*beta*distr.xs[i]^2))^2 for i in 1:length(distr.xs)),
		)
	optimize!(model)
	beta = value(beta)

	# Plotting both the actual histogram and the fitted Maxwell-Boltzmann distribution
	ns_boltz = zeros(length(distr.xs))
	for i in 1:length(ns_boltz)
		ns_boltz[i] = m*beta*distr.xs[i]*exp(-0.5*m*beta*distr.xs[i]^2)
	end
	@show(beta)
	T = 1/(beta*kB)
	@show(T)

	p = plot(distr.xs, [distr.ys ns_boltz], label = ["data" "Maxwell-Boltzmann, T="*string(T)])
	savefig(p, name)
	return T
end

#"""
#	plot_velocity_distr(path::String, name:: String)
#
#Plot velocity distribution of particles loaded from a VTK file.
#"""
#function plot_velocity_distr(path::String, name::String)
#	@error "Not working, sorry"
#	domain = Rectangle(-box_width, -box_width, 2*box_width, 3*box_height)
#	sys = ParticleSystem(Particle, domain, h)
#	read_vtk!(sys, path, x -> Particle(x = x, type = 0.0))
#	plot_velocity_distr(sys, name)  #not yet working
#end

end


text:
struct BASH <: AbstractShell
    path::AbstractString
    prompt::AbstractString
end

BASH(path::AbstractString) = BASH(path, "\\e[0;35m\\u@\\h:\\W (playground)> \\e[m")
BASH() = BASH(strip(readstring(`which bash`)))

function Base.run(shell::BASH, env::Environment)
    init(env)
    prompt = getprompt(shell, env)
    ENV["PS1"] = prompt

    usr_rc = join(home(), ".bashrc")
    pg_rc = join(parent(Path(ENV["JULIA_PKGDIR"])), ".bashrc")

    if !exists(pg_rc)
        exists(usr_rc) ? cp(usr_rc, pg_rc, follow_symlinks=true) : touch(pg_rc)

        content = string(
            "\nexport PATH=\"", ENV["PATH"], "\"\n",
            "export PS1=\"$prompt\"\n",
            "export JULIA_PKGDIR=\"", ENV["JULIA_PKGDIR"], "\"\n",
            "export JULIA_DEPOT_PATH=\"", ENV["JULIA_DEPOT_PATH"], "\"\n",
        )

        if haskey(ENV, "HISTFILE")
            content = string(content, "export HISTFILE=\"", ENV["HISTFILE"], "\"\n")
        end

        write(pg_rc, content, "a")
    end

    @mock run(`$(shell.path) --rcfile $pg_rc -i`)
end


text:
module HearThereTests

using HearThere
using DataFrame
using Base.Test

function generate_rawdata(n::Integer)
    generate_random_rawdata(n, 100, 0, 100)
end

function generate_rawdata(n::Integer, μ::Real, σ::Real)
    generate_random_rawdata(n, fill(μ, 9), fill(σ, 9))
end

function generate_rawdata(n::Integer, μ::Array{Real}, σ::Array{Real})
    DataFrame(
        timestamp=linspace(0, 10, n),
        accel_x=convert(Array{Int32}, map(round, randn(n) * σ[1] + μ[1])),
        accel_y=convert(Array{Int32}, map(round, randn(n) * σ[2] + μ[2])),
        accel_z=convert(Array{Int32}, map(round, randn(n) * σ[3] + μ[3])),
        gyro_x=convert(Array{Int32}, map(round, randn(n) * σ[4] + μ[4])),
        gyro_y=convert(Array{Int32}, map(round, randn(n) * σ[5] + μ[5])),
        gyro_z=convert(Array{Int32}, map(round, randn(n) * σ[6] + μ[6])),
        mag_x=convert(Array{Int32}, map(round, randn(n) * σ[7] + μ[7])),
        mag_y=convert(Array{Int32}, map(round, randn(n) * σ[8] + μ[8])),
        mag_z=convert(Array{Int32}, map(round, randn(n) * σ[9] + μ[9])))
end

function test_gyro_bias()
    n = 1000
    df = generate_rawdata(n, 2, 100)
end

test_gyro_bias()

end


text:
@testset "781.rabbits-in-forest.jl" begin
    @test num_rabbits([1,0,1,0,0]) == 5
    @test num_rabbits([1, 1, 2]) == 5
    @test num_rabbits([10, 10, 10]) == 11
    @test num_rabbits([0,0,1,1,1]) == 6
end


text:
function initialize_tpu!(sess)
    topo = run(sess, TensorFlow.Ops.configure_distributed_tpu())
    readproto(IOBuffer(topo), TopologyProto())
end

topo_string(topo::TopologyProto) = join(topo.mesh_shape, "x")

function permuted_coords(topo)
    permutedims(reshape(topo.device_coordinates, (length(topo.mesh_shape),
        topo.num_tpu_devices_per_task, topo.num_tasks)), (3,2,1))
end

function Base.show(io::IO, topo::TopologyProto)
    coordinates = permuted_coords(topo)
    print(io, topo_string(topo), " TPU Topology\n")
    for i = 1:topo.num_tasks
        println(" Task $(i-1):  ( x, y, core)")
        for j = 1:topo.num_tpu_devices_per_task
            println("  Device $(j-1): ", join(coordinates[i, j, :], "  "))
        end
    end
end

struct MeshIndex
    x::Int
    y::Int
    core::Int
end

struct DeviceIndex
    task::Int
    index::Int
end

function Base.getindex(topo::TopologyProto, idx::DeviceIndex)
    MeshIndex(permuted_coords(topo)[idx.task+1, idx.index+1, :]...)
end

function Base.getindex(topo::TopologyProto, idxs::Vector{DeviceIndex})
    coords = permuted_coords(topo)
    [MeshIndex(coords[idx.task+1, idx.index+1, :]...) for idx in idxs]
end

job_name(sess) = "tpu_worker"

function tf_tpu_device(sess, device::DeviceIndex)
    TensorFlow.Device("/job:$(job_name(sess))/replica:1/task:$(device.task+1)/device:TPU:$(device.index+1)")
end

function tf_host_cpu_device(sess, device::DeviceIndex)
    TensorFlow.Device("/job:$(job_name(sess))/replica:1/task:$(device.task+1)/device:CPU:1")
end


text:
using WAV
using ArgParse

push!(LOAD_PATH, pwd())

using Decoder

# Main entry
function main(wav_path)

	# Read the WAV file
	s, fs = wavread(wav_path)

	sequence = decode(s, fs)

	print("The dialing code is: ")
	for s in sequence
		print(s)
	end
	println()
end

if PROGRAM_FILE == "SignalDecoding.jl"

	# Parse the argument
	s = ArgParseSettings("Audio dialing decoder")

	@add_arg_table s begin
		"--audio"
		help = "The audio to be decoded"
		required = true
	end

	parsed_args = parse_args(ARGS, s)

	main(parsed_args["audio"])
end

max_stars_repo_name: PatricioIribarneCatella/TP-Telefono

max_issues_repo_name: PatricioIribarneCatella/TP-Telefono

max_forks_repo_name: PatricioIribarneCatella/TP-Telefono