text:
using Ju
using BenchmarkTools
using StatsBase:sample

buffer = CircularArrayBuffer{Array{Float64, 3}}(10, (80, 80, 1));

println("\n", repeat('=', 50))
println("\npush! into CircularArrayBuffer\n")
display(@benchmark push!($buffer, $(rand(80, 80, 1))))

println("\n", repeat('=', 50))
println("\ngetindex of CircularArrayBuffer\n")
display(@benchmark $buffer[$(rand(1:length(buffer)))])

println("\n", repeat('=', 50))
println("\nview element of CircularArrayBuffer\n")
display(@benchmark view($buffer, $(rand(1:length(buffer)))))

n = 16
println("\n", repeat('=', 50))
println("\nview $n elements of CircularArrayBuffer\n")
display(@benchmark view($buffer, $(sample(1:length(buffer), n))))

n = 32
println("\n", repeat('=', 50))
println("\nview $n elements of CircularArrayBuffer\n")
display(@benchmark view($buffer, $(sample(1:length(buffer), n))))

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

# push! into CircularArrayBuffer

# BenchmarkTools.Trial:
#   memory estimate:  0 bytes
#   allocs estimate:  0
#   --------------
#   minimum time:     2.411 μs (0.00% GC)
#   median time:      2.456 μs (0.00% GC)
#   mean time:        2.464 μs (0.00% GC)
#   maximum time:     5.133 μs (0.00% GC)
#   --------------
#   samples:          10000
#   evals/sample:     9

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

# getindex of CircularArrayBuffer

# BenchmarkTools.Trial:
#   memory estimate:  50.08 KiB
#   allocs estimate:  2
#   --------------
#   minimum time:     5.540 μs (0.00% GC)
#   median time:      15.940 μs (0.00% GC)
#   mean time:        16.852 μs (27.09% GC)
#   maximum time:     8.680 ms (99.80% GC)
#   --------------
#   samples:          10000
#   evals/sample:     5

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

# view element of CircularArrayBuffer

# BenchmarkTools.Trial:
#   memory estimate:  64 bytes
#   allocs estimate:  1
#   --------------
#   minimum time:     12.913 ns (0.00% GC)
#   median time:      14.916 ns (0.00% GC)
#   mean time:        23.723 ns (30.27% GC)
#   maximum time:     45.953 μs (99.93% GC)
#   --------------
#   samples:          10000
#   evals/sample:     999

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

# view 16 elements of CircularArrayBuffer

# BenchmarkTools.Trial:
#   memory estimate:  400 bytes
#   allocs estimate:  6
#   --------------
#   minimum time:     96.007 ns (0.00% GC)
#   median time:      101.996 ns (0.00% GC)
#   mean time:        130.636 ns (17.69% GC)
#   maximum time:     45.819 μs (99.68% GC)
#   --------------
#   samples:          10000
#   evals/sample:     952

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

# view 32 elements of CircularArrayBuffer

# BenchmarkTools.Trial:
#   memory estimate:  528 bytes
#   allocs estimate:  6
#   --------------
#   minimum time:     109.849 ns (0.00% GC)
#   median time:      123.340 ns (0.00% GC)
#   mean time:        166.186 ns (18.24% GC)
#   maximum time:     51.037 μs (99.33% GC)
#   --------------
#   samples:          10000
#   evals/sample:     934


text:
# module TestSlice

using Test, Dates, Temporal

@testset "Filtering" begin
    @testset "Head/Tail" begin
        X = TS(cumsum(randn(100,5), dims=1)) + 100
        n = 5
        @test head(X, n).values == X.values[1:n,:]
        @test tail(X, n).values == X.values[end-5+1:end,:]
    end
    @testset "Missing Data Removal" begin
        @testset "Rows" begin
            X = TS(cumsum(randn(100,5), dims=1)) + 100
            X[1,1] = NaN
            @test size(dropnan(X, dim=1, fun=any), 1) == size(X, 1) - 1
            @test size(dropnan(X, dim=1, fun=all), 1) == size(X, 1)
            X[1] = NaN
            @test size(dropnan(X, dim=1, fun=any), 1) == size(X, 1) - 1
            @test size(dropnan(X, dim=1, fun=all), 1) == size(X, 1) - 1
        end
        @testset "Columns" begin
            X = TS(cumsum(randn(100,5), dims=1)) + 100
            X[1,1] = NaN
            @test size(dropnan(X, dim=2, fun=any), 2) == size(X, 2) - 1
            @test size(dropnan(X, dim=2, fun=all), 2) == size(X, 2)
            X[:,1] = NaN
            @test size(dropnan(X, dim=2, fun=any), 2) == size(X, 2) - 1
            @test size(dropnan(X, dim=2, fun=all), 2) == size(X, 2) - 1
        end
    end
    @testset "Missing Data Fills" begin
        X = TS(cumsum(randn(100), dims=1)) + 100
        X[3:5] = NaN
        @test size(dropnan(fillnan(X, :ffill)), 1) == size(X, 1)
        @test size(dropnan(fillnan(X, :bfill)), 1) == size(X, 1)
        @test size(dropnan(fillnan(X, :linear)), 1) == size(X, 1)
        fillnan!(X, :bfill)
        @test all(.!nanrows(X))
        X[3:5] = NaN
        fillnan!(X, :ffill)
        @test all(.!nanrows(X))
        X[3:5] = NaN
        fillnan!(X, :linear)
        @test all(.!nanrows(X))
    end
end

# end

ext: jl


text:
module PlotsExample

using Matte

const title = "Matte Example :: Using Plots with Matte"

function ui()
    sidebar_layout(
        side_panel(
            selector("xvar", "X variable", "['SepalLength', 'SepalWidth', 'PetalLength', 'PetalWidth', 'Species']"),
            br(),
            selector("yvar", "Y variable", "['SepalLength', 'SepalWidth', 'PetalLength', 'PetalWidth', 'Species']"),
        ),
        main_panel(
            plots_output("output_plot")
        )
    )
end

module Server

using Plots
using RDatasets
iris = dataset("datasets","iris")

function output_plot(xvar, yvar)
    if typeof(xvar) <: AbstractString && typeof(yvar) <: AbstractString
        plot(iris[!, Symbol(xvar)], iris[!, Symbol(yvar)], seriestype=:scatter)
    end
end

end

end


text:
using VegaDatasets
import IteratorInterfaceExtensions, TableTraits, TableShowUtils
using Test

@testset "VegaDatasets" begin

for filename in readdir(joinpath(@__DIR__,"..","data","data"))
    if splitext(filename)[2]==".json" || splitext(filename)[2]==".csv"
        if !in(splitdir(filename)[2], ["earthquakes.json", "graticule.json", "londonBoroughs.json", "londonTubeLines.json", "miserables.json", "sf-temps.csv", "us-10m.json", "world-110m.json"])
            vd = dataset(filename)

            @test IteratorInterfaceExtensions.isiterable(vd) == true
            @test TableTraits.isiterabletable(vd) == true

            cvd = collect(vd)

            @test isa(cvd, Array)
        end
    end
end

df = dataset("iris")

@test isa(df, VegaDatasets.VegaDataset)

df2 = dataset("airports")

@test isa(df2, VegaDatasets.VegaDataset)

@test_throws ErrorException dataset("doesnotexist")

@test_throws ErrorException dataset("7zip.png")

df3 = dataset("weather.csv")

@test isa(df3, VegaDatasets.VegaDataset)

df4 = dataset("weather.json")

@test isa(df4, VegaDatasets.VegaDataset)

@test sprint(show, dataset("lookup_groups")) == """
9x2 Vega dataset
group │ person
──────┼───────
1     │ Alan
1     │ George
1     │ Fred
2     │ Steve
2     │ Nick
2     │ Will
3     │ Cole
3     │ Rick
3     │ Tom   """

end


text:
using BEAST
using Test
using LinearAlgebra

width = 4.0
delay = 6.0
g = BEAST.creategaussian(width, delay)
G = BEAST.fouriertransform(g)

dx = width / 15
x0 = 0.0
x = x0 : dx : 3*delay
n = length(x)
y = g.(x)

Y, dω, ω0 = BEAST.fouriertransform(y, dx, x0)
ω = collect(ω0 .+ (0:n-1)*dω)
Z = G.(ω)

# using Plots
# plot(ω, real(Y))
# scatter!(ω, real(Z))
# plot!(ω, imag(Y))
# scatter!(ω, imag(Z))

@test norm(Y-Z) < 1.e-10


text:
import Literate
using Literate: script, markdown, notebook

src_dirs = [joinpath(@__DIR__, "..", "src", "models")]

trgt_dir_test = joinpath(@__DIR__, "..", "test", "models")
trgt_dir = joinpath(@__DIR__, "src", "models")
mkpath(trgt_dir)

for src_dir in src_dirs
    for dir in readdir(src_dir)
        if endswith(dir, ".jl")
            #@warn "ignoring $src_dir/$file"
        else
            src_path = abspath(joinpath(src_dir, dir, dir*".jl"))
            text = script(src_path, trgt_dir_test, credit=false)
            code = strip(read(text, String))
            mdpost(str) = replace(str, "@__CODE__" => code)
            markdown(src_path, trgt_dir, postprocess=mdpost, credit=false)
        end
    end
end

max_stars_repo_name: mforets/ReachabilityLibrary

max_issues_repo_name: mforets/ReachabilityLibrary

max_forks_repo_name: mforets/ReachabilityLibrary


text:
# functions related to F distribution

# R implementations
# For pdf and logpdf we use the Julia implementation
using .RFunctions:
    fdistcdf,
    fdistccdf,
    fdistlogcdf,
    fdistlogccdf,
    fdistinvcdf,
    fdistinvccdf,
    fdistinvlogcdf,
    fdistinvlogccdf

# Julia implementations
fdistpdf(ν1::Real, ν2::Real, x::Real) = exp(fdistlogpdf(ν1, ν2, x))

fdistlogpdf(ν1::Real, ν2::Real, x::Real) = fdistlogpdf(promote(ν1, ν2, x)...)
function fdistlogpdf(ν1::T, ν2::T, x::T) where {T<:Real}
    # we ensure that `log(x)` does not error if `x < 0`
    ν1ν2 = ν1 / ν2
    y = max(x, 0)
    val = (xlogy(ν1, ν1ν2) + xlogy(ν1 - 2, y) - xlogy(ν1 + ν2, 1 + ν1ν2 * y)) / 2 - logbeta(ν1 / 2, ν2 / 2)
    return x < 0 ? oftype(val, -Inf) : val
end


text:
# compute reduced s bubble
function compute_channel_s_reduced!(
    Λ      :: Float64,
    kernel :: Int64,
    w1     :: Int64,
    w3     :: Int64,
    r      :: Reduced_lattice,
    m      :: Mesh,
    a1     :: Action,
    a2     :: Action,
    tbuff  :: NTuple{3, Matrix{Float64}},
    temp   :: Array{Float64, 3},
    corrs  :: Array{Float64, 3},
    eval   :: Int64,
    Γ_tol  :: NTuple{2, Float64}
    )      :: Nothing

    # reset buffer
    tbuff[1] .= 0.0

    # get frequency arguments
    s, vsp = 0.0, 0.0

    if kernel == 1
        s, vsp = m.Ωs[w1], Inf
    else
        s, vsp = m.Ωs[w1], m.νs[w3]
    end

    # define integrand
    integrand!(b, v, dv) = compute_s_reduced!(Λ, b, v, dv, s, vsp, r, m, a1, temp)

    # compute integrals
    ref = Λ + 0.5 * s
    val = m.Ωs[end] + m.νs[end]
    integrate_log!((b, v, dv) -> integrand!(b, v, dv), tbuff,  2.0 * ref, 1.0 * val, eval, Γ_tol[1], Γ_tol[2], sgn = -1.0)
    integrate_lin!((b, v, dv) -> integrand!(b, v, dv), tbuff, -2.0 * ref, 0.0 * ref, eval, Γ_tol[1], Γ_tol[2])
    integrate_lin!((b, v, dv) -> integrand!(b, v, dv), tbuff,  0.0 * ref, 2.0 * ref, eval, Γ_tol[1], Γ_tol[2])
    integrate_log!((b, v, dv) -> integrand!(b, v, dv), tbuff,  2.0 * ref, 1.0 * val, eval, Γ_tol[1], Γ_tol[2])

    # correct boundaries
    integrand!(tbuff[1], -val, corrs[1, 1, w1])
    integrand!(tbuff[1],  val, corrs[2, 1, w1])

    # parse result
    for i in eachindex(a2.Γ)
        if kernel == 1
            @turbo a2.Γ[i].ch_s.q1[:, w1] .= view(tbuff[1], i, :)
        else
            @turbo a2.Γ[i].ch_s.q2_2[:, w1, w3] .= view(tbuff[1], i, :)
        end
    end

    return nothing
end

# compute self energy from SDE
function compute_Σ!(
    Λ      :: Float64,
    r      :: Reduced_lattice,
    m      :: Mesh,
    a1     :: Action,
    a2     :: Action,
    tbuffs :: Vector{NTuple{3, Matrix{Float64}}},
    temps  :: Vector{Array{Float64, 3}},
    corrs  :: Array{Float64, 3},
    eval   :: Int64,
    Γ_tol  :: NTuple{2, Float64},
    Σ_tol  :: NTuple{2, Float64}
    )      :: Nothing

    # compute boundary corrections
    compute_corrs!(Λ, m, a1, corrs, Γ_tol)

    # compute reduced s bubble in q1 and q2_2
    @sync begin
        for w1 in 1 : m.num_Ω
            # compute q1
            Threads.@spawn compute_channel_s_reduced!(Λ, 1, w1, -1, r, m, a1, a2, tbuffs[Threads.threadid()], temps[Threads.threadid()], corrs, eval, Γ_tol)

            for w3 in 1 : m.num_ν
                # compute q2_2
                Threads.@spawn compute_channel_s_reduced!(Λ, 3, w1, w3, r, m, a1, a2, tbuffs[Threads.threadid()], temps[Threads.threadid()], corrs, eval, Γ_tol)
            end
        end
    end

    # compute self energy for all frequencies
    @sync for i in 2 : length(m.σ)
        Threads.@spawn begin
            integrand = v -> compute_Σ_kernel(Λ, m.σ[i], v, r, m, a1, a2)
            a2.Σ[i]   = quadgk(integrand, -Inf, -2.0 * Λ, 0.0, 2.0 * Λ, Inf, atol = Σ_tol[1], rtol = Σ_tol[2])[1]
        end
    end

    return nothing
end


text:
include("plain.jl")
include("jld.jl")


text:
using DReal
using Base.Test
reset_ctx!()

a = Var(Bool,"a")
b = Var(Bool,"b")
add!(a → b) # \rightarrow - equivalent to implies(a,b)
a,b =  model(a,b)
@show a,b
@test !a | b


text:
function l2_squared_norm(f; tol=1e-12)
    return pcubature(x -> norm(f(x))^2, [0., 0.], [1., 1.], abstol=tol)[1]
end

function l2_squared_distance(f, g; tol=1e-12)
    return pcubature(x -> norm(f(x) - g(x))^2, [0., 0.], [1., 1.], abstol=tol)[1] #, rtol=tol, atol=tol)[1]
end

# Calculate the L2-innerproduct between two function f, g: R^2 -> R^3
function l2_inner_product(f, g; tol=1e-12)
    function euclidean_inner_product(x)
        return f(x) ⋅ g(x)
    end
    return pcubature(euclidean_inner_product, [0., 0.], [1., 1.], abstol=tol)[1] #, rtol=tol, atol=tol)[1]
end

function tangent_vector(q, r)
    scaling = √l2_squared_distance(q, r)
    return x -> (q(x) - r(x) ) / scaling
end

function unscaled_tangent_vector(q, r)
    return x -> q(x) - r(x)
end


text:
infile = open("file1", "r")
outfile = open("file2", "w")
write(outfile, read(infile, String))
close(outfile)
close(infile)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
export OrdLoc, OrdLocElem

###############################################################################
#
#   Declaration types
#   OrdLoc / OrdLocElem
#
###############################################################################

mutable struct OrdLoc{T<:nf_elem} <: Hecke.Ring
   OK::NfAbsOrd{AnticNumberField,T}
   prime::NfAbsOrdIdl{AnticNumberField,T}
   comp::Bool

   function OrdLoc{T}(OK::NfAbsOrd{AnticNumberField,T}, prime::NfAbsOrdIdl{AnticNumberField,T}, cached::Bool = true, comp::Bool = false) where {T <: nf_elem}
      if cached && haskey(OrdLocDict, (OK, prime, comp))
         return OrdLocDict[OK, prime, comp]::OrdLoc{T}
      else
         z = new(OK, prime, comp)
         if cached
            OrdLocDict[OK, prime, comp] = z
         end
         return z
      end
   end
end

function ppio(a::NfOrdIdl, b::NfOrdIdl)
   c = gcd(a, b)
   n = divexact(a, c)
   g = gcd(c, n)
   while !isone(g)
      c *= g
      n = divexact(n, g)
      g = gcd(c, n)
   end
   return c, n
end

function isin(a::nf_elem, L::OrdLoc{<:nf_elem})
  iszero(a) && return true
  n, d = integral_split(a*L.OK)
  L.comp || (!isone(gcd(d, L.prime)) && return false)
  L.comp && ppio(d, L.prime)[1] != d && return false
  return true
end

#=
s an ideal (or element = principal ideal)
  S = { x in R | v_p(x) = 0 for all p | s}
  L = { x in K | gcd(den(<x>), s) = 1} = {a/b | v_p(a) - v_p(b) >= 0 for all p|s}
    N(x) = norm(num(<x>) + s^infty)
         = prod N(p)^v_p(x)    p|s
    K -> prod R_(p)  for p | s where R_(p) = {x | v_p(x) >=0 } subset R_p = completion at p
      find q_p, r_p s.th.
      a = q_p b + r_p and v_p(r_p) < v_p(b) (either a or 0)
      any r in K s.th. v_p(r) = v_p(r_p) should be a valid remainder, q = (a-r)/b is the quotient
        c = num(<a>) + s^infty
        d = num(<b>) + s^infty
        den(c/d) has support all p where v_p(c) < v_p(d)
        r = a mod (den(c/d)^2)
          = 0 mod s/gcd(s, den(c,d))
        should work...

some more detail in divrem

  S = { x in R | v_p(x) = 0 for all p not | s}
  L = { x in K | gcd(den(<x>), s^infty) = den(<x>) } = {a/b | v_p(a) - v_p(b) >= 0 for all p not | s}
    N(x) = norm(num(<x>):s^infty)
         = prod N(p)^v_p(x)    p not | s
    in general not euclidean since not PID (S = 1 => R = R).
    if S large enough (gen class group), then PID

=#

OrdLocDict = Dict{Tuple{NfAbsOrd{AnticNumberField,nf_elem}, NfAbsOrdIdl{AnticNumberField,nf_elem}, Bool}, Hecke.Ring}()

mutable struct OrdLocElem{T<:nf_elem} <: RingElem
   data::T
   parent::OrdLoc{T}

   function OrdLocElem{T}(data::T, par::OrdLoc, checked::Bool = true) where {T <:nf_elem}
      data == zero(parent(data)) && return new{T}(data,par)
      checked && (!isin(data, par)) && error("No valid element of localization")
      return new{T}(data,par)
   end
end

###############################################################################
#
#   Unsafe operators and functions
#
###############################################################################

add!(c::OrdLocElem, a::OrdLocElem, b::OrdLocElem) = a + b

mul!(c::OrdLocElem, a::OrdLocElem, b::OrdLocElem) = a * b

addeq!(a::OrdLocElem, b::OrdLocElem) = a + b

###############################################################################
#
#   Data type and parent object methods
#
###############################################################################

elem_type(::Type{OrdLoc{T}}) where {T <: nf_elem} = OrdLocElem{T}

parent_type(::Type{OrdLocElem{T}}) where {T <: nf_elem} = OrdLoc{T}

order(L::OrdLoc{T}) where {T <: nf_elem}  = L.OK

order(a::OrdLocElem{T}) where {T <: nf_elem}  = order(parent(a))

nf(L::OrdLoc{T}) where {T <: nf_elem}  = nf(L.OK)::parent_type(T)

nf(a::OrdLocElem{T}) where {T <: nf_elem} = nf(parent(a))

parent(a::OrdLocElem{T})  where {T <: nf_elem} = a.parent

function check_parent(a::OrdLocElem{T}, b::OrdLocElem{T})  where {T <: nf_elem}
    parent(a) != parent(b) && error("Parent objects do not match")
end

###############################################################################
#
#   Basic manipulation
#
###############################################################################

data(a::OrdLocElem{T}) where {T <: nf_elem} = a.data

numerator(a::OrdLocElem{T}) where {T <: nf_elem} = numerator(data(a))

denominator(a::OrdLocElem{T}) where {T <: nf_elem} = denominator(data(a))

prime(L::OrdLoc{T}) where {T <: nf_elem} = L.prime

prime(a::OrdLocElem{T}) where {T <: nf_elem} = prime(parent(a))

zero(L::OrdLoc{T}) where {T <: nf_elem} = L(0)

one(L::OrdLoc{T}) where {T <: nf_elem} = L(1)

iszero(a::OrdLocElem{T}) where {T <: nf_elem} = iszero(data(a))

isone(a::OrdLocElem{T}) where {T <: nf_elem} = isone(data(a))

function in(x::nf_elem, L::OrdLoc)
   iszero(x) ? true :
   return isin(x, L)
end

function isunit(a::OrdLocElem{T})  where {T <: nf_elem}
   iszero(a) ? false :
    return isin(inv(a.data), parent(a))
end

deepcopy_internal(a::OrdLocElem{T}, dict::IdDict) where {T <: nf_elem} = parent(a)(deepcopy(data(a)))

###############################################################################
#
#   AbstractString I/O
#
###############################################################################

function show(io::IO, a::OrdLocElem{T}) where {T <: nf_elem}
   print(io, data(a))
end

function show(io::IO, L::OrdLoc{T}) where {T <: nf_elem}
  if L.comp
    print(io, "Localization of ", order(L), " at complement of ", prime(L))
  else
    print(io, "Localization of ", order(L), " at ", prime(L))
  end
end

needs_parentheses(x::OrdLocElem{T})  where {T <: nf_elem} = needs_parentheses(data(x))

displayed_with_minus_in_front(x::OrdLocElem{T})  where {T <: nf_elem} = displayed_with_minus_in_front(data(x))

show_minus_one(::Type{OrdLocElem{T}}) where {T <: nf_elem} = true

##############################################################################
#
#   Unary operations
#
##############################################################################

function -(a::OrdLocElem{T})  where {T <: nf_elem}
   parent(a)(-data(a))
end

###############################################################################
#
#   Binary operators
#
###############################################################################

function +(a::OrdLocElem{T}, b::OrdLocElem{T})  where {T <: nf_elem}
   check_parent(a,b)
   return parent(a)(data(a) + data(b), false)
end

function -(a::OrdLocElem{T}, b::OrdLocElem{T})  where {T <: nf_elem}
   check_parent(a,b)
   return parent(a)(data(a) - data(b), false)
end

function *(a::OrdLocElem{T}, b::OrdLocElem{T})  where {T <: nf_elem}
   check_parent(a,b)
   return parent(a)(data(a) * data(b), false)
end

###############################################################################
#
#   Comparison
#
###############################################################################

function ==(a::OrdLocElem{T}, b::OrdLocElem{T}) where {T <: nf_elem}
   check_parent(a, b)
   return data(a) == data(b)
end

##############################################################################
#
#  Inversion
#
##############################################################################

@doc Markdown.doc"""
     inv(a::OrdLocElem{T}, checked::Bool = true)  where {T <: nf_elem}
> Returns the inverse element of $a$ if $a$ is a unit.
> If 'checked = false' the invertibility of $a$ is not checked and the corresponding inverse element
> of the numberfield is returned.
"""
function inv(a::OrdLocElem{T}, checked::Bool = true)  where {T <: nf_elem}
   b = inv(data(a))
   checked && (!isin(b, parent(a))) && error("$a not invertible in given localization")
   return parent(a)(b, false)
end

##############################################################################
#
#  Exact division
#
##############################################################################

@doc Markdown.doc"""
     divides(a::OrdLocElem{T}, b::OrdLocElem{T}, checked::Bool = true) where {T <: nf_elem}
> Returns tuple (`true`,`c`) if $b$ divides $a$ where `c`*$b$ = $a$.
> If 'checked = false' the corresponding element of the numberfield is returned and it is not
> checked whether it is an element of the given localization.
"""
function divides(a::OrdLocElem{T}, b::OrdLocElem{T}, checked::Bool = true) where {T <: nf_elem}
   check_parent(a,b)

   if iszero(b)
     if iszero(a)
       return true, parent(a)()
     else
       return false, parent(a)()
     end
   end

   elem = divexact(data(a), data(b))
   if !checked
      return true, parent(a)(elem, checked)
   elseif checked && in(elem,parent(a))
      return true, parent(a)(elem)
   else
      return false, parent(a)()
   end
end

@doc Markdown.doc"""
     divexact(a::OrdLocElem{T}, b::OrdLocElem{T}, checked::Bool = true)  where {T <: nf_elem}
> Returns element 'c' of given localization s.th. `c`*$b$ = $a$ if such element exists.
> If 'checked = false' the corresponding element of the numberfield is returned and it is not
> checked whether it is an element of the given localization.
"""
function divexact(a::OrdLocElem{T}, b::OrdLocElem{T}, checked::Bool = true)  where {T <: nf_elem}
   d = divides(a, b, checked)
   d[1] ? d[2] : error("$a not divisible by $b in the given localization")
end

@doc Markdown.doc"""
     div(a::OrdLocElem{T}, b::OrdLocElem{T}, checked::Bool = true)  where {T <: nf_elem}
"""
function _make_legal(a::nf_elem, S::NfOrdIdl)
  d = denominator(a, order(S))
  n = order(S)(d*a)
  b, _ = ppio(d*order(S), S)
  bi = inv(b)
  x = bi.num.gen_two.elem_in_nf//bi.den
  # a in OrdLoc(S) implies the true den of a is coprime to S
  # return a = n/d with n, d in the order and d coprime to S
  R = order(S)
  @assert (n*x) in R
  @assert (d*x) in R
  @assert isone(R(d*x)*R + S)
  return n*x, d*x
end

function divrem(a::OrdLocElem, b::OrdLocElem, checked::Bool = true)
  L = parent(a)
  L.comp && error("ring not known to be useful euclidean")

  na, da = _make_legal(a.data, L.prime)
  nb, db = _make_legal(b.data, L.prime)

  A = na*db
  B = nb*da
  #Plan: A = qB + R
  #then a = q/da b + R/db/da
  #moving to the residue ring (B+S^infty)*S helps
  # the euclidean function is the same, and division works in there
  # Plan: get the reminder from there, then lift and compute the quotient

  C = ppio(L.OK(B)*L.OK, L.prime)[1]*L.prime
  Q, mQ = quo(L.OK, C)
  R = L(preimage(mQ, divrem(mQ(L.OK(A)), mQ(L.OK(B)))[2]))

  r = divexact(R, L(da*db))
  q = divexact(a-r, b)
  @assert iszero(r) || euclid(r) < euclid(b)
  return q, r
end
#=
#XXX: those do not seem to be there for AbsOrdQuoRingElems
function div(a::OrdLocElem{T}, b::OrdLocElem{T}, checked::Bool = true)  where {T <: nf_elem}
   return divrem(a, b)[1]
end

function rem(a::OrdLocElem{T}, b::OrdLocElem{T}, checked::Bool = true)  where {T <: nf_elem}
   return divrem(a, b)[2]
end
=#

function euclid(a::OrdLocElem)
  iszero(a) && return fmpz(0)
  L = parent(a)
  L.comp && error("ring not known to be useful euclidean")
  N, _ = integral_split(a.data * L.OK)
  N, _ = ppio(N, L.prime)
  return norm(N)
end

###############################################################################
#
#   GCD
#
###############################################################################

@doc Markdown.doc"""
    gcd(a::OrdLocElem{T}, b::OrdLocElem{T}) where {T <: nf_elem}
> Returns gcd of $a$ and $b$ in canonical representation.
"""
function gcd(a::OrdLocElem{T}, b::OrdLocElem{T}) where {T <: nf_elem}
  L = parent(a)
  L.comp && error("ring not known to be useful euclidean")
   check_parent(a,b)
   iszero(a) && return inv(canonical_unit(b)) * b
   iszero(b) && return inv(canonical_unit(a)) * a
   na, _ = _make_legal(a.data, L.prime)
   nb, _ = _make_legal(b.data, L.prime)
   g = L.OK(na) * L.OK + L.OK(nb) * L.OK + L.prime
   assure_2_normal(g)
   elem = L(g.gen_two)
   return inv(canonical_unit(elem)) * (elem)
end

###############################################################################
#
#   GCDX
#
###############################################################################

@doc Markdown.doc"""
    gcdx(a::OrdLocElem{T}, b::OrdLocElem{T}) where {T <: nf_elem}
> Returns tuple `(g,u,v)` s.th. `g` = gcd($a$,$b$) and `g` = `u` * $a$ + `v` * $b$.
"""
function gcdx(a::OrdLocElem{T}, b::OrdLocElem{T}) where {T <: nf_elem}
  L = parent(a)
  L.comp && error("ring not known to be useful euclidean")
   check_parent(a,b)
   L = parent(a)
   M = [L(1) L(0); L(0) L(1)]
   while !iszero(b)
     q,r = divrem(a, b)
     M = [L(0) L(1); L(1) L(-q)] * M
     a, b = b, r
   end
   c = inv(canonical_unit(a))
   return a*c, M[1,1]*c, M[1,2]*c
end

###############################################################################
#
#   PID
#
###############################################################################

function principal_generator(L::OrdLoc{T}, I::NfAbsOrdIdl{AnticNumberField,T}) where {T <: nf_elem}
  #possible for !L.comp due to semi local
  #theoretical for L.comp if L.prime large enough...
   valuation(L(I.gen_one)) >= valuation(L(I.gen_two)) ? L(I.gen_two) : L(I.gen_one)
end

###############################################################################
#
#   Powering
#
###############################################################################

function ^(a::OrdLocElem{T}, b::Int) where {T <: nf_elem}
   return parent(a)(data(a)^b, false)
end

###############################################################################
#
#   Random Functions
#
###############################################################################

#mainly for testing
function rand(L::OrdLoc{T}, scale = (-100:100)) where {T <: nf_elem}#rand
   Qx,x = FlintQQ["x"]
   K = nf(L)
   d = degree(K)
   while true
      temp = K(rand(Qx, 0:d-1, scale))
      try
         temp = L(temp)
         return temp
      catch
      end
   end
end

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

promote_rule(::Type{OrdLocElem{T}}, ::Type{OrdLocElem{T}}) where {T <: nf_elem} = OrdLocElem{T}

###############################################################################
#
#   Parent object call overloading
#
###############################################################################

(L::OrdLoc{T})() where {T <: nf_elem} = L(zero(nf(L)))

(L::OrdLoc{T})(a::Integer)  where {T <: nf_elem} = L(nf(L)(a))

function (L::OrdLoc{T})(data::T, checked::Bool = true) where {T <: nf_elem}
   return OrdLocElem{T}(data,L,checked)
end

function (L::OrdLoc{T})(data::NfAbsOrdElem{AnticNumberField,T}, checked::Bool = true) where {T <: nf_elem}
   return OrdLocElem{T}(nf(parent(data))(data),L,checked)
end

function (L::OrdLoc{T})(data::Rational{<: Integer}, checked::Bool = true) where {T <: nf_elem}
   return OrdLocElem{T}(nf(L)(numerator(data)) // nf(L)(denominator(data)),L,checked)
end

function (L::OrdLoc{T})(data::fmpz, checked::Bool = true) where {T <: nf_elem}
   return OrdLocElem{T}(nf(L)(data),L,checked)
end

function (L::OrdLoc{T})(a::OrdLocElem{T}) where {T <: nf_elem}
   L != parent(a) && error("No element of $L")
   return a
end

################################################################################
#
#   Valuation
#
################################################################################

@doc Markdown.doc"""
    valuation(a::OrdLocElem{T}, prime::NfAbsOrdIdl{AnticNumberField,T}) where {T <: nf_elem}
> Returns the valuation `n` of $a$ at $P$.
"""
valuation(a::OrdLocElem{T}, prime::NfAbsOrdIdl{AnticNumberField,T}) where {T <: nf_elem} = valuation(data(a), prime)

###############################################################################
#
#   Canonicalisation
#
###############################################################################

@doc Markdown.doc"""
    canonical_unit(a::OrdLocElem{T}) where {T <: nf_elem}
> Returns unit `b`::OrdLocElem{T} s.th. ($a$ * inv(`b`)) is hopefully nicer.
"""
function canonical_unit(a::OrdLocElem{T}) where {T <: nf_elem}
   iszero(a) && return parent(a)(1)
   isunit(a) && return a
   L = parent(a)
   if L.comp
     return L(1)
   else
     na, _ = _make_legal(a.data, L.prime)
     A = ppio(L.OK(na)*L.OK, L.prime)[1] * L.prime
     z = mod(L.OK(na), A)
     #L(z) should be equivalent to a, since v_p(A) = v_p(a) + 1`for all p | prime
     # u = a//z
     u = divexact(a, L(z))
     @assert isunit(u)
     return divexact(a, L(z))
   end

end

###############################################################################
#
#   Constructors
#
###############################################################################

@doc Markdown.doc"""
    Localization(OK::NfAbsOrd{AnticNumberField,T}, S::NfAbsOrdIdl{AnticNumberField,T}; cached=true, comp = false) where {T <: nf_elem}
> Returns the localization of the order $OK$ at the ideal $S$.
> If `cached == true` (the default) then the resulting
> localization parent object is cached and returned for any subsequent calls
> to the constructor with the same order $OK$ and ideal $S$.
> `comp == false` means primes dividing $S$ are invertable,
> `comp == true` means all primes not dividing $S$ become units.
"""
function Localization(OK::NfAbsOrd{AnticNumberField,T}, S::NfAbsOrdIdl{AnticNumberField,T}; cached=true, comp::Bool = false) where {T <: nf_elem}
   return OrdLoc{T}(OK, S, cached, comp)
end


text:
module TestInternalBags

using ArgCheck: @check
using ProgressLogging: Progress
using Reagents.Internal: Internal, Bag
using Test
using UUIDs: uuid4
using ..Utils: random_sleep, run_concurrently

struct Positive{T}
    value::T
end

Internal.isdeleted(p::Positive) = p.value <= 0

function test_concurrent_push_iterate(; nrepeat = 10, kwargs...)
    id = uuid4()
    try
        for i in 1:nrepeat
            @debug Progress(
                id,
                (i - 1) / nrepeat;
                name = "`test_concurrent_push_iterate`: trial $i",
            )
            check_concurrent_push_iterate(; kwargs...)
        end
    finally
        @debug Progress(id; done = true)
    end
    @test true
end

function check_concurrent_push_iterate(; randomize = true, nitems = 100)
    items = Bag{Positive{Int}}()
    sentinel = typemax(Int)

    receivers = map(1:5) do i
        function receiver_task()
            while true
                for p in items
                    if p.value == sentinel
                        @debug "`check_concurrent_push_iterate`: Iterator $i done"
                        return
                    end
                    randomize && random_sleep()
                end
                yield()
            end
        end
    end

    nsenders = 5
    nfinished = Threads.Atomic{Int}(0)
    senders = map(1:nsenders) do i
        function sender_task()
            for j in 1:nitems
                push!(items, Positive(nsenders * (j - 1) + i))
                push!(items, Positive(-1))
                randomize && random_sleep()
            end
            if Threads.atomic_add!(nfinished, 1) + 1 == nsenders
                # Note: using `nfinished` instead of using a task that waits for
                # all senders, so that we can use a single `run_concurrently`
                # call (to make debugging easier).

                push!(items, Positive(sentinel))
                @debug "`check_concurrent_push_iterate`: All sent"
            end
            @debug "`check_concurrent_push_iterate`: Sender $i done"
        end
    end
    allsent = 1:nsenders*nitems

    run_concurrently([senders; receivers]; spawn = true)
    @debug "`check_concurrent_push_iterate`: All sendres and receivers done"

    allreceived = sort!(Int[p.value for p in items])
    if allreceived[end] == sentinel
        pop!(allreceived)
    end
    @check length(allreceived) == length(allsent)
    @check sort!(allreceived) == allsent
end

end  # module


text:
#=
Implements the model described in Fogli and Veldkamp, "Germs, Social Networks, and Growth."
@author: Arnav Sood
@date: 2017-09-27
=#

# Name module.
module networksBackend

# Export types.
export Node, Model, ModelState, DataHolder, ModelState, NodeDist

# Export functions.
export τ, isSick, Node, ourMod, runModel, update!, initialize, nodes

# Load dependencies. Alias some objects with our notation.
using Distributions
using LightGraphs
const Probability = Real
const nodes = vertices
import Base.rand

# Define the objects used in the model.
type Node
    θ::Bool # Collectivist flag.
    nf::Integer # Degree. Specifically, twice the invariant number of outgoing (left) bonds.
    m::Probability # Mobility probability.
    A::Number # Technology.
    mobilityFlag::Bool # Whether the agent has a mobility link or doesn't.
end

type Model
    ξ::Distribution # Shock distribution.
    T::Integer # Number of time periods.
    N::Integer # Number of agents.
    F1::Bool # Innovation feedback.
    F2::Bool # Infection feedback.
    λ::Real # Innovation Poisson parameter.
    δ::Real # Innovation gain.
    A₀::Real # Initial technology level.
    π::Probability # Disease transmission probability.
    ϕ::Probability # Technology transmission probability.
    α::Distribution # Distribution of initial node types.
    Δ::Integer # Spacing for individualist nodes.
    Λ::Real # Feedback 1 parameter.
    Π::Real # Feedback 2 parameter.
end

type DataHolder
    typeChanges::Integer # Total number of type changes.
    maxTech::Real # Tech of node with the highest tech.
    frontierVec::Array # Vector of times when maxTech is pushed forward.
end

# Create α.
type NodeDist <: Distribution{Multivariate, Discrete}
    θ::Distribution # Bernoulli type.
    nf # Degree.
    m # Mobility.
end

type ModelState
    m::Model
    g::DiGraph # Model state represented as a directed graph; left arrows are outgoing, and right arrows are incoming.
    nodeList::Array{Node} # Vector of nodes corresponding to each index.
    runtimeData::DataHolder # Captures data we need to track as the state evolves.
end

# Create accessors and constructors for these objects.
function rand(α::NodeDist, n::Int64=1) # For τ = rand(m.α) in fillNode.
    rtn = [rand(α.θ), rand(α.nf), rand(α.m)]
    return rtn
end

function τ(n::Node)
    return [n.θ, n.nf, n.m]
end

function isSick(n::Node)
    return iszero(n.A)
end

function Node(θ, nf, m)
    return Node(θ, nf, m, 0) # If unset, nodes initialize with A = 0.
end

function fillNode(m::Model)
    τ = rand(m.α) # Draw a node type.
    return Node(τ..., m.A₀, false) # ... unpacks f([x,y,z],γ) into f(x,y,z,γ).
end

function ModelState(m::Model, g::DiGraph, nodeList)
    return ModelState(m, g, nodeList, DataHolder(0, m.A₀, []))
end

# Modular arithmetic on [1...N] (instead of the default [0...N-1]).
function ourMod(n::Number, N::Integer)
    if n <= N && 0 < n
        return n
    elseif n >= N
        return ourMod(n - N, N)
    elseif n <= 0
        return ourMod(n + N, N)
    end
end

# Draw bonds based on type.
function setBonds(g::DiGraph, m::Model, indexSet::Array,
nodeList::Array{Node})
    N = m.N
    for index in indexSet
        nf = nodeList[index].nf
        if nodeList[index].θ == true
            # Collectivist.
            friends = collect(1:nf/2)
            for f in friends
                add_edge!(g, index, ourMod(index + f, N))
            end
        elseif nodeList[index].θ == false
            # Individualist.
            friends = [1]
            append!(friends, collect(m.Δ : m.Δ -2 + nf/2))
            for f in friends
                add_edge!(g, index, ourMod(index + f, N))
            end
        end
    end
end

# Initialize model.
function initialize(m::Model)
    # Create objects.
    nodeList = Array{Node}(m.N)
    g = DiGraph(m.N)
    # Fill the objects.
    map!(n -> fillNode(m), nodeList, nodeList)
    setBonds(g, m, collect(1:m.N), nodeList)
    # Create a ModelState around them.
    return ModelState(m, g, nodeList)
end

# Do technology transmission.
function updateTech!(s::ModelState, ϕ::Distribution)
    # Create a changes list, so that we run each node's changes independently.
    # deepcopy() ensures that we're actually creating an independent object.
    newList = deepcopy(s.nodeList)
    for nodeIndex in nodes(s.g)
        friendsList = all_neighbors(s.g, nodeIndex)
        techList = [s.nodeList[f].A for f in friendsList]
        # Do the shock for each.
        map!(n -> rand(ϕ) == 1 ? n : 0, techList, techList)
        bestNeighbor = maximum(techList)
        # Make sure the new tech is actually better than the current one.
        if bestNeighbor >= s.nodeList[nodeIndex].A
            newList[nodeIndex].A = bestNeighbor
        end
    end
    # Replace.
    s.nodeList = newList
end

# Do disease transmission.
function updateDisease!(s::ModelState, π::Distribution)
    newList = deepcopy(s.nodeList)
    for nodeIndex in nodes(s.g)
        numSick = count(isSick, [s.nodeList[n] for n in all_neighbors(s.g, nodeIndex)])
        draw = rand(π, numSick)
        # If at least one shock is positive, we infect the current node.
        if 1 ∈ draw
            newList[nodeIndex].A = 0
        end
    end
    s.nodeList = newList
end

# Tech innovation.
function updateInnovation!(s::ModelState, λ::Real)
    shock = Poisson(λ)
    # No need to store changes in a copy, since these events are individual.
    for node in s.nodeList
        draw = rand(shock)
        newTech = (1 + draw * s.m.δ) * node.A
        node.A = newTech
    end
end

# Check technology adoption lag.
function lagCheck!(s::ModelState, t::Integer)
    oldMax = s.runtimeData.maxTech;
    newMax = maximum(n.A for n in s.nodeList);
    if newMax > oldMax
        # print_with_color(32, "New tech max attained at time $(t)\n")
        s.runtimeData.maxTech = newMax
        push!(s.runtimeData.frontierVec, t)
    end
end

# Perform tech frontier-related reporting.
function checkFrontiers!(s::ModelState, t::Integer)
    lagCheck!(s::ModelState, t::Integer)
end

# Update types.
function updateTypes!(s::ModelState, ξ::Array)
    newList = deepcopy(s.nodeList)
    changedIndices = find(ξ)
    for nodeIndex in changedIndices
        friendsList = all_neighbors(s.g, nodeIndex)
        techList = [s.nodeList[f].A for f in friendsList]
        A, loc = findmax(techList)
        # Increment type-change flag as necessary.
        if τ(s.nodeList[loc]) != τ(s.nodeList[nodeIndex])
            s.runtimeData.typeChanges += 1
        end
        # Replace type.
        newList[nodeIndex] = Node(τ(s.nodeList[loc])..., A, false) # false indicates that nodes start without mobility links.
    end
    s.nodeList = newList
end

# Type-changed agents update bonds.
function updateBonds!(s::ModelState, ξ::Array)
    changedIndices = find(ξ)
    for nodeIndex in changedIndices
        # Break all outgoing bonds.
        neighbors = deepcopy(out_neighbors(s.g, nodeIndex))
        for neighborIndex in neighbors
            rem_edge!(s.g, nodeIndex, neighborIndex)
        end
    end
    # Redraw new outgoing bonds, based on the node's type.
    setBonds(s.g, s.m, changedIndices, s.nodeList)
end

# Redraw bonds for mobile agents.
function redrawBonds!(s::ModelState, mobList::Array)
    newGraph = deepcopy(s.g)
    for nodeIndex in mobList
        shock = Bernoulli(s.nodeList[nodeIndex].m)
        draw = rand(shock)
        if s.nodeList[nodeIndex].mobilityFlag == false && draw == 1 # No existing mobility links, and a valid shock.
            neighbors = deepcopy(out_neighbors(s.g, nodeIndex)) # Neighbors from old graph.
            neighborIndex = rand(neighbors) # Pick one.
            rem_edge!(newGraph, nodeIndex, neighborIndex) # Break the corresponding bond on the new graph.
            unconnected = setdiff(collect(1:s.m.N), all_neighbors(newGraph, nodeIndex))
            # Remove the actual node.
            unconnected = setdiff(unconnected, [nodeIndex])
            newNeighbor = rand(unconnected)
            add_edge!(newGraph, nodeIndex, newNeighbor) # Add the new edge on the new graph.
            # println("Replaced edge $(neighborIndex) with edge $(newNeighbor) for $(nodeIndex)")
            s.nodeList[nodeIndex].mobilityFlag = true; # Update mobility flag.
        end
    end
    # Replace the graph.
    s.g = newGraph
end

# Mobile agents redraw bonds.
function updateMobile!(s::ModelState)
    mobilities = [n.m for n in s.nodeList]
    mobileAgents = find(mobilities)
    # println("Mobile agents:")
    # println(mobileAgents)
    redrawBonds!(s, mobileAgents)
end

function update!(s::ModelState, t::Integer)
    # Unpack parameters, draw ξ shock vector, etc.
    m = s.m
    g = s.g
    nodeList = s.nodeList
    π = Bernoulli(m.π)
    Π = m.Π
    ϕ = Bernoulli(m.ϕ)
    ξ = rand(m.ξ, m.N) # This is an Array, instead of a Distribution!
    λ = m.λ
    Λ = m.Λ
    N = m.N
    # Update parameters for any feedbacks.
    if m.F1 == true
        numSick = count(isSick, nodeList)
        avgSick = numSick/N
        λ = 1 - cdf(Normal(0, 1), avgSick/Λ)
    end
    if m.F2 == true
        techList = [n.A for n in nodeList]
        meanTech = mean(techList)
        ϕ = Bernoulli(1-cdf(Normal(0, 1), meanTech/Π))
    end
    updateTech!(s,ϕ)
    updateDisease!(s,π)
    updateInnovation!(s, λ)
    checkFrontiers!(s, t::Integer)
    updateTypes!(s, ξ)
    updateBonds!(s, ξ)
    updateMobile!(s)
end

# Run model.
function runModel(m::Model)
    state = initialize(m)
    for t in 1:m.T # So initialization is at time 0.
        # print_with_color(32, "t = $(t) \n")
        update!(state, t)
    end
    print_with_color(32, "Run complete \n")
    return state
end

# End module.
end

max_stars_repo_name: ArnavSood/networks-paper

max_issues_repo_name: arnavs/networks-paper

max_forks_repo_name: arnavs/networks-paper


text:
using Test
import ArchGDAL; const AG = ArchGDAL

@testset "Testing ConfigOptions" begin
    @test AG.getconfigoption("GDAL_CACHEMAX") == ""
    AG.setconfigoption("GDAL_CACHEMAX", "64")
    @test AG.getconfigoption("GDAL_CACHEMAX") == "64"
    AG.clearconfigoption("GDAL_CACHEMAX")
    @test AG.getconfigoption("GDAL_CACHEMAX", "128") == "128"

    @test AG.getthreadconfigoption("GDAL_CACHEMAX") == ""
    AG.setthreadconfigoption("GDAL_CACHEMAX","32")
    @test AG.getthreadconfigoption("GDAL_CACHEMAX") == "32"
    AG.clearthreadconfigoption("GDAL_CACHEMAX")
    @test AG.getthreadconfigoption("GDAL_CACHEMAX", "128") == "128"

    @test AG.getconfigoption("GDAL_CACHEMAX") == ""
    @test AG.getconfigoption("CPL_LOG_ERRORS") == ""
    @test AG.getthreadconfigoption("GDAL_CACHEMAX") == ""
    @test AG.getthreadconfigoption("CPL_LOG_ERRORS") == ""

    AG.environment(globalconfig=[("GDAL_CACHEMAX","64"),
                                 ("CPL_LOG_ERRORS","ON")],
                   threadconfig=[("GDAL_CACHEMAX","32"),
                                 ("CPL_LOG_ERRORS","OFF")]) do
        # it seems that thread settings overwrites global settings?
        @test AG.getconfigoption("GDAL_CACHEMAX") == "32"
        @test AG.getconfigoption("CPL_LOG_ERRORS") == "OFF"
        @test AG.getthreadconfigoption("GDAL_CACHEMAX") == "32"
        @test AG.getthreadconfigoption("CPL_LOG_ERRORS") == "OFF"
    end

    AG.environment(globalconfig=[("GDAL_CACHEMAX","64"),
                                 ("CPL_LOG_ERRORS","ON")]) do
        # everything normal here
        @test AG.getconfigoption("GDAL_CACHEMAX") == "64"
        @test AG.getconfigoption("CPL_LOG_ERRORS") == "ON"
        @test AG.getthreadconfigoption("GDAL_CACHEMAX") == ""
        @test AG.getthreadconfigoption("CPL_LOG_ERRORS") == ""
    end
end

@testset "Test Driver Capabilities" begin
    drivers = AG.listdrivers()
    @test drivers["GTiff"] == "GeoTIFF"
    @test length(AG.driveroptions("GTiff")) > 100
    @test sprint(print, AG.identifydriver("data/point.geojson")) == "Driver: GeoJSON/GeoJSON"
    @test sprint(print, AG.identifydriver("data/utmsmall.tif")) == "Driver: GTiff/GeoTIFF"

    @test AG.validate(AG.getdriver("GTiff"), ["COMPRESS=LZW", "INTERLEAVE=PIXEL"]) == true
    @test AG.validate(AG.getdriver("GTiff"), ["COMPRESS=LZW"]) == true
    @test AG.validate(AG.getdriver("GTiff"), ["INTERLEAVE=PIXEL"]) == true
    AG.read("data/point.geojson") do dataset
        @test AG.listcapability(dataset) == Dict(
            "CreateLayer"=>false,
            "DeleteLayer"=>false,
            "CreateGeomFieldAfterCreateLayer"=>false,
            "CurveGeometries"=>false,
            "Transactions"=>false,
            "EmulatedTransactions"=>false
        )
        @test AG.listcapability(AG.getlayer(dataset,0)) == Dict(
            "SequentialWrite"=>false,    "DeleteField"=>false,
            "IgnoreFields"=>false,       "FastSpatialFilter"=>false,
            "DeleteFeature"=>false,      "FastFeatureCount"=>true,
            "StringsAsUTF8"=>true,       "CreateGeomField"=>false,
            "ReorderFields"=>false,      "MeasuredGeometries"=>true,
            "FastSetNextByIndex"=>true, "CreateField"=>false,
            "RandomWrite"=>false,        "RandomRead"=>true,
            "CurveGeometries"=>false,    "FastGetExtent"=>false,
            "Transactions"=>false,       "AlterFieldDefn"=>false
        )
    end
end


text:
"""
    struct DirichletFESpace <: SingleFieldFESpace
      space::SingleFieldFESpace
    end
"""
struct DirichletFESpace{S<:SingleFieldFESpace} <: SingleFieldFESpace
  space::S
  function DirichletFESpace(space::SingleFieldFESpace)
    new{typeof(space)}(space)
  end
end

ConstraintStyle(::Type{DirichletFESpace{B}}) where B = ConstraintStyle(B)

CellField(t::DirichletFESpace,cell_vals) = CellField(t.space,cell_vals)

get_cell_isconstrained(f::DirichletFESpace) = get_cell_isconstrained(f.space)

get_cell_constraints(f::DirichletFESpace) = get_cell_constraints(f.space)

function get_free_dof_ids(f::DirichletFESpace)
  get_dirichlet_dof_ids(f.space)
end

function get_vector_type(f::DirichletFESpace)
  get_vector_type(f.space)
end

function get_cell_dof_ids(f::DirichletFESpace)
  lazy_map(Broadcasting(-),get_cell_dof_ids(f.space))
end

function get_dirichlet_dof_ids(f::DirichletFESpace)
  get_free_dof_ids(f.space)
end

function num_dirichlet_tags(f::DirichletFESpace)
  1
end

function get_dirichlet_dof_tag(f::DirichletFESpace)
  ones(Int8,num_dirichlet_dofs(f))
end

function scatter_free_and_dirichlet_values(f::DirichletFESpace,fv,dv)
  scatter_free_and_dirichlet_values(f.space,dv,fv)
end

function gather_free_and_dirichlet_values!(fv,dv,f::DirichletFESpace,cv)
  gather_free_and_dirichlet_values!(dv,fv,f.space,cv)
  (fv,dv)
end

function TrialFESpace(f::DirichletFESpace)
  U = TrialFESpace(f.space)
  DirichletFESpace(U)
end

function get_cell_shapefuns(f::DirichletFESpace)
  get_cell_shapefuns(f.space)
end

function get_cell_shapefuns_trial(f::DirichletFESpace)
  get_cell_shapefuns_trial(f.space)
end

function get_cell_dof_basis(f::DirichletFESpace)
  get_cell_dof_basis(f.space)
end

function get_triangulation(f::DirichletFESpace)
  get_triangulation(f.space)
end


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/FEMBeam.jl/blob/master/LICENSE

# Read ABAQUS .mtx files

s = """
** Assembled nodal loads
*CLOAD, REAL
2,3,  5.400000000000027e+04
2,5,  9.000000000000071e+03
"""

t = """
-1,1, -1,1,  3.333333333333317e+01
-1,2, -1,1, -8.333333333333492e+00
1,1, -1,1,  2.500000000000002e+01
2,1, -1,1, -2.500000000000002e+01
-1,2, -1,2,  3.333333333333316e+01
1,1, -1,2,  2.500000000000003e+01
2,1, -1,2, -2.500000000000003e+01
1,1, 1,1,  1.000000000000000e+36
2,1, 1,1, -3.000000000000020e+02
1,2, 1,2,  1.000000000000000e+36
1,6, 1,2,  5.999999999999996e+00
2,2, 1,2, -1.199999999999999e+01
2,6, 1,2,  5.999999999999996e+00
1,3, 1,3,  1.000000000000000e+36
1,5, 1,3, -8.999999999999996e+00
2,3, 1,3, -1.799999999999999e+01
2,5, 1,3, -8.999999999999996e+00
1,4, 1,4,  1.000000000000000e+36
2,4, 1,4, -2.400000000000012e-01
1,5, 1,5,  1.000000000000000e+36
2,3, 1,5,  8.999999999999996e+00
2,5, 1,5,  2.999999999999990e+00
1,6, 1,6,  1.000000000000000e+36
2,2, 1,6, -5.999999999999996e+00
2,6, 1,6,  1.999999999999993e+00
2,1, 2,1,  3.000000000000020e+02
2,2, 2,2,  1.199999999999999e+01
2,6, 2,2, -5.999999999999996e+00
2,3, 2,3,  1.799999999999999e+01
2,5, 2,3,  8.999999999999996e+00
2,4, 2,4,  2.400000000000012e-01
2,5, 2,5,  6.000000000000006e+00
2,6, 2,6,  4.000000000000004e+00
"""

function read_mtx(data; dim=0)
    if dim == 0
        for ln in eachline(copy(data))
            i,idof,j,jdof,value = map(Meta.parse, split(ln, ','))
            dim = max(dim, idof, jdof)
        end
    end
    I = Int64[]
    J = Int64[]
    V = Float64[]
    for ln in eachline(data)
        i,idof,j,jdof,value = map(Meta.parse, split(ln, ','))
        if i < 1 || j < 1
            continue
        end
        push!(I, (i-1)*dim+idof)
        push!(J, (j-1)*dim+jdof)
        push!(V, value)
    end
    A = Matrix(sparse(I, J, V))
    A += transpose(tril(A,-1))
    return A
end

function read_mtx_from_str(data)
    read_mtx(IOBuffer(data))
end

#A = read_mtx(IOBuffer(t))
#display(A)


text:
##############################################################################
#
# A macro for doing delegation
#
# This macro call
#
#     @delegate MyContainer.elems [:size, :length, :ndims, :endof]
#
# produces this block of expressions
#
#     size(a::MyContainer) = size(a.elems)
#     length(a::MyContainer) = length(a.elems)
#     ndims(a::MyContainer) = ndims(a.elems)
#     endof(a::MyContainer) = endof(a.elems)
#
##############################################################################

macro delegate(source, targets)
    typename = esc(source.args[1])
    sourceargs2 = source.args[2]
    fieldname = Expr(:quote, isa(sourceargs2, QuoteNode) ? sourceargs2.value : sourceargs2.args[1])
    funcnames = targets.args
    n = length(funcnames)
    result = quote begin end end
    for i in 1:n
        funcname = esc(funcnames[i])
        f = quote
            ($funcname)(a::($typename), args...; kwargs...) = ($funcname)(getfield(a, $fieldname), args...; kwargs...)
        end
        push!(result.args[2].args, f)
    end
    return result
end

"""
Wrapper for a `StatisticalModel` that has been fit from a `@formula` and tabular
data.

Most functions from the StatsBase API are simply delegated to the wrapped model,
with the exception of functions like `fit`, `predict`, and `coefnames` where the
tabular nature of the data means that additional processing is required or
information provided by the formula.

# Fields
* `model::M` the wrapped `StatisticalModel`.
* `mf::ModelFrame` encapsulates the formula, schema, and model type.
* `mm::ModelMatrix{T}` the model matrix that the model was fit from.
"""
struct TableStatisticalModel{M,T} <: StatisticalModel
    model::M
    mf::ModelFrame
    mm::ModelMatrix{T}
end

"""
Wrapper for a `RegressionModel` that has been fit from a `@formula` and tabular
data.

Most functions from the StatsBase API are simply delegated to the wrapped model,
with the exception of functions like `fit`, `predict`, and `coefnames` where the
tabular nature of the data means that additional processing is required or
information provided by the formula.

# Fields
* `model::M` the wrapped `RegressioModel`.
* `mf::ModelFrame` encapsulates the formula, schema, and model type.
* `mm::ModelMatrix{T}` the model matrix that the model was fit from.
"""
struct TableRegressionModel{M,T} <: RegressionModel
    model::M
    mf::ModelFrame
    mm::ModelMatrix{T}
end

for (modeltype, dfmodeltype) in ((:StatisticalModel, TableStatisticalModel),
                                 (:RegressionModel, TableRegressionModel))
    @eval begin
        function StatsBase.fit(::Type{T}, f::FormulaTerm, data, args...;
                               contrasts::Dict{Symbol,<:Any} = Dict{Symbol,Any}(),
                               kwargs...) where T<:$modeltype

            Tables.istable(data) || throw(ArgumentError("expected data in a Table, got $(typeof(data))"))
            cols = columntable(data)

            mf = ModelFrame(f, cols, model=T, contrasts=contrasts)
            mm = ModelMatrix(mf)
            y = response(mf)
            $dfmodeltype(fit(T, mm.m, y, args...; kwargs...), mf, mm)

            ## TODO: consider doing this manually, without the ModelFrame/ModelMatrix
            # schema = schema(data, cols, contrasts)
            # f = apply_schema(f, schema, T)
            # y, X = modelcols(f, cols)
            # $dfmodeltype(fit(T, X, y, args...; kwargs...))
        end
    end
end

"""
    formula(model)

Retrieve formula from a fitted or specified model
"""
function formula end

formula(m::TableStatisticalModel) = m.mf.f
formula(m::TableRegressionModel) = m.mf.f

@doc """
    fit(Mod::Type{<:StatisticalModel}, f::FormulaTerm, data, args...;
        contrasts::Dict{Symbol}, kwargs...)

Convert tabular data into a numeric response vector and predictor matrix using
the formula `f`, and then `fit` the specified model type, wrapping the result in
a [`TableRegressionModel`](@ref) or [`TableStatisticalModel`](@ref) (as
appropriate).

This is intended as a backstop for modeling packages that implement model types
that are subtypes of `StatsBase.StatisticalModel` but do not explicitly support
the full StatsModels terms-based interface.  Currently this works by creating a
[`ModelFrame`](@ref) from the formula and data, and then converting this to a
[`ModelMatrix`](@ref), but this is an internal implementation detail which may
change in the near future.
""" fit

# Delegate functions from StatsBase that use our new types
const TableModels = Union{TableStatisticalModel, TableRegressionModel}
@delegate TableModels.model [StatsBase.coef, StatsBase.confint,
                             StatsBase.deviance, StatsBase.nulldeviance,
                             StatsBase.loglikelihood, StatsBase.nullloglikelihood,
                             StatsBase.dof, StatsBase.dof_residual, StatsBase.nobs,
                             StatsBase.stderror, StatsBase.vcov, StatsBase.fitted]
@delegate TableRegressionModel.model [StatsBase.modelmatrix,
                                      StatsBase.residuals, StatsBase.response,
                                      StatsBase.predict, StatsBase.predict!,
                                      StatsBase.cooksdistance]
StatsBase.predict(m::TableRegressionModel, new_x::AbstractMatrix; kwargs...) =
    predict(m.model, new_x; kwargs...)
# Need to define these manually because of ambiguity using @delegate

StatsBase.r2(mm::TableRegressionModel) = r2(mm.model)
StatsBase.adjr2(mm::TableRegressionModel) = adjr2(mm.model)
StatsBase.r2(mm::TableRegressionModel, variant::Symbol) = r2(mm.model, variant)
StatsBase.adjr2(mm::TableRegressionModel, variant::Symbol) = adjr2(mm.model, variant)
StatsBase.loglikelihood(mm::TableModels, c::Colon) = loglikelihood(mm.model, c)

isnested(m1::TableModels, m2::TableModels; kwargs...) = isnested(m1.model, m2.model; kwargs...)

function _return_predictions(T, yp::AbstractVector, nonmissings, len)
    out = Vector{Union{eltype(yp),Missing}}(missing, len)
    out[nonmissings] = yp
    out
end

function _return_predictions(T, yp::AbstractMatrix, nonmissings, len)
    out = Matrix{Union{eltype(yp),Missing}}(missing, len, 3)
    out[nonmissings, :] = yp
    T((prediction = out[:,1], lower = out[:,2], upper = out[:,3]))
end

function _return_predictions(T, yp::NamedTuple, nonmissings, len)
    y = Vector{Union{eltype(yp.prediction),Missing}}(missing, len)
    l, h = similar(y), similar(y)
    out = (prediction = y, lower = l, upper = h)
    for key in (:prediction, :lower, :upper)
        out[key][nonmissings] = yp[key]
    end
    T(out)
end

# Predict function that takes data table as predictor instead of matrix
function StatsBase.predict(mm::TableRegressionModel, data; kwargs...)
    Tables.istable(data) ||
        throw(ArgumentError("expected data in a Table, got $(typeof(data))"))

    f = mm.mf.f
    cols, nonmissings = missing_omit(columntable(data), f.rhs)
    new_x = modelcols(f.rhs, cols)
    y_pred = predict(mm.model, reshape(new_x, size(new_x, 1), :);
                     kwargs...)
    _return_predictions(Tables.materializer(data), y_pred, nonmissings, length(nonmissings))
end

StatsBase.coefnames(model::TableModels) = coefnames(model.mf)

# coeftable implementation
function StatsBase.coeftable(model::TableModels; kwargs...)
    ct = coeftable(model.model, kwargs...)
    cfnames = coefnames(model.mf)
    if length(ct.rownms) == length(cfnames)
        ct.rownms = cfnames
    end
    ct
end

# show function that delegates to coeftable
function Base.show(io::IO, model::TableModels)
    try
        ct = coeftable(model)
        println(io, typeof(model))
        println(io)
        println(io, model.mf.f)
        println(io)
        println(io,"Coefficients:")
        show(io, ct)
    catch e
        if isa(e, ErrorException) && occursin("coeftable is not defined", e.msg)
            show(io, model.model)
        else
            rethrow(e)
        end
    end
end


text:
module DiffEqFlux

using DiffEqBase, Zygote, DiffResults, DiffEqSensitivity, ForwardDiff,
      Flux, Requires, Adapt, LinearAlgebra

gpu_or_cpu(x) = Array
function __init__()
    @require CuArrays="3a865a2d-5b23-5a0f-bc46-62713ec82fae" begin
        gpu_or_cpu(x::CuArrays.CuArray) = CuArrays.CuArray
        gpu_or_cpu(x::Transpose{<:Any,<:CuArrays.CuArray}) = CuArrays.CuArray
        gpu_or_cpu(x::Adjoint{<:Any,<:CuArrays.CuArray}) = CuArrays.CuArray
    end
end

include("Flux/layers.jl")
include("Flux/neural_de.jl")
include("Flux/utils.jl")

export diffeq_adjoint
#= export diffeq_fd, diffeq_rd, diffeq_adjoint =#
#= export neural_ode, neural_ode_rd =#
export neural_ode
#= export neural_dmsde =#
end


text:
"""
    Multiply(input=inputDimension, output=outputDimension, winit=xavier, atype=KnetLayers.arrtype)

Creates a matrix multiplication layer based on `inputDimension` and `outputDimension`.
    (m::Multiply)(x::AbstractArray{T,N}; keepsize=true)= m.w * x

By default parameters initialized with xavier, you may change it with `winit` argument.

# Keywords
* `input=inputDimension`: input size
* `output=outputDimension`: output size
* `winit=xavier`: weight initialization distribution
* `atype=KnetLayers.arrtype` : underlying array type for parameters.
   Default arrtype=KnetArray{Float32} if running on a GPU device, otherwise Array{Float32}.
   see `KnetLayers.settype!`
* `keepsize=true`: When `N=3` or higher dimensional arrays supplied, `Multiply` squeezes the
                   dimensions except the first dimension in order to perform multiplication.
                   It reshapes the output to recover input sizeality `N`.
                   However, one may obtain 2D vesion by setting keepsize=false

"""
mutable struct Multiply{P} <: Layer
    weight::P
end
Multiply(;input::Int, output::Int, winit=xavier, atype=arrtype) = Multiply(param(output, input; init=winit, atype=atype))
@inline (m::Multiply)(x::AbstractArray{<:Integer}) = m.weight[:,x] # Lookup (EmbedLayer)

# TODO: Find a faster (or compound) way for tensor-product
function (m::Multiply)(x; keepsize=true)
    if ndims(x) > 2
        s = size(x)
        y = m.weight * reshape(x, s[1], prod(s[2:end]))
        return (keepsize ? reshape(y, size(y, 1), s[2:end]...) : y)
    else
        return m.weight * x
    end
end
Base.show(io::IO,m::Multiply{P}) where P = print(io,Multiply{P},"(input=",size(m.weight,2)," output=",size(m.weight,1),")")
@inline Base.size(m::Multiply) = size(m.weight)

"""
    Embed(input=inputSize, output=embedSize, winit=xavier, atype=KnetLayers.arrtype)
Creates an embedding layer according to given `inputSize` and `embedSize`.
By default parameters initialized with xavier, change it with `winit` argument.

    (m::Embed)(x::Array{T}) where T<:Integer
    (m::Embed)(x::AbstractArray{T,N}; keepsize=true) # see `Multiply`

Embed objects are callable with an input which is either and integer array
(one hot encoding) or an AbstractArray{T,N}
For N-dimensional array, `size(x,1)==inputSize`, see `Multiply`

# Keywords

* `input=inputDimension`: input size
* `output=embeddingDimension`: output size
* `winit=xavier`: weight initialization distribution
* `atype=KnetLayers.arrtype` : underlying array type for parameters.
   Default arrtype=KnetArray{Float32} if running on a GPU device, otherwise Array{Float32}.
   see `KnetLayers.settype!`
"""
const Embed = Multiply

"""
    Bias(sizes..., atype=KnetLayers.arrtype, winit=zeros)
Creates a bias layer according to given `sizes`
By default parameters initialized with zeros, change it with `winit` argument.
`Bias` layer can be initialized with Bias(nothing) to make it ineffective.

    (m::Bias)(x) = m.b .+ x
    (m::Bias{Nothing})(x) = x
# Keywords

* `winit=zeros`: weight initialization distribution
* `atype=KnetLayers.arrtype` : underlying array type for parameters.
   Default arrtype=KnetArray{Float32} if running on a GPU device, otherwise Array{Float32}.
   see `KnetLayers.settype!`
"""
struct Bias{T}
    b::T
end
@inline (m::Bias)(x) = m.b .+ x
@inline (m::Bias{Nothing})(x) = x
#Bias(size::Int;atype=arrtype, winit=zeros, o...) = Bias(param(size;atype=atype, init=winit,o...))
Bias(sizes::Int...;atype=arrtype, winit=zeros, o...) = Bias(param(sizes...;atype=atype, init=winit, o...))
Bias() = Bias(nothing)
Base.show(io::IO,m::Bias{T}) where T = print(io,Bias{T},"(length=",length(m.b),")")
@inline Base.size(m::Bias) = size(b)
"""
    Linear(input=inputSize, output=outputSize, winit=xavier, binit=zeros, atype=KnetLayers.arrtype)

Creates and linear layer according to given `inputSize` and `outputSize`.
`Linear` consists of a `Multiply` and a `Bias` layer.

    (m::Linear)(x) = m.bias(m.mult(x))

# Keywords
* `input=inputSize`   input size
* `output=outputSize` output size
* `winit=xavier`: weight initialization distribution
* `bias=zeros`: bias initialization distribution. It can be `nothing` for unbiased `Linear` layer.
* `atype=KnetLayers.arrtype` : underlying array type for parameters.
   Default arrtype=KnetArray{Float32} if running on a GPU device, otherwise Array{Float32}.
   see `KnetLayers.settype!`
"""
mutable struct Linear <: Layer
    mult::Multiply
    bias::Bias
end
function Linear(;input::Int, output::Int, winit=xavier, binit=zeros, atype=arrtype)
    Linear(Multiply(input=input, output=output, winit=winit, atype=atype),Bias(output, winit=binit, atype=atype))
end
@inline (m::Linear)(x) = m.bias(m.mult(x))
@inline Base.size(m::Linear) = size(m.mult)
"""
    Dense(input=inputSize, output=outputSize, activation=ReLU(), winit=xavier, binit=zeros, atype=KnetLayers.arrtype)
Creates a dense layer according to given `input` and `output` sizes.
If activation is `nothing`, it acts like a `Linear` Layer. see `Linear`

    (m::Dense{Nothing})(x) = m.linear(x)
    (m::Dense{<:Activation})(x)= m.activation(m.linear(x))

# Keywords
* `input=inputSize`   input size
* `output=outputSize` output size
* `winit=xaiver`: weight initialization distribution
* `bias=zeros`:   bias initialization distribution
* `activation=ReLU()`  an  activation layer, or nothing for a `Linear` equivalent layer.
* `atype=KnetLayers.arrtype` : underlying array type for parameters.
   Default arrtype=KnetArray{Float32} if running on a GPU device, otherwise Array{Float32}.
   see `KnetLayers.settype!`
"""
mutable struct Dense{T<:ActOrNothing} <: Layer
    linear::Linear
    activation::T
end

function Dense(;input::Int, output::Int, activation::ActOrNothing=ReLU(), winit=xavier, binit=zeros, atype=arrtype)
    Dense(Linear(input=input, output=output, winit=winit, binit=binit, atype=atype), activation)
end
@inline (m::Dense{Nothing})(x) = m.linear(x)
@inline (m::Dense{<:Activation})(x)= m.activation(m.linear(x))

Base.show(io::IO, x::Dense) = print(io,typeof(x),"(",x.linear,")")
@inline Base.size(m::Dense) = size(m.linear)
#TO-DO: Remove after the issue is resolved:
#https://github.com/denizyuret/Knet.jl/issues/418
"""
    BatchNorm(channels:Int;options...)
    (m::BatchNorm)(x;training=false) #forward run
# Options
* `momentum=0.1`: A real number between 0 and 1 to be used as the scale of
 last mean and variance. The existing running mean or variance is multiplied by
 (1-momentum).
* `mean=nothing': The running mean.
* `var=nothing`: The running variance.
* `meaninit=zeros`: The function used for initialize the running mean.
* `varinit=ones`: The function used for initialize the running variance of batchnorm
* `dataType=eltype(KnetLayers.arrtype)` : element type ∈ {Float32,Float64} for parameters.
* `usegpu=KnetLayers.arrtype <: KnetArray`
# Keywords
* `training`=nothing: When training is true, the mean and variance of x are used and moments
 argument is modified if it is provided. When training is false, mean and variance
 stored in the moments argument are used. Default value is true when at least one
 of x and params is AutoGrad.Value, false otherwise.
"""
mutable struct BatchNorm <: Layer
    params
    moments::BNMoments
end

function BatchNorm(channels::Int;  usegpu = arrtype <: KnetArray, dataType=eltype(arrtype), o...)
    w = bnparams(dataType,channels)
    m = bnmoments(;o...)
    p = usegpu ? Param(KnetArray(w)) : Param(w)
    BatchNorm(p,m)
end
@inline (m::BatchNorm)(x;o...) = batchnorm(x,m.moments,m.params;o...)
Base.show(io::IO,x::BatchNorm) where P = print(io,BatchNorm,"(",x.params,", ",x.moments,")")

"""
    Diagonal(D::Integer)
Creates an element-wise linear transformation layer with learnable
vectors `w` and `b`:
    y = w .* x .+ b
The input `x` must be a array where `size(x, 1) == D`.
"""
struct Diagonal{T} <: Layer
  w::T
  b::T
end

Diagonal(D::Integer; winit = ones, binit = zeros, atype=arrtype) =
  Diagonal(param(D; init=winit, atype=atype),
           param(D; init=binit, atype=atype))

function (a::Diagonal)(x)
  w, b = a.w, a.b
  return a.w.*x .+ a.b
end

Base.length(l::Diagonal) = length(l.w)

function Base.show(io::IO, l::Diagonal)
  print(io, "Diagonal(", length(l), ")")
end

"""
    LayerNorm(h::Integer)
A [normalisation layer](https://arxiv.org/pdf/1607.06450.pdf) designed to be
used with recurrent hidden states of size `h`. Normalises the mean/stddev of
each input before applying a per-neuron gain/bias.
"""
struct LayerNorm{T} <: Layer
  diag::Diagonal{T}
end

LayerNorm(h::Integer) =
  LayerNorm(Diagonal(h))

(a::LayerNorm)(x) = a.diag(normalise(x))

function Base.show(io::IO, l::LayerNorm)
  print(io, "LayerNorm(", length(l.diag), ")")
end


text:
using {PKGNAME}
using Test