text:
const SimpleGraphEdge = SimpleEdge

"""
    SimpleGraph{T}

A type representing an undirected graph.
"""
mutable struct SimpleGraph{T <: Integer} <: AbstractSimpleGraph{T}
    ne::Int
    fadjlist::Vector{Vector{T}} # [src]: (dst, dst, dst)
end

eltype(x::SimpleGraph{T}) where T = T

# Graph{UInt8}(6), Graph{Int16}(7), Graph{UInt8}()
function SimpleGraph{T}(n::Integer=0) where T <: Integer
    fadjlist = [Vector{T}() for _ = one(T):n]
    return SimpleGraph{T}(0, fadjlist)
end

# SimpleGraph()
SimpleGraph() = SimpleGraph{Int}()

# SimpleGraph of a SimpleGraph
SimpleGraph(g::SimpleGraph) = copy(g)

# SimpleGraph(6), SimpleGraph(0x5)
SimpleGraph(n::T) where T <: Integer = SimpleGraph{T}(n)

# SimpleGraph(UInt8)
SimpleGraph(::Type{T}) where T <: Integer = SimpleGraph{T}(zero(T))

# Graph{UInt8}(adjmx)
function SimpleGraph{T}(adjmx::AbstractMatrix) where T <: Integer
    dima, dimb = size(adjmx)
    isequal(dima, dimb) || throw(ArgumentError("Adjacency / distance matrices must be square"))
    issymmetric(adjmx) || throw(ArgumentError("Adjacency / distance matrices must be symmetric"))

    g = SimpleGraph(T(dima))
    @inbounds for i in findall(triu(adjmx) .!= 0)
        add_edge!(g, i[1], i[2])
    end
    return g
end

# converts Graph{Int} to Graph{Int32}
function SimpleGraph{T}(g::SimpleGraph) where T <: Integer
    h_fadj = [Vector{T}(x) for x in fadj(g)]
    return SimpleGraph(ne(g), h_fadj)
end

# SimpleGraph(adjmx)
SimpleGraph(adjmx::AbstractMatrix) = SimpleGraph{Int}(adjmx)

# SimpleGraph(digraph)
function SimpleGraph(g::SimpleDiGraph)
    gnv = nv(g)
    edgect = 0
    newfadj = deepcopy(g.fadjlist)
    @inbounds for i in vertices(g)
        for j in badj(g, i)
            index = searchsortedfirst(newfadj[i], j)
            if index <= length(newfadj[i]) && newfadj[i][index] == j
                edgect += 1     # this is an existing edge - we already have it
                if i == j
                    edgect += 1 # need to count self loops
                end
            else
                insert!(newfadj[i], index, j)
                edgect += 2      # this is a new edge only in badjlist
            end
        end
    end
    iseven(edgect) || throw(AssertionError("invalid edgect in graph creation - please file bug report"))
    return SimpleGraph(edgect ÷ 2, newfadj)
end

@inbounds function cleanupedges!(fadjlist::Vector{Vector{T}}) where T <: Integer
    neg = 0
    for v in 1:length(fadjlist)
        if !issorted(fadjlist[v])
            sort!(fadjlist[v])
        end
        unique!(fadjlist[v])
        neg += length(fadjlist[v])
        # self-loops should count as one edge
        for w in fadjlist[v]
            if w == v
                neg += 1
                break
            end
        end
    end
    return neg ÷ 2
end

function SimpleGraph(edge_list::Vector{SimpleGraphEdge{T}}) where T <: Integer
    nvg = zero(T)
    @inbounds(
    for e in edge_list
        nvg = max(nvg, src(e), dst(e))
    end)

    list_sizes = ones(Int, nvg)
    degs = zeros(Int, nvg)
    @inbounds(
    for e in edge_list
        s, d = src(e), dst(e)
        (s >= 1 && d >= 1) || continue
        degs[s] += 1
        if s != d
            degs[d] += 1
        end
    end)

    fadjlist = Vector{Vector{T}}(undef, nvg)
    @inbounds(
    for v in 1:nvg
        fadjlist[v] = Vector{T}(undef, degs[v])
    end)

    @inbounds(
    for e in edge_list
        s, d = src(e), dst(e)
        (s >= 1 && d >= 1) || continue
        fadjlist[s][list_sizes[s]] = d
        list_sizes[s] += 1
        if s != d
            fadjlist[d][list_sizes[d]] = s
            list_sizes[d] += 1
        end
    end)

    neg = cleanupedges!(fadjlist)
    g = SimpleGraph{T}()
    g.fadjlist = fadjlist
    g.ne = neg

    return g
end

@inbounds function add_to_fadjlist!(fadjlist::Vector{Vector{T}}, s::T, d::T) where T <: Integer
    nvg = length(fadjlist)
    nvg_new = max(nvg, s, d)
    for v = (nvg + 1):nvg_new
        push!(fadjlist, Vector{T}())
    end

    push!(fadjlist[s], d)
    if s != d
        push!(fadjlist[d], s)
    end
end

function _SimpleGraphFromIterator(iter)::SimpleGraph
    T = Union{}
    fadjlist = Vector{Vector{T}}()
    @inbounds(
    for e in iter
        typeof(e) <: SimpleGraphEdge ||
                        throw(ArgumentError("iter must be an iterator over SimpleEdge"))
        s, d = src(e), dst(e)
        (s >= 1 && d >= 1) || continue
        if T != eltype(e)
            T = typejoin(T, eltype(e))
            fadjlist = convert(Vector{Vector{T}}, fadjlist)
        end
        add_to_fadjlist!(fadjlist, s, d)
    end)

    T == Union{} && return SimpleGraph(0)
    neg  = cleanupedges!(fadjlist)
    g = SimpleGraph{T}()
    g.fadjlist = fadjlist
    g.ne = neg

    return g
end

function _SimpleGraphFromIterator(iter, ::Type{SimpleGraphEdge{T}}) where T <: Integer
    fadjlist = Vector{Vector{T}}()
    @inbounds(
    for e in iter
        s, d = src(e), dst(e)
        (s >= 1 && d >= 1) || continue
        add_to_fadjlist!(fadjlist, s, d)
    end)

    neg  = cleanupedges!(fadjlist)
    g = SimpleGraph{T}()
    g.fadjlist = fadjlist
    g.ne = neg

    return g
end

"""
    SimpleGraphFromIterator(iter)

Creates a SimpleGraph from an iterator iter. The elements in iter must
be of type <: SimpleEdge.
"""
function SimpleGraphFromIterator(iter)::SimpleGraph
    if Base.IteratorEltype(iter) == Base.EltypeUnknown()
        return _SimpleGraphFromIterator(iter)
    end
    if eltype(iter) <: SimpleGraphEdge && isconcretetype(eltype(iter))
        return _SimpleGraphFromIterator(iter, eltype(iter))
    end
    return _SimpleGraphFromIterator(iter)
end

edgetype(::SimpleGraph{T}) where T <: Integer = SimpleGraphEdge{T}

"""
    badj(g::SimpleGraph[, v::Integer])

Return the backwards adjacency list of a graph. If `v` is specified,
return only the adjacency list for that vertex.

###Implementation Notes
Returns a reference, not a copy. Do not modify result.
"""
badj(g::SimpleGraph) = fadj(g)
badj(g::SimpleGraph, v::Integer) = fadj(g, v)

"""
    adj(g[, v])

Return the adjacency list of a graph. If `v` is specified, return only the
adjacency list for that vertex.

### Implementation Notes
Returns a reference, not a copy. Do not modify result.
"""
adj(g::SimpleGraph) = fadj(g)
adj(g::SimpleGraph, v::Integer) = fadj(g, v)

copy(g::SimpleGraph) =  SimpleGraph(g.ne, deepcopy(g.fadjlist))

==(g::SimpleGraph, h::SimpleGraph) =
vertices(g) == vertices(h) &&
ne(g) == ne(h) &&
fadj(g) == fadj(h)

"""
    is_directed(g)

Return `true` if `g` is a directed graph.
"""
is_directed(::Type{SimpleGraph}) = false
is_directed(::Type{SimpleGraph{T}}) where T = false
is_directed(g::SimpleGraph) = false

function has_edge(g::SimpleGraph{T}, s, d) where T
    verts = vertices(g)
    (s in verts && d in verts) || return false  # edge out of bounds
    @inbounds list_s = g.fadjlist[s]
    @inbounds list_d = g.fadjlist[d]
    if length(list_s) > length(list_d)
        d = s
        list_s = list_d
    end
    return insorted(d, list_s)
end

function has_edge(g::SimpleGraph{T}, e::SimpleGraphEdge{T}) where T
    s, d = T.(Tuple(e))
    return has_edge(g, s, d)
end

"""
    add_edge!(g, e)

Add an edge `e` to graph `g`. Return `true` if edge was added successfully,
otherwise return `false`.
"""
function add_edge!(g::SimpleGraph{T}, e::SimpleGraphEdge{T}) where T
    s, d = T.(Tuple(e))
    verts = vertices(g)
    (s in verts && d in verts) || return false  # edge out of bounds
    @inbounds list = g.fadjlist[s]
    index = searchsortedfirst(list, d)
    @inbounds (index <= length(list) && list[index] == d) && return false  # edge already in graph
    insert!(list, index, d)

    g.ne += 1
    s == d && return true  # selfloop

    @inbounds list = g.fadjlist[d]
    index = searchsortedfirst(list, s)
    insert!(list, index, s)
    return true  # edge successfully added
end

"""
    rem_edge!(g, e)

Remove an edge `e` from graph `g`. Return `true` if edge was removed successfully,
otherwise return `false`.

### Implementation Notes
If `rem_edge!` returns `false`, the graph may be in an indeterminate state, as
there are multiple points where the function can exit with `false`.
"""
function rem_edge!(g::SimpleGraph{T}, e::SimpleGraphEdge{T}) where T
    s, d = T.(Tuple(e))
    verts = vertices(g)
    (s in verts && d in verts) || return false  # edge out of bounds
    @inbounds list = g.fadjlist[s]
    index = searchsortedfirst(list, d)
    @inbounds (index <= length(list) && list[index] == d) || return false  # edge not in graph
    deleteat!(list, index)

    g.ne -= 1
    s == d && return true  # selfloop

    @inbounds list = g.fadjlist[d]
    index = searchsortedfirst(list, s)
    deleteat!(list, index)
    return true  # edge successfully removed
end

"""
    add_vertex!(g)

Add a new vertex to the graph `g`. Return `true` if addition was successful.
"""
function add_vertex!(g::SimpleGraph{T}) where T
    (nv(g) + one(T) <= nv(g)) && return false       # test for overflow
    push!(g.fadjlist, Vector{T}())
    return true
end


text:
using BinDeps
@BinDeps.setup

dl_dir = joinpath(dirname(dirname(@__FILE__)), "deps", "downloads")
deps_dir = joinpath(dirname(dirname(@__FILE__)), "deps")
usr_dir = joinpath(deps_dir, "usr")
lib_dir = joinpath(deps_dir, "usr", "lib")
src_dir = joinpath(deps_dir, "src")

# Cleanup
isdir(dl_dir) && rm(dl_dir, recursive=true)
isdir(usr_dir) && rm(usr_dir, recursive=true)
isdir(src_dir) && rm(src_dir, recursive=true)

bit = (Int==Int32) ? "32" : "64"

# Verions of libraries
pathlib_v = "4.7.03"
pathjulia_v = "0.1.2"

# The main dependency
libpath47julia = library_dependency("libpath47julia")

# File paths for linux
libpath47_so = joinpath(src_dir, "pathlib-$pathlib_v", "lib", "linux$bit", "libpath47.so")
path47julia_c = joinpath(src_dir, "PathJulia-$pathjulia_v", "src", "pathjulia.c")

# File paths for windows
libpath47_dll = joinpath(src_dir, "pathlib-$pathlib_v", "lib", "win$bit", "path47.dll")
libpath47_lib = joinpath(src_dir, "pathlib-$pathlib_v", "lib", "win$bit", "path47.lib")
libpath47julia_dll = joinpath(src_dir, "PathJulia-$pathjulia_v", "lib", "win$bit", "libpath47julia.dll")

provides(Sources, URI("https://github.com/chkwon/PathJulia/archive/$pathjulia_v.tar.gz"),
  libpath47julia, unpacked_dir="PathJulia-$pathjulia_v", os = :Linux)

# Mac OS X
provides(BuildProcess,
    (@build_steps begin
        GetSources(libpath47julia)
        CreateDirectory(lib_dir, true)
        @build_steps begin
            ChangeDirectory(joinpath(src_dir, "PathJulia-$pathjulia_v", "src"))
            # cp("../lib/osx/libpath47julia.dylib", lib_dir)
            `cp -f ../lib/osx/libpath47julia.dylib $lib_dir`
        end
    end), libpath47julia, os = :Darwin)
# provides(Binaries, URI("https://github.com/chkwon/PathJulia/archive/$pathjulia_v.zip"),
#          libpath47julia, unpacked_dir="PathJulia-$pathjulia_v/lib/osx", os = :Darwin)

# Linux 32/64
provides(BuildProcess,
    (@build_steps begin
        GetSources(libpath47julia)
        CreateDirectory(lib_dir, true)
        @build_steps begin
            ChangeDirectory(joinpath(src_dir, "PathJulia-$pathjulia_v", "src"))
            `make linux$bit`
            `cp -f ../lib/linux$bit/libpath47julia.so $lib_dir`
        end
    end), libpath47julia, os = :Linux)

# Windows 32/64
provides(BuildProcess,
    (@build_steps begin
        CreateDirectory(lib_dir, true)
        CreateDirectory(src_dir, true)
        @build_steps begin
            ChangeDirectory(src_dir)
            FileDownloader("https://github.com/ampl/pathlib/archive/$pathlib_v.zip", joinpath(dl_dir, "pathlib.zip"))
            FileUnpacker(joinpath(dl_dir, "pathlib.zip"), src_dir, libpath47_dll)
            `powershell -NoProfile -Command "Copy-Item -Path '$libpath47_dll' -Destination '$lib_dir' -force"`
            `powershell -NoProfile -Command "Copy-Item -Path '$libpath47_lib' -Destination '$lib_dir' -force"`
        end
        @build_steps begin
            FileDownloader("https://github.com/chkwon/PathJulia/archive/$pathjulia_v.zip", joinpath(dl_dir, "PathJulia.zip"))
            FileUnpacker(joinpath(dl_dir, "PathJulia.zip"), src_dir, libpath47julia_dll)
            `powershell -NoProfile -Command "Copy-Item -Path '$libpath47julia_dll' -Destination '$lib_dir' -force"`
        end
    end), libpath47julia, os = :Windows)

@static if Sys.iswindows()
    push!(BinDeps.defaults, BuildProcess)
end

@BinDeps.install Dict(:libpath47julia => :libpath47julia)

@static if Sys.iswindows()
    pop!(BinDeps.defaults)
end

ext: jl


text:
@inline function predictive_stepsize_controller!(integrator, alg)
  # Gustafsson predictive stepsize controller
  if iszero(integrator.EEst)
    q = inv(integrator.opts.qmax)
  else
    gamma = integrator.opts.gamma
    if typeof(alg) <: Union{RKC,IRKC}
      fac = gamma
    else
      if alg isa RadauIIA5
        @unpack nl_iters = integrator.cache
        @unpack max_iter = alg
      else
        @unpack nl_iters, max_iter = integrator.cache.nlsolver
      end
      fac = min(gamma,(1+2*max_iter)*gamma/(nl_iters+2*max_iter))
    end
    expo = 1/(get_current_adaptive_order(integrator.alg,integrator.cache)+1)
    qtmp = (integrator.EEst^expo)/fac
    @fastmath q = max(inv(integrator.opts.qmax),min(inv(integrator.opts.qmin),qtmp))
    integrator.qold = q
  end
  q
end

@inline function standard_stepsize_controller!(integrator, alg)
  # Standard stepsize controller
  if iszero(integrator.EEst)
    q = inv(integrator.opts.qmax)
  else
    qtmp = integrator.EEst^(1/(get_current_adaptive_order(integrator.alg,integrator.cache)+1))/integrator.opts.gamma
    @fastmath q = max(inv(integrator.opts.qmax),min(inv(integrator.opts.qmin),qtmp))
    integrator.qold = integrator.dt/q
  end
  q
end

@inline function PI_stepsize_controller!(integrator, alg)
  # PI-controller
  EEst,beta1,q11,qold,beta2 = integrator.EEst, integrator.opts.beta1, integrator.q11,integrator.qold,integrator.opts.beta2
  if iszero(EEst)
    q = inv(integrator.opts.qmax)
  else
    @fastmath q11 = EEst^beta1
    @fastmath q = q11/(qold^beta2)
    integrator.q11 = q11
    @fastmath q = max(inv(integrator.opts.qmax),min(inv(integrator.opts.qmin),q/integrator.opts.gamma))
  end
  q
end

function stepsize_controller!(integrator,alg)
  if ispredictive(alg)
    predictive_stepsize_controller!(integrator, alg)
  elseif isstandard(alg)
    standard_stepsize_controller!(integrator, alg)
  else # Default is PI-controller
    PI_stepsize_controller!(integrator, alg)
  end
end

@inline function predictive_step_accept_controller!(integrator, alg, q)
  if q <= integrator.opts.qsteady_max && q >= integrator.opts.qsteady_min
    q = one(q)
  end
  if integrator.success_iter > 0
    expo = 1/(get_current_adaptive_order(integrator.alg,integrator.cache)+1)
    qgus=(integrator.dtacc/integrator.dt)*(((integrator.EEst^2)/integrator.erracc)^expo)
    qgus = max(inv(integrator.opts.qmax),min(inv(integrator.opts.qmin),qgus/integrator.opts.gamma))
    qacc=max(q,qgus)
  else
    qacc = q
  end
  integrator.dtacc = integrator.dt
  integrator.erracc = max(1e-2,integrator.EEst)
  return integrator.dt/qacc
end

standard_step_accept_controller!(integrator, alg, q) = integrator.dt/q

@inline function PI_step_accept_controller!(integrator, alg, q)
  if q <= integrator.opts.qsteady_max && q >= integrator.opts.qsteady_min
    q = one(q)
  end
  integrator.qold = max(integrator.EEst,integrator.opts.qoldinit)
  return integrator.dt/q # dtnew
end

function step_accept_controller!(integrator, alg, q)
  if ispredictive(alg)
    predictive_step_accept_controller!(integrator, alg, q)
  elseif isstandard(alg)
    standard_step_accept_controller!(integrator, alg, q)
  else
    PI_step_accept_controller!(integrator, alg, q)
  end
end

function step_reject_controller!(integrator,alg)
  if ispredictive(alg)
    integrator.dt = integrator.success_iter == 0 ? 0.1integrator.dt : integrator.dt/integrator.qold
  elseif isstandard(alg)
    integrator.dt = integrator.qold
  else #PI
    integrator.dt = integrator.dt/min(inv(integrator.opts.qmin),integrator.q11/integrator.opts.gamma)
  end
  return
end

function stepsize_controller!(integrator, alg::JVODE)
  if iszero(integrator.EEst)
    η = integrator.opts.qmax
  else
    η = integrator.cache.η
    integrator.qold = η
  end
  η
end

function step_accept_controller!(integrator,alg::JVODE,η)
  return η * integrator.dt  # dtnew
end

function step_reject_controller!(integrator,alg::JVODE)
  integrator.dt *= integrator.qold
end

function stepsize_controller!(integrator, alg::QNDF)
  cnt = integrator.iter
  if cnt <= 3
    q = standard_stepsize_controller!(integrator, alg)
  else
    q = integrator.dt/integrator.cache.h
    integrator.qold = integrator.dt/q
  end
  q
end


text:
using BinDeps

@BinDeps.setup

nettle = library_dependency("nettle", aliases = ["libnettle"], runtime = false)
gnutls = library_dependency("gnutls", aliases = ["libgnutls.so.28","libgnutls","libgnutls28"], depends = [nettle], validate = function(p,h)
	if !haskey(ENV,"GNUTLS_VERSION")
		return true
	end
	ccall(dlsym(h,:gnutls_check_version),Ptr{Uint8},(Ptr{Uint8},),ENV["GNUTLS_VERSION"]) != C_NULL
end)

provides(Sources,{
	URI("http://www.lysator.liu.se/~nisse/archive/nettle-2.7.tar.gz") => nettle,
	URI("ftp://ftp.gnutls.org/gcrypt/gnutls/v3.2/gnutls-3.2.1.tar.xz") => gnutls})

provides(Binaries,URI("ftp://ftp.gnutls.org/gcrypt/gnutls/w32/gnutls-3.2.1-w32.zip"),gnutls,os = :Windows)

if haskey(ENV,"GNUTLS_VERSION")
	requested_version = convert(VersionNumber,ENV["GNUTLS_VERSION"])
	pkgmanager_validate	= (p,dep)->available_version(p) >= requested_version
else
	pkgmanager_validate = true
end
## Pending successfull bottling, we're just going to fall back on Source for now
#provides(Homebrew,"gnutls",gnutls,validate = pkgmanager_validate)
provides(AptGet,"libgnutls28",gnutls,validate = pkgmanager_validate) # Yes, this is the most current version, I guess they broke binary compatibility in v2.8?
provides(Yum,"libgnutls",gnutls,validate = pkgmanager_validate)

julia_usrdir = normpath(JULIA_HOME*"/../") # This is a stopgap, we need a better builtin solution to get the included libraries
libdirs = String["$(julia_usrdir)/lib"]
includedirs = String["$(julia_usrdir)/include"]

env = {"HOGWEED_LIBS" => "-L$(libdirs[1]) -L$(BinDeps.libdir(nettle)) -lhogweed -lgmp",
		"NETTLE_LIBS" => "-L$(libdirs[1]) -L$(BinDeps.libdir(nettle)) -lnettle -lgmp", "LIBS" => "-lgmp ","LD_LIBRARY_PATH"=>join([libdirs[1];BinDeps.libdir(nettle)],":")}

provides(BuildProcess,Autotools(lib_dirs = libdirs, include_dirs = includedirs, env = env),nettle)

# If we're installing gnutls from source we better also installl nettle from source, otherwise we end up with a giant mess.
provides(BuildProcess,Autotools(libtarget = "lib/libgnutls.la", lib_dirs = libdirs, include_dirs = includedirs, env = env),gnutls,force_depends = {BuildProcess => nettle})

@BinDeps.install


text:
let
    global gammln, gammq

    FPMIN::Float64 = 2.22507e-308 / eps()
    ngau::Int = 18

    y::Vector{Float64} = [0.0021695375159141994,
    0.011413521097787704,0.027972308950302116,0.051727015600492421,
    0.082502225484340941, 0.12007019910960293,0.16415283300752470,
    0.21442376986779355, 0.27051082840644336, 0.33199876341447887,
    0.39843234186401943, 0.46931971407375483, 0.54413605556657973,
    0.62232745288031077, 0.70331500465597174, 0.78649910768313447,
    0.87126389619061517, 0.95698180152629142]

    w::Vector{Float64} = [0.0055657196642445571,
    0.012915947284065419,0.020181515297735382,0.027298621498568734,
    0.034213810770299537,0.040875750923643261,0.047235083490265582,
    0.053244713977759692,0.058860144245324798,0.064039797355015485,
    0.068745323835736408,0.072941885005653087,0.076598410645870640,
    0.079687828912071670,0.082187266704339706,0.084078218979661945,
    0.085346685739338721,0.085983275670394821]

    cof::Vector{Float64} = [57.1562356658629235,-59.5979603554754912,
        14.1360979747417471,-0.491913816097620199,.339946499848118887e-4,
        .465236289270485756e-4,-.983744753048795646e-4,.158088703224912494e-3,
        -.210264441724104883e-3,.217439618115212643e-3,-.164318106536763890e-3,
        .844182239838527433e-4,-.261908384015814087e-4,.368991826595316234e-5]

    function gammpapprox(a::Float64, x::Float64, psig::Int)
        gln::Float64    = gammln(a)
        a1::Float64     = a - 1.0
        lna1::Float64   = log(a1)
        sqrta1::Float64 = sqrt(a1)

        xu::Float64 = 0.0
        t::Float64 = 0.0
        sum__::Float64 = 0.0
        ans::Float64 = 0.0

        if (x > a1)
            xu = max(a1 + 11.5 * sqrta1, x + 6.0 * sqrta1)
        else
            xu = max(0.0, min(a1 - 7.5 * sqrta1, x - 5.0 * sqrta1))
        end

        for j in 1:ngau
            t = x + (xu - x) * y[j]
            sum__ += w[j] * exp(-(t - a1) + a1*(log(t) - lna1))
        end

        ans = sum*(xu - x)*exp(a1*(lna1 - 1.0) - gln)

        return (psig ? (ans > 0.0 ? 1.0 - ans : -ans)
                    : (ans >= 0.0 ? ans : 1.0 + ans))
    end

    function gser(a::Float64, x::Float64)
        gln::Float64 = gammln(a)

        ap::Float64  = a
        sum__::Float64 = 1.0 / a
        del::Float64 = sum__

        while true
            ap += 1
            del *= x/ap
            sum__ += del
            if (abs(del) < abs(sum__) * eps())
                return sum__ * exp(-x + a * log(x) - gln)
            end
        end
    end

    function gcf(a::Float64, x::Float64)
        gln::Float64 = gammln(a)

        b::Float64 = x + 1.0 - a
        c::Float64 = 1.0/FPMIN
        d::Float64 = 1.0/b
        h::Float64 = d

        i::Int = 1
        while i < typemax(Int)
            an::Float64 = -i * (i - a)

            b += 2.0
            d = an * d + b
            if (abs(d) < FPMIN)
                d = FPMIN
            end

            c = b + an/c

            if (abs(c) < FPMIN)
                c = FPMIN
            end

            d = 1.0/d
            del::Float64 = d*c
            h *= del

            if (abs(del - 1.0) <= eps())
                break
            end
            i += 1
        end

        return h * exp( -x + a * log(x) - gln)
    end

    function gammq(a::Float64, x::Float64)
        """
        Implements the q-gamma function which seemingly isn't available
        in SpecialFunctions.jl.
        """
        ASWITCH::Int = 100

        @assert (x >= 0.0 && a > 0.0) "Negative numbers are not valid for gamma function!"

        if (x == 0.0)
            return 1.0
        elseif (Int(round(a)) >= ASWITCH)
            return gammpapprox(a, x, 0)
        elseif (x < a + 1.0)
            return 1.0 - gser(a,x)
        else
            return gcf(a,x)
        end
    end

    function gammln(xx::Float64)
        """
        Log gamma function
        """
        @assert xx > 0.0 "Input must be greater than 0.0."

        x::Float64 = xx
        z::Float64 = xx

        tmp::Float64 = x + 5.24218750000000000
        tmp = (x + 0.5) * log(tmp) - tmp

        ser = 0.999999999999997092
        for j in 1:length(cof)
            z += 1
            ser += cof[j] / z
        end

        return tmp + log(2.5066282746310005 * ser / x)
    end

end

max_stars_repo_name: heindelj/MolecularPotentials

max_issues_repo_name: heindelj/MolecularPotentials

max_forks_repo_name: heindelj/MolecularPotentials


text:
module TestUtils
using Test
using BenchmarkTools
import Statistics: mean
import Markdown: Table, MD

export btest

function init()
    BenchmarkTools.DEFAULT_PARAMETERS.seconds = 1
    global table = [["Title", "Mem [MiB]", "Allocs", "Time [ms]"]]

    if "AOC_PERF" in keys(ENV)
        println(join(table[1], "\t"))
    end
end

function btest(fun, title)
    @test fun()
    if "AOC_PERF" in keys(ENV)
        r = @benchmark $fun()
        row = string.([title,
                       round(r.memory/(1024^2); digits=2),
                       r.allocs,
                       round(mean(r.times)/1000000; digits=2)])
        push!(table, row)
        println(join(row, "\t"))
    end
end

function output(fname)
    if "AOC_PERF" in keys(ENV)
        open(fname, "w") do f
            write(f, string(MD(Table(table, [:l,:r,:r,:r]))))
        end
    end
end

end #module


text:
using BlochSim, Test, MAT

function testA5b()

    answer = matread("matlabtestdata/testA5b.mat")

    T1 = 600 # ms
    T2 = 100 # ms
    Δf = 10 # Hz
    M = [1.0, 0, 0]
    dt = 1 # ms
    time = 0:dt:1000 # ms
    spin = Spin(M, 1, T1, T2, Δf, [0,0,0])
    mag = zeros(3,length(time))
    mag[:,1] = M
    (A, B) = freeprecess(spin, dt)
    for t = 2:length(time)
      applydynamics!(spin, A, B)
      mag[:,t] = spin.M
    end

    return mag ≈ answer["M"]

end

function testF1a()

    answer = matread("matlabtestdata/testF1a.mat")

    t = 2.3 # ms
    α = π/4 # rad
    θ = 0 # rad
    Δf = -500:5:500 # Hz
    T1 = 600 # ms
    T2 = 100 # ms
    spins = map(Δf -> Spin([0,0,1.0], 1, T1, T2, Δf, [0,0,0]), Δf)
    for f = 1:length(Δf)
      (Ae, _) = excitation(spins[f], θ, α)
      (Af, Bf) = freeprecess(spins[f], t)
      applydynamics!(spins[f], Ae)
      applydynamics!(spins[f], Af, Bf)
      applydynamics!(spins[f], Ae)
    end
    sig = map(spin -> spin.signal, spins)

    return sig ≈ vec(answer["sig"])

end

function testF1b()

    answer = matread("matlabtestdata/testF1b.mat")

    t = 2.3 # ms
    α = π/4 # rad
    θ = 0 # rad
    Δf = 0 # Hz
    T1 = 600 # ms
    T2 = 100 # ms
    grad = [0.1, 0, 0] # G/cm
    xpos = -2:0.01:2 # cm
    pos = [[x, 0, 0] for x in xpos]
    spins = map(pos -> Spin([0,0,1.0], 1, T1, T2, Δf, pos), pos)
    for i = 1:length(xpos)
      (Ae, _) = excitation(spins[i], θ, α)
      (Af, Bf) = freeprecess(spins[i], t, grad)
      applydynamics!(spins[i], Ae)
      applydynamics!(spins[i], Af, Bf)
      applydynamics!(spins[i], Ae)
    end
    sig = map(spin -> spin.signal, spins)

    return sig ≈ vec(answer["sig"])

end

function testF1c()

    answer = matread("matlabtestdata/testF1c.mat")

    t = 2.3 # ms
    α = π/4 # rad
    θ = 0 # rad
    Δf = 0 # Hz
    T1 = 600 # ms
    T2 = 100 # ms
    grad = [0.1, 0, 0] # G/cm
    xpos = -2:0.01:2 # cm
    pos = [[x, 0, 0] for x in xpos]
    spins = map(pos -> Spin([0,0,1.0], 1, T1, T2, Δf, pos), pos)
    for i = 1:length(xpos)
      (Ae, _) = excitation(spins[i], θ, α)
      (Af, Bf) = freeprecess(spins[i], t, grad)
      (Ar, Br) = freeprecess(spins[i], t, -grad/2)
      applydynamics!(spins[i], Ae)
      applydynamics!(spins[i], Af, Bf)
      applydynamics!(spins[i], Ae)
      applydynamics!(spins[i], Ar, Br)
    end
    sig = map(spin -> spin.signal, spins)

    return sig ≈ vec(answer["sig"])

end

function testF2c()

    answer = matread("matlabtestdata/testF2c.mat")

    dt = 0.04 # ms
    t = 0:dt:6 # ms
    rf = 0.05 * sinc.(t .- 3) # G
    Δf = -1000:20:1000 # Hz
    spins = map(Δf -> Spin([0,0,1.0], 1, 600, 100, Δf, [0,0,0]), Δf)
    map(spin -> excitation!(spin, rf, 0, zeros(3), dt), spins)
    sig = map(spin -> spin.signal, spins)

    return sig ≈ vec(answer["sig"])

end

function testF3a()

    answer = matread("matlabtestdata/testF3a.mat")

    dt = 0.1 # ms
    t = 0:dt:6 # ms
    T = length(t)
    xpos = -2:0.01:2 # cm
    rf = 0.05 * sinc.(t .- 3) # G
    grad = [0.1, 0, 0]
    T1 = 600 # ms
    T2 = 100 # ms
    Δf = 0 # Hz
    pos = [[x, 0, 0] for x in xpos]
    spins = map(pos -> Spin([0,0,1.0], 1, T1, T2, Δf, pos), pos)
    for i = 1:length(xpos)
      (Ae, Be) = excitation(spins[i], rf, 0, grad, dt)
      applydynamics!(spins[i], Ae, Be)
    end
    sig = map(spin -> spin.signal, spins)

    return sig ≈ vec(answer["sig"])

end

function testF3c()

    answer = matread("matlabtestdata/testF3c.mat")

    dt = 0.05 # ms
    t = 0.1:dt:6 # ms
    T = length(t)
    xpos = -2:0.1:2 # cm
    rf = 0.05 * sinc.(t .- 3) # G
    grad = [0.1ones(T)'; zeros(T)'; zeros(T)']
    T1 = 600 # ms
    T2 = 100 # ms
    Δf = 0 # Hz
    pos = [[x, 0, 0] for x in xpos]
    spins = map(pos -> Spin([0,0,1.0], 1, T1, T2, Δf, pos), pos)
    for i = 1:length(xpos)
        (Ae, Be) = excitation(spins[i], [rf; 0rf], 0, [grad -0.52grad], dt)
        applydynamics!(spins[i], Ae, Be)
    end
    sig = map(spin -> spin.signal, spins)
    Mz = map(spin -> spin.M[3], spins)

    return sig ≈ vec(answer["sig"]) && Mz ≈ vec(answer["mz"])

end

function testF3d()

    answer = matread("matlabtestdata/testF3d.mat")

    dt = 0.05 # ms
    t = 0.1:dt:6 # ms
    T = length(t)
    xpos = -2:0.1:2 # cm
    rf = 0.05 * sinc.(t .- 3) # G
    grad = [0.1ones(T)'; zeros(T)'; zeros(T)']
    T1 = 600 # ms
    T2 = 100 # ms
    Δf = 100 # Hz
    pos = [[x, 0, 0] for x in xpos]
    spins = map(pos -> Spin([0,0,1.0], 1, T1, T2, Δf, pos), pos)
    map(spin -> excitation!(spin, [rf; 0rf], 0, [grad -0.52grad], dt), spins)
    sig = map(spin -> spin.signal, spins)
    Mz = map(spin -> spin.M[3], spins)

    return sig ≈ vec(answer["sig"]) && Mz ≈ vec(answer["mz"])

end

function testF3f()

    answer = matread("matlabtestdata/testF3f.mat")

    dt = 0.05 # ms
    t = 0.1:dt:6 # ms
    T = length(t)
    xpos = -5:0.1:5 # cm
    rf = 0.05 * sinc.(t .- 3) .* (exp.(im * 2π * 900 * t/1000) + exp.(-im * 2π * 900 * t/1000)) # G
    grad = [0.1ones(T)'; zeros(T)'; zeros(T)']
    T1 = 600 # ms
    T2 = 100 # ms
    Δf = 0 # Hz
    pos = [[x, 0, 0] for x in xpos]
    spins = map(pos -> Spin([0,0,1.0], 1, T1, T2, Δf, pos), pos)
    map(spin -> excitation!(spin, [rf; 0rf], 0, [grad -0.52grad], dt), spins)
    sig = map(spin -> spin.signal, spins)
    Mz = map(spin -> spin.M[3], spins)

    return sig ≈ vec(answer["sig"]) && Mz ≈ vec(answer["mz"])

end

function testA5bMC()

    answer = matread("matlabtestdata/testA5b.mat")

    T1 = 600 # ms
    T2 = 100 # ms
    Δf = 10 # Hz
    M = [1.0, 0, 0]
    dt = 1 # ms
    time = 0:dt:1000 # ms
    spin = SpinMC(M, 1, [1], [T1], [T2], [Δf], Vector{Int}(), [0,0,0])
    mag = zeros(3,length(time))
    mag[:,1] = M
    for t = 2:length(time)
      freeprecess!(spin, dt)
      mag[:,t] = spin.M
    end

    return mag ≈ answer["M"]

end

function testF1cMC()

    answer = matread("matlabtestdata/testF1c.mat")

    t = 2.3 # ms
    α = π/4 # rad
    θ = 0 # rad
    Δf = 0 # Hz
    T1 = 600 # ms
    T2 = 100 # ms
    grad = [0.1, 0, 0] # G/cm
    xpos = -2:0.01:2 # cm
    pos = [[x, 0, 0] for x in xpos]
    spins = map(pos -> SpinMC([0,0,1.0], 1, [1], [T1], [T2], [Δf], Vector{Int}(), pos), pos)
    excitation!.(spins, θ, α)
    map(spin -> freeprecess!(spin, t, grad), spins)
    excitation!.(spins, θ, α)
    map(spin -> freeprecess!(spin, t, -grad/2), spins)
    sig = map(spin -> spin.signal, spins)

    return sig ≈ vec(answer["sig"])

end

function testF2cMC()

    answer = matread("matlabtestdata/testF2c.mat")

    dt = 0.04 # ms
    t = 0:dt:6 # ms
    rf = 0.05 * sinc.(t .- 3) # G
    Δf = -1000:20:1000 # Hz
    spins = map(Δf -> SpinMC([0,0,1.0], 1, [1], [600], [100], [Δf], Vector{Int}(), [0,0,0]), Δf)
    map(spin -> excitation!(spin, rf, 0, zeros(3,length(rf)), dt), spins)
    sig = map(spin -> spin.signal, spins)

    return sig ≈ vec(answer["sig"])

end

function testF3cMC()

    answer = matread("matlabtestdata/testF3c.mat")

    dt = 0.05 # ms
    t = 0.1:dt:6 # ms
    T = length(t)
    xpos = -2:0.1:2 # cm
    rf = 0.05 * sinc.(t .- 3) # G
    grad = [0.1ones(T)'; zeros(T)'; zeros(T)']
    T1 = 600 # ms
    T2 = 100 # ms
    Δf = 0 # Hz
    pos = [[x, 0, 0] for x in xpos]
    spins = map(pos -> SpinMC([0,0,1.0], 1, [1], [T1], [T2], [Δf], Vector{Int}(), pos), pos)
    map(spin -> excitation!(spin, [rf; 0rf], 0, [grad -0.52grad], dt), spins)
    sig = map(spin -> spin.signal, spins)
    Mz = map(spin -> spin.M[3], spins)

    return sig ≈ vec(answer["sig"]) && Mz ≈ vec(answer["mz"])

end

function testF3dMC()

    answer = matread("matlabtestdata/testF3d.mat")

    dt = 0.05 # ms
    t = 0.1:dt:6 # ms
    T = length(t)
    xpos = -2:0.1:2 # cm
    rf = 0.05 * sinc.(t .- 3) # G
    grad = [0.1ones(T)'; zeros(T)'; zeros(T)']
    T1 = 600 # ms
    T2 = 100 # ms
    Δf = 100 # Hz
    pos = [[x, 0, 0] for x in xpos]
    spins = map(pos -> SpinMC([0,0,1.0], 1, [1], [T1], [T2], [Δf], Vector{Int}(), pos), pos)
    map(spin -> excitation!(spin, [rf; 0rf], 0, [grad -0.52grad], dt), spins)
    sig = map(spin -> spin.signal, spins)
    Mz = map(spin -> spin.M[3], spins)

    return sig ≈ vec(answer["sig"]) && Mz ≈ vec(answer["mz"])

end

function testF3fMC()

    answer = matread("matlabtestdata/testF3f.mat")

    dt = 0.05 # ms
    t = 0.1:dt:6 # ms
    T = length(t)
    xpos = -5:0.1:5 # cm
    rf = 0.05 * sinc.(t .- 3) .* (exp.(im * 2π * 900 * t/1000) + exp.(-im * 2π * 900 * t/1000)) # G
    grad = [0.1ones(T)'; zeros(T)'; zeros(T)']
    T1 = 600 # ms
    T2 = 100 # ms
    Δf = 0 # Hz
    pos = [[x, 0, 0] for x in xpos]
    spins = map(pos -> SpinMC([0,0,1.0], 1, [1], [T1], [T2], [Δf], Vector{Int}(), pos), pos)
    map(spin -> excitation!(spin, [rf; 0rf], 0, [grad -0.52grad], dt), spins)
    sig = map(spin -> spin.signal, spins)
    Mz = map(spin -> spin.M[3], spins)

    return sig ≈ vec(answer["sig"]) && Mz ≈ vec(answer["mz"])

end

@testset "Spin" begin

    @testset "Single Compartment" begin

        @test testA5b()
        @test testF1a()
        @test testF1b()
        @test testF1c()
        @test testF2c()
        @test testF3a()
        @test testF3c()
        @test testF3d()
        @test testF3f()

    end

    @testset "Multicompartment" begin

        @test testA5bMC()
        @test testF1cMC()
        @test testF2cMC()
        @test testF3cMC()
        @test testF3dMC()
        @test testF3fMC()

    end

end


text:
cd("/Users/florian.oswald/Dropbox/teaching/ScPo/ScPo-CompEcon/ScPoExample.jl")
include("src/ScPoExample.jl")
include("test/runtests.jl")

@step ScPoExample.dbg_fun(rand(10))


text:
abstract type AbstractBernstein{T,X} <: AbstractSpecialPolynomial{T,X} end

"""

    Bernstein{N, T, X}

A [Bernstein  polynomial](https://en.wikipedia.org/wiki/Bernstein_polynomial) is a polynomial expressed in terms of
Bernstein basic polynomials. For each degree, `𝐍`, this is a set of `𝐍+1` degree `𝐍` polynomials of the form:
`β_{𝐍,ν} =  (ν choose 𝐍) x^ν  (1-x)^{𝐍-ν}`, `0 ≤ x ≤ 1.`

The `Bernstein{𝐍,T}` type represents a polynomial of degree `𝐍` or
less with a linear combination of the basis vectors using coefficients
of type `T`.

```jldoctest Bernstein
julia> using Polynomials, SpecialPolynomials

julia> p = basis(Bernstein{3},  2)
Bernstein(1⋅β₂,₂(x))

julia> convert(Polynomial, p)
Polynomials.Polynomial(1.0*x^2)
```

!!! note
    [StaticUnivariatePolynomials](https://github.com/tkoolen/StaticUnivariatePolynomials.jl) Offers a  more  performant version.

"""
struct Bernstein{𝐍,T,X} <: AbstractBernstein{T,X}
    coeffs::Vector{T}
    function Bernstein{𝐍,T,X}(coeffs::AbstractVector{T}) where {𝐍,T,X}
        N = findlast(!iszero, coeffs)
        N == nothing && return new{𝐍,T,X}(zeros(T, 0))
        (N > 𝐍 + 1) && throw(ArgumentError("Wrong length for coefficents"))
        return new{𝐍,T,X}(coeffs[1:N])
    end
end

export Bernstein
Polynomials.@registerN Bernstein 𝐍
Polynomials.:⟒(::Type{<:Bernstein{𝐍}}) where {𝐍} = Bernstein{𝐍}

function Polynomials.showterm(
    io::IO,
    ::Type{Bernstein{N,T,X}},
    pj::T,
    var,
    j,
    first::Bool,
    mimetype,
) where {N,T,X}
    iszero(pj) && return false
    !first && print(io, " ")
    print(io, Polynomials.hasneg(T) && Polynomials.isneg(pj) ? "- " : (!first ? "+ " : ""))
    print(io, "$(abs.(pj))⋅β")
    Polynomials.unicode_subscript(io, N)
    print(io, ",")
    Polynomials.unicode_subscript(io, j)
    print(io, "($var)")
    return true
end

# add default case with no N specified
function Bernstein(coeffs::AbstractVector{T}, var::Symbol=:x) where {T}
    l = findlast(!iszero, coeffs)
    𝐍 = l == nothing ? -1 : l - 1
    Bernstein{𝐍,T,Symbol(var)}(coeffs[1:(𝐍 + 1)])
end

# promotion between family
function Base.promote(
    p::P,
    q::Q,
) where {𝐍,T,X,P<:Bernstein{𝐍,T,X},𝐌,S,Y,Q<:Bernstein{𝐌,S,Y}}
    X == Y || throw(ArgumentError("different indeterminates"))
    MN = max(𝐌, 𝐍)
    R = promote_type(T, S)
    PQ = Bernstein{MN,R}
    convert(PQ, p), convert(PQ, q)
end

# promote up
function Base.convert(P::Type{<:Bernstein{𝐍}}, p::Bernstein{𝐌,S}) where {𝐍,𝐌,S}
    @assert 𝐍 >= 𝐌
    𝐑 = 𝐍 - 𝐌
    R = eltype(one(promote_type(eltype(P), S)) / 1)
    cs = zeros(R, 𝐍 + 1)
    for j in 0:𝐌
        pⱼ = p[j]
        iszero(pⱼ) && continue
        for i in j:(j + 𝐑)
            cs[1 + i] += pⱼ * binomial(𝐌, j) * binomial(𝐑, i - j) / binomial(𝐍, i)
        end
    end
    X = Polynomials.indeterminate(P, p)
    Bernstein{𝐍,R,X}(cs)
end

function Base.convert(
    P::Type{<:Polynomials.StandardBasisPolynomial},
    p::Bernstein{N,T},
) where {N,T}
    x = variable(P)
    p(x)
end

function Base.convert(C::Type{<:Bernstein{𝐍,T}}, p::Polynomial) where {𝐍,T}
    @assert degree(p) <= 𝐍

    R = eltype(one(T) / one(Int))
    cs = zeros(R, 𝐍 + 1)

    for (i, a_nu) in enumerate(coeffs(p))
        k = i - 1
        nk = binomial(𝐍, k)
        for j in k:𝐍
            cs[j + 1] += a_nu / nk * binomial(j, k)
        end
    end

    X = Polynomials.indeterminate(C, p)
    Bernstein{𝐍,R,X}(cs)
end

function Base.convert(P::Type{<:AbstractBernstein}, p::Polynomial{T}) where {T}
    𝐍 = degree(p)
    X = Polynomials.indeterminate(P, p)
    convert(Bernstein{𝐍,T,X}, p)
end

Polynomials.domain(::Type{<:AbstractBernstein}) = Polynomials.Interval(0, 1)
Polynomials.degree(p::Bernstein{N,T,X}) where {N,T,X} = degree(convert(Polynomial{T,X}, p))

Polynomials.constantterm(p::Bernstein) = p[0] # β.,ᵢ = 0 + ... when i > 0

## need to do one  and variable differently,as coefficients have length N+1
Base.one(P::Type{Bernstein{N,T,X}}) where {N,T,X} = Bernstein{N,T,X}(ones(T, N + 1))
Base.one(P::Type{Bernstein{N,T}}, var::Polynomials.SymbolLike) where {N,T} =
    Bernstein{N,T,Symbol(var)}(ones(T, N + 1))
function Base.one(P::Type{<:Bernstein{N}}, var::Polynomials.SymbolLike=:x) where {N}
    N′ = max(N, 0)
    one(Bernstein{N′,eltype(P),Symbol(var)})
end

function Polynomials.variable(P::Type{Bernstein{𝐍,T,X}}) where {𝐍,T,X}
    𝐍 <= 0 && throw(ArgumentError("Need  𝐍 ≥ 1"))
    R = typeof(one(T) / 1)
    Bernstein{𝐍,R,X}([1 - (i * one(T)) / 𝐍 for i in 𝐍:-1:0])
end

function Polynomials.variable(
    P::Type{<:Bernstein{𝐍}},
    var::Polynomials.SymbolLike=:x,
) where {𝐍}
    S = eltype(P)
    variable(Bernstein{𝐍,S,Symbol(var)})
end
Polynomials.variable(P::Type{Bernstein}, var::Polynomials.SymbolLike=:x) =
    variable(Bernstein{1,eltype(P),Symbol(var)})

function basis(::Type{P}, k::Int) where {𝐍,T,X,P<:Bernstein{𝐍,T,X}}
    zs = zeros(Int, k + 1)
    zs[end] = 1
    P(zs)
end

function basis(P::Type{<:Bernstein}, k::Int, _var::Polynomials.SymbolLike=:x; var=_var)
    zs = zeros(Int, k + 1)
    zs[end] = 1
    Bernstein(zs, var)
end

# poly evaluation
Polynomials.evalpoly(x, p::Bernstein) = simple_deCasteljau_polyval(p, x)

# we could  use the compensated one
# from https://doi.org/10.1016/j.camwa.2010.05.021
# Accurate evaluation of a polynomial and its derivative in Bernstein form
# Hao Jiang, Shengguo Li, Lizhi Cheng, and Fang Su
#
# for increased accuracy, but is it really necessary here?
function simple_deCasteljau_polyval(p::Bernstein{𝐍,T,X}, t::S) where {𝐍,T,X,S}
    p₀ = p[0]
    R = eltype(p₀ * t)
    zR = p₀ * zero(t)

    R = promote_type(T, S)
    𝐍 == -1 && return zR
    𝐍 == 0 && return p[0] * one(t)
    n = length(coeffs(p))
    iszero(n) && return zR

    bs = [zR for _ in 1:(𝐍 + 1)]
    for i in eachindex(coeffs(p))
        bs[i] = p[i - 1]
    end

    r = 1 - t
    for j in 1:𝐍
        for i in 0:(𝐍 - j)
            ii = i + 1
            bs[ii] = bs[ii] * r + bs[ii + 1] * t
        end
    end
    bs[1]
end

# function deCasteljau_polyval(p::Bernstein{N,T}, t) where {N, T}
#     b̂ = coeffs(p) .* one(t)
#     errb̂ = 0 .* b̂ #zeros(eltype(b̂), N+1)
#     r̂,ρ = twoSum(1, -t)
#     for j in 1:N
#         for ii in 0:(N-j)
#             i = ii+1
#             s, π̂ = twoProd(b̂[i], r̂)
#             v, σ̂ = twoProd(b̂[i+1], t)
#             b̂[i],β̂ = twoSum(s, v)
#             ŵ = π̂ .+ σ̂ .+ β̂ .+ (b̂[i] .* ρ)
#             errb̂[i] = errb̂[i] * r̂ + (errb̂[i+1] * t + ŵ)
#         end
#     end
#     b̂[1]+ errb̂[1]

# end

# # compensated form to  find p'(t)
# function deCasteljau_derivative_val(p::Bernstein{N,T}, t) where {N,T}
#     R = eltype(one(T)*one(eltype(t)))
#     r̂,ρ  = twoSum(one(R), -t)
#     Δb = zeros(R, N).*one(t)
#     errΔb = similar(Δb) #zeros(R, N).*one(t)
#     for ii=0:N-1
#         i = ii+1
#         Δb[i], errΔb[i] = twoSum(p[ii+1], -p[ii])
#     end
#     for j in 1:N-1
#         for ii in 0:N-1-j
#             i = ii+1
#             s, π̂ = twoProd(Δb[i], r̂)
#             v, σ̂ = twoProd(Δb[i+1], t)
#             Δb[i], β̂ = twoSum(s, v)
#             ŵ = π̂ + σ̂ + β̂ + (Δb[i] * ρ)
#             errΔb[i] = errΔb[i] * r̂ + (errΔb[i+1] * t + ŵ)
#         end
#     end
#     (Δb[1] + errΔb[1]) * N
# end

# # Moller-Knuth from Exact computations with approximate arithmetic
# # http://perso.ens-lyon.fr/jean-michel.muller/
# function twoSum(a::Number, b::Number)
#     s = a + b
#     â = s - b
#     b̂ = s - â
#     δ_a = a - â
#     δ_b = b - b̂
#     s, δ_a + δ_b
# end
# twoSum(x, y) = (x + y, zero(x)+zero(y))

# function twoProd(a::Number, b::Number)
#     ab = a * b
#     ab, fma(a,b,-ab)
# end
# twoProd(x, y)  = x*y, 0

# not needed, but speeds things along
function Base.:+(p1::P, p2::Q) where {𝐍,T,X,P<:Bernstein{𝐍,T,X},𝐌,S,Q<:Bernstein{𝐌,S,X}}
    p, q = promote(p1, p2)
    𝐎, R = length(p.coeffs), eltype(p)
    P′ = Bernstein{𝐎,R,X}
    return Bernstein([p[i] + q[i] for i in 0:𝐎], X)
end

# no promote(p1,p2)  called here
function Base.:*(p::P, q::Q) where {𝐍,T,X,P<:Bernstein{𝐍,T,X},𝐌,S,Y,Q<:Bernstein{𝐌,S,Y}}
    ## use b(n,k) * b(m,j) = choose(n,k)choose(m,j)/choose(n+m,k+j) b(n+m, k+j)

    isconstant(p) && return q * constantterm(p)
    isconstant(q) && return p * constantterm(q)
    Polynomials.assert_same_variable(p, q) ||
        throw(ArgumentError("p1 and p2 must have the same symbol"))

    R = typeof(one(promote_type(T, S)) / 1)
    x = variable(Polynomial{R})
    return convert(Bernstein, p(x) * q(x))

    ## multiplication  through conversion  is as  efficient
    # n,  m =  degree(p), degree(q)
    # 𝐍 = n + m
    # cs = zeros(R, 𝐍 + 1)

    # for i in 0:𝐍
    #     for j in max(0,i-m):min(n, i)
    #         cs[1+i] += p[j] * q[i-j] * (one(R) * binomial(n,j)  *  binomial(m,i-j))
    #     end
    #     cs[1+i] /= binomial(m+n, i)
    # end

    # Bernstein{𝐍,R}(cs, p.var)

end

#  use p'(t) = n ∑ (b_{i+1} - b_i) B_i,n-1
function Polynomials.derivative(p::Bernstein{𝐍,T,X}, order::Integer=1) where {𝐍,T,X}
    cs = coeffs(p)
    𝐍 < -1 && return p
    iszero(𝐍) && return zero(p)
    order == 0 && return p

    cs = zeros(T, 𝐍)
    dp = Bernstein{𝐍 - 1,T,X}(𝐍 * diff(coeffs(p)))

    order > 1 ? derivative(dp, order - 1) : dp
end

# use ∫b_n,nu  = 1/(n+1)  * ∑_{nu+1}^{n+1} b_n+1, j
function Polynomials.integrate(p::Bernstein{𝐍,T,X}) where {𝐍,T,X}
    R = eltype(one(T) / 1)
    𝐍𝐍 = 𝐍 + 1

    cs = zeros(R, 𝐍𝐍 + 1)

    @inbounds for ν in eachindex(p)
        for j in (ν + 1):𝐍𝐍
            cs[1 + j] += p[ν] / 𝐍𝐍
        end
    end

    ∫p = Bernstein{𝐍𝐍,R,X}(cs)
    return ∫p
end

#  could do p 33 of http://mae.engr.ucdavis.edu/~farouki/bernstein.pdf
function Base.divrem(num::Bernstein{𝐍,T}, den::Bernstein{𝐌,S}) where {𝐍,T,𝐌,S}
    R = eltype((one(T) + one(S)) / 1)
    p1 = convert(Polynomial{R}, num)
    p2 = convert(Polynomial{R}, den)
    q, r = divrem(p1, p2)
    convert.(Bernstein, (q, r))
end

# Replace me. Doesn't seem necessary and poor stability of conversion.
# cf https://arxiv.org/pdf/1910.01998.pdf for one possible way, though for approximate GCD
function Base.gcd(num::Bernstein{𝐍,T}, den::Bernstein{𝐌,S}) where {𝐍,T,𝐌,S}
    ps = convert.(Polynomial, (num, den))
    qs = gcd(ps...)
    convert.(Bernstein, qs)
end

function Polynomials.vander(
    P::Type{<:AbstractBernstein},
    xs::AbstractVector{T},
    n::Integer,
) where {T<:Number}
    N = length(xs) - 1 # xs = [x0, x1, ...]
    R = typeof(one(T) / one(T))
    V = zeros(R, N + 1, n + 1)
    for j in 0:n
        bnj = binomial(n, j)
        for i in 0:N
            x = xs[i + 1]
            V[i + 1, j + 1] = bnj * x^j * (1 - x)^(n - j)     #beta(n,j)(xj)
        end
    end
    V
end

# direct compuation of roots using numerical stable method of y Guðbjörn and Jónsson as
# noted by Corless, Sevyeri  in §2.1 of
# https://arxiv.org/pdf/1910.01998.pdf
function Polynomials.roots(p::Bernstein{𝐍,T}) where {𝐍,T}
    R = eltype(one(T) / 1)

    as = zeros(R, 𝐍 + 1)
    for (i, pᵢ) in pairs(coeffs(p))
        as[i] = pᵢ
    end

    Ap = diagm(-1 => ones(R, 𝐍 - 1))
    for j in 1:𝐍
        Ap[1, j] = -as[end - j]
    end

    Bp = diagm(-1 => ones(R, 𝐍 - 1))
    for j in 1:𝐍
        Bp[1, j] = -as[end - j]
    end

    for j in 2:𝐍
        Bp[j, j] = j / (𝐍 - j + 1)
    end
    Bp[1, 1] += as[end] / 𝐍

    eigvals(Ap, Bp)
end

# http://www.cecm.sfu.ca/personal/pborwein/MITACS/papers/LinMatPol.pdf
# Sec 4.1

function comrade_matrix(p::P) where {𝐍,T,P<:Bernstein{𝐍,T}}
    C₀, C₁ = comrade_pencil(p)
    C₀ * inv(C₁)
end

function comrade_pencil(p::P) where {𝐍,T,P<:Bernstein{𝐍,T}}
    n = 𝐍
    a, b = 0, 1
    C₀ = diagm(
        0 => a / (b - a) .* ((n:-1:1) ./ (1:n)) .* ones(T, n),
        -1 => b / (b - a) .* ones(T, n - 1),
    )
    C₁ = diagm(
        0 => 1 / (b - a) .* ((n:-1:1) ./ (1:n)) .* ones(T, n),
        -1 => 1 / (b - a) .* ones(T, n - 1),
    )

    C₀[end, end] *= p[end]
    C₁[end, end] *= p[end]

    for i in 0:(n - 1)
        C₀[1 + i, end] -= b / (b - a) * p[i]
        C₁[1 + i, end] -= 1 / (b - a) * p[i]
    end
    #    C₁[1,end] = -(n-1)/(b-a) * p[0]

    C₀, C₁
end

# block version
function comrade_pencil(p::P) where {𝐍,T,M<:AbstractMatrix{T},P<:Bernstein{𝐍,M}}
    m, m′ = size(p[0])
    @assert m == m′

    Iₘ = I(m)
    n = 𝐍

    a, b = 0, 1 # could generalize with Bernsteing
    R = eltype(one(T) / 1)
    C₀ = zeros(R, n * m, n * m)
    C₁ = zeros(R, n * m, n * m)
    Δ = 1:m

    for i in 1:(n - 1)
        C₀[(i - 1) * m .+ Δ, (i - 1) * m .+ Δ] .= (n + 1 - i) / i * a / (b - a) .* Iₘ
        C₀[i * m .+ Δ, (i - 1) * m .+ Δ] .= b / (b - a) .* Iₘ

        C₁[(i - 1) * m .+ Δ, (i - 1) * m .+ Δ] .= (n + 1 - i) / i * 1 / (b - a) .* Iₘ
        C₁[i * m .+ Δ, (i - 1) * m .+ Δ] .= 1 / (b - a) .* Iₘ
    end

    C₀[(n - 1) * m .+ Δ, (n - 1) * m .+ Δ] .= 1 / n * a / (b - a) .* p[n]
    C₁[(n - 1) * m .+ Δ, (n - 1) * m .+ Δ] .= 1 / n * 1 / (b - a) .* p[n]

    for i in 0:(n - 1)
        C₀[i * m .+ Δ, (n - 1) * m .+ Δ] .-= b / (b - a) .* p[i]
        C₁[i * m .+ Δ, (n - 1) * m .+ Δ] .-= 1 / (b - a) .* p[i]
    end

    #    C₁[Δ, (n-1)*Δ] .= -(n-1) * 1 / (b-a) .* p[0]

    C₀, C₁
end


text:
using Test

@testset "Registrator" begin

include("server.jl")

# Travis CI gets rate limited easily unless we have access to an API key.
if get(ENV, "TRAVIS", "") == "true" && !haskey(ENV, "GITHUB_API_TOKEN")
    @info "Skipping web tests on Travis CI (credentials are unavailable)"
else
    include("webui.jl")
end
include("webui/gitutils.jl")

end


text:
using Test
using Nash
using Distributions
using Polyhedra
using CDDLib
using LinearAlgebra
using Combinatorics
using IterTools
using QuantEcon
using LightGraphs
using SymPy

game = random_2players_game(Normal(0,2),2,3)
@testset "games" begin
    ggame = generate_game([1 0; 0 1], [1 0; 0 1])
    @test typeof(ggame) == Dict{String,Array{Int64,2}}
    @test typeof(ggame["player1"]) == Array{Int64,2}

    p2game = random_2players_game(Normal(0,2),2,3)
    @test typeof(p2game["player1"]) == Array{Float64,2}
    @test typeof(p2game["player2"]) == Array{Float64,2}

    npgame = random_nplayers_game(Binomial(5,0.5),[2,3,3]);
    @test typeof(npgame["player1"]) == Array{Int64,3}
end

@testset "utilities" begin
    s = [[0.5,0.5],[0.75,0.25]]
    out = outer(s)
    @test typeof(out) == Array{Float64,2}
    @test out[1][1] == s[1][1]*s[2][1]

    s = [Sym.([0.5,0.5]),Sym.([0.75,0.25])]
    out = outer(s)
    @test typeof(out) == Array{Float64,2}
    @test out[1][1] == s[1][1]*s[2][1]
end

@testset "equilibira" begin
    s = [[0.5,0.5],[0.75,0.25]]
    game = random_nplayers_game(Binomial(5,0.5),[2,2]);
    payoffs = get_payoff(game,s)
    @test typeof(payoffs) == Dict{String,Float64}
    @test payoffs["player2"] == s[1]'*game["player2"]*s[2]

    bp1 = best_reply(generate_game([1 0; 0 1], [1 0; 0 1]), [[1, 0], [1, 0]], 1)
    bp2 = best_reply(generate_game([1 0; 0 1], [1 0; 0 1]), [[1, 0], [1, 0]], 2)
    bp1p = best_reply(generate_game([1 0; 0 1], [1 0; 0 1]), [[1, 0], [1, 0]], 1, 1/5)
    # bpch = best_reply(generate_game(Matrix(I,3,3), Matrix(I,3,3)), [[1/2,1/2,0],[1/3,1/3,1/3]], 1, return_val = "chull")
    @test bp1 == bp2
    @test bp1 == [1 0]
    @test bp2 isa Array{<:Real,2}
    @test bp1p isa Array{<:Real,2}
    # @test bpch isa CDDLib.Polyhedron{<:Real}

    # nashq = is_nash_q(generate_game([1 0 ; 0 1], [1 0; 0 1]), [[1, 0], [1, 0]])
    # @test typeof(nashq) == Dict{String,Bool}
    # @test nashq["Nash"] == true

    # game_hist = iterate_best_reply(game, s, 0.5, 7)
    # @test length(game_hist) == 8
    # @test game_hist isa Array{Any,1} #TODO in future not Any type
    # @test game_hist[1] isa Array{Array{Float64,1},1}
end

@testset "markov" begin
    game = generate_game(Matrix(I, 3, 3), Matrix(I, 3, 3))
    s = [[1, 0, 0],[0, 1, 0]]
    mc = game2markov(game, s)
    @test mc isa MarkovChain{Int64,Array{Int64,2},Array{Int64,1}}
    @test all(simulate(mc, 4, init = 1) .== 1)
    @test is_irreducible(mc) == false
    @test period(mc) == 2
    @test is_aperiodic(mc) == false
end

@testset "symmetric" begin
    game = random_symmetric_2players_game(Binomial(10,1/2),2)
    @test typeof(game) == Dict{String,AbstractArray{Int64,2}}
    @test game["player1"] == game["player2"]'

    game = generate_game([1 0; 0 1], [1 0; 1 1])
    x = symbols("x", real=true)
    res = find_symmetric_nash_equilibrium_2players_game(game, [x, 1-x])
    @test res == [0.0, 0.5, 1.0]

    actions_no_vector = [2, 2];
    graph = create_symmetries_graph(actions_no_vector)
    @test length(vertices(graph.graph)) == 8
    @test length(edges(graph.graph)) == 4
    @test check_equality_condition(graph,((1,1),1), ((1,1),2)) == true
    @test check_equality_condition(graph,((1,1),1), ((1,2),2)) == false
    @test length(find_all_equalities(graph)) == 4
end


text:
export AbstractTracerProvider,
    global_tracer_provider,
    global_tracer_provider!,
    provider,
    tracer,
    resource,
    Tracer,
    AbstractSpan,
    with_span,
    span_context,
    span_kind,
    parent_span_context,
    attributes,
    is_recording,
    span_status!,
    end_span!,
    span_name!,
    span_name,
    span_links,
    span_events,
    span_status

"""
A tracer provider is a part of an [`Tracer`](@ref). For each concrete tracer
provider, `resource` and `OpenTelemetryAPI.create_span(name::String, tracer::Tracer{<:YourCustomProvider})` should also be implemented.
"""
abstract type AbstractTracerProvider end

struct DummyTracerProvider <: AbstractTracerProvider end

resource(::DummyTracerProvider) = Resource()

const DUMMY_TRACER_PROVIDER = DummyTracerProvider()

const GLOBAL_TRACER_PROVIDER = Ref{AbstractTracerProvider}(DUMMY_TRACER_PROVIDER)

"""
    global_tracer_provider()

Get the global tracer provider.
"""
global_tracer_provider() = GLOBAL_TRACER_PROVIDER[]

"""
    global_tracer_provider!(p)

Set the global tracer provider to `p`.
"""
global_tracer_provider!(p) = GLOBAL_TRACER_PROVIDER[] = p

"""
    Tracer(;instrumentation_info=InstrumentationInfo(), provider=global_meter_provider())
"""
Base.@kwdef struct Tracer{P<:AbstractTracerProvider}
    instrumentation_info::InstrumentationInfo = InstrumentationInfo()
    provider::P = global_tracer_provider()
end

provider(t::Tracer) = t.provider

"""
Each concrete span should have the following interfaces implemented.

  - [`tracer`](@ref)
  - [`span_context`](@ref)
  - [`span_kind`](@ref)
  - [`parent_span_context`](@ref)
  - [`attributes`](@ref)
  - [`is_recording`](@ref)
  - [`start_time`](@ref)
  - [`end_time`](@ref)
  - [`end_span!`](@ref)
  - [`span_status!`](@ref)
  - [`span_status`](@ref)
  - [`span_name!`](@ref)
  - [`span_name`](@ref)
  - [`span_links`](@ref)
  - [`span_events`](@ref)
"""
abstract type AbstractSpan end

struct NonRecordingSpan <: AbstractSpan
    name::String
    span_context::SpanContext
    parent_span_context::Union{Nothing,SpanContext}
end

"""
    tracer(s::AbstractSpan)

Get the [`Tracer`](@ref) which generates the span `s`.
"""
tracer(::NonRecordingSpan) = Tracer(; provider = DUMMY_TRACER_PROVIDER)

"""
    provider(s::AbstractSpan)

Get the [`AbstractTracerProvider`](@ref) which generates the tracer that the
span `s` resides in.
"""
provider(s::AbstractSpan) = provider(tracer(s))

"""
    resource(s::AbstractSpan)

Get the associated resource of the span `s`. Fall back to `resource(provider(::AbstractSpan))`
"""
resource(s::AbstractSpan) = resource(provider(s))

"""
    is_recording([current_span()])

Returns `true` if this span `s` is recording information like [`Event`](@ref) operations, attribute modification using [`setindex!`](@ref), etc.
"""
is_recording() = is_recording(current_span())
is_recording(s::NonRecordingSpan) = false

"""
    attributes(s::AbstractSpan)

Return either [`StaticAttrs`](@ref) or [`DynamicAttrs`](@ref) in the span `s`.

!!! warning

    Do not modify the returned attributes directly. Users should always
    modify attributes through `(s::AbstractSpan)[key]=val` because it will check
    the `is_recording(s)` first.
"""
attributes(s::NonRecordingSpan) = StaticAttrs()

"""
    start_time(s::AbstractSpan)

Get the start time of span `s` in nanoseconds.
"""
start_time(::NonRecordingSpan) = UInt(0)

"""
    end_time(s::AbstractSpan)

Get the end time of span `s` in nanoseconds.
"""
end_time(::NonRecordingSpan) = UInt(0)

"""
    (s::AbstractSpan)[key] = val

Set the attributes in span `s`. Only valid when the span is not ended yet.
"""
function Base.setindex!(s::AbstractSpan, val, key)
    if is_recording(s)
        setindex!(attributes(s), val, key)
    else
        @debug "the span is not recording."
    end
end

"""
    Base.getindex(s::AbstractSpan, key)

Look up `key` in the attributes of the span `s`.
"""
Base.getindex(s::AbstractSpan, key) = getindex(attributes(s), key)

"""
    Base.haskey(s::AbstractSpan, key)

Check if the span `s` has the key in its attributes.
"""
Base.haskey(s::AbstractSpan, key) = haskey(attributes(s), key)

"""
    Base.push!([s::AbstractSpan], event::Event)

Add an [`Event`](@ref) into the span `s`.
"""
Base.push!(event::Event) = push!(current_span(), event)

function Base.push!(s::AbstractSpan, event::Event)
    if is_recording(s)
        push!(span_events(s), event)
    else
        @warn "the span is not recording."
    end
end

"""
    span_events(s::AbstractSpan)

Get the recorded events in a span. A [`Limited`](@ref) is expected.
"""
span_events() = span_events(current_span())
span_events(s::NonRecordingSpan) = Limited(Event[])

"""
    Base.push!([current_span()], link::Link)

Add a [`Link`](@ref) into the span `s`.
"""
Base.push!(link::Link) = push!(current_span(), link)

function Base.push!(s::AbstractSpan, link::Link)
    if is_recording(s)
        push!(span_links(s), link)
    else
        @warn "the span is not recording."
    end
end

"""
    span_links(s::AbstractSpan)

Get the recorded links in a span. A [`Limited`](@ref) is expected.
"""
span_links() = span_links(current_span())
span_links(s::NonRecordingSpan) = Limited(Link[])

"""
    span_status!([current_span()], code::SpanStatusCode, description=nothing)

Update the status of span `s` by following the original specification. `description` is only considered when the `code` is `SPAN_STATUS_ERROR`. Only valid when the span is not ended yet.
"""
span_status!(code::SpanStatusCode, description) =
    span_status!(current_span(), code, description)
span_status!(s::NonRecordingSpan, code::SpanStatusCode, description = nothing) =
    @warn "the span is not recording."

"""
    span_status([current_span()])

Get status of the span. A [`SpanStatusCode`](@ref) is returned.
"""
span_status() = span_status(current_span())
span_status(s::NonRecordingSpan) = SpanStatus(SPAN_STATUS_UNSET)

"""
    end_span!([s=current_span()], [t=UInt(time()*10^9)])

Set the end time of the span and trigger span processors. Note `t` is the nanoseconds.
"""
end_span!() = end_span!(current_span())
end_span!(s::AbstractSpan) = end_span!(s, UInt(time() * 10^9))
end_span!(s::NonRecordingSpan) = @debug "the span is not recording."
end_span!(s::NonRecordingSpan, t) = @debug "the span is not recording."

"""
    Base.push!(s::AbstractSpan, ex::Exception; is_rethrow_followed = false)

A specialized variant of [`Event`](@ref) to record exceptions. Usually used in a `try... catch...end` to capture the backtrace. If the `ex` is `rethrow`ed in the `catch...end`, `is_rethrow_followed` should be set to `true`.
"""
Base.push!(s::NonRecordingSpan, ex::Exception; is_rethrow_followed = false) =
    @warn "the span is not recording."

"""
    span_name([current_span()])
"""
span_name() = span_name(current_span())
span_name(s::NonRecordingSpan) = s.name

"""
    span_name!([current_span()], name::String)
"""
span_name!(name::String) = span_name!(current_span(), name)
span_name!(s::NonRecordingSpan, name::String) = @warn "the span is not recording."

#####

"""
    span_context([s::AbstractSpan])

Get the [`SpanContext`](@ref) from a span `s`. If `s` is not specified,
[`current_span()`](@ref) will be used. `nothing` is returned if no span context
found.
"""
span_context() = span_context(current_context())
span_context(ctx::Context) = span_context(current_span(ctx))
span_context(::Nothing) = nothing
span_context(s::NonRecordingSpan) = s.span_context

"""
    span_kind([current_span()])

Return [`SpanKind`](@ref)
"""
span_kind() = span_kind(current_span())
span_kind(::NonRecordingSpan) = SPAN_KIND_INTERNAL

"""
    parent_span_context(s::AbstractSpan)

Get the [`SpanContext`](@ref) from parent span.
"""
parent_span_context() = parent_span_context(current_span())
parent_span_context(s::NonRecordingSpan) = s.parent_span_context

"""
Here we follow the [Behavior of the API in the absence of an installed SDK](https://github.com/open-telemetry/opentelemetry-specification/blob/main/specification/trace/api.md#behavior-of-the-api-in-the-absence-of-an-installed-sdk).
"""
function create_span(
    name::String,
    tracer::Tracer{DummyTracerProvider};
    context = current_context(),
    kw...,
)
    parent_span = current_span(context)
    if isnothing(parent_span)
        NonRecordingSpan(name, INVALID_SPAN_CONTEXT, nothing)
    elseif is_recording(parent_span)
        parent_span_ctx = span_context(parent_span)
        NonRecordingSpan(name, parent_span_ctx, parent_span_ctx)
    else
        parent_span
    end
end

"""
    with_span(f, name::String, [tracer=Tracer()];kw...)

Call function `f` with the current span set a newly created one of `name` with `tracer`.

# Keyword arguments

  - `end_on_exit=true`, controls whether to call [`end_span!`](@ref) after `f` or not.
  - `record_exception=true`, controls whether to record the exception.
  - `set_status_on_exception=true`, decides whether to set status to [`SPAN_STATUS_ERROR`](@ref) automatically when an exception is caught.
  - The rest keyword arguments are forwarded to [`create_span`](@ref).
"""
function with_span(
    f,
    name::String,
    tracer::Tracer = Tracer();
    end_on_exit = true,
    record_exception = true,
    set_status_on_exception = true,
    kw...,
)
    s = create_span(name, tracer; kw...)
    with_context(; SPAN_KEY_IN_CONTEXT => s) do
        try
            f()
        catch ex
            if is_recording(s)
                if record_exception
                    push!(s, ex; is_rethrow_followed = true)
                end
                if set_status_on_exception
                    span_status!(s, SPAN_STATUS_ERROR, string(ex))
                end
            end
            rethrow(ex)
        finally
            if end_on_exit
                end_span!(s)
            end
        end
    end
end


text:
include("../src/UNSflow.jl")
using UNSflow

cref = 1.
bref = 10.
sref = 10.

patch1 = patch(0., -5., 0., 0., "FlatPlate", 1., 0., 0.11, 10)
patch2 = patch(0., 0., 0., 0., "FlatPlate", 1., 0., 0.11, 10)
patch3 = patch(0., 5., 0., 0., "FlatPlate", 1., 0., 0.11, 5)
patchdata = [patch1; patch2; patch3]

alphadef = ConstDef(15.*pi/180)
hdef = ConstDef(0.)
udef = ConstDef(1.)
kin = KinemDef3D(alphadef, hdef, udef)

surf = ThreeDSurf(cref, bref, sref, patchdata, kin, 1., 70, 35, 21)
field = ThreeDFlowField()
dtstar = 0.015
nsteps =round(Int,5./dtstar) + 1
@time mat3d, surf3d, field3d = QSLLT_ldvm(surf, field, nsteps, dtstar)


text:
## We need to use a library for Symbolic calculations, a very common one is SymPy from the python
## language. Luckily Julia has excellent support for it.
using SymPy

#jl Remove the follow `#jl` to use this file standalone
#jl include("lv_model.jl")

# We want to deal with all of these as symbolic variables, which are not the default type in Julia,
# unlike Mathematica.
@vars u1 u2
@vars r1 r2 a12 a21 K1 K2
# We need to make a symbolic parameter list, as `LVPar` is numeric
spar = Dict(
    :α12 => a12,
    :α21 => a21,
    :r1 => r1,
    :r2 => r2,
    :K1 => K1,
    :K2 => K2);

# we need to make symbolic versions of the model equations. We do this by calling the function with
# the symbolic parameters. the last parameter could be anythign as the time (`t`) argument is not
# used. I have set it to `NaN` which is a name for not a number.
f1, f2 = lv_comp([u1, u2], spar, NaN)
#-
SymPy.solve(f1, u1)
#-
SymPy.solve(f1, u2)
#-
SymPy.solve(f2, u1)
#-
SymPy.solve(f2, u2)
#-

#### Another version of above (does the same thing but not as elegant)

# x, y, r1, r2, k1, k2, a12, a21 = symbols("x, y, r1, r2, k1, k2, a12, a21", real = true)
#
# f(x, y) = x * r1 * (1 - (x + a12 * y) / k1)
# g(x ,y) = y * r2 * (1 - (y + a21 * x) / k2)
#
# SymPy.solve(f(x,y),x)
# SymPy.solve(f(x,y),y)
# SymPy.solve(g(x,y),x)
# SymPy.solve(g(x,y),y)

max_stars_repo_name: McCannLab/TheoryClub

max_issues_repo_name: McCannLab/TheoryClub

max_forks_repo_name: McCannLab/TheoryClub


text:
# Toy array type for simple tests
# Actual applications should use ROCArray

mutable struct HSAArray{T,N} <: AbstractArray{T,N}
    size::Dims{N}
    buffer::Mem.Buffer
end

# TODO: Support non-isbitstype allocations
function HSAArray(agent::HSAAgent, ::Type{T}, size::NTuple{N,Int}) where {T,N}
    @assert isprimitivetype(T) "$T is not a primitive type"
    @assert all(x->x>0, size) "Invalid array size: $size"

    nbytes = sizeof(T) * prod(size)

    # Always allocated as coherent memory since fill(), getindex(), and setindex!()
    # access the buffer pointer directly.
    buffer = Mem.alloc(agent, nbytes; coherent=true)

    arr = HSAArray{T,N}(size, buffer)

    hsaref!()
    finalizer(arr) do arr
        Mem.free(buffer)
        hsaunref!()
    end

    return arr
end

HSAArray(::Type{T}, size::NTuple{N,Int}) where {T,N} =
    HSAArray(DEFAULT_AGENT[], T, size)

### HSAArray Methods ###

HSAArray(arr::Array{T,N}) where {T,N} =
    HSAArray(DEFAULT_AGENT[], arr)

function HSAArray(agent::HSAAgent, arr::Array{T,N}) where {T,N}
    harr = HSAArray(agent, T, size(arr))
    Mem.upload!(harr.buffer, arr)
    return harr
end

function Array(harr::HSAArray{T,N}) where {T,N}
    array = Array{T}(undef, size(harr))
    Mem.download!(array, harr.buffer)
    array
end

Base.pointer(arr::HSAArray{T,N}) where {T, N} = Ptr{T}(arr.buffer.ptr)
Base.IndexStyle(::Type{<:HSAArray}) = Base.IndexLinear()
Base.IndexStyle(::HSAArray) = Base.IndexLinear()

function Base.iterate(A::HSAArray, i=1) # copy-pasta from Base
    Base.@_inline_meta
    (i % UInt) - 1 < length(A) ? (@inbounds A[i], i + 1) : nothing
end

Base.similar(arr::HSAArray{T,N}) where {T,N} =
    HSAArray(T, size(arr))
Base.similar(agent, arr::HSAArray{T,N}) where {T,N} =
    HSAArray(agent, T, size(arr))

Base.size(arr::HSAArray) = arr.size
Base.length(arr::HSAArray) = sum(size(arr))

function Base.fill!(arr::HSAArray{T,N}, value::T) where {T,N}
    for idx in 1:length(arr)
        arr[idx] = value
    end
end

@inline function Base.getindex(arr::HSAArray{T,N}, idx) where {T,N}
    @boundscheck checkbounds(arr, idx)
    Base.unsafe_load(pointer(arr), idx)::T
end

@inline function Base.setindex!(arr::HSAArray{T,N}, value, idx) where {T,N}
    @boundscheck checkbounds(arr, idx)
    Base.unsafe_store!(pointer(arr), value, idx)
end


text:
export ASB07_decomposed

struct ASB07_decomposed <: ContinuousPost
    options::TwoLayerOptions

    function ASB07_decomposed(𝑂::Options)
        𝑂new = validate_and_wrap_options(𝑂, options_ASB07_decomposed())
        return new(𝑂new)
    end
end

# convenience constructor from pairs of symbols
ASB07_decomposed(𝑂::Pair{Symbol,<:Any}...) =
    ASB07_decomposed(Options(Dict{Symbol,Any}(𝑂)))

# default options (they are added in the function validate_and_wrap_options)
ASB07_decomposed() = ASB07_decomposed(Options())

include("options.jl")
include("init.jl")
include("post.jl")
include("reach.jl")


text:
using ITensors
using PastaQ
using Test

@testset "depracated" begin
  N = 5

  sites = siteinds("Qubit", N)

  q = qubits(N)
  @test q isa MPS
  @test length(q) == N

  q = qubits(sites)
  @test q isa MPS
  @test length(q) == N
  @test q ≈ productstate(sites)

  ψ = randomstate(sites)
  q = qubits(ψ)
  @test q isa MPS
  @test length(q) == N
  @test siteinds(ψ) == siteinds(q)

  q = qubits(sites, ["Z+","Z-","X+","Y-","X+"])
  @test q isa MPS
  @test length(q) == N
  @test q ≈ productstate(sites, ["Z+","Z-","X+","Y-","X+"])
end


text:
using OptimalTransport
using Pkg: Pkg
using SafeTestsets

using Test

const GROUP = get(ENV, "GROUP", "All")

@testset "OptimalTransport" begin
    if GROUP == "All" || GROUP == "OptimalTransport"
        @safetestset "Utilities" begin
            include("utils.jl")
        end

        @testset "Entropically regularized OT" begin
            @safetestset "Sinkhorn Gibbs" begin
                include(joinpath("entropic", "sinkhorn_gibbs.jl"))
            end
            @safetestset "Stabilized Sinkhorn" begin
                include(joinpath("entropic", "sinkhorn_stabilized.jl"))
            end
            @safetestset "Sinkhorn with ε-scaling" begin
                include(joinpath("entropic", "sinkhorn_epsscaling.jl"))
            end
            @safetestset "Unbalanced Sinkhorn" begin
                include(joinpath("entropic", "sinkhorn_unbalanced.jl"))
            end
            @safetestset "Sinkhorn barycenter" begin
                include(joinpath("entropic", "sinkhorn_barycenter.jl"))
            end
            @safetestset "Sinkhorn divergence" begin
                include(joinpath("entropic", "sinkhorn_divergence.jl"))
            end
        end

        @safetestset "Quadratically regularized OT" begin
            include("quadratic.jl")
        end
    end

    # CUDA requires Julia >= 1.6
    if (GROUP == "All" || GROUP == "GPU") && VERSION >= v"1.6"
        # activate separate environment: CUDA can't be added to test/Project.toml since it
        # is not available on older Julia versions
        Pkg.activate("gpu")
        Pkg.instantiate()

        @safetestset "Simple GPU" begin
            include(joinpath("gpu/simple_gpu.jl"))
        end
    end
end


text:
export iauAper
"""
In the star-independent astrometry parameters, update only the
Earth rotation angle, supplied by the caller explicitly.

This function is part of the International Astronomical Union's
SOFA (Standards of Fundamental Astronomy) software collection.

Status:  support function.

Given:
   theta   double      Earth rotation angle (radians, Note 2)
   astrom  iauASTROM*  star-independent astrometry parameters:
      pmt    double       not used
      eb     double[3]    not used
      eh     double[3]    not used
      em     double       not used
      v      double[3]    not used
      bm1    double       not used
      bpn    double[3][3] not used
      along  double       longitude + s' (radians)
      xpl    double       not used
      ypl    double       not used
      sphi   double       not used
      cphi   double       not used
      diurab double       not used
      eral   double       not used
      refa   double       not used
      refb   double       not used

Returned:
   astrom  iauASTROM*  star-independent astrometry parameters:
      pmt    double       unchanged
      eb     double[3]    unchanged
      eh     double[3]    unchanged
      em     double       unchanged
      v      double[3]    unchanged
      bm1    double       unchanged
      bpn    double[3][3] unchanged
      along  double       unchanged
      xpl    double       unchanged
      ypl    double       unchanged
      sphi   double       unchanged
      cphi   double       unchanged
      diurab double       unchanged
      eral   double       "local" Earth rotation angle (radians)
      refa   double       unchanged
      refb   double       unchanged

Notes:

   1. This function exists to enable sidereal-tracking applications to
   avoid wasteful recomputation of the bulk of the astrometry
   parameters:  only the Earth rotation is updated.

   2. For targets expressed as equinox based positions, such as
   classical geocentric apparent (RA,Dec), the supplied theta can be
   Greenwich apparent sidereal time rather than Earth rotation
   angle.

   3. The function iauAper13 can be used instead of the present
   function, and starts from UT1 rather than ERA itself.

   4. This is one of several functions that inserts into the astrom
   structure star-independent parameters needed for the chain of
   astrometric transformations ICRS <-> GCRS <-> CIRS <-> observed.

   The various functions support different classes of observer and
   portions of the transformation chain:

         functions         observer        transformation

      iauApcg iauApcg13    geocentric      ICRS <-> GCRS
      iauApci iauApci13    terrestrial     ICRS <-> CIRS
      iauApco iauApco13    terrestrial     ICRS <-> observed
      iauApcs iauApcs13    space           ICRS <-> GCRS
      iauAper iauAper13    terrestrial     update Earth rotation
      iauApio iauApio13    terrestrial     CIRS <-> observed

   Those with names ending in "13" use contemporary SOFA models to
   compute the various ephemerides.  The others accept ephemerides
   supplied by the caller.

   The transformation from ICRS to GCRS covers space motion,
   parallax, light deflection, and aberration.  From GCRS to CIRS
   comprises frame bias and precession-nutation.  From CIRS to
   observed takes account of Earth rotation, polar motion, diurnal
   aberration and parallax (unless subsumed into the ICRS <-> GCRS
   transformation), and atmospheric refraction.

This revision:   2013 September 25

SOFA release 2018-01-30

Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""
function iauAper(theta::Real, astrom::iauASTROM)
   ref_astrom = Ref{iauASTROM}(astrom)
   ccall((:iauAper, libsofa_c), Cvoid,
            (Cdouble, Ref{iauASTROM}),
            convert(Float64, theta), ref_astrom)

   return ref_astrom[]
end


text:
module Interpolations

export linerinterp

function linerinterp(grid,vals)

    function funfun(x)

        if x < grid[1]

            return (vals[2]-vals[1])/(grid[2]-grid[1])*(x-grid[1])+vals[1]

        end

          if  grid[end] <= x

               return (vals[end]-vals[end-1])/(grid[end]-grid[end-1])*(x-grid[end-1])+vals[end-1]

        end

        if grid[1] <= x < grid[end]

            i = searchsortedlast(grid,x)

            y = (vals[i+1]-vals[i])/(grid[i+1]-grid[i])*(x-grid[i])+vals[i]

            return y

        end

    end

    return funfun

end

end


text:
module PProfTest

using PProf

using Test
using Profile
using ProtoBuf

const out = tempname() * ".pb.gz"

function foo(n, a, out=[])
    # make this expensive to ensure it's sampled
    for i in 1:n
        push!(out, i*a)
    end
end

@testset "empty output" begin
    Profile.clear()
    # Function doesn't error if no Profile recorded
    rm(out, force=true)
    @test (pprof(out=out, web=false); true)
    @test ispath(out)
end

@testset "export basic profile" begin
    Profile.clear()

    let x = 1
        @profile for _ in 1:10000; x += 1; end
        sleep(2)
    end

    # Cache profile output to test that it isn't changed
    _prior_profile_output = Profile.retrieve()

    # Write the profile
    outf = pprof(out=out, web=false)

    # Read the exported profile
    prof = open(io->readproto(io, PProf.perftools.profiles.Profile()), outf, "r")

    # Verify that we exported stack trace samples:
    @test length(prof.sample) > 0
    # Verify that we exported frame information
    @test length(prof.location) > 0
    @test length(prof._function) > 0

    # Test that we didn't modify the Profile output
    Profile.retrieve() == _prior_profile_output
end

function load_prof_proto(file)
    @show file
    open(io->readproto(io, PProf.perftools.profiles.Profile()), file, "r")
end

@testset "with_c" begin
    Profile.clear()

    let arr = []
        @profile foo(1000000, 5, arr)
        sleep(2)
    end
    for i in 1:2
        if i == 1
            data = Profile.fetch()
            args = (data,)
        else
            data,lidict = Profile.retrieve()
            args = (data, lidict)
        end

        # Write a profile that includes C function frames
        with_c = load_prof_proto(pprof(args..., out=tempname(), web=false, from_c = true))

        # Write a profile that excludes C function frames
        without_c = load_prof_proto(pprof(args..., out=tempname(), web=false, from_c = false))

        # Test that C frames were excluded
        @test length(with_c.sample) == length(without_c.sample)
        @test length(with_c.location) > length(without_c.location)
        @test length(with_c._function) > length(without_c._function)
    end
end

@testset "full_signatures" begin
    Profile.clear()
    @profile foo(1000000, 5, [])
    @test "foo" in load_prof_proto(pprof(out=tempname(), web=false, full_signatures = false)).string_table
    @test any(occursin.(Regex("^foo\\(::$Int, ::$Int, ::(Vector|Array)"),
                load_prof_proto(pprof(out=tempname(), web=false, full_signatures = true)).string_table))
end

@testset "drop_frames/keep_frames" begin
    @test load_prof_proto(pprof(out=tempname(), web=false, drop_frames = "foo")).drop_frames != 0
    @test load_prof_proto(pprof(out=tempname(), web=false, keep_frames = "foo")).keep_frames != 0
end

@testset "@pprof macro" begin

    @pprof foo(10000, 5, [])

    @test PProf.proc[] !== nothing
    @test process_running(PProf.proc[])

    PProf.kill()
    @test PProf.proc[] === nothing

    @test isfile("profile.pb.gz")
    rm("profile.pb.gz")
end

@testset "subprocess refresh" begin

    @pprof foo(10000, 5, [])

    current_proc = PProf.proc[]
    @test process_running(current_proc)

    PProf.refresh()
    sleep(1)

    @test process_running(PProf.proc[])
    @test process_exited(current_proc)
end

@testset "subprocess kill" begin

    @pprof foo(10000, 5, [])

    current_proc = PProf.proc[]
    @test process_running(current_proc)

    PProf.kill()
    sleep(1)

    @test process_exited(current_proc)
    @test PProf.proc[] === nothing
    @test isfile("profile.pb.gz")
    rm("profile.pb.gz")
end

@testset "enforce correct output file extension (.pb.gz)" begin
    dir = mktempdir()
    @test basename(pprof(out="$dir/test")) == "test.pb.gz"
    @test basename(pprof(out="$dir/test.pb.gz")) == "test.pb.gz"

    @test_throws AssertionError pprof(out="$dir/")  # directory path with no file throws err
end

end # module


text:
@doc "Free Space path loss Model"
type FreeSpaceModel<:Model
  freq::Float64
  FreeSpaceModel() = new(800)
end

@doc "Calculate the path loss with the FreeSpace Model for ERP"
function pathloss_op(model::FreeSpaceModel,dist::AbstractFloat)
       return 32.44 + 20*log10(dist) + 20*log10(model.freq);
end


text:
pwd() ∉ LOAD_PATH && push!(LOAD_PATH, pwd())
"C:\\Users\\Jones\\Desktop\\home\\Code\\julia\\RCS-Invert\\src" ∉ LOAD_PATH && push!(LOAD_PATH, "C:\\Users\\Jones\\Desktop\\home\\Code\\julia\\RCS-Invert\\src")
"C:\\Users\\Jones\\Desktop\\home\\Code\\julia\\RCS-Invert\\src\\TMFDTD" ∉ LOAD_PATH && push!(LOAD_PATH, "C:\\Users\\Jones\\Desktop\\home\\Code\\julia\\RCS-Invert\\src\\TMFDTD")

"C:\\Users\\Jones\\Desktop\\home\\Code\\julia\\RCS-Invert\\tests" ∉ LOAD_PATH && push!(LOAD_PATH, "C:\\Users\\Jones\\Desktop\\home\\Code\\julia\\RCS-Invert\\tests")
#"/home/zander/Code/julia/RCS-Invert/src" ∉ LOAD_PATH && push!(LOAD_PATH, "/home/zander/Code/julia/RCS-Invert/src")
#"/home/zander/Code/julia/RCS-Invert/src" ∉ LOAD_PATH && push!(LOAD_PATH, "/home/zander/Code/julia/RCS-Invert/src")
using FDTD_2D
using DiscreteAxis

function epsilon(r::AbstractVector{T}) where T
    d = [0.305, 0.4]
    if all(0.0.-d./2 .<= r .<= d./2)
        return convert(T, 2.3)
    else
        return one(T)
    end
end
dB(x::Number) = 20*log10(abs(x))

function Medium_Test2D(θ=0.0)
    fmax = 0.5*10^9
    f₀ = fmax/2
    nPML = 14
    space, time = setup_spacetime(fmax, 1.5, nPML;
                                    xlims = (-1.0,1.0),
                                    ylims = (-1.0,1.0))
    println("$(space.x.N), $(space.y.N), $time.N")

    println("Δt is $(time.Δ)")

    ###############################################################################
    # MEDIUM DEFINITION
    ###############################################################################

    @time medium = Medium(space; θ = θ, ε = epsilon) #freespace medium
    Esource = sinc.(4π.*fmax.*(time.-maximum(time.pts)/4))

    Hsource = zeros(time.N)
    interior = interior_range(space, nPML)
    sourceloc = [nPML+2, interior[2]]
    sourceindex = (2, sourceloc)
    detectorindex = (:,[:,:,:])

    @time field = FDTD_propagate(space, time, f₀, nPML; #Try with a PML
                                    source = (Esource, Hsource),
                                    medium = medium,
                                    sourceindex = sourceindex,
                                    detectorindex = detectorindex
                                    )
    #s_size = size(space)

    return field

end
field = Medium_Test2D(0.1)

using Plots
plotly()
heatmap(dB.(field[200][1][3]))


text:
@testset "Sudoku" begin

@testset "Sudoku from opensourc.es" begin
    grid = [0 2 1 0 7 9 0 8 5;
            0 4 5 3 1 0 0 0 9;
            0 7 0 0 4 0 0 1 0;
            0 0 0 1 0 8 0 3 6;
            0 6 0 0 0 0 2 0 8;
            0 0 0 0 0 3 0 0 4;
            6 0 8 0 0 0 0 0 0;
            0 9 4 0 0 7 8 0 0;
            2 0 0 5 0 0 0 4 0]

    m = Model(with_optimizer(CS.Optimizer))
    @variable(m, 1 <= x[1:9,1:9] <= 9, Int)
    # set variables
    for r=1:9, c=1:9
        if grid[r,c] != 0
            @constraint(m, x[r,c] == grid[r,c])
        end
    end
    # sudoku constraints
    jump_add_sudoku_constr!(m, x)

    optimize!(m)
    @test JuMP.termination_status(m) == MOI.OPTIMAL
    @test jump_fulfills_sudoku_constr(JuMP.value.(x))
end

@testset "Hard sudoku" begin
    com = CS.init()

    grid = zeros(Int,(9,9))
    grid[1,:] = [0 0 0 5 4 6 0 0 9]
    grid[2,:] = [0 2 0 0 0 0 0 0 7]
    grid[3,:] = [0 0 3 9 0 0 0 0 4]
    grid[4,:] = [9 0 5 0 0 0 0 7 0]
    grid[5,:] = [7 0 0 0 0 0 0 2 0]
    grid[6,:] = [0 0 0 0 9 3 0 0 0]
    grid[7,:] = [0 5 6 0 0 8 0 0 0]
    grid[8,:] = [0 1 0 0 3 9 0 0 0]
    grid[9,:] = [0 0 0 0 0 0 8 0 6]

    com_grid = create_sudoku_grid!(com, grid)
    add_sudoku_constr!(com, com_grid)

    @test CS.solve!(com, CS.get_default_options()) == :Solved
    @test fulfills_sudoku_constr(com_grid)
end

@testset "Hard sudoku infeasible" begin
    grid = [0 0 0 5 4 6 0 0 9;
            0 2 0 0 0 0 0 0 7;
            0 0 3 9 0 0 0 0 4;
            9 0 5 0 0 0 0 7 3;
            7 0 0 0 0 0 0 2 0;
            0 0 0 0 9 3 0 0 0;
            0 5 6 0 0 8 0 0 0;
            0 1 0 0 3 9 0 0 0;
            0 0 0 0 0 0 8 0 6]

    m = Model(with_optimizer(CS.Optimizer))
    @variable(m, 1 <= x[1:9,1:9] <= 9, Int)
    # set variables
    for r=1:9, c=1:9
        if grid[r,c] != 0
            @constraint(m, x[r,c] == grid[r,c])
        end
    end
    # sudoku constraints
    jump_add_sudoku_constr!(m, x)

    optimize!(m)
    @test JuMP.termination_status(m) == MOI.INFEASIBLE
    @test !jump_fulfills_sudoku_constr(JuMP.value.(x))
end

@testset "Hard fsudoku repo" begin
    com = CS.init()

    grid = zeros(Int,(9,9))
    grid[1,:] = [0 0 0 0 0 0 0 0 0]
    grid[2,:] = [0 1 0 6 2 0 0 9 0]
    grid[3,:] = [0 0 2 0 0 9 3 1 0]
    grid[4,:] = [0 0 4 0 0 6 0 8 0]
    grid[5,:] = [0 0 8 7 0 2 1 0 0]
    grid[6,:] = [0 3 0 8 0 0 5 0 0]
    grid[7,:] = [0 6 9 1 0 0 4 0 0]
    grid[8,:] = [0 8 0 0 7 3 0 5 0]
    grid[9,:] = [0 0 0 0 0 0 0 0 0]

    m = Model(with_optimizer(CS.Optimizer))
    @variable(m, 1 <= x[1:9,1:9] <= 9, Int)
    # set variables
    nvars_set = 0
    for r=1:9, c=1:9
        if grid[r,c] != 0
            @constraint(m, x[r,c] == grid[r,c])
            nvars_set += 1
        end
    end

    @test nvars_set == length(filter(n -> n != 0, grid))

    # sudoku constraints
    jump_add_sudoku_constr!(m, x)

    optimize!(m)

    @test JuMP.termination_status(m) == MOI.OPTIMAL

    # check that it actually solves the given sudoku
    for r=1:9, c=1:9
        if grid[r,c] != 0
            @test JuMP.value(x[r,c]) == grid[r,c]
        end
    end

    @test jump_fulfills_sudoku_constr(JuMP.value.(x))
end

@testset "Hard fsudoku repo 0-8" begin
    com = CS.init()

    grid = zeros(Int,(9,9))
    grid[1,:] = [0 0 0 0 0 0 0 0 0]
    grid[2,:] = [0 1 0 6 2 0 0 9 0]
    grid[3,:] = [0 0 2 0 0 9 3 1 0]
    grid[4,:] = [0 0 4 0 0 6 0 8 0]
    grid[5,:] = [0 0 8 7 0 2 1 0 0]
    grid[6,:] = [0 3 0 8 0 0 5 0 0]
    grid[7,:] = [0 6 9 1 0 0 4 0 0]
    grid[8,:] = [0 8 0 0 7 3 0 5 0]
    grid[9,:] = [0 0 0 0 0 0 0 0 0]
    grid .-= 1

    com_grid = Array{CS.Variable, 2}(undef, 9, 9)
    for (ind,val) in enumerate(grid)
        if val == -1
            if ind == 81 # bottom right
                # some other values are possible there
                com_grid[ind] = CS.add_var!(com, 9, 11)
            elseif ind == 80 # one above (will be 9 in the end)
                com_grid[ind] = CS.add_var!(com, 7, 11)
            else
                com_grid[ind] = CS.add_var!(com, 0, 8)
            end
        else
            com_grid[ind] = CS.add_var!(com, 0, 8; fix=val)
        end
    end

    add_sudoku_constr!(com, com_grid)

    @test CS.solve!(com, CS.get_default_options()) == :Solved
    @test fulfills_sudoku_constr(com_grid)
end

@testset "top95 some use backtracking" begin
    grids = sudokus_from_file("data/top95")
    c = 0
    for grid in grids
        m = Model(with_optimizer(CS.Optimizer))

        @variable(m, 1 <= x[1:9,1:9] <= 9, Int)
        # set variables
        for r=1:9, c=1:9
            if grid[r,c] != 0
                @constraint(m, x[r,c] == grid[r,c])
            end
        end

        # sudoku constraints
        jump_add_sudoku_constr!(m, x)

        optimize!(m)

        @test JuMP.termination_status(m) == MOI.OPTIMAL
        @test jump_fulfills_sudoku_constr(JuMP.value.(x))
        c += 1
    end
    # check that actually all 95 problems were tested
    @test c == 95
end

@testset "Number 7 in top95.txt w/o backtracking" begin
    m = Model(with_optimizer(CS.Optimizer, backtrack=false))

    grid = Int[6,0,2,0,5,0,0,0,0,0,0,0,0,0,3,0,4,0,0,0,0,0,0,0,0,0,0,4,3,0,0,0,8,0,
              0,0,0,1,0,0,0,0,2,0,0,0,0,0,0,0,0,7,0,0,5,0,0,2,7,0,0,0,0,0,0,0,0,0,
              0,0,8,1,0,0,0,6,0,0,0,0,0]
    grid = transpose(reshape(grid, (9,9)))

    @variable(m, 1 <= x[1:9,1:9] <= 9, Int)
    # set variables
    for r=1:9, c=1:9
        if grid[r,c] != 0
            @constraint(m, x[r,c] == grid[r,c])
        end
    end

    # sudoku constraints
    jump_add_sudoku_constr!(m, x)

    optimize!(m)

    com = JuMP.backend(m).optimizer.model.inner
    @test JuMP.termination_status(m) == MOI.OTHER_LIMIT

    @test !com.info.backtracked
    @test com.info.backtrack_fixes == 0
    @test com.info.in_backtrack_calls == 0
    @show com.info
end

end


text:
using SparseArrays
using LinearAlgebra
using MatrixDepot
using BenchmarkTools
using UnicodePlots
using PrettyTables
using ChainPartitioners

for mtx in [
            "Boeing/ct20stif",
            "DIMACS10/chesapeake",
            #"Schmid/thermal1",
            "Rothberg/3dtube",
           ]
    A = permutedims(1.0 * sparse(mdopen(mtx).A))

    mdl = BlockComponentCostModel{Int64}((8, 8), 0, 0, (1, identity), (4, identity))

    println()
    println()
    println(spy(A, maxwidth=50, maxheight=50, title="$mtx"))
    rows = []
    for (key, method) in [
        ("original", nothing),
        ("1D", AlternatingPacker(EquiChunker(1), DynamicTotalChunker(mdl, 8))),
        ("strict", AlternatingPacker(StrictChunker(8), StrictChunker(8))),
        ("overlap", AlternatingPacker(OverlapChunker(0.9, 8), OverlapChunker(0.9, 8))),
        ("dynamic", AlternatingPacker(DynamicTotalChunker(AffineFillConnectivityModel(0, 0, 1, 4), 8), DynamicTotalChunker(mdl, 8))),
        ("dynamic2", AlternatingPacker(
            DynamicTotalChunker(AffineFillConnectivityModel(0, 0, 1, 4), 8),
            DynamicTotalChunker(mdl, 8),
            DynamicTotalChunker(mdl, 8),
            DynamicTotalChunker(mdl, 8),
            DynamicTotalChunker(mdl, 8),
        )),
        ("dynamic3", AlternatingPacker(
            DynamicTotalChunker(AffineFillConnectivityModel(0, 0, 1, 1), 8),
            DynamicTotalChunker(BlockComponentCostModel{Int64}((8, 8), 0, 0, (1, identity), (2, identity)), 8),
            DynamicTotalChunker(BlockComponentCostModel{Int64}((8, 8), 0, 0, (1, identity), (3, identity)), 8),
            DynamicTotalChunker(mdl, 8),
            DynamicTotalChunker(mdl, 8)
        )),
        ("dynamic4", AlternatingPacker(
            DynamicTotalChunker(AffineFillConnectivityModel(0, 0, 1, 7), 8),
            DynamicTotalChunker(BlockComponentCostModel{Int64}((8, 8), 0, 0, (1, identity), (6, identity)), 8),
            DynamicTotalChunker(BlockComponentCostModel{Int64}((8, 8), 0, 0, (1, identity), (5, identity)), 8),
            DynamicTotalChunker(mdl, 8),
            DynamicTotalChunker(mdl, 8)
        )),
        #=
        ("strict_sym", SymmetricPacker(StrictChunker(8))),
        ("overlap_sym", SymmetricPacker(OverlapChunker(0.9, 8))),
        ("dynamic_sym", SymmetricPacker(DynamicTotalChunker(AffineFillConnectivityModel(0, 0, 1, 4), 8), DynamicTotalChunker(mdl, 8))),
        ("dynamic2_sym", SymmetricPacker(
            DynamicTotalChunker(AffineFillConnectivityModel(0, 0, 1, 4), 8),
            DynamicTotalChunker(BlockComponentCostModel{Int64}((8, 8), 0, 0, (1, identity), (4, identity)), 8),
            DynamicTotalChunker(BlockComponentCostModel{Int64}((8, 8), 0, 0, (1, identity), (4, identity)), 8),
            DynamicTotalChunker(mdl, 8)
        )),
        ("dynamic3_sym", SymmetricPacker(
            DynamicTotalChunker(AffineFillConnectivityModel(0, 0, 1, 1), 8),
            DynamicTotalChunker(BlockComponentCostModel{Int64}((8, 8), 0, 0, (1, identity), (2, identity)), 8),
            DynamicTotalChunker(BlockComponentCostModel{Int64}((8, 8), 0, 0, (1, identity), (3, identity)), 8),
            DynamicTotalChunker(mdl, 8)
        )),
        =#
    ]
        if method === nothing
            setup_time = 0
            memory = 5 * nnz(A)
        else
            setup_time = time(@benchmark (Π, Φ) = pack_plaid($A, $method))
            (Π, Φ) = pack_plaid(A, method)
            memory = total_value(A, Π, Φ, mdl)
        end

        push!(rows, [key setup_time memory])
    end
    pretty_table(vcat(rows...), ["method", "setuptime", "memory"])
end


text:
"""
    apply_objects!(g, objects)

Internal function to apply a vector of `Object` contained in an `Axes` container in a GLE context.
"""
function apply_objects!(g::GLE, objects::Vector{<:Object}, figid::String)
    foreach(o -> apply_object!(g, o, figid), objects)
    return nothing
end

function apply_object!(g::GLE, obj::Text2D, ::String)
    "\ngsave"    |> g
    "\nset just $(obj.position)"                        |> g
    apply_textstyle!(g, obj.textstyle, addset=true)
    "\namove xg($(obj.anchor[1])) yg($(obj.anchor[2]))" |> g
    "\nwrite \"$(obj.text)\""                           |> g
    "\ngrestore" |> g
    return nothing
end

function apply_object!(g::GLE, obj::StraightLine2D, ::String)
    "\ngsave"    |> g
    apply_linestyle!(g, obj.linestyle; addset=true)
    if obj.horiz
        "\namove xg(xgmin) yg($(obj.anchor))" |> g
        "\naline xg(xgmax) yg($(obj.anchor))" |> g
    else
        "\namove xg($(obj.anchor)) yg(ygmin)" |> g
        "\naline xg($(obj.anchor)) yg(ygmax)" |> g
    end
    "\ngrestore" |> g
    return nothing
end

function apply_object!(g::GLE, obj::Box2D, ::String)
    "\ngsave" |> g
    obj.position == "bl" || "\nset just $(obj.position)"      |> g
    apply_linestyle!(g, obj.linestyle; addset=true)
    "\namove xg($(obj.anchor[1])) yg($(obj.anchor[2]))"       |> g
    "\nbox xg($(obj.size[1]))-xg(0) yg($(obj.size[2]))-yg(0)" |> g
    obj.nobox && "nobox"                    |> g
    fs = obj.fillstyle
    isdef(fs) && "fill $(col2str(fs.fill))" |> g
    "\ngrestore" |> g
    return nothing
end

function apply_object!(g::GLE, obj::Colorbar, figid::String)

    add_sub_palette!(Figure(figid; _noreset=true), obj.cmap)

    dx, dy = obj.offset

    width, height = 0.0, 0.0
    if isdef(obj.size)
        width, height = obj.size
    else
        if obj.position ∈ ["left", "right"]
            width, height = "0.25", "abs(yg(ygmax)-yg(ygmin))"
        else
            width, height = "abs(xg(xgmax)-xg(xgmin))", "0.25"
        end
    end

    if obj.position == "right"
        "\namove xg(xgmax)+$dx yg(ygmin)+$dy"             |> g
    elseif obj.position == "left"
        "\namove xg(xgmin)-$dx-0.3-$width yg(ygmin)+$dy"  |> g
    elseif obj.position == "bottom"
        "\namove xg(xgmin)+$dx yg(ygmin)-$dy-0.3-$height" |> g
    else
        "\namove xg(xgmin)+$dx yg(ygmax)+$dy"             |> g
    end

    "\nbegin box name cmap" |> g
    obj.nobox && "nobox"    |> g

    #
    # colormap "y" 0 1 0 1 1 pixels width height palette palette$
    #
    "\n\tcolormap" |> g
    if obj.position ∈ ["right", "left"]
        "\"y\" 0 1 0 1 1 $(obj.pixels)"     |> g
    else
        "\"x\" 0 1 0 1 $(obj.pixels) 1"     |> g
    end
    "$width $height palette cmap_$(hash(obj.cmap))" |> g
    "\nend box" |> g

    # ticks (not axis elem this time)
    # map the ticks from 0.0 to 1.0
    ticks = collect(obj.ticks.places)
    ticks .-= obj.zmin
    ticks ./= (obj.zmax - obj.zmin) # now ticks on [0.0, 1.0]

    if obj.position == "right"
        if !(obj.ticks.off)
            tlength = isdef(obj.ticks.length) ? obj.ticks.length : "$width/3"
            "\ngsave"    |> g
            apply_linestyle!(g, obj.ticks.linestyle; nosmooth=true, addset=true)
            for (i, tick) ∈ enumerate(ticks)
                "\namove xg(xgmax)+$dx+$width yg(ygmin)+$dy+$height*$tick" |> g
                "\nrline $tlength 0" |> g # draw tick
            end
            "\ngrestore" |> g
        end
        if !(obj.ticks.labels.off)
            offset = isdef(obj.ticks.labels.dist) ? obj.ticks.labels.dist : "$width/3"
            shift  = isdef(obj.ticks.labels.shift) ? obj.ticks.labels.shift : 0
            "\ngsave"       |> g
            apply_textstyle!(g, obj.ticks.labels.textstyle; addset=true)
            "\nset just lc" |> g # justify center wrt anchor
            # retrieve the labels
            labels = obj.ticks.labels.names
            if isempty(labels)
                labels = ["$(round(ticks[i]*(obj.zmax-obj.zmin)+obj.zmin, digits=1))"
                                for i ∈ 1:length(ticks)]
            end
            # write them
            for (i, tick) ∈ enumerate(ticks)
                "\namove xg(xgmax)+$dx+$width yg(ygmin)+$dy+$height*$tick" |> g
                "\nrmove $offset*1.3 0"    |> g # move a bit more to write the label
                iszero(shift) || "\nrmove 0 $shift" |> g # shift vertical
                "\nwrite $(labels[i])" |> g # write label
            end
            "\ngrestore" |> g
        end
    elseif obj.position == "left"
        if !(obj.ticks.off)
            tlength = isdef(obj.ticks.length) ? obj.ticks.length : "$width/3"
            "\ngsave"    |> g
            apply_linestyle!(g, obj.ticks.linestyle; nosmooth=true, addset=true)
            for (i, tick) ∈ enumerate(ticks)
                "\namove xg(xgmin)-$dx yg(ygmin)+$dy" |> g
                "\nrline -$tlength 0" |> g # draw tick
            end
            "\ngrestore" |> g
        end
        if !(obj.ticks.labels.off)
            offset = isdef(obj.ticks.labels.dist) ? obj.ticks.labels.dist : "$width/3"
            shift  = isdef(obj.ticks.labels.shift) ? obj.ticks.labels.shift : 0
            "\ngsave"       |> g
            "\nset just lc" |> g # justify center wrt anchor
            apply_textstyle!(g, obj.ticks.labels.textstyle; addset=true)
            # retrieve the labels
            labels = obj.ticks.labels.names
            if isempty(labels)
                labels = ["$(round(obj.ticks.places[i], digits=1))" for i ∈ 1:length(ticks)]
            end
            # write them
            for (i, tick) ∈ enumerate(ticks)
                "\namove xg(xgmin)-$dx yg(ygmin)+$dy" |> g
                "\nrmove -$offset*1.3 0"    |> g # move a bit more to write the label
                iszero(shift) || "\nrmove 0 $shift"   |> g # shift vertical
                "\nset just lc"         |> g # justify center wrt anchor
                "\nwrite $(labels[i])"  |> g # write label
            end
            "\ngrestore" |> g
        end
    elseif obj.position == "bottom"
        if !(obj.ticks.off)
            tlength = isdef(obj.ticks.length) ? obj.ticks.length : "$height/3"
            "\ngsave"    |> g
            apply_linestyle!(g, obj.ticks.linestyle; nosmooth=true, addset=true)
            for (i, tick) ∈ enumerate(ticks)
                "\namove xg(xgmin)+$dx yg(ygmin)-$dy-0.3-$height" |> g
                "\nrline 0 -$tlength" |> g # draw tick
            end
            "\ngrestore" |> g
        end
        if !(obj.ticks.labels.off)
            offset = isdef(obj.ticks.labels.dist) ? obj.ticks.labels.dist : "$width/3"
            shift  = isdef(obj.ticks.labels.shift) ? obj.ticks.labels.shift : 0
            "\ngsave"       |> g
            apply_textstyle!(g, obj.ticks.labels.textstyle; addset=true)
            "\nset just lc" |> g # justify center wrt anchor
            # retrieve the labels
            labels = obj.ticks.labels.names
            if isempty(labels)
                labels = ["$(round(obj.ticks.places[i], digits=1))" for i ∈ 1:length(ticks)]
            end
            # write them
            for (i, tick) ∈ enumerate(ticks)
                "\namove xg(xgmin)+$dx yg(ygmin)-$dy-0.3-$height" |> g
                "\nrmove 0 -$offset*1.3"    |> g # move a bit more to write the label
                iszero(shift) || "\nrmove $shift 0" |> g # shift horizontal
                "\nset just lc"         |> g # justify center wrt anchor
                "\nwrite $(labels[i])"  |> g # write label
            end
            "\ngrestore" |> g
        end
    else
        if !(obj.ticks.off)
            tlength = isdef(obj.ticks.length) ? obj.ticks.length : "$height/3"
            "\ngsave"    |> g
            apply_linestyle!(g, obj.ticks.linestyle; nosmooth=true, addset=true)
            for (i, tick) ∈ enumerate(ticks)
                "\namove xg(xgmin)+$dx yg(ygmax)+$dy" |> g
                "\nrline 0 $tlength" |> g # draw tick
            end
            "\ngrestore" |> g
        end
        if !(obj.ticks.labels.off)
            offset = isdef(obj.ticks.labels.dist) ? obj.ticks.labels.dist : "$width/3"
            shift  = isdef(obj.ticks.labels.shift) ? obj.ticks.labels.shift : 0
            "\ngsave"       |> g
            apply_textstyle!(g, obj.ticks.labels.textstyle; addset=true)
            "\nset just lc" |> g # justify center wrt anchor
            # retrieve the labels
            labels = obj.ticks.labels.names
            if isempty(labels)
                labels = ["$(round(obj.ticks.places[i], digits=1))" for i ∈ 1:length(ticks)]
            end
            # write them
            for (i, tick) ∈ enumerate(ticks)
                "\namove xg(xgmin)+$dx yg(ygmax)+$dy" |> g
                "\nrmove 0 $offset*1.3"     |> g # move a bit more to write the label
                iszero(shift) || "\nrmove $shift 0"   |> g # shift horizontal
                "\nset just lc"         |> g # justify center wrt anchor
                "\nwrite $(labels[i])"  |> g # write label
            end
            "\ngrestore" |> g
        end
    end
    return nothing
end