text:
## GFF3 Writer

struct Writer{S <: TranscodingStream} <: BioGenerics.IO.AbstractWriter
    output::S
end

function BioGenerics.IO.stream(writer::Writer)
    return writer.output
end

function Writer(output::IO)
    if output isa TranscodingStream
        return Writer{typeof(output)}(stream)
    else
        stream = TranscodingStreams.NoopStream(output)
        return Writer{typeof(stream)}(stream)
    end
end

function Base.write(writer::Writer, record::Record)
    printgff(writer.output, record)
end

function gffstring(gene::Gene)
    buf = IOBuffer()
    firstattribute = true
    for field in names(parent(gene).genedata)
        field in [:source, :score, :phase] && continue
        v = parent(gene).genedata[index(gene), field]
        if !ismissing(v)
            if firstattribute
                firstattribute = false
            else
                print(buf, ";")
            end
            if v isa AbstractVector
                print(buf, field, "=")
                for i in eachindex(v)
                    print(buf, oneline(v[i]))
                    i == lastindex(v) ? print(buf, ";") : print(buf, ",")
                end
            else
                print(buf, field, "=", oneline(v))
            end
        end
    end
    join([parent(gene).name,
        get(gene, :source, "."),
        feature(gene),
        locus(gene).position.start,
        locus(gene).position.stop,
        get(gene, :score, "."),
        locus(gene).strand,
        get(gene, :phase, "."),
        String(take!(buf))], '\t')
end

"""
    printgff(io::IO, chr)
    printgff(path::AbstractString, chr)

Print `chr` in GFF3 format.
"""
printgff(filepath::AbstractString, chrs) = printgff(open(filepath, "w"), chrs)
printgff(io::IO, chr::Record) = printgff(io, [chr])
function printgff(io::IO, chrs::AbstractVector{Record{G}}) where G <: AbstractGene
    iobuffer = IOBuffer()
    ### Header
    if chrs[1].header[1:15] == "#gff-version 3\n"
        print(iobuffer, chrs[1].header)
    else
        println(iobuffer, "##gff-version 3")
    end
    ### Body
    for chr in chrs
        for gene in chr.genes
            println(iobuffer, gffstring(gene))
        end
    end
    ### Footer
    if !all(isempty(chr.sequence) for chr in chrs)
        println(iobuffer, "##FASTA")
        for chr in chrs
            println(iobuffer, ">", chr.name)
            for s in Iterators.partition(chr.sequence, 80)
                println(iobuffer, join(s))
            end
        end
    end
    print(io, String(take!(iobuffer)))
end


text:
export hessian_sparse
function hessian_sparse(nlp, x)
    n = length(x)
    H = hess(nlp, x)
    tempH = (H + tril(H,-1)')
    H = tempH
    return H
end

ext: jl


text:
abstract type Unit end

mutable struct Goblin <: Unit
    hp::Int64
    ap::Int64
end

mutable struct Elf <: Unit
    hp::Int64
    ap::Int64
end

mutable struct Empty <: Unit end
mutable struct Wall <: Unit end

mutable struct Node{T}
    loc::Tuple{Int64,Int64}
    data::T
    neighbors::Dict{String,Node{T}}
end

mutable struct Graph{T}
    data::Dict{Tuple{Int64,Int64},Node{T}}
    width::Int64
    height::Int64
end

function make_grid()
    f = open("inputs/15_test.txt", "r")
    lines = readlines(f)
    ysize = length(lines)
    xsize = maximum(length.(lines))
    d::Dict{Tuple{Int64,Int64},Node{Unit}} = Dict()
    graph = Graph(d,xsize,ysize)
    for row in 1:ysize
        for col in 1:xsize
            neighbors::Dict{String,Node{Unit}} = Dict()
            newnode = Node((col,row),Empty(),neighbors)
            graph.data[(col,row)] = newnode
            if row > 1
                # connect above
                nodeabove = graph.data[(col,row-1)]
                nodeabove.neighbors["south"] = newnode
                newnode.neighbors["north"] = nodeabove
            end
            if col > 1
                # connect left
                nodeleft = graph.data[(col-1,row)]
                nodeleft.neighbors["east"] = newnode
                newnode.neighbors["west"] = nodeleft
            end
        end
    end
    for y in 1:length(lines)
        for x in 1:length(lines[y])
            c = lines[y][x]
            if c == '#'
                # remove edges
                node = graph.data[(x,y)]
                try
                    n = node.neighbors["north"]
                    delete!(n.neighbors, "south")
                    delete!(node.neighbors, "north")
                catch
                end
                try
                    n = node.neighbors["south"]
                    delete!(n.neighbors, "north")
                    delete!(node.neighbors, "south")
                catch
                end
                try
                    n = node.neighbors["west"]
                    delete!(n.neighbors, "east")
                    delete!(node.neighbors, "west")
                catch
                end
                try
                    n = node.neighbors["east"]
                    delete!(n.neighbors, "west")
                    delete!(node.neighbors, "east")
                catch
                end
                delete!(graph.data, (x,y))
            elseif c == 'G'
                goblin = Goblin(200,3)
                graph.data[(x,y)].data = goblin
            elseif c == 'E'
                elf = Elf(200,3)
                graph.data[(x,y)].data = elf
            end
        end
    end

    return graph
end

function path(currentloc::Tuple{Int64,Int64},visited::Set{Tuple{Int64,Int64}},target::Tuple{Int64,Int64},grid::Graph{Unit})::Array{Tuple{Int64,Int64}}
    options = []
    try
        north = grid.data[currentloc].neighbors["north"]
        if north.loc == target
            return [north.loc]
        elseif !(in(north.loc, visited))
            v2 = copy(visited)
            push!(v2,currentloc)
            push!(vcat([currentloc], path(north.loc, v2, target, grid)), options)
        end
    catch
    end
    try
        west = grid.data[currentloc].neighbors["west"]
        if west.loc == target
            return [west.loc]
        elseif !(in(west.loc, visited))
            v2 = copy(visited)
            push!(v2,currentloc)
            push!(vcat([currentloc], path(west.loc, v2, target, grid)), options)
        end
    catch
    end
    try
        east = grid.data[currentloc].neighbors["east"]
        if east.loc == target
            return [east.loc]
        elseif !(in(east.loc, visited))
            v2 = copy(visited)
            push!(v2,currentloc)
            push!(vcat([currentloc], path(east.loc, v2, target, grid)), options)
        end
    catch
    end
    try
        south = grid.data[currentloc].neighbors["south"]
        if south.loc == target
            return [south.loc]
        elseif !(in(south.loc, visited))
            v2 = copy(visited)
            push!(v2,currentloc)
            push!(vcat([currentloc], path(south.loc, v2, target, grid)), options)
        end
    catch
    end
    shortest = sort(map(x->(length(x),x), options), by=x->x[1])[1]
    return shortest
end

function _empty(x)
    if x == (0,0)
        return true
    else
        return false
    end
end

function update_grid(direction::String,coords::Tuple{Int64,Int64},paths::Array{Tuple{Int64,Int64},2},grid::Graph, target::Tuple{Int64,Int64})
    try
        node = grid.data[coords].neighbors[direction]
        if paths[node.loc...] == (0,0) && (node.loc == target || typeof(grid.data[node.loc...].data) == Empty)
            paths[node.loc...] = coords
            return true
        end
        return false
    catch
        return false
    end
end

function _path(target::Tuple{Int64,Int64},start::Tuple{Int64,Int64},grid::Graph{Unit})
    paths::Array{Tuple{Int64,Int64},2} = map(x->(0,0),zeros(grid.height,grid.width))
    paths[start...] = (-1,-1)
    active = [start]
    while length(active) > 0
        n = active[1]
        deleteat!(active,1)
        for direction in ["west", "north", "south", "east"]
            result = update_grid(direction, n, paths, grid, target)
            if result
                push!(active, grid.data[n].neighbors[direction].loc)
            end
        end
    end
    if paths[target...] == (0,0)
        return []
    end
    retval = []
    push!(retval, target)
    while retval[end] != (-1,-1)
        push!(retval, paths[retval[end]...])
    end
    return retval[2:end-1]
end

function nearest_target(x::Tuple{Int64,Int64},grid::Graph)
    target = Elf
    if typeof(grid.data[x].data) == Elf
        target = Goblin
    end
    println(target)
    targets = []
    for pair in grid.data
        if typeof(pair[2].data) == target
            push!(targets, pair[1])
        end
    end
    paths = [_path(x,y,grid) for y in targets]
    distances = [length(x) for x in paths]
    if length(distances) == 0
        return x
    end
    attack_positions = []
    if minimum(distances) == 1
        # skip move
        return x
    else
        # find closest attack position
        println("closing in")
        for target in targets
            push!(attack_positions, (target[1]+1,target[2]))
            push!(attack_positions, (target[1],target[2]+1))
            push!(attack_positions, (target[1],target[2]-1))
            push!(attack_positions, (target[1]-1,target[2]))
        end
        attack_positions = filter(y->length(_path(x,y,grid)) > 0, attack_positions)
        println(attack_positions)
    end
    closest_distance = minimum(map(y->length(_path(x,y,grid)), attack_positions))
    println(closest_distance)
    positions = filter(y->length(_path(x,y,grid))==closest_distance, attack_positions)
    println(positions)
    if length(positions) == 0
        return x
    end
    position_to_go_to = positions[end]
    return position_to_go_to
end

function move(x::Tuple{Int64,Int64}, grid::Graph)
    println("moving " * string(x))
    nt = nearest_target(x,grid)
    draw_grid(grid)
    if x != nt
        println("Actually moving")
        _p = _path(x,nt,grid)
        println(_p)
        p = _p[1]
        u = grid.data[x].data
        grid.data[p].data = u
        grid.data[x].data = Empty()
    end
    println("...")
    draw_grid(grid)
end

function _get_units(grid::Graph)
    units = []
    for row in 1:grid.height
        for col in 1:grid.width
            try
                n = grid.data[(row,col)]
                if typeof(n.data) == Elf || typeof(n.data) == Goblin
                    push!(units, (row,col))
                end
            catch
            end
        end
    end
    return units
end

function move_units(grid::Graph)
    types_left = Set()
    for unit in _get_units(grid)
        push!(types_left, typeof(grid.data[unit].data))
    end
    if length(collect(types_left)) > 1
        for unit in _get_units(grid)
            move(unit, grid)
        end
    end
end

function draw_grid(grid::Graph)
    retval = ""
    for row in 1:grid.height
        for col in 1:grid.width
            try
                if typeof(grid.data[(row,col)].data) == Wall
                    retval *= '#'
                end
                if typeof(grid.data[(row,col)].data) == Empty
                    retval *= '.'
                end
                if typeof(grid.data[(row,col)].data) == Goblin
                    retval *= 'G'
                end
                if typeof(grid.data[(row,col)].data) == Elf
                    retval *= 'E'
                end
            catch
                retval *= '#'
            end
        end
        retval *= '\n'
    end
    print(retval)
end

grid = make_grid();
move_units(grid)
draw_grid(grid)


text:
struct SSHAgent
    handle::Ptr{Cvoid}
    session::SSHSession
end

Base.convert(::Type{Ptr{LIBSSH2_AGENT}}, ag::SSHAgent) = ag.handle
Base.unsafe_convert(::Type{Ptr{LIBSSH2_AGENT}}, ag::SSHAgent) = convert(Ptr{LIBSSH2_AGENT}, ag.handle)

SSHAgent(session::SSHSession) = SSHAgent(libssh2_agent_init(session), session)

Sockets.connect(agent::SSHAgent) = libssh2_agent_connect(agent) == 0
disconnect(agent::SSHAgent) = libssh2_agent_disconnect(agent) == 0
free(agent::SSHAgent) = libssh2_agent_free(agent)


text:
using LaserTag
using POMDPToolbox
using ParticleFilters
using ContinuousPOMDPTreeSearchExperiments
using QMDP
using JLD
using BasicPOMCP
using StaticArrays
using POMDPs
using Reel
using ARDESPOT

f = ARGS[1]
k = ARGS[2]
i = parse(Int, ARGS[3])

rdict = load(f, "rdict")

n = 100_000

pomdp = gen_lasertag(rng=MersenneTwister(i+600_000))

sol = DESPOTSolver(lambda=0.1,
                   T_max=5.0,
                   max_trials=1_000_000,
                   bounds=LaserBounds{typeof(pomdp)}(),
                   rng=MersenneTwister(4))

if isa(sol,Solver)
    p = solve(deepcopy(sol), pomdp)
else
    p = sol
end
hr = HistoryRecorder(max_steps=100, rng=MersenneTwister(i))
up_rng = MersenneTwister(i+100_000)
# up = ObsAdaptiveParticleFilter(deepcopy(pomdp), LowVarianceResampler(10_000), 0.05, up_rng)
up = ObsAdaptiveParticleFilter(deepcopy(pomdp), LowVarianceResampler(100_000), 0.05, up_rng)

frames = Frames(MIME("image/png"), fps=2)
disc = 1.0
rew = 0.0
states = []
beliefs = []
for (s, b, a, r, sp, o, bp) in stepthrough(pomdp, p, up, "sbarspobp", max_steps=100, rng=MersenneTwister(i))
    rew += disc*r
    disc*=discount(pomdp)
    show(STDOUT, MIME("text/plain"), LaserTagVis(pomdp, s=s))
    println()
    push!(frames, LaserTagVis(pomdp, s=sp, b=bp, o=o, a=a))
    push!(states, s)
    push!(beliefs, b)
end

println("Total discounted reward for this simulation: $rew")
println("Total discounted reward from original: $(rdict[k][i])")

gifname = "out"*randstring()*".gif"
println("Writing $gifname...")
write(gifname, frames)
println("done.")


text:
using ITensors
using OptimKit
using Zygote
using Random: seed!

include("circuit.jl")

function gate(::Val{:Ry}; θ)
  return [
    cos(θ / 2) -sin(θ / 2)
    sin(θ / 2) cos(θ / 2)
  ]
end

function gate(::Val{:CX})
  return [
    1 0 0 0
    0 1 0 0
    0 0 0 1
    0 0 1 0
  ]
end

function Rylayer(N, θ⃗)
  return [("Ry", (n,), (θ=θ⃗[n],)) for n in 1:N]
end

function CXlayer(N)
  return [("CX", (n, n + 1)) for n in 1:2:(N - 1)]
end

layer(N, θ⃗) = vcat(Rylayer(N, θ⃗), CXlayer(N))

# The variational circuit we want to optimize
function variational_circuit(N, nlayers, θ⃗)
  range = 1:N
  circuit = layer(N, θ⃗[range])
  for n in 1:(nlayers - 1)
    circuit = vcat(circuit, layer(N, θ⃗[range .+ n * N]))
  end
  return circuit
end

N = 4
s = siteinds("Qubit", N)

h = 1.3
ℋ = OpSum()
for j in 1:(N - 1)
  ℋ .+= -1, "Z", j, "Z", j + 1
end
for j in 1:N
  ℋ .+= h, "X", j
end
H = prod(MPO(ℋ, s))
ψ = prod(MPS(s, "0"))

N = 4
nlayers = 5

function loss(θ⃗)
  gates = variational_circuit(N, nlayers, θ⃗)
  U = buildcircuit(gates, s)
  return rayleigh_quotient(H, (U, ψ))
end

seed!(1234)
θ⃗₀ = 2π * rand(N * nlayers)
@show θ⃗₀
@show loss(θ⃗₀)
@show loss'(θ⃗₀)

fg(x) = (loss(x), convert(Vector, loss'(x)))
θ⃗ₒₚₜ, fₒₚₜ, gₒₚₜ, numfg, normgradhistory = optimize(fg, θ⃗₀, GradientDescent())
@show θ⃗ₒₚₜ
@show fₒₚₜ
@show gₒₚₜ
@show numfg
println("normgradhistory = ")
display(normgradhistory)

Uₒₚₜ = buildcircuit(variational_circuit(N, nlayers, θ⃗ₒₚₜ), s)
Uₒₚₜψ = apply(Uₒₚₜ, ψ)
@show (Uₒₚₜψ' * H * Uₒₚₜψ)[]

nothing


text:
@testset "Statistics" begin
  # half sample mode
  d = LogNormal(0,1)
  rs = rand(d,1000)
  @test GeoStatsBase.hsm_mode([1,2,2,3]) == 2.0
  @test GeoStatsBase.hsm_mode([1,2,2,3,5]) == 2.0
  @test GeoStatsBase.hsm_mode(rs) < mean(rs)
  @test GeoStatsBase.hsm_mode(rs) < median(rs)
  d = MixtureModel([Normal(), Normal(3, 0.2)], [0.7, 0.3])
  rs = rand(d, 1000)
  @test GeoStatsBase.hsm_mode(rs) < mean(rs)
  @test GeoStatsBase.hsm_mode(rs) < median(rs)

  # load data with bias towards large values (gold mine)
  sdata = georef(CSV.File(joinpath(datadir,"clustered.csv")), (:x,:y))

  # spatial mean
  μn = mean(sdata[:Au])
  μs = mean(sdata, :Au)
  @test abs(μn - 0.5) > abs(μs - 0.5)
  @test mean(sdata)[:Au] ≈ μs

  # spatial variance
  σn = var(sdata[:Au])
  σs = var(sdata, :Au)
  @test isapprox(σn, σs, atol=1e-2)
  @test var(sdata)[:Au] ≈ σs

  # spatial quantile
  qn = quantile(sdata[:Au], 0.5)
  qs = quantile(sdata, :Au, 0.5)
  @test qn ≥ qs
  @test quantile(sdata, 0.5)[:Au] ≈ qs
end


text:
# Bifurcation diagram for 1D climate model
using DrWatson
@quickactivate "NonlinearDynamicsTextbook"
include(srcdir("style.jl"))
using DynamicalSystems, PyPlot

include(srcdir("simple_continuation.jl"))

using ForwardDiff, PyPlot
αtan(T) = 0.5 - 0.2*tanh((T-263)/4)
dTdt(T, ε = 0.65, α=αtan, s= 1.0) = s*(1 - α(T)) - 1.6e-10 * ε * T^4
dTdt(T; ε = 0.65, α=αtan, s = 1.0) = dTdt(T, ε, α, s)
d²Tdt²(T, ε) = ForwardDiff.derivative(T -> dTdt(T, ε), T)
climate_f(u, p, t) = SVector(dTdt(u[1], p))
climate_J(u, p, t) = SMatrix{1,1}(d²Tdt²(u[1], p))

u0 = SVector(225.0)
p0 = 0.9
ds = ContinuousDynamicalSystem(climate_f, u0, p0, climate_J)

f, J = bifurcation_rule_form(ds)

pmin = 0.3
pmax = 1

dp0 = 0.05
dx0 = SVector(0.1)

xs, ps, stability = continuation(f, J, u0, p0;
    pmin, pmax, dp0, dx0, N = 1500
)

xs = [x[1] for x in xs]

# %% make figure
stable1 = findfirst(!, stability)
unstable1 = findlast(!, stability)

fig = figure(figsize = (figx/2, figy))
ax = gca()

plot(ps[1:stable1], xs[1:stable1]; c = "C2", ls = "solid", label="stable")
plot(ps[stable1+1:unstable1], xs[stable1+1:unstable1]; c = "C2", ls = "dashed", label="unstable")
plot(ps[unstable1+1:end], xs[unstable1+1:end]; c = "C2", ls = "solid")

# scatter bifurcation points
ax.plot([ps[stable1], ps[unstable1]], [xs[stable1], xs[unstable1]], ls = "None",
		marker="o", ms = 12, color = "k", mew = 1, mec = "k",
        fillstyle = "top", markerfacecoloralt= "w"
)

# plot arrows
_a = 300
arrow1ε = ps[_a:stable1] .- 0.03
arrow1T = xs[_a:stable1] .- 5
plot(arrow1ε, arrow1T; color = "C1", lw = 4.0)
arrow(arrow1ε[end], arrow1T[end], 0, 50;
width = 0.004, head_length = 5.0, color = "C1", zorder = 99)

_b = 200
arrow2ε = ps[unstable1:unstable1+_b] .+ 0.03
arrow2T = xs[unstable1:unstable1+_b] .+ 5

plot(arrow2ε, arrow2T; color = "C0", lw = 4.0)
arrow(arrow2ε[1], arrow2T[1], 0, -50;
width = 0.004, head_length = 5.0, color = "C0", zorder = 99)

# scatter(ps[stability], xs[stability]; c = "C3", legend = "stable")
# scatter(ps[stability], xs[stability]; c = "C3", legend = "stable")
legend()
xticks(0.3:0.2:0.9)
yticks(200:50:350)
ylabel("\$T^*\$", labelpad = -20)
xlabel("\$\\varepsilon\$"; labelpad=-20)
# yticks(230:50:310)
# ylim(210, 320)
# ax.xaxis.set_label_coords(0.4, -0.025)
fig.tight_layout(pad=0.3)
wsave(plotsdir("4", "bif_example"), fig)

max_stars_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_issues_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_forks_repo_name: JuliaDynamics/NonlinearDynamicsTextbook


text:
function train_MNIST_model(model, trainX, trainY, valX, valY)
    model(trainX[1]);

    loss(x, y) = my_xent(model(x), y)
    accuracy(x, y) = mean(onecold(model(x)) .== onecold(y))

    evalcb = throttle(() -> @show(loss(valX[1], valY[1])) , 5);
    opt = ADAM(1e-3);
    @epochs 10 mytrain!(loss, Flux.params(model), zip(trainX, trainY), opt, cb=evalcb)
end

function advtrain!(model, attack_fn, loss, ps, data, opt; cb = () -> ())
    ps = Flux.Tracker.Params(ps)
    cb = runall(cb)
    m_cleanloss = MeanMetric()
    m_cleanacc = MeanMetric()
    m_advloss = MeanMetric()
    m_advacc = MeanMetric()
    step = 0
    @showprogress 0.1 "Training..." for d in data
        # FIXME can I use x,y in for variable?
        x, y = d
        x_adv = attack_fn(model, d...)

        # train xent loss using adv data
        gs = Flux.Tracker.gradient(ps) do
            clean_logits = model(x)
            clean_loss = my_xent(clean_logits, y)
            adv_logits = model(x_adv)
            adv_loss = my_xent(adv_logits, y)

            add!(m_cleanloss, clean_loss.data)
            add!(m_cleanacc, accuracy_with_logits(clean_logits.data, y))
            add!(m_advloss, adv_loss.data)
            add!(m_advacc, accuracy_with_logits(adv_logits.data, y))

            l = adv_loss
            # l = adv_loss + clean_loss
            l
        end
        Flux.Tracker.update!(opt, ps, gs)

        step += 1

        if step % 40 == 0
            println()
            @show get!(m_cleanloss)
            @show get!(m_cleanacc)
            @show get!(m_advloss)
            @show get!(m_advacc)
        end
        # cb(step, total_loss)
        cb()
    end
end

function advtrain(model, attack_fn, trainX, trainY, valX, valY)
    model(trainX[1]);

    # evalcb = throttle(cb_fn , 5);
    evalcb = () -> ()

    # train
    opt = ADAM(1e-4);
    @epochs 3 advtrain!(model, attack_fn, Flux.params(model), zip(trainX, trainY), opt, cb=evalcb)
end

function evaluate_attack(model, attack_fn, trainX, trainY, testX, testY)
    accuracy(x, y) = mean(onecold(model(x)) .== onecold(y))
    println("Clean accuracy:")
    # FIXME this may need to be evaluated on CPU
    @show accuracy(trainX[1], trainY[1])
    @show accuracy(testX[1], testY[1])

    sample_and_view(testX[1], testY[1])

    x = testX[1]
    y = testY[1]

    @info "performing attack .."
    x_adv = attack_fn(model, x, y)
    @info "attack done."

    @show accuracy(x_adv, y)

    sample_and_view(x_adv, model(x_adv))

    # all test data
    m_acc = MeanMetric()
    @showprogress 0.1 "testing all data .." for d in zip(testX, testY)
        x, y = d
        x_adv = attack_fn(model, x, y)
        acc = accuracy(x_adv, y)
        add!(m_acc, acc)
    end
    @show get(m_acc)
    nothing
end

function evaluate_AE(ae, cnn, ds)
    # test clean cnn accuracy
    accuracy(x, y) = mean(onecold(cnn(x)) .== onecold(y))
    println("Computing clean accuracy ..")
    @show accuracy(trainX[1], trainY[1])
    @show accuracy(testX[1], testY[1])
    sample_and_view(testX[1], testY[1])
    # use ae and visualize the images
    println("Visualizing decoded images ..")
    decoded = ae(testX[1]);
    maximum(decoded)
    minimum(decoded)
    # FIXME use .data because decoded is TrackedArray, why???
    sample_and_view(decoded.data, testY[1])
    sample_and_view(testX[1], testY[1])
    # add ae and test all accuracy
    ae_accuracy(x, y) = mean(onecold(cnn(ae(x))) .== onecold(y))
    println("Computing AE accuracy ..")
    @show ae_accuracy(trainX[1], trainY[1])
    @show ae_accuracy(testX[1], testY[1])
end

function upsample(x)
    ratio = (2,2,1,1)
    y = similar(x, (size(x) .* ratio)...)
    for i in Iterators.product(Base.OneTo.(ratio)...)
        loc = map((i,r,s)->range(i, stop = s, step = r), i, ratio, size(y))
        @inbounds y[loc...] = x
    end
    y
end

function AdvAE_train(ae, cnn, attack_fn, trainX, trainY, valX, valY)
    model = Chain(ae, cnn)
    model(trainX[1]);

    loss(x, y) = crossentropy(model(x), y)
    # how about try add regularizer
    # loss(x, y) = crossentropy(cnn(ae(x)), y) + mymse(ae(x), x)

    accuracy(x, y) = mean(onecold(model(x)) .== onecold(y))
    adv_accuracy(x, y) = begin
        x_adv = attack_fn(model, loss, x, y)
        accuracy(x_adv, y)
    end
    adv_loss(x, y) = begin
        x_adv = attack_fn(model, loss, x, y)
        loss(x_adv, y)
    end
    cb_fn() = begin
        # add a new line so that it plays nicely with progress bar
        @time begin
            println("")
            @show loss(valX[1], valY[1])
            @show adv_loss(valX[1], valY[1])
            @show accuracy(valX[1], valY[1])
            @show adv_accuracy(valX[1], valY[1])
            @show accuracy(trainX[1], trainY[1])
            @show adv_accuracy(trainX[1], trainY[1])
        end
    end
    evalcb = throttle(cb_fn, 10)

    # train
    # opt = Flux.Optimiser(Flux.ExpDecay(0.001, 0.5, 1000, 1e-4), ADAM(0.001))
    opt = ADAM(1e-3);
    @epochs 3 advtrain!(model, attack_fn, loss, Flux.params(ae), zip(trainX, trainY), opt, cb=evalcb)

    # DEBUG testing gc time
    # @time advtrain!(model, attack_fn, loss, Flux.params(ae), zip(trainX, trainY), opt, cb=evalcb)

    # DEBUG testing the trainning of AE using clean image CNN
    # @epochs 3 mytrain!(loss, Flux.params(ae), zip(trainX, trainY), opt, cb=evalcb)
end

function AdvAE()
    (trainX, trainY), (valX, valY), (testX, testY) = load_MNIST(batch_size=64);

    # cnn = get_LeNet5()
    cnn = get_Madry_model()

    # FIXME this is in adv.jl
    train_MNIST_model(cnn, trainX, trainY, valX, valY)

    # ae
    ae = CNN_AE()
    train_AE(ae, trainX, trainY, valX, valY)

    # test accuracy
    evaluate_AE(ae, cnn, trainX, trainY, testX, testY)

    # adv train of ae
    AdvAE_train(ae, cnn, attack_PGD_k(7), trainX, trainY, valX, valY)
    AdvAE_train(ae, cnn, attack_PGD_k(20), trainX, trainY, valX, valY)
    AdvAE_train(ae, cnn, attack_PGD_k(40), trainX, trainY, valX, valY)

    # evaluate attack
    evaluate_attack(Chain(ae, cnn), attack_FGSM, trainX, trainY, testX, testY)
    evaluate_attack(Chain(ae, cnn), attack_PGD_k(20), trainX, trainY, testX, testY)
    evaluate_attack(Chain(ae, cnn), attack_PGD_k(40), trainX, trainY, testX, testY)
    evaluate_attack(cnn, attack_PGD_k(40), trainX, trainY, testX, testY)
end

# Extend distributions slightly to have a numerically stable logpdf for `p` close to 1 or 0.
logpdf(b::Bernoulli, y::Bool) = y * log(b.p + eps(Float32)) + (1f0 - y) * log(1 - b.p + eps(Float32))

function vae()
    # Latent dimensionality, # hidden units.
    Dz, Dh = 5, 500
    # Components of recognition model / "encoder" MLP.
    A, μ, logσ = Dense(28^2, Dh, tanh), Dense(Dh, Dz), Dense(Dh, Dz)
    g(X) = (h = A(X); (μ(h), logσ(h)))
    z(μ, logσ) = μ + exp(logσ) * randn(Float32)
    # Generative model / "decoder" MLP.
    f = Chain(Dense(Dz, Dh, tanh), Dense(Dh, 28^2, σ))

    # KL-divergence between approximation posterior and N(0, 1) prior.
    kl_q_p(μ, logσ) = 0.5f0 * sum(exp.(2f0 .* logσ) + μ.^2 .- 1f0 .- (2 .* logσ))
    # logp(x|z) - conditional probability of data given latents.
    logp_x_z(x, z) = sum(logpdf.(Bernoulli.(f(z)), x))
    # Monte Carlo estimator of mean ELBO using M samples.
    L̄(X) = ((μ̂, logσ̂) = g(X); (logp_x_z(X, z.(μ̂, logσ̂)) - kl_q_p(μ̂, logσ̂)) * 1 // M)

    loss(X) = -L̄(X) + 0.01f0 * sum(x->sum(x.^2), Flux.params(f))

    # Sample from the learned model.
    modelsample() = rand.(Bernoulli.(f(z.(zeros(Dz), zeros(Dz)))))

    evalcb = throttle(() -> @show(-L̄(X[:, rand(1:N, M)])), 30)
    opt = ADAM()
    ps = Flux.params(A, μ, logσ, f)

    @progress for i = 1:20
        @info "Epoch $i"
        Flux.train!(loss, ps, zip(data), opt, cb=evalcb)
    end
end

function make_minibatch(X, Y, idxs)
    X_batch = Array{Float32}(undef, size(X[1])..., 1, length(idxs))
    for i in 1:length(idxs)
        X_batch[:, :, :, i] = Float32.(X[idxs[i]])
    end
    Y_batch = onehotbatch(Y[idxs], 0:9)
    return (X_batch, Y_batch)
end

"""
FIXME batch_size cannot equally divide
FIXME move to GPU on demand, otherwise GPU memory may exaust for large dataset
FIXME shuffle
FIXME minibatch
"""
function load_MNIST(; batch_size=32, val_split=0.1)
    imgs = Flux.Data.MNIST.images();
    labels = Flux.Data.MNIST.labels();

    # I don't want to reshape the data. I'll just add reshape layer
    # X = hcat(float.(reshape.(imgs, :))...);
    # Do I need to do the partition here?
    # N = length(imgs)
    # X = [float(hcat(vec.(imgs)...)) for imgs in partition(imgs, 1000)];

    function toXY(imgs, labels)
        let N = length(imgs),
            X = float.(imgs),
            Y = Flux.onehotbatch(labels, 0:9)

            # size(X)                     # (60000,)
            # size(X[1])                  # (28, 28)
            # size(Y)                     # (10, 60000)

            X, Y
        end
    end

    # batch_size = 128
    # mb_idxs = partition(1:length(train_imgs), batch_size)
    # train_set = [make_minibatch(train_imgs, train_labels, i) for i in mb_idxs];

    # Prepare test set as one giant minibatch:
    # test_imgs = MNIST.images(:test);
    # test_labels = MNIST.labels(:test);
    # test_set = make_minibatch(test_imgs, test_labels, 1:length(test_imgs));

    X, Y = toXY(Flux.Data.MNIST.images(:train), Flux.Data.MNIST.labels(:train))
    X2, Y2 = toXY(Flux.Data.MNIST.images(:test), Flux.Data.MNIST.labels(:test))

    N = length(X)
    # val_split = 0.1
    mid = round(Int, N * (1 - val_split))
    N2 = length(X2)

    trainX = gpu.([cat(X[i]..., dims = 4) for i in partition(1:mid, batch_size)]);
    trainY = gpu.([Y[:,i] for i in partition(1:mid, batch_size)]);

    valX = [cat(X[i]..., dims = 4) for i in partition(mid+1:N, batch_size)] .|> gpu;
    valY = [Y[:,i] for i in partition(mid+1:N, batch_size)] .|> gpu;

    testX = [cat(X2[i]..., dims = 4) for i in partition(1:N2, batch_size)] .|> gpu;
    testY = [Y2[:,i] for i in partition(1:N2, batch_size)] .|> gpu;

    return (trainX, trainY), (valX, valY), (testX, testY)

end

function test_load_MNIST()
    (trainX, trainY), (valX, valY), (testX, testY) = load_MNIST();
    sample_and_view(trainX[3])
    size(trainX)
    size(trainX[1])
    size(testX)
end

"""
    load_CIFAR10()

Return trainX, trainY, valX, valY, testX, testY

TODO I'm not going to move it onto GPU yet.

"""
function load_CIFAR10(; batch_size=32, val_split=0.1)
    getarray(X) = Float32.(permutedims(channelview(X), (2, 3, 1)))

    function toXY(imgs)
        let N = length(imgs),
            X = [getarray(imgs[i].img) for i in 1:N],
            # FIXME this has to be 1:10, not 0:9. Take care about the
            # consistency
            Y = onehotbatch([imgs[i].ground_truth.class for i in 1:N],1:10)

            X, Y
        end
    end

    X, Y = toXY(trainimgs(Metalhead.CIFAR10));

    N = length(X)
    # TODO use random index
    mid = round(Int, N * (1 - val_split))
    # println("N=$(N), mid=$(mid)")

    # There is no testimgs. the valimgs is the test.
    # imgs3 = testimgs(CIFAR10);
    X2, Y2 = toXY(valimgs(Metalhead.CIFAR10));
    N2 = length(X2);

    trainX = gpu.([cat(X[i]..., dims = 4) for i in partition(1:mid, batch_size)]);
    trainY = gpu.([Y[:,i] for i in partition(1:mid, batch_size)]);

    valX = gpu.([cat(X[i]..., dims = 4) for i in partition(mid+1:N, batch_size)]);
    valY = gpu.([Y[:,i] for i in partition(mid+1:N, batch_size)]);

    testX = gpu.([cat(X2[i]..., dims = 4) for i in partition(1:N2, batch_size)]);
    testY = gpu.([Y2[:,i] for i in partition(1:N2, batch_size)]);

    return (trainX, trainY), (valX, valY), (testX, testY)
end

function test_load_CIFAR10()
    data = load_CIFAR10();
    (trainX, trainY), (valX, valY), (testX, testY) = data;

    (trainX, trainY), (valX, valY), (testX, testY) = load_CIFAR10();
    (trainX, trainY), (valX, valY), (testX, testY) = load_CIFAR10(batch_size=50);
    size(trainX)
    size(trainY)
    size(testY)
end

call(f, xs...) = f(xs...)
runall(f) = f
runall(fs::AbstractVector) = () -> foreach(call, fs)
"""FIXME replace Tracker with Zygote
"""
function mytrain!(loss, ps, data, opt; cb = () -> ())
    # FIXME maybe this is used to reset the gradients?
    ps = Flux.Tracker.Params(ps)
    cb = runall(cb)
    @showprogress 0.1 "Training..." for d in data
        gs = Flux.Tracker.gradient(ps) do
            loss(d...)
        end
        Flux.Tracker.update!(opt, ps, gs)
        cb()
    end
end

function get_MNIST_FC_model()
    model = Chain(
        # (28,28,N)
        x -> reshape(x, :, size(x, 4)),
        # TODO regularizers?
        # TODO l2 weight decay?
        Dense(28^2, 32, relu),
        Dense(32, 10),
        softmax) |> gpu;
    return model
end

function get_MNIST_CNN_model()
    model = Chain(
        # First convolution, operating upon a 28x28 image
        Conv((3, 3), 1=>16, pad=(1,1), relu),
        # x -> maxpool(x, (2,2)),
        MaxPool((2,2)),

        # Second convolution, operating upon a 14x14 image
        Conv((3, 3), 16=>32, pad=(1,1), relu),
        # x -> maxpool(x, (2,2)),
        MaxPool((2,2)),

        # Third convolution, operating upon a 7x7 image
        Conv((3, 3), 32=>32, pad=(1,1), relu),
        # x -> maxpool(x, (2,2)),
        MaxPool((2,2)),

        # Reshape 3d tensor into a 2d one, at this point it should be (3, 3, 32, N)
        # which is where we get the 288 in the `Dense` layer below:
        x -> reshape(x, :, size(x, 4)),
        Dense(3*3*32, 10),

        # Finally, softmax to get nice probabilities
        softmax,
    ) |> gpu;
    return model
end

function test_tracker()
    (trainX, trainY), (valX, valY), (testX, testY) = load_MNIST();
    cnn = get_MNIST_CNN_model()

    size(trainX[1])
    # param(trainX[1]).grad[:,:,:,1]

    model = cnn
    x = trainX[1];
    y = trainY[1];
    # but this will clear it as well
    px = Flux.param(x);

    sum(px.grad)
    sum(px.tracker.grad)

    typeof(px.grad)
    typeof(px.data)

    loss(x, y) = crossentropy(model(x), y)

    grad = Flux.Tracker.extract_grad!(px)

    typeof(grad[1])
    grad.grad

    sign(grad.data)
    px.grad .= 0;

    theta = Flux.params(model);
    typeof(theta)
    typeof(px)

    newtheta = params(theta..., px)

    g = Flux.Tracker.gradient(() -> loss(px, y));
    g = Flux.Tracker.gradient(() -> loss(px, y), theta);
    g = Flux.Tracker.gradient(() -> loss(px, y), newtheta);

    for t in theta
        # @show typeof(t)
        @show sum(t.grad)
        @show sum(g[t])
    end

    sum(g[px])
end

function CNN3_AE()
    encoder = Chain(Conv((3,3), 1=>16, pad=(1,1), relu),
                    MaxPool((2,2)),
                    Conv((3,3), 16=>8, pad=(1,1), relu),
                    MaxPool((2,2)),
                    # FIMXE this is from 7 to 3, not good
                    Conv((3,3), 8=>8, pad=(1,1), relu),
                    MaxPool((2,2)))
    decoder = Chain(Conv((3,3), 8=>8, pad=(1,1), relu),
                    upsample,
                    Conv((3,3), 8=>8, pad=(1,1), relu),
                    upsample,
                    Conv((3,3), 8=>16, pad=(1,1), relu),
                    upsample,
                    Conv((3,3), 16=>1, pad=(1,1)),
                    x -> σ.(x))
    Chain(encoder, decoder) |> gpu
end

function test_MNIST_model(model)
    (trainX, trainY), (valX, valY), (testX, testY) = load_MNIST();

    model(trainX[1]);

    function alternative_loss(x, y)
        # We augment `x` a little bit here, adding in random noise
        x_aug = x .+ 0.1f0*gpu(randn(eltype(x), size(x)))
        y_hat = model(x_aug)
        return crossentropy(y_hat, y)
    end

    loss(x, y) = crossentropy(model(x), y)
    accuracy(x, y) = mean(onecold(model(x)) .== onecold(y))

    evalcb = throttle(() -> @show(loss(valX[1], valY[1])) , 5);
    opt = ADAM(0.001);

    # TODO early stopping
    # TODO learning rate decay

    @epochs 10 mytrain!(loss, Flux.params(model), zip(trainX, trainY), opt, cb=evalcb)

    # print out training details, e.g. accuracy
    @show accuracy(trainX[1], trainY[1])
    @show accuracy(valX[1], valY[1])
    # Test set accuracy
    @show accuracy(testX[1], testY[1])

    sample_and_view(trainX[1], trainY[1])
    sample_and_view(testX[1], testY[1])
end

function test_MNIST()
    fc = get_MNIST_FC_model()
    test_MNIST_model(fc)
    cnn = get_MNIST_CNN_model()
    test_MNIST_model(cnn)
end

"""
CIFAR raw CNN models. Training acc 0.43, val 0.55, testing 0.4
"""
function test_CIFAR_model(model)
    (trainX, trainY), (valX, valY), (testX, testY) = load_CIFAR10();
    model(trainX[1]);

    function loss(x, y)
        # We augment `x` a little bit here, adding in random noise
        x_aug = x .+ 0.1f0*gpu(randn(eltype(x), size(x)))
        y_hat = model(x_aug)
        return crossentropy(y_hat, y)
    end

    accuracy(x, y) = mean(onecold(model(x)) .== onecold(y))

    evalcb = throttle(() -> @show(loss(valX[1], valY[1])) , 5);
    opt = ADAM(0.001)

    # training
    @epochs 10 mytrain!(loss, Flux.params(model), zip(trainX, trainY), opt, cb=evalcb)

    # into testing mode
    Flux.testmode!(model)

    # test accuracy
    @show accuracy(trainX[1], trainY[1])
    @show accuracy(valX[1], valY[1])
    @show accuracy(testX[1], testY[1])

    # visualize the dataset
    sample_and_view(trainX[1], trainY[1])
    sample_and_view(testX[1], testY[1])

    # back into training node
    Flux.testmode!(model, false)
end

function test_CIFAR()
    accuracy_with_logits(model(testX[1]), testY[1])

    model = get_CIFAR_CNN_model()
    test_CIFAR_model(model)
    resnet20 = resnet(3) |> gpu;
    # 0.87, 0.74, 0.7
    resnet32 = resnet(5) |> gpu;
    resnet56 = resnet(9) |> gpu;
    resnet68 = resnet(11) |> gpu;
    test_CIFAR_model(resnet32)
end

function test_nested_meter()
    @showprogress 0.1 "Level 0 " for i in 1:10
        sleep(0.5)
        if i % 2 == 0
            @showprogress 0.1 " Level 1 " 1 for i2 in 1:10
                sleep(0.5)
                @showprogress 0.1 "  Level 2 " 2 for i3 in 1:10
                    sleep(0.5)
                end
            end
        end
        sleep(0.5)
    end
end

function test_logger()
    name="$(now())"
    TeeLogger(global_logger(),
              TBLogger("tensorboard_logs/exp-$name",
                       min_level=Logging.Info))
    TBLogger("tensorboard_logs/exp-$name", min_level=Logging.Info)

    # 7937 bytes
    @save "test.bson" weights=Tracker.data.(Flux.params(cpu(model)))
    # 37482 bytes, so I'm going to just save the model
    @save "test.bson" model
    with_logger(ConsoleLogger(stderr, Logging.Debug)) do
        @info "hello"
        @debug "world"
    end
end

"""
Several problems before going forward:

1. the augmentation should happen on-the-fly, so that each batch receives new
kinds of augmentation

2. I need a batched version of this for performance

"""
function test_augment()
    # the padding is easy, but there does not seem to be a random crop, I'll
    # need to write a simple random slicing algorithm
    #
    # In the meanwhile, there are more advanced augmentation techniques
    # using Augmentor
    pl = ElasticDistortion(6, scale=0.3, border=true) |>
        Rotate([10, -5, -3, 0, 3, 5, 10]) |>
        FlipX(0.5) |>
        ShearX(-10:10) * ShearY(-10:10) |>
        CropSize(32, 32) |>
        Zoom(0.9:0.1:1.2)
    augment(img, pl)
end

function evaluate(model, ds; attack_fn=(m,x,y)->x)
    xx,yy = next_batch!(ds) |> gpu

    @info "Sampling BEFORE images .."
    sample_and_view(xx, yy, model)
    @info "Sampling AFTER images .."
    adv = attack_fn(model, xx, yy)
    sample_and_view(adv, yy, model)

    @info "Testing multiple batches .."
    acc_metric = MeanMetric()
    @showprogress 0.1 "Testing..." for step in 1:10
        x,y = next_batch!(ds) |> gpu
        adv = attack_fn(model, x, y)
        acc = accuracy_with_logits(model(adv), y)
        add!(acc_metric, acc)
    end
    @show get!(acc_metric)
end

function test()
    Random.seed!(1234);

    ds, test_ds = load_MNIST_ds(batch_size=50);
    x, y = next_batch!(ds) |> gpu;

    model = get_Madry_model()[1:end-1]
    # model = get_LeNet5()[1:end-1]

    # Warm up the model
    model(x)

    # FIXME would this be 0.001 * 0.001?
    # FIXME decay on pleau
    # FIXME print out information when decayed
    # opt = Flux.Optimiser(Flux.ExpDecay(0.001, 0.5, 1000, 1e-4), ADAM(0.001))
    opt = ADAM(1e-4);

    logger = create_logger()
    with_logger(logger) do
        train!(model, opt, ds, print_steps=20)
    end

    @info "Adv training .."
    # custom_train!(model, opt, train_ds)

    with_logger(logger) do
        @epochs 2 advtrain!(model, opt, attack_PGD_k(40), ds, print_steps=20)
    end

    @info "Evaluating clean accuracy .."
    evaluate(model, test_ds)
    evaluate(model, test_ds, attack_fn=attack_FGSM)
    evaluate(model, test_ds, attack_fn=attack_PGD_k(40))
end

# FIXME model cannot be easily saved and loaded, weights can, but
# needs to get rid of CuArrays and TrackedArrays
#
@time @save model_file weights=Tracker.data.(Flux.params(cpu(model))) from_steps=step

max_stars_repo_name: lihebi/AdvAE

max_issues_repo_name: lihebi/AdvAE

max_forks_repo_name: lihebi/AdvAE


text:
module example

export func

"""
    func(x)

Returns double the number `x` plus `1`.
"""
func(x) = 2x + 1

end


text:
using ..Fed: STATS

"""
    VanillaPayloadSerde() <: PayloadSerde

Basic serialization where no compression is applied. This `PayloadSerde` is
only a wrapper around the data serialization in order to match the interface to
implement.
"""
struct VanillaPayloadSerde <: PayloadSerde end

"""
    serialize_payload(::VanillaPayloadSerde, weights)

Serialize `weights` with the `VanillaPayloadSerde` where only basic
serialization is applied.
"""
function serialize_payload(::VanillaPayloadSerde, weights::Vector{Float32})::Vector{UInt8}
    STATS.base.req_data = weights
    return pack(weights)
end

"""
    deserialize_payload(::VanillaPayloadSerde, data, from)

Deserialize `data` with the `VanillaPayloadSerde` where only basic
deserialization is applied.
"""
function deserialize_payload(::VanillaPayloadSerde, data::Vector{UInt8}, from::String)::Vector{Float32}
    weights = unpack(data)
    push!(STATS.base.res_data, weights)
    return weights
end


text:
mutable struct NatureHook10 <: AbstractHook
    step::Int
    max_steps::Int
    food_counts::Vector{Float32}
    food_balance::Vector{Float32}
    act_counts::Vector{Int}
    act_probs::Vector{Float32}
    player_acts::Vector{Vector{Int}}
    total_rewards::Vector{Float32}
    trial_id::String
end
function NatureHook10(env::NatureEnv, trial_id::String, max_steps::Int)
    NatureHook(0,
        max_steps,
        zeros(Float32, env.food_types),
        [],
        zeros(length(action_space(env, 1))),
        Vector{Float32}(),
        [zeros(length(action_space(env, 1))) for i in 1:length(env.players)],
        zeros(length(env.players)),
        trial_id,
       )
end
NatureHook = NatureHook10

function (hook::NatureHook)(::PreActStage, policy, env, action)
    for p in 1:length(env.players)
        hook.act_counts[action[p].action] += 1
        push!(hook.act_probs, exp(action[p].meta.action_log_prob[1]))
        hook.player_acts[p][action[p].action] += 1
    end
end

function (hook::NatureHook)(::PostActStage, policy, env)
    for p in 1:length(env.players)
        hook.total_rewards[p] += reward(env, p)
        hook.food_counts = hook.food_counts .+ env.players[p].food_counts
        push!(hook.food_balance, max(env.players[p].food_counts...) / min(env.players[p].food_counts...))
    end
    hook.step += 1
    steps_per_checkpoint = 1_000
    if hook.step % steps_per_checkpoint == 0
        log_step = if hasproperty(global_logger(), :global_step)
            global_logger().global_step
        else -1 end
        serialize("checkpoints/$(hook.trial_id)/$(Int(hook.step / steps_per_checkpoint)).jls",
                  Dict(:policy => policy,
                       :tb_step => log_step,
                       :hook_step => hook.step))
    end

end

function (hook::NatureHook)(::PostEpisodeStage, policy, env)

    @info "episode"     len=env.step                         log_step_increment=0
    @info "total_food"  total_food=Tuple(hook.food_counts)   log_step_increment=0
    @info "actions"     act=Tuple(hook.act_counts)           log_step_increment=0
    @info "rewards"     reward=Tuple(hook.total_rewards)     log_step_increment=0
    @info "rewards"     total_reward=sum(hook.total_rewards) log_step_increment=0
    @info "act_prob"    act_prob=mean(hook.act_probs)        log_step_increment=0
    @info "player_acts" player_acts=Tuple([Tuple(hook.player_acts[p]) for p in 1:length(env.players)]) log_step_increment=0
    @info "exchanges"   ex=Tuple(env.exchanges)              log_step_increment=0
    @info "food_balance" food_balance=mean(hook.food_balance)

    hook.act_counts    = zeros(length(action_space(env, 1)))
    hook.total_rewards = zeros(length(env.players))
    hook.food_counts   = zeros(env.food_types)
    hook.act_probs     = Vector{Float32}()
    hook.food_balance  = []

    hook.player_acts = [zeros(length(action_space(env, 1))) for _ in env.players]
end


text:
#=
Given an array of integers and a number k, where 1 <= k <= length of the array, compute the maximum values of each subarray of length k.

For example, given array = [10, 5, 2, 7, 8, 7] and k = 3, we should get: [10, 7, 8, 8], since:

    10 = max(10, 5, 2)
    7 = max(5, 2, 7)
    8 = max(2, 7, 8)
    8 = max(7, 8, 7)

Do this in O(n) time and O(k) space. You can modify the input array in-place and you do not need to store the results. You can simply print them out as you compute them.
=#

function max_in_subarrays(X::Array{Int64,1}, k::Int64; store_results = false)
    if store_results
        maximums = []
    end
    for index in 1:length(X)-k+1
        subarray = X[1:k]
        popfirst!(X)
        if store_results
            append!(maximums, maximum(subarray))
        else
            println(maximum(subarray))
        end
    end
    if store_results
        return maximums
    end
end

max_stars_repo_name: DominiqueCaron/daily-coding-problem

max_issues_repo_name: DominiqueCaron/daily-coding-problem

max_forks_repo_name: DominiqueCaron/daily-coding-problem


text:
using NNlib

# Activation functions
@cufunc softplus(x::Real) = ifelse(x > 0, x + log1p(exp(-x)), log1p(exp(x)))

@cufunc logσ(x::Real) = -softplus(-x)

@cufunc function gelu(x::Real)
    p = oftype(x / 1, π)
    λ = oftype(x / 1, √(2 / p))
    α = oftype(x / 1, 0.044715)
    h = oftype(x / 1, 0.5)
    h * x * (one(x) + tanh(λ * (x + α * x^3)))
end

@cufunc lisht(x::Real) = x * tanh(x)

@cufunc logcosh(x::Real) = x + softplus(-2x) - log(oftype(x, 2))

@cufunc mish(x::Real) = x * tanh(softplus(x))

@cufunc tanhshrink(x::Real) = x - tanh(x)

# Batched matrix multiplication
# 1st argument is produced by NNlib.storage_type(A)
NNlib._batched_gemm!(::Type{<:CuArray}, transA::Char, transB::Char, α::Number, A, B, β::Number, C) =
     CUBLAS.gemm_strided_batched!(transA, transB, α, A, B, β, C)

Base.unsafe_convert(::Type{CuPtr{T}}, A::NNlib.BatchedAdjOrTrans{T}) where {T} =
    Base.unsafe_convert(CuPtr{T}, parent(A))

#
# 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.

function upsample_bilinear_whcn_kernel!(n_elem, rheight, rwidth, x, y)
    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
        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
        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, rheight, rwidth, Δ, dx)
    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
        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
        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]
                @atomic dx[ow0, oh0, c, n] += h0lambda * w0lambda * val
                @atomic dx[ow1, oh0, c, n] += h0lambda * w1lambda * val
                @atomic dx[ow0, oh1, c, n] += h1lambda * w0lambda * val
                @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}) where T
    w,h,c,n = size(x)
    out_w, out_h = (size(y,1), size(y,2))

    out_size = out_h*out_w
    rheight = T((h-1)/(out_h-1))
    rwidth  = T((w-1)/(out_w-1))

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

function NNlib.∇upsample_bilinear_whcn!(dx::CuArray{T,4}, Δ::CuArray{T,4}) where T
    w,h,c,n = Base.size(Δ)
    out_w, out_h = (size(dx, 1), size(dx, 2))
    in_size = h*w
    rheight = T((out_h-1)/(h-1)) # reversed compared to forward pass
    rwidth  = T((out_w-1)/(w-1))

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


text:
abstract PSMMALAState <: MCSamplerState

abstract MuvPSMMALAState <: PSMMALAState

abstract PSMMALA <: LMCSampler

mala_only_update!(sstate::MuvPSMMALAState, pstate::ParameterState{Continuous, Multivariate}, i::Integer, tot::Integer) =
  sstate.presentupdatetensor = false

smmala_only_update!(sstate::MuvPSMMALAState, pstate::ParameterState{Continuous, Multivariate}, i::Integer, tot::Integer) =
  sstate.presentupdatetensor = true

mod_update!(
  sstate::MuvPSMMALAState,
  pstate::ParameterState{Continuous, Multivariate},
  i::Integer,
  tot::Integer,
  n::Integer=10
) =
  sstate.presentupdatetensor = mod(i, n) == 0 ? true : false

cos_update!(
  sstate::MuvPSMMALAState,
  pstate::ParameterState{Continuous, Multivariate},
  i::Integer,
  tot::Integer,
  a::Real=10.,
  b::Real=0.2,
  c::Real=0.2
) =
  sstate.presentupdatetensor = rand(Bernoulli(b*cos(a*i*pi/tot)+c))

rand_update!(
  sstate::MuvPSMMALAState,
  pstate::ParameterState{Continuous, Multivariate},
  i::Integer,
  tot::Integer,
  p::Real=0.5
) =
  sstate.presentupdatetensor = rand(Bernoulli(p))

rand_exp_decay_update!(
  sstate::MuvPSMMALAState,
  pstate::ParameterState{Continuous, Multivariate},
  i::Integer,
  tot::Integer,
  a::Real=10.,
  b::Real=0.
) =
  sstate.presentupdatetensor = rand(Bernoulli(exp_decay(i, tot, a, b)))

rand_pow_decay_update!(
  sstate::MuvPSMMALAState,
  pstate::ParameterState{Continuous, Multivariate},
  i::Integer,
  tot::Integer,
  a::Real=1e-3,
  b::Real=0.
) =
  sstate.presentupdatetensor = rand(Bernoulli(pow_decay(i, tot, a, b)))

rand_lin_decay_update!(
  sstate::MuvPSMMALAState,
  pstate::ParameterState{Continuous, Multivariate},
  i::Integer,
  tot::Integer,
  a::Real=1e-2,
  b::Real=0.
) =
  sstate.presentupdatetensor = rand(Bernoulli(lin_decay(i, tot, a, b)))

rand_quad_decay_update!(
  sstate::MuvPSMMALAState,
  pstate::ParameterState{Continuous, Multivariate},
  i::Integer,
  tot::Integer,
  a::Real=1e-5,
  b::Real=0.
) =
  sstate.presentupdatetensor = rand(Bernoulli(quad_decay(i, tot, a, b)))


text:
using ReachabilityAnalysis
using ReachabilityAnalysis: solve
using JuMP, Ipopt
using Plots
using LaTeXStrings

# IVP parameters
A = [0.0 1; -4π 0]
X = Universe(2)

# discretization parameter
δ = 0.005

# varying initial states
X0s = [
       (Singleton([0.0, 10.0]), "□(x₀, 0)"),
       (BallInf([0.0, 10.0], 0.01), "□(x₀, 0.01)"),
       (BallInf([0.0, 10.0], 0.1), "□(x₀, 0.1)"),
       (BallInf([0.0, 10.0], 1.0), "□(x₀, 1)")
      ]

# methods
algs = [
        (FirstOrderddt(), "d/dt"),
        (FirstOrderZonotope(), "Zonotope"),
        (CorrectionHull(order=4), "Correction hull"),
        (FirstOrder(), "First order"),
        (ForwardBackward(solver=Ipopt.Optimizer()), "Forward/backward"),
        (Forward(), "Forward")
       ]

# results
X02Ω0s = []
for (X0, name_X0) in X0s
    local Ω0s = []
    local P = @ivp(x' = A*x, x ∈ X, x(0) ∈ X0)
    for (alg, name_alg) in algs
        local Pd = discretize(P, δ, alg)
        local Ω0 = Pd.x0
        push!(Ω0s, Ω0)
    end
    push!(X02Ω0s, Ω0s)
end

# plot sets
colors = [:orange, :green, :blue, :purple, :red, :black]
plots = []
for (i, (X0, name_X0)) in enumerate(X0s)
    local fig = plot(leg=:bottomleft, xlab="X₀ = $name_X0")
    for (j, (alg, name_alg)) in enumerate(algs)
        local Ω0 = X02Ω0s[i][j]
        if name_alg == "Forward/backward"
            Ω0 = overapproximate(Ω0, PolarDirections(30))  # workaround because σ is not available
        end
        plot!(Ω0, alpha=0, color=colors[j], linecolor=colors[j], linealpha=1, lw=2, lab="Ω₀ $(name_alg)")
    end

    # reachable states at time t
    R(t) = linear_map(exp(A * t), X0)

    [plot!(R(δi), color=:gray, alpha=0.05, linealpha=0, lab="") for δi in range(0.0; stop=δ, length=20)]

    savefig("experiment_X0_$i.pdf")
    push!(plots, fig)
end

return plots


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct AzureReachabilityReportItem <: SwaggerModel
    provider::Any # spec type: Union{ Nothing, String } # spec name: provider
    azureLocation::Any # spec type: Union{ Nothing, String } # spec name: azureLocation
    latencies::Any # spec type: Union{ Nothing, Vector{AzureReachabilityReportLatencyInfo} } # spec name: latencies

    function AzureReachabilityReportItem(;provider=nothing, azureLocation=nothing, latencies=nothing)
        o = new()
        validate_property(AzureReachabilityReportItem, Symbol("provider"), provider)
        setfield!(o, Symbol("provider"), provider)
        validate_property(AzureReachabilityReportItem, Symbol("azureLocation"), azureLocation)
        setfield!(o, Symbol("azureLocation"), azureLocation)
        validate_property(AzureReachabilityReportItem, Symbol("latencies"), latencies)
        setfield!(o, Symbol("latencies"), latencies)
        o
    end
end # type AzureReachabilityReportItem

const _property_map_AzureReachabilityReportItem = Dict{Symbol,Symbol}(Symbol("provider")=>Symbol("provider"), Symbol("azureLocation")=>Symbol("azureLocation"), Symbol("latencies")=>Symbol("latencies"))
const _property_types_AzureReachabilityReportItem = Dict{Symbol,String}(Symbol("provider")=>"String", Symbol("azureLocation")=>"String", Symbol("latencies")=>"Vector{AzureReachabilityReportLatencyInfo}")
Base.propertynames(::Type{ AzureReachabilityReportItem }) = collect(keys(_property_map_AzureReachabilityReportItem))
Swagger.property_type(::Type{ AzureReachabilityReportItem }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_AzureReachabilityReportItem[name]))}
Swagger.field_name(::Type{ AzureReachabilityReportItem }, property_name::Symbol) =  _property_map_AzureReachabilityReportItem[property_name]

function check_required(o::AzureReachabilityReportItem)
    true
end

function validate_property(::Type{ AzureReachabilityReportItem }, name::Symbol, val)
end


text:
using Test

include("test_one_dimensional_time_series.jl")
# include("test_forward_map.jl")


text:
using SolverLogging
using Test
using Crayons
using SolverLogging: EntrySpec
# using BenchmarkTools # using Formatting
## Formatting defaults
lg = SolverLogging.Logger()
def = SolverLogging._default_formats()
@test !haskey(def, Float64)

# Add a new entry
@test SolverLogging.default_format(lg, Float64) == def[AbstractFloat]
@test haskey(lg.defaults, Float64)

# Set default
@test !haskey(lg.defaults, Int)
SolverLogging.set_default_format(lg, Int, "%3d")
@test haskey(lg.defaults, Int)
@test SolverLogging.default_format(lg, Int) !==
    SolverLogging.default_format(lg, Int32)
@test haskey(lg.defaults, Int32)

# Insert fields
@test length(lg.data) == 0
SolverLogging.setentry(lg, "alpha", Float64)
@test length(lg.data) == 1
@test lg.fmt["alpha"] == EntrySpec(Float64, def[AbstractFloat],1, 1, SolverLogging.DEFAULT_WIDTH)
@test lg.idx[1] == 1
@test SolverLogging.getidx(lg, "alpha") == 1
lg.data[1] = "alpha"

# Add an integer entry
SolverLogging.setentry(lg, "iter", Int, fmt="%10d")
@test lg.fmt["iter"] == EntrySpec(Int,"%10d",2,1,11)
@test length(lg.data) == 2
@test SolverLogging.getidx(lg, "iter") == 2
lg.data[2] = "iter"

# Insert twice
SolverLogging.setentry(lg, "iter", Int, fmt="%3d", index=2)
@test length(lg.data) == 2
@test lg.data[2] == "iter"
@test length(keys(lg.fmt)) == 2
@test lg.data == ["alpha","iter"]
@test SolverLogging.getidx(lg, "iter") == 2

# New string field
SolverLogging.setentry(lg, "info", String, index=1, width=20)
@test lg.fmt["info"] == EntrySpec(String,"%s",3,1,20)
@test SolverLogging.getidx(lg, "info") == 1
@test SolverLogging.getidx(lg, "alpha") == 2
@test SolverLogging.getidx(lg, "iter") == 3
lg.data[1] = "info"
@test lg.data == ["info", "alpha","iter"]

# Move an existing field (1 to 3)
SolverLogging.setentry(lg, "info", String, index=-1)
@test lg.fmt["info"] == EntrySpec(String,"%s",3,1,20)
@test lg.data == ["alpha","iter","info"]
@test lg.idx == [1,2,3]

# Change the width (not specifying type)
@test lg.fmt["alpha"].width == SolverLogging.DEFAULT_WIDTH
SolverLogging.setentry(lg, "alpha", width=12)
@test lg.fmt["alpha"].width == 12

# Change the ccrayon and test
SolverLogging.setentry(lg, "alpha", ccrayon=ConditionalCrayon(1, 10))
entry = lg.fmt["alpha"]
@test entry.ccrayon(5) == Crayon(reset=true)
@test entry.ccrayon(0) == crayon"green"
@test entry.ccrayon(11) == crayon"red"

# Make sure it doesn't erase the previous crayon if another field is modified
SolverLogging.setentry(lg, "alpha", width=9)
entry = lg.fmt["alpha"]
@test entry.ccrayon(5) == Crayon(reset=true)
@test entry.ccrayon(0) == crayon"green"
@test entry.ccrayon(11) == crayon"red"

# # Try adding new field without the type
# @test_throws AssertionError SolverLogging.setentry(lg, "ϕ")

# Add to 2nd to last field
SolverLogging.setentry(lg, "ϕ", Int32, index=-2)
@test lg.fmt["ϕ"] == EntrySpec(Int32,def[Integer],4,1, 10)
lg.data[3] = "phi"
@test lg.data == ["alpha","iter","phi","info"]
@test lg.idx == [1,2,4,3]

# move and use new format (3 to 2)
SolverLogging.setentry(lg, "ϕ", Int32, index=-3, fmt="%5d")
@test lg.fmt["ϕ"] == EntrySpec(Int32,"%5d", 4, 1, 10)
@test lg.data == ["alpha","phi","iter","info"]
@test lg.idx == [1,3,4,2]

# Change level
SolverLogging.setentry(lg, "ϕ", Int32, level=2)
@test lg.fmt["ϕ"] == EntrySpec(Int32,"%5d", 4, 2, 10)

# Change Width
SolverLogging.setentry(lg, "ϕ", Int32, width=12)
@test lg.fmt["ϕ"] == EntrySpec(Int32,"%5d", 4, 2, 12)

# Test clear
SolverLogging.clear!(lg)
for i in (1,3,4)
    @test all(isspace, lg.data[1])
end
@test isempty(lg.data[2])

#############################################
## Log values
#############################################
SolverLogging._log!(lg, "iter", 1)
@test lg.data[3] == "  1        "
@test parse(Int,lg.data[3]) == 1
SolverLogging._log!(lg, "iter", 200)
@test parse(Int,lg.data[3]) == 200

SolverLogging._log!(lg, "info", "hi there")
@test lg.data[4] == rpad("hi there", 20)

SolverLogging.setentry(lg, "alpha", Float64, width=6)
lg.opts.autosize = false
@test_logs (:warn,) SolverLogging._log!(lg, "alpha", 1.234567e-3)
@test lg.data[1] == "1.23e-03"

SolverLogging.setentry(lg, "alpha", Float64, width=10)
SolverLogging._log!(lg, "alpha", 1e-3)
@test lg.data[1] == "1.00e-03  "
@test lg.crayons[1] == crayon"green"

# Move alpha and make sure the crayon moves too
SolverLogging.setentry(lg, "alpha", index=2)
@test lg.data[2] == "1.00e-03  "
@test lg.crayons[2] == crayon"green"

SolverLogging.setentry(lg, "alpha", index=1)

# Test append operation
lg.opts.autosize = false
SolverLogging._log!(lg, "info", "hi there")
info = SolverLogging._getdata(lg, "info")
length(info) == 20
SolverLogging._log!(lg, "info", "Something", :append)
newinfo = SolverLogging._getdata(lg, "info")
@test newinfo == rpad("hi there. Something", 20)
@test_logs (:warn,) SolverLogging._log!(lg, "info", "new", :append)
info = SolverLogging._getdata(lg, "info")
@test length(info) > 20

lg.opts.autosize = true
SolverLogging._log!(lg, "info", "hi there")
info = SolverLogging._getdata(lg, "info")
length(info) == 20
SolverLogging._log!(lg, "info", "Something", :append)
newinfo = SolverLogging._getdata(lg, "info")
@test newinfo == rpad("hi there. Something", 20)
@test_nowarn SolverLogging._log!(lg, "info", "new", :append)
info = SolverLogging._getdata(lg, "info")
@test length(info) > 20
@test Int(lg.fmt["info"].width) == length(info)

@test_logs (:warn,) SolverLogging._log!(lg, "iter", 1, :append)
SolverLogging._log!(lg, "iter", "new", :append)
newinfo = SolverLogging._getdata(lg, "info")

# Test add operation
SolverLogging._log!(lg, "iter", 11)
iter = parse(Int,SolverLogging._getdata(lg, "iter"))
SolverLogging._log!(lg, "iter", 2, :add)
newiter = parse(Int,SolverLogging._getdata(lg, "iter"))
@test newiter == iter + 2

SolverLogging._log!(lg, "info", "hi there")
@test_logs (:warn,r"Cannot add*") SolverLogging._log!(lg, "info", "a", :add)

## Test printing and verbosity
setentry(lg, "info", width=20)
SolverLogging._log!(lg, "info", "")
SolverLogging.setlevel!(lg, 1)
@test lg.data[2] == ""
SolverLogging._log!(lg, "ϕ", 3)
@test lg.data[2] == ""
@test !occursin("ϕ", SolverLogging.formheader(lg))
@test length(SolverLogging.formrow(lg)) == 41
setentry(lg, "info", index=2)

setentry(lg, "info", width=20)
SolverLogging._log!(lg, "info", "hi there")
SolverLogging.printheader(lg)
for i = 1:10
    SolverLogging._log!(lg, "alpha", 2i-5)
    SolverLogging.printrow(lg)
end

# Try expanding column
setentry(lg, "info", width=20)
SolverLogging._log!(lg, "info", "hi there")
lg.opts.freq = 10
SolverLogging.resetcount!(lg)
for i = 1:10
    SolverLogging._log!(lg, "alpha", 2i-5)
    SolverLogging._log!(lg, "info", "$i", :append)
    SolverLogging.printlog(lg)
end

@test SolverLogging.setlevel!(lg, 2) == 1
SolverLogging._log!(lg, "ϕ", 3)
@test lg.data[3] == rpad("    3", 12)
@test occursin("ϕ", SolverLogging.formheader(lg))
@test length(SolverLogging.formrow(lg)) == 85

begin
    SolverLogging.printheader(lg)
    for i = 1:3
        SolverLogging._log!(lg, "iter", i)
        SolverLogging.printrow(lg)
    end
end

# Try append operation for strings
setentry(lg, "info", width=25)
SolverLogging._log!(lg, "info", "more info", :append)
@test occursin("more info", lg.data[SolverLogging.getidx(lg, "info")]) &&
    occursin("hi there.", lg.data[SolverLogging.getidx(lg, "info")])
printlog(lg)

# Print a lower verbosity level and make sure there
# aren't any extra entries
@test SolverLogging.setlevel!(lg, 1) == 2
begin
    SolverLogging.printheader(lg)
    for i = 1:3
        SolverLogging._log!(lg, "alpha", 10(i-1)-5)
        SolverLogging._log!(lg, "iter", i)
        SolverLogging.printrow(lg)
    end
end
SolverLogging._log!(lg, "ϕ", 21.2)
header = SolverLogging.formheader(lg)
@test !occursin("ϕ", header)
row = SolverLogging.formrow(lg)
@test !occursin("21.5", row)


text:
##############################################################################
##
## AbstractDataArray is a type of AbstractArray that can handle NA's
##
##############################################################################

abstract AbstractDataArray{T, N} <: AbstractArray{T, N}

##############################################################################
##
## DataArray type definition
##
##############################################################################

type DataArray{T, N} <: AbstractDataArray{T, N}
    data::Array{T, N}
    na::BitArray{N}

    # Sanity check that new data values and missingness metadata match
    function DataArray(d::Array{T, N}, m::BitArray{N})
        if size(d) != size(m)
            msg = "Data and missingness arrays must be the same size"
            throw(ArgumentError(msg))
        end
        new(d, m)
    end
end

##############################################################################
##
## DataVector and DataMatrix are typealiases for 1D and 2D DataArray's
##
##############################################################################

typealias AbstractDataVector{T} AbstractDataArray{T, 1}
typealias AbstractDataMatrix{T} AbstractDataArray{T, 2}
typealias DataVector{T} DataArray{T, 1}
typealias DataMatrix{T} DataArray{T, 2}

##############################################################################
##
## DataArray constructors
##
##############################################################################

# Need to redefine inner constructor as outer constuctor
DataArray{T, N}(d::Array{T, N}, m::BitArray{N}) = DataArray{T, N}(d, m)

# Convert an Array into a DataArray w/o NA's
DataArray(d::Array) = DataArray(d, falses(size(d)))

# A no-op constructor
DataArray(d::DataArray) = d

# Convert Array{Bool}'s to BitArray's to save space
DataArray(d::Array, m::Array{Bool}) = DataArray(d, bitpack(m))

# Convert a BitArray into a DataArray w/ NA's
DataArray(d::BitArray, m::BitArray) = DataArray(convert(Array{Bool}, d), m)

# Convert a BitArray into a DataArray w/o NA's
DataArray(d::BitArray) = DataArray(convert(Array{Bool}, d), falses(size(d)))

# Convert a Ranges object into a DataVector
DataArray(r::Ranges) = DataArray([r], falses(length(r)))

# Construct an all-NA DataArray of a specific type
DataArray(t::Type, dims::Integer...) = DataArray(Array(t, dims...),
                                                 trues(dims...))

# Wrap a scalar in a DataArray w/ repetition
function DataArray(val::Any, dims::Integer...)
    vals = Array(typeof(val), dims...)
    for i in 1:length(vals)
        vals[i] = val
    end
    DataArray(vals, falses(dims...))
end

# Wrap a scalar in a DataArray w/o repetition
DataArray(val::Any) = DataArray([val], falses(1))

##############################################################################
##
## Initialized constructors
##
##############################################################################

# Initialized constructors with 0's, 1's
for (f, basef) in ((:datazeros, :zeros), (:dataones, :ones))
    @eval begin
        ($f)(dims::Int...) = DataArray(($basef)(dims...), falses(dims...))
        ($f)(t::Type, dims::Int...) = DataArray(($basef)(t, dims...),
                                                falses(dims...))
    end
end

# Initialized constructors with false's or true's
for (f, basef) in ((:datafalses, :falses), (:datatrues, :trues))
    @eval begin
        ($f)(dims::Int...) = DataArray(($basef)(dims...), falses(dims...))
    end
end

##############################################################################
##
## Copying
##
##############################################################################

copy(d::DataArray) = DataArray(copy(d.data), copy(d.na))
deepcopy(d::DataArray) = DataArray(deepcopy(d.data), deepcopy(d.na))
function copy_to(dest::DataArray, src::Any)
    i = 1
    for x in src
        dest[i] = x
        i += 1
    end
    return dest
end

##############################################################################
##
## similar()
##
##############################################################################

similar(d::DataArray) = d

function similar{T}(d::DataArray{T}, dims::Int...)
    DataArray(Array(T, dims...), trues(dims...))
end

function similar{T}(d::DataArray{T}, dims::Dims)
    DataArray(Array(T, dims), trues(dims))
end

##############################################################################
##
## Size information
##
##############################################################################

size(d::DataArray) = size(d.data)
ndims(d::DataArray) = ndims(d.data)
length(d::DataArray) = length(d.data)
endof(d::DataArray) = endof(d.data)
eltype{T, N}(d::DataArray{T, N}) = T

##############################################################################
##
## Generic Strategies for dealing with NA's
##
## Editing Functions:
##
## * failNA: Operations should die on the presence of NA's.
## * removeNA: What was once called FILTER.
## * replaceNA: What was once called REPLACE.
##
## Iterator Functions:
##
## * each_failNA: Operations should die on the presence of NA's.
## * each_removeNA: What was once called FILTER.
## * each_replaceNA: What was once called REPLACE.
##
## v = failNA(dv)
##
## for v in each_failNA(dv)
##     do_something_with_value(v)
## end
##
##############################################################################

function failNA(da::DataArray)
    if anyna(da)
        throw(NAException())
    else
        return copy(da.data)
    end
end

# Can do strange things on DataArray of rank > 1
function removeNA(da::DataArray)
    return copy(da.data[!da.na])
end

function replaceNA(da::DataArray, replacement_val::Any)
    res = copy(da.data)
    for i in 1:length(da)
        if da.na[i]
            res[i] = replacement_val
        end
    end
    return res
end

# TODO: Re-implement these methods for PooledDataArray's
function failNA{T}(da::AbstractDataArray{T})
    if anyna(da)
        throw(NAException())
    else
        res = Array(T, size(da))
        for i in 1:length(da)
            res[i] = da[i]
        end
        return res
    end
end

# TODO: Figure out how to make this work for Array's
function removeNA{T}(da::AbstractDataVector{T})
    n = length(da)
    res = Array(T, n)
    total = 0
    for i in 1:n
        if !isna(da[i])
            total += 1
            res[total] = convert(T, da[i])
        end
    end
    return res[1:total]
end

function replaceNA{S, T}(da::AbstractDataArray{S}, replacement_val::T)
    res = Array(S, size(da))
    for i in 1:length(da)
        if isna(da[i])
            res[i] = replacement_val
        else
            res[i] = da[i]
        end
    end
    return res
end

##
## NA-aware iterators
##

type EachFailNA{T}
    da::AbstractDataArray{T}
end
each_failNA{T}(da::AbstractDataArray{T}) = EachFailNA(da)
start(itr::EachFailNA) = 1
function done(itr::EachFailNA, ind::Integer)
    return ind > length(itr.da)
end
function next(itr::EachFailNA, ind::Integer)
    if isna(itr.da[ind])
        throw(NAException())
    else
        (itr.da[ind], ind + 1)
    end
end

type EachRemoveNA{T}
    da::AbstractDataArray{T}
end
each_removeNA{T}(da::AbstractDataArray{T}) = EachRemoveNA(da)
start(itr::EachRemoveNA) = 1
function done(itr::EachRemoveNA, ind::Integer)
    return ind > length(itr.da)
end
function next(itr::EachRemoveNA, ind::Integer)
    while ind <= length(itr.da) && isna(itr.da[ind])
        ind += 1
    end
    (itr.da[ind], ind + 1)
end

type EachReplaceNA{S, T}
    da::AbstractDataArray{S}
    replacement_val::T
end
function each_replaceNA(da::AbstractDataArray, val::Any)
    EachReplaceNA(da, convert(eltype(da), val))
end
start(itr::EachReplaceNA) = 1
function done(itr::EachReplaceNA, ind::Integer)
    return ind > length(itr.da)
end
function next(itr::EachReplaceNA, ind::Integer)
    if isna(itr.da[ind])
        (itr.replacement_val, ind + 1)
    else
        (itr.da[ind], ind + 1)
    end
end

##############################################################################
##
## ref()
##
##############################################################################

# v[dv]
function ref(x::Vector, inds::AbstractDataVector{Bool})
    return x[find(replaceNA(inds, false))]
end
function ref{S, T}(x::Vector{S}, inds::AbstractDataVector{T})
    return x[removeNA(inds)]
end

# d[SingleItemIndex]
function ref(d::DataArray, i::Real)
	if d.na[i]
		return NA
	else
		return d.data[i]
	end
end

# d[MultiItemIndex]
# TODO: Return SubDataArray
function ref(d::DataArray, inds::AbstractDataVector{Bool})
    inds = find(replaceNA(inds, false))
    return DataArray(d.data[inds], d.na[inds])
end
# TODO: Return SubDataArray
function ref(d::DataArray, inds::AbstractDataVector)
    inds = removeNA(inds)
    return DataArray(d.data[inds], d.na[inds])
end
# TODO: Return SubDataArray
# TODO: Make inds::AbstractVector
function ref(d::DataArray, inds::Union(Vector, BitVector, Ranges))
    return DataArray(d.data[inds], d.na[inds])
end

# dm[SingleItemIndex, SingleItemIndex)
function ref(d::DataMatrix, i::Real, j::Real)
    if d.na[i, j]
        return NA
    else
        return d.data[i, j]
    end
end

# dm[SingleItemIndex, MultiItemIndex]
function ref(x::DataMatrix, i::Real, col_inds::AbstractDataVector{Bool})
    ref(x, i, find(replaceNA(col_inds, false)))
end
function ref(x::DataMatrix, i::Real, col_inds::AbstractDataVector)
    ref(x, i, removeNA(col_inds))
end
# TODO: Make inds::AbstractVector
function ref(x::DataMatrix,
             i::Real,
             col_inds::Union(Vector, BitVector, Ranges))
    DataArray(x.data[i, col_inds], x.na[i, col_inds])
end

# dm[MultiItemIndex, SingleItemIndex]
function ref(x::DataMatrix, row_inds::AbstractDataVector{Bool}, j::Real)
    ref(x, find(replaceNA(row_inds, false)), j)
end
function ref(x::DataMatrix, row_inds::AbstractVector, j::Real)
    ref(x, removeNA(row_inds), j)
end
# TODO: Make inds::AbstractVector
function ref(x::DataMatrix,
             row_inds::Union(Vector, BitVector, Ranges),
             j::Real)
    DataArray(x.data[row_inds, j], x.na[row_inds, j])
end

# dm[MultiItemIndex, MultiItemIndex]
function ref(x::DataMatrix,
             row_inds::AbstractDataVector{Bool},
             col_inds::AbstractDataVector{Bool})
    ref(x, find(replaceNA(row_inds, false)), find(replaceNA(col_inds, false)))
end
function ref(x::DataMatrix,
             row_inds::AbstractDataVector{Bool},
             col_inds::AbstractDataVector)
    ref(x, find(replaceNA(row_inds, false)), removeNA(col_inds))
end
# TODO: Make inds::AbstractVector
function ref(x::DataMatrix,
             row_inds::AbstractDataVector{Bool},
             col_inds::Union(Vector, BitVector, Ranges))
    ref(x, find(replaceNA(row_inds, false)), col_inds)
end
function ref(x::DataMatrix,
             row_inds::AbstractDataVector,
             col_inds::AbstractDataVector{Bool})
    ref(x, removeNA(row_inds), find(replaceNA(col_inds, false)))
end
function ref(x::DataMatrix,
             row_inds::AbstractDataVector,
             col_inds::AbstractDataVector)
    ref(x, removeNA(row_inds), removeNA(col_inds))
end

# TODO: Make inds::AbstractVector
function ref(x::DataMatrix,
             row_inds::AbstractDataVector,
             col_inds::Union(Vector, BitVector, Ranges))
    ref(x, removeNA(row_inds), col_inds)
end

# TODO: Make inds::AbstractVector
function ref(x::DataMatrix,
             row_inds::Union(Vector, BitVector, Ranges),
             col_inds::AbstractDataVector{Bool})
    ref(x, row_inds, find(replaceNA(col_inds, false)))
end

# TODO: Make inds::AbstractVector
function ref(x::DataMatrix,
             row_inds::Union(Vector, BitVector, Ranges),
             col_inds::AbstractDataVector)
    ref(x, row_inds, removeNA(col_inds))
end

# TODO: Make inds::AbstractVector
function ref(x::DataMatrix,
             row_inds::Union(Vector, BitVector, Ranges),
             col_inds::Union(Vector, BitVector, Ranges))
    DataArray(x.data[row_inds, col_inds], x.na[row_inds, col_inds])
end

##############################################################################
##
## assign()
##
##############################################################################

# d[SingleItemIndex] = NA
function assign(da::DataArray, val::NAtype, i::Real)
	da.na[i] = true
end

# d[SingleItemIndex] = Single Item
function assign(da::DataArray, val::Any, i::Real)
	da.data[i] = val
	da.na[i] = false
end

# d[MultiIndex] = NA
function assign(da::DataArray{NAtype}, val::NAtype, inds::AbstractVector{Bool})
    throw(ArgumentError("DataArray{NAtype} is incoherent"))
end
function assign(da::DataArray{NAtype}, val::NAtype, inds::AbstractVector)
    throw(ArgumentError("DataArray{NAtype} is incoherent"))
end
function assign(da::DataArray, val::NAtype, inds::AbstractVector{Bool})
    da.na[find(inds)] = true
    return NA
end
function assign(da::DataArray, val::NAtype, inds::AbstractVector)
    da.na[inds] = true
    return NA
end

# d[MultiIndex] = Multiple Values
function assign(da::AbstractDataArray,
                vals::AbstractVector,
                inds::AbstractVector{Bool})
    assign(da, vals, find(inds))
end
function assign(da::AbstractDataArray,
                vals::AbstractVector,
                inds::AbstractVector)
    for (val, ind) in zip(vals, inds)
        da[ind] = val
    end
    return vals
end

# x[MultiIndex] = Single Item
function assign{T}(da::AbstractDataArray{T},
                   val::Union(Number, String, T),
                   inds::AbstractVector{Bool})
    assign(da, val, find(inds))
end
function assign{T}(da::AbstractDataArray{T},
                   val::Union(Number, String, T),
                   inds::AbstractVector)
    val = convert(T, val)
    for ind in inds
        da[ind] = val
    end
    return val
end
function assign(da::AbstractDataArray,
                val::Any,
                inds::AbstractVector{Bool})
    assign(da, val, find(inds))
end
function assign{T}(da::AbstractDataArray{T},
                   val::Any,
                   inds::AbstractVector)
    val = convert(T, val)
    for ind in inds
        da[ind] = val
    end
    return val
end

# dm[SingleItemIndex, SingleItemIndex] = NA
function assign(dm::DataMatrix, val::NAtype, i::Real, j::Real)
    dm.na[i, j] = true
    return NA
end

# dm[SingleItemIndex, SingleItemIndex] = Single Item
function assign(dm::DataMatrix, val::Any, i::Real, j::Real)
    dm.data[i, j] = val
    dm.na[i, j] = false
    return val
end

# dm[MultiItemIndex, SingleItemIndex] = NA
function assign(dm::DataMatrix,
                val::NAtype,
                row_inds::Union(Vector, BitVector, Ranges),
                j::Real)
    dm.na[row_inds, j] = true
    return NA
end

# dm[MultiItemIndex, SingleItemIndex] = Multiple Items
function assign{S, T}(dm::DataMatrix{S},
                      vals::Vector{T},
                      row_inds::Union(Vector, BitVector, Ranges),
                      j::Real)
    dm.data[row_inds, j] = vals
    dm.na[row_inds, j] = false
    return val
end

# dm[MultiItemIndex, SingleItemIndex] = Single Item
function assign(dm::DataMatrix,
                val::Any,
                row_inds::Union(Vector, BitVector, Ranges),
                j::Real)
    dm.data[row_inds, j] = val
    dm.na[row_inds, j] = false
    return val
end

# dm[SingleItemIndex, MultiItemIndex] = NA
function assign(dm::DataMatrix,
                val::NAtype,
                i::Real,
                col_inds::Union(Vector, BitVector, Ranges))
    dm.na[i, col_inds] = true
    return NA
end

# dm[SingleItemIndex, MultiItemIndex] = Multiple Items
function assign{S, T}(dm::DataMatrix{S},
                      vals::Vector{T},
                      i::Real,
                      col_inds::Union(Vector, BitVector, Ranges))
    dm.data[i, col_inds] = vals
    dm.na[i, col_inds] = false
    return val
end

# dm[SingleItemIndex, MultiItemIndex] = Single Item
function assign(dm::DataMatrix,
                val::Any,
                i::Real,
                col_inds::Union(Vector, BitVector, Ranges))
    dm.data[i, col_inds] = val
    dm.na[i, col_inds] = false
    return val
end

# dm[MultiItemIndex, MultiItemIndex] = NA
function assign(dm::DataMatrix,
                val::NAtype,
                row_inds::Union(Vector, BitVector, Ranges),
                col_inds::Union(Vector, BitVector, Ranges))
    dm.na[row_inds, col_inds] = true
    return NA
end

# dm[MultiIndex, MultiIndex] = Multiple Items
function assign{S, T}(dm::DataMatrix{S},
                      vals::Vector{T},
                      row_inds::Union(Vector, BitVector, Ranges),
                      col_inds::Union(Vector, BitVector, Ranges))
    dm.data[row_inds, col_inds] = vals
    dm.na[row_inds, col_inds] = false
    return val
end

# dm[MultiItemIndex, MultiItemIndex] = Single Item
function assign(dm::DataMatrix,
                val::Any,
                row_inds::Union(Vector, BitVector, Ranges),
                col_inds::Union(Vector, BitVector, Ranges))
    dm.data[row_inds, col_inds] = val
    dm.na[row_inds, col_inds] = false
    return val
end

##############################################################################
##
## Predicates
##
##############################################################################

isna(da::DataArray) = copy(da.na)
isnan(da::DataArray) = DataArray(isnan(da.data), copy(da.na))
isfinite(da::DataArray) = DataArray(isfinite(da.data), copy(da.na))
anyna(d::AbstractDataArray) = anyp(isna, d)
allna(d::AbstractDataArray) = allp(isna, d)

isna(a::AbstractArray) = falses(size(a))
anyna(a::AbstractArray) = false
allna(a::AbstractArray) = false

##############################################################################
##
## Generic iteration over AbstractDataArray's
##
##############################################################################

start(x::AbstractDataArray) = 1

function next(x::AbstractDataArray, state::Integer)
    return (x[state], state + 1)
end

function done(x::AbstractDataArray, state::Integer)
    return state > length(x)
end

##############################################################################
##
## Promotion rules
##
##############################################################################

promote_rule{T, T}(::Type{AbstractDataArray{T}},
                   ::Type{T}) = promote_rule(T, T)
promote_rule{S, T}(::Type{AbstractDataArray{S}},
                   ::Type{T}) = promote_rule(S, T)
promote_rule{T}(::Type{AbstractDataArray{T}}, ::Type{T}) = T

##############################################################################
##
## Conversion rules
##
##############################################################################

function convert{N}(::Type{BitArray{N}}, d::DataArray{BitArray{N}, N})
    throw(ArgumentError("Can't convert to BitArray"))
end

function convert{T, N}(::Type{BitArray{N}}, d::DataArray{T, N})
    throw(ArgumentError("Can't convert to BitArray"))
end

function convert{T, N}(::Type{Array{T, N}}, x::DataArray{T, N})
    if anyna(x)
        err = "Cannot convert DataArray with NA's to base type"
        throw(NAException(err))
    else
        return x.data
    end
end

function convert{S, T, N}(::Type{Array{S, N}}, x::DataArray{T, N})
    if anyna(x)
        err = "Cannot convert DataArray with NA's to desired type"
        throw(NAException(err))
    else
        return convert(S, x.data)
    end
end

##############################################################################
##
## Conversion convenience functions
##
##############################################################################

for f in (:int, :float, :bool)
    @eval begin
        function ($f){T}(da::DataArray{T})
            if anyna(da)
                err = "Cannot convert DataArray with NA's to desired type"
                throw(NAException(err))
            else
                ($f)(dv.data)
            end
        end
    end
end
for (f, basef) in ((:dataint, :int), (:datafloat, :float64), (:databool, :bool))
    @eval begin
        function ($f){T}(a::Array{T})
            DataArray(($basef)(a))
        end
        function ($f){T}(dv::DataArray{T})
            DataArray(($basef)(dv.data), copy(dv.na))
        end
    end
end

##############################################################################
##
## padNA
##
##############################################################################

function padNA(dv::AbstractDataVector, front::Int, back::Int)
    n = length(dv)
    res = similar(dv, front + n + back)
    for i in 1:n
        res[i + front] = dv[i]
    end
    return res
end

##############################################################################
##
## Conversion
##
##############################################################################

function vector(adv::AbstractDataVector, t::Type, replacement_val::Any)
    n = length(adv)
    res = Array(t, n)
    for i in 1:n
        if isna(adv[i])
            res[i] = replacement_val
        else
            res[i] = adv[i]
        end
    end
    return res
end
function vector(adv::AbstractDataVector, t::Type)
    n = length(adv)
    res = Array(t, n)
    for i in 1:n
        res[i] = adv[i]
    end
    return res
end
vector{T}(adv::AbstractDataVector{T}) = vector(adv, T)

function matrix(adm::AbstractDataMatrix, t::Type, replacement_val::Any)
    n, p = size(adm)
    res = Array(t, n, p)
    for i in 1:n
        for j in 1:p
            if isna(adm[i, j])
                res[i, j] = replacement_val
            else
                res[i, j] = adm[i, j]
            end
        end
    end
    return res
end
function matrix(adm::AbstractDataMatrix, t::Type)
    n, p = size(adm)
    res = Array(t, n, p)
    for i in 1:n
        for j in 1:p
            res[i, j] = adm[i, j]
        end
    end
    return res
end
matrix{T}(adm::AbstractDataMatrix{T}) = matrix(adm, T)

##############################################################################
##
## Hashing
##
## Make sure this agrees with is_equals()
##
##############################################################################

function hash(a::AbstractDataArray)
    h = hash(size(a)) + 1
    for i in 1:length(a)
        h = bitmix(h, int(hash(a[i])))
    end
    return uint(h)
end