text:
module DependencyWalker

using ObjectFile, Crayons
import Libdl

export Library

struct Library{OH<:Union{ObjectHandle,Missing,Nothing}}
    path::String
    handle_type::Type{OH}
    level::Int
    deps::Vector{Library}
    # Whether to show libraries whose metadata can't be read.  They're mostly
    # noisy.
    shownothing::Bool
end

function Library(path::String, level::Int = 0; shownothing::Bool=false)
    if !isfile(path)
        # TODO: should check also if it can be opened?
        return Library(path, Missing, level, Library[], shownothing)
    end
    io = open(path, "r")
    oh, deps_names = try
        oh = readmeta(io)
        # Get the list of needed libraries
        oh, keys(find_libraries(oh))
    catch
        nothing, []
    end
    close(io)
    deps = dependency_tree(deps_names, level; shownothing=shownothing)
    return Library(path, typeof(oh), level, deps, shownothing)
end

Library(path::String, oht::Type{T}, level::Int; shownothing::Bool=false) where {T<:Union{Missing,Nothing}} =
    Library{T}(path, oht, level, Library[], shownothing)

Library(path::AbstractString, oht, level; shownothing::Bool=false) =
    Library(String(path), oht, level; shownothing=shownothing)

# Check if `lib` matches `open_lib`, one of the libraries currently loaded in
# the system.  The condition to ignore the case is not quite accurate as this is
# a file-system property rather than an OS one, but in most cases this is
# correct.
@static if Sys.iswindows() || Sys.isapple()
    is_library_open(lib, open_lib) = occursin(lowercase(lib), lowercase(open_lib))
else
    is_library_open(lib, open_lib) = occursin(lib, open_lib)
end

function dependency_tree(deps_names, level::Int; dlext::String = Libdl.dlext, shownothing::Bool=false)
    # Initialise list of dependencies
    deps = Library[]
    # Get list of already dlopen'ed libraries
    open_dls = Libdl.dllist()
    for dep in deps_names
        if Sys.iswindows() && occursin(r"^api-ms-win-.*\.dll"i, dep)
            # Skipp all "api-ms-win-*" libraries on Windows
            continue
        end
        split_dep = split(dep, '.')
        # Get rid of the soversion.  TODO: this is only for GNU/Linux and
        # FreeBSD, no idea what we have to do for the other operating systems.
        idx = findlast(isequal(dlext), split_dep)
        if !isnothing(idx)
            dep = join(split_dep[1:idx], '.')
        end
        # Get the first dlopen'ed library matching the needed library, if any
        idx = findfirst(d -> is_library_open(dep, d), open_dls)
        if isnothing(idx)
            # Push a missing library to the list
            push!(deps, Library(dep, Missing, level + 1; shownothing=shownothing))
        else
            # Push the found library to the list
            push!(deps, Library(open_dls[idx], level + 1; shownothing=shownothing))
        end
    end
    return deps
end

reduce_hash(x::UInt64) = Base.hash_64_32(x)
reduce_hash(x::UInt32) = x

function Base.show(io::IO, lib::Library{T}; shownothing::Bool=false) where {T}
    if T === Missing
        symbol = "✗"
        extra_info = " (NOT FOUND)"
    elseif T === Nothing
        symbol = "❓"
        extra_info = " (COULD NOT READ METADATA)"
    else
        symbol = "◼"
        extra_info = ""
    end
    if !(T === Nothing && !shownothing)
        if lib.level > 0
            println(io)
        end
        print(io, repeat("  ", lib.level), Crayon(foreground = reduce_hash(hash(lib.path))), symbol,
              Crayon(reset=true), " ", lib.path, Crayon(foreground = :yellow), "$(extra_info)")
        for dep in lib.deps
            show(io, dep)
        end
    end
end

"""
DependencyWalker lets you walk through the dependencies of a shared library loaded in Julia.

The package exports a single function, `Library`, which takes as only argument the path to the shared library:

```
julia> using DependencyWalker, LibSSH2_jll

julia> LibSSH2_jll.libssh2_path # Path to the libssh2 library
"/home/user/.julia/artifacts/26c7d3a6c17151277018b133ab0034e93ddc3d1e/lib/libssh2.so"

julia> Library(LibSSH2_jll.libssh2_path)
◼ /home/user/.julia/artifacts/26c7d3a6c17151277018b133ab0034e93ddc3d1e/lib/libssh2.so
  ◼ /usr/bin/../lib/julia/libmbedtls.so.12
    ◼ /usr/bin/../lib/julia/libmbedx509.so.0
...
```
"""
DependencyWalker

end # module


text:
function FieldShortPeriodicsInterpolatedCoefficient(arg0::jint)
    return FieldShortPeriodicsInterpolatedCoefficient((jint,), arg0)
end

function add_grid_point(obj::FieldShortPeriodicsInterpolatedCoefficient, arg0::FieldAbsoluteDate, arg1::Vector{RealFieldElement})
    return jcall(obj, "addGridPoint", void, (FieldAbsoluteDate, Vector{RealFieldElement}), arg0, arg1)
end

function clear_history(obj::FieldShortPeriodicsInterpolatedCoefficient)
    return jcall(obj, "clearHistory", void, ())
end

function value(obj::FieldShortPeriodicsInterpolatedCoefficient, arg0::FieldAbsoluteDate)
    return jcall(obj, "value", Vector{RealFieldElement}, (FieldAbsoluteDate,), arg0)
end

ext: jl


text:
@testset "bayes.jl" begin

	x = dataset("datasets","faithful")  # RDatasets
	y = convert(Vector{Float64},x.Waiting)

	# test standard run
	g = GMM(2,y)
	grps = sortperm(g.w)
	@test all( isapprox.(g.w[grps], [0.36; 0.639],atol=0.01) )
	@test all( isapprox.(g.μ[grps], [54.61; 80.09],atol=0.01) )
	@test all( isapprox.(sqrt.(g.Σ[grps]), [5.867; 5.871],atol=0.01) )

	## only do k-means
	x = convert(Matrix,x)
	g = GMM(8, x, kind=:full, nIter=0)
	p = GMMprior(g.d, 1.0, 1.0)
	v = VGMM(g, p)
	em!(v, x)
	## show the reduction in Gaussians
	println(history(v))
	## we have no pretty printing yet...
	for f in [:α, :β, :m, :ν, :W]
	    println(f, " = ", getfield(v, f))
	end
end


text:
"""
Abstract type that describes Empirical Bayes estimands (which we want to estimate or conduct inference for).
"""
abstract type EBayesTarget end

Base.broadcastable(target::EBayesTarget) = Ref(target)

function (targets::AbstractArray{<:EBayesTarget})(prior)
    [target(prior) for target in targets]
end

abstract type AbstractPosteriorTarget <: EBayesTarget end
abstract type BasicPosteriorTarget <: AbstractPosteriorTarget end

"""
    LinearEBayesTarget <: EBayesTarget

Abstract type that describes Empirical Bayes estimands that are linear functionals of the
prior `G`.
"""
abstract type LinearEBayesTarget <: EBayesTarget end

Base.extrema(::EBayesTarget) = (-Inf, +Inf) # allow distribution-dependent choice?
_support(::LinearEBayesTarget) = Interval(nothing, nothing)

abstract type AbstractTargetComputation end

struct Conjugate <: AbstractTargetComputation end
struct NumeratorOfConjugate <: AbstractTargetComputation end
struct QuadgkQuadrature <: AbstractTargetComputation end

struct LinearOverLinear{N,D} <: AbstractTargetComputation
    num::N
    denom::D
end

Base.numerator(lin::LinearOverLinear) = lin.num
Base.denominator(lin::LinearOverLinear) = lin.denom

function (target::EBayesTarget)(prior)
    _loc = location(target) #TODO What about targets w/o location -> make it a trait!
    _comp = default_target_computation(target, _loc, prior)
    compute_target(_comp, target, _loc, prior)
end

function compute_target(comp::AbstractTargetComputation, target, prior)
    compute_target(comp, target, location(target), prior)
end

function default_target_computation(target::LinearEBayesTarget, sample, prior)
    QuadgkQuadrature()
end

function default_target_computation(target::AbstractPosteriorTarget, sample, prior)
    LinearOverLinear(nothing, nothing)
end

function compute_target(lin::LinearOverLinear, target::AbstractPosteriorTarget, sample, prior)
    # TODO: Actually take advantage of options in LinearOverLinear fields
    _num = numerator(target)(prior)
    _denom = denominator(target)(prior)
    _num/_denom
end

# TODO: Allow setting tolerances.
function compute_target(lin::QuadgkQuadrature, target::LinearEBayesTarget, sample, prior::ContinuousUnivariateDistribution)
   _interval = intersect(_support(target), _support(prior))
   if isbounded(_interval)
        _lower = first(_interval)
        _upper = last(_interval)
   elseif isa(_interval, Interval{T, Unbounded, Unbounded} where {T})
        _lower = -Inf
        _upper = +Inf
   elseif isa(_interval, Interval{T, Unbounded} where {T})
        _lower = -Inf
        _upper = last(_interval)
   else
        _lower = first(_interval)
        _upper = +Inf
   end
   quadgk( μ -> target(μ)*pdf(prior,μ), _lower, _upper)[1]
end

Base.@kwdef struct AffineTransformedLinearTarget{T, I, S} <: LinearEBayesTarget
    a::I = 0
    b::S = 1
    target::T
end

import Base.:*

function Base.:*(s::Number, t::LinearEBayesTarget)
    AffineTransformedLinearTarget(;b=s, target=t)
end

Base.:*(t::LinearEBayesTarget, s::Number) = s*t

(t::AffineTransformedLinearTarget)(prior::Union{<:Number, <:Distribution}) = t.b * t.target(prior) + t.a

"""
	cf(::LinearEBayesTarget, t)

The characteristic function of ``L(\\cdot)``, a `LinearEBayesTarget`, which we define as follows:

For ``L(\\cdot)`` which may be written as ``L(G) = \\int \\psi(\\mu)dG\\mu``
(for a measurable function ``\\psi``) this returns the Fourier transform of ``\\psi``
evaluated at t, i.e., ``\\psi^*(t) = \\int \\exp(it x)\\psi(x)dx``.
Note that ``\\psi^*(t)`` is such that for distributions ``G`` with density ``g``
(and ``g^*`` the Fourier Transform of ``g``) the following holds:
```math
L(G) = \\frac{1}{2\\pi}\\int g^*(\\mu)\\psi^*(\\mu) d\\mu
```
"""
function Distributions.cf(::LinearEBayesTarget, t) end

"""
	PriorDensity(z::Float64) <: LinearEBayesTarget
## Example call
```jldoctest
julia> PriorDensity(2.0)
PriorDensity{Float64}(2.0)
```
## Description
This is the evaluation functional of the density of ``G`` at `z`, i.e.,
``L(G) = G'(z) = g(z)`` or in Julia code `L(G) = pdf(G, z)`.
"""
struct PriorDensity{T} <: LinearEBayesTarget
    μ::T
end

location(target::PriorDensity) = target.μ

function Distributions.cf(target::PriorDensity{<:Real}, t)
    exp(im * location(target) * t)
end

#TODO: Not sure this is the right dispatch one(\mu) instead of e.g. one(Float64)
function (target::PriorDensity{<:Real})(μ::Number)
    location(target) == μ ? one(μ) : zero(μ)
end

function (target::PriorDensity{<:Interval})(μ::Number)
    in(μ, location(target)) ? one(μ) : zero(μ)
end

function (target::PriorDensity)(prior::Distribution)
    _pdf(prior, location(target))
end

Base.extrema(target::PriorDensity) = (0, Inf)
Base.extrema(target::PriorDensity{<:AbstractInterval}) = (0, 1)

"""
	MarginalDensity(Z::EBayesSample) <: LinearEBayesTarget
## Example call
```julia
MarginalDensity(StandardNormalSample(2.0))
```
## Description
Describes the marginal density evaluated at ``Z=z``  (e.g. ``Z=2`` in the example above).
In the example above the sample is drawn from the hierarchical model
```math
\\mu \\sim G, Z \\sim \\mathcal{N}(0,1)
```
In other words, letting ``\\varphi`` the Standard Normal pdf
```math
L(G) = \\varhi \\star dG(z)
```
Note that `2.0` has to be wrapped inside `StandardNormalSample(2.0)` since this target
depends not only on `G` and the location, but also on the likelihood.
"""
struct MarginalDensity{T} <: LinearEBayesTarget
    Z::T
end

location(target::MarginalDensity) = target.Z

function (target::MarginalDensity)(μ::Number)
    likelihood(location(target), μ)
end

function (target::MarginalDensity)(prior::Distribution)
    pdf(prior, location(target))
end

location(target::AbstractPosteriorTarget) = target.Z

"""
    Base.denominator(target::AbstractPosteriorTarget)

Suppose a posterior target ``\\theta_G(z)``, such as the posterior mean can be written as:
```math
\\theta_G(z) = \\frac{ a_G(z)}{f_G(z)} = \\frac{ \\int h(\\mu)dG(\\mu)}{\\int p(z \\mid \\mu)dG(\\mu)}.
```

For example, for the posterior mean ``h(\\mu) =  \\mu \\cdot p(z \\mid \\mu)``. Then `Base.denominator`
returns the linear functional representing ``G \\mapsto f_G(z)`` (i.e., typically the marginal density).
Also see [`Base.numerator(::AbstractPosteriorTarget)`](@ref).
"""
Base.denominator(target::AbstractPosteriorTarget) = MarginalDensity(location(target))

struct PosteriorTargetNumerator{T} <: LinearEBayesTarget
    posterior_target::T
end

location(target::PosteriorTargetNumerator) = location(target.posterior_target)

function default_target_computation(target::PosteriorTargetNumerator{<:BasicPosteriorTarget}, sample, prior)
    _post_target = default_target_computation(target.posterior_target, sample, prior)
    _post_target == Conjugate() ? NumeratorOfConjugate() : QuadgkQuadrature()
end

function compute_target(::NumeratorOfConjugate, post_numerator::PosteriorTargetNumerator, sample, prior)
    _post = post_numerator.posterior_target
    post_numerator.posterior_target(prior) * denominator(_post)(prior)
end

function (post_numerator::PosteriorTargetNumerator)(μ::Number)
    _post = post_numerator.posterior_target
    post_numerator.posterior_target(μ) * denominator(_post)(μ)
end

"""
    Base.numerator(target::AbstractPosteriorTarget)

Suppose a posterior target ``\\theta_G(z)``, such as the posterior mean can be written as:
```math
\\theta_G(z) = \\frac{ a_G(z)}{f_G(z)} = \\frac{ \\int h(\\mu)dG(\\mu)}{\\int p(z \\mid \\mu)dG(\\mu)}.
```

For example, for the posterior mean ``h(\\mu) =  \\mu \\cdot p(z \\mid \\mu)``. Then `Base.numerator`
returns the linear functional representing ``G \\mapsto a_G(z)``.
"""
Base.numerator(target::AbstractPosteriorTarget) = PosteriorTargetNumerator(target)

struct PosteriorTargetNullHypothesis{T,S} <: LinearEBayesTarget
    posterior_target::T
    c::S
end

location(target::PosteriorTargetNullHypothesis) = location(target.posterior_target)

function (post_null::PosteriorTargetNullHypothesis)(prior::T) where {T<:Union{Distribution,Number}}
    c = post_null.c
    _post = post_null.posterior_target
    numerator(_post)(prior) - c * denominator(_post)(prior)
end

"""
    PosteriorMean(Z::EBayesSample) <: AbstractPosteriorTarget

Type representing the posterior mean, i.e.,

```math
E_G[\\mu_i \\mid Z_i = z]
```
"""
struct PosteriorMean{T} <: BasicPosteriorTarget
    Z::T
end
PosteriorMean() = PosteriorMean(missing)

function compute_target(::Conjugate, postmean::PosteriorMean, Z::EBayesSample, prior)
    mean(posterior(Z, prior))
end

function (postmean::PosteriorMean)(μ::Number)
    μ
end

"""
    PosteriorVariance(Z::EBayesSample) <: AbstractPosteriorTarget

Type representing the posterior variance, i.e.,

```math
V_G[\\mu_i \\mid Z_i = z]
```

"""
struct PosteriorVariance{T} <: BasicPosteriorTarget
    Z::T
end
PosteriorVariance() = PosteriorVariance(missing)

function compute_target(::Conjugate, postmean::PosteriorVariance, Z::EBayesSample, prior)
    var(posterior(Z, prior))
end

Base.numerator(::PosteriorVariance) = throw("Posterior Variance is not fractional.")
Base.denominator(::PosteriorVariance) = throw("Posterior Variance is not fractional.")

"""
    PosteriorProbability(Z::EBayesSample, s) <: AbstractPosteriorTarget

Type representing the posterior probability, i.e.,

```math
\\Prob_G[\\mu_i \\in s \\mid Z_i = z]
```

"""
struct PosteriorProbability{T,S} <: BasicPosteriorTarget
    Z::T
    s::S
end

function compute_target(::Conjugate, postprob::PosteriorProbability, Z::EBayesSample, prior)
    _pdf(posterior(Z, prior), postprob.s)
end

function (postprob::PosteriorProbability{T, <:Interval})(μ::Number) where {T}
    _interval = postprob.s
    eltype(_interval)(μ in postprob.s)
end

function Base.extrema(target::PosteriorProbability)
    (0.0, 1.0)
end

# TODO: once we define support for posterior targets as well?
function _support(target::PosteriorTargetNumerator{<:PosteriorProbability})
    target.posterior_target.s
end

# Some additional linear targets that we keep unexported for now.

struct PriorMean <: LinearEBayesTarget end

(target::PriorMean)(μ::Number) = μ
(target::PriorMean)(prior::Distribution) = mean(prior)

struct PriorSecondMoment <: LinearEBayesTarget end

(target::PriorSecondMoment)(μ::Number) = abs2(μ)
(target::PriorSecondMoment)(prior::Distribution) = abs2(mean(prior)) + var(prior)

# Plotting code

@recipe function f(targets::AbstractVector{<:EBayesTarget}, g)
    length(unique(typeof.(targets))) == 1 || error("Expected homogeneous targets")
    xs = Float64.(location.(targets))
    ys = targets.(g)

    background_color_legend --> :transparent
    foreground_color_legend --> :transparent

    seriestype --> :path
    seriescolor --> "#550133"

    xs, ys
end


text:
@testset "softmax" begin
    X = rand(100,5)
    max_ = maximum(X, dims=2)
    exp_ = exp.(X .- max_)
    @test isapprox(MLJMultivariateStatsInterface.softmax(X), exp_ ./ sum(exp_, dims=2))
    X_ = MLJMultivariateStatsInterface.softmax!(X)
    @test X_ === X
    @test isapprox(X_, exp_ ./ sum(exp_, dims=2))
end


text:
# Return start of commn substring in s1, start of common substring in s2, and length of substring
# Indexes refer to character number, not index (differ for Unicode strings)
function longest_common_substring(s1::AbstractStringorGraphemeIterator, s2::AbstractStringorGraphemeIterator)
    if length(s1) > length(s2)
        start2, start1, size= longest_common_substring(s2, s1)
    else
        start1, start2, size = 0, 0, 0
        p = zeros(Int, length(s2))
        i1 = 0
        for ch1 in s1
            i1 += 1
            i2 = 0
            oldp = 0
            for ch2 in s2
                i2 += 1
                newp = 0
                if ch1 == ch2
                    newp = oldp > 0 ? oldp : i2
                    currentlength = (i2 - newp + 1)
                    if currentlength > size
                        start1, start2, size = i1 - currentlength + 1, newp, currentlength
                    end
                end
                p[i2], oldp = newp, p[i2]
            end
        end
    end
    return start1, start2, size
end

function matching_blocks!(x::Set{Tuple{Int, Int, Int}}, s1::AbstractStringorGraphemeIterator, s2::AbstractStringorGraphemeIterator, start1::Integer, start2::Integer)
    a = longest_common_substring(s1, s2)
    if a[3] > 0
        push!(x, (a[1] + start1 - 1, a[2] + start2 - 1, a[3]))
        s1before = SubString(s1, start(s1), chr2ind(s1, a[1]) - 1)
        s2before = SubString(s2, start(s2), chr2ind(s2, a[2]) - 1)
        matching_blocks!(x, s1before, s2before, start1, start2)
        if (a[1] + a[3]) <= endof(s1) && (a[2] + a[3]) <= endof(s2)
            s1after = SubString(s1, chr2ind(s1, a[1] + a[3]), endof(s1))
            s2after = SubString(s2, chr2ind(s2, a[2] + a[3]), endof(s2))
            matching_blocks!(x, s1after, s2after, start1 + a[1] + a[3] - 1, start2 + a[2] + a[3] - 1)
        end
    end
end

function matching_blocks(s1::AbstractStringorGraphemeIterator, s2::AbstractStringorGraphemeIterator)
    x = Set{Tuple{Int, Int, Int}}()
    matching_blocks!(x, s1, s2, 1, 1)
    return x
end

type RatcliffObershelp <: PreMetric end
function evaluate(dist::RatcliffObershelp, s1::AbstractStringorGraphemeIterator, s2::AbstractStringorGraphemeIterator, len1::Integer, len2::Integer)
    result = matching_blocks(s1, s2)
    matched = 0
    for x in result
        matched += x[3]
    end
    1.0 - 2 * matched / (len1 + len2)
end


text:
using RunTests
using Base.Test

@testmodule HouseHoldersTest begin
  using MLLibrary
  function test_QRDecompose()
    A=[2. 2 4. 18.;
       1. 3. -2. 1.;
       3. 1. 3. 14.]

    Q,R = MLLibrary.LinearAlgebra.HouseHolders.QRDecompose(A)
    @test round(Q,2) ==[-.53 -.22 .82;-.27 -.87 -.41; -.8 .44 -.41]
    @test round(R,2)==[-3.74 -2.67 -4.01 -21.11;0 -2.62 2.18 1.31; 0 0 2.86 8.57]
    @test round(norm(Q)) ==1.

    B=[1. -1. -1.;
       1. 2. 3.;
       2. 1. 1.;
       2. -2. 1.;
       3. 2. 1.]
     Q,R = MLLibrary.LinearAlgebra.HouseHolders.QRDecompose(B)
     @test round(Q,2) ==[-.23 -.23 -.46 -.46 -.69;.35 -.49 -.13 .71 -.34;-.37 .7 -.06 .47 -.38; -.67 -.42 .59 .11 -.1; -.49 -.19 -.65 .23 .5]'
     @test round(R,2) ==[-4.36 -1.15 -2.06; 0 -3.56 -1.58; 0 0 2.5; 0 0 0; 0 0 0]
     @test round(norm(Q)) ==1.  end

  function test_decompose()
    A=[2. 2 4. 18.;
       1. 3. -2. 1.;
       3. 1. 3. 14.]

    Q,R = MLLibrary.LinearAlgebra.HouseHolders.decompose(A)
    @test round(Q,2) ==[-.53 -.22 .82;-.27 -.87 -.41; -.8 .44 -.41]'
    @test round(R,2)==[-3.74 -2.67 -4.01 -21.11;0 -2.62 2.18 1.31; 0 0 2.86 8.57]
    @test round(norm(Q)) ==1.

    B=[1. -1. -1.;
       1. 2. 3.;
       2. 1. 1.;
       2. -2. 1.;
       3. 2. 1.]
     Q,R = MLLibrary.LinearAlgebra.HouseHolders.decompose(B)
     @test round(Q,2) ==[-.23 -.23 -.46 -.46 -.69;.35 -.49 -.13 .71 -.34;-.37 .7 -.06 .47 -.38; -.67 -.42 .59 .11 -.1; -.49 -.19 -.65 .23 .5]
     @test round(R,2) ==[-4.36 -1.15 -2.06; 0 -3.56 -1.58; 0 0 2.5; 0 0 0; 0 0 0]
     @test round(norm(Q)) ==1.

  end

end

max_stars_repo_name: ajermaky/CustomMLLibrary

max_issues_repo_name: ajermaky/CustomMLLibrary

max_forks_repo_name: ajermaky/CustomMLLibrary


text:
using Pkg

Pkg.add("ArgParse")
Pkg.add("DifferentialEquations")
Pkg.add("DiffEqCallbacks")
Pkg.add("Random")
Pkg.add("SparseArrays")
Pkg.add("StatsBase")
Pkg.add("Distributions")


text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: sagemaker_runtime
using AWS.Compat
using AWS.UUIDs

"""
    invoke_endpoint(body, endpoint_name)
    invoke_endpoint(body, endpoint_name, params::Dict{String,<:Any})

After you deploy a model into production using Amazon SageMaker hosting services, your
client applications use this API to get inferences from the model hosted at the specified
endpoint.  For an overview of Amazon SageMaker, see How It Works.  Amazon SageMaker strips
all POST headers except those supported by the API. Amazon SageMaker might add additional
headers. You should not rely on the behavior of headers outside those enumerated in the
request syntax.  Calls to InvokeEndpoint are authenticated by using AWS Signature Version
4. For information, see Authenticating Requests (AWS Signature Version 4) in the Amazon S3
API Reference. A customer's model containers must respond to requests within 60 seconds.
The model itself can have a maximum processing time of 60 seconds before responding to
invocations. If your model is going to take 50-60 seconds of processing time, the SDK
socket timeout should be set to be 70 seconds.  Endpoints are scoped to an individual
account, and are not public. The URL does not contain the account ID, but Amazon SageMaker
determines the account ID from the authentication token that is supplied by the caller.

# Arguments
- `body`: Provides input data, in the format specified in the ContentType request header.
  Amazon SageMaker passes all of the data in the body to the model.  For information about
  the format of the request body, see Common Data Formats-Inference.
- `endpoint_name`: The name of the endpoint that you specified when you created the
  endpoint using the CreateEndpoint API.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"Accept"`: The desired MIME type of the inference in the response.
- `"Content-Type"`: The MIME type of the input data in the request body.
- `"X-Amzn-SageMaker-Custom-Attributes"`: Provides additional information about a request
  for an inference submitted to a model hosted at an Amazon SageMaker endpoint. The
  information is an opaque value that is forwarded verbatim. You could use this value, for
  example, to provide an ID that you can use to track a request or to provide other metadata
  that a service endpoint was programmed to process. The value must consist of no more than
  1024 visible US-ASCII characters as specified in Section 3.3.6. Field Value Components of
  the Hypertext Transfer Protocol (HTTP/1.1).  The code in your model is responsible for
  setting or updating any custom attributes in the response. If your code does not set this
  value in the response, an empty value is returned. For example, if a custom attribute
  represents the trace ID, your model can prepend the custom attribute with Trace ID: in your
  post-processing function. This feature is currently supported in the AWS SDKs but not in
  the Amazon SageMaker Python SDK.
- `"X-Amzn-SageMaker-Inference-Id"`: If you provide a value, it is added to the captured
  data when you enable data capture on the endpoint. For information about data capture, see
  Capture Data.
- `"X-Amzn-SageMaker-Target-Container-Hostname"`: If the endpoint hosts multiple containers
  and is configured to use direct invocation, this parameter specifies the host name of the
  container to invoke.
- `"X-Amzn-SageMaker-Target-Model"`: The model to request for inference when invoking a
  multi-model endpoint.
- `"X-Amzn-SageMaker-Target-Variant"`: Specify the production variant to send the inference
  request to when invoking an endpoint that is running two or more variants. Note that this
  parameter overrides the default behavior for the endpoint, which is to distribute the
  invocation traffic based on the variant weights. For information about how to use variant
  targeting to perform a/b testing, see Test models in production
"""
invoke_endpoint(Body, EndpointName; aws_config::AbstractAWSConfig=global_aws_config()) = sagemaker_runtime("POST", "/endpoints/$(EndpointName)/invocations", Dict{String, Any}("Body"=>Body); aws_config=aws_config)
invoke_endpoint(Body, EndpointName, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = sagemaker_runtime("POST", "/endpoints/$(EndpointName)/invocations", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("Body"=>Body), params)); aws_config=aws_config)


text:
using Test
using SystemBenchmark

@testset "SystemBenchmark" begin
    res = runbenchmark();
    path, io = mktemp()
    savebenchmark(path, res)
    res2 = readbenchmark(path);
    comp = comparetoref(res);
    show(comp, allrows=true, allcols=true)

    crowd = getsubmittedbenchmarks()
    println("") #tidy test printing
end


text:
function cost_fn_sir(β, γ; k=1)
    return 1_000_000_000_000 * logistic(k*(10*base_reproduction_rate(β, γ) - 30))
end

const COST_FUNCTIONS = (
    cost_fn_sir = cost_fn_sir,
)


text:
using Stan

## testing accessors to CMDSTAN_HOME
let oldpath = Stan.CMDSTAN_HOME
    newpath = Stan.CMDSTAN_HOME * "##test##"
    set_cmdstan_home!(newpath)
    @test Stan.CMDSTAN_HOME == newpath
    set_cmdstan_home!(oldpath)
    @test Stan.CMDSTAN_HOME == oldpath
end


text:
println("Loading required packages...")
using Plots, JSON

function plot_results()
    println("Generating plots...")
    for filename in readdir()
        endswith(filename, ".json") || continue
        plot_results(replace(filename, ".json" => ""))
    end
end

function plot_results(filename)
    results = JSON.parsefile(filename*".json")
    if get(results, "type", nothing) == nothing
        @warn "Badly formatted file: `$filename.json'"
        return
    end

    println("-> $filename")
    plot_results(Val(Symbol(results["type"])), filename, results)
end

function plot_results(::Val{:accuracy}, filename, results)
    title = results["title"]
    labels = results["labels"]
    data = results["data"]

    scatter(title=title,
            xscale=:log10, yscale=:log10,
            xlabel="Condition number",
            ylabel="Relative error")

    markers = Symbol[:circle, :+, :rect, :x, :utriangle]

    for i in 1:length(labels)
        scatter!(Float64.(data[1]), Float64.(data[i+1]), label=labels[i], markershape=markers[i])
    end

    savefig(filename*".pdf")
    savefig(filename*".svg")
end

function plot_results(::Val{:ushift}, filename, results)
    title  = results["title"]
    labels = results["labels"]
    data   = results["data"]

    p = plot(title=title,
             xlabel="log2(U)",
             ylabel="Time [ns/elem]")

    for i in 1:length(labels)
        plot!(Float64.(data[1]), Float64.(data[i+1]), label="2^$(Int(round(log2(labels[i])))) elems")
    end

    savefig(filename*".pdf")
    savefig(filename*".svg")
end

function plot_results(::Val{:performance}, filename, results)
    title  = results["title"]
    labels = results["labels"]
    data   = results["data"]

    p = plot(title=title,
             xscale=:log10,
             xlabel="Vector size",
             ylabel="Time [ns/elem]")

    for i in 1:length(labels)
        plot!(Float64.(data[1]), Float64.(data[i+1]), label=labels[i])
    end

    savefig(filename*".pdf")
    savefig(filename*".svg")
end


text:
#######################################
#
# FTest
#
#######################################

function ftest(s::SSR; freq_of_interest::Union{Real, AbstractArray}=modulationrate(s), side_freq::Number=0.5,
    ID::AbstractString="", spill_bins::Int=2, results_key::AbstractString="statistics", kwargs... )

    # Do calculation here once, instead of in each low level call
    spectrum    = EEG._ftest_spectrum(s.processing["sweeps"])
    spectrum    = compensate_for_filter(s.processing, spectrum, samplingrate(s))
    frequencies = range(0, stop=1, length=Int(size(spectrum, 1)))*samplingrate(s)/2

    for freq in freq_of_interest

        snrDb, phase, signal, noise, statistic = ftest(spectrum, frequencies, freq, side_freq, spill_bins)

        result = DataFrame(ID                 = vec(repmat([ID], length(channelnames(s)), 1)),
                           Channel            = copy(channelnames(s)),
                           ModulationRate     = copy(modulationrate(s)),
                           AnalysisType       = vec(repmat(["F-test"], length(channelnames(s)))),
                           AnalysisFrequency  = vec(repmat([freq], length(channelnames(s)))),
                           SignalAmplitude    = vec(sqrt.(signal)),
                           SignalPhase        = vec(phase),
                           NoiseAmplitude     = vec(sqrt.(noise)),
                           SNRdB              = vec(snrDb),
                           Statistic          = vec(statistic))

        result = add_dataframe_static_rows(result, kwargs)

        if haskey(s.processing, results_key)
            s.processing[results_key] = vcat(s.processing[results_key], result)
        else
            s.processing[results_key] = result
        end

    end

    return s
end

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

# Save ftest results to file
function save_results(a::SSR; name_extension::AbstractString="", results_key::AbstractString="statistics", kwargs...)

    file_name = string(a.file_name, name_extension, ".csv")

    # Rename to save space
    results = a.processing

    # Index of keys to be exported
    result_idx = find_keys_containing(results, results_key)

    if length(result_idx) > 0

        to_save = get(results, collect(keys(results))[result_idx[1]], 0)

        CSV.write(file_name, to_save)
    end

    Logging.info("File saved to $file_name")

    return a
end


text:
export Cells, update_cells!

import CUDA
CUDA.allowscalar(false)

mutable struct Cells
    M::Int32
    cutoff::Float32
    action_cells::CUDA.CuArray{Int32,2}
    reaction_cells::CUDA.CuArray{Int32,2}
    head::CUDA.CuArray{Int32,1}
    next::CUDA.CuArray{Int32,1}
    index::CUDA.CuArray{Int32,1}
    population::CUDA.CuArray{Int32,1}
    collected::CUDA.CuArray{Int32,1}
    num_threads::Int32
    num_baskets::Int32
    basket_head::CUDA.CuArray{Int32,1}
    basket_count::CUDA.CuArray{Int32,1}
end

function index2voxel(index, M)
    k, l = divrem(index-1, M*M)
    j, i = divrem(l, M)
    return [i, j, k]
end

function stencil_vectors(rc, action)
    nmax = ceil(Int, rc)
    M = 1+2*nmax
    range = action ? (1:M^3÷2) : (M^3÷2+2:M^3)
    vectors = map(i->index2voxel(i, M) .- nmax, range)
    return filter((x->sum(x.^2) < rc^2) ∘ (x->abs.(x).-1), vectors)
end

cells_per_dimension(L, cutoff, ndiv) = floor(Int32, ndiv*L/cutoff)

function surrounding_cells(L, cutoff, M, action)
    pbc(x) = x < 0 ? x + M : (x >= M ? x - M : x)
    nearby_cell(index, vector) = 1 + (Array{Int32}([1, M, M^2]')*pbc.(index2voxel(index, M) + vector))[1]
    vectors = stencil_vectors(M*cutoff/L, action)
    nearby_cells = [Int32(nearby_cell(i, v)) for v in vectors, i in 1:M^3]
    return nearby_cells
end

function distribute!(head, next, index, population)
    icell = (CUDA.blockIdx().x - 1) * CUDA.blockDim().x + CUDA.threadIdx().x
    if icell <= length(head)
        head[icell] = 0
        population[icell] = 0
        for i = 1:length(index)
            if icell == index[i]
                next[i] = head[icell]
                head[icell] = i
                population[icell] += 1
            end
        end
    end
    return nothing
end

function clean_cells!(head, next, index, population, basket_head, basket_count, r, L, M)
    num_threads = CUDA.blockDim().x
    nbytes = sizeof(Int32)*num_threads
    thread_head = CUDA.@cuDynamicSharedMem(Int32, num_threads)
    thread_tail = CUDA.@cuDynamicSharedMem(Int32, num_threads, offset=nbytes)
    thread_count = CUDA.@cuDynamicSharedMem(Int32, num_threads, offset=2*nbytes)
    tid = CUDA.threadIdx().x
    bid = CUDA.blockIdx().x
    thread_head[tid] = 0
    thread_tail[tid] = 0
    thread_count[tid] = 0
    icell = (bid - 1)*num_threads + tid
    if icell <= length(head)
        if head[icell] != 0
            previous = head[icell]
            current = next[previous]
            while current != 0
                s1 = r[1, current]/L
                s2 = r[2, current]/L
                s3 = r[3, current]/L
                v1 = floor(Int32, M*(s1 - floor(s1)))
                v2 = floor(Int32, M*(s2 - floor(s2)))
                v3 = floor(Int32, M*(s3 - floor(s3)))
                index[current] = 1 + v1 + (v2 + v3*M)*M
                if index[current] == icell
                    previous = current
                    current = next[current]
                else
                    removed = current
                    current = next[current]
                    next[previous] = current
                    next[removed] = thread_head[tid]
                    (thread_head[tid] == 0) && (thread_tail[tid] = removed)
                    thread_head[tid] = removed
                    thread_count[tid] += 1
                    population[icell] -= 1
                end
            end
            current = head[icell]
            s1 = r[1, current]/L
            s2 = r[2, current]/L
            s3 = r[3, current]/L
            v1 = floor(Int32, M*(s1 - floor(s1)))
            v2 = floor(Int32, M*(s2 - floor(s2)))
            v3 = floor(Int32, M*(s3 - floor(s3)))
            index[current] = 1 + v1 + (v2 + v3*M)*M
            if index[current] != icell
                removed = current
                head[icell] = next[current]
                next[removed] = thread_head[tid]
                (thread_head[tid] == 0) && (thread_tail[tid] = removed)
                thread_head[tid] = removed
                thread_count[tid] += 1
                population[icell] -= 1
            end
        end
    end
    CUDA.sync_threads()

    s = 1
    while s < num_threads
        if tid % (2*s) == 1
            icell = tid + s
            if thread_head[icell] != 0
                next[thread_tail[icell]] = thread_head[tid]
                (thread_head[tid] == 0) && (thread_tail[tid] = thread_tail[icell])
                thread_head[tid] = thread_head[icell]
            end
            thread_count[tid] += thread_count[icell]
        end
        CUDA.sync_threads()
        s *= 2
    end

    if tid == 1
        basket_head[bid] = thread_head[1]
        basket_count[bid] = thread_count[1]
    end
    return nothing
end

function collect_baskets!(collected, basket_count, basket_head, next)
    ibasket = (CUDA.blockIdx().x - 1) * CUDA.blockDim().x + CUDA.threadIdx().x
    if ibasket <= length(basket_head)
        position = 0
        for i = 1:ibasket-1
            position += basket_count[i]
        end
        current = basket_head[ibasket]
        while current != 0
            position += 1
            collected[position] = current
            current = next[current]
        end
        basket_count[ibasket] = position
    end
    return nothing
end

function renew_cells!(head, next, population, collected, count, index)
    icell = (CUDA.blockIdx().x - 1) * CUDA.blockDim().x + CUDA.threadIdx().x
    if icell <= length(head)
        for position = 1:last(count)
            i = collected[position]
            if index[i] == icell
                next[i] = head[icell]
                head[icell] = i
                population[icell] += 1
            end
        end
    end
    return nothing
end

function Cells(r::CUDA.CuArray{T,2}, L, cutoff; ndiv=2, num_threads=256) where {T<:Number}
    M = cells_per_dimension(L, cutoff, ndiv)
    action_cells = surrounding_cells(L, cutoff, M, true)
    reaction_cells = surrounding_cells(L, cutoff, M, false)
    s = r'/L
    index = map(x->floor(Int32, x), M*(s .- floor.(s)))*CUDA.cu([1, M, M^2]) .+ 1
    num_cells = M^3
    num_particles = size(r, 2)
    num_baskets = ceil(Int, num_cells/num_threads)
    head = CUDA.zeros(Int32, num_cells)
    next = CUDA.zeros(Int32, num_particles)
    population = CUDA.zeros(Int32, num_cells)
    collected = CUDA.zeros(Int32, num_particles)
    basket_head = CUDA.zeros(Int32, num_baskets)
    basket_count = CUDA.zeros(Int32, num_baskets)
    CUDA.@cuda threads=num_threads blocks=num_baskets distribute!(head, next, index, population)
    return Cells(M, cutoff, action_cells, reaction_cells, head, next, index, population, collected,
                 num_threads, num_baskets, basket_head, basket_count)
end

function update_cells!(cells, r, L)
    CUDA.@cuda(
        threads=cells.num_threads,
        blocks=cells.num_baskets,
        shmem=3*cells.num_threads*sizeof(Int32),
        clean_cells!(
            cells.head, cells.next, cells.index, cells.population,
            cells.basket_head, cells.basket_count,
            r, L, cells.M
        )
    )
    CUDA.@cuda(
        threads=cells.num_threads,
        blocks=ceil(Int, cells.num_baskets/cells.num_threads),
        collect_baskets!(
            cells.collected, cells.basket_count, cells.basket_head, cells.next
        )
    )
    CUDA.@cuda(
        threads=cells.num_threads,
        blocks=cells.num_baskets,
        renew_cells!(
            cells.head, cells.next, cells.population, cells.collected, cells.basket_count, cells.index
        )
    )
    return nothing
end

function find_action_partners1!(action_partner, action_previous, action_count, pair_count,
                               head, next, index, action_cells,
                               rc, s, max_pairs_per_block, max_neighbors_per_atom)
    num_threads = CUDA.blockDim().x
    tid = CUDA.threadIdx().x
    bid = CUDA.blockIdx().x
    overall_previous = CUDA.@cuDynamicSharedMem(Int32, 1)
    tid == 1 && (overall_previous[1] = (bid - 1)*max_pairs_per_block)
    counter = CUDA.@cuDynamicSharedMem(Int32, num_threads, offset=sizeof(Int32))
    num_particles = length(index)
    i = (bid - 1)*num_threads + tid
    if i <= num_particles
        neighbor = CUDA.@cuDynamicSharedMem(Int32, max_neighbors_per_atom, offset=sizeof(Int32)+(num_threads+(tid-1)*max_neighbors_per_atom)*sizeof(Int32))
        s1i = s[1,i]
        s2i = s[2,i]
        s3i = s[3,i]
        rcsq = rc^2
        inside(dx, dy, dz) = (dx-round(dx))^2 + (dy-round(dy))^2 + (dz-round(dz))^2  <= rcsq
        n = 0
        icell = index[i]
        j = head[icell]
        while j != 0
            if j > i && inside(s[1,j]-s1i, s[2,j]-s2i, s[3,j]-s3i)
                n += 1
                if n <= max_neighbors_per_atom
                    neighbor[n] = j
                else
                    # distribute()
                end
            end
            j = next[j]
        end
        for k = 1:size(action_cells, 1)
            jcell = action_cells[k,icell]
            j = head[jcell]
            while j != 0
                if inside(s[1,j]-s1i, s[2,j]-s2i, s[3,j]-s3i)
                    n += 1
                    if n <= max_neighbors_per_atom
                        neighbor[n] = j
                    else
                        # distribute()
                    end
                end
                j = next[j]
            end
        end
        counter[tid] = min(n, max_neighbors_per_atom)
        action_count[i] = n
    end
    CUDA.sync_threads()

    if i <= num_particles
        previous_in_block = 0
        for t = 1:tid-1
            previous_in_block += counter[t]
        end
        previous = (bid - 1)*max_pairs_per_block + previous_in_block
        upper_limit = bid*max_pairs_per_block
        for n = 1:counter[tid]
            previous+n <= upper_limit && (action_partner[previous+n] = neighbor[n])
        end
        action_previous[i] = previous
        if tid == num_threads || i == num_particles
            pair_count[bid] = previous_in_block + counter[tid]
        end
    end
    CUDA.sync_threads()
    previous = (bid - 1)*max_pairs_per_block
    for n = pair_count[bid]+tid:num_threads:max_pairs_per_block
        action_partner[previous+n] = 0
    end
    return nothing
end

N = 2000
L = 10.0f0
cutoff = 2.0f0
r = L*CUDA.rand(Float32, 3, N)
cells = Cells(r, L, cutoff)

density = N/L^3
half_sphere_volume = 2π*cells.cutoff^3/3

max_neighbors_per_atom = ceil(Int, density*half_sphere_volume)  # Number of neighbors per particle

# max_neighbors_per_atom*(cells.num_threads + 1)*sizeof(Int32) = 49152
# max_max_neighbors_per_atom = 49152÷(cells.num_threads*sizeof(UInt8)) - 1
max_pairs_per_block = cells.num_threads
num_blocks = ceil(Int32, N/cells.num_threads)
action_partner = CUDA.zeros(Int32, max_pairs_per_block*num_blocks)
action_previous = CUDA.zeros(Int32, N)
action_count = CUDA.zeros(Int32, N)
action_loc = CUDA.zeros(Int32, N)
reaction_previous = CUDA.zeros(Int32, N)
reaction_count = CUDA.zeros(Int32, N)
pair_count = CUDA.zeros(Int32, num_blocks)

# max_neighbors_per_atom = (49152 - cells.num_threads*sizeof(Bool))÷(cells.num_threads*sizeof(Int32)) - 1
CUDA.@cuda(
    threads=cells.num_threads,
    blocks=num_blocks,
    shmem=(1 + max_neighbors_per_atom)*cells.num_threads*sizeof(Int32),
    find_action_partners1!(
        action_partner, action_previous, action_count, pair_count,
        cells.head, cells.next, cells.index, cells.action_cells,
        cutoff/L, r/L, max_pairs_per_block, max_neighbors_per_atom
    )
)
max_neighbors_per_atom


text:
using Cairo

minx(c) = min(c.p1[1],c.p2[1],c.p3[1])
maxx(c) = max(c.p1[1],c.p2[1],c.p3[1])
miny(c) = min(c.p1[2],c.p2[2],c.p3[2])
maxy(c) = max(c.p1[2],c.p2[2],c.p3[2])

function transformf(m,w,h)
    xs = map(v->v[1],m.vertices)
    ys = map(v->v[2],m.vertices)
    mminx, mmaxx = extrema(xs)
    mminy, mmaxy = extrema(ys)

    function transform(v)
        typeof(v)(tocr(v[1],w,mminx,mmaxx),tocr(v[2],h,mminy,mmaxy))
    end
end

export vertices

vertices(c::Cell2D) = [c.p1,c.p2,c.p3]

function _drawVertices(cr,vs,transform,colors)
    if colors != nothing
        @assert length(colors) == length(vs)
    end

    for i in 1:length(vs)
        v = vs[i]
        colors != nothing && set_source(cr,colors[i])
        vʹ = transform(v)
        circle(cr,vʹ[1],vʹ[2],2.0)
        fill(cr)
    end
end

function drawVertices(m::Meshes.Mesh{Vertex2},vs::Vector{Vertex2},w,h; colors = nothing, draw_mesh = false)
    c,cr = setupDrawing(w,h)

    transform = transformf(m,w,h)

    _drawVertices(cr,vs,transform,colors)

    set_source_rgb(cr,0.0,0.0,0.0);    # white

    draw_mesh && _drawMesh(cr,m,w,h,transform=transform)

    c
end

function _drawFaces(cr,faces,tvs; color = color("black"))
    set_source(cr,color);
    set_line_width (cr, 1.0);

    for f in faces
        s = tvs[f.v1][1:2]
        move_to(cr,s...)
        line_to(cr,tvs[f.v2][1:2]...)
        line_to(cr,tvs[f.v3][1:2]...)
        line_to(cr,s...)
    end
    stroke(cr);
end

function drawFaces(m,faces,w,h)
    transform = transformf(m,w,h)
    tvs = map(transform,m.vertices)
    c,cr = setupDrawing(w,h)
    _drawFaces(c,cr,faces,tvs)
    c
end

function setupDrawing(w,h)
    c = CairoRGBSurface(w,h);
    cr = CairoContext(c);

    save(cr);
    set_source_rgb(cr,1.0,1.0,1.0);    # white
    rectangle(cr,0.0,0.0,w,h); # background
    fill(cr);
    restore(cr);

    scale(cr,1.0,-1.0)
    translate(cr,0.0,-h)

    (c,cr)
end

function _drawMesh(cr::CairoContext, m::Meshes.Mesh{Vertex2}, w, h; transform = transformf(m,w,h))
    tvs = map(transform,m.vertices)

    set_source_rgb(cr,0.0,0.0,0.0);    # black
    set_line_width (cr, 1.0);
    for f in m.faces
        s = tvs[f.v1][1:2]
        move_to(cr,s...)
        line_to(cr,tvs[f.v2][1:2]...)
        line_to(cr,tvs[f.v3][1:2]...)
        line_to(cr,s...)
    end
    stroke(cr);
end
function drawMesh(m::Meshes.Mesh{Vertex2}, w, h)
    c,cr = setupDrawing(w,h)
    _drawMesh(cr,m,w,h)
    c
end

function drawNormals(m::Meshes.Mesh{Vertex2},edges::Vector{Edge},w,h; inward=false)
    c,cr = setupDrawing(w,h)

    transform = transformf(m,w,h)

    _drawMesh(cr,m,w,h,transform=transform)

    set_source_rgb(cr,1.0,0.0,0.0)
    for e in edges
        m = e.p1 + (e.p2-e.p1)./2.
        n = inward ? -n⃗(e)/20 : n⃗(e)/20
        endp = m+n
        endp = transform(Vertex2(endp[1],endp[2]))
        m = transform(m)
        move_to(cr,m[1],m[2])
        line_to(cr,endp[1],endp[2])
    end
    stroke(cr)

    c
end

using Winston

#colorbar
function colorbar(dmin, dmax; orientation="horizontal", colormap=Winston._current_colormap, kvs...)

    if orientation == "vertical"
        p=FramedPlot()
        setattr(p.x, draw_ticks=false)
        setattr(p.y1, draw_ticks=false)
        setattr(p.x1, draw_ticklabels=false)
        setattr(p.y1, draw_ticklabels=false)
        setattr(p.y2, draw_ticklabels=true)

        xr=(1,2)
        yr=(dmin,dmax)

        setattr(p, :xrange, xr)
        setattr(p.y1, range = yr)
        setattr(p.y2, range = yr)

        y=linspace(dmin, dmax, 256)*1.0
        data=[y y]
    elseif orientation == "horizontal"
        p=FramedPlot()
        setattr(p.y, draw_ticks=false)
        setattr(p.x1, draw_ticks=false)
        setattr(p.y1, draw_ticklabels=false)
        setattr(p.x1, draw_ticklabels=false)
        setattr(p.x2, draw_ticklabels=true)

        yr=(1,2)
        xr=(dmin, dmax)

        setattr(p, :yrange, yr)
        setattr(p.x1, range = xr)
        setattr(p.x2, range = xr)

        x=linspace(dmin,dmax,256)*1.0
        data=[x', x']
    end

    #setattr(p; kvs...)

    clims = (minimum(data),maximum(data))

    img = Winston.data2rgb(data, clims, colormap)
    add(p, Image(xr, yr, img))
    p
end

function default_colorize(data)
    colors = Winston.data2rgb(data, (minimum(data),maximum(data)), Winston._current_colormap)
    for x=1:size(data,1), y=1:size(data,2)
        if isnan(data[x,y])
            colors[x,y] = 0xD3D3D3
        end
    end

    c = colorbar(maximum(data),minimum(data), orientation = "vertical")

    colors, c
end

immutable BBox
    mminx::Float64
    mmaxx::Float64
    mminy::Float64
    mmaxy::Float64
end

tocr(x,s,mmin,mmax) = s*(x-mmin)/(mmax-mmin)
tom(x,s,mmin,mmax) = mmin+(x/s)*(mmax-mmin)

tocrx(x,s,bbox::BBox) = tocr(x,s,bbox.mminx,bbox.mmaxx)
tocry(x,s,bbox::BBox) = tocr(x,s,bbox.mminy,bbox.mmaxy)

tomx(x,s,bbox::BBox) = tom(x,s,bbox.mminx,bbox.mmaxx)
tomy(x,s,bbox::BBox) = tom(x,s,bbox.mminy,bbox.mmaxy)

tom(c::Cell2D,x,y,w,h,bbox::BBox) = Vertex2(tomx(x,w,bbox),tomy(y,h,bbox))
tocr(c::Cell2D,x,y,w,h,bbox::BBox)= Vertex2(tocrx(x,w,bbox),tocry(y,h,bbox))

function computebbox(m)
    xs = map(v->v[1],m.vertices)
    ys = map(v->v[2],m.vertices)
    mminx, mmaxx = extrema(xs)
    mminy, mmaxy = extrema(ys)
    BBox(mminx,mmaxx,mminy,mmaxy)
end

function pixelwise!(data,m::Meshes.Mesh{Vertex2},w,h,func; bbox=computebbox(m), naval=NaN)
    ns = computeNeighbors(m)
    for i in 1:length(m.faces)
        f = m.faces[i]
        # Find all the pixels whose center is in this triangle.
        # There's probably a smarter way to do this, but this is fine
        # for now.
        #@show i
        cell = cell2d(m,f,i)
        for x=floor(tocrx(minx(cell),w,bbox)):ceil(tocrx(maxx(cell),w,bbox)),
            y=floor(tocry(miny(cell),h,bbox)):ceil(tocry(maxy(cell),h,bbox))
            p = Vertex2(tomx(x,w,bbox),tomy(y,h,bbox))
            ff, dist = facingface(cell,p)
            if ff == 0
                data[clamp(x,1,w),clamp(y,1,h)] = func(cell,p)
            else
#                #@show dist
#                mindist = Inf
#                for face in 1:3
#                    has_neighbor(m,cell,face;ns=ns) || continue
#                    nn = neighbor(m,cell,face;ns=ns)
#                    _, ndist = facingface(nn,p)
#                    mindist = min(ndist,mindist)
#                    #@show mindist
#                    contains(nn,p) && assert(mindist == 0.)
#                    mindist == 0. && break
#                end
#                if mindist >= dist-w*h*eps() && data[clamp(x,1,w),clamp(y,1,h)] === naval
#                    data[clamp(x,1,w),clamp(y,1,h)] = func(cell,p)
#                end
            end
        end
    end
end

function plotMesh(m::Meshes.Mesh{Vertex2},w,h; func=nothing, colorize=identity, drawMesh=true)
    c,cr = setupDrawing(w+(func !== nothing ? 10+ceil(h/10) : 0),h)

    if func !== nothing
        #typeof(func(first(m),Vertex2(0,0)))
        data = Array(Float64,w,h)
        fill!(data,NaN)
        pixelwise!(data,m,w,h,func)
        colors, colorbar = colorize(data)
        ccr = Winston.CairoRenderer(c)
        Winston.compose(colorbar, ccr, Base.Graphics.BoundingBox(w+10,w+10+ceil(h/10),0,h))
        Cairo.image(cr,CairoImageSurface(colors,Cairo.FORMAT_RGB24,flipxy=false),0,0,w,h)
    end

    drawMesh && _drawMesh(cr,m,w,h)

    c
end

function drawArrow(cr, from, to; arrow_length = 15.0, arrow_degrees = 0.5)
    angle = atan2(to[2] - from[2], to[1] - from[1]) + pi;
    p1 = to + Vertex2(arrow_length * cos(angle - arrow_degrees),
                arrow_length * sin(angle - arrow_degrees))
    p2 = to + Vertex2(arrow_length * cos(angle + arrow_degrees),
                arrow_length * sin(angle + arrow_degrees))
    set_source_rgb(cr,1.0,0.0,0.0)
    move_to(cr,from[1],from[2])
    line_to(cr,to[1],to[2])
    line_to(cr,p1[1],p1[2])
    move_to(cr,to[1],to[2])
    line_to(cr,p2[1],p2[2])
    stroke(cr)
end

norm(v::Vertex2) = norm(v.coords)

function draw_vector_field(m,Q,w,h;nump=10)
    c,cr = setupDrawing(w,h)

    #_drawMesh(cr,m,w,h)

    bbox = computebbox(m)

    xps = map(round,linspace(0,w,nump))
    yps = map(round,linspace(0,h,nump))

    data = zeros(Int64,w,h)
    pixelwise!(data,m,w,h,(cell,p)->cell.cid; bbox=bbox)

    for x in xps, y in yps
        cid = data[clamp(x,1,w),clamp(y,1,h)]
        if cid==0
            continue
        end
        p = tom(m[cid],x,y,w,h,bbox)

        gf = ∇I(m[cid],Q[cid],p)

        #v = tocr(m[cid],gf[1],gf[2],w,h,bbox)
        v = (Ak(m[cid])')\gf
        #v = gf
        v ./= norm(v)
        v = 20.0*v
        #
        drawArrow(cr,Vertex2(x,y),Vertex2(v[1]+x,v[2]+y); arrow_length = 5)
    end

    c
end

##
# f = a f_1 + b f_2 + c f_3
# g = f(chi(p))
# g' = f'(chi(p))*chi'(p)
##

function drawref(w,h,m,cell,coeffs,func; colorize=identity, nump=10)
    c,cr = setupDrawing(w,h)
    data = Array(Float64,w,h)
    fill!(data,NaN)
    bbox = computebbox(m)
    for x=1:w, y=1:h
        point = cell.p1 + (x/w)*(cell.p2-cell.p1) + (y/h)*(cell.p3-cell.p1)
        if contains(cell,point)
            data[clamp(x,1,w),clamp(y,1,h)] = func(cell,point)
        end
    end

    Cairo.image(cr,CairoImageSurface(colorize(data)[1],Cairo.FORMAT_RGB24,flipxy=false),0,0,w,h)

    xps = map(round,linspace(0,w,nump))
    yps = map(round,linspace(0,h,nump))

    for x in xps, y in yps
        point = cell.p1 + (x/w)*(cell.p2-cell.p1) + (y/h)*(cell.p3-cell.p1)
        if contains(cell,point)
            gf = ∇I(cell,coeffs,point)

            #v = tocr(m[cid],gf[1],gf[2],w,h,bbox)
            #v = 𝜒⁻¹(m[cid],Vertex2(gf[1],gf[2]))-m[cid].p1
            v = gf
            v ./= norm(v)
            v = 20.0*v
            #
            drawArrow(cr,Vertex2(x,y),Vertex2(v[1]+x,v[2]+y); arrow_length = 5)
        end
    end

    c
end

drawSupport(pcg::CG2D,k::Int64,w=256,h=256;drawMesh=false) = drawSupport(pcg,[k],w,h;drawMesh=drawMesh)

function drawSupport(pcg::CG2D,k::Array{Int64},w=256,h=256;drawMesh=false)
    m = mesh(pcg)
    r = Face[]
    for c in m
        for i = 1:𝖓(pcg)
            if ℳ(pcg,c,porder(pcg),i) in k
                push!(r,m.faces[cid(c)])
            end
        end
    end
    c,cr = AC274.setupDrawing(w,h)
    transform = AC274.transformf(m,w,h)
    tvs = map(transform,m.vertices)
    drawMesh && AC274._drawMesh(cr,m,w,h,transform=transform)
    AC274._drawFaces(cr,r,tvs; color = drawMesh ? color("blue") : color("black"))
    v = [i > length(m.vertices) ? pcg.dualmesh.vertices[i-length(m.vertices)] : m.vertices[i] for i in k]
    AC274._drawVertices(cr,v,transform,[color("red") for _ in 1:length(v)])
    c
end

function plotSolution{N,T}(p::DG2D, Q::Array{Coeffs{N,T},2}, t)
    coeffs = Q[t,:]
    plotMesh(p.mesh,256,256; func=function (c,point)
        evaluate(coeffs[c.cid],AC274.getPhi(porder(p)),chi(c,point))
    end,
    colorize=AC274.default_colorize,
    drawMesh=false)
end

function plotSolution{T}(p::DG2D, QQ::Array{T,3}, t; tfunc = identity, w =256, h=256)
    coeffs = QQ[t,:,:]
    plotMesh(p.mesh,w,h; func=function (c,point)
        evaluate(tfunc(coeffs[1,c.cid,:]),AC274.getPhi(porder(p)),chi(c,point))
    end,
    colorize=AC274.default_colorize,
    drawMesh=false)
end

dof(m::Meshes.Mesh,c::AC274.Cell2D) = m.faces[cid(c)]
function ℳ(p::CG2D,c,porder,i)
    @assert i <= 3 || porder != 1
    m = mesh(p)
    dm = p.dualmesh
    f = dof(m,c)
    df = dof(dm,c)
    i == 1 ? f.v1 : i == 2 ? f.v2 : i == 3 ? f.v3 : (length(m.vertices) + (
    i == 4 ? df.v1 : i == 5 ? df.v2 : i == 6 ? df.v3 : error(BoundsError())))
end
ℳ(m::CG2D,c,i) = ℳ(m,c,1,i)

function AC274.plotSolution{T}(p::CG2D, QQ::Array{T,1}; tfunc = identity, w =256, h=256)
    plotMesh(p.mesh,w,h; func=function (c,point)
        coeffs = [QQ[ℳ(p,c,porder(p),i)] for i = 1:𝖓(p)]
        evaluate(tfunc(coeffs),AC274.getPhi(AC274.porder(p)),chi(c,point))
    end,
    colorize=AC274.default_colorize,
    drawMesh=false)
end

# Plotting Streamlines for (Navier-)Stokes
import PyPlot

function plotStreamLines(pu,Fsol)
    m = mesh(pu)
    bbox = AC274.computebbox(m)

    w,h = 256, 256

    data = zeros(Int64,w,h)
    AC274.pixelwise!(data,m,w,h,(cell,p)->cell.cid; bbox=bbox)
    data

    xps = linspace(bbox.mminx,bbox.mmaxx,w)
    yps = linspace(bbox.mminy,bbox.mmaxy,h)

    𝒩vel = 𝒩(pu)

    u = Array(Float64,w,h)
    v = Array(Float64,w,h)
    for i = 1:w, j=1:h
        x = xps[i]
        y = yps[j]
        cid = data[i,j]
        if cid == 0
            u[i,j] = 0
            v[i,j] = 0
            continue
        end
        c = m[cid]
        point = Vertex2(x,y)
        coeffsu = [Fsol[ℳ(pu,c,AC274.porder(pu),k)] for k = 1:𝖓(pu)]
        coeffsv = [Fsol[𝒩vel+ℳ(pu,c,AC274.porder(pu),k)] for k = 1:𝖓(pu)]
        u[i,j] = evaluate(coeffsu,AC274.getPhi(AC274.porder(pu)),chi(c,point))
        v[i,j] = evaluate(coeffsv,AC274.getPhi(AC274.porder(pu)),chi(c,point))
    end
    PyPlot.axes(aspect=1)
    PyPlot.streamplot(xps,yps,u',v',density=5)
end


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

using BAT
using Test
using Statistics, StatsBase, Distributions

@testset "Multimodal Tests" begin
    cm = MixtureModel([Distributions.Cauchy(), Distributions.Cauchy()])
    prod = Distributions.Product([Distributions.Cauchy(), Distributions.Cauchy()])

    #Instantiates multimodal Cauchy with Mixture Model and Product
    @test @inferred(BAT.MultimodalCauchy(cm, 0.1, 2, prod)) isa BAT.MultimodalCauchy
    #Instantiates multimodal Cauchy with μ, σ and n
    @test @inferred(BAT.MultimodalCauchy(μ = 1., σ = 0.1, n=4)) isa BAT.MultimodalCauchy
    #Error if n = 1
    @test_throws ArgumentError  BAT.MultimodalCauchy(μ = 1., σ = 0.1, n=1)

    mmc = BAT.MultimodalCauchy(μ = 1., σ = 0.1, n=4)

    @test @inferred(Base.size(mmc)) == (4,)
    @test @inferred(Base.length(mmc)) == 4
    @test @inferred(Base.eltype(mmc)) == Float64

    params = @inferred(StatsBase.params(mmc))

    #Parameters
    @test params[1] == [-1., 1., 0., 0.]
    @test params[2] == [0.1, 0.1, 0.1, 0.1]
    @test params[3] == 4

    #logpdf
    @test @inferred(Distributions._logpdf(mmc, [-1., 0., 1., 2.])) == -11.283438151658

    #If μ = 0 the Distribution should be Cauchy-like in every dimension
    mmc = BAT.MultimodalCauchy(μ = 0., σ = 0.1, n=2)
    ks_test = HypothesisTests.ExactOneSampleKSTest(rand(mmc, 10^6)[1,:], Cauchy(0., 0.1))#First dimension
    @test pvalue(ks_test) > 0.05
    ks_test = HypothesisTests.ExactOneSampleKSTest(rand(mmc, 10^6)[2,:], Cauchy(0., 0.1))#Second dimension
    @test pvalue(ks_test) > 0.05
end


text:
#####
##### Profiling batched inference
#####

using AlphaZero
using Setfield
import CUDA

using Plots

# using Revise; Revise.includet("scripts/profile/inference.jl")
# profile_inference(on_gpu=true, batch_size=128, num_filters=128)
# With 128 filters: 27μs per sample
# With  64 filters: 17μs per sample

"""
Estimate inference time in μs/sample.
"""
function profile_inference(
    exp::Experiment = Examples.experiments["connect-four"];
    on_gpu=true,
    nrep=100,
    batch_size=128,
    num_filters=64)

  params = exp.netparams
  params = @set params.num_filters = num_filters
  net = exp.mknet(exp.gspec, params)
  net = Network.copy(net; on_gpu, test_mode=true)
  state = GI.current_state(GI.init(exp.gspec))
  batch = [state for _ in 1:batch_size]
  Network.evaluate_batch(net, batch) # To compile everything
  GC.gc(true)
  info = CUDA.@timed for _ in 1:nrep Network.evaluate_batch(net, batch) end
  return info.time / nrep / batch_size * 1_000_000
end

"""
Plot the inference time speedup (in μs/sample) as a function of batch-size.
"""
function plot_inference_speedup(
  file::String,
  exp::Experiment = Examples.experiments["connect-four"];
  on_gpu=true,
  num_filters=64)

  N = 8
  batch_sizes = [2^i for i in 0:N]
  ts = [
    profile_inference(exp;
      nrep=10,
      num_filters, batch_size, on_gpu)
    for batch_size in batch_sizes ]

  xticks = (0:N, map(string, batch_sizes))
  plot(0:N, ts[1] ./ ts,
    title="Inference Speedup (baseline: $(Int(ceil(ts[1]))) μs)",
    ylabel="Speedup",
    xlabel="Batch size",
    ylims=(0, Inf),
    legend=:none,
    xticks=xticks)
  savefig(file)
end

function all_plots()
  plot_inference_speedup("inference-gpu-64.png",  on_gpu=true,  num_filters=64)
  plot_inference_speedup("inference-cpu-64.png",  on_gpu=false, num_filters=64)
  plot_inference_speedup("inference-gpu-128.png", on_gpu=true,  num_filters=128)
  plot_inference_speedup("inference-cpu-128.png", on_gpu=false, num_filters=128)
end


text:
"""
    add_crowding_constraints(weights::Dict{Int64, Float64})

Adds a molecular crowding constraint to the optimization problem: `∑ wᵢ × vᵢ ≤ 1` where `wᵢ`
is a weight and `vᵢ` is a flux index in the model's reactions specified in `weights` as `vᵢ
=> wᵢ` pairs.

See Beg, Qasim K., et al. "Intracellular crowding defines the mode and sequence of
substrate uptake by Escherichia coli and constrains its metabolic activity." Proceedings of
the National Academy of Sciences 104.31 (2007) for more details.
"""
add_crowding_constraints(weights::Dict{Int64,Float64}) =
    (model, opt_model) -> begin
        idxs = collect(keys(weights)) # order of keys and values is the same
        ws = values(weights)
        # since fluxes can be positive or negative, need absolute value: ∑ wᵢ × |vᵢ| ≤ 1
        # introduce slack variables to handle this
        @variable(opt_model, crowding_slack[1:length(weights)])
        @constraint(opt_model, crowding_slack .>= opt_model[:x][idxs])
        @constraint(opt_model, crowding_slack .>= -opt_model[:x][idxs])
        @constraint(opt_model, sum(w * crowding_slack[i] for (i, w) in enumerate(ws)) <= 1)
    end

"""
    add_crowding_constraints(weights::Dict{String, Float64})

Variant of [`add_crowding_constraints`](@ref) that takes a dictinary of reactions `ids`
instead of reaction indices mapped to weights.
"""
add_crowding_constraints(weights::Dict{String,Float64}) =
    (model, opt_model) -> begin
        idxs = indexin(keys(weights), reactions(model))
        nothing in idxs && throw(ArgumentError("Reaction id not found in model."))
        add_crowding_constraints(Dict(zip(Int.(idxs), values(weights))))(model, opt_model)
    end


text:
function store_int!(X::Array{Int64,1}, buf::Array{UInt8,1}, i::Int64, n::Int64)
  n == 0 && return
  nx = length(X)
  if n > 0
    (n > nx) && append!(X, zeros(Int64, n-nx))
    X[n] = buf_to_int(buf, i)
  end
  return nothing
end

function store_dbl!(X::Array{Float64,1}, buf::Array{UInt8,1}, i::Int64, n::Int64)
  n == 0 && return
  nx = length(X)
  if n > 0
    (n > nx) && append!(X, zeros(Float64, n-nx))
    X[n] = buf_to_double(buf, i)
  end
  return nothing
end

function get_coeff_n(io::IO, c::UInt8, buf::Array{UInt8,1})
  i = 1
  while c != 0x20
    buf[i] = c
    i += 1
    c = fastread(io)
  end
  return buf_to_int(buf, i-1)+1
end

function to_newline(io::IO, c::UInt8)
  while c != 0x0a
    c = fastread(io)
  end
  return c
end

function skip_whitespace(io::IO, c::UInt8)
  c = fastread(io)
  while c == 0x20
    c = fastread(io)
  end
  return c
end

function parse_resp_date(io::IO, T::Array{UInt16,1})
  fill!(T, zero(UInt16))
  i = mark(io)
  k = 1
  while true
    c = fastread(io)
    if c in (0x2c, 0x2e, 0x3a)
      L = fastpos(io)-i-1
      reset(io)
      T[k] = stream_int(io, L)
      fastskip(io, 1)
      i = mark(io)
      k += 1
    elseif c == 0x7e
      L = fastpos(io)-i-1
      if L > 0
        reset(io)
        T[k] = stream_int(io, L)
        fastskip(io, 1)
      end
      break
    end
  end
  return mktime(T)
end

function string_field(sio::IO)
  m = mark(sio)
  c = fastread(sio)
  while c != 0x7e
    c = fastread(sio)
  end
  p = fastpos(sio)-m
  reset(sio)
  return String(fastread(sio, p))[1:p-1]
end

function skip_string!(sio::IO)
  c = fastread(sio)
  while c != 0x7e
    c = fastread(sio)
  end
  return nothing
end

function blk_string_read(io::IO, nb::Int64, v::Integer)
  checkbuf!(BUF.buf, nb)
  if nb + BUF.k > 4096
    part1 = 4096-BUF.k
    fast_readbytes!(io, BUF.buf, part1)
    fastread!(io, BUF.seq)
    nb -= part1
    while nb > 4088
      bv = pointer(BUF.buf, part1+1)
      fast_unsafe_read(io, bv, 4088)
      fastread!(io, BUF.seq)
      part1 += 4088
      nb -= 4088
    end
    bv = pointer(BUF.buf, part1+1)
    fast_unsafe_read(io, bv, nb)
    BUF.k = 8+nb
  else
    fast_readbytes!(io, BUF.buf, nb)
    BUF.k += nb
  end
  sio = IOBuffer(BUF.buf)
  return sio
end

seed_time(u16::Array{UInt16, 1}, hh::UInt8, mm::UInt8, ss::UInt8, δt::Int64) =
  y2μs(u16[1]) + Int64(u16[2] - one(UInt16))*86400000000 + Int64(u16[3])*100 +
  Int64(hh)*3600000000 + Int64(mm)*60000000 + Int64(ss)*1000000 + δt


text:
ENV["PLOTS_TEST"] = "true"
ENV["GKSwstype"] = "100"
using TrackRoots
import TrackRoots:Mark
using Base.Test

using DataDeps

ENV["DATADEPS_ALWAY_ACCEPT"]=true
register(DataDep("test",
                "These are two folders with their `nd` files and multiple dark and light timelapse 16 bit TIF images (1.9 GB total).",
                "https://s3.eu-central-1.amazonaws.com/yakirgagnon/roots/test.zip",
                "d99b8f2edce5e72104ba1cfed02798a6683bc22f4b76a5cb823c80257bc4bf48",
                post_fetch_method=unpack))

# # incase there are old test files, remove them just in case they are old/bad
rm(joinpath(datadep"test", "1", "d2"), recursive=true, force=true)

@testset "all" begin
    startpoints = [[Mark[[774.3011654713115, 911.7642161885246]]], [Mark[[761.334080430328, 413.93003970286884]], Mark[[811.9541095671107, 510.3906089907787]]]]
    endpoints = [[Mark[[905.2760077185163, 918.9115390162102]]], [Mark[[1016.9281490727263, 209.57046484226038]], Mark[[1017.350592982095, 315.11227808371984]]]]
    for (i, dataset) in enumerate(["1","5"])
        files = readdir(joinpath(datadep"test", dataset))
        ind = findfirst(x -> last(splitext(x)) == ".nd", files)
        ndfile = joinpath(datadep"test", dataset, files[ind])
        stages = TrackRoots.nd2stages(ndfile)
        Δx = TrackRoots.pixel_width(stages[1])
        calibstages = TrackRoots.stages2calib(stages, Δx)
        tracks = TrackRoots.trackroot.(calibstages, startpoints[i])
        endpoint = map(tracks) do t
            map(t) do r
                r.coordinates[end]
            end
        end
        @test endpoints[i] == endpoint
    end
end

@testset "utils" begin
    a = TrackRoots.disk(1)
    @test all(CartesianIndex(i) ∈ a for i in [(-1,0), (0,-1), (0,0), (0,1), (1,0)])
    @test TrackRoots.isnd("name.nd")
    @test !TrackRoots.isnd("name.ndd")
    ndfiles = TrackRoots.findall_nd(datadep"test")
    @test length(ndfiles) == 2
end

@testset "startpoint" begin

    @test TrackRoots.get_point(GtkReactive.XY{GtkReactive.UserUnit}(.2,.4)) == Mark(.4,.2)

    ndfile = joinpath(datadep"test", "1", "d2.nd")
    stages = TrackRoots.nd2stages(ndfile)
    file = stages[1].timelapse[1].dark
    img, imgtmp, ms = TrackRoots.fetch_image(file)
    @test img == imgtmp
    @test ms == [0.0038452735179674985, 0.005874723430228122]

    M = 0.5
    message = "test"
    s, b = TrackRoots.build_slider(M, message)
    @test GtkReactive.value(s) == M
    @test b isa Gtk.GtkBoxLeaf

end

@testset "plotting" begin

    startpoint = [Mark[[774.3011654713115, 911.7642161885246]]]
    ndfile = joinpath(datadep"test", "1", "d2.nd")
    main(ndfile, startpoint, Base.DevNullStream())

    @test isfile(joinpath(datadep"test", "1", "d2", "stage 1", "roots.png"))
    @test isfile(joinpath(datadep"test", "1", "d2", "stage 1", "root 1", "coordinates.csv"))
    @test isfile(joinpath(datadep"test", "1", "d2", "stage 1", "root 1", "intensities.csv"))
    @test isfile(joinpath(datadep"test", "1", "d2", "stage 1", "root 1", "summary.mp4"))

    #clean up
    #=rm(joinpath(first(DataDeps.default_loadpath), "test", "1", "d2"), recursive=true, force=true)=#
end