text:
module SnoopCompile

using Serialization

export @snoopc
if VERSION >= v"1.2.0-DEV.573"
    export @snoopi

    const __inf_timing__ = Tuple{Float64,Core.MethodInstance}[]

    function typeinf_ext_timed(linfo::Core.MethodInstance, params::Core.Compiler.Params)
        tstart = time()
        ret = Core.Compiler.typeinf_ext(linfo, params)
        tstop = time()
        push!(__inf_timing__, (tstop-tstart, linfo))
        return ret
    end
    function typeinf_ext_timed(linfo::Core.MethodInstance, world::UInt)
        tstart = time()
        ret = Core.Compiler.typeinf_ext(linfo, world)
        tstop = time()
        push!(__inf_timing__, (tstop-tstart, linfo))
        return ret
    end

    function sort_timed_inf(tmin)
        data = __inf_timing__
        if tmin > 0
            data = filter(tl->tl[1] >= tmin, data)
        end
        return sort(data; by=tl->tl[1])
    end

    """
        inf_timing = @snoopi commands
        inf_timing = @snoopi tmin=0.0 commands

    Execute `commands` while snooping on inference. Returns an array of `(t, linfo)`
    tuples, where `t` is the amount of time spent infering `linfo` (a `MethodInstance`).

    Methods that take less time than `tmin` will not be reported.
    """
    macro snoopi(args...)
        tmin = 0.0
        if length(args) == 1
            cmd = args[1]
        elseif length(args) == 2
            a = args[1]
            if isa(a, Expr) && a.head == :(=) && a.args[1] == :tmin
                tmin = a.args[2]
                cmd = args[2]
            else
                error("unrecognized input ", a)
            end
        else
            error("at most two arguments are supported")
        end
        quote
            empty!($__inf_timing__)
            ccall(:jl_set_typeinf_func, Cvoid, (Any,), $typeinf_ext_timed)
            try
                $(esc(cmd))
            finally
                ccall(:jl_set_typeinf_func, Cvoid, (Any,), Core.Compiler.typeinf_ext)
            end
            $sort_timed_inf($tmin)
        end
    end
end

"""
```
@snoopc "compiledata.csv" begin
    # Commands to execute, in a new process
end
```
causes the julia compiler to log all functions compiled in the course
of executing the commands to the file "compiledata.csv". This file
can be used for the input to `SnoopCompile.read`.
"""
macro snoopc(flags, filename, commands)
    return :(snoopc($(esc(flags)), $(esc(filename)), $(QuoteNode(commands))))
end
macro snoopc(filename, commands)
    return :(snoopc(String[], $(esc(filename)), $(QuoteNode(commands))))
end

function snoopc(flags, filename, commands)
    println("Launching new julia process to run commands...")
    # addprocs will run the unmodified version of julia, so we
    # launch it as a command.
    code_object = """
            using Serialization
            while !eof(stdin)
                Core.eval(Main, deserialize(stdin))
            end
            """
    process = open(`$(Base.julia_cmd()) $flags --eval $code_object`, stdout, write=true)
    serialize(process, quote
        let io = open($filename, "w")
            ccall(:jl_dump_compiles, Nothing, (Ptr{Nothing},), io.handle)
            try
                $commands
            finally
                ccall(:jl_dump_compiles, Nothing, (Ptr{Nothing},), C_NULL)
                close(io)
            end
        end
        exit()
    end)
    wait(process)
    println("done.")
    nothing
end

function split2(str, on)
    i = findfirst(isequal(on), str)
    i === nothing && return str, ""
    return (SubString(str, firstindex(str), prevind(str, first(i))),
            SubString(str, nextind(str, last(i))))
end

"""
`SnoopCompile.read("compiledata.csv")` reads the log file produced by the compiler and returns the functions as a pair of arrays. The first array is the amount of time required to compile each function, the second is the corresponding function + types. The functions are sorted in order of increasing compilation time. (The time does not include the cost of nested compiles.)
"""
function read(filename)
    times = Vector{UInt64}()
    data = Vector{String}()
    toplevel = false # workaround for Julia#22538
    for line in eachline(filename)
        if toplevel
            if endswith(line, '"')
                toplevel = false
            end
            continue
        end
        time, str = split2(line, '\t')
        length(str) < 2 && continue
        (str[1] == '"' && str[end] == '"') || continue
        if startswith(str, """"<toplevel thunk> -> """)
            # consume lines until we find the terminating " character
            toplevel = true
            continue
        end
        tm = tryparse(UInt64, time)
        tm === nothing && continue
        push!(times, tm)
        push!(data, str[2:prevind(str, lastindex(str))])
    end
    # Save the most costly for last
    p = sortperm(times)
    return (permute!(times, p), permute!(data, p))
end

# pattern match on the known output of jl_static_show
extract_topmod(e::QuoteNode) = extract_topmod(e.value)
function extract_topmod(e)
    Meta.isexpr(e, :.) &&
        return extract_topmod(e.args[1])
    Meta.isexpr(e, :call) && length(e.args) == 3 && e.args[1] == :getfield &&
        return extract_topmod(e.args[2])
    Meta.isexpr(e, :call) && length(e.args) == 2 && e.args[1] == :typeof &&
        return extract_topmod(e.args[2])
    #Meta.isexpr(e, :call) && length(e.args) == 2 && e.args[1] == :Symbol &&
    #    return Symbol(e.args[2])
    # parametrized anonymous functions
    Meta.isexpr(e, :call) && e.args[1].args[1] == :getfield &&
        return extract_topmod(e.args[1].args[2].args[2])
    isa(e, Symbol) &&
        return e
    return :unknown
end

const anonrex = r"##\d*#{1,2}\d*"   # detect anonymous functions

function parse_call(line; subst=Vector{Pair{String, String}}(), blacklist=String[])
    match(anonrex, line) === nothing || return false, line, :unknown
    for (k, v) in subst
        line = replace(line, k=>v)
    end
    if any(b -> occursin(b, line), blacklist)
        println(line, " contains a blacklisted substring")
        return false, line, :unknown
    end

    curly = ex = Meta.parse(line, raise=false)
    while Meta.isexpr(curly, :where)
        curly = curly.args[1]
    end
    if !Meta.isexpr(curly, :curly)
        @warn("failed parse of line: ", line)
        return false, line, :unknown
    end
    func = curly.args[2]
    topmod = (func isa Expr ? extract_topmod(func) : :Main)

    # make some substitutions to try to form a leaf types tuple
    changed = false
    for i in 3:length(curly.args)
        e = curly.args[i]
        if Meta.isexpr(e, :where) && length(e.args) == 2 && Meta.isexpr(e.args[1], :curly) && length(e.args[1].args) == 3 && e.args[1].args[1] === :Vararg
            # Unwrap a Varargs argument, make it a single instance of the allowed Varargs type.
            e = e.args[1].args[2]
            curly.args[i] = e
            changed = true
        end
        if e === :Function
            # `Function` is not an allowable type for `precompile()`
            e = :(typeof(identity))
        elseif e === :Any
            # `Any` is an abstract type, and can't be precompiled for. Replace it with something that can.
            e = :Int
        elseif e == :(Type{T} where T)
            # Type is an abstract type; replace it with a concrete type.
            e = :(Type{Int})
        else
            continue
        end
        curly.args[i] = e
        changed = true
    end
    # In Julia 0.6, functions with symbolic names like `Base.:(+)` don't output correctly.
    # So if we didn't change the arg list above, use the original string.
    if changed
        line = string(ex)
    end
    return true, line, topmod
end

"""
`pc = parcel(calls; subst=[], blacklist=[])` assigns each compile statement to the module that owns the function. Perform string substitution via `subst=["Module1"=>"Module2"]`, and omit functions in particular modules with `blacklist=["Module3"]`. On output, `pc[:Module2]` contains all the precompiles assigned to `Module2`.

Use `SnoopCompile.write(prefix, pc)` to generate a series of files in directory `prefix`, one file per module.
"""
function parcel(calls; subst=Vector{Pair{String, String}}(), blacklist=String[])
    pc = Dict{Symbol, Vector{String}}()
    for c in calls
        local keep, pcstring, topmod
        keep, pcstring, topmod = parse_call(c, subst=subst, blacklist=blacklist)
        keep || continue
        # Add to the appropriate dictionary
        if !haskey(pc, topmod)
            pc[topmod] = String[]
        end
        push!(pc[topmod], "precompile($pcstring)")
    end
    return pc
end

"""
`pc = format_userimg(calls; subst=[], blacklist=[])` generates precompile directives intended for your base/userimg.jl script. Use `SnoopCompile.write(filename, pc)` to create a file that you can `include` into `userimg.jl`.
"""
function format_userimg(calls; subst=Vector{Pair{String, String}}(), blacklist=String[])
    pc = Vector{String}()
    for c in calls
        keep, pcstring, topmod = parse_call(c, subst=subst, blacklist=blacklist)
        keep || continue
        push!(pc, "precompile($pcstring)")
    end
    return pc
end

# Write precompiles for userimg.jl
function write(io::IO, pc::Vector)
    for ln in pc
        println(io, ln)
    end
end

function write(filename::AbstractString, pc::Vector)
    path, fn = splitdir(filename)
    if !isdir(path)
        mkpath(path)
    end
    open(filename, "w") do io
        write(io, pc)
    end
    nothing
end

"""
    write(prefix::AbstractString, pc::Dict; always::Bool = false)

Write each modules' precompiles to a separate file.  If `always` is
true, the generated function will always run the precompile statements
when called, otherwise the statements will only be called during
package precompilation.
"""
function write(prefix::AbstractString, pc::Dict; always::Bool = false)
    if !isdir(prefix)
        mkpath(prefix)
    end
    for (k, v) in pc
        open(joinpath(prefix, "precompile_$k.jl"), "w") do io
            println(io, "function _precompile_()")
            !always && println(io, "    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing")
            for ln in v
                println(io, "    ", ln)
            end
            println(io, "end")
        end
    end
end

macro snoop(args...)
    @warn "@snoop is deprecated, use @snoopc instead"
    :(@snoopc $(args...))
end

if VERSION >= v"1.2.0-DEV.573"
    function __init__()
        # typeinf_ext_timed must be compiled before it gets run
        # We do this in __init__ to make sure it gets compiled to native code
        # (the *.ji file stores only the inferred code)
        precompile(typeinf_ext_timed, (Core.MethodInstance, Core.Compiler.Params))
        precompile(typeinf_ext_timed, (Core.MethodInstance, UInt))
        nothing
    end
end

end # module


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

struct PrivateLinkResourcesApi <: SwaggerApi
    client::Swagger.Client
end

"""

Gets the private link resources that need to be created for a storage account.
Param: resourceGroupName::String (required)
Param: accountName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: PrivateLinkResourceListResult
"""
function _swaggerinternal_privateLinkResourcesListByStorageAccount(_api::PrivateLinkResourcesApi, resourceGroupName::String, accountName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    Swagger.validate_param("resourceGroupName", "privateLinkResourcesListByStorageAccount", :maxLength, resourceGroupName, 90)
    Swagger.validate_param("resourceGroupName", "privateLinkResourcesListByStorageAccount", :minLength, resourceGroupName, 1)

    Swagger.validate_param("accountName", "privateLinkResourcesListByStorageAccount", :maxLength, accountName, 24)
    Swagger.validate_param("accountName", "privateLinkResourcesListByStorageAccount", :minLength, accountName, 3)

    Swagger.validate_param("api_version", "privateLinkResourcesListByStorageAccount", :minLength, api_version, 1)

    Swagger.validate_param("subscriptionId", "privateLinkResourcesListByStorageAccount", :minLength, subscriptionId, 1)

    _ctx = Swagger.Ctx(_api.client, "GET", PrivateLinkResourceListResult, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Storage/storageAccounts/{accountName}/privateLinkResources", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "accountName", accountName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function privateLinkResourcesListByStorageAccount(_api::PrivateLinkResourcesApi, resourceGroupName::String, accountName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_privateLinkResourcesListByStorageAccount(_api, resourceGroupName, accountName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function privateLinkResourcesListByStorageAccount(_api::PrivateLinkResourcesApi, response_stream::Channel, resourceGroupName::String, accountName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_privateLinkResourcesListByStorageAccount(_api, resourceGroupName, accountName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

export privateLinkResourcesListByStorageAccount

ext: jl


text:
"""
$(TYPEDEF)

Store both grid and density for KDE over ``ℝ²``.

Reading the fields directly is part of the API, and

```julia
sum(density) * step(x) ≈ 1
```

# Fields

$(FIELDS)
"""
mutable struct UnivariateKDE{R<:AbstractRange} <: AbstractKDE
    "Gridpoints for evaluating the density."
    x::R
    "Kernel density at corresponding gridpoints `x`."
    density::Vector{Float64}
end

# construct kernel from bandwidth
kernel_dist(::Type{Normal},w::Real) = Normal(0.0,w)
kernel_dist(::Type{Uniform},w::Real) = (s = 1.7320508075688772*w; Uniform(-s,s))

const LocationScale = Union{Laplace,Logistic,SymTriangularDist}
kernel_dist(::Type{D},w::Real) where {D} = (s = w/std(D(0.0,1.0)); D(0.0,s))

# Silverman's rule of thumb for KDE bandwidth selection
function default_bandwidth(data::RealVector, alpha::Float64 = 0.9)
    # Determine length of data
    ndata = length(data)
    ndata <= 1 && return alpha

    # Calculate width using variance and IQR
    var_width = std(data)
    q25, q75 = quantile(data, [0.25, 0.75])
    quantile_width = (q75 - q25) / 1.34

    # Deal with edge cases with 0 IQR or variance
    width = min(var_width, quantile_width)
    if width == 0.0
        if var_width == 0.0
            width = 1.0
        else
            width = var_width
        end
    end

    # Set bandwidth using Silverman's rule of thumb
    return alpha * width * ndata^(-0.2)
end

function default_weights(data::RealVector)
    UniformWeights(length(data))
end

# Roughly based on:
#   B. W. Silverman (1982) "Algorithm AS 176: Kernel Density Estimation Using
#   the Fast Fourier Transform", Journal of the Royal Statistical
#   Society. Series C (Applied Statistics) , Vol. 31, No. 1, pp. 93-99
#   URL: http://www.jstor.org/stable/2347084
# and:
#   M. C. Jones and H. W. Lotwick (1984) "Remark AS R50: A Remark on Algorithm
#   AS 176. Kernal Density Estimation Using the Fast Fourier Transform",
#   Journal of the Royal Statistical Society. Series C (Applied Statistics) ,
#   Vol. 33, No. 1, pp. 120-122
#   URL: http://www.jstor.org/stable/2347674

# default kde range
# Should extend enough beyond the data range to avoid cyclic correlation from the FFT
function kde_boundary(data::RealVector, bandwidth::Real)
    lo, hi = extrema(data)
    lo - 4.0*bandwidth, hi + 4.0*bandwidth
end

# convert boundary and npoints to Range object
function kde_range(boundary::Tuple{Real,Real}, npoints::Int)
    lo, hi = boundary
    lo < hi || error("boundary (a,b) must have a < b")

    range(lo, stop=hi, length=npoints)
end

struct UniformWeights{N} end

UniformWeights(n) = UniformWeights{n}()

Base.sum(x::UniformWeights) = 1.
Base.getindex(x::UniformWeights{N}, i) where {N} = 1/N

const Weights = Union{UniformWeights, RealVector, StatsBase.Weights}

# tabulate data for kde
function tabulate(data::RealVector, midpoints::R, weights::Weights=default_weights(data)) where R<:AbstractRange
    npoints = length(midpoints)
    s = step(midpoints)

    # Set up a grid for discretized data
    grid = zeros(Float64, npoints)
    ainc = 1.0 / (sum(weights)*s*s)

    # weighted discretization (cf. Jones and Lotwick)
    for (i,x) in enumerate(data)
        k = searchsortedfirst(midpoints,x)
        j = k-1
        if 1 <= j <= npoints-1
            grid[j] += (midpoints[k]-x)*ainc*weights[i]
            grid[k] += (x-midpoints[j])*ainc*weights[i]
        end
    end

    # returns an un-convolved KDE
    UnivariateKDE(midpoints, grid)
end

# convolve raw KDE with kernel
# TODO: use in-place fft
function conv(k::UnivariateKDE, dist::UnivariateDistribution)
    # Transform to Fourier basis
    K = length(k.density)
    ft = rfft(k.density)

    # Convolve fft with characteristic function of kernel
    # empirical cf
    #  = \sum_{n=1}^N e^{i*t*X_n} / N
    #  = \sum_{k=0}^K e^{i*t*(a+k*s)} N_k / N
    #  = e^{i*t*a} \sum_{k=0}^K e^{-2pi*i*k*(-t*s*K/2pi)/K} N_k / N
    #  = A * fft(N_k/N)[-t*s*K/2pi + 1]
    c = -twoπ/(step(k.x)*K)
    for j = 0:length(ft)-1
        ft[j+1] *= cf(dist,j*c)
    end

    dens = irfft(ft, K)
    # fix rounding error.
    for i = 1:K
        dens[i] = max(0.0,dens[i])
    end

    # Invert the Fourier transform to get the KDE
    UnivariateKDE(k.x, dens)
end

# main kde interface methods
function kde(data::RealVector, weights::Weights, midpoints::R, dist::UnivariateDistribution) where R<:AbstractRange
    k = tabulate(data, midpoints, weights)
    conv(k,dist)
end

function kde(data::RealVector, dist::UnivariateDistribution;
             boundary::Tuple{Real,Real}=kde_boundary(data,std(dist)), npoints::Int=2048, weights=default_weights(data))

    midpoints = kde_range(boundary,npoints)
    kde(data,weights,midpoints,dist)
end

function kde(data::RealVector, midpoints::R;
             bandwidth=default_bandwidth(data), kernel=Normal, weights=default_weights(data)) where R<:AbstractRange
    bandwidth > 0.0 || error("Bandwidth must be positive")
    dist = kernel_dist(kernel,bandwidth)
    kde(data,weights,midpoints,dist)
end

function kde(data::RealVector; bandwidth=default_bandwidth(data), kernel=Normal,
             npoints::Int=2048, boundary::Tuple{Real,Real}=kde_boundary(data,bandwidth), weights=default_weights(data))
    bandwidth > 0.0 || error("Bandwidth must be positive")
    dist = kernel_dist(kernel,bandwidth)
    kde(data,dist;boundary=boundary,npoints=npoints,weights=weights)
end

# Select bandwidth using least-squares cross validation, from:
#   Density Estimation for Statistics and Data Analysis
#   B. W. Silverman (1986)
#   sections 3.4.3 (pp. 48-52) and 3.5 (pp. 61-66)

function kde_lscv(data::RealVector, midpoints::R;
                  kernel=Normal,
                  bandwidth_range::Tuple{Real,Real}=(h=default_bandwidth(data); (0.25*h,1.5*h)),
                  weights=default_weights(data)) where R<:AbstractRange

    ndata = length(data)
    k = tabulate(data, midpoints, weights)

    # the ft here is K/ba*sqrt(2pi) * u(s), it is K times the Yl in Silverman's book
    K = length(k.density)
    ft = rfft(k.density)

    ft2 = abs2.(ft)
    c = -twoπ/(step(k.x)*K)
    hlb, hub = bandwidth_range

    opt = Optim.optimize(hlb, hub) do h
        dist = kernel_dist(kernel, h)
        ψ = 0.0
        for j = 1:length(ft2)-1
            ks = real(cf(dist, j*c))
            ψ += ft2[j+1]*(ks-2.0)*ks
        end
        ψ*step(k.x)/K + pdf(dist,0.0)/ndata
    end

    dist = kernel_dist(kernel, Optim.minimizer(opt))
    for j = 0:length(ft)-1
        ft[j+1] *= cf(dist, j*c)
    end

    dens = irfft(ft, K)
    # fix rounding error.
    for i = 1:K
        dens[i] = max(0.0,dens[i])
    end

    # Invert the Fourier transform to get the KDE
    UnivariateKDE(k.x, dens)
end

function kde_lscv(data::RealVector;
                  boundary::Tuple{Real,Real}=kde_boundary(data,default_bandwidth(data)),
                  npoints::Int=2048,
                  kernel=Normal,
                  bandwidth_range::Tuple{Real,Real}=(h=default_bandwidth(data); (0.25*h,1.5*h)),
                  weights::Weights = default_weights(data))

    midpoints = kde_range(boundary,npoints)
    kde_lscv(data,midpoints; kernel=kernel, bandwidth_range=bandwidth_range, weights=weights)
end


text:
using Soss, Stheno
using DataFrames, CSV

births1 = CSV.read("/home/chad/git/fivethirtyeight_data/births/US_births_1994-2003_CDC_NCHS.csv")
births2 = CSV.read("/home/chad/git/fivethirtyeight_data/births/US_births_2000-2014_SSA.csv")

using Plots

gpc = GPC()
k7 = kernel(EQ(); l=7, s=0.9)

m = @model x begin
    α ~ Normal(9,5)
    β ~ Normal()
    ϕ ~ Exponential() # overdispersion
    f7 = GP(x -> α + β*x, k7, gpc)
    f = f7
    μ ~ f(x, ϕ)
    y .~ Poisson.(exp.(μ))
end

x = [1:50;];

rand(m(x=x))

# logdensityof(m(x=x), (;α=9.0, β=0.1, y=births2.births))

post = dynamicHMC(m(x=x), (;y=births2.births[1:50]))

k = kernel(EQ(); l=7, s=0.9);
f = GP(k, gpc)([1:7;])

f([0.1],4)


text:
using LinearAlgebra

import Tulip

include("LPBenchmark.jl")

"""
    run_tulip

Run benchmark for Tulip.
"""
function run_tulip(fname::String, fsol::String="")

    BLAS.set_num_threads(1)

    tulip = Tulip.Optimizer()

    # Read file
    load_problem!(tulip, fname)
    MOI.set(tulip, MOI.Silent(), false)
    MOI.set(tulip, MOI.NumberOfThreads(), 1)
    MOI.set(tulip, MOI.TimeLimitSec(), 10_000)
    MOI.set(tulip, MOI.RawParameter("BarrierIterationsLimit"), 500)

    # Solve instance
    run(tulip)
    return nothing
end

if abspath(PROGRAM_FILE) == @__FILE__
    run_tulip(ARGS[1])
end

max_stars_repo_name: mtanneau/LPBenchmarks

max_issues_repo_name: mtanneau/LPBenchmarks

max_forks_repo_name: mtanneau/LPBenchmarks


text:
using RandomNumbers.Xorshifts
using Distributions

@testset "bernoulli_process.jl" begin
    p = [0.1, 0.1]
    n = 10

    rng = Xorshift128Plus(1)

    bs = PMDPs.CountingProcesses.BernoulliScheme(n, p)

    @test params(bs) == (10, [0.1,0.1])
    @test params(bs[3])==([0.1,0.1,0.8],)

    samples = Matrix{Float64}(undef, 2,1000)
    for i in 1:1000
        samples[:,i] = rand(rng, bs)
    end
    max_shift, max_outcome =  maximum(samples; dims=2)
    min_shift, min_outcome =  minimum(samples; dims=2)
    @test min_outcome == 1
    @test max_outcome == 2
    @test min_shift == 1

    nhbs = PMDPs.CountingProcesses.NonhomogenousBernoulliScheme([0.1 0.1;0.1 0.1])[1]
end


text:
struct PDESystem <: ModelingToolkit.AbstractSystem
  eq
  bcs
  domain
  indvars
  depvars
end


text:
@testset "unname" begin
    for orig in ([1 2; 3 4], spzeros(2, 2))
        @test unname(NamedDimsArray(orig, (:x, :y))) === orig
        @test unname(orig) === orig
    end
end

@testset "dimnames" begin
    nda = NamedDimsArray([10 20; 30 40], (:x, :y))

    @test dimnames(nda) === (:x, :y)
    @test dimnames(nda, 2) === :y
    @test dimnames(nda, 3) === :_

    @test dimnames([10 20; 30 40]) === (:_, :_)
    @test dimnames([10 20; 30 40], 2) === :_
    @test dimnames([10 20; 30 40], 3) === :_

    v = NamedDimsArray(1:2, :a)
    @test dimnames(v, 2) == dimnames(permutedims(v), 1) # That's why :_ for d > ndims

    @test_throws Exception dimnames(nda, 0)
end

@testset "refine_names" begin
    @testset "Named Array into a NamedDimsArray" begin
        nda = refine_names(ones(3, 4, 5), (:a, :b, :c))
        @test dimnames(nda) == (:a, :b, :c)
        @test nda isa NamedDimsArray
    end

    @testset "Functioning on a fully named NamedDimsArray" begin
        orig_full = NamedDimsArray(ones(3, 4, 5), (:a, :b, :c))
        @test dimnames(refine_names(orig_full, (:a, :b, :c))) == (:a, :b, :c)
        @test dimnames(refine_names(orig_full, (:a, :b, :_))) == (:a, :b, :c)
        @test_throws DimensionMismatch refine_names(orig_full, (:a, :b, :wrong))
        @test_throws DimensionMismatch refine_names(orig_full, (:c, :a, :b))
    end

    @testset "Functioning on a partially named  NamedDimsArray" begin
        orig_partial = NamedDimsArray(ones(3, 4, 5), (:a, :_, :c))
        @test dimnames(refine_names(orig_partial, (:a, :b, :c))) == (:a, :b, :c)
        @test dimnames(refine_names(orig_partial, (:a, :_, :c))) == (:a, :_, :c)
    end
end

@testset "rename" begin
    nda = NamedDimsArray{(:a, :b, :c, :d)}(ones(10, 1, 1, 20))
    new_nda = rename(nda, (:w, :x, :y, :z))

    @test dimnames(new_nda) === (:w, :x, :y, :z)
    @test parent(new_nda) === parent(nda)
end


text:
using ASCIIPlots

function labelnames()
  # Set File Name
  DataDir = joinpath(dirname(@__FILE__),"data","bin")
  LabelNameFile = joinpath(DataDir,"batches.meta.txt")

  # Settings
  Labels = 10

  # Open file
  h_LabelNameFile = open(LabelNameFile)

  # Read Labels
  LabelNames = readlines(h_LabelNameFile)

  # Close File
  close(h_LabelNameFile)

  # Strip newlines
  for i=1:length(LabelNames)
    LabelNames[i] = LabelNames[i][1:end-1]
  end
  LabelNames=LabelNames[1:end-1]

  return LabelNames
end

function normalize(X)
    X = convert(Array{Float64},X)

    samples = size(X,2)

    for i=1:samples
      x = X[:,i]
      minx = minimum(x)
      maxx = maximum(x)
      ranx = maxx-minx

      X[:,i] = (x - minx) ./ ranx
    end

    return X
end

function toLuminance(dataset)
  features = size(dataset,1)
  pixels = convert(Int,features/3)
  # Get Color Channels
  r = dataset[1:pixels,:]
  g = dataset[pixels+1:2*pixels,:]
  b = dataset[2*pixels+1:end,:]
  # Convert to floating point
  r = convert(Array{Float64},r)
  g = convert(Array{Float64},g)
  b = convert(Array{Float64},b)
  # Apply Luminance Projection
  lum = 0.2126*r .+ 0.7152*g .+ 0.0722*b;

  return lum
end

function traindata(batch_number::Int; normalize_images=false, grey=false)
  # Assert to the number of batches
  @assert minimum(batch_number) >= 1 && maximum(batch_number) <= 5

  # Set File Name
  DataDir = joinpath(dirname(@__FILE__),"data","bin")
  BatchFile = @sprintf "data_batch_%d.bin" batch_number
  BatchFile = joinpath(DataDir,BatchFile)
  LabelNameFile = joinpath(DataDir,"batches.meta.txt")

  # Fixed dataset values
  Pixels = 1024
  Features = 3*Pixels
  Samples = 10000
  Labels = 10

  # Initialize Dataset Structures
  Dataset = Array(UInt8,Features,Samples)
  Labels = Array(UInt8,Samples,1)
  LabelNames = Array(AbstractString,10,1)

  # Open Batch File
  h_BatchFile = open(BatchFile)

  # Loop Over Samples
  r = Array(UInt8,Pixels,1)
  g = Array(UInt8,Pixels,1)
  b = Array(UInt8,Pixels,1)
  for sampleIdx = 1:Samples
    # Read Label Byte
    Labels[sampleIdx] = read(h_BatchFile,UInt8)+1
    # Read Color Byte Arrays
    read!(h_BatchFile,r)
    read!(h_BatchFile,g)
    read!(h_BatchFile,b)

    Dataset[1:Pixels,sampleIdx] = r
    Dataset[Pixels+1:2*Pixels,sampleIdx] = g
    Dataset[2*Pixels+1:end,sampleIdx] = b
  end

  # Close File
  close(h_BatchFile)

  # Get Label Names
  LabelNames = labelnames()

  # Conversion to Int type
  Dataset = convert(Array{Int},Dataset)
  Labels = convert(Array{Int},Labels)

  # Check for grey
  if grey
    Dataset = toLuminance(Dataset)
  end

  # Apply Normalization
  if normalize_images
    Dataset = normalize(Dataset)
  end

  return Dataset, Labels, LabelNames
end

function testdata(;normalize_images=false,grey=false)
  DataDir = joinpath(dirname(@__FILE__),"data","bin")
  BatchFile = "test_batch.bin"
  BatchFile = joinpath(DataDir,BatchFile)

  # Fixed dataset values
  Pixels = 1024
  Features = 3*Pixels
  Samples = 10000
  Labels = 10

  # Initialize Dataset Structures
  Dataset = Array(UInt8,Features,Samples)
  Labels = Array(UInt8,Samples,1)
  LabelNames = Array(AbstractString,10,1)

  # Open Batch File
  h_BatchFile = open(BatchFile)

  # Loop Over Samples
  r = Array(UInt8,Pixels,1)
  g = Array(UInt8,Pixels,1)
  b = Array(UInt8,Pixels,1)
  for sampleIdx = 1:Samples
    # Read Label Byte
    Labels[sampleIdx] = read(h_BatchFile,UInt8)+1
    # Read Color Byte Arrays
    read!(h_BatchFile,r)
    read!(h_BatchFile,g)
    read!(h_BatchFile,b)
    # Set Features in Dataset
    Dataset[1:Pixels,sampleIdx] = r
    Dataset[Pixels+1:2*Pixels,sampleIdx] = g
    Dataset[2*Pixels+1:end,sampleIdx] = b
  end

  # Close File
  close(h_BatchFile)

  # Get Label Names
  LabelNames = labelnames()

  # Conversion to Int type
  Dataset = convert(Array{Int},Dataset)
  Labels = convert(Array{Int},Labels)

  # Check for grey
  if grey
    Dataset = toLuminance(Dataset)
  end

  # Apply Normalization
  if normalize_images
    Dataset = normalize(Dataset)
  end

  return Dataset, Labels, LabelNames
end


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

struct DndApi <: SwaggerApi
    client::Swagger.Client
end

function _swaggerinternal_dndEndDnd(_api::DndApi; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "POST", Dict{String, Any}, "/dnd.endDnd", ["slackAuth"])
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/x-www-form-urlencoded", "application/json"] : [_mediaType])
    return _ctx
end

"""

Ends the current user's Do Not Disturb session immediately.
Return: Dict{String, Any}
"""
function dndEndDnd(_api::DndApi; _mediaType=nothing)
    _ctx = _swaggerinternal_dndEndDnd(_api; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function dndEndDnd(_api::DndApi, response_stream::Channel; _mediaType=nothing)
    _ctx = _swaggerinternal_dndEndDnd(_api; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

function _swaggerinternal_dndEndSnooze(_api::DndApi; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "POST", Dict{String, Any}, "/dnd.endSnooze", ["slackAuth"])
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/x-www-form-urlencoded", "application/json"] : [_mediaType])
    return _ctx
end

"""

Ends the current user's snooze mode immediately.
Return: Dict{String, Any}
"""
function dndEndSnooze(_api::DndApi; _mediaType=nothing)
    _ctx = _swaggerinternal_dndEndSnooze(_api; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function dndEndSnooze(_api::DndApi, response_stream::Channel; _mediaType=nothing)
    _ctx = _swaggerinternal_dndEndSnooze(_api; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

function _swaggerinternal_dndInfo(_api::DndApi; user=nothing, _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", Dict{String, Any}, "/dnd.info", ["slackAuth"])
    Swagger.set_param(_ctx.query, "user", user)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/x-www-form-urlencoded"] : [_mediaType])
    return _ctx
end

"""

Retrieves a user's current Do Not Disturb status.
Param: user::String
Return: Dict{String, Any}
"""
function dndInfo(_api::DndApi; user=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_dndInfo(_api; user=user, _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function dndInfo(_api::DndApi, response_stream::Channel; user=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_dndInfo(_api; user=user, _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

function _swaggerinternal_dndSetSnooze(_api::DndApi, num_minutes::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "POST", Dict{String, Any}, "/dnd.setSnooze", ["slackAuth"])
    Swagger.set_param(_ctx.form, "num_minutes", num_minutes)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/x-www-form-urlencoded"] : [_mediaType])
    return _ctx
end

"""

Turns on Do Not Disturb mode for the current user, or changes its duration.
Param: num_minutes::String (required)
Return: Dict{String, Any}
"""
function dndSetSnooze(_api::DndApi, num_minutes::String; _mediaType=nothing)
    _ctx = _swaggerinternal_dndSetSnooze(_api, num_minutes; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function dndSetSnooze(_api::DndApi, response_stream::Channel, num_minutes::String; _mediaType=nothing)
    _ctx = _swaggerinternal_dndSetSnooze(_api, num_minutes; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

function _swaggerinternal_dndTeamInfo(_api::DndApi; users=nothing, _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", Dict{String, Any}, "/dnd.teamInfo", ["slackAuth"])
    Swagger.set_param(_ctx.query, "users", users)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/x-www-form-urlencoded"] : [_mediaType])
    return _ctx
end

"""

Retrieves the Do Not Disturb status for up to 50 users on a team.
Param: users::String
Return: Dict{String, Any}
"""
function dndTeamInfo(_api::DndApi; users=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_dndTeamInfo(_api; users=users, _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function dndTeamInfo(_api::DndApi, response_stream::Channel; users=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_dndTeamInfo(_api; users=users, _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

export dndEndDnd, dndEndSnooze, dndInfo, dndSetSnooze, dndTeamInfo


text:
module LeafAreaIndex

using Optim: optimize, minimizer
using StatsBase: mean, midpoints
import FileIO, JLD2, Netpbm, Parameters, StatsBase, Optim

# import StatsBase #for `midpoints`
import ColorTypes #in crops.jl for AbstractRGB

# in calibration.jl
using Graphics: Point, norm
using DataFrames

export rawblueread, rawcolourread,
    CameraLens, CameraLensParams, PolarImage,
    Slope, SlopeParams, Mask, MaskParams,
    RedMax,#pixels, gapfraction, contactfreqs,
    threshold, #EdgeDetection, MinimumThreshold, RidlerCalvard,
    inverse, #Zenith57, Miller, Lang, EllipsLUT, EllipsOpt,
    langxiang #chencihlar,
    #calibrate_center, calibrate_projfun

include("Types.jl")
include("thresholding.jl")
include("segments.jl")
include("gapfraction.jl")
include("inversion.jl")
include("clumping.jl")
include("crops.jl")
include("ChenCihlar.jl")
include("calibration.jl")
include("io.jl")

# reoccuring argument checks
function checkθ1θ2(θ1, θ2)
    θ1 < 0   && throw(DomainError())
    θ2 > π/2 && throw(DomainError())
    @assert θ1 < θ2
end

#TODO replace with OnlineStats.jl
"Specialed type for immutable streaming sum, based on PR #18 from StreamStats.jl"
struct StreamMean
    streamsum::Float64
    len::Int
end
StreamMean() = StreamMean(0.0, 0)
update(sm::StreamMean, term) = StreamMean(sm.streamsum + term, sm.len + 1)
StatsBase.mean(sm::StreamMean) = sm.streamsum / sm.len
Base.empty!(sm::StreamMean) = StreamMean()
Base.length(sm::StreamMean) = sm.len

# A fast histogram method derived from julialang PR #8952. Used in thresholding
# and slope adjustment.
function fasthist(img::AbstractVector, edg::AbstractRange)
    n = length(edg) - 1
    histcount = zeros(Int, n)
    step(edg) <= 0 && error("step(edg) must be positive")
    for pixel in img
        f = (pixel - first(edg)) / step(edg)
        if 0 < f <= n
            histcount[ceil(Int, f)] += 1
        end
    end
    histcount
end

end # module


text:
#=
Created on 30/01/2021 16:39:07
Last update: -

@author: Michiel Stock
michielfmstock@gmail.com

Generates the figures for flatland.
=#

include("examples/flatland.jl")

# examples of shapes

square = Square((1, 1), l=2)
rectangle = Rectangle((2, 6), l=1, w=2.5)
circle = Circle((-2, 1), R=1)
pentagon = RegularPolygon((-3, 4), 5, R=1.2)
hexagon = RegularPolygon((-3.4, 7), 5, R=1.5, θ=0.7)
triangle = Triangle((-4, 1), (-1, -1), (-8, 1) )

plot(square, label="square", legend=:topleft)
plot!(rectangle, label="rectangle")
plot!(circle, label="circle")
plot!(pentagon, label="pentagon")
plot!(hexagon, label="hexagon")
plot!(triangle, label="triangle")
title!("Some shapes")
savefig("../img/boundingbox.svg")
savefig("../img/boundingbox.png")

# outline of regular polygon

plot(hexagon, legend=:topleft)
corn_hex = corners(hexagon)
scatter!(first.(corn_hex), last.(corn_hex), label="corners")
plot!([first(center(hexagon)), corn_hex[end][1]], [last(center(hexagon)), corn_hex[end][2]],
                label="radius R", lw=2, color="orange")
plot!([first(center(hexagon)), corn_hex[1][1]], [last(center(hexagon)), corn_hex[1][2]],
    label="angle theta", lw=2, color="red", ls=:dash)
hline!([7], lw=2, color="red", ls=:dash, label="")

# bounding box

plot(triangle)
vline!(xlim(triangle), label="xlims")
hline!(ylim(triangle), label="ylims")
plot!(boundingbox(triangle), ls=:dash, color="red", label="bounding box", fillcolor=nothing)

# in/out

points = [8rand(2) for i in 1:500]
pinrect = filter(p->p ∈ rectangle, points)
pnotinrect = filter(p->p ∉ rectangle, points)
plot(rectangle)
scatter!(first.(pinrect), last.(pinrect), label="point in shape")
scatter!(first.(pnotinrect), last.(pnotinrect), label="point not in shape")

# overlap

t1 = Triangle((0,0), (1, 0.4), (0.5, 0.7))
t2 = deepcopy(t1)
rotate!(t2, π/2)
move!(t2, (1, 0.2))
t3 = deepcopy(t1)
rotate!(t3, π/3)
move!(t3, (0.4, 0.2))

plot(t1, alpha=0.5)
plot!(t2, alpha=0.5)
plot!(t3, alpha=0.5)
plot!(boundingbox(t1), ls=:dash, color="red", fillcolor=nothing)
plot!(boundingbox(t2), ls=:dash, color="red", fillcolor=nothing)
plot!(boundingbox(t3), ls=:dash, color="red", fillcolor=nothing)
title!("Overlaping bounding boxes is a required,
    but not sufficient condition for overlap")
#title!("Overlaping bounding boxes is a required, but not sufficient condition for overlap")
savefig("../img/boundingbox.svg")
savefig("../img/boundingbox.png")

# examples

shapes1, trials1 = rejection_sampling(triangle, 50, (0, 100), (0, 100))
plotshapes(shapes1, title="50 triangles\n $trials1 trials")

small_large_squares = [Rectangle((0, 0), l=10, w=15) for i in 1:5]
append!(small_large_squares, [Rectangle((0, 0), l=1, w=1.5) for i in 1:25])
shapes2, trials2 = rejection_sampling!(small_large_squares, (0, 100), (0, 100))
plotshapes(small_large_squares, title="5 large and 25 small squares\n $trials2 trials")


text:
"""
	Reduces the set of possible children of a node using the grammar's constraints
"""
function propagate_constraints(grammar::ContextSensitiveGrammar, context::GrammarContext, child_rules::Vector{Int})
	domain = child_rules

	for propagator in grammar.constraints
		domain = propagate(propagator, context, domain)
	end

	return domain
end

"""
	reimplementation of the ExprRules._next_state!
	Change: child expressions are filtered so that the constraints are not violated
"""
function _next_state!(node::RuleNode, grammar::ContextSensitiveGrammar, max_depth::Int, context::GrammarContext)

	if max_depth < 1
	    return (node, false) # did not work
	elseif isterminal(grammar, node)
	    # do nothing
	    if iseval(grammar, node.ind) && (node._val === nothing)  # evaluate the rule
		node._val = eval(grammar.rules[node.ind].args[2])
	    end
	    return (node, false) # cannot change leaves
	else # !isterminal
	    # if node is not terminal and doesn't have children, expand every child
	    if isempty(node.children)
		if max_depth ≤ 1
		    return (node,false) # cannot expand
		end

		child_index = 1  # keepa track of which child we are processing now (needed for context)

		# build out the node
		for c in child_types(grammar, node)
		    worked = false
		    i = 0
		    child = RuleNode(0)

		    new_context = GrammarContext(context.originalExpr, deepcopy(context.nodeLocation))
		    push!(new_context.nodeLocation, child_index)

		    child_rules = [x for x in grammar[c]]  # select all applicable rules
		    child_rules = propagate_constraints(grammar, new_context, child_rules)  # filter out those that violate constraints

		    while !worked && i < length(child_rules)
			i += 1
			child = RuleNode(child_rules[i])

			if iseval(grammar, child.ind) # if rule needs to be evaluated (_())
			    child._val = eval(grammar.rules[child.ind].args[2])
			end

			worked = true
			if !isterminal(grammar, child)
			    child, worked = _next_state!(child, grammar, max_depth-1, new_context)
			end
		    end
		    if !worked
			return (node, false) # did not work
		    end
		    push!(node.children, child)

		    child_index += 1
		end

		return (node, true)
	    else # not empty
		# make one change, starting with rightmost child
		worked = false
		child_index = length(node.children) + 1
		while !worked && child_index > 1
		    child_index -= 1
		    child = node.children[child_index]

		    new_context = GrammarContext(context.originalExpr, deepcopy(context.nodeLocation))
		    push!(new_context.nodeLocation, child_index)

		    child, child_worked = _next_state!(child, grammar, max_depth-1, new_context)
		    while !child_worked
			child_type = return_type(grammar, child)

			child_rules = [x for x in grammar[child_type]]  # get all applicable rules
			child_rules = propagate_constraints(grammar, new_context, child_rules)  # filter ones that violate constraints

			i = something(findfirst(isequal(child.ind), child_rules), 0)
			if i < length(child_rules)
			    child_worked = true
			    child = RuleNode(child_rules[i+1])

			    # node needs to be evaluated
			    if iseval(grammar, child.ind)
				child._val = eval(grammar.rules[child.ind].args[2])
			    end

			    if !isterminal(grammar, child)
				child, child_worked = _next_state!(child, grammar, max_depth-1, new_context)
			    end
			    node.children[child_index] = child
			else
			    break
			end
		    end

		    if child_worked
			worked = true

			# reset remaining children
			for child_index2 in child_index+1 : length(node.children)
			    c = child_types(grammar, node)[child_index2]
			    worked = false
			    i = 0
			    child = RuleNode(0)

			    new_context = GrammarContext(context.originalExpr, deepcopy(context.nodeLocation))
			    push!(new_context.nodeLocation, child_index2)

			    child_rules = [x for x in grammar[c]]  # take all applicable rules
			    child_rules = propagate_constraints(grammar, new_context, child_rules)  # reomove ones that violate constraints

			    while !worked && i < length(child_rules)
				i += 1
				child = RuleNode(child_rules[i])

				if iseval(grammar, child.ind)
				    child._val = eval(grammar.rules[child.ind].args[2])
				end

				worked = true
				if !isterminal(grammar, child)
				    child, worked = _next_state!(child, grammar, max_depth-1, new_context)
				end
			    end
			    if !worked
				break
			    end
			    node.children[child_index2] = child
			end
		    end
		end

		return (node, worked)
	    end
	end
    end

mutable struct ContextSensitiveIterator <: ExpressionIterator
	grammar::ContextSensitiveGrammar
	max_depth::Int
	sym::Symbol
end

function Base.iterate(iter::ContextSensitiveIterator)
	init_node = RuleNode(0)  # needed for propagating constraints on the root node
	init_context = GrammarContext(init_node)

	grammar, sym, max_depth = iter.grammar, iter.sym, iter.max_depth

	# propagate constraints on the root node
	sym_rules = [x for x in grammar[sym]]
	sym_rules = propagate_constraints(grammar, init_context, sym_rules)
	#node = RuleNode(grammar[sym][1])
	node = RuleNode(sym_rules[1])

	if isterminal(grammar, node)
	    return (deepcopy(node), node)
	else
	    context = GrammarContext(node)
	    node, worked =  _next_state!(node, grammar, max_depth, context)
	    while !worked
		# increment root's rule
		rules = [x for x in grammar[sym]]
		rules = propagate_constraints(grammar, init_context, rules) # propagate constraints on the root node

		i = something(findfirst(isequal(node.ind), rules), 0)
		if i < length(rules)
		    node, worked = RuleNode(rules[i+1]), true
		    if !isterminal(grammar, node)
			node, worked = _next_state!(node, grammar, max_depth, context)
		    end
		else
		    break
		end
	    end
	    return worked ? (deepcopy(node), node) : nothing
	end
end

function Base.iterate(iter::ContextSensitiveIterator, state::RuleNode)
	grammar, max_depth = iter.grammar, iter.max_depth
	context = GrammarContext(state)
	node, worked = _next_state!(state, grammar, max_depth, context)

	while !worked
	    # increment root's rule
	    init_node = RuleNode(0)  # needed for propagating constraints on the root node
    	    init_context = GrammarContext(init_node)

	    rules = [x for x in grammar[iter.sym]]
	    rules = propagate_constraints(grammar, init_context, rules)

	    i = something(findfirst(isequal(node.ind), rules), 0)
	    if i < length(rules)
		node, worked = RuleNode(rules[i+1]), true
		if !isterminal(grammar, node)
		    context = GrammarContext(node)
		    node, worked = _next_state!(node, grammar, max_depth, context)
		end
	    else
		break
	    end
	end
	return worked ? (deepcopy(node), node) : nothing
end


text:
using AutoClutch
using Test

@testset "AutoClutch.jl" begin
    # Write your tests here.
end

max_stars_repo_name: powerzf/AutoClutch

max_issues_repo_name: powerzf/AutoClutch

max_forks_repo_name: powerzf/AutoClutch


text:
readstring("/devel/myfile.txt") # read file into a string
read("/devel/myfile.txt", filesize("/devel/myfile.txt")) # read file into an array of bytes

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# # Tutorials for Trixi.jl

# The tutorial section for [Trixi.jl](https://github.com/trixi-framework/Trixi.jl) also contains
# interactive step-by-step explanations via [Binder](https://mybinder.org).

# Right now, you are using the classic documentation. The corresponding interactive notebooks can
# be opened in [Binder](https://mybinder.org/) and viewed in [nbviewer](https://nbviewer.jupyter.org/)
# via the icons ![](https://mybinder.org/badge_logo.svg) and ![](https://raw.githubusercontent.com/jupyter/design/master/logos/Badges/nbviewer_badge.svg)
# in the respective tutorial.
# You can download the raw notebooks from GitHub via ![](https://camo.githubusercontent.com/aea75103f6d9f690a19cb0e17c06f984ab0f472d9e6fe4eadaa0cc438ba88ada/68747470733a2f2f696d672e736869656c64732e696f2f62616467652f646f776e6c6f61642d6e6f7465626f6f6b2d627269676874677265656e).

# **Note:** To improve responsiveness via caching, the notebooks are updated only once a week. They are only
# available for the latest stable release of Trixi at the time of caching.

# There are tutorials for the following topics:

# ### [1 Introduction to DG methods](@ref scalar_linear_advection_1d)
#-
# This tutorial gives an introduction to discontinuous Galerkin (DG) methods with the example of the
# scalar linear advection equation in 1D. Starting with some theoretical explanations, we first implement
# a raw version of a discontinuous Galerkin spectral element method (DGSEM). Then, we will show how
# to use features of Trixi.jl to achieve the same result.

# ### [2 DGSEM with flux differencing](@ref DGSEM_FluxDiff)
#-
# To improve stability often the flux differencing formulation of the DGSEM (split form) is used.
# We want to present the idea and formulation on a basic 1D level. Then, we show how this formulation
# can be implemented in Trixi and analyse entropy conservation for two different flux combinations.

# ### [3 Shock capturing with flux differencing and stage limiter](@ref shock_capturing)
#-
# Using the flux differencing formulation, a simple procedure to capture shocks is a hybrid blending
# of a high-order DG method and a low-order subcell finite volume (FV) method. We present the idea on a
# very basic level and show the implementation in Trixi. Then, a positivity preserving limiter is
# explained and added to an exemplary simulation of the Sedov blast wave with the 2D compressible Euler
# equations.

# ### 4 Adding a new equation
# #### [4.1 Scalar conservation law](@ref cubic_conservation_law)
#-
# This tutorial explains how to add a new physics model using the example of the cubic conservation
# law. First, we define the equation using a `struct` `CubicEquation` and the physical flux. Then,
# the corresponding standard setup in Trixi.jl (`mesh`, `solver`, `semi` and `ode`) is implemented
# and the ODE problem is solved by OrdinaryDiffEq's `solve` method.

# #### [4.2 Nonconservative advection](@ref nonconservative_advection)
#-
# In this part, another physics model is implemented, the nonconservative linear advection equation.
# We run two different simulations with different levels of refinement and compare the resulting errors.

# ### [5 Differentiable programming](@ref differentiable_programming)
#-
# This part deals with some basic differentiable programming topics. For example, a Jacobian, its
# eigenvalues and a curve of total energy (through the simulation) are calculated and plotted for
# a few semidiscretizations. Moreover, we calculate an example for propagating errors with Measurement.jl
# at the end.

# ### [6 Unstructured meshes with HOHQMesh.jl](@ref hohqmesh_tutorial)
#-
# The purpose of this tutorial is to demonstrate how to use the [`UnstructuredMesh2D`](@ref)
# functionality of Trixi.jl. This begins by running and visualizing an available unstructured
# quadrilateral mesh example. Then, the tutorial will demonstrate how to conceptualize a problem
# with curved boundaries, generate a curvilinear mesh using the available [HOHQMesh](https://github.com/trixi-framework/HOHQMesh)
# software in the Trixi.jl ecosystem, and then run a simulation using Trixi.jl on said mesh.

# ## Examples in Trixi
# Trixi already contains several more coding examples, the so-called `elixirs`. You can find them
# in the folder [`examples`](https://github.com/trixi-framework/Trixi.jl/blob/main/examples/).
# They are structured by the underlying mesh type and the respective number of spatial dimensions.
# The name of an elixir is composed of the underlying system of conservation equations (for instance
# `advection` or `euler`) and other special characteristics like the initial condition
# (e.g. `gauss`, `astro_jet`, `blast_wave`) or the included simulation features
# (e.g. `amr`, `shockcapturing`).


text:
"""
    Server(project, builder)

A server manages interactively running a `Project`, coordinating events
to update project state efficiently.
"""
mutable struct Server
    project
    builder
    updates
    lock
end

Server(project, builder = FileBuilder(HTML(), Path(mktempdir()))) = Server(project, builder, Updates(), ReentrantLock())

"""
    abstract type ServerMode

Mode that a `Server` can run in. Controls how events affect the project
state.
"""
abstract type ServerMode end

geteventsource(::ServerMode, server, ch) = nothing
geteventsource(::Rewriter, ch) = nothing

initialize(::ServerMode, server) = return
"""
    handle(server, mode, event)
"""
function handle end

"""
    runserver(server, mode)

Run a `server` in `mode`. Handles start and cleanup of event sources
from rewriters and `mode`. Synchronizes updates to project state.
"""
function runserver(server, mode; dt=1 / 60)
    eventch = Channel()
    eventsources = servereventsources(server, mode, eventch)
    initialize(mode, server)
    @async begin
        for event in eventch
            @debug "Received $(typeof(event))"
            try
            handle(server, mode, event)
            catch e
                @error "Error was thrown during handling of event!" event=event error=e
            end
        end
    end
    try
        @info "Starting server..."
        start(eventsources)
        while true
            updates = server.updates
            lock(server.lock) do
                server.updates = Updates()
            end
            applyupdates!(server.project, server.builder, updates)
            if !all(isempty.((updates.sources, updates.torewrite, updates.torebuild)))
                @info "Rebuild complete."
            end
            sleep(dt)
        end
    catch e
        if e isa InterruptException
            @info "Shutting down server..."
        else
            rethrow(e)
        end
    finally
        stop(eventsources)
        close(eventch)
    end
end

#=
[`Updates`](#) is a container that asynchronously collects updates to
the project state which are then applied snychronously. For example,
a file watcher can update source documents using [`addsource!`](#) and
trigger rewrites and builds using [`addrewrite!`](#) and [`addbuild!`].
Only in [`applyupdates!`](#) are the changes actually applied to a project. =#

"""
    struct Updates

Tracks updates to project state. Used internally by [`runserver`](#).
"""
struct Updates
    sources
    torewrite
    torebuild
end

Updates() = Updates(Dict{AbstractPath,XTree}(), Set{AbstractPath}(), Set{AbstractPath}())

"""
    applyupdates!(project, builder, updates::Updates)

Apply `updates` to a project by updating sources, rewriting and building
documents.
"""
function applyupdates!(project, builder, updates::Updates)
    for (p, doc) in updates.sources
        project.sources[p] = doc
    end

    dirtypaths = if !isempty(updates.torewrite)
        rewritesources!(project, updates.torewrite)
    else
        Set()
    end
    dirtypaths = union(dirtypaths, updates.torebuild)

    if !isempty(dirtypaths)
        build(builder, project, dirtypaths)
    end

    return project
end

applyupdates!(server::Server) = applyupdates!(server.project, server.builder, server.updates)

function addsource!(server, path, doc)
    lock(server.lock) do
        server.updates.sources[path] = doc
    end
end

function addrewrite!(server, path)
    lock(server.lock) do
        if haskey(server.project.sources, path) || haskey(server.updates.sources, path)
            push!(server.updates.torewrite, path)
        else
            error("Cannot find source document $path to rewrite!")
        end
    end
end

function addbuild!(server, path)
    lock(server.lock) do
        if haskey(server.project.outputs, path) || path in server.updates.torewrite
            push!(server.updates.torebuild, path)
        else
            error("Cannot find output document $path to build!")
        end
    end
end

function servereventsources(server, mode, ch)
    eventsources = []
    for rewriter in server.project.rewriters
        push!(eventsources, geteventsource(rewriter, ch))
    end
    push!(eventsources, geteventsource(mode, server, ch))
    return collect(filter(!isnothing, eventsources))
end


text:
module cindex

export parse_header, cu_type, ty_kind, name, spelling,
       value, children, cu_file, resolve_type, return_type,
       tokenize, pointee_type, typedef_type
export CLType, CLCursor, CXString, CXTypeKind, CursorList, TokenList
export getindex, start, next, done, search, show, endof

export function_return_modifiers, function_arg_modifiers

import Base.getindex, Base.start, Base.next, Base.done, Base.search, Base.show
import Base.endof, Base.length

using Compat

###############################################################################

# Name of the helper library
import ..libwci

###############################################################################

include("cindex/defs.jl")
include("cindex/types.jl")
include("cindex/base.jl")

###############################################################################

# Main entry point for parsing
# Returns root CXCursor in the TranslationUnit for a given header
#
# Required argument:
#   "header.h"          header file to parse
#
# Optional (keyword) arguments:
#   index:              CXIndex pointer (pass to avoid re-allocation)
#   diagnostics:        Display Clang diagnostics
#   cplusplus:          Parse as C++
#   includes:           Vector of extra include directories to search
#   args:               Compiler switches as string array, eg: ["-x", "c++", "-fno-elide-type"]
#   flags:              Bitwise OR of TranslationUnitFlags
#
function parse_header(header::String;
                index                           = None,
                diagnostics::Bool               = true,
                cplusplus::Bool                 = false,
                args                            = ASCIIString[""],
                includes                        = ASCIIString[],
                flags                           = TranslationUnit_Flags.None)
    if !isfile(header)
        error(header, " not found")
    end
    if (index == None)
        index = idx_create(0, (diagnostics ? 1 : 0))
    end
    if (cplusplus)
        append!(args, ["-x", "c++"])
    end
    if (length(includes) > 0)
        args = vcat(args, [["-I",x] for x in includes]...)
    end

    tu = tu_parse(index, header, args, length(args),
                  C_NULL, 0, flags)
    if (tu == C_NULL)
        error("ParseTranslationUnit returned NULL; unable to create TranslationUnit")
    end

    return tu_cursor(tu)
end

###############################################################################
# Search function for CursorList
# Returns vector of CXCursors in CursorList matching predicate
#
# Required arguments:
#   CursorList      List to search
#   IsMatch(CXCursor)
#                   Predicate Function, accepting a CXCursor argument
#
function search(cl::CursorList, ismatch::Function)
    ret = CLCursor[]
    for cu in cl
        ismatch(cu) && push!(ret, cu)
    end
    ret
end
search(cu::CLCursor, ismatch::Function) = search(children(cu), ismatch)
search(cu::CLCursor, T::DataType) = search(cu, x->isa(x, T))
search(cu::CLCursor, name::ASCIIString) = search(cu, x->(cindex.spelling(x) == name))

###############################################################################
# Extended search function
# Returns a Dict{ DataType => CLCursor

function matchchildren(cu::CLCursor, types::Array{DataType,1})
    ret = @compat Dict{Any,Any}([ t => CLCursor[] for t in types])
    for child in children(cu)
        for t in types
            isa(child, t) && push!(ret[t], child)
        end
    end
    return ret
end

###############################################################################

# TODO: macro version should be more efficient.
anymatch(first, args...) = any([==(first, a) for a in args])

cu_type(c::CLCursor) = getCursorType(c)
ty_kind(c::CLType) = convert(Int, c.data[1].kind)
name(c::CLCursor) = getCursorDisplayName(c)
spelling(c::CLType) = getTypeKindSpelling(convert(Int32, ty_kind(c)))
spelling(c::CLCursor) = getCursorSpelling(c)
is_null(c::CLCursor) = (Cursor_isNull(c) != 0)

function resolve_type(rt::CLType)
    # This helper attempts to work around some limitations of the
    # current libclang API.
    if ty_kind(rt) == cindex.TypeKind.Unexposed
        # try to resolve Unexposed type to cursor definition.
        rtdef_cu = cindex.getTypeDeclaration(rt)
        if (!is_null(rtdef_cu) && !isa(rtdef_cu, NoDeclFound))
            return cu_type(rtdef_cu)
        end
    end
    # otherwise, this will either be a builtin or unexposed
    # client needs to sort out.
    return rt
end

function return_type(c::Union(FunctionDecl, CXXMethod), resolve::Bool)
    if (resolve)
        return resolve_type( getCursorResultType(c) )
    else
        return getCursorResultType(c)
    end
end
return_type(c::CLCursor) = return_type(c, true)

function pointee_type(t::Pointer)
    return cindex.getPointeeType(t)
end
pointee_type(cu::CLCursor) = error("pointee_type(CLCursor) is discontinued, please use pointee_type(cindex.cu_type(cu))")

function typedef_type(c::TypedefDecl)
    t = cindex.getTypedefDeclUnderlyingType(c)
    if isa(t, Unexposed)
        t = cu_type(children(c)[1])
    end
    return t
end

function value(c::EnumConstantDecl)
    t = cu_type(c)
    if isa(t, IntType) || isa(t, Long) || isa(t, LongLong)
            return getEnumConstantDeclValue(c)
    elseif isa(t, UInt) || isa(t, ULong) || isa(t, ULongLong)
            return getEnumConstantDeclUnsignedValue(c)
    end
    error("Unknown EnumConstantDecl type: ", t, " cursor: ", c)
end

function getindex(cl::CursorList, clid::Int, default::UnionType)
    try
        getindex(cl, clid)
    catch
        return default
    end
end
function getindex(cl::CursorList, clid::Int)
    if (clid < 1 || clid > cl.size) error("Index out of range or empty list") end
    cu = TmpCursor()
    ccall( (:wci_getCLCursor, libwci),
            Void,
            (Ptr{Void}, Ptr{Void}, Int),
            cu.data, cl.ptr, clid-1)
    return CXCursor(cu)
end

function children(cu::CLCursor)
    cl = cl_create()
    ccall(  (:wci_getChildren, libwci),
            Ptr{Void},
            (Ptr{CXCursor}, Ptr{Void}),
            cu.data, cl.ptr)
    size = cl_size(cl.ptr)
    return CursorList(cl.ptr,size)
end

function cu_file(cu::CLCursor)
    str = CXString()
    ccall( (:wci_getCursorFile, libwci),
            Void,
            (Ptr{Void}, Ptr{Void}),
            cu.data, str.data)
    return get_string(str)
end

start(cl::CursorList) = 1
done(cl::CursorList, i) = (i > cl.size)
next(cl::CursorList, i) = (cl[i], i+1)
length(cl::CursorList) = cl.size

################################################################################
# Tokenizer access
################################################################################

# Returns TokenList
function tokenize(cursor::CLCursor)
    tu = Cursor_getTranslationUnit(cursor)
    sourcerange = getCursorExtent(cursor)
    return cindex.tokenize(tu, sourcerange)
end

length(tl::TokenList) = tl.size
start(tl::TokenList) = 1
next(tl::TokenList, i) = (tl[i], i+1)
endof(tl::TokenList) = length(tl)
done(tl::TokenList, i) = (i > length(tl))

function getindex(tl::TokenList, i::Int)
    if (i < 1 || i > length(tl)) throw(BoundsError()) end

    c = CXToken(unsafe_load(tl.ptr, i))
    kind = c.data[1].int_data1
    spelling = cindex.getTokenSpelling(tl.tunit, c)

    if (kind == TokenKind.Punctuation)
        return Punctuation(spelling)
    elseif (kind == TokenKind.Keyword)
        return Keyword(spelling)
    elseif (kind == TokenKind.Identifier)
        return Identifier(spelling)
    elseif (kind == TokenKind.Literal)
        return Literal(spelling)
    elseif (kind == TokenKind.Comment)
        return Comment(spelling)
    end
end

################################################################################
# Helper functions
################################################################################

# Retrieve function arguments for a given cursor
function function_args(cursor::Union(FunctionDecl, CXXMethod))
    search(cursor, ParmDecl)
end

#returns a tuple with the default argument values for a C++ function
#only seems to work if cplusplus=true in parse_header
function function_arg_defaults(method::Union(cindex.CXXMethod, cindex.FunctionDecl, cindex.Constructor))
    defvals = Any[]
    for c in children(method)
        if isa(c,cindex.ParmDecl)

            #value in case default was not specified
            lit = nothing
            toks = tokenize(c)

            n = 1
            #get index of '='
            for tok in toks
                if isa(tok, cindex.Punctuation) && tok.text == "="
                    break
                end
                n += 1
            end

            ts = collect(toks)

            #default value is the one after '='
            if n<length(toks)
                lit = join(map(x->x.text, ts[n+1:end]), "")
            end

            #parse default value as Julia code
            pl = lit != nothing ? parse(lit) : lit

            #for empty string, need to add quotes
            if pl == ""
                pl = "\"\""
            end
            push!(defvals, pl)
        end
    end
    return tuple(defvals...)
end

#returns an array of the function return type modifiers
function function_return_modifiers(f::Union(FunctionDecl, CXXMethod))
    modifs = Any[]
    for tok in tokenize(f)
        if isa(tok, cindex.Keyword)
            if tok.text in ["const"]
                push!(modifs, tok)
            end
        end
        isa(tok, cindex.Identifier) && tok.text==spelling(f) && break
    end
    return modifs
end

function function_arg_modifiers(p::ParmDecl)
    modifs = cindex.Keyword[]
    for tok in tokenize(p)

        #only read keywords
        isa(tok, cindex.Keyword) || continue

        #read up to variable identifier
        isa(tok, cindex.Identifier) && break

        const tt = tok.text
        if tt in ["const", "virtual"]
            push!(modifs, tok)
        end
    end
    return modifs
end

################################################################################
# Display overrides
################################################################################

show(io::IO, tk::CLToken)   = print(io, typeof(tk), "(\"", tk.text, "\")")
show(io::IO, ty::CLType)    = print(io, "CLType (", typeof(ty), ") ")
show(io::IO, cu::CLCursor)    = print(io, "CLCursor (", typeof(cu), ") ", name(cu))

###############################################################################
# Internal functions
###############################################################################

tu_init(hdrfile::Any) = tu_init(hdrfile, 0, false, 0)
function tu_init(hdrfile::Any, diagnostics, cpp::Bool, opts::Int)
    idx = idx_create(0,diagnostics)
    tu = tu_parse(idx, hdrfile, (cpp ? ["-x", "c++"] : [""]), opts)
    return tu
end

tu_dispose(tu::CXTranslationUnit) = ccall( (:clang_disposeTranslationUnit, "libclang"), Void, (Ptr{Void},), tu)

function tu_cursor(tu::CXTranslationUnit)
    if (tu == C_NULL)
        error("Invalid TranslationUnit!")
    end
    getTranslationUnitCursor(tu)
end

function tu_parse(CXIndex,
                  source_filename::ASCIIString,
                  cl_args::Array{ASCIIString,1},
                  num_clargs,
                  unsaved_files::CXUnsavedFile,
                  num_unsaved_files,
                  options)

    ccall(  (:clang_parseTranslationUnit, "libclang"),
            CXTranslationUnit,
            (Ptr{Void}, Ptr{Uint8}, Ptr{Ptr{Uint8}}, Uint32, Ptr{Void}, Uint32, Uint32),
                CXIndex,
                source_filename,
                cl_args,
                num_clargs,
                unsaved_files,
                num_unsaved_files,
                options)
end

function idx_create(excludeDeclsFromPCH::Int, displayDiagnostics::Int)
    ccall(  (:clang_createIndex, "libclang"),
            CXTranslationUnit,
            (Int32, Int32),
                excludeDeclsFromPCH,
                displayDiagnostics)
end
idx_create() = idx_create(0,0)

function getFile(tu::CXTranslationUnit, file::ASCIIString)
    ccall( (:clang_getFile, "libclang"),
            CXFile,
            (Ptr{Void}, Ptr{Uint8}),
                tu,
                file)
end

function cl_create()
    ptr = ccall( (:wci_createCursorList, libwci),
                Ptr{Void},
                () )
    return CursorList(ptr,0)
end

function cl_dispose(cl::CursorList)
    ccall( (:wci_disposeCursorList, libwci),
            None,
            (Ptr{Void},),
                cl.ptr)
end

function cl_size(clptr::Ptr{Void})
    ccall( (:wci_sizeofCursorList, libwci),
            Int,
            (Ptr{Void},),
                clptr)
end
cl_size(cl::CursorList) = cl.size

end # module


text:
"""
    struct ProfilePoint
        value::Float64               # x value of profile point
        loss::Float64                # y value of profile point (loss function value)
        params::Array{Float64, 1}    # vector of optimal values of `loss_func` arguments
        ret::Symbol                  # `NLOpt.optimize()` return value
        counter::Union{Int, Nothing} # number of `loss_func` evaluations
    end

Structure storing one point from profile function.
ret values: `:FORCED_STOP`, `:MAXEVAL_REACHED`, `:FTOL_REACHED`

"""
struct ProfilePoint
    value::Float64
    loss::Float64
    params::Array{Float64, 1}
    ret::Symbol
    counter::Union{Int, Nothing}
end

"""
    struct EndPoint
        value::Union{Real, Nothing}           # value of endpoint or `nothing`
        profilePoints::Array{ProfilePoint, 1} # vector of profile points
        status::Symbol                        # status
        direction::Symbol                     # `:right` or `:left`
        counter::Int                          # number of `loss_func` evaluations
        supreme::Union{Real, Nothing}         # maximal value inside profile interval
    end

Structure storing one endpoint of confidence interval.

status values: `:BORDER_FOUND_BY_SCAN_TOL`, `:BORDER_FOUND_BY_LOSS_TOL`,
 `:SCAN_BOUND_REACHED`, `:MAX_ITER_STOP`, `:LOSS_ERROR_STOP`
"""
struct EndPoint
    value::Union{Real, Nothing}
    profilePoints::Array{ProfilePoint, 1}
    status::Symbol
    direction::Symbol
    counter::Int
    supreme::Union{Real, Nothing}
end


text:
using MathGL
mgl = MathGL

gr = mgl.Graph(800, 500)
mgl.xlabel(gr, "x")
mgl.ylabel(gr, "y")
mgl.box(gr)
mgl.axis(gr)
mgl.grid(gr)
mgl.plot(gr, 0.8*sin(linspace(-4pi, 4pi, 200)))
mgl.write(gr, "../graphs/sin_julia.png")


text:
#################################################################################################
# AIM 5: two-sample comparison performance
#################################################################################################
## Deps
using StatsPlots
using Distributed
using Combinatorics
using Distributions
using DelimitedFiles
@everywhere using Pkg
@everywhere Pkg.activate("./") # <=
@everywhere using CpelNano

## Source code
root_dir = dirname(dirname(pathof(CpelNano)))
include("$(root_dir)/calls/PaperSimulations/Aim-5/Aim-5-Source.jl")

## Constants

# IO
const data_dir = "/Users/jordiabante/OneDrive - Johns Hopkins/CpelNano/Data/Simulations/Aim-5/" # <=
const blind_friend_col = ["#999999","#E69F00","#56B4E9","#009E73","#F0E442","#0072B2","#D55E00","#CC79A7"]

# Set up
@everywhere const M = 5
@everywhere const POBS = 1.0
@everywhere const Δ_LEVEL_MEAN = 5.0
@everywhere const LEVEL_STD_M = 3.0 / sqrt(4000 / 450)
@everywhere const LEVEL_STD_U = 3.0 / sqrt(4000 / 450)
@everywhere const LEVEL_MEAN_M = Δ_LEVEL_MEAN / 2.0
@everywhere const LEVEL_MEAN_U = -Δ_LEVEL_MEAN / 2.0
@everywhere const EF_PERC_RNG = 0.05:0.05:0.25

## Default attributes
default(titlefont=(14, "arial"),guidefont=(16, "arial"),tickfont=(12, "arial"))

#################################################################################################
# Calls
#################################################################################################

# Arguments
eff_type = ARGS[1]
n_sim_reps = parse(Int64, ARGS[2])
CpelNano.print_log("Running ROC simulations for two-sample test $(eff_type) w/ $(n_sim_reps) reps")

# Run
eff_size_arr, tpr, fpr = run_sim(n_sim_reps, eff_type)

# Plot
plt = plt_roc(eff_size_arr, tpr, fpr, "ROC - $(eff_type) - two-sample test ")
savefig(plt, "$(data_dir)/ROC-Two-Sample-Test-$(eff_type).pdf")


text:
"Parameters for right-hand side function of Variational Partitioned Runge-Kutta methods."
const ParametersVPRKpSymmetric = AbstractParametersVPRK{:vprk_psymmetric}

"Variational partitioned Runge-Kutta integrator."
struct IntegratorVPRKpSymmetric{DT, TT, D, S,
                PT <: ParametersVPRKpSymmetric{DT,TT},
                ST <: NonlinearSolver{DT},
                IT <: InitialGuessIODE{TT}} <: AbstractIntegratorVPRKwProjection{DT,TT,D,S}

    params::PT
    solver::ST
    iguess::IT
    caches::CacheDict{PT}

    function IntegratorVPRKpSymmetric(params::ParametersVPRKpSymmetric{DT,TT,D,S}, solver::ST, iguess::IT, caches) where {DT,TT,D,S,ST,IT}
        new{DT, TT, D, S, typeof(params), ST, IT}(params, solver, iguess, caches)
    end

    function IntegratorVPRKpSymmetric{DT,D}(equations::NamedTuple, tableau::TableauVPRK{TT}, Δt::TT) where {DT,TT,D}
        # get number of stages
        S = tableau.s

        # create params
        R = convert(Vector{TT}, [1, tableau.R∞])
        params = ParametersVPRKpSymmetric{DT,D}(equations, tableau, Δt, NamedTuple{(:R,)}((R,)))

        # create cache dict
        caches = CacheDict(params)

        # create solver
        solver = create_nonlinear_solver(DT, D*(S+1), params, caches)

        # create initial guess
        iguess = InitialGuessIODE(get_config(:ig_extrapolation), equations[:v̄], equations[:f̄], Δt)

        # create integrator
        IntegratorVPRKpSymmetric(params, solver, iguess, caches)
    end

    function IntegratorVPRKpSymmetric(equation::Union{IODE{DT},LODE{DT}}, tableau, Δt; kwargs...) where {DT}
        IntegratorVPRKpSymmetric{DT, ndims(equation)}(get_functions(equation), tableau, Δt; kwargs...)
    end
end

function Integrators.get_internal_variables(int::IntegratorVPRKpSymmetric{DT,TT,D,S}) where {DT, TT, D, S}
    Q = create_internal_stage_vector(DT, D, S)
    P = create_internal_stage_vector(DT, D, S)
    V = create_internal_stage_vector(DT, D, S)
    F = create_internal_stage_vector(DT, D, S)
    λ = zeros(DT,D)

    solver = get_solver_status(int.solver)

    (Q=Q, P=P, V=V, F=F, λ=λ, solver=solver)
end

function initial_guess!(int::IntegratorVPRKpSymmetric{DT,TT}, sol::AtomicSolutionPODE{DT,TT},
                        cache::IntegratorCacheVPRK{DT}=int.caches[DT]) where {DT,TT}
    for i in eachstage(int)
        evaluate!(int.iguess, sol.q̄, sol.p̄, sol.v̄, sol.f̄,
                              sol.q, sol.p, sol.v, sol.f,
                              cache.q̃, cache.ṽ,
                              tableau(int).q.c[i])

        for k in eachdim(int)
            int.solver.x[ndims(int)*(i-1)+k] = cache.ṽ[k]
        end
    end

    for k in eachdim(int)
        int.solver.x[ndims(int)*nstages(int)+k] = 0
    end
end

function compute_projection_vprk!(x::Vector{ST},
                q::SolutionVector{ST}, p::SolutionVector{ST}, λ::SolutionVector{ST},
                Q::Vector{Vector{ST}}, V::Vector{Vector{ST}}, U::Vector{Vector{ST}}, G::Vector{Vector{ST}},
                params::ParametersVPRKpSymmetric{DT,TT,D,S}) where {ST,DT,TT,D,S}

    local t₀::TT = params.t̄
    local t₁::TT = params.t̄ + params.Δt
    local y1::ST
    local y2::ST

    # copy x to λ and q̄
    for k in 1:D
        λ[k] = x[D*S+k]
    end

    # compute U=λ
    U[1] .= λ
    U[2] .= λ

    # compute G=g(q,λ)
    for k in 1:D
        y1 = 0
        y2 = 0
        for j in 1:S
            y1 += params.tab.q.b[j] * V[j][k]
            y2 += params.tab.q.b̂[j] * V[j][k]
        end
        q[k] = params.q̄[k] + params.Δt * (y1 + y2) + params.Δt * (params.pparams[:R][1] * U[1][k] + params.pparams[:R][2] * U[2][k])
    end

    # compute G=g(q,λ)
    params.equ[:g](t₀, params.q̄, λ, G[1])
    params.equ[:g](t₁, q, λ, G[2])

    # compute p=ϑ(q)⋅p
    params.equ[:ϑ](t₁, q, λ, p)
end

"Compute stages of variational partitioned Runge-Kutta methods."
function Integrators.function_stages!(x::Vector{ST}, b::Vector{ST},
                params::ParametersVPRKpSymmetric{DT,TT,D,S},
                caches::CacheDict) where {ST,DT,TT,D,S}

    # get cache for internal stages
    cache = caches[ST]

    compute_stages!(x, cache.q̃, cache.p̃, cache.λ, cache.Q, cache.V, cache.U, cache.P, cache.F, cache.G, params)

    # compute b = - [P-AF-U]
    compute_rhs_vprk!(b, cache.P, cache.F, cache.G, params)

    # compute b = - [p-bF-G]
    compute_rhs_vprk_projection_p!(b, cache.p̃, cache.F, cache.G, D*(S+0), params)

    compute_rhs_vprk_correction!(b, cache.V, params)
end

function Integrators.integrate_step!(int::IntegratorVPRKpSymmetric{DT,TT}, sol::AtomicSolutionPODE{DT,TT},
                                     cache::IntegratorCacheVPRK{DT}=int.caches[DT]) where {DT,TT}
    # update nonlinear solver parameters from cache
    update_params!(int.params, sol)

    # compute initial guess
    initial_guess!(int, sol, cache)

    # reset solution
    reset!(sol, timestep(int))

    # call nonlinear solver
    solve!(int.solver)

    # print solver status
    print_solver_status(int.solver.status, int.solver.params)

    # check if solution contains NaNs or error bounds are violated
    check_solver_status(int.solver.status, int.solver.params)

    # compute vector fields at internal stages and projection vector fields
    compute_stages!(int.solver.x,
                    cache.q̃, cache.p̃, cache.λ,
                    cache.Q, cache.V, cache.U,
                    cache.P, cache.F, cache.G, int.params)

    # compute unprojected solution
    update_solution!(int, sol, cache)

    # add projection to solution
    project_solution!(int, sol, int.params.pparams[:R], cache)

    # copy solution to initial guess
    update_vector_fields!(int.iguess, sol.t, sol.q, sol.p, sol.v, sol.f)

    # copy internal stage variables
    sol.internal.Q .= cache.Q
    sol.internal.P .= cache.P
    sol.internal.V .= cache.V
    sol.internal.F .= cache.F
    sol.internal.λ .= cache.λ

    # copy solver status
    get_solver_status!(int.solver, sol.internal[:solver])
end


text:
"""
Unit tests for QuNet

Hudson's style convention for testing:

+ Redundency is good.
+ Don't aim to be redundant, but don't aim for tight scripting either.
+ If you want to manipulate a network from the network-library, (I.E. refreshing,
  adding temporal links, percolating etc.) do it on a copy.
+ One test set per file. Test every function in a given file thoroughly
+ Each test or subset of tests within the testset should begin with a comment
  briefly explaining what it's testing. Ideally, the test should work completely independently
  from other tests.
+ If you find something here that doesn't follow the style guide but it works regardless,
  don't change it.
+ Don't follow this style guide if you have a better one
+ Redundency is good

"""

using QuNet
using Test
using LightGraphs
using SimpleWeightedGraphs

include_networks = ["barbell", "simple_network", "simple_satnet",
"small_square", "shortest_path_test", "smalltemp", "greedy_test", "bridge"]

for N in include_networks
    include("network-library/$N.jl")
end

@testset "Network.jl" begin
    # Test: QNetwork is correctly initialised
    @test(typeof(barbell) == QNetwork)

    # Test: Nodes are correctly added
    @test length(barbell.nodes) == 2

    # Test: Channels are correctly added
    @test length(barbell.channels) == 1

    # Test: getnode works for id
    newnode = QuNet.getnode(barbell, 1)
    @test newnode == barbell.nodes[1]

    # Test: getnode works for name
    newnode = QuNet.getnode(barbell, "A")
    @test newnode == barbell.nodes[1]

    # Test: getchannel works for id
    newchannel = QuNet.getchannel(barbell, 1, 2)
    @test newchannel == barbell.channels[1]

    # Test: getchannel works for string
    newerchannel = QuNet.getchannel(barbell, "A", "B")
    @test newerchannel == barbell.channels[1]

    # Test: Update a sat network and check that the costs have changed
    AS = QuNet.getchannel(simple_satnet, "A", "S")
    old_costs = AS.costs
    update(simple_satnet)
    new_costs = AS.costs
    for key in keys(old_costs)
        @test old_costs[key] != new_costs[key]
    end

    # Test: Reset the network back to t=0 and check position goes back to init.
    S = QuNet.getnode(simple_satnet, "S")
    update(simple_satnet, 0.0)
    @test S.location.x == 500

    # Test: Check that deepcopy can clone network structure
    Q = deepcopy(barbell)
    @test all(Q.nodes[i] != barbell.nodes[i] for i in 1:length(Q.nodes))
    @test cmp(string(Q), string(barbell)) == 0

    # Test: refresh_graph! creates SimpleWeightedGraph copies of Network for all costs
    Q = deepcopy(barbell)
    QuNet.refresh_graph!(Q)
    @test length(Q.graph) == length(zero_costvector())
    g = Q.graph["Z"]
    @test nv(g) == 2
    @test ne(g) == 2
    @test g.weights[1, 2] == 0.5
    @test g.weights[2, 1] == 0.5

    # Test refresh_graph! on satellite network
    Q = deepcopy(simple_satnet)
    QuNet.refresh_graph!(Q)
    g = Q.graph["loss"]
    @test nv(g) == 3
    @test ne(g) == 4

    # Test 12: Test that update works on copied graph
    Q = deepcopy(barbell)
    update(Q)
    @test cmp(string(Q), string(barbell)) != 0

    # Test 13 / 14: Test that getchannel fetches the right channel in
    # a copied graph
    Q = deepcopy(barbell)
    AB = QuNet.getchannel(barbell, "A", "B")
    CAB = QuNet.getchannel(Q, "A", "B")
    @test (AB in barbell.channels) && (CAB in Q.channels)
    @test !(CAB in barbell.channels) && !(AB in Q.channels)
end

@testset "Channel.jl" begin
    #Test Basic Channel is initalised with costs
    AB = BasicChannel(A, B)
    @test typeof(AB) == BasicChannel

    # Test AirChannel is a subtype of QChannel
    AS = AirChannel(A, S)
    @test typeof(AS) <: QuNet.QChannel
end

@testset "Node.jl" begin
    # Test Basic node is initalised
    A = BasicNode("A")
    @test A.name == "A"

    # Test node properties can be updated
    B = BasicNode("B")
    B.location = Coords(100, 0, 0)
    @test isequal(B.location.x, 100)

    # Test PlanSatNode is initialised
    S = PlanSatNode("S")
    S.location = Coords(0, 0, 1000)
    S.velocity = Velocity(1000, 0, 0)
end

@testset "Utilities.jl" begin
    # Test maximal coherence
    Z = 1.0
    dB = Z_to_dB(Z)
    @test dB == 0
    Zn = dB_to_Z(dB)
    @test Zn == 1

    # Test maximal decoherence
    Z = 0.5
    dB = Z_to_dB(Z)
    @test dB == Inf

    Z = dB_to_Z(dB)
    @test Z == 0.5

    # Test decoherence < 0.5
    Z = 0.3
    @test_throws DomainError Z_to_dB(Z)

    # Test Purification
    cv1 = Dict()
    cv2 = Dict()
    cv1["loss"] = 1.0
    cv2["loss"] = 2.0
    cv1["Z"] = 1.0
    cv2["Z"] = 2.0
    result = purify([cv1, cv2])
    # TODO: Compare with hand calc
end

@testset "CostVector.jl" begin
    # Test: convert_costs to decibelic form
    cost_vector = Dict("loss"=>0.5, "Z"=>0.5)
    dB_cv = convert_costs(cost_vector, false)
    @test (cost_vector["loss"] == 3.0102999566398116 &&
    cost_vector["Z"]== Inf)

    # Test: convert_cost from decibelic form
    metric_cv = convert_costs(cost_vector, true)
    @test cost_vector == metric_cv

    # Test: get_pathcv for a vector of QChannels
    Q = QNetwork()
    A = BasicNode("A")
    S = PlanSatNode("S")
    S.location = Coords(0,0,1000)
    C = AirChannel(A, S)
    for i in [A, S, C]
        add(Q, i)
    end
    cost_vector = get_pathcv([C])
    @test(cost_vector["loss"] == 0.4970454276591223 &&
    cost_vector["Z"]== 0.49694603901995044)

    # Test: get_pathcv for a single QChannel
    Q = deepcopy(barbell)
    channel = Q.channels[1]
    cv = get_pathcv(channel)
    @test(cv["loss"] == 1 && cv["Z"] == 0.5)

    # Test: get_pathcv for a QNetwork, and a path of Int tuples
    Q = deepcopy(simple_network)
    path = [(1,2),(2,3)]
    pathcv = get_pathcv(Q, path)
    @test pathcv["loss"] == 2 && pathcv["Z"] == 2

    # Test: get_pathcv for a TemporalGraph and a path of Int tuples
    T = deepcopy(smalltemp)
    path = [(1,5)]
    pathcv = get_pathcv(T, path)
    @test pathcv["loss"] == 1.0 && pathcv["Z"] == 1.0
end

@testset "Routing.jl" begin
    # Test: shortest_path for a SimpleWeightedDiGraph
    g = SimpleWeightedDiGraph(2)
    add_edge!(g, 1, 2, 3.14)
    short_path = shortest_path(g, 1, 2)
    @test(short_path[1] == SimpleWeightedEdge(1, 2, 3.14))

    # Test: path_length
    @test(3.14 == QuNet.path_length(g, short_path))

    # Test: remove_shortest_path! for a SimpleWeightedDiGraph
    removed_path, removed_path_cost = QuNet.remove_shortest_path!(g, 1, 2)
    @test(removed_path_cost == 3.14)
    @test(length(edges(g)) == 0)

    # Test: remove_shortest_path! for a QNetwork
    Q = deepcopy(shortest_path_test)
    removed_path, removed_cv = QuNet.remove_shortest_path!(Q, "loss", 1, 2)
    # Test removed path is correct
    shortestpath = [(1,3),(3,2)]
    shortestpath = QuNet.int_to_simpleedge(shortestpath)
    @test(shortestpath == removed_path)

    # Test removed path costs are correct
    @test(removed_cv["Z"] == 10 && removed_cv["loss"] == 2)

    # Test that the correct path was removed for both weighted graphs of the QNetwork
    # And check that edges are removed in both directions
    G = Q.graph["loss"]
    @test has_edge(G, 1, 3) == false && has_edge(G, 3, 1) == false &&
    has_edge(G, 3, 2) == false && has_edge(G, 2, 3) == false

    # Check that no nodes and no other edges were removed
    @test nv(G) == length(Q.nodes)
    @test ne(G) == (length(Q.channels) - 2) * 2

    G = Q.graph["Z"]
    @test has_edge(G, 1, 3) == false && has_edge(G, 3, 1) == false &&
    has_edge(G, 3, 2) == false && has_edge(G, 2, 3) == false

    # Check that no nodes and no other edges were removed
    @test nv(G) == length(Q.nodes)
    @test ne(G) == (length(Q.channels) - 2) * 2

    # Test: remove_shortest_path! for a TemporalGraph returning cost vector
    # and using asynchronus nodes
    T = deepcopy(smalltemp)
    src = 1+8
    dst = 2+8
    QuNet.add_async_nodes!(T, [(src, dst)])
    removed_path, removed_cv = QuNet.remove_shortest_path!(T, "loss", src, dst)
    # Test removed path is correct
    shortestpath = [(1,2)]
    shortestpath = QuNet.int_to_simpleedge(shortestpath)
    @test(shortestpath == removed_path)
    # Test removed path cost vectors is correct.
    @test removed_cv["loss"] == 1 && removed_cv["Z"] == 1

    # Remove shortest path again, and test that the path in the next temporal layer was removed
    removed_path, removed_cv = QuNet.remove_shortest_path!(T, "loss", src, dst)
    shortestpath = [(5,6)]
    shortestpath = QuNet.int_to_simpleedge(shortestpath)
    @test(shortestpath == removed_path)

    # Remove shortest path one final time:
    # test that the path now takes a longer route in first temporal layer
    removed_path, removed_cv = QuNet.remove_shortest_path!(T, "loss", src, dst)
    shortestpath = [(1,3),(3,4),(4,2)]
    shortestpath = QuNet.int_to_simpleedge(shortestpath)
    @test(shortestpath == removed_path)
    # test that the costs of this path are correct
    @test(removed_cv["loss"] == 3.0 && removed_cv["Z"] == 3.0)

    # Test that the correct path was removed for both weighted graphs of the TempNet
    # And check that edges are removed in both directions
    G = T.graph["loss"]
    @test has_edge(G, 1, 3) == false && has_edge(G, 3, 1) == false &&
    has_edge(G, 3, 4) == false && has_edge(G, 4, 3) == false &&
    has_edge(G, 4, 2) == false && has_edge(G, 2, 4) == false

    G = T.graph["loss"]
    @test has_edge(G, 1, 3) == false && has_edge(G, 3, 1) == false &&
    has_edge(G, 3, 4) == false && has_edge(G, 4, 3) == false &&
    has_edge(G, 4, 2) == false && has_edge(G, 2, 4) == false

    # Test: Check that no nodes and no other channels are removed for TempNet
    T = deepcopy(smalltemp)
    QuNet.add_async_nodes!(T, [(1, 2)])
    removed_path, removed_cv = QuNet.remove_shortest_path!(T, "loss", 1, 2)
    QuNet.remove_async_nodes!(T)

    G = T.graph["loss"]
    @test nv(G) == 8
    @test ne(G) == 18

    G = T.graph["Z"]
    @test nv(G) == 8
    @test ne(G) == 18

    #Test: greedy_multi_path
    Q = deepcopy(greedy_test)
    pathset, pur_paths, pathuse_count = QuNet.greedy_multi_path!(Q, purify, [(1, 2)])
    # Test that all three paths were purified together.
    @test pathuse_count == [0, 0, 0, 1]

    # Test that the pathset contains all the paths from "A" to "B"
    all_paths = [[(1,3),(3,2)], [(1,4), (4,2)], [(1,5), (5,2)]]
    new_pathset = []
    for path in pathset[1]
        path = QuNet.simpleedge_to_int(path)
        push!(new_pathset, path)
    end
    @test new_pathset == all_paths

    # Test: greedy_multi_path on a TemporalGraph with 2 asynchronus end-users
    T = deepcopy(smalltemp)
    offset = smalltemp.nv * smalltemp.steps
    src1 = 1 + offset
    dst1 = 4 + offset
    src2 = 2 + offset
    dst2 = 3 + offset
    QuNet.add_async_nodes!(T, [(src1, dst1), (src2, dst2)])
    pathset, pur_paths, pathuse_count = QuNet.greedy_multi_path!(T, purify, [(src1, dst1), (src2, dst2)])
    # Test that 2 edge-disjoint paths were found for each usepair
    @test(pathuse_count[3] == 2)
    # Test that the pathsets have the same destinations
    for bundle in pathset
        dsts = []
        for path in bundle
            push!(dsts, last(path).dst)
        end
        @test(all(i == dsts[1] for i in dsts))
    end

    # # Test: A grid lattice completely saturated with end-users has no purification
    # # i.e. no end-user used 2 or more paths
    # NOTE: I actually proved myself wrong on this one. Purification of 2 and even 3
    # Paths is possible depending on the choice of userpair. Will leave the code here
    # As proof:
    # Q = GridNetwork(10, 10)
    # userpairs = make_user_pairs(Q, 50)
    # pathset, purpaths, pathuse_count = QuNet.greedy_multi_path!(Q, purify, userpairs)
    # @test pathuse_count[3] == 0
    # @test pathuse_count[4] == 0

    # # Test that Greedy_multi_path and other routing algorithms can distinguish between
    # # The asynchronous edges of a temporal graph (i.e. prefering short times over long ones)
    # Q = deepcopy(bridge)
    # T = QuNet.TemporalGraph(Q, 2)
    # # Make asynchronus userpairs
    # QuNet.add_async_nodes!(T, [(13,17), (14,18)])
    # pathset, pur_paths, pathuse_count = QuNet.greedy_multi_path!(T, purify, [(13,17), (14,18)])
    # # Expected output: [[(1-3),(3-4),(4-5)],(2-8),(8-9),(9-10),(10-12)]
    # println(pathset)
end

@testset "TemporalGraphs.jl" begin
    # Test: TemporalGraph initialises correctly
    G = GridNetwork(2, 2)
    T = QuNet.TemporalGraph(G, 2)
    @test nv(T.graph["loss"]) == 8
    @test T.nv == 4

    # Test: TemporalGraph attribute memory_prob works
    G = GridNetwork(10, 10)
    T = QuNet.TemporalGraph(G, 2, memory_prob=0.0)
    @test ne(T.graph["Z"]) == 720

    # Test: TemporalGraph adds 100 more edges when memory_prob = 100%
    T = QuNet.TemporalGraph(G, 2, memory_prob=1.0)
    @test ne(T.graph["Z"]) == 820

    # Test TemporalGraph with a default memory_cost, check that temporal edges have that cost.
    T = QuNet.TemporalGraph(G, 2, memory_prob=1.0, memory_costs=Dict("Z"=>3, "loss"=>4))
    graph = T.graph["Z"]
    @test graph.weights[101, 1] == 3
    graph = T.graph["loss"]
    @test graph.weights[104, 4] == 4

    # Test add_async_nodes! for endusers that are both temporal
    G = deepcopy(barbell)
    T = QuNet.TemporalGraph(G, 3)
    @test T.has_async_nodes == false
    async_pairs = [(7,8)]
    QuNet.add_async_nodes!(T, async_pairs)
    @test nv(T.graph["loss"]) == 8
    @test nv(T.graph["Z"]) == 8
    @test T.nv == 2
    @test T.has_async_nodes == true

    # Test add_async_nodes! won't let you double dip. Uncomment for good warning
    # QuNet.add_async_nodes!(T, async_pairs)

    # Test that async edges were added to the correct nodes and in the proper directions
    g = T.graph["loss"]
    # async_pairs = [(7,8)]
    # usage: g.weights[dst, src]
    @test g.weights[1,7] != 0 && g.weights[7,1] == 0
    @test g.weights[2,8] == 0 && g.weights[8,2] != 0
    @test g.weights[3,7] != 0 && g.weights[7,3] == 0
    @test g.weights[4,8] == 0 && g.weights[8,4] != 0
    @test g.weights[5,7] != 0 && g.weights[7,5] == 0
    @test g.weights[6,8] == 0 && g.weights[8,6] != 0

    # Test add_async_nodes! for endusers where src is fixed on the top plane and dst is asynchronus
    G = deepcopy(barbell)
    T = QuNet.TemporalGraph(G, 2)
    async_pairs = [(1, 6)]
    QuNet.add_async_nodes!(T, async_pairs)
    g = T.graph["Z"]
    # Test that no async links were added to src
    @test g.weights[1,5] == 0 && g.weights[5,1] == 0

    # # Test that async links were added to dst and in proper direction
    # # usage: g.weights[dst, src]
    @test g.weights[6, 2] != 0 && g.weights[2, 6] == 0
    @test g.weights[6, 4] != 0 && g.weights[4, 6] == 0

    # Test remove_async_nodes!
    G = deepcopy(barbell)
    T = QuNet.TemporalGraph(G, 2)
    async_pairs = make_user_pairs(T, 1)
    QuNet.add_async_nodes!(T, async_pairs)
    QuNet.remove_async_nodes!(T)
    @test nv(T.graph["loss"]) == 4
    @test T.has_async_nodes == false

    # Test that remove_async_nodes! won't let you double dip
    # uncomment for a useful warning
    # QuNet.remove_async_nodes!(T)

    # Test fix async_nodes_in_time
    Q = deepcopy(barbell)
    depth = 5
    T = QuNet.TemporalGraph(Q, depth)
    src = 1 + 2*depth
    dst = 2 + 2*depth
    user_pairs = [(src, dst)]
    QuNet.add_async_nodes!(T, user_pairs)
    # Fix the third node
    QuNet.fix_async_nodes_in_time!(T, [3])

    g = T.graph["loss"]
    users = filter(isodd, 1:depth)
    for user in users
        if user == 3
            @test has_edge(g, src, user) == true
        else
            @test has_edge(g, src, user) == false
        end
    end

    # TODO: Not working
    # # Test remove_async_edges! for specific time layers
    # Q = deepcopy(barbell)
    # T = QuNet.TemporalGraph(Q, 5)
    # user_pairs = [(11, 12)]
    # QuNet.add_async_nodes!(T, user_pairs)
    # QuNet.remove_async_edges!(T, 3)
    # # Async_edges should have been removed at t = 3
    # g = T.graph["loss"]
    # src = 11; dst = 12
    # #usage g.weights[dst, src]
    # @test g.weights[1, src] != 0
    # @test g.weights[3, src] != 0
    # @test g.weights[5, src] == 0
    # @test g.weights[7, src] != 0

    # @test(has_path(g, src, 1) == true)
    # @test(has_path(g, src, 3) == true)
    # @test(has_path(g, src, 5) == false)
    # @test(has_path(g, src, 7) == true)

end

@testset "Benchmarking.jl" begin
    # TODO Test dict_average
    # TODO Test dict_err

    # Test make_user_pairs for QNetwork.

    # NOTE: One important detail I noticed:
    # make_user_pairs ordinarily generates random user pairs. I verified this in
    # A seperate test file. Impressively, the @testset seems to fix the seed, such
    # that no matter how many times you roll the dice, you get the same outcome.
    # Convenient for unit testing, not so much for verifying randomness.

    # Test make_user_pairs for a regular QNetwork
    Q = GridNetwork(2,2)
    user_pair = [(4,3),(1,2)]
    @test user_pair == [(4,3),(1,2)]

    # Test make_user_pairs for a TemporalGraph, specifying that pairs should be asynchronus
    G = deepcopy(barbell)
    T = QuNet.TemporalGraph(G, 2)
    async_pairs = QuNet.make_user_pairs(T, 1, src_layer=-1, dst_layer=-1)
    @test async_pairs == [(5,6)]

    # Test make_user_pairs for a TemporalGraph, specifying that src should be asynchronus, and dst on layer 1
    G = deepcopy(barbell)
    T = QuNet.TemporalGraph(G, 3)
    lopsided_pair = QuNet.make_user_pairs(T, 1, src_layer=-1, dst_layer=1)
    @test lopsided_pair[1][1] > T.nv * T.steps && lopsided_pair[1][2] <= T.nv

    # Test make_user_pairs for a TemporalGraph for many pairs, specifying dst should be asynchronus and src on arbitrary layer,
    G = GridNetwork(5,5)
    T = QuNet.TemporalGraph(G, 5)
    lopsided_pairs = QuNet.make_user_pairs(T, 10, src_layer=3, dst_layer=-1)
    @test all(T.nv * 2 < lopsided_pairs[i][1] <= T.nv * 3 for i in 1:10)
    @test all(T.nv * T.steps < lopsided_pairs[i][2] for i in 1:10)

    # Test ave_paths_used
    pathuse_count = [5.0,1.0,2.0,4.0]
    # Expected answer:
    # (0*5 + 1*1 + 2*2 + 3*4) / (5 + 1 + 2 + 4) == 17/12
    ave_pathuse = QuNet.ave_paths_used(pathuse_count)
    @test(ave_pathuse == 17/12)

    # TODO Test net_performance
    # Test for 1 trial of barbell network
    # usage: (network::QNetwork, num_trials::Int64, num_pairs::Int64; max_paths=3)
    Q = deepcopy(barbell)
    performance, performance_err, ave_pathcounts, ave_pathcounts_err = net_performance(barbell, 1, 1)
    # Check ave_pathcounts is correct
    @test ave_pathcounts == [0, 1, 0, 0]
    # Only one sample taken, so no variance in ave_pathcounts:
    @test all(isnan(i) == true for i in ave_pathcounts_err)

    # Test net_performance for many trials of barbell network
    Q = deepcopy(barbell)
    performance, performance_err, ave_pathcounts, ave_pathcounts_err = net_performance(barbell, 100, 1)
    @test ave_pathcounts == [0, 1, 0, 0]
    # Only one path possible, so no error
    @test ave_pathcounts_err == [0, 0, 0, 0]

    # Test net_performance for TemporalGraph
    T = deepcopy(smalltemp)
    performance, performance_err, ave_pathcounts, ave_pathcounts_err = net_performance(T, 100, 2)
    # Check that the average number of paths used is 2
    @test ave_pathcounts[3] == 2.0
    @test ave_pathcounts_err == [0.0, 0.0, 0.0, 0.0]

    # Test make_user_pairs_wbuffer
    G = GridNetwork(5,5)
    pairs = QuNet.make_user_pairs_wbuffer(G, 9, buffer = 1)
    println("H")
    println(pairs)
end

max_stars_repo_name: FalafelGood/QuNet

max_issues_repo_name: FalafelGood/QuNet

max_forks_repo_name: FalafelGood/QuNet