text:
# overview

# so this isn't a neural net.  roughly, it's a series of matrix
# operations, where the "matrix" values are the same size as the
# playing board.

# the final result, after evaluating the operations, is treated as a
# probability distribution for where to make the next move.

# the input matrices are taken from the position itself - so, for
# example, a matrix where spaces are 0, black is 1 and white is -1.
# there are also two contant matrices, with the values 0 and 1 at
# every point.

# the operations available are:
# - convolution with a kernel
# - addition of two values scaled by some factor
# - product of two values raised to some power

# implementation

# to make meta-heuristics as effective as possible we want the
# representation to be very compact (everything counts), and for all
# values to be normalised (everything is useful).

# so we use bytes for values, but evaluate expressins using floats.
# the conversion from byte to float depends on the context.

# byte/serialisation format

# the idea here is to make this as flexible as possible, so that we
# can move data between different experiments.  its quite possible
# that any particular run / optimisation will be much more
# restrictive.

# header: "go" + "xp" (expression) + byte for version (version 0 below)
# length: unsigned two bytes (for following data, excluding header + length)
#         (any final operation that is incomplete because of the length
#         is simply discarded)
# operations: repeated until end (discarding any incomplete at end)
#   bit 6+7 of first byte gives type
#   00: jump
#     byte 0:
#       bit 6-7 (delegate tag)
#       bit 0-5 index to jump to (-n-1 so always -ve, relative)
#   01: kernel
#     byte 0:
#       bit 6-7: 0 (kernel tag)
#       bit 0-5: size in x and y (see code)
#     byte 1: index into input
#     byte 2: byte used for input outside board edges
#     byte 3: (mod size of grid) location of output within the kernel
#     byte 4: (mod number of combines) how to combine transforms
#     remaining bytes as kernel (see below for byte -> float)
#     all normalised by sqrt(n entries)
#   10, 11: arithmetic
#     byte 0:
#       bit 7: 1 (arithmetic tag)
#       bit 6:
#         0: additiion/subtraction
#         1: multiplication/division
#       bit 5-2: 4 flags (after first) mult/div and sine
#       bit 1-0: (+1) number of terms
#     remaining bytes as pairs:
#       byte 0: index into input
#       byte 1: scale (see below for byte -> float)
# indexing into input is
#   positive (mod available): direct index into input
#   negative (mod available): indirect index into input (relative to "here")
# general encoding from byte <-> float is
#   float(b) = signed(b) / scale (see below) with hard clipping

# chunking

# new fragments start every chunk (see below) bytes after the header.
# this allows changes without 'breaking' all following fragments.

# evaluation order

# originally fragments were evaluated in the order they appeared in the
# expression, with the final value being taken from the last expression.
# but this evaluated unused expressions, to the evaluation was made lazy.
# and then the ordering was reversed, so that the final value is taken
# from the first fragment defined, and references are to later fragments.
# this gives more persistence across merges, because early fragments
# are not affected by re-interpretation of modified later fragments.

# --- data structuring

const kernel_extra = 5
const arith_extra = 1
unpack_kernel_size(n::UInt8) = ((n & 0x38) >> 3) + 1, (n & 0x07) + 1
unpack_arithmetic_size(n::UInt8) = (n & 0x03) + 1

const scale = 64.0
b2f(b::Int8) = Float32(reinterpret(Int8, b) / scale)
f2b(f::Number) = Int8(min(127, max(-128, round(f * scale))))

const chunk = 8
function chunkend(s)
    n = (s.state - lheader - 1) % chunk
    if n > 0
        read(s, UInt8, min(available(s), chunk - n))
    else
        UInt8[]
    end
end

const given_atom = 16
const given_kern = 10
const header = map(UInt8, collect("goxp"))
const lheader = 7   # 4 chars, 1 version, 2 length

const transforms = ([1 0; 0 1], [0 1; -1 0], [-1 0; 0 -1], [0 -1; 1 0],
                    [1 0; 0 -1], [0 1; 1 0], [-1 0; 0 1], [0 -1; -1 0])

const combines = Dict{Int, Function}(0 => x(a, d) = sum(a, d) / length(transforms),
                                     1 => x(a, d) = prod(a, d),
                                     2 => x(a, d) = maximum(a, d) - minimum(a, d),
                                     3 => minimum,
                                     4 => maximum,
                                     5 => x(a, d) = minimum(abs(a), d),
                                     6 => x(p...) = maximum(abs(a), d))
const combine_names = Dict{Int, ASCIIString}(0 => "average",
                                             1 => "product",
                                             2 => "range",
                                             3 => "max",
                                             4 => "min",
                                             5 => "maxabs",
                                             6 => "minabs")

@enum Operation jump=0 kernel=1 addition=2 product=3 junk=4

@auto_hash_equals immutable Fragment
    # we only parse to this level because we need to preserve all bits
    # for evolutionary algorithms.  parsing to semantical level
    # discards some bits and life gets complicated.
    operation::Operation
    data::Vector{UInt8}
    padding::Vector{UInt8}
end

@auto_hash_equals type Expression
    version::UInt8
    # must be maintained by mutating operations
    length::UInt16
    fragment::Vector{Fragment}
    Expression(version::UInt8) = new(version, 7, Fragment[])
    Expression() = Expression(0x00)
end

immutable Index
    # indexing is simple for atoms - at level n there are
    # available+given_atom indices that can be used, where available
    # is i-1 (the number of entries below i) - see read_input.
    # for kernels, however, we are restricted to kernel enties.  so
    # 'available' depends on the number of kernels below, and indexing
    # is not linear.  to handle this, we construct these index tables,
    # which give the number of kernels below a given point and the
    # translation from input to direct index (both these include
    # given_kern).
    available::Vector{Int}
    index::Vector{Int}
end

function build_index(e::Expression)
    n = length(e.fragment)
    ok = zeros(Bool, n+given_atom)
    ok[1:given_kern] = true
    index = Index(zeros(Int, n), Int[1:given_kern;])
    for i = n:-1:1   # work through fragments in reverse
        index.available[i] = length(index.index) - given_kern
        if e.fragment[i].operation == kernel
            push!(index.index, given_kern + (n-i+1))
            ok[n-i+1 + given_atom] = true
        elseif e.fragment[i].operation == jump
            if ok[unpack_jump(e.fragment[i], n-i)]
                push!(index.index, given_kern + (n-i+1))
                ok[n-i+1 + given_atom] = true
            end
        elseif e.fragment[i].operation in (addition, product)
            for (input, scale, change) in unpack_arithmetic(e.fragment[i], n-i)
                if ok[input]
                    push!(index.index, given_kern + (n-i+1))
                    ok[n-i+1 + given_atom] = true
                    break
                end
            end
        end
    end
    index
end

function read_input(e::StatefulIterator{Vector{UInt8},Int}, available, given)
    # signed value
    input = read(e, Int8)
    n = given + available
    # abs(-128) = -128 for Int8
    index = abs(Int(input)) % n
    if input < 0
        # relative indexing for -ve values, -1 means next
        n - index
    else
        # absolute indexing for +ve values, 0 means bottom
        1 + index
    end
end

function unpack_jump(f::Fragment, available)
    masked = reinterpret(Int8, f.data[1] & 0x3f)
    # clunky, but allows same routine to be called
    # zero maps to -1
    read_input(StatefulIterator([reinterpret(UInt8, -masked-one(Int8))]),
               available, given_atom)
end

# positive input is not handled correctly
function pack_jump(input)
    # -1 maps to zero
    n = UInt8((-input-1) & 0x3f)
    Fragment(jump, [n], [])
end

function unpack_kernel(f::Fragment, i, index)
    e = StatefulIterator(f.data)
    nx, ny = unpack_kernel_size(read(e))
    idx = read_input(e, index.available[i], given_kern)
    input = index.index[idx]
    edge = b2f(read(e, Int8))
    cy, cx = map(Int, [divrem(read(e) % (nx*ny), nx)...] + [1,1])
    combine = read(e) % length(combines)
    kernel = reshape([b2f(read(e, Int8)) for i in 1:nx*ny], nx, ny)
    (input, edge, (cx, cy), combine, kernel)
end

# input is not handled correctly
function pack_kernel(input, edge, c, combine, coeffs)
    cx, cy = c
    nx, ny = size(coeffs)
    data = [UInt8(nx-1) << 3 | UInt8(ny-1)]
    push!(data, UInt8(input-1), f2b(edge), UInt8((cy-1)*nx + cx-1), UInt8(combine))
    append!(data, reinterpret(UInt8, map(f2b, reshape(coeffs, nx*ny))))
    Fragment(kernel, data, [])
end

function unpack_arithmetic(f::Fragment, available)
    e = StatefulIterator(f.data)
    tag = read(e)
    n = unpack_arithmetic_size(tag)
    flags = ((tag >> 2) & 0x0f) << 1
    function term(e)
        input = read_input(e, available, given_atom)
        scale = b2f(read(e, Int8))
        change = flags & 0x01 == 0x01
        flags = flags >> 1
        (input, scale, change)
    end
    [term(e) for i in 1:n]
end

unpack_product = unpack_arithmetic
unpack_addition = unpack_arithmetic

# -ve input is not handled correctly
function pack_arithmetic(coeffs...)
    n = length(coeffs)
    @assert 1 <= n <= 5
    data = [UInt8(n-1) | 0x80]
    flags, mask = 0x00, 0x01
    for (input, scale, change) in coeffs
        push!(data, UInt8(input-1), f2b(scale))
        flags = flags | (change ? mask : 0x00)
        mask = mask << 1
    end
    data[1] = data[1] | ((flags >> 1) << 2)
    data
end

pack_addition(coeffs...) = Fragment(addition, pack_arithmetic(coeffs...), [])

function pack_product(coeffs...)
    data = pack_arithmetic(coeffs...)
    data[1] = data[1] | 0x40
    Fragment(product, data, [])
end

function Base.push!(e::Expression, f::Fragment)
    # by appending junk to the previous padding we don't have to worry
    # about junk during evaluation.
    if f.operation == junk && length(e.fragment) > 0
        append!(e.fragment[end].padding, f.data)
    else
        push!(e.fragment, f)
    end
    e.length += length(f.data) + length(f.padding)
end

function unpack_expression(data::Vector{UInt8})
    d = StatefulIterator(data)
    @assert read(d, 4) == header
    e = Expression(read(d))
    n_bytes = read(d, UInt16)

    while !done(d)
        tag = peek(d)
        if tag & 0x80 == 0
            if tag & 0x40 == 0
                op = jump
                n = 1
            else
                op = kernel
                n = prod(unpack_kernel_size(tag)) + kernel_extra
            end
        else
            if tag & 0x40 == 0
                op = addition
                n = unpack_arithmetic_size(tag) * 2 + arith_extra
            else
                op = product
                n = unpack_arithmetic_size(tag) * 2 + arith_extra
            end
        end
        if n > available(d)
            f = Fragment(junk, read(d, available(d)), [])
        else
            f = Fragment(op, read(d, n), chunkend(d))
        end
        push!(e, f)
    end
    # this is all we support right now
    @assert e.version == 0x00
    @assert e.length == n_bytes
    @assert (e.length - lheader) % chunk == 0
    e
end

Expression(data::Vector{UInt8}) = unpack_expression(data)

function pack_expression(e::Expression)
    data = vcat(header, [0x00], [0x00, 0x00])
    for f in e.fragment
        fragment = vcat(f.data, f.padding)
        append!(data, fragment)
        # this is not needed when packing previously unpacked data,
        # but is when assembling a hand-written fragment
        n = length(data) - lheader
        if n % chunk != 0
            m = chunk - (n % chunk)
            append!(data, zeros(UInt8, m))
            append!(f.padding, zeros(UInt8, m))
            e.length += m
        end
    end
    # set length at end, with extra padding
    data[lheader-1:lheader] = reinterpret(UInt8, [e.length])
    @assert length(data) == e.length
    @assert (e.length - lheader) % chunk == 0
    data
end

# --- lazy evaluation

function lookup{N}(f, x, y, e, d::Array{Int8, 3}, input, edge, p::Position{N}, t::Point, index::Index)

    if input == 11
        if p.stats.moves == 0
            Float32(0.5)
        elseif isnull(p.stats.prev)
            zero(Float32)
        else
            one(Float32)
        end
    elseif input == 12
        zero(Float32)
    elseif input == 13
        one(Float32)
    elseif input == 14
        Float32(p.stats.score * Int(t))
    elseif input == 15
        Float32(p.stats.owned / (N*N - p.stats.stones))
    elseif input == 16
        Float32(p.stats.stones)
    elseif 1 <= x <= N && 1 <= y <= N
        if input == 1
            # -1 for "our" colour, 0 for space, 1 for opponent's colour
            Float32(-1 * Int(t) * Int(point(p, x, y)))
        elseif input == 2
            # distance to nearest stone (-ve for opponent's colour)
            Float32(Int(t) * p.flood.distance[x, y])
        elseif input == 3
            # -1 if owned by us, 1 owned by opponent, 0 otherwise
            b = p.space.border[x, y]
            if b == 0 || b == 3
                zero(Float32)
            elseif b == border_mask(t)
                -one(Float32)
            else
                one(Float32)
            end
        elseif input == 4
            i = p.groups.index[x, y]
            if i > 0
                Float32(p.groups.size[i])
            else
                edge
            end
        elseif input == 5
            i = p.groups.index[x, y]
            if i > 0
                Float32(p.groups.lives[i])
            else
                edge
            end
        elseif input == 6
            m = (N+1) / 2
            dx = abs(m - x)
            dy = abs(m - y)
            (dx + dy) / N
        elseif input == 7
            dx = min(x, N+1-x)
            dy = min(y, N+1-y)
            2 * (min(dx, dy) - 1) / N
        elseif input == 8
            2 * min(abs(x-y), abs(N+1-x-y)) / N
        elseif input == 9
            x / N
        elseif input == 10
            y / N
        else
            i = input-given_atom
            e[i] || evaluate(f, length(f)-i+1, p, t, e, d, index)
            b2f(d[x, y, i])
        end
    elseif 0 <= x <= N+1 && 0 <= y <= N+1
        edge
    else
        0
    end
end

function evaluate_jump{N}(f::Vector{Fragment}, idx, p::Position{N}, t::Point, e, d::Array{Int8, 3}, index::Index)
    n = length(f)
    input = unpack_jump(f[idx], n-idx)
    @forall i j N begin
        d[i, j, n-idx+1] = f2b(lookup(f, i, j, e, d, input, 0.0, p, t, index))
    end
    e[n-idx+1] = true
end

function evaluate_kernel{N}(f::Vector{Fragment}, idx, p::Position{N}, t::Point, e, d::Array{Int8, 3}, index::Index)
    n, nt = length(f), length(transforms)
    input, edge, (cx, cy), combine, coeffs = unpack_kernel(f[idx], idx, index)
    nx, ny = size(coeffs)
    my_acc = zeros(Float32, N, N, nt)
    @forall i j N begin
        for k in 1:nt
            for di in 1:nx
                for dj in 1:ny
                    ddi, ddj = [di-cx dj-cy] * transforms[k]
                    # the best way to understand this is to draw a picture.
                    # it's basically a coord change from one frame (coeffs)
                    # to the other (data).
                    my_acc[i,j,k] += coeffs[di, dj] * lookup(f, i+ddi, j+ddj, e, d, input, edge, p, t, index)
                end
            end
        end
    end
    d[1:end, 1:end, n-idx+1] = map(f2b, combines[combine](my_acc, 3))
    e[n-idx+1] = true
end

function evaluate_arithmetic{N}(f::Vector{Fragment}, idx, p::Position{N}, t::Point, e, d::Array{Int8, 3}, g, index::Index)
    n = length(f)
    coeffs = unpack_addition(f[idx], n-idx)
    my_acc = zeros(Float32, N, N)
    for (input, scale, change) in coeffs
        @forall i j N begin
            value = lookup(f, i, j, e, d, input, 0.0, p, t, index)
            my_acc[i, j] += scale * (change ? g(value) : value)
        end
    end
    @forall i j N begin
        d[i, j, n-idx+1] = f2b(my_acc[i, j] / sqrt(length(coeffs)))
    end
    e[n-idx+1] = true
end

function evaluate(f::Vector{Fragment}, i, p::Position, t::Point, e, d::Array{Int8, 3}, index::Index)
    fragment = f[i]
    if fragment.operation == jump
        evaluate_jump(f, i, p, t, e, d, index)
    elseif fragment.operation == kernel
        evaluate_kernel(f, i, p, t, e, d, index)
    elseif fragment.operation == addition
        evaluate_arithmetic(f, i, p, t, e, d, x -> sin(pi * x / 128), index)
    elseif fragment.operation == product
        evaluate_arithmetic(f, i, p, t, e, d, x -> 1/x, index)
    end
end

function evaluate{N}(e::Expression, p::Position{N}, t::Point)
    index = build_index(e)
    n = length(e.fragment)
    evaluated = zeros(Bool, n)
    output_data = zeros(Int8, N, N, n)
    evaluate(e.fragment, 1, p, t, evaluated, output_data, index)
    # output is final frame (most complex)
    map(b2f, output_data[:,:,n])
end

function evaluate(e::Array{UInt8, 1}, p::Position, t::Point)
    evaluate(Expression(e), p, t)
end

# --- move extraction

# the above generates a matrix of 'values' that can be positive or
# negative.  the best interpretation i can give to that is log(prob).
# and in a sense i get to choose, because if i decide that and evolve
# on that basis, then that is what i get (in an ideal world).

# however, it's damn spiky.  so any random sampling will be dominated
# by the largest numbers.  so instead, we just pick off the largest
# scoring locations (and then exclude invalid moves etc).

function moves{N}(e::Array{UInt8, 1}, h::Vector{Position{N}}, t::Point)
    p = h[end]
    logp = evaluate(e, p, t)
    indexed = reshape([(logp[i, j], (i, j)) for i in 1:N, j in 1:N], N*N)
    positive = filter(x -> x[1] > 0, indexed)
    possible = filter(x -> valid(h, t, x[2]...), positive)
    m = map(x -> x[2], sort(possible; by=x -> x[1], rev=true))
    if p.stats.moves == 0
        r = hex2bytes(name(e))[1]
        t = transforms[1 + (r % length(transforms))]
        o = [Int((N+1) / 2) Int((N+1) / 2)]
        function transform(xy)
            xy2 = hcat(xy...)
            xy3 = o + ((xy2 - o) * t)
            tuple(xy3...)
        end
        m = map(transform, m)
    end
    m
end

# --- analysis

const lookup_dict = Dict{Int,ASCIIString}(1 => "point",
                                          2 => "flood",
                                          3 => "owner",
                                          4 => "size",
                                          5 => "lives",
                                          6 => "centre",
                                          7 => "edge",
                                          8 => "diag",
                                          9 => "x",
                                          10 => "y",
                                          11 => "prev",
                                          12 => "zero",
                                          13 => "one",
                                          14 => "score",
                                          15 => "owned%",
                                          16 => "stones%")

function lookup_name(input, n, given)
    if input > given
        n-(input - given)+1
    else
        lookup_dict[input]
    end
end

int(x) = typeof(x) <: Integer

function dump_jump(i, f)
    n = length(f)
    input = lookup_name(unpack_jump(f[i], n-i), n, given_atom)
    println(" $i jump to $(input)")
    filter(int, [input])
end

function cell(value, lo, hi, symbols)
    if value <= lo
        symbols[1]
    elseif value >= hi
        symbols[3]
    else
        symbols[2]
    end
end

function sketch_kernel(kernel, cx, cy)
    nx, ny = size(kernel)
    ord = sort(reshape(kernel, nx*ny))
    lo = nx*ny >=9 ? ord[cld(length(ord), 3)] : ord[1]
    hi = nx*ny >=9 ? ord[length(ord) - fld(length(ord), 3)] : ord[end]
    for j in 1:ny
        for i in 1:nx
            if i == cx && j == cy
                print(cell(kernel[i,j], lo, hi, "O.@"))
            else
                print(cell(kernel[i,j], lo, hi, "- +"))
            end
        end
        println()
    end
end

function dump_kernel(i, f, index)
    (input, edge, (cx, cy), combine, kernel) = unpack_kernel(f[i], i, index)
    input = lookup_name(input, length(f), given_kern)
    println(" $i kernel($(combine_names[combine]), $(input))")
    sketch_kernel(kernel, cx, cy)
    filter(int, [input])
end

term2str(n, input, scale, change) = string("$(scale)*$(lookup_name(input, n, given_atom))", change ? "!" : "")

function dump_arithmetic(i, f, op)
    n = length(f)
    terms = unpack_arithmetic(f[i], n-i)
    expr = join(map(x -> term2str(n, x...), terms), " $op ")
    println(" $i $expr")
    filter(int, map(x -> lookup_name(x[1], n, given_atom), terms))
end

function dump_fragment(i, f, used, index)
    if f[i].operation == jump
        refs = dump_jump(i, f)
    elseif f[i].operation == kernel
        refs = dump_kernel(i, f, index)
    elseif f[i].operation == addition
        refs = dump_arithmetic(i, f, "+")
    elseif f[i].operation == product
        refs = dump_arithmetic(i, f, "*")
    else
        refs = Int[]
    end
    for ref in refs
        push!(used, ref)
    end
end

function dump_expression(e::Expression)
    index = build_index(e)
    used = Set([1])
    for (i, fragment) in collect(enumerate(e.fragment))
        if i in used
            dump_fragment(i, e.fragment, used, index)
        end
    end
end

dump_expression(e::Vector{UInt8}) = dump_expression(unpack_expression(e))


text:
using Test
using StaticArrays
using FiniteVolumes
using FiniteVolumes.CartesianMeshes: Half, dx, _direction

@testset "Mesh" begin
    @test dx(CartesianMesh(10)) == SVector(0.1)
    @test dx(CartesianMesh(10, 10)) == SVector(0.1, 0.1)
    @test dx(CartesianMesh(SVector(1.0, 0.0), SVector(2.0, 2.0), SVector(10, 20))) == SVector(0.1, 0.1)

    @test FiniteVolumes.nb_cells(CartesianMesh(10)) == 10
    @test FiniteVolumes.nb_cells(CartesianMesh(5, 5)) == 25
    @test FiniteVolumes.nb_cells(CartesianMesh(2, 5)) == 10

    @test CartesianIndex(3, 3) in FiniteVolumes.all_cells(CartesianMesh(10, 10))
    @test FiniteVolumes.all_cells(CartesianMesh(2)) |> length == 2
    @test FiniteVolumes.all_cells(CartesianMesh(2, 2)) |> length == 4

    @test FiniteVolumes.inner_faces(CartesianMesh(2)) |> collect |> length == 1
    @test FiniteVolumes.inner_faces(PeriodicCartesianMesh(10, 10)) |> collect |> length == 200

    for m in [CartesianMesh(10), PeriodicCartesianMesh(3), CartesianMesh(5, 2), PeriodicCartesianMesh(10, 4)]
        @test FiniteVolumes.inner_faces(m) |> collect |> length == FiniteVolumes.nb_inner_faces(m)
        @test FiniteVolumes.boundary_faces(m) |> collect |> length == FiniteVolumes.nb_boundary_faces(m)
    end

    @test FiniteVolumes.boundary_faces(CartesianMesh(2)) |> collect |> length == 2
    # @test FiniteVolumes.boundary_faces(CartesianMesh(5, 5)) |> collect |> length == 20

    @test _direction((Half(3), Half(2))) == 1
    @test _direction((Half(2), Half(3))) == 2

    @test FiniteVolumes.cells_next_to_inner_face(PeriodicCartesianMesh(10), (Half(19),)) == (CartesianIndex(9), CartesianIndex(10))
    @test FiniteVolumes.cells_next_to_inner_face(PeriodicCartesianMesh(10), (Half(19),)) == (CartesianIndex(9), CartesianIndex(10))
    @test FiniteVolumes.cells_next_to_inner_face(PeriodicCartesianMesh(10), (Half(1),)) == (CartesianIndex(10), CartesianIndex(1))

    @test FiniteVolumes.cells_next_to_inner_face(CartesianMesh(10, 10), (Half(3), Half(2))) == (CartesianIndex(1, 1), CartesianIndex(2, 1))
    @test FiniteVolumes.cells_next_to_inner_face(CartesianMesh(10, 10), (Half(2), Half(3))) == (CartesianIndex(1, 1), CartesianIndex(1, 2))
    # @test FiniteVolumes.cell_next_to_boundary_face(CartesianMesh(10, 10), (Half(1), Half(2))) == CartesianIndex(1, 1)
    # @test FiniteVolumes.cell_next_to_boundary_face(CartesianMesh(10, 10), (Half(21), Half(2))) == CartesianIndex(10, 1)

    @test FiniteVolumes.cell_volume(CartesianMesh(10), CartesianIndex(3)) == 0.1
    @test FiniteVolumes.cell_volume(CartesianMesh(2, 5), CartesianIndex(1, 1)) == 0.1

    @test FiniteVolumes.face_area(CartesianMesh(5, 2), (Half(2), Half(1))) == 0.2
    @test FiniteVolumes.face_area(CartesianMesh(5, 2), (Half(1), Half(2))) == 0.5

    @test FiniteVolumes.cell_center(CartesianMesh(5, 5), CartesianIndex(1, 1)) ≈ [0.1, 0.1]
    @test FiniteVolumes.cell_center(CartesianMesh(5, 5), CartesianIndex(2, 1)) ≈ [0.3, 0.1]
    @test FiniteVolumes.face_center(CartesianMesh(5, 5), (Half(3), Half(2))) == [0.2, 0.1]
    @test FiniteVolumes.face_center(CartesianMesh(5, 5), (Half(2), Half(3))) == [0.1, 0.2]

    # @test FiniteVolumes.cell_corners(CartesianMesh(5, 5), ).bottom_left == @SVector [0.0, 0.0]
    # @test FiniteVolumes.cell_corners(CartesianMesh(5, 5), ).top_right == @SVector [1.0, 1.0]
end

# @btime cells_next_to_inner_face(mesh, i_face) setup=(mesh=CartesianMesh(10); i_face=(Half(19),))
# @btime cells_next_to_inner_face(mesh, i_face) setup=(mesh=PeriodicCartesianMesh(10); i_face=(Half(19),))
# @btime cell_center(mesh, (3,)) setup=(mesh=CartesianMesh(10))
# @btime face_center(mesh, (Half(3),)) setup=(mesh=CartesianMesh(10))

ext: jl


text:
# export play_game_with_policy

# using Flux:softmax
# using StatsBase: wsample

# const EPSILON = 0.1

# function play_game_with_policy(policy; verbose=false, goal=8, epsilon=EPSILON)
#     board = initboard()

#     histories = Array{Int8, 2}[]
#     dir_histories = Dirs[]
#     while true
#         res = move.(Ref(board), DIRS)
#         # possible = [p for (_, (_, p)) in res]
#         possible::Vector{Bool} = [p for (_, (_, p)) in res]

#         if all(!, possible)
#             return histories, dir_histories, board
#         end

#         choices = (1:4)[possible]

#         if rand() < epsilon
#             epsed = true
#             i = rand(choices)
#         else
#             epsed = false
#             tmp = policy(board)
#             prob = reshape(tmp, 4)
#             if any(isnan, prob)
#                 error("wtf")
#             end
#             # sample a direction from policy
#             i = wsample(choices, prob)
#         end

#         dir = DIRS[i]
#         if verbose
#             println("****************************begin: one move")
#             display(board)
#             println(dir)
#             println(DIRS[choices])
#             if epsed
#                 "epsilon'ed"
#             else
#                 println(prob)
#             end
#         end
#         push!(histories, board)
#         push!(dir_histories, dir)
#         board = res[i][1]

#         if verbose
#             display(board)
#             println("****************************end: one move")
#         end
#         gen_new_tile!(board)
#     end
# end


text:
include("logger.jl")
x = include("load_data.jl")
include("k-means.jl")

## visualize raw data
#using Gadfly
#plot(x=x[1, :], y=x[2, :])

function compute_final_error(centroids, nearest_centroids)
    sum = 0.0
    for (i, c) in nearest_centroids
        real_centroid = get_real_centroid(i)
        diff = real_centroid - centroids[c]
        sum += dot(diff, diff)
    end

    sum / length(nearest_centroids)
end

function cluster_stat(k, nearest_centroids)
    counts = zeros(k)
    for (i, c) in nearest_centroids
        counts[c] += 1
    end
    counts
end

function get_real_centroid(i)
    if i <= 200
        [1., -1.]
    elseif i <= 400
        [5.5, -4.5]
    elseif i <= 600
        [1., 4.]
    elseif i <= 800
        [6., 4.5]
    else
        [9., 0.0]
    end
end

function main()
    inputs = Vector{Float64}[x[:, i] for i = 1:size(x)[2]]

    centroids, nearest_centroids = kmeans(5, inputs, lowerbound=0.000001)

    println("=====\nfinal=$(formatnum(centroids))")
    error = compute_final_error(centroids, nearest_centroids)
    dbglog("error=$(formatnum(error))")
    stat = cluster_stat(5, nearest_centroids)
    dbglog("stat=$(formatnum(stat))")

    return centroids, nearest_centroids
end

if !isinteractive()
    main()
end

max_stars_repo_name: zxteloiv/julia-ex

max_issues_repo_name: zxteloiv/julia-ex

max_forks_repo_name: zxteloiv/julia-ex


text:
# -*- coding: utf-8 -*-
# + {}
module PetriModels
import Base: show
using Petri
using MacroTools
import MacroTools: prewalk
using Catlab.WiringDiagrams
import Catlab.Graphics.Graphviz: Graph

using SemanticModels
using SemanticModels.CategoryTheory
import SemanticModels.CategoryTheory: ⊔, FinSetMorph, pushout
import SemanticModels: model

export PetriModel, model, rewrite!, PetriSpan, DPOProblem, solve, dropdown, equnion, ⊔

struct PetriModel <: AbstractModel
  model::Petri.Model
end

function model(::Type{PetriModel}, m::Petri.Model)
    return PetriModel(m)
end

function equnion(a::Vector, b::Vector)
    x = copy(a)
    for item in b
        if !any(item2 -> isequal(item2, item), x)
            push!(x, item)
        end
    end
    return x
end

dictReplace(item, dict) = prewalk(i -> i in keys(dict) ? dict[i] : i, item)

function ⊔(gModel::PetriModel, hModel::PetriModel)
    g = gModel.model
    h = hModel.model
    dict = Dict(i => any(isequal.(i, g.S)) ? i : Symbol(i, "′") for i in h.S)
    newS = [dictReplace(n, dict) for n in h.S]
    newΔ = [(dictReplace(first(n), dict), dictReplace(last(n), dict)) for n in h.Δ]
    PetriModel(Petri.Model(union(g.S, newS),
                           equnion(g.Δ, newΔ)))
end

function (f::FinSetMorph)(gModel::G) where G <: PetriModel
    g = gModel.model
    length(dom(f)) == length(g.S) || throw(DomainError(g.S))
    ϕ = func(f)
    outS = 1:length(Set(ϕ.m.fun))
    outΔ = Vector{Tuple{Vector{Int}, Vector{Int}}}()
    for t in g.Δ
        ins = ϕ.(t[1])
        outs = ϕ.(t[2])
        push!(outΔ, (ins,outs))
     end
    out = deepcopy(g)
    PetriModel(Petri.Model(outS, outΔ))
end

# TODO: this version uses model toolkit variables, we want to delete it
# function (f::FinSetMorph)(gModel::G) where G <: PetriModel
#     g = gModel.model
#     length(dom(f)) == length(g.S) || throw(DomainError(g.S))
#     ϕ = func(f)
#     outS = Array{Operation}(undef, length(Set(ϕ.m.fun)))
#     for i in dom(f)
#         outS[ϕ(i)] = g.S[i]
#     end
#     out = deepcopy(g)
#     PetriModel(Petri.Model(outS, out.Δ, out.Λ, out.Φ))
# end

# TODO: this looks like coproduct not pushout.

"""    pushout(m1::Model, m2::Model)

compute the CT pushout of two models.
"""
function pushout(pModel::PetriModel, pModel2::PetriModel)
    pm = pModel.model
    pm2 = pModel2.model
    states = equnion(pm.S, pm2.S)
    Δ = equnion(pm.Δ, pm2.Δ)
    return model(PetriModel, Petri.Model(states, Δ))
end

"""
    dropdown(l::Model, c::Model, l′::Model)

compute c′ given l, c, l′, the formula for petrinets is

T_{c'} = T_{l'} \\setdiff f(T_l) \\cup a(T_c)
"""
function dropdown(pL::PetriModel, pC::PetriModel, pL′::PetriModel)
    pl = pL.model
    pc = pC.model
    pl′ = pL′.model
    states = union(pl′.S, pc.S)
    Δ = union(setdiff(pl′.Δ, pl.Δ), pc.Δ)
    return model(PetriModel, Petri.Model(states, Δ))
end

end


text:
using NeuralQuantum, Random, Plots

Random.seed!(Random.GLOBAL_RNG, 1234)

# Parameters of the transverse field ising model
dims = [5,5]
h = 3.0
J = 1.0

# Constructs the Hilbert space for N 1//2 spins.
N    = prod(dims)
hilb = HomogeneousSpin(N, 1//2)

coord(i,j) = i + (j-1)*dims[2]

# Builds the hamiltonian
H = LocalOperator(hilb)
for i=1:dims[1]
    for j=1:dims[2]
        p = coord(i,j)
        global H  -= h * sigmax(hilb, p)
        hop = J * sigmaz(hilb, p) * sigmaz(hilb, coord(mod(i, dims[1])+1,j)) +
                J * sigmaz(hilb, p) * sigmaz(hilb, coord(i, mod(j, dims[2])+1))
        global H += hop
    end
end

net  = RBM(Float32, N, 1, af_logcosh)
init_random_pars!(net, sigma=0.01)

sampl = MetropolisSampler(LocalRule(), 75, N, burn=100)
algo  = SR(ϵ=(0.1), algorithm=sr_qlp)

is = BatchedSampler(net, sampl, H, algo;
                    batch_sz=16)

optimizer = Optimisers.Descent(0.1)

Evalues = Float64[];
Eerr = Float64[];
for i=1:1000
    ldata, prec = sample!(is)
    ob = compute_observables(is)

    println("$i - $ldata")

    push!(Evalues, real(ldata.mean))
    push!(Eerr, ldata.error)
    grad = precondition!(prec, algo, i)
    Optimisers.update!(optimizer, net, grad)
end

plot(Evalues, yerr=Eerr)

exact = -80.13310152422413
hline!([exact, exact])


text:
# partial constraint development
using Test
using IncrementalInference
# going to introduce two new constraint mutable structs
using Statistics
using TensorCast
using Manifolds

import IncrementalInference: getManifold, getSample

##

mutable struct DevelopPartial{P <: Tuple} <: AbstractPrior
  x::Distribution
  partial::P
end
getSample(cf::CalcFactor{<:DevelopPartial}) = rand(cf.factor.x, 1)

#
mutable struct DevelopDim2 <: AbstractPrior
  x::Distribution
end
getSample(cf::CalcFactor{<:DevelopDim2}) = rand(cf.factor.x, 1)

mutable struct DevelopPartialPairwise <: AbstractRelativeMinimize
  x::Distribution
  partial::Tuple
  DevelopPartialPairwise(x::Distribution) = new(x, (2,))
end
getManifold(::IIF.InstanceType{DevelopPartialPairwise}) = TranslationGroup(1)

getSample(cf::CalcFactor{<:DevelopPartialPairwise}) = rand(cf.factor.x, 1)

function (dp::CalcFactor{<:DevelopPartialPairwise})(meas,
                                                    x1,
                                                    x2  )
  #
  return meas[1] - (x2[2]-x1[2])
end

## prep

N=100

pr = DevelopDim2(MvNormal([0.0;0.0], diagm([0.01;0.01])))
dp = DevelopPartial(Normal(2.0, 1.0),(1,))

## build graph

fg = initfg()

v1 = addVariable!(fg,:x1,ContinuousEuclid{2}(),N=N)
f1  = addFactor!(fg,[:x1],pr)
f2  = addFactor!(fg,[:x1],dp, graphinit=false)

doautoinit!(fg, :x1)

##

@testset "test evaluation of full constraint prior" begin

pts_, _ = evalFactor(fg, f1, v1.label, N=N)
@cast pts[i,j] := pts_[j][i]
@test size(pts,1) == 2
@test size(pts,2) == N
@test norm(Statistics.mean(pts,dims=2)[1] .- [0.0]) < 0.3

end

##

@testset "test evaluation of partial constraint prior" begin
##

memcheck_ = getVal(v1)
@cast memcheck[i,j] := memcheck_[j][i]

X1pts_ = getVal(v1)
@cast X1pts[i,j] := X1pts_[j][i]
pts_ = approxConv(fg, getLabel(f2), :x1, N=N)
@cast pts[i,j] := pts_[j][i]

@test size(pts, 1) == 2
@test size(pts, 2) == N
@test norm(Statistics.mean(pts,dims=2)[1] .- [2.0]) < 0.75
# ensure the correct response from
@test norm(X1pts[1,:] - pts[1,:]) > 2.0
@test norm(X1pts[2,:] - pts[2,:]) < 1e-10
@test norm(X1pts - memcheck) < 1e-10

##
end

@testset "check that partials are received through convolutions of prior" begin
##

# check that a partial belief is obtained
X1_ = approxConvBelief(fg, :x1f2, :x1)

@test isPartial(X1_)

##
end

@testset "test solving of factor graph" begin
##

getSolverParams(fg).N = N
tree = solveTree!(fg)
pts_ = getVal(fg, :x1)
@cast pts[i,j] := pts_[j][i]

@test norm(Statistics.mean(pts,dims=2)[1] .- [0.0]) < 0.4
@test norm(Statistics.mean(pts,dims=2)[2] .- [0.0]) < 0.4

##
end
# plotKDE(getBelief(fg, :x1),levels=3)

## partial relative gradient and graph

v2 = addVariable!(fg,:x2,ContinuousEuclid{2},N=N)

dpp = DevelopPartialPairwise(Normal(10.0, 1.0))
f3  = addFactor!(fg,[:x1;:x2],dpp)

dp2 = DevelopPartial( Normal(-20.0, 1.0), (1,) )
f4  = addFactor!(fg,[:x2;], dp2, graphinit=false)

doautoinit!(fg, :x2)

##

@testset "test partial info per coord through relative convolution (conditional)" begin
##

one_meas = [10.0;]
pts = ([0;0.0], [0;10.0])
gradients = FactorGradientsCached!(dpp, (ContinuousEuclid{2}, ContinuousEuclid{2}), one_meas, pts);

##

# check that the gradients can be calculated
J = gradients(one_meas, pts...)

@test size(J) == (4,4)
@test norm(J - [0 0 0 0; 0 0 0 1; 0 0 0 0; 0 1 0 0] ) < 1e-4

## check perturbation logic

prtb = calcPerturbationFromVariable(gradients, [1=>[1;1]])

# self variation is taken as 0 at this time
@test isapprox( prtb[1], [0;0] )
# variable 1 influences 2 only through partial dimension 2 (as per DevelopPartialPairwise)
@test isapprox( prtb[2], [0;1] )

##  test evaluation through the convolution operation withing a factor graph

# add relative IPC calculation inside evalFactor
bel = approxConvBelief(fg, getLabel(f3), :x2)
@test isPartial(bel)

##
end

##

@testset "test evaluation of multiple simultaneous partial constraints" begin
global fg

##

initAll!(fg)
valx2_ = getVal(fg, :x2)
@cast valx2[i,j] := valx2_[j][i]
pts_ = approxConv(fg, :x1x2f1, :x2, N=N)
@cast pts[i,j] := pts_[j][i]
@test size(pts,1) == 2
@test norm(Statistics.mean(pts,dims=2)[2] .- [10.0]) < 3.0
@test norm(valx2[1,:] - pts[1,:]) < 1e-5

pts_ = approxConv(fg, :x2f1, :x2, N=N)
@cast pts[i,j] := pts_[j][i]
@test size(pts,1) == 2
@test norm(Statistics.mean(pts,dims=2)[1] .- [-20.0]) < 0.75
@test (Statistics.std(pts,dims=2)[1] .- 1.0) < 0.4

##

end

##

# keep previous values to ensure funciton evaluation is modifying correct data fields

@warn "restore findRelatedFromPotential as testset!"
# @testset "test findRelatedFromPotential..." begin
# global v2, fg, f3, f4, N

thefac = getFactor(fg, :x1x2f1)

X2lpts_ = getVal(getVariable(fg, :x2))
@cast X2lpts[i,j] := X2lpts_[j][i]
keepaside, = findRelatedFromPotential(fg, thefac, :x2, N=N)
@test Ndim(keepaside) == 2
lpts_ = getPoints(keepaside, false)
@cast lpts[i,j] := lpts_[j][i]
@test length(lpts_) == N

@show X2lpts[2,95:100]
@show lpts[2,95:100]
@show getPoints(keepaside)

##

# DevelopPartialPairwise must only modify the second dimension of proposal distribution on X2
@test norm(X2lpts[1,:] - lpts[1,:]) < 1e-10
# @test norm(X2lpts[2,:] - lpts[2,:]) > 1e-10 # 10*N*0.5 # why so big?
memcheck_ = getVal(v2)
@cast memcheck[i,j] := memcheck_[j][i]
@test norm(X2lpts - memcheck) < 1e-10

X2lpts_ = getVal(v2)
@cast X2lpts[i,j] := X2lpts_[j][i]
p4, = findRelatedFromPotential(fg, f4, v2.label, N=N)
@test Ndim(p4) == 2
lpts_ = getPoints(keepaside, false)
@cast lpts[i,j] := lpts_[j][i]
@test length(lpts_) == N

# DevelopPartialPairwise must only modify the second dimension of proposal distribution on X2
@test norm(X2lpts[1,:] - lpts[1,:]) < 1e-10
@test norm(X2lpts[2,:] - lpts[2,:]) > 1e-10 # 10*N*0.5 # why so big?
memcheck_ = getVal(v2)
@cast memcheck[i,j] := memcheck_[j][i]
@test norm(X2lpts - memcheck) < 1e-10

# end

##

@testset "test belief prediction with partials..." begin

##

global v2, fg

# partial prior
X2pts_ = getVal(v2)
@cast X2pts[i,j] := X2pts_[j][i]
# NOTE, SUPER IMPORTANT, predictbelief returns full dimension points (even if only partials are sent in for proposals)
val_, = predictbelief(fg, v2, [f4], N=N)
@cast val[i,j] := val_[j][i]
@show X2pts_[1]';
@show val_[1]';
@test norm(X2pts[2,:] - val[2,:]) < 1e-10
@test 0.0 < norm(X2pts[1,:] - val[1,:])
@test norm(Statistics.mean(val[1,:]) .+ 20.0) < 0.75

# partial pairwise
X2pts_ = getVal(v2)
@cast X2pts[i,j] := X2pts_[j][i]
val_, = predictbelief(fg, v2, [f3], N=N)
@cast val[i,j] := val_[j][i]
@test norm(X2pts[1,:] - val[1,:]) < 1e-10
@test 0.0 < norm(X2pts[2,:] - val[2,:])
val2_ = getVal(v1)
@cast val2[i,j] := val2_[j][i]
@test abs(Statistics.mean(val[2,:] - val2[2,:]) .- 10.0) < 0.75

##

# combination of partials
val_, = predictbelief(fg, v2, [f3;f4], N=N)
@cast val[i,j] := val_[j][i]
# plotKDE(kde!(val),levels=3)
@test norm(Statistics.mean(val,dims=2)[1] .- [-20.0]) < 1
@test_broken norm(Statistics.mean(val,dims=2)[2] .- [10.0]) < 0.01
@test (Statistics.std(val,dims=2)[1] .- 1.0) < 3.0
@test (Statistics.std(val,dims=2)[2] .- 1.0) < 3.0

##

getSolverParams(fg).N = N
tree = solveTree!(fg)

pts_ = getVal(fg, :x1)
@cast pts[i,j] := pts_[j][i]
@test norm(Statistics.mean(pts,dims=2)[1] .- [0.0]) < 0.5
@test norm(Statistics.mean(pts,dims=2)[2] .- [0.0]) < 0.5

pts_ = getVal(fg, :x2)

ppe = getPPE(fg, :x2).mean

X2 = getBelief(fg, :x2)

# check mean is close
@test_broken isapprox(mean(X2), [-20;10], atol=0.01)

# check covariance is close too
@test 0 < AMP.calcCovarianceBasic(getManifold(X2), getPoints(X2))

@cast pts[i,j] := pts_[j][i]
@test (Statistics.std(pts,dims=2)[1]-1.0) < 3.0
@test (Statistics.std(pts,dims=2)[2]-1.0) < 3.0

##

end

# plotKDE(getBelief(fg, :x2),levels=3)

# spyCliqMat(tree, :x2)

#


text:
using JLD
using Mustache

# counts = load("data/compare_Friday_24_Mar_00_52.jld", "counts")
# steps = load("data/compare_Friday_24_Mar_00_52.jld", "steps")

# without the
# counts = load("data/compare_Sunday_26_Mar_17_56.jld", "counts")
# steps = load("data/compare_Sunday_26_Mar_17_56.jld", "steps")

# solver_keys = load("data/compare_Wednesday_29_Mar_10_05.jld", "solver_keys")
counts = load("data/compare_Wednesday_29_Mar_10_05.jld", "counts")
rewards = load("data/compare_Wednesday_29_Mar_10_05.jld", "rewards")
steps = load("data/compare_Wednesday_29_Mar_10_05.jld", "steps")

for k in keys(counts)
    # cts_per_step[k] = sum(counts[k])/sum(steps[k])
    rew = rewards[k]
    sem = std(rew)/sqrt(length(rew))
    cps = sum(counts[k])/sum(steps[k])
    println("$k & $(mean(rew)) \\pm $sem & $cps")
end

counts = load("data/compare_Wednesday_29_Mar_19_04.jld", "counts")
rewards = load("data/compare_Wednesday_29_Mar_19_04.jld", "rewards")
steps = load("data/compare_Wednesday_29_Mar_19_04.jld", "steps")

for k in keys(counts)
    # cts_per_step[k] = sum(counts[k])/sum(steps[k])
    rew = rewards[k]
    sem = std(rew)/sqrt(length(rew))
    cps = sum(counts[k])/sum(steps[k])
    println("$k & $(mean(rew)) \\pm $sem & $cps")
end

times = load("data/compare_Friday_31_Mar_20_23.jld", "times")
steps = load("data/compare_Friday_31_Mar_20_23.jld", "steps")

for k in keys(times)
    tps = sum(times[k])/sum(steps[k])
    println("$k: $tps")
end

# @show 50_000*mean(cts_per_step["pomcpow"])/mean(cts_per_step["modified_pomcp"])

# mt"""
# \begin{tabular}{l}
# \toprule
# Algorithm & \$n\$ & Avg.
# \midrule
#
# \bottomrule
# \end{tabular}
# """


text:
using DataStructures
using Base.Test

# construction

@test typeof(OrderedDict()) == OrderedDict{Any,Any}
@test typeof(OrderedDict([(1,2.0)])) == OrderedDict{Int,Float64}
@test typeof(OrderedDict([("a",1),("b",2)])) == OrderedDict{ASCIIString,Int}
if VERSION >= v"0.4.0-dev+980"
    @test typeof(OrderedDict(Pair(1, 1.0))) == OrderedDict{Int,Float64}
    @test typeof(OrderedDict(Pair(1, 1.0), Pair(2, 2.0))) == OrderedDict{Int,Float64}
    @test typeof(OrderedDict(Pair(1, 1.0), Pair(2, 2.0), Pair(3, 3.0))) == OrderedDict{Int,Float64}
end

# empty dictionary
d = OrderedDict{Char, Int}()
@test length(d) == 0
@test isempty(d)
@test_throws KeyError d['c'] == 1
d['c'] = 1
@test !isempty(d)
empty!(d)
@test isempty(d)

# access, modification

for c in 'a':'z'
    d[c] = c-'a'+1
end

@test (d['a'] += 1) == 2
@test 'a' in keys(d)
@test haskey(d, 'a')
@test get(d, 'B', 0) == 0
@test !('B' in keys(d))
@test !haskey(d, 'B')
@test pop!(d, 'a') == 2

@test collect(keys(d)) == collect('b':'z')
@test collect(values(d)) == collect(2:26)
if VERSION >= v"0.4.0-dev+980"
    @test collect(d) == [Pair(a,i) for (a,i) in zip('b':'z', 2:26)]
else
    @test collect(d) == [(a,i) for (a,i) in zip('b':'z', 2:26)]
end

# Test for #60

od60 = OrderedDict{Int,Int}()
od60[1] = 2

ranges = Range[2:5,6:9,10:13]
for range in ranges
    for i = range
        od60[i] = i+1
    end
    for i = range
        delete!( od60, i )
    end
end
od60[14]=15

@test od60[14] == 15

##############################
# Copied and modified from Base/test/dict.jl

# OrderedDict
h = OrderedDict{Int,Int}()
for i=1:10000
    h[i] = i+1
end

if VERSION >= v"0.4.0-dev+980"
    @test collect(h) == [Pair(x,y) for (x,y) in zip(1:10000, 2:10001)]
else
    @test collect(h) == collect(zip(1:10000, 2:10001))
end

for i=1:2:10000
    delete!(h, i)
end
for i=1:2:10000
    h[i] = i+1
end

for i=1:10000
    @test h[i]==i+1
end

for i=1:10000
    delete!(h, i)
end
@test isempty(h)

h[77] = 100
@test h[77]==100

for i=1:10000
    h[i] = i+1
end

for i=1:2:10000
    delete!(h, i)
end

for i=10001:20000
    h[i] = i+1
end

for i=2:2:10000
    @test h[i]==i+1
end

for i=10000:20000
    @test h[i]==i+1
end

h = OrderedDict{Any,Any}([("a", 3)])
@test h["a"] == 3
h["a","b"] = 4
@test h["a","b"] == h[("a","b")] == 4
h["a","b","c"] = 4
@test h["a","b","c"] == h[("a","b","c")] == 4

let
    z = OrderedDict()
    get_KeyError = false
    try
        z["a"]
    catch _e123_
        get_KeyError = isa(_e123_,KeyError)
    end
    @test get_KeyError
end

_d = OrderedDict([("a", 0)])
@test isa([k for k in filter(x->length(x)==1, collect(keys(_d)))], Vector{Any})

let
    ## TODO: this should work, but inference seems to be working incorrectly
    #d = OrderedDict(((1, 2), (3, 4)))
    d = OrderedDict([(1, 2), (3, 4)])
    @test d[1] === 2
    @test d[3] === 4
    ## TODO: @compat only rewrites Dict, not OrderedDict
    # d2 = OrderedDict(1 => 2, 3 => 4)
    # d3 = OrderedDict((1 => 2, 3 => 4))
    # @test d == d2 == d3
    # @test typeof(d) == typeof(d2) == typeof(d3) == OrderedDict{Int,Int}
    @test typeof(d) == OrderedDict{Int,Int}

    #d = OrderedDict(((1, 2), (3, "b")))
    d = OrderedDict([(1, 2), (3, "b")])
    @test d[1] === 2
    @test d[3] == "b"
    # d2 = OrderedDict(1 => 2, 3 => "b")
    # d3 = OrderedDict((1 => 2, 3 => "b"))
    # @test d == d2 == d3
    # @test typeof(d) == typeof(d2) == typeof(d3) == OrderedDict{Int,Any}
    @test typeof(d) == OrderedDict{Int,Any}

    #d = OrderedDict(((1, 2), ("a", 4)))
    d = OrderedDict([(1, 2), ("a", 4)])
    @test d[1] === 2
    @test d["a"] === 4
    # d2 = OrderedDict(1 => 2, "a" => 4)
    # d3 = OrderedDict((1 => 2, "a" => 4))
    # @test d == d2 == d3
    # @test typeof(d) == typeof(d2) == typeof(d3) == OrderedDict{Any,Int}
    @test typeof(d) == OrderedDict{Any,Int}

    #d = OrderedDict(((1, 2), ("a", "b")))
    d = OrderedDict([(1, 2), ("a", "b")])
    @test d[1] === 2
    @test d["a"] == "b"
    # d2 = OrderedDict(1 => 2, "a" => "b")
    # d3 = OrderedDict((1 => 2, "a" => "b"))
    # @test d == d2 == d3
    # @test typeof(d) == typeof(d2) == typeof(d3) == OrderedDict{Any,Any}
    @test typeof(d) == OrderedDict{Any,Any}
end

# TODO: this is a BoundsError on v0.3, ArgumentError on v0.4
#@test_throws ArgumentError first(OrderedDict())

if VERSION >= v"0.4.0-dev+980"
    @test first(OrderedDict([(:f, 2)])) == Pair(:f,2)
else
    @test first(OrderedDict([(:f, 2)])) == (:f,2)
end

# issue #1821
let
    d = OrderedDict{UTF8String, Vector{Int}}()
    d["a"] = [1, 2]
    @test_throws MethodError d["b"] = 1
    @test isa(repr(d), AbstractString)  # check that printable without error
end

# issue #2344
let
    local bar
    bestkey(d, key) = key
    bestkey{K<:AbstractString,V}(d::Associative{K,V}, key) = string(key)
    bar(x) = bestkey(x, :y)
    @test bar(OrderedDict([(:x, [1,2,5])])) == :y
    @test bar(OrderedDict([("x", [1,2,5])])) == "y"
end

@test  isequal(OrderedDict(), OrderedDict())
@test  isequal(OrderedDict([(1, 1)]), OrderedDict([(1, 1)]))
@test !isequal(OrderedDict([(1, 1)]), OrderedDict())
@test !isequal(OrderedDict([(1, 1)]), OrderedDict([(1, 2)]))
@test !isequal(OrderedDict([(1, 1)]), OrderedDict([(2, 1)]))

# Generate some data to populate dicts to be compared
data_in = [ (rand(1:1000), randstring(2)) for _ in 1:1001 ]

# Populate the first dict
d1 = OrderedDict{Int, ASCIIString}()
for (k,v) in data_in
    d1[k] = v
end
data_in = collect(d1)
# shuffle the data
for i in 1:length(data_in)
    j = rand(1:length(data_in))
    data_in[i], data_in[j] = data_in[j], data_in[i]
end
# Inserting data in different (shuffled) order should result in
# equivalent dict.
d2 = OrderedDict{Int, AbstractString}()
for (k,v) in data_in
    d2[k] = v
end

@test  isequal(d1, d2)
d3 = copy(d2)
d4 = copy(d2)
# Removing an item gives different dict
delete!(d1, data_in[rand(1:length(data_in))][1])
@test !isequal(d1, d2)
# Changing a value gives different dict
d3[data_in[rand(1:length(data_in))][1]] = randstring(3)
!isequal(d1, d3)
# Adding a pair gives different dict
d4[1001] = randstring(3)
@test !isequal(d1, d4)

@test isequal(OrderedDict(), sizehint(OrderedDict(),96))

# Here is what currently happens when dictionaries of different types
# are compared. This is not necessarily desirable. These tests are
# descriptive rather than proscriptive.
@test !isequal(OrderedDict([(1, 2)]), OrderedDict([("dog", "bone")]))
@test isequal(OrderedDict{Int,Int}(), OrderedDict{AbstractString,AbstractString}())

# get! (get with default values assigned to the given location)

## TODO: get! not implemented for OrderedDict
# let f(x) = x^2, d = OrderedDict((8, 19))

#     @test get!(d, 8, 5) == 19
#     @test get!(d, 19, 2) == 2

#     @test get!(d, 42) do  # d is updated with f(2)
#         f(2)
#     end == 4

#     @test get!(d, 42) do  # d is not updated
#         f(200)
#     end == 4

#     @test get(d, 13) do   # d is not updated
#         f(4)
#     end == 16

#     @test d == OrderedDict((8, 19), (19, 2), (42, 4))
# end

# issue #5886
d5886 = OrderedDict()
for k5886 in 1:11
   d5886[k5886] = 1
end
for k5886 in keys(d5886)
   # undefined ref if not fixed
   d5886[k5886] += 1
end

# issue #8877
## TODO: merge not implemented for OrderedDict
# let
#     a = OrderedDict("foo"  => 0.0, "bar" => 42.0)
#     b = OrderedDict("フー" => 17, "バー" => 4711)
#     @test is(typeof(merge(a, b)), OrderedDict{UTF8String,Float64})
# end

# issue 9295
let
    d = OrderedDict()
    @test is(push!(d, ('a', 1)), d)
    @test d['a'] == 1
    @test is(push!(d, ('b', 2), ('c', 3)), d)
    @test d['b'] == 2
    @test d['c'] == 3
    @test is(push!(d, ('d', 4), ('e', 5), ('f', 6)), d)
    @test d['d'] == 4
    @test d['e'] == 5
    @test d['f'] == 6
    @test length(d) == 6
end


text:
using Documenter, Ekztazy

makedocs(;
    modules=[Ekztazy],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
        "Types" => "types.md",
        "Client" => "client.md",
        "REST API" => "rest.md",
        "Helpers" => "helpers.md",
        "Events" => "events.md"
    ],
    repo="https://github.com/Humans-of-Julia/Ekztazy.jl/blob/{commit}{path}#L{line}",
    sitename="Ekztazy.jl",
    authors="Xh4H <sindur.esl@gmail.com>, christopher-dG <chrisadegraaf@gmail.com> Kyando <izawa.iori.tan@gmail.com>",
)

deploydocs(;
    repo="github.com/Humans-of-Julia/Ekztazy.jl",
)


text:
"""
    calc_deviance_resids(;model, df, y)

Return deviance residuals (unstandardised) from a fitted linear model.

...
# Arguments
- `model`: Variable containing a fitted linear model (e.g. general linear model [GLM]).
- `df::DataFrame=df`: DataFrame array containing response and predictor variable(s)/covariates.
- `y`: Column name from `df` containing the response variable.
...

"""
function calc_deviance_resids(;model, df, y)

    ###################################################
    # Section #1: Extract response vector and model residuals
    ###################################################

    # Extract fitted/predicted values from model object
    pred = GLM.predict(model)

    # Extract vector containing response variable
    response_y = df[!, y]

    ###################################################
    # Section #2: Calculate deviance residuals
    #             - requires conditional statements to
    #               loop through different data distributions
    ###################################################

    # Extract statistical distribution specified in model
    distType = typeof(model).parameters[1].parameters[1].parameters[2]
    distType = string.(distType)

    # Use if-elseif-else statements to calculate deviance residuals depending
    # on which errors were specified during model specification

        # If model was specified with Poisson errors
        if distType == "Poisson{Float64}"
            sign.(response_y .- pred) .* sqrt.(GLM.devresid.(Distributions.Poisson(), response_y, pred))

        # If model was specified with negative binomial errors
        elseif distType == "NegativeBinomial{Float64}"
            sign.(response_y .- pred) .* sqrt.(GLM.devresid.(Distributions.NegativeBinomial(), response_y, pred))

        elseif distType == "Normal{Float64}"
            sign.(response_y .- pred) .* sqrt.(GLM.devresid.(Distributions.Normal(), response_y, pred))

        # If model was specified with Gamma errors
        elseif distType == "Gamma{Float64}"
            sign.(response_y .- pred) .* sqrt.(GLM.devresid.(Distributions.Gamma(), response_y, pred))

        # If model was specified with Bernoulli errors
        elseif distType == "Bernoulli{Float64}"
            sign.(response_y .- pred) .* sqrt.(GLM.devresid.(Distributions.Bernoulli(), response_y, pred))

        # If model was specified with Binomial errors
        elseif distType == "Binomial{Float64}"
            sign.(response_y .- pred) .* sqrt.(GLM.devresid.(Distributions.Binomial(), response_y, pred))

        # If model was specified with unsupported error distribution (e.g. tweedie)
        else
            println("ERROR: Unsupported error distribution specified.
                    Please refer to GLM.jl documentation for supported error distributions.")

        end
end


text:
"""
Alias for union of type `T` or `nothing`.

$(TYPEDEF)
"""
const Maybe{T} = Union{T, Nothing}

"""
Alias for Dict with Symbol as key.

$(TYPEDEF)
"""
const PnmlDict = Dict{Symbol, Any}

pid(pdict::PnmlDict)::Symbol = pdict[:id]
tag(pdict::PnmlDict)::Symbol = pdict[:tag]

#has_xml(pdict::PnmlDict) = false
xmlnode(pdict::PnmlDict) = pdict[:xml]

has_labels(pdict::PnmlDict) = haskey(pdict, :labels)
has_label(d::PnmlDict, tagvalue::Symbol) = if has_labels(d)
    has_label(labels(d), tagvalue)
else
    false
end

labels(pdict::PnmlDict) = pdict[:labels]

get_label(d::PnmlDict, tagvalue::Symbol) = has_labels(d) ? get_label(labels(d), tagvalue) : nothing
get_labels(d::PnmlDict, tagvalue::Symbol) = get_labels(labels(d), tagvalue)

"""
$(TYPEDEF)
$(TYPEDFIELDS)

Wrap `PnmlDict` holding well-formed XML.
See [`ToolInfo`](@ref) and [`PnmlLabel`](@ref).
"""
@auto_hash_equals struct AnyElement
    tag::Symbol
    dict::PnmlDict
    xml::XMLNode
end

AnyElement(p::Pair{Symbol,PnmlDict}, xml::XMLNode) = AnyElement(p.first, p.second, xml)

tag(a::AnyElement) = a.tag
dict(a::AnyElement) = a.dict
xmlnode(a::AnyElement) = a.xml

"""
$(TYPEDEF)
Labels are attached to the Petri Net Graph object subtypes. See [`PnmlObject`](@ref).
"""
abstract type AbstractLabel end

xmlnode(::AbstractLabel) = nothing

# TODO: Doc HLLabel interface in manual.
"Return `true` if label has `text` field."
has_text(l::L) where {L <: AbstractLabel} = hasproperty(l, :text) && !isnothing(l.text)
"Return `text` field"
text(l::L) where {L <: AbstractLabel} = l.text

"Return `true` if label has a `structure` field."
has_structure(l::L) where {L <: AbstractLabel} = hasproperty(l, :structure) && !isnothing(l.structure)
"Return `structure` field."
structure(l::L) where {L <: AbstractLabel} = has_structure(l) ? l.structure : nothing

has_graphics(l::L) where {L <: AbstractLabel} =
    hasproperty(l, :graphics) && !isnothing(l.graphics)
graphics(l::L) where {L <: AbstractLabel} =
    has_graphics(l) ? l.graphics : nothing

has_tools(l::L) where {L <: AbstractLabel} = has_tools(l.com)#! && !isempty(tools(l))
#    hasproperty(l, :tools) && !isnothing(l.tools)
tools(l::L) where {L <: AbstractLabel} = tools(l.com)

has_labels(l::L) where {L <: AbstractLabel} = has_labels(l.com)#! && !isempty(labels(l))
#    hasproperty(l, :labels) && !isnothing(l.labels)
labels(l::L) where {L <: AbstractLabel} = labels(l.com)

has_label(l::L, tagvalue::Symbol)  where {L <: AbstractLabel}= if has_labels(l)
    has_label(labels(l), tagvalue)
else
    false
end

"""
$(TYPEDEF)
Label that may be displayed.
It differs from an Attribute Label by possibly having a [`Graphics`](@ref) field.
"""
abstract type Annotation <: AbstractLabel end

"""
$(TYPEDEF)
$(TYPEDFIELDS)

Wrap a `PnmlDict` that may be the root of an XML-tree.

Used for labels that do not have, or we choose not to use, a dedicated parse method.
Claimed labels will have a type defined to make use of the structure
defined by the pntd schema. See [`Name`](@ref), the only label defined in pnmlcore
and [`HLLabel`](@ref) for similat treatment of "unclaimed" high-level labels.
"""
@auto_hash_equals struct PnmlLabel <: Annotation
    tag::Symbol
    dict::PnmlDict
    xml::XMLNode
end

PnmlLabel(node::XMLNode; kw...) = PnmlLabel(unclaimed_label(node; kw...), node)
PnmlLabel(p::Pair{Symbol,PnmlDict}, node::XMLNode; kw...) =
        PnmlLabel(p.first, p.second, node; kw...)

tag(label::PnmlLabel) = label.tag
dict(label::PnmlLabel) = label.dict
xmlnode(label::PnmlLabel) = label.xml

function has_label(v::Vector{PnmlLabel}, tagvalue::Symbol)
    any(label->tag(label) === tagvalue, v)
end

function get_label(v::Vector{PnmlLabel}, tagvalue::Symbol)
    getfirst(l->tag(l) === tagvalue, v)
end

function get_labels(v::Vector{PnmlLabel}, tagvalue::Symbol)
    filter(l -> tag(l) === tagvalue, v)
end

"""
$(TYPEDEF)

Annotation label that uses <text> and <structure>.
"""
abstract type HLAnnotation <: AbstractLabel end

"""
$(TYPEDEF)
Objects of a Petri Net Graph are pages, arcs, nodes.
"""
abstract type PnmlObject end

"PnmlObjects are exected to have unique pnml ids."
pid(object::PnmlObject) = object.id

has_name(o::T) where {T <: PnmlObject} = o.name !== nothing
name(o::T) where {T <: PnmlObject} = has_name(o) ? o.name.text : ""

has_labels(o::T) where {T <: PnmlObject} = has_labels(o.com)
labels(o::T) where {T <: PnmlObject} = labels(o.com)

has_label(o::T, tagvalue::Symbol) where {T <: PnmlObject} =
    if has_labels(o)
        l = labels(o)
        l !== nothing ? has_label(l, tagvalue) : false
    else
        false
    end
get_label(o::T, tagvalue::Symbol) where {T <: PnmlObject} =
    if has_labels(o)
        l = labels(o)
        l !== nothing ? get_label(l, tagvalue) : nothing
    else
        nothing
    end

has_tools(o::T) where {T <: PnmlObject} = has_tools(o.com) && !isnothing(tools(o.com))
tools(o::T) where {T <: PnmlObject} = tools(o.com)
#TODO has_tool, get_tool

"""
$(TYPEDEF)
Petri Net Graph nodes are places, transitions.
"""
abstract type PnmlNode <: PnmlObject end

xmlnode(node::PnmlNode) = node.xml

"""
For common behavior shared by [`RefPlace`](@ref), [`RefTransition`](@ref).
"""
abstract type ReferenceNode <: PnmlNode end

refid(reference::ReferenceNode) = reference.ref

"""
$(TYPEDEF)

Tool specific objects can be attached to `PnmlObject`s and `AbstractLabel`s subtypes.
"""
abstract type AbstractPnmlTool end #TODO see ToolInfo

xmlnode(tool::AbstractPnmlTool) = tool.xml


text:
function Ne_to_fpe(Ne)
    # convert electron concentration to electron plasma frequency
    return 8978.18235285877*sqrt(abs(Ne)) # модуль нужен, чтобы небыло предупреждений при построении графиков из-за взятия корня от -0.
end

function fpe_to_Ne(fpe)
# convert electron plasma frequency to electron concentration
    return 1.24057537429491*fpe^2*1e-8
end

function heaviside(t)
    if t==0
        return 0
    end
    return 0.5 * (sign(t) + 1)
 end

function N0_fun(z, Nmax, z1, z2, zmax)
    N0 = Nmax*(z-z1)*(z-z2)/(zmax-z1)/(zmax-z2)*(heaviside(z-z1)-heaviside(z-z2))
    return N0
end

function dN_max_fun(z, dN, z0, d_z)
    return -dN*exp(-(z-z0)^2/2/d_z^2)
end

function func_f(t,t_r,t_d,tau_r,tau_d)
    return (1-exp(-(t-t_r)/tau_r))*(heaviside(t-t_r)-(heaviside(t-t_d)))+(1-exp(-(t_d-t_r)/tau_r))*exp(-(t-t_d)/tau_d)*heaviside(t-t_d)
end

function N_model_fun(z, t, z1, z2, zmax, dN, z0, d_z, t_r, t_d, tau_r, tau_d)
    N0 = N0_fun(z, Nmax, z1, z2, zmax)
    dN_max = dN_max_fun(z, dN, z0, d_z)
    f = func_f(t,t_r,t_d,tau_r,tau_d)
    return N0 + dN_max*f
end

# z_axe = z1:z2
# t_axe = 0:2:14

# N_model = zeros(length(z_axe),length(t_axe))
# for i=1:length(z_axe)
#     for j=1:length(t_axe)
#         N_model[i,j] = N_model_fun(z_axe[i], t_axe[j], z1, z2, zmax, dN, z0, d_z, t_r, t_d, tau_r, tau_d)
#     end
# end

# using Plots
# plot(N_model,z_axe/1000,lw=3)
# ylims!((z0/1000-10,z0/1000+10))
# xlims!((0.85*N_pump,1.15*N_pump))


text:
input = [(parse.(Int, split(pair, ','))..., ) for pair in readlines("06.input")]

manhattan_dist(a,b) = sum(abs.(a.-b))

function findmin_noties(xs, tie_sentinel=0)
    minx = first(xs)
    mini = 0
    for (i,x) in enumerate(xs)
        mini = x > minx ? mini : x == minx ? tie_sentinel : i
        minx = min(x, minx)
    end
    mini
end

function count_map(xs, n)
    counts = zeros(Int, n)
    for x in xs
        x > 0 && (counts[x] += 1)
    end
    counts
end

function star1(input)
    xs = 1:maximum(first.(input))+50
    ys = 1:maximum(last.(input))+50
    closest = [findmin_noties(manhattan_dist((x,y), i) for i in input) for x in xs, y in ys]
    dq = union(Set(closest[1:end, 1]),
               Set(closest[1:end, end]),
               Set(closest[1, 1:end]),
               Set(closest[end, 1:end]))
    delete!(dq, 0)
    counts[collect(dq)] .= 0
    findmax(counts)
end

function star2(input)
    xs = 1:maximum(first.(input))+50
    ys = 1:maximum(last.(input))+50
    sum(sum(manhattan_dist((x,y), i) for i in input)<10000 for x in xs, y in ys)
end


text:
@testset "arithmetic.jl" begin

    # Sampling based validation

    Nsamps = 10^5

    a = interval(0,2) ∪ interval(3,4)
    b = interval(-10,-9) ∪ interval(20, 25)

    aSamps = zeros(Nsamps); bSamps = zeros(Nsamps)

    # Sample inside set
    for i = 1: Nsamps
        if rand() < 0.5
            aSamps[i] = rand() * 2
        else
            aSamps[i] = rand() + 3
        end

        if rand() < 0.5
            bSamps[i] = (rand() * 5) + 20
        else
            bSamps[i] = rand() - 10
        end
    end

    @test all(aSamps .∈ a)
    @test all(bSamps .∈ b)

    BivOps = [+, -, *, / , min, max, ^] # Error with /
    UniOps = [-, sin, cos, tan, exp, log, sqrt]

    for op in BivOps

        cInt = op(a, b)
        cSamps = op.(aSamps, bSamps)
        @test all( cSamps .∈ cInt)

    end

    for op in UniOps

        cInt = op(a)
        cSamps = op.(aSamps)
        @test all( cSamps .∈ cInt)

    end

    # With 0 ∈ b

    a = interval(0,2) ∪ interval(3,4)
    b = interval(-5,5) ∪ interval(20, 25)

    for i = 1: Nsamps
        if rand() < 0.5
            aSamps[i] = rand() * 2
        else
            aSamps[i] = rand() + 3
        end

        if rand() < 0.5
            bSamps[i] = (rand() * 5) + 20
        else
            bSamps[i] = rand() * 10  - 5
        end
    end

    @test all(aSamps .∈ a)
    @test all(bSamps .∈ b)

    BivOps = [+, -, *, / , min, max] # Error with /

    for op in BivOps

        cInt = op(a, b)
        cSamps = op.(aSamps, bSamps)
        @test all( cSamps .∈ cInt)

    end

end


text:
mutable struct Target
    node::Union{Nothing,Int}
    path::String
    extensions::Union{Nothing,Dict}
    extras::Union{Nothing,Dict}
    function Target()
        obj = new()
        obj.node = nothing
        obj.extensions = nothing
        obj.extras = nothing
        obj
    end
end

function Base.setproperty!(obj::Target, sym::Symbol, x)
    if sym === :node && x !== nothing
        x ≥ 0 || throw(ArgumentError("the index of the node to target should be ≥ 0"))
    elseif sym === :path
        x === "translation" || x === "rotation" || x === "scale" || x === "weights" ||
            throw(ArgumentError("""path should be one of: "translation", "rotation", "scale", "weights" """))
    end
    setfield!(obj, sym, x)
end

JSON3.StructType(::Type{Target}) = JSON3.Mutable()
JSON3.omitempties(::Type{Target}) = (:node, :extensions, :extras)

mutable struct Channel
    sampler::Int
    target::Target
    extensions::Union{Nothing,Dict}
    extras::Union{Nothing,Dict}
    function Channel()
        obj = new()
        obj.extensions = nothing
        obj.extras = nothing
        obj
    end
end

function Base.setproperty!(obj::Channel, sym::Symbol, x)
    if sym === :sampler
        x ≥ 0 || throw(ArgumentError("the index of a sampler used to compute the value for the target should be ≥ 0"))
    end
    setfield!(obj, sym, x)
end

JSON3.StructType(::Type{Channel}) = JSON3.Mutable()
JSON3.omitempties(::Type{Channel}) = (:extensions, :extras)

mutable struct AnimationSampler
    input::Int
    interpolation::Union{Nothing,String}
    output::Int
    extensions::Union{Nothing,Dict}
    extras::Union{Nothing,Dict}
    function AnimationSampler()
        obj = new()
        obj.interpolation = nothing
        obj.extensions = nothing
        obj.extras = nothing
        obj
    end
end

function Base.getproperty(obj::AnimationSampler, sym::Symbol)
    x = getfield(obj, sym)
    sym === :interpolation && x === nothing && return "LINEAR"
    return x
end

function Base.setproperty!(obj::AnimationSampler, sym::Symbol, x)
    if sym === :input
        x ≥ 0 || throw(ArgumentError("the index of an accessor containing keyframe input values should be ≥ 0"))
    elseif sym === :interpolation && x !== nothing
        x === "LINEAR" || x === "STEP" || x === "CUBICSPLINE" ||
            throw(ArgumentError("""interpolation should be one of: "LINEAR", "STEP", "CUBICSPLINE" """))
    elseif sym === :output
        x ≥ 0 || throw(ArgumentError("the index of an accessor containing keyframe output values should be ≥ 0"))
    end
    setfield!(obj, sym, x)
end

JSON3.StructType(::Type{AnimationSampler}) = JSON3.Mutable()
JSON3.omitempties(::Type{AnimationSampler}) = (:interpolation, :extensions, :extras)

mutable struct Animation
    channels::ZVector{Channel}
    samplers::ZVector{AnimationSampler}
    name::Union{Nothing,String}
    extensions::Union{Nothing,Dict}
    extras::Union{Nothing,Dict}
    function Animation()
        obj = new()
        obj.name = nothing
        obj.extensions = nothing
        obj.extras = nothing
        obj
    end
end

JSON3.StructType(::Type{Animation}) = JSON3.Mutable()
JSON3.omitempties(::Type{Animation}) = (:name, :extensions, :extras)


text:
using Revise

include("../src/signalRepresentations.jl")
include("../src/types.jl")

"Update buffer using audio stream every Stream and track must have the same number of channels."
function updateBufferContinuously(stream::PortAudioStream, buffer::Observable, done::Observable{Bool}; window=10ms)
  while done[] == false {
    buffer[] = read(stream, window)
    sleep_time = uconvert(s, window).val
    sleep(sleep_time)
  }
  end
end

convertToMonoThroughAvg(sample_buffer) = (sample_buffer[:, 1] + sample_buffer[:, 2]) / 2

done = Observable(false)
stream = PortAudioStream("PulseAudio JACK Sink", 2, 0)
input_buffer = Observable(read(stream, 5ms))
mono_input_buffer = map(convertToMonoThroughAvg, input_buffer)
representations = SignalRepresentations(mono_input_buffer)

@async updateBufferContinuously(stream, input_buffer, done)

# done[] = false
# close(stream)

# s1 = slider(LinRange(1, 100, 100), raw = true, camera = campixel!, start = 10)
# buffer_size = s1[end][:value]

# scene = hbox(visualize(obs_buffer), s1)

# @async begin
#   while !isopen(scene) # wait for screen to be open
#     sleep(0.01)
#   end
#   while isopen(scene)
#     obs_buffer[] = read(stream, buffer_size[]ms)
#     sleep(buffer_size[] / 1000)
#   end

#   close(stream)
# end


text:
using PiCraft

function buildArena(pos::Tuple{Real, Real, Real}, length, width, floornumber)
    if length < 10 || width < 10 || floornumber < 1
        post("The size doesn't make sense!!")
        return
    end

    initialBuild(pos, length, width, floornumber)
    fillArena(pos, length, width, floornumber)
end

function initialBuild(pos::Tuple{Real, Real, Real}, length, width, floornumber)
    height = (floornumber + 1) * 10

    #clean the area
    setBlocks(pos .+ (0, -1, 0), pos .+ (length, height - 1, width), PiCraft.AIR)

    #build the floor
    setBlocks(pos .+ (0, -1, 0), pos .+ (length, -1, width), PiCraft.STONE)
    #build the roof
    setBlocks(pos .+ (0, height - 1, 0), pos .+ (length, height - 1, width), PiCraft.STONE)

    #4 corners are needed to build the walls
    cor1 = pos
    cor2 = pos .+ (length, 0, width)
    #and upper corners:
    cor3 = pos .+ (length, height - 1, 0)
    cor4 = pos .+ (0, height - 1, width)

    #build the walls
    setBlocks(cor1, cor3, PiCraft.GLASS)
    setBlocks(cor1, cor4, PiCraft.GLASS)
    setBlocks(cor2, cor3, PiCraft.GLASS)
    setBlocks(cor2, cor4, PiCraft.GLASS)

    #place one bedrock so the arena can be easily refilled
    setBlock(pos .+ (-1, -1, -1), PiCraft.BEDROCK)
end

function fillArena(pos::Tuple{Real, Real, Real}, length, width, floornumber)
    for i = 1:floornumber
        cor1 = pos .+ (1, i * 10, 1)
        cor2 = pos .+ (length - 1, i * 10, width - 1)

        setBlocks(cor1, cor2, PiCraft.TNT)
        setBlocks(cor1 .+ (0, 1, 0), cor2 .+ (0, 1, 0), PiCraft.SAND)
    end
end


text:
### A Pluto.jl notebook ###
# v0.12.3

using Markdown
using InteractiveUtils

# ╔═╡ bd3bbd0c-056e-11eb-05c3-1dae71fbcfd3
using DataFrames

# ╔═╡ f5b05dfa-056e-11eb-3e75-49e213daff01
using CSV

# ╔═╡ fc1efae4-0577-11eb-17df-ddcf53b193cb
using DelimitedFiles

# ╔═╡ 0a043ca0-0578-11eb-151b-2d706034d2cc
using BenchmarkTools

# ╔═╡ 98ec3f42-0581-11eb-112d-973e42d69028
using Queryverse

# ╔═╡ 7e79203c-0591-11eb-183b-3d463db683a1
using Dates

# ╔═╡ 0890fdf0-0593-11eb-2617-8993a9014069
using Plots

# ╔═╡ 7e184ee0-084a-11eb-26d1-3f27e600df7d
using Flux

# ╔═╡ 9834bb80-084a-11eb-37de-412fbca3cee7
using Flux: @epochs

# ╔═╡ a375f890-084a-11eb-1ec9-ddea0b709d31
using IterTools: ncycle

# ╔═╡ b0ff871a-084a-11eb-11f3-55f601c5f63f
using Parameters: @with_kw

# ╔═╡ c02c25c2-084a-11eb-15a6-cda4fb28c236
using WebIO

# ╔═╡ 950cbb44-0562-11eb-1949-e3d9165a2436
md" # Covid-19, Analysis, Visualization, Prediction"

# ╔═╡ 42e72a4c-0563-11eb-10f8-b5db40ef30a8
md"**Author: Zhuofei, Zhou**

***"

# ╔═╡ 774be886-0563-11eb-0e13-31a4fe1d9271
md"## Introduction
### Covid-19

![Covid-19](https://ss0.bdstatic.com/94oJfD_bAAcT8t7mm9GUKT-xh_/timg?image&quality=100&size=b4000_4000&sec=1601719355&di=00632f6975b3b3d60acdef92bcd8fd25&src=http://www.yxppt.com/wp-content/uploads/2020/02/24/334627/yxppt-2020-02-24_12-48-48_812888.jpg)

Novel Coronavirus 2019, on 12 January 2020, WHO officially named it 2019-NCOV. Coronaviruses are a large family of viruses known to cause colds and more serious illnesses such as Middle East Respiratory syndrome (MERS) and severe acute respiratory syndrome (SARS). Novel Coronavirus is a novel coronavirus strain that has never been found in humans before. More information is available.[百度百科](https://baike.baidu.com/item/2019新型冠状病毒/24267858?fr=aladdin)

### data

Data from COVID-19 Data Repository by the Center for Systems Science and Engineering (CSSE) at Johns Hopkins University.[link](https://github.com/CSSEGISandData/COVID-19)
- time series covid19 confirmed global.csv.([github link](https://github.com/CSSEGISandData/COVID-19/blob/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_confirmed_global.csv)).
- time series covid19 deaths global.csv.([github link](https://github.com/CSSEGISandData/COVID-19/blob/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_deaths_global.csv))
- time serise covid19 recoverd global.csv([github link](https://raw.githubusercontent.com/CSSEGISandData/COVID-19/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_recovered_global.csv))

- cases contry.csv.([github link](https://raw.githubusercontent.com/CSSEGISandData/COVID-19/web-data/data/cases_country.csv))

- cases time.csv. ([github link](https://raw.githubusercontent.com/CSSEGISandData/COVID-19/web-data/data/cases_time.csv))

download here:

"

# ╔═╡ 27e37c60-0582-11eb-1d76-05359f6465a4
md"
***
import some packages in julia:"

# ╔═╡ 885c8cee-0602-11eb-078c-bf0c90d22d00
pgfplotsx()

# ╔═╡ 5056ae0c-084a-11eb-2940-297a3406887e
import Plots.plot

# ╔═╡ 5a93b1bc-084a-11eb-1e37-d3665d2c30a8
import PlotlyJS

# ╔═╡ ef325640-059b-11eb-02e1-85c67610f89b
md"***"

# ╔═╡ 817a46da-056e-11eb-0d48-97ee88055ce0
comfirmed = load("../datasets/time_series_covid19_comfirmed_global.csv") |> DataFrame

# ╔═╡ 6ad7e4be-058c-11eb-3fac-35f5e9f14513
dt, Header = readdlm("../datasets/time_series_covid19_comfirmed_global.csv", ',', header=true)

# ╔═╡ f88187b8-0594-11eb-3181-41bdee367c4a
dt_deaths, Header_deaths = readdlm("../datasets/time_series_covid19_deaths_global.csv", ',', header=true)

# ╔═╡ 1be24c26-05fc-11eb-15e8-81d8904a38ab
dt_recover,Header_recov = readdlm("../datasets/time_series_covid19_recovered_global.csv",',',header=true)

# ╔═╡ 15ed28e6-059c-11eb-190c-2512c8cc9d8e
md"## Analysis
### Analysis Global Data
Now, we just use the data to rebuild a new dataframe to store the Global  data.
I create a new dataframe named as **Global_num** to store the number of deaths and comfirmed numbers，time from 2020-01-22 to 2020-10-02.

#### Data Operate
operating as folows:"

# ╔═╡ 1037f522-0590-11eb-2a1e-318b94418794
begin
	confirmed_num = []
	deaths_num = []
	recovered_num = []
	for i in 5:259
		push!(confirmed_num, sum(dt[:, i]))
		push!(deaths_num, sum(dt_deaths[:, i]))
		push!(recovered_num, sum(dt_recover[:,i]))
	end
	dates = Date(2020, 1, 22):Day(1):Date(2020, 10, 2)
end

# ╔═╡ b8eb5ea0-0592-11eb-2240-ff69390ac636
Global_num = DataFrame(Date=dates, confirmed_num=confirmed_num, deaths_num=deaths_num, recovered_num=recovered_num)

# ╔═╡ a1753130-0601-11eb-1371-ab5cf2290e5b
begin
	n, = size(confirmed_num)
	new_confirmed = []
	new_death = []
	new_recovered = []
	death_rate = deaths_num ./ confirmed_num
	recover_rate = recovered_num ./ confirmed_num
	for i in 2:n
		push!(new_confirmed, confirmed_num[i]-confirmed_num[i-1])
		push!(new_death, deaths_num[i] - deaths_num[i-1])
		push!(new_recovered, recovered_num[i]-recovered_num[i-1])
	end
	Global_num.death_rate = death_rate
	Global_num.recover_rate = recover_rate
end

# ╔═╡ 40d40f2e-059d-11eb-05a2-47b885c74f74
md"
#### Plot show
Then i plot the number of comfirmed cases and deaths cases over time.

As we can see, the number of comfirmed cases continues to increase exponentially. But the number of deaths cases seems to increase linearly.
"

# ╔═╡ 7cb20d98-0593-11eb-1ad9-37117ba3e5a0
p1 = plot(dates, confirmed_num, legend=:none, ylabel="Numbers", xlabel="Date", title="Global comfirmed numbers")

# ╔═╡ 56f43390-0595-11eb-116d-e10adb862fea
p2 = plot(dates, deaths_num, legend=:none, linecolor=:red, ylabel="Numbers", xlabel="Date", title="Global Deaths numbers")

# ╔═╡ 990e195a-05fc-11eb-37c9-d15175e25493
p_re = plot(dates, recovered_num, legend=:none, linecolor=:green, ylabel="Numbers", xlabel="Date", title="Global recovered numbers")

# ╔═╡ 6b000e4a-0595-11eb-1d6c-6fb54a74b93f
begin
	p3 = plot(dates, [confirmed_num deaths_num recovered_num], label=["Confirmed" "Deaths" "recovered"], legend=:topleft)
	l = @layout([a [b;c;d]])
	plot(p3, p1, p2,p_re, layout=l,xticks = [Date(2020, 1, 22), Date(2020,10,2)],size=(700,400), ylabel="Numbers", xlabel="Date", titlefontsize=10, xaxis=(font(8)))
end

# ╔═╡ f7e7f280-05f1-11eb-1d21-71674a8adbc0
md"Now, start calculate the New cases every day，and the deaths rate."

# ╔═╡ 8f2f1d20-05f3-11eb-25b0-1784220d2499
p4 = plot(Date(2020, 1, 23):Day(1):Date(2020, 10, 2),new_confirmed,shape=:circle, legend=:topleft, label="new case", xlabel="Date", ylabel="Numbers", title="new cases every day")

# ╔═╡ c6d41bda-05f9-11eb-347b-6f367c5b70fb
p5 = plot(Date(2020,1,23):Day(1):Date(2020,10,2), new_death,linecolor=:red,shape=:circle, legend=:topleft, label="new death", xlabel="Date", ylabel="Numbers", title="New death every day")

# ╔═╡ 116a25a6-05fd-11eb-2969-c37e38fab7ae
p6 = plot(Date(2020,1,23):Day(1):Date(2020,10,2), new_recovered,shape=:circle,linecolor=:green, legend=:topleft, label="new recovered", xlabel="Date", ylabel="Numbers", title="New recoverd every day")

# ╔═╡ 0e76b13c-05fa-11eb-3ec2-d76b973dbf5a
Global_num[Global_num.death_rate .== maximum(Global_num.death_rate), :]

# ╔═╡ 1ac7c3b8-05f5-11eb-048c-47188f25dd78
begin
	plot(Date(2020, 1, 22):Day(1):Date(2020, 10, 2),death_rate, linecolor=:red, title="Death rate", xlabel="date", ylabel="Numbers", legend=:none, ylims=(0,0.08))
	plot!([Date(2020,1,22),Date(2020,4,29),Date(2020,4,29)], [maximum(Global_num.death_rate),maximum(Global_num.death_rate), 0],
    line=:dot, linecolor=:black, annotations=[(Date(2020,4,29), 0.075, text("max death rate: 0.0723", 10))])
end

# ╔═╡ ba3e86c2-05fd-11eb-1431-354e2a567533
plot(Date(2020, 1, 22):Day(1):Date(2020, 10, 2),[death_rate recover_rate],xlabel="date", ylabel="Numbers", legend=:outerright, linecolor=[:red :green], title="death rate and recover rate", label=["death" "recover"])

# ╔═╡ b0ece2e4-05f8-11eb-30e4-39e726ac635c
md"the maximum death rate in 2020-4-29

So, the new comfirmed cases continues to rise, the death rate is already falling."

# ╔═╡ a41cc41c-08a5-11eb-04ed-43bbcf5b8bc8
md"
#### Prediction
**Global Trend:**

It is useful to understand the global trend of an increase in the number of cases over time. There is always a pattern in any data, but the concern is how strongly data follows a pattern. COVID-19 spreads exponentially.

i think of some ways to estimate the curve:

- Numerical Analysis: ❎
- linear regression, after transmission: ✅
- neural network: ✅

So, we first focus on comfirmed case:

do $y'=\sqrt{y}$ when i use to transition. maybe will get a linea function.:
"

# ╔═╡ 074a5f34-08a8-11eb-2e6f-3fe69c3be774
begin
	#transition
	x = 0:n-1
	y = Global_num.confirmed_num
	y₁ = sqrt.(y)
	X = [x ones(n)]
	β = inv(X'*X)*X'*y₁
	ŷ = (X*β).^2
	plot(dates,Global_num.confirmed_num, label="real")
	plot!(dates, ŷ, label="predict", annotations=[(Date(2020,4,22), 30000000, text("predict function: y=(24.77t - 584.29)²", 10))])
end

# ╔═╡ 690f7ea2-08a8-11eb-3a77-c32e93b0b950
md"""
Other ways: _**Neural Network**_

`import Flux` to do meachine learning.

the all code here:
```julia
using Flux, Statistics
using Flux.Data: DataLoader
using Flux: throttle
using Parameters: @with_kw
using DelimitedFiles
using IterTools: ncycle
using Dates
using DataFrames
using CSV
using Plots

@with_kw mutable struct Args
    η::Float64 = 0.001
    batchsize::Int = 1
    epochs::Int = 1000
end
cd(@__DIR__)
pwd()
#read data:
dt, Header = readdlm("time_series_covid19_comfirmed_global.csv", ',', header=true)
dt_deaths, Header_deaths = readdlm("time_series_covid19_deaths_global.csv", ',', header=true)
dt_recover,Header_recov = readdlm("time_series_covid19_recovered_global.csv",',',header=true)

#create a new datafram to store datum
comfirmed_num = []
deaths_num = []
recovered_num = []
for i in 5:259
    push!(comfirmed_num, sum(dt[:, i]))
    push!(deaths_num, sum(dt_deaths[:, i]))
    push!(recovered_num, sum(dt_recover[:,i]))
end
dates = Date(2020, 1, 22):Day(1):Date(2020, 10, 2)
Global_num = DataFrame(Date=dates,
    comfirmed_num=comfirmed_num,
    deaths_num=deaths_num,
    recovered_num=recovered_num)

#read train_data
n, = size(comfirmed_num)
x = 0:n-1
y = Global_num.comfirmed_num
args = Args()
train_data = DataLoader((Array(x), Float64.(y)), batchsize=args.batchsize)

# define leaky relu
Lelu(x, α=100) = (x ≥ 0 ? x : x/α)

#define Model
m = Chain(
    Dense(1, 40, Lelu),
    Dense(40, 40, Lelu),
    Dense(40, 40, Lelu),
    Dense(40, 1, Lelu)
)

#define loss function
loss(x, y) = Flux.mse(m(x), y)

#define parameters
ps = Flux.params(m)

#define Opt
opt = ADAM(args.η)

#train model
Flux.train!(loss, ps, ncycle(train_data, args.epochs), opt)

#visualize
flux_y = []
for i in Array(x)
    push!(flux_y, Array(m([i]))[1])
end
flux_y = Float64.(flux_y)
plot(dates, [flux_y Float64.(y)],
    label=["predict" "real"], xlabel="date", ylabel="Numbers", size=(900, 600))
```

"""

# ╔═╡ 87b583e0-08aa-11eb-0cd9-396999de8d5b
#write a function
function predict_curve(x, y, η::Float64=0.001, epochs::Int=1000, batchsize::Int=1)
    dates = Date(2020, 1, 22):Day(1):Date(2020, 10, 2)
	#prepare data
    train_data = Flux.Data.DataLoader((x, y), batchsize=batchsize)
	#define leaky relu
    Lelu(x, α=100) = (x ≥ 0 ? x : x/α)

	#define model
    m = Chain(
    Dense(1, 40, Lelu),
    Dense(40, 40, Lelu),
    Dense(40, 40, Lelu),
    Dense(40, 1, Lelu))

	#define mse
    loss(x, y) = Flux.mse(m(x), y)

	#parameters
    ps = Flux.params(m)

    opt = ADAM(η)
    Flux.train!(loss, ps, ncycle(train_data, epochs), opt)
    flux_y = []
    for i in Array(x)
        push!(flux_y, Array(m([i]))[1])
    end
    plot(dates, [Float64.(flux_y) Float64.(y)],label=["predict" "real"], xlabel="date", ylabel="Numbers", legend=:topleft)
end

# ╔═╡ 18d4a4fa-08ab-11eb-2912-79d0c89add25
predict_curve(Array(0:254), Float64.(Global_num.confirmed_num))

# ╔═╡ a1cc6374-08ab-11eb-2b6e-e1e9d6848603
predict_curve(Array(0:254), Float64.(Global_num.deaths_num))

# ╔═╡ f499186a-08ab-11eb-271d-277d68f0c0da
predict_curve(Array(0:254), Float64.(Global_num.recovered_num))

# ╔═╡ 30d23aa2-08ad-11eb-174d-bdaab2b5bb8d
md"
### Analysis country data

consider analyzing the situation in different countries.
"

# ╔═╡ ccaa2948-08b3-11eb-1848-c9006ef99374
dt_country = CSV.read("../datasets/cases_country.csv")

# ╔═╡ 597c6202-08ad-11eb-241f-75f2794fe332
sort(dt_country, [:Confirmed], rev=true)

# ╔═╡ a4862d54-08b3-11eb-1887-ab22f0db7600
dt_country

# ╔═╡ Cell order:
# ╟─950cbb44-0562-11eb-1949-e3d9165a2436
# ╟─42e72a4c-0563-11eb-10f8-b5db40ef30a8
# ╟─774be886-0563-11eb-0e13-31a4fe1d9271
# ╟─27e37c60-0582-11eb-1d76-05359f6465a4
# ╠═bd3bbd0c-056e-11eb-05c3-1dae71fbcfd3
# ╠═f5b05dfa-056e-11eb-3e75-49e213daff01
# ╠═fc1efae4-0577-11eb-17df-ddcf53b193cb
# ╠═0a043ca0-0578-11eb-151b-2d706034d2cc
# ╠═98ec3f42-0581-11eb-112d-973e42d69028
# ╠═7e79203c-0591-11eb-183b-3d463db683a1
# ╠═0890fdf0-0593-11eb-2617-8993a9014069
# ╠═885c8cee-0602-11eb-078c-bf0c90d22d00
# ╠═5056ae0c-084a-11eb-2940-297a3406887e
# ╠═5a93b1bc-084a-11eb-1e37-d3665d2c30a8
# ╠═7e184ee0-084a-11eb-26d1-3f27e600df7d
# ╠═9834bb80-084a-11eb-37de-412fbca3cee7
# ╠═a375f890-084a-11eb-1ec9-ddea0b709d31
# ╠═b0ff871a-084a-11eb-11f3-55f601c5f63f
# ╠═c02c25c2-084a-11eb-15a6-cda4fb28c236
# ╟─ef325640-059b-11eb-02e1-85c67610f89b
# ╠═817a46da-056e-11eb-0d48-97ee88055ce0
# ╠═6ad7e4be-058c-11eb-3fac-35f5e9f14513
# ╠═f88187b8-0594-11eb-3181-41bdee367c4a
# ╠═1be24c26-05fc-11eb-15e8-81d8904a38ab
# ╟─15ed28e6-059c-11eb-190c-2512c8cc9d8e
# ╠═1037f522-0590-11eb-2a1e-318b94418794
# ╠═b8eb5ea0-0592-11eb-2240-ff69390ac636
# ╠═a1753130-0601-11eb-1371-ab5cf2290e5b
# ╟─40d40f2e-059d-11eb-05a2-47b885c74f74
# ╠═7cb20d98-0593-11eb-1ad9-37117ba3e5a0
# ╠═56f43390-0595-11eb-116d-e10adb862fea
# ╟─990e195a-05fc-11eb-37c9-d15175e25493
# ╠═6b000e4a-0595-11eb-1d6c-6fb54a74b93f
# ╟─f7e7f280-05f1-11eb-1d21-71674a8adbc0
# ╠═8f2f1d20-05f3-11eb-25b0-1784220d2499
# ╠═c6d41bda-05f9-11eb-347b-6f367c5b70fb
# ╠═116a25a6-05fd-11eb-2969-c37e38fab7ae
# ╠═0e76b13c-05fa-11eb-3ec2-d76b973dbf5a
# ╠═1ac7c3b8-05f5-11eb-048c-47188f25dd78
# ╟─ba3e86c2-05fd-11eb-1431-354e2a567533
# ╟─b0ece2e4-05f8-11eb-30e4-39e726ac635c
# ╟─a41cc41c-08a5-11eb-04ed-43bbcf5b8bc8
# ╠═074a5f34-08a8-11eb-2e6f-3fe69c3be774
# ╟─690f7ea2-08a8-11eb-3a77-c32e93b0b950
# ╟─87b583e0-08aa-11eb-0cd9-396999de8d5b
# ╠═18d4a4fa-08ab-11eb-2912-79d0c89add25
# ╠═a1cc6374-08ab-11eb-2b6e-e1e9d6848603
# ╠═f499186a-08ab-11eb-271d-277d68f0c0da
# ╟─30d23aa2-08ad-11eb-174d-bdaab2b5bb8d
# ╠═ccaa2948-08b3-11eb-1848-c9006ef99374
# ╠═597c6202-08ad-11eb-241f-75f2794fe332
# ╠═a4862d54-08b3-11eb-1887-ab22f0db7600

max_stars_repo_name: ZhouZhuofei/juliaDL

max_issues_repo_name: ZhouZhuofei/juliaDL

max_forks_repo_name: ZhouZhuofei/juliaDL


text:
module test_primitive_cell

using Test, Xtals, XtalsPyTools

@testset "pymatgen dep" begin
    pymatgen = rc[:pymatgen]
    rc[:pymatgen] = nothing
    @test_throws ErrorException primitive_cell(Crystal("IRMOF-1.cif"))
    rc[:pymatgen] = pymatgen
end

@testset "primitive_cell" begin
    xtal = Crystal("str_m1_o3_o20_pcu_sym.22.cif")
    prim = primitive_cell(xtal)
    @test prim.atoms.n == 996
    @test isnothing(assert_P1_symmetry(prim))
end

end

max_stars_repo_name: SimonEnsemble/XtalsPyTools

max_issues_repo_name: SimonEnsemble/XtalsPyTools


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule diffutils_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("diffutils")
JLLWrappers.@generate_main_file("diffutils", UUID("11f5bf75-d780-529b-b59b-a3a04c974070"))
end  # module diffutils_jll


text:
function FGradDiv2Vec!(F,v1CG,v2CG,CG,Param)
OP=CG.OrdPoly+1;
NF=Param.Grid.NumFaces;
nz=Param.Grid.nz;

vCon = Param.CacheC1
DvCon = Param.CacheC2
vC1 = Param.CacheC3
D1vC1 = vCon
D2vC1 = DvCon
JC = Param.cache.JC

vCon .= v1CG.*Param.dXdxIC11 .+ v2CG.*Param.dXdxIC12;
mul!(reshape(vC1,OP,OP*NF*nz),CG.DS,reshape(vCon,OP,OP*nz*NF))
vCon .= v1CG.*Param.dXdxIC21 .+ v2CG.*Param.dXdxIC22;
mul!(reshape(PermutedDimsArray(DvCon,(2,1,3,4)),OP,OP*NF*nz),CG.DS,reshape(PermutedDimsArray(vCon,(2,1,3,4)),OP,OP*nz*NF))
vC1 .= vC1 .+ DvCon

mul!(reshape(D1vC1,OP,OP*NF*nz),CG.DW,reshape(vC1,OP,OP*nz*NF))
mul!(reshape(PermutedDimsArray(D2vC1,(2,1,3,4)),OP,OP*NF*nz),CG.DW,reshape(PermutedDimsArray(vC1,(2,1,3,4)),OP,OP*nz*NF))

@views F[:,:,:,:,1] .-= Param.HyperDGrad .* (Param.dXdxIC11.*D1vC1 .+ Param.dXdxIC21.*D2vC1)./JC;
@views F[:,:,:,:,2] .-= Param.HyperDGrad .* (Param.dXdxIC12.*D1vC1 .+ Param.dXdxIC22.*D2vC1)./JC;

end

function FGradDiv2Vec(v1CG,v2CG,CG,Param)
OP=CG.OrdPoly+1;
NF=Param.Grid.NumFaces;
nz=Param.Grid.nz;
dXdxIC = Param.cache.dXdxIC
JC = Param.cache.JC
vC1=reshape(
  CG.DS*(reshape(v1CG.*dXdxIC[:,:,:,:,1,1] +
  v2CG.*dXdxIC[:,:,:,:,1,2]
  ,OP,OP*NF*nz))
  ,OP,OP,NF,nz) +
  permute(
  reshape(
  CG.DS*reshape(
  permute(
   v1CG.*dXdxIC[:,:,:,:,2,1] +
   v2CG.*dXdxIC[:,:,:,:,2,2]
  ,[2 1 3 4])
  ,OP,OP*NF*nz)
  ,OP,OP,NF,nz)
  ,[2 1 3 4]);

D1cCG=reshape(
  CG.DW*reshape(vC1,OP,OP*NF*nz)
  ,OP,OP,NF,nz);
D2cCG=permute(reshape(
  CG.DW*reshape(
  permute(
  reshape(vC1,OP,OP,NF,nz)
  ,[2 1 3 4])
  ,OP,OP*NF*nz)
  ,OP,OP,NF,nz)
  ,[2 1 3 4]);

gradCG=zeros(OP,OP,NF,nz,2);
gradCG[:,:,:,:,1]=(dXdxIC[:,:,:,:,1,1].*D1cCG +
  dXdxIC[:,:,:,:,2,1].*D2cCG)./JC;
gradCG[:,:,:,:,2]=(dXdxIC[:,:,:,:,1,2].*D1cCG +
  dXdxIC[:,:,:,:,2,2].*D2cCG)./JC;
return gradCG
end


text:
"""
"""
abstract type AbstractGF end

"""
"""
struct TimeOrderedGF <: AbstractGF
    L   # Lesser
    G   # Greater
    R   # Retarded
    ts  # time-grid

    TimeOrderedGF(L, G, ts) = new(L, G, LowerTriangular(G - L), ts)
end

Base.:*(a::Number, g::TimeOrderedGF) = TimeOrderedGF(a * lesser(g), a * greater(g), g.ts)
Base.:+(g1::TimeOrderedGF, g2::TimeOrderedGF) = TimeOrderedGF(lesser(g1) + lesser(g2), greater(g1) + greater(g2), g1.ts)

"""
"""
struct TimeOrderedConvolution <: AbstractGF
    A::AbstractGF # Left GF
    B::AbstractGF # Right GF
    dts::AbstractMatrix # Trapezoidal rule is achieved by a multiplication with this matrix
end

"""
"""
function conv(A::TimeOrderedGF, B::TimeOrderedGF;
    dts = reduce(hcat, ([calculate_weights(A.ts[1:i], ones(Int64, i-1)); zeros(length(A.ts)-i)] for i in eachindex(A.ts)))
    )
    c = TimeOrderedConvolution(A, B, dts)
    return TimeOrderedGF(lesser(c), greater(c), B.ts) # no need for lazy `TimeOrderedConvolution`
end
const ⋆ = conv # left-associative operator


text:
module GeoMakie

using Reexport
using GeometryBasics, EarCut, Colors, AbstractPlotting.MakieLayout, AbstractPlotting, ImageMagick, GeoJSON, Glob

using AbstractPlotting.MakieLayout: Formatting

using AbstractPlotting.MakieLayout.GridLayoutBase
using AbstractPlotting.MakieLayout.GridLayoutBase: Side
using GeoInterface: GeoInterface, coordinates, AbstractPolygon, AbstractMultiPolygon, features, geometry

@reexport using Colors, AbstractPlotting
@reexport using Proj4, NASAEarthObservations

import AbstractPlotting: convert_arguments, convert_attribute
using AbstractPlotting.DocStringExtensions

# fix conflicts
const AbstractGeometry = GeometryBasics.AbstractGeometry
const Point = AbstractPlotting.Point
const attributes = AbstractPlotting.attributes
const transform! = Proj4.transform!
const volume = AbstractPlotting.volume
const Mesh = GeometryBasics.Mesh
const Text = AbstractPlotting.Text

include("conversions.jl")
include("projections.jl")
include("utils.jl")
include("data.jl")

include("recipes/geoaxis.jl")
include("recipes/stock.jl")

# export recipes

export geoaxis, geoaxis!

# export projections

export LonLat, PlateCarree, Mercator, TransverseMercator, UTM, AlbersEqualArea, Robinson, WinkelTripel

# export utilities

export rec_project, rec_point, xygrid

# function __init__()
#     Base.delete_method(methods(coastlines!).ms[6])
#     Base.delete_method(methods(earth!).ms[6])
# end

end # module

## TODOS
# TODO note you need to return a tuple rom conv-arg


text:
#using ExoplanetsSysSim
#using StatsFuns
if !isdefined(:JLD) using JLD end
if !isdefined(:DataFrames) using DataFrames end
if !isdefined(:CSV) using CSV end
#import DataFrames.DataFrame, DataFrames.isna
import DataFrames.DataFrame, DataArrays.ismissing
#import ExoplanetsSysSim.StellarTable.df

## Old code for generating stellar properties # TODO: WARNING: Should eventually use version in main branch to make sure have recent improvements

## stellar_table
function setup_star_table_christiansen(sim_param::SimParam; force_reread::Bool = false)
  #global df
  df = ExoplanetsSysSim.StellarTable.df
  if haskey(sim_param,"read_stellar_catalog") && !force_reread
     return df
     #return data
  end
  stellar_catalog_filename = convert(String,joinpath(Pkg.dir("ExoplanetsSysSim"), "data", convert(String,get(sim_param,"stellar_catalog","q1_q17_dr25_stellar.csv")) ) )
  df = setup_star_table_christiansen(stellar_catalog_filename)
  add_param_fixed(sim_param,"read_stellar_catalog",true)
  ExoplanetsSysSim.StellarTable.set_star_table(df)
  return df
end

function setup_star_table_christiansen(filename::String; force_reread::Bool = false)
  #global df
  wf = WindowFunction.setup_window_function(sim_param)
  df = ExoplanetsSysSim.StellarTable.df
  if ismatch(r".jld$",filename)
  try
    data = JLD.load(filename)
    df::DataFrame = data["stellar_catalog"]
    ExoplanetsSysSim.StellarTable.set_star_table(df)
  catch
    error(string("# Failed to read stellar catalog >",filename,"< in jld format."))
  end

  else
  try
    df = CSV.read(filename,nullable=true)
  catch
    error(string("# Failed to read stellar catalog >",filename,"< in ascii format."))
  end
  end # if ismatch

  has_mass = .! (ismissing.(df[:mass]) .| ismissing.(df[:mass_err1]) .| ismissing.(df[:mass_err2]))
  has_radius = .! (ismissing.(df[:radius]) .| ismissing.(df[:radius_err1]) .| ismissing.(df[:radius_err2]))
  has_dens = .! (ismissing.(df[:dens]) .| ismissing.(df[:dens_err1]) .| ismissing.(df[:dens_err2]))
  has_rest = .! (ismissing.(df[:rrmscdpp04p5]) .| ismissing.(df[:dataspan]) .| ismissing.(df[:dutycycle]))

#=
  in_Q1Q12 = []
  for x in df[:st_quarters]
    subx = string(x)
    subx = ("0"^(17-length(subx)))*subx
    indQ = search(subx, '1')
    if ((indQ < 1) | (indQ > 12))
      push!(in_Q1Q12, false)
    else
      push!(in_Q1Q12, true)
    end
  end
=#

#Comment out following block for 'is_FGK' if using stellar catalogs with cuts already made:
#=
  is_FGK = []
  for x in 1:length(df[:teff])
    if ((df[x,:teff] > 4000.0) & (df[x,:teff] < 7000.0) & (df[x,:logg] > 4.0))
#println("logg?: ", df[x,:logg])
      push!(is_FGK, true)
    else
      push!(is_FGK, false)
    end
  end
=#

  is_usable = has_radius .& has_mass .& has_rest #.& in_Q1Q12 #.& has_dens #.& is_FGK
  println("Total number of stars in our stellar catalog: ", sum(is_usable))
  #if contains(filename,"q1_q12_christiansen.jld")
#if contains(filename,"q1_q12_christiansen")   # TODO: Ask Danely what he's trying to do here.
#is_usable = is_usable #& in_Q1Q12
#end
  # See options at: http://exoplanetarchive.ipac.caltech.edu/docs/API_keplerstellar_columns.html
  # TODO SCI DETAIL or IMPORTANT?: Read in all CDPP's, so can interpolate?
  symbols_to_keep = [ :kepid, :mass, :mass_err1, :mass_err2, :radius, :radius_err1, :radius_err2, :dens, :dens_err1, :dens_err2, :rrmscdpp04p5, :dataspan, :dutycycle ]
  delete!(df, [~(x in symbols_to_keep) for x in names(df)])    # delete columns that we won't be using anyway
  usable = find(is_usable)
  df = df[usable, symbols_to_keep]
  ExoplanetsSysSim.StellarTable.set_star_table(df)
#end
  return df
end

function test_stellar_table() # TODO: Write test for stellar catalog functions
end


text:
module TestNEOSServer

using NEOSServer
using Test

const EMAIL = "odow@users.noreply.github.com"
const SERVER = NEOSServer.Server(EMAIL)

function test_help()
    @test occursin("class NeosServer", neos_help(SERVER))
end

function test_welcome()
    @test occursin("WELCOME TO NEOS!", neos_welcome(SERVER))
end

function test_version()
    @test occursin("neos version", neos_version(SERVER))
end

function test_ping()
    @test "NeosServer is alive\n" == neos_ping(SERVER)
end

function test_printQueue()
    @test occursin("Running:\n", neos_printQueue(SERVER))
end

function test_listAllSolvers()
    @test length(neos_listAllSolvers(SERVER)) > 200
end

function test_listCategories()
    categories = neos_listCategories(SERVER)
    @test length(categories) > 20
    @test length(neos_listSolversInCategory(SERVER, categories[1])) > 0
end

function test_Optimizer_no_email()
    @test_throws UndefVarError Optimizer(solver = "Ipopt")
end

function test_Optimizer()
    io = IOBuffer()
    model = NEOSServer.Optimizer(email = EMAIL, solver = "Ipopt", stdout = io)
    NEOSServer.AmplNLWriter.MOI.optimize!(model)
    seekstart(io)
    ret = String(take!(io))
    m = match(r"Job ([0-9]+) dispatched\npassword: ([a-zA-Z]+)\n"i, ret)
    job = NEOSServer.Job(parse(Int, m[1]), m[2])
    server = model.solver_command.server
    @test neos_getCompletionCode(server, job) == "Normal"
    @test neos_getJobInfo(server, job) == Any["nco", "Ipopt", "NL", "Done"]
    @test neos_killJob(server, job) == Any["Job #$(job.number) is finished"]
    ret, offset = neos_getIntermediateResults(server, job, 0)
    @test occursin("dispatched", ret)
    @test offset > 0
    @test occursin("ipopt", neos_getFinalResults(server, job))
    @test occursin("ipopt", neos_getFinalResultsNonBlocking(server, job))
end

function runtests()
    for name in names(@__MODULE__, all = true)
        if startswith("$(name)", "test_")
            @testset "$(name)" begin
                getfield(@__MODULE__, name)()
            end
        end
    end
end

end

TestNEOSServer.runtests()


text:
# Hmm, mutates active_context
function create_term_definition(
    active_context::Dict{String, Any},
    local_context::Dict{String, Any},
    term::String,
    defined::Dict{String, Bool} = Dict{String, Bool}())
  # 4.2.2 -> 1
  if term in keys(defined)
    if defined[term]
      # Already defined, don't define again
      return nothing
    else
      # In the process of being defined
      throw("Cyclic IRI mapping")
    end
  end
  # 4.2.2 -> 2
  # Set the term as in the process of being defined
  defined[term] = false
  # 4.2.2 -> 3
  # Throw an exception if redefining a a keyword
  if term in jsonld_keywords
    throw("keyword redefinition error")
  end
  # 4.2.2 -> 4
  # Remove any existing term definition in active_context
  delete!(active_context, term)

  # 4.2.2 -> 5
  # Initialize value to a copy of term from the local_context
  value = copy(local_context[term])

  # 4.2.2 -> 6
  # If value is null or it's @id is null, set the term as null, then return
  if value == nothing || (isa(value, Dict) && value["@id"] == nothing)
    active_context[term] = nothing
    defined[term] = true
    return;
  end

  # 4.2.2 -> 7
  # Normalize value to a Dict
  simple_term = isa(value, String)
  if simple_term
    value = Dict{String, String}("@id" => value)
  end

  # 4.2.2 -> 8
  # Throw an error if value isn't a Dict by now
  if !isa(value, Dict)
    throw("Invalid term definition")
  end

  # 4.2.2 -> 9
  # Create a new term definition, definition
  definition = TermDefinition()

  # 4.2.2 -> 10
  if "@type" in keys(value)
    # 4.2.2 -> 10.1
    # Initialize type to @type
    ttype = value["@type"]
    isa(value["@type"], String) || throw("Invalid type mapping error")

    # 4.2.2 -> 10.2
    # TODO Use IRI expansion

    definition["@type"] = ttype
  end

  # 4.2.2 -> 11
  if "@reverse" in keys(value)
    # 4.2.2 -> 11.1
    ("@id" in keys(value) || "@nest" in keys(value)) && throw("Invalid reverse property")
    isa(value["@reverse"], String) || throw("Invalid IRI mapping")
    # TODO set IRI mapping of definition, not sure what the IRI mapping was
    definition[]
  end

end


text:
using PackageCompiler
# Designed to be launched from the test directory
create_sysimage(["TMLEEpistasis"];
                 sysimage_path="TMLEEpistasisSysimage.so",
                 precompile_execution_file="example_script.jl")