text:
import Pkg
import UUIDs

function maketempdir()::String
    path::String = mktempdir()
    atexit(() -> _force_remove_path(path))
    return path
end

function _force_remove_path(path)::Nothing
    try
        rm(path; force = true, recursive = true)
    catch
    end
    return nothing
end

original_directory = pwd()
project_root = joinpath(splitpath(@__DIR__)...)
cd(project_root)

if length(ARGS) < 3
    @error(
        string(
            "syntax: julia add_from_external_registry.jl ",
            "\"EXTERNAL_REGISTRY_NAME\" ",
            "\"EXTERNAL_REGISTRY_UUID\" ",
            "\"EXTERNAL_REGISTRY_URL\" ",
            "[list of packages]",
            )
        )
    @error(
        string(
            "example usage #1: julia add_from_external_registry.jl ",
            "\"General\" ",
            "\"23338594-aafe-5451-b93e-139f81909106\" ",
            "\"https://github.com/JuliaRegistries/General.git\" ",
            "\"Foo\"",
            )
        )
    @error(
        string(
            "example usage #2: julia add_from_external_registry.jl ",
            "\"General\" ",
            "\"23338594-aafe-5451-b93e-139f81909106\" ",
            "\"https://github.com/JuliaRegistries/General.git\" ",
            "\"Foo\" \"Bar\" \"Baz\"",
            )
        )
    error("Not enough command-line arguments")
end

external_registry = Pkg.RegistrySpec(
    name = strip(ARGS[1]),
    uuid = UUIDs.UUID(strip(ARGS[2])),
    url = strip(ARGS[3]),
    )

this_registry_toml_path= joinpath(project_root, "Registry.toml")
this_registry_configuration = Pkg.TOML.parsefile(this_registry_toml_path)
this_registry_name_to_path = Dict{String, String}()
this_registry_name_to_uuid = Dict{String, String}()
this_registry_all_packages = String[]
if !haskey(this_registry_configuration, "packages")
    this_registry_configuration["packages"] = Dict{String, Any}()
end
for pair in this_registry_configuration["packages"]
    uuid = pair[1]
    name = pair[2]["name"]
    path = pair[2]["path"]
    push!(this_registry_all_packages, name,)
    this_registry_name_to_path[name] = path
    this_registry_name_to_uuid[name] = uuid
end

unique!(this_registry_all_packages)
sort!(this_registry_all_packages)
n = length(this_registry_all_packages)
packages_to_add_from_external_registry = String[
    strip(s) for s in ARGS[4:end]
    ]
unique!(packages_to_add_from_external_registry)
sort!(packages_to_add_from_external_registry)
@debug(
    "packages_to_add_from_external_registry",
    length(packages_to_add_from_external_registry),
    packages_to_add_from_external_registry,
    repr(packages_to_add_from_external_registry),
    )

my_depot = joinpath(maketempdir(), "depot",)
my_environment = joinpath(maketempdir(), "depot",)
rm(my_depot; force = true, recursive = true,)
rm(my_environment;force = true,recursive = true,)
mkpath(my_depot)
mkpath(my_environment)
original_depot_path = [x for x in Base.DEPOT_PATH]
empty!(Base.DEPOT_PATH)
pushfirst!(Base.DEPOT_PATH, my_depot,)
Pkg.activate(my_environment)
Pkg.Registry.add(external_registry)
external_registry_afteradding = first(
    Pkg.Types.collect_registries()
    )
external_registry_configuration = Pkg.TOML.parsefile(
    joinpath(external_registry_afteradding.path,"Registry.toml",)
    )
external_registry_name_to_path = Dict{String, String}()
external_registry_name_to_uuid = Dict{String, String}()
for pair in external_registry_configuration["packages"]
    uuid = pair[1]
    name = pair[2]["name"]
    path = pair[2]["path"]
    external_registry_name_to_path[name] = path
    external_registry_name_to_uuid[name] = uuid
end
external_registry_all_packages = Set(
    collect(keys(external_registry_name_to_path))
    )

info_messages = Vector{String}(undef, 0)
n = length(packages_to_add_from_external_registry)
for i = 1:n
    name = packages_to_add_from_external_registry[i]
    @debug("Adding \"$(name)\" (package $(i) of $(n))")
    if name in external_registry_all_packages
        if name in this_registry_all_packages
            error(
                string(
                    "We already have the package \"$(name)\" ",
                    "in our registry. ",
                    "Please delete the \"$(name)\" package ",
                    "and then try again.",
                    )
                )
        else
            their_uuid = external_registry_name_to_uuid[name]
            their_path = external_registry_name_to_path[name]
            their_full_path = joinpath(
                external_registry_afteradding.path,
                their_path,
                )
            our_path = joinpath(
                "packages",
                "$(uppercase(name[1:1]))",
                "$(name)",
                )
            our_full_path = joinpath(project_root, our_path)
            if isdir(our_full_path) || ispath(our_full_path)
                error(
                    string(
                        "The directory $(our_path) already exists at ",
                        "$(our_full_path). ",
                        "Please delete that directory and try again.",
                        )
                    )
            else
                rm(
                    our_full_path;
                    force = true,
                    recursive = true,
                    )
                mkpath(our_full_path)
                our_compat_toml =joinpath(our_full_path, "Compat.toml")
                our_deps_toml=joinpath(our_full_path, "Deps.toml")
                our_package_toml=joinpath(our_full_path, "Package.toml")
                our_versions_toml=joinpath(our_full_path, "Versions.toml")
                their_compat_toml=joinpath(their_full_path, "Compat.toml")
                their_deps_toml=joinpath(their_full_path, "Deps.toml")
                their_package_toml=joinpath(their_full_path, "Package.toml")
                their_versions_toml=joinpath(their_full_path, "Versions.toml")
                if isfile(their_compat_toml)
                    cp(
                        their_compat_toml,
                        our_compat_toml;
                        force = true,
                        )
                end
                if isfile(their_deps_toml)
                    cp(
                        their_deps_toml,
                        our_deps_toml;
                        force = true,
                        )
                end
                if isfile(their_package_toml)
                    cp(
                        their_package_toml,
                        our_package_toml;
                        force = true,
                        )
                end
                if isfile(their_versions_toml)
                    cp(
                        their_versions_toml,
                        our_versions_toml;
                        force = true,
                        )
                end
                this_registry_configuration["packages"][their_uuid] = Dict(
                    "name" => name,
                    "path" => our_path,
                    )
                push!(info_messages, "Added $(name)")
            end
        end
    else
        error(
            string(
                "The external registry does not have the package ",
                "\"$(name)\" ",
                "in their registry.",
                )
            )
    end
end

empty!(Base.DEPOT_PATH)
for x in original_depot_path
    push!(Base.DEPOT_PATH, x,)
end
unique!(Base.DEPOT_PATH)

rm(this_registry_toml_path; force = true, recursive = true)
open(this_registry_toml_path, "w") do f
    Pkg.TOML.print(f, this_registry_configuration)
end
temp_registry_toml = read(this_registry_toml_path, String)
temp_registry_toml_parsed = Pkg.TOML.parse(temp_registry_toml)
if !haskey(temp_registry_toml_parsed, "packages")
    temp_output = string(temp_registry_toml,
                         "\n\n\n\n",
                         "[packages]",
                         "\n\n\n\n")
    rm(this_registry_toml_path; force = true, recursive = true)
    write(this_registry_toml_path, temp_output)
end
push!(info_messages, "Wrote Registry.toml file to $(this_registry_toml_path)")

for message in info_messages
    @info(message)
end

cd(original_directory)


text:
# File: riskreturn.jl
# Author:
# Date: 29-June-2020

# import XL data
[names,~,returns] = getXLData('hfma_matlab_data','10 Hedge Funds');

am1 = fMean(returns,12); # annual mean
am2 = fStd(returns,0,12); # annual std. dev.

# settings for scatter plot
scatter(am2,am1,'b*')
xlabel('Risk (# )');
ylabel('Return (# )');

ext: jl


text:
"Create a new field definition."
createfielddefn(name::AbstractString, etype::GDAL.OGRFieldType) =
    GDAL.fld_create(name, etype)

"Destroy a field definition."
destroy(fielddefn::FieldDefn) = GDAL.destroy(fielddefn)

"Set the name of this field."
setname(fielddefn::FieldDefn, name::AbstractString) =
    GDAL.setname(fielddefn, name)

"Fetch name of this field."
getname(fielddefn::FieldDefn) = GDAL.getnameref(fielddefn)

"Fetch type of this field."
gettype(fielddefn::FieldDefn) = GDAL.gettype(fielddefn)

"Set the type of this field."
settype(fielddefn::FieldDefn, etype::GDAL.OGRFieldType) =
    GDAL.settype(fielddefn, etype)

# """
#     OGR_Fld_GetSubType(OGRFieldDefnH hDefn) -> OGRFieldSubType
# Fetch subtype of this field.
# ### Parameters
# * **hDefn**: handle to the field definition to get subtype from.
# ### Returns
# field subtype.
# """
# function getsubtype(arg1::Ptr{OGRFieldDefnH})
#     ccall((:OGR_Fld_GetSubType,libgdal),OGRFieldSubType,(Ptr{OGRFieldDefnH},),arg1)
# end

# """
#     OGR_Fld_SetSubType(OGRFieldDefnH hDefn,
#                        OGRFieldSubType eSubType) -> void
# Set the subtype of this field.
# ### Parameters
# * **hDefn**: handle to the field definition to set type to.
# * **eSubType**: the new field subtype.
# """
# function setsubtype(arg1::Ptr{OGRFieldDefnH},arg2::OGRFieldSubType)
#     ccall((:OGR_Fld_SetSubType,libgdal),Void,(Ptr{OGRFieldDefnH},OGRFieldSubType),arg1,arg2)
# end

"Get the justification for this field."
getjustify(fielddefn::FieldDefn) = GDAL.getjustify(fielddefn)

"Set the justification for this field."
setjustify(fielddefn::FieldDefn, ejustify::GDAL.OGRJustification) =
    GDAL.setjustify(fielddefn, ejustify)

"""Get the formatting width for this field.

### Returns
the width, zero means no specified width.
"""
getwidth(fielddefn::FieldDefn) = GDAL.getwidth(fielddefn)

"Set the formatting width for this field in characters."
setwidth(fielddefn::FieldDefn, width::Integer) =
    GDAL.setwidth(fielddefn, width)

"Get the formatting precision for this field."
getprecision(fielddefn::FieldDefn) = GDAL.getprecision(fielddefn)

"Set the formatting precision for this field in characters."
setprecision(fielddefn::FieldDefn, precision::Integer) =
    GDAL.setprecision(fielddefn, precision)

# """
#     OGR_Fld_Set(OGRFieldDefnH hDefn,
#                 const char * pszNameIn,
#                 OGRFieldType eTypeIn,
#                 int nWidthIn,
#                 int nPrecisionIn,
#                 OGRJustification eJustifyIn) -> void
# Set defining parameters for a field in one call.
# ### Parameters
# * **hDefn**: handle to the field definition to set to.
# * **pszNameIn**: the new name to assign.
# * **eTypeIn**: the new type (one of the OFT values like OFTInteger).
# * **nWidthIn**: the preferred formatting width. Defaults to zero indicating undefined.
# * **nPrecisionIn**: number of decimals places for formatting, defaults to zero indicating undefined.
# * **eJustifyIn**: the formatting justification (OJLeft or OJRight), defaults to OJUndefined.
# """
# function set(arg1::Ptr{OGRFieldDefnH},arg2,arg3::OGRFieldType,arg4::Integer,arg5::Integer,arg6::OGRJustification)
#     ccall((:OGR_Fld_Set,libgdal),Void,(Ptr{OGRFieldDefnH},Cstring,OGRFieldType,Cint,Cint,OGRJustification),arg1,arg2,arg3,arg4,arg5,arg6)
# end

"Return whether this field should be omitted when fetching features."
isignored(fielddefn::FieldDefn) = Bool(GDAL.isignored(fielddefn))

"Set whether this field should be omitted when fetching features."
setignored(fielddefn::FieldDefn, ignore::Bool) =
    GDAL.setignored(fielddefn, ignore)

"Return whether this field can receive null values."
isnullable(fielddefn::FieldDefn) = Bool(GDAL.isnullable(fielddefn))

"Set whether this field can receive null values."
setnullable(fielddefn::FieldDefn, nullable::Bool) =
    GDAL.setnullable(fielddefn, nullable)

"Get default field value."
getdefault(fielddefn::FieldDefn) = GDAL.getdefault(fielddefn)

"Set default field value."
setdefault(fielddefn::FieldDefn, default) =
    GDAL.setdefault(fielddefn, default)

"Returns whether the default value is driver specific."
isdefaultdriverspecific(fielddefn::FieldDefn) =
    Bool(GDAL.isdefaultdriverspecific(fielddefn))

"Create a new field geometry definition."
creategeomfielddefn(name::AbstractString, etype::GDAL.OGRwkbGeometryType) =
    GDAL.gfld_create(name, etype)

"Destroy a geometry field definition."
destroy(gfd::GeomFieldDefn) = GDAL.destroy(gfd)

"Set the name of this field."
setname(gfd::GeomFieldDefn, name::AbstractString) =
    GDAL.setname(gfd, name)

"Fetch name of this field."
getname(gfd::GeomFieldDefn) = GDAL.getnameref(gfd)

"Fetch geometry type of this field."
gettype(gfd::GeomFieldDefn) = GDAL.gettype(gfd)

"Set the geometry type of this field."
settype(gfd::GeomFieldDefn, etype::GDAL.OGRwkbGeometryType) =
    GDAL.settype(gfd, etype)

"Fetch spatial reference system of this field."
getspatialref(gfd::GeomFieldDefn) = SpatialRef(GDAL.getspatialref(gfd))

"Set the spatial reference of this field."
setspatialref(gfd::GeomFieldDefn, spatialref::SpatialRef) =
    GDAL.setspatialref(gfd, spatialref)

"Return whether this geometry field can receive null values."
isnullable(gfd::GeomFieldDefn) = Bool(GDAL.isnullable(gfd))

"Set whether this geometry field can receive null values."
setnullable(gfd::GeomFieldDefn, nullable::Bool) =
    GDAL.setnullable(gfd, nullable)

"Return whether this field should be omitted when fetching features."
isignored(gfd::GeomFieldDefn) = Bool(GDAL.isignored(gfd))

"Set whether this field should be omitted when fetching features."
setignored(gfd::GeomFieldDefn, ignore::Bool) =
    GDAL.setignored(gfd, ignore)


text:
@testset "Tuple" begin

s = """
    (a,b) = (1,2,3)
"""
msgs = lintstr(s)
@test msgs[1].code == :W546
@test occursin("implicitly discarding values, 2 of 3 used", msgs[1].message)

s = """
    a, = (1,2,3)
"""
msgs = lintstr(s)
@test msgs[1].code == :W546
@test occursin("implicitly discarding values, 1 of 3 used", msgs[1].message)

@test isempty(lintstr("a = (1, 2, 3)"))

s = """
function f()
    (a,b,c) = (1,2,3)
    return (b,c,a)
end
"""
msgs = lintstr(s)
@test isempty(msgs)

s = """
    (a,b,c) = (1,2)
"""
msgs = lintstr(s)
@test msgs[1].code == :E418
@test occursin("RHS is a tuple, 3 of 2 variables used", msgs[1].message)

end


text:
@testset "Univariate leading coefficient" begin
    R, (x, y, z) = Nemo.PolynomialRing(Nemo.QQ, ["x", "y", "z"])
    p = x^2 * y + x^2 * (z + z^3) + y - 5
    @test lc_univariate(p, x) == y + z + z^3
    @test lc_univariate(p, z) == x^2
    @test lc_univariate(p, y) == 1 + x^2
    @test lc_univariate(x + y, z) == x + y
end


text:
#=
    The file contains rational reconstrction and CRT functions
=#

#------------------------------------------------------------------------------

# Rational number reconstruction implementation borrowed from CLUE
# and modified a bit to suit the 'Modern Computer Algebra' definitions
# Returns a rational r // h of QQ field in a canonical form such that
#   r // h ≡ a (mod m)
#
# let n = max( λ(a), λ(m) ) , where λ(x) is a number of bits for x
# O(n^2)
function rational_reconstruction(a::I, m::I) where {I<:Union{Int, BigInt}}
    a = mod(a, m)
    if a == 0 || m == 0
        return QQ(0, 1)
    end
    if m < 0
        m = -m
    end
    if a < 0
        a = m - a
    end
    if a == 1
        return QQ(1, 1)
    end
    bnd = sqrt(float(m) / 2)

    @debug "" bnd
    U = (I(1), I(0), m)
    V = (I(0), I(1), a)
    while abs(V[3]) >= bnd
        q = div(U[3], V[3])
        T = U .- q .* V
        U = V
        V = T
    end

    t = abs(V[2])
    r = V[3] * sign(V[2])

    @debug "" r t

    # changed from `<= bnd` to `<= m / bnd`
    # we can speed up this !
    # if t <= bnd && gcd(r, t) == 1
    #    return QQ(r, t)
    # end

    return QQ(r, t)

    # throw(DomainError(
    #    :($a//$m), "rational reconstruction of $a (mod $m) does not exist"
    # ))

    # return QQ(0, 1)
end

# returns the bound for rational reconstruction
# based on the current modulo
#
# The bound for rational reconstrction:
# as soon as the numerator in rational reconstruction
# exceeds this bound, the gcd iteration is stopped
function rational_reconstruction_bound(modulo::BigInt)
    setprecision(2*Base.GMP.MPZ.sizeinbase(modulo, 2)) do
        ceil(BigInt, sqrt(BigFloat(modulo) / 2))
    end
end

# 0 allocations!
#=
    Computes the rational reconstruction of `a` mod `m`.
    Namely, a pair of numbers num, den , such that
        num//den ≡ a (mod m)

    Additional params:
        bnd  = stores stopping criterion threshold
        buf, buf1, buf2, buf3  = buffers
        u1, u2, u3 = buffers
        v1, v2, v3 = buffers
=#
function rational_reconstruction!(
            num::BigInt, den::BigInt, bnd::BigInt, buf::BigInt,
            buf1::BigInt, buf2::BigInt, buf3::BigInt,
            u1::BigInt, u2::BigInt, u3::BigInt,
            v1::BigInt, v2::BigInt, v3::BigInt,
            a::BigInt, m::BigInt)

    if Base.GMP.MPZ.cmp_ui(a, 0) == 0
        Base.GMP.MPZ.set_ui!(num, 0)
        Base.GMP.MPZ.set_ui!(den, 1)
        return true
    end

    # assumes input is nonnegative
    @assert Base.GMP.MPZ.cmp_ui(a, 0) > 0

    if Base.GMP.MPZ.cmp_ui(a, 1) == 0
        Base.GMP.MPZ.set_ui!(num, 1)
        Base.GMP.MPZ.set_ui!(den, 1)
        return true
    end

    Base.GMP.MPZ.set_ui!(u1, 1)
    Base.GMP.MPZ.set_ui!(u2, 0)
    Base.GMP.MPZ.set!(u3, m)
    Base.GMP.MPZ.set_ui!(v1, 0)
    Base.GMP.MPZ.set_ui!(v2, 1)
    Base.GMP.MPZ.set!(v3, a)

    while true
        if Base.GMP.MPZ.cmp(v2, bnd) > 0
            # @error "" v2 bnd
            return false
        end

        Base.GMP.MPZ.set!(buf, v3)
        if Base.GMP.MPZ.cmp_ui(buf, 0) < 0
            Base.GMP.MPZ.neg!(buf)
        end

        if Base.GMP.MPZ.cmp(buf, bnd) < 0
            break
        end

        Base.GMP.MPZ.tdiv_q!(buf, u3, v3)

        Base.GMP.MPZ.mul!(buf1, buf, v1)
        Base.GMP.MPZ.mul!(buf2, buf, v2)
        Base.GMP.MPZ.mul!(buf3, buf, v3)

        Base.GMP.MPZ.sub!(buf1, u1, buf1)
        Base.GMP.MPZ.sub!(buf2, u2, buf2)
        Base.GMP.MPZ.sub!(buf3, u3, buf3)

        Base.GMP.MPZ.set!(u1, v1)
        Base.GMP.MPZ.set!(u2, v2)
        Base.GMP.MPZ.set!(u3, v3)

        Base.GMP.MPZ.set!(v1, buf1)
        Base.GMP.MPZ.set!(v2, buf2)
        Base.GMP.MPZ.set!(v3, buf3)
    end

    Base.GMP.MPZ.set!(den, v2)
    Base.GMP.MPZ.set!(num, v3)

    #=
    num and den are obviously coprime,
    gcd is not nessesary

    Base.GMP.MPZ.gcd!(buf, den, num)
    Base.GMP.MPZ.tdiv_q!(den, buf)
    Base.GMP.MPZ.tdiv_q!(num, buf)
    =#

    if Base.GMP.MPZ.cmp_ui(den, 0) < 0
        Base.GMP.MPZ.neg!(den)
        Base.GMP.MPZ.neg!(num)
    end

    true
end

#------------------------------------------------------------------------------

# 0 allocations!
#=
    Computes the unique x s.t
        x ≡ a1 mod m1
        x ≡ a2 mod m2

    Additional params:
        M   = m1*m2
        buf = buffer
        n1  = buffer
        n2  = buffer
=#
function CRT!(
            M::BigInt, buf::BigInt, n1::BigInt, n2::BigInt,
            a1::BigInt, minv1::BigInt, a2::UInt, minv2::BigInt,
            m1::BigInt, m2::BigInt)

    Base.GMP.MPZ.mul!(buf, m1, minv1)
    Base.GMP.MPZ.mul_ui!(n1, buf, a2)

    Base.GMP.MPZ.mul!(buf, m2, minv2)
    Base.GMP.MPZ.mul!(n2, buf, a1)

    Base.GMP.MPZ.add!(buf, n1, n2)
    Base.GMP.MPZ.fdiv_r!(buf, M)
end

#------------------------------------------------------------------------------


text:
import Dates
import Pkg
import Tar

function my_exit(process::Base.Process)
    wait(process)

    @info(
        "",
        process.exitcode,
        process.termsignal,
    )

    # Pass the exit code back up
    if process.termsignal != 0
        ccall(:raise, Cvoid, (Cint,), process.termsignal)

        # If for some reason the signal did not cause an exit, we'll exit manually.
        # We need to make sure that we exit with a non-zero exit code.
        if process.exitcode != 0
            exit(process.exitcode)
        else
            exit(1)
        end
    end
    exit(process.exitcode)
end

function get_from_env(name::AbstractString)
    value = ENV[name]
    result = convert(String, strip(value))::String
    return result
end

cleanup_string(str::AbstractString) = replace(str, r"[^A-Za-z0-9_]" => "_")

if Base.VERSION < v"1.6"
    throw(ErrorException("The `$(basename(@__FILE__))` script requires Julia 1.6 or greater"))
end

if length(ARGS) < 1
    throw(ErrorException("Usage: julia $(basename(@__FILE__)) [command...]"))
end

const build_number                      = get_from_env("BUILDKITE_BUILD_NUMBER") |> cleanup_string
const commit_full                       = get_from_env("BUILDKITE_COMMIT")       |> cleanup_string
const job_name                          = get_from_env("BUILDKITE_STEP_KEY")     |> cleanup_string
const buildkite_timeout_minutes_string  = get_from_env("BUILDKITE_TIMEOUT")

const commit_short = first(commit_full, 10)
const buildkite_timeout_minutes = parse(Int, buildkite_timeout_minutes_string)::Int
const cleanup_minutes = 15
const ctest_timeout_minutes = buildkite_timeout_minutes - cleanup_minutes
if ctest_timeout_minutes < 1
    msg = "ctest_timeout_minutes must be strictly positive"
    @error(
        msg,
        ctest_timeout_minutes,
        buildkite_timeout_minutes,
        cleanup_minutes,
    )
    throw(ErrorException(msg))
end

@info(
    "",
    build_number,
    job_name,
    commit_full,
    commit_short,
    ctest_timeout_minutes,
    buildkite_timeout_minutes,
    cleanup_minutes,
)

const my_archives_dir    = joinpath(pwd(), "my_archives_dir")
const my_temp_parent_dir = joinpath(pwd(), "my_temp_parent_dir")

mkpath(my_archives_dir)
mkpath(my_temp_parent_dir)

const TMPDIR = mktempdir(my_temp_parent_dir)

proc = nothing

mktempdir(my_temp_parent_dir) do dir
    Pkg.activate(dir)
    Pkg.add("Zstd_jll")
    zstd_jll = Base.require(Base.PkgId(Base.UUID("3161d3a3-bdf6-5164-811a-617609db77b4"), "Zstd_jll"))
    # zstdmt(func) = Base.invokelatest(zstd_jll.zstdmt, func; adjust_LIBPATH=false)
    zstdmt(func) = Base.invokelatest(zstd_jll.zstdmt, func)

    new_env = copy(ENV)
    new_env["TMPDIR"] = TMPDIR
    command = setenv(`$ARGS`, new_env)
    global proc = run(command, (stdin, stdout, stderr); wait = false)

    # Start asynchronous timer that will kill `ctest`
    @async begin
        sleep(ctest_timeout_minutes * 60)

        # If we've exceeded the timeout and `ctest` is still running, kill it
        if isopen(proc)
            @error(
                string(
                    "Process timed out ",
                    "(with a timeout of $(ctest_timeout_minutes) minutes). ",
                    "Killing with SIGTERM.",

                )
            )
            kill(proc, Base.SIGTERM)
        end
    end

    # Wait for `ctest` to finish, either through naturally finishing its run, or `SIGTERM`
    wait(proc)

    if proc.termsignal != 0
        @info "Command signalled $(proc.termsignal)"
    else
        @info "Command returned $(proc.exitcode)"
    end
    date_str = Dates.format(Dates.now(), Dates.dateformat"yyyy_mm_dd_HH_MM_SS")
    artifact_specifications = [
        ("TMPDIR", TMPDIR),
    ]
    for (artifact_name, artifact_input_dir) in artifact_specifications
        dst_file_name = string(
            artifact_name,
            "--build_$(build_number)",
            "--$(job_name)",
            "--commit_$(commit_short)",
            "--$(date_str)",
            ".tar.zst",
        )
        dst_full_path = joinpath(my_archives_dir, dst_file_name)
        run(`find . -type p`) # list the named pipes before we delete them
        run(`find . -type p -delete`)
        run(`find . -type s`) # list the sockets before we delete them
        run(`find . -type s -delete`)
        zstdmt() do zstdp
            tarproc = open(`$(zstdp) -o $(dst_full_path)`, "w")
            Tar.create(artifact_input_dir, tarproc)
            close(tarproc.in)
        end

        buildkite_upload_cmd = `buildkite-agent artifact upload $(dst_file_name)`
        if !success(proc)
            run(setenv(buildkite_upload_cmd; dir = my_archives_dir))
        end
    end
end

my_exit(proc)

max_stars_repo_name: mozilla/rr

max_issues_repo_name: rocallahan/rr

max_forks_repo_name: rocallahan/rr


text:
Polish = Command(
   "polish",
   "pangraph polish <options> [pangraph.json]",
   "realigns pancontigs of multiple sequence alignment graph",
   """zero or one pangraph file (native json)
      if no file, reads from stdin
      stream can be optionally gzipped.""",
   [
    Arg(
        Int,
        "maximum length",
        (short="-l", long="--length"),
        "cutoff above which we won't realign",
        typemax(Int),
    ),
    Arg(
        Bool,
        "preserve case",
        (short="-c", long="--preserve-case"),
        "ensure case (upper/lower) is preserved after realignment",
        false,
    ),
   ],
   function(args)
       path = parse(Polish, args)
       path = if (path === nothing || length(path) == 0)
           nothing
       elseif length(path) == 1
           path
       else
           return 2
       end

       graph = load(path, Polish)
       if !Shell.havecommand("mafft")
           panic("external command mafft not found. please install before running polish step\n")
       end

       case = arg(Polish, "-c")
       accept = function(blk)
           length(blk) ≤ arg(Polish, "-l") && Graphs.depth(blk) > 1
       end
       Graphs.realign!(graph; accept=accept, case=case)

       marshal(stdout, graph; fmt=:json)
       return 0
   end
)

max_stars_repo_name: neherlab/pangraph

max_issues_repo_name: neherlab/pangraph

max_forks_repo_name: neherlab/pangraph


text:
using Pkg
Pkg.activate(joinpath(@__DIR__, ".."))
Pkg.instantiate()

# ## setup
using Dojo
using IterativeLQR
using LinearAlgebra

# ## system
timestep = 0.05
gravity=-9.81
env = get_environment(:raiberthopper,
    representation=:minimal,
    timestep=timestep,
    gravity=gravity);

# ## visualizer
open(env.vis)

# ## dimensions
n = env.num_states
m = env.num_inputs

# ## states
z1 = maximal_to_minimal(env.mechanism, Dojo.raiberthopper_nominal_max())
zM = maximal_to_minimal(env.mechanism, Dojo.raiberthopper_offset_max(0.5, 0.5, 0.5))
zT = maximal_to_minimal(env.mechanism, Dojo.raiberthopper_offset_max(0.5, 0.5, 0.0))

# ## nominal control
u_control = [0.0; 0.0; env.mechanism.bodies[1].mass * env.mechanism.gravity[3] * env.mechanism.timestep]

# ## horizon
T = 21
Tm = convert(Int, floor((T - 1) / 2))

# ## model
dyn = IterativeLQR.Dynamics(
    (y, x, u, w) -> dynamics(y, env, x, u, w),
    (dx, x, u, w) -> dynamics_jacobian_state(dx, env, x, u, w),
    (du, x, u, w) -> dynamics_jacobian_input(du, env, x, u, w),
    n, n, m)

model = [dyn for t = 1:T-1]

# ## rollout
ū = [[0.0; 0.0; env.mechanism.bodies[1].mass * env.mechanism.gravity[3] * env.mechanism.timestep + 0.0 * randn(1)[1]] for t = 1:T-1]
x̄ = IterativeLQR.rollout(model, z1, ū)
visualize(env, x̄)

# ## objective
ot1 = (x, u, w) -> 1 * (transpose(x - zM) * Diagonal([1.0 * ones(3); 0.01 * ones(3); 0.1 * ones(3); 0.01 * ones(3); 1.0; 0.01]) * (x - zM) + transpose(u) * Diagonal(1.0e-2 * [1.0; 1.0; 1.0]) * u)
ot2 = (x, u, w) -> 1 * (transpose(x - zT) * Diagonal([1.0 * ones(3); 0.01 * ones(3); 0.1 * ones(3); 0.01 * ones(3); 1.0; 0.01]) * (x - zT) + transpose(u) * Diagonal(1.0e-2 * [1.0; 1.0; 1.0]) * u)
oT = (x, u, w) -> transpose(x - zT) * Diagonal([1.0 * ones(3); 0.01 * ones(3); 0.1 * ones(3); 0.01 * ones(3); 1.0; 0.01]) * (x - zT)

ct1 = IterativeLQR.Cost(ot1, n, m)
ct2 = IterativeLQR.Cost(ot2, n, m)
cT = IterativeLQR.Cost(oT, n, 0)
obj = [[ct1 for t = 1:Tm]..., [ct2 for t = 1:Tm]..., cT]

# ## constraints
function goal(x, u, w)
    Δ = x - zT
    return [Δ[collect(1:3)]; Δ[collect(7:9)]; Δ[13]]
end

cont = IterativeLQR.Constraint()
conT = IterativeLQR.Constraint(goal, n, 0)
cons = [[cont for t = 1:T-1]..., conT]

# ## solver
s = IterativeLQR.solver(model, obj, cons,
    opts=Options(
        linesearch=:armijo,
        α_min=1.0e-5,
        obj_tol=1.0e-3,
        grad_tol=1.0e-3,
        con_tol=0.005,
        max_iter=100,
        max_al_iter=10,
        ρ_init=1.0,
        ρ_scale=10.0,
        verbose=true))
IterativeLQR.initialize_controls!(s, ū)
IterativeLQR.initialize_states!(s, x̄)

# ## solve
@time IterativeLQR.solve!(s)

# ## solution
x_sol, u_sol = IterativeLQR.get_trajectory(s)
@show IterativeLQR.eval_obj(s.m_data.obj.costs, s.m_data.x, s.m_data.u, s.m_data.w)
@show s.s_data.iter[1]
@show norm(goal(s.m_data.x[T], zeros(0), zeros(0)), Inf)

# ## visualize
visualize(env, x_sol)


text:
# We wish to extend operations to identically named methods dispatched
# on `Machine`s and `NodalMachine`s. For example, we have from the model API
#
# `predict(model::M, fitresult, X) where M<:Supervised`
#
# but want also want
#
# `predict(machine::Machine, X)` where `X` is data
#
# and "networks.jl" requires us to define
#
# `predict(machine::NodalMachine, X)` where `X` is data
#
# and we would like the syntactic sugar (for `X` a node):
#
# `predict(machine::NodalMachine, X::Node)=node(predict, machine, X)`
#
# (If an operation has zero arguments, we cannot achieve the last
# desire because of ambiguity with the preceding one.)
#
# The following macros are for this purpose.

## TODO: need to add checks on the arguments of
## predict(::AbstractMachine, ) and transform(::AbstractMachine, )
## Refactor using simpler non-macro code generation

macro extend_to_machines(operation)
    quote

        # with arguments specified:
        function $(esc(operation))(machine::AbstractMachine, args...)
            if isdefined(machine, :fitresult)
                return $(esc(operation))(machine.model,
                                         machine.fitresult,
                                         args...)
            else
                throw(error("$machine has not been trained."))
            end
        end
        $(esc(operation))(machine::Machine; rows=:) =
            $(esc(operation))(machine, selectrows(machine.args[1], rows))
        $(esc(operation))(machine::Machine, task::MLJTask) =
            $(esc(operation))(machine, task.X)

    end
end

macro sugar(operation)
    quote
        $(esc(operation))(machine::NodalMachine, args::AbstractNode...) =
            node($(esc(operation)), machine, args...)
    end
end

@extend_to_machines predict
@extend_to_machines predict_mode
@extend_to_machines predict_mean
@extend_to_machines predict_median
@extend_to_machines transform
@extend_to_machines inverse_transform
@extend_to_machines se

@sugar predict
@sugar predict_mode
@sugar predict_mean
@sugar predict_median
@sugar transform
@sugar inverse_transform
@sugar se

# the zero argument special cases:
function evaluate(machine::AbstractMachine)
    if isdefined(machine, :fitresult)
        return evaluate(machine.model, machine.fitresult)
    else
        throw(error("$machine has not been trained."))
    end
end
function fitted_params(machine::AbstractMachine)
    if isdefined(machine, :fitresult)
        return fitted_params(machine.model, machine.fitresult)
    else
        throw(error("$machine has not been trained."))
    end
end

getindex(n::Node{<:NodalMachine{<:Model}}, s::Symbol) = getproperty(n.machine.model, s)
setindex!(n::Node{<:NodalMachine{<:Model}}, v, s::Symbol) = setproperty!(n.machine.model, s, v)


text:
include("utils.jl")

binary = false

binarize_data(x) = round(x) <= 0 ? -1 : 1

binarize_state(x) = x <= 0 ? -1 : 1

##

mutable struct RBM

    visibles::Array{Float32}
    hiddens::Array{Float32}

    weights::Array{Float32}

    RBM(in_size, hidden_size) = new(

        [0 for _ in 1:in_size],
        [0 for _ in 1:hidden_size],

        randn(in_size, hidden_size),

    )

end

(rbm::RBM)(input) =
begin

    rbm.hiddens = tanh.(input * rbm.weights)

    binary ? rbm.hiddens = binarize_state.(rbm.hiddens) : ()

end

(rbm::RBM)() =
begin

    rbm.visibles = tanh.(rbm.hiddens * transpose(rbm.weights))

    binary ? rbm.visibles = binarize_state.(rbm.visibles) : ()

end

##

energy(rbm) = -(rbm.visibles * rbm.weights * rbm.hiddens')

propogate_until_convergence!(rbm, input; k=1_000) =
begin

    prev_visibles = nothing
    prev_hiddens = nothing

    ctr = 0

    while (rbm.visibles != prev_visibles || rbm.hiddens != prev_hiddens) && ctr < k

        prev_visibles = rbm.visibles
        prev_hiddens = rbm.hiddens

        ctr +=1

        rbm(input)
        rbm()

        input = rbm.visibles

    end

input
end

##

alternating_gibbs_grads!(rbm, input; converge=false, k=1) =
begin

    rbm(input)

    pos_hiddens = rbm.hiddens

    if converge

        rbm()

        propogate_until_convergence!(rbm, rbm.visibles)

    else

        for _ in 1:k

            rbm()

            rbm(rbm.visibles)

        end

    end

    grads = (transpose(input) * pos_hiddens) .- (transpose(rbm.visibles) * rbm.hiddens)

grads
end

batch_grads(rbm, batch) = threadpool(sum, args->alternating_gibbs_grads!(args...), [[deepcopy(rbm), input] for input in batch]) ./ length(batch)

update_weights!(rbm, grads, lr) =
begin

    rbm.weights += lr .* grads

end

##


text:
@testset "ROC" begin
    y = [  0   0   0   1   0   1   1   0] |> vec |> categorical
    s = [0.0 0.1 0.1 0.1 0.2 0.2 0.5 0.5] |> vec
    ŷ = UnivariateFinite([0, 1], s, augment=true, pool=y)

    fprs, tprs, ts = roc(ŷ, y)

    sk_fprs = [0. , 0.2, 0.4, 0.8, 1. ]
    sk_tprs = [0. , 0.33333333, 0.66666667, 1., 1.]

    @test fprs ≈ sk_fprs
    @test tprs ≈ sk_tprs
end


text:
using SuperResolution
using Test

@testset "evalSR.jl" begin
    @test true
    # #my_f(x,y)=2x+3y
    # @test my_f(2,1)==7
    # @test my_f(2,3)==13
    # # Write your tests here.
end

max_stars_repo_name: PGSch/SuperResolution

max_issues_repo_name: PGSch/SuperResolution

max_forks_repo_name: PGSch/SuperResolution


text:
#################################################################################################
# Functions
#################################################################################################
function get_ϕs_eff_rng(eff_type, config)

    # Pick max values
    if eff_type == "mml"
        ϕ1_max = [-1.0,-10.0,2.0]
        ϕ2_max = [+1.0,+10.0,2.0]
    elseif eff_type == "nme"
        ϕ1_max = [0.0,0.0,0.0]
        ϕ2_max = [0.0,0.0,5.0]
    elseif eff_type == "cmd"
        ϕ1_max = [+1.0,+10.0,2.0]
        ϕ2_max = [+1.0,+10.0,2.0]
    else
        CpelNano.print_log("Effect type $(eff_type) type not recognized.)")
    end

    # Loop over percentages
    i = 1
    eff_perc = 0.0
    eff_size_arr = Vector{Float64}()
    ϕ1_arr = Vector{Vector{Float64}}()
    ϕ2_arr = Vector{Vector{Float64}}()
    next_eff_size = minimum(EF_PERC_RNG)
    eff_size_check = fill(false, length(EF_PERC_RNG))
    while ! all(eff_size_check)

        # Scale vector
        if eff_type == "mml"
            ϕ1 = [eff_perc * ϕ1_max[1],eff_perc * ϕ1_max[2],ϕ1_max[3]]
            ϕ2 = [eff_perc * ϕ2_max[1],eff_perc * ϕ2_max[2],ϕ2_max[3]]
        elseif eff_type == "nme"
            ϕ1 = eff_perc * ϕ1_max
            ϕ2 = eff_perc * ϕ2_max
        else
            ϕ1 = ϕ1_max * ( 1.0 .+ eff_perc)
            ϕ2 = ϕ2_max / ( 1.0 .+ eff_perc)
        end

        # Baseline struct
        rs_g1 = create_baseline_struct(config)
        rs_g2 = create_baseline_struct(config)

        # Set vector
        rs_g1.ϕhat = ϕ1
        rs_g2.ϕhat = ϕ2

        # Compute stats
        proc_rep_rs!(rs_g1, config)
        proc_rep_rs!(rs_g2, config)

        # Non-empty analysis regions
        kp_in = [int != 0:0 for int in rs_g1.nls_rgs.cpg_ind]

        # Effect size
        if eff_type == "mml"
            eff_size = mean(abs.(rs_g1.mml[kp_in] - rs_g2.mml[kp_in]))
        elseif eff_type == "nme"
            eff_size = mean(abs.(rs_g1.nme[kp_in] - rs_g2.nme[kp_in]))
        elseif eff_type == "cmd"
            eff_size = mean(CpelNano.comp_cmd(rs_g1, rs_g2)[kp_in])
        end

        println("eff_size=$(eff_size)")

        # Record vector
        if eff_size >= next_eff_size

            # Record vector
            push!(ϕ1_arr, ϕ1)
            push!(ϕ2_arr, ϕ2)
            push!(eff_size_arr, eff_size)

            # Break
            i == length(EF_PERC_RNG) && break

            # Next effect size
            next_eff_size = EF_PERC_RNG[i + 1]

            # Increase counter
            i += 1

        end

        # Increase effect percentage
        eff_perc += 0.0005

    end

    # Return vector
    return eff_size_arr, ϕ1_arr, ϕ2_arr

end
@everywhere function create_baseline_struct(config::CpelNano.CpelNanoConfig)

    # Initialize
    rs = CpelNano.RegStruct()

    # Read in file
    rs_data = readdlm("$(aim_dir)/estimation_region_chr_chr1_42900254_42903250.txt")

    # Get values
    rs.chrst = 42900254
    rs.chrend = 42903250
    rs.cpg_pos = Int64.(rs_data[:,1])
    rs.ρn = Float64.(rs_data[:,2])
    rs.dn = Float64.(rs_data[:,3])

    # Get rest of values
    rs.N = length(rs.cpg_pos)
    rs.cpg_grps = CpelNano.get_grps_from_cgs(rs.cpg_pos, config.min_grp_dist)
    rs.Nl = map(grp -> length(grp.cpg_ind), rs.cpg_grps)
    rs.ρl = CpelNano.get_ρl(rs.ρn, rs.cpg_grps)
    rs.dl = CpelNano.get_dl(rs.cpg_pos, rs.cpg_grps)
    rs.L = length(rs.cpg_grps)

    # Get analysis region info
    CpelNano.get_nls_reg_info!(rs, config)

    # Return struct
    return rs

end
@everywhere function cp_rs_flds!(rs::CpelNano.RegStruct, rs_aux::CpelNano.RegStruct)

    # Fill fields
    rs_aux.L = rs.L
    rs_aux.Nl = rs.Nl
    rs_aux.ρl = rs.ρl
    rs_aux.dl = rs.dl
    rs_aux.N = rs.N
    rs_aux.ρn = rs.ρn
    rs_aux.dn = rs.dn
    rs_aux.cpg_pos = rs.cpg_pos
    rs_aux.nls_rgs = rs.nls_rgs

    # Return nothing
    return nothing

end
@everywhere function gen_pair_data(ϕ1::Vector{Float64}, ϕ2::Vector{Float64}, config::CpelNano.CpelNanoConfig)

    ## Set baseline struct
    rs = create_baseline_struct(config)
    αs1, βs1 = CpelNano.get_αβ_from_ϕ(ϕ1, rs)
    αs2, βs2 = CpelNano.get_αβ_from_ϕ(ϕ2, rs)

    ## Group 1

    # Generate calls
    ms_g1 = CpelNano.cpel_samp_ont(M, αs1, βs1, POBS, LEVEL_MEAN_M, LEVEL_STD_M, LEVEL_MEAN_U, LEVEL_STD_U)

    # Fill fields
    cp_rs_flds!(rs, ms_g1)

    ## Group 2

    # Generate data
    ms_g2 = CpelNano.cpel_samp_ont(M, αs2, βs2, POBS, LEVEL_MEAN_M, LEVEL_STD_M, LEVEL_MEAN_U, LEVEL_STD_U)

    # Fill fields
    cp_rs_flds!(rs, ms_g2)

    # Returns group data
    return ms_g1, ms_g2

end
@everywhere function proc_rep_rs!(rs::CpelNano.RegStruct, config::CpelNano.CpelNanoConfig)

    # Estimate parameter vector ϕ
    if length(rs.ϕhat) == 0
        CpelNano.get_ϕhat!(rs, config)
        length(rs.ϕhat) > 0 || return nothing
    end

    # Re-scale to per cpg site resolution
    CpelNano.rscle_grp_mod!(rs)

    # Get matrices
    CpelNano.get_rs_lgtrck_mats!(rs)

    # Partition function
    CpelNano.get_rs_logZ!(rs)

    # Get them all
    CpelNano.get_rs_log_gs!(rs)

    # Compute E[log g(X)]
    CpelNano.get_rs_exp_log_g1!(rs)
    CpelNano.get_rs_exp_log_g2!(rs)

    # Compute E[X] and E[XX]
    CpelNano.get_rs_exps!(rs)

    # Compute μ(X)
    CpelNano.comp_mml!(rs)

    # Compute h(X)
    CpelNano.comp_nme!(rs)

    # Set processed
    rs.proc = true

    # Return nothing
    return nothing

end
function get_pvals(test_out::CpelNano.RegStatTestStruct, eff_type::String)

    # Get
    if eff_type == "mml"
        pvals = [x[2] for x in test_out.tests.tmml_test]
    elseif eff_type == "nme"
        pvals = [x[2] for x in test_out.tests.tnme_test]
    elseif eff_type == "cmd"
        pvals = [x[2] for x in test_out.tests.tcmd_test]
    else
        CpelNano.print_log("Wrong effect type $(eff_type)")
    end

    # Return p-values
    return pvals

end
@everywhere function pmap_two_samp_null_stats(rs_ref::CpelNano.RegStruct, calls_s1::Vector{Vector{CpelNano.MethCallCpgGrp}}, calls_s2::Vector{Vector{CpelNano.MethCallCpgGrp}}, perm_ids::Vector{Int64}, config::CpelNano.CpelNanoConfig)::NTuple{3,Vector{Float64}}

    # Init output
    tmml = fill(NaN, rs_ref.nls_rgs.num)
    tnme = fill(NaN, rs_ref.nls_rgs.num)
    tcmd = fill(NaN, rs_ref.nls_rgs.num)

    # Create two aux structs
    aux_rs_s1 = CpelNano.RegStruct()
    aux_rs_s2 = CpelNano.RegStruct()
    cp_rs_flds!(rs_ref, aux_rs_s1)
    cp_rs_flds!(rs_ref, aux_rs_s2)
    aux_rs_s1.m = length(calls_s1)
    aux_rs_s2.m = length(calls_s2)

    # Assign calls
    aux_calls = vcat(calls_s1, calls_s2)
    aux_rs_s1.calls = aux_calls[perm_ids]
    deleteat!(aux_calls, perm_ids)
    aux_rs_s2.calls = aux_calls

    # Do estimation
    proc_rep_rs!(aux_rs_s1, config)
    proc_rep_rs!(aux_rs_s2, config)
    (aux_rs_s1.proc && aux_rs_s2.proc) || return tmml, tnme, tcmd

    # Compute stats
    tmml = abs.(aux_rs_s1.mml - aux_rs_s2.mml)
    tnme = abs.(aux_rs_s1.nme - aux_rs_s2.nme)
    tcmd = CpelNano.comp_cmd(aux_rs_s1, aux_rs_s2)

    # Return tuple
    return tmml, tnme, tcmd

end
function pmap_diff_two_samp_comp(mod_s1::CpelNano.RegStruct, mod_s2::CpelNano.RegStruct, config::CpelNano.CpelNanoConfig)::CpelNano.RegStatTestStruct

    # Init output struct
    test_struct = CpelNano.RegStatTestStruct(mod_s1)

    # Observed statistics
    tmml_obs = abs.(mod_s1.mml - mod_s2.mml)
    tnme_obs = abs.(mod_s1.nme - mod_s2.nme)
    tcmd_obs = CpelNano.comp_cmd(mod_s1, mod_s2)

    # Init p-values
    tmml_pvals = fill(NaN, length(tmml_obs))
    tnme_pvals = fill(NaN, length(tnme_obs))
    tcmd_pvals = fill(NaN, length(tcmd_obs))

    # Compute number of possible randomizations
    L = binomial(length(mod_s1.calls) + length(mod_s2.calls), length(mod_s1.calls))

    # If enough permutations
    if L > 20

        #  Check if exact test
        exact = L < LMAX_TWO_SAMP_AIM_5

        # Store calls in aux vars
        calls_s1 = deepcopy(mod_s1.calls)
        calls_s2 = deepcopy(mod_s2.calls)

        # Create iteratable object with all combinations
        comb_iter = combinations(1:(length(calls_s1) + length(calls_s2)), length(calls_s1))

        # Get sample label combinations to use
        comb_iter_used = []
        if exact
            # Use all sample assignments
            comb_iter_used = comb_iter
        else
            # Use Lmax group assignments
            ind_subset = rand(1:L, LMAX_TWO_SAMP_AIM_5)
            @inbounds for (ind, comb) in enumerate(comb_iter)
                (ind in ind_subset) && push!(comb_iter_used, comb)
            end
        end

        ## Null statistics

        # Compute null statistics
        pmap_out = pmap(perm -> pmap_two_samp_null_stats(mod_s1, calls_s1, calls_s2, perm, config), comb_iter_used)

        # Distribute statistics
        tmml_perms = [x[1] for x in pmap_out]
        tnme_perms = [x[2] for x in pmap_out]
        tcmd_perms = [x[3] for x in pmap_out]

        ## P-value computation

        # Loop over analysis regions
        @inbounds for k in 1:length(tmml_obs)

            # Check if data
            isnan(tmml_obs[k]) && continue

            # Get permutation stats from k-th analysis region
            tmml_perms_k = [perm[k] for perm in tmml_perms]
            tnme_perms_k = [perm[k] for perm in tnme_perms]
            tcmd_perms_k = [perm[k] for perm in tcmd_perms]

            # Clean NaNs in null stats
            tmml_perms_k = tmml_perms_k[.!isnan.(tmml_perms_k)]
            tnme_perms_k = tnme_perms_k[.!isnan.(tnme_perms_k)]
            tcmd_perms_k = tcmd_perms_k[.!isnan.(tcmd_perms_k)]

            # Check enough null stats after filtering NaNs
            # length(tmml_perms_k) > 20 || continue

            # Get number of permutation stats equal or above observed
            tmml_pval_k = sum(abs.(tmml_perms_k) .>= abs(tmml_obs[k]))
            tnme_pval_k = sum(abs.(tnme_perms_k) .>= abs(tnme_obs[k]))
            tcmd_pval_k = sum(tcmd_perms_k .>= tcmd_obs[k])

            # Compute p-values
            tmml_pvals[k] = exact ? tmml_pval_k / length(tmml_perms_k) : (1.0 + tmml_pval_k) / (1.0 + length(tmml_perms_k))
            tnme_pvals[k] = exact ? tnme_pval_k / length(tnme_perms_k) : (1.0 + tnme_pval_k) / (1.0 + length(tnme_perms_k))
            tcmd_pvals[k] = exact ? tcmd_pval_k / length(tcmd_perms_k) : (1.0 + tcmd_pval_k) / (1.0 + length(tcmd_perms_k))

        end

    end

    # Fill return object
    @inbounds for k in 1:length(tmml_obs)
        test_struct.tests.tmml_test[k] = (tmml_obs[k], tmml_pvals[k])
        test_struct.tests.tnme_test[k] = (tnme_obs[k], tnme_pvals[k])
        test_struct.tests.tcmd_test[k] = (tcmd_obs[k], tcmd_pvals[k])
    end

    # Return test struct
    return test_struct

end
@everywhere function map_run_sim(ϕ1::Vector{Float64}, ϕ2::Vector{Float64}, config::CpelNano.CpelNanoConfig)

    # Choose whether effect or not (50% chance)
    eff_bool_j = rand() > 0.5

    # Generate data
    if eff_bool_j
        # Introduce effect
        m_g1, m_g2 = gen_pair_data(ϕ1, ϕ2, config)
    else
    # Do not introduce effect
        m_g1, m_g2 = gen_pair_data(ϕ1, ϕ1, config)
    end

    # Estimate parameters
    proc_rep_rs!(m_g1, config)
    proc_rep_rs!(m_g2, config)

    # Find valid reps
    m_g1.proc && m_g2.proc || return nothing

    # Perform test
    test_out = pmap_diff_two_samp_comp(m_g1, m_g2, config)

    # Return test out and effect bool
    return eff_bool_j, test_out

end
function map_run_null_sim(ϕ::Vector{Float64}, config::CpelNano.CpelNanoConfig)

    CpelNano.print_log("Generating set of p-values...")

    # Generate data
    m_g1, m_g2 = gen_pair_data(ϕ, ϕ, config)

    # Estimate parameters
    proc_rep_rs!(m_g1, config)
    proc_rep_rs!(m_g2, config)

    # Find valid reps
    m_g1.proc && m_g2.proc || return nothing

    # Perform test
    test_out = pmap_diff_two_samp_comp(m_g1, m_g2, config)

    # Return test out
    return test_out

end
function run_sim(n_sim_reps::Int64, eff_type::String)

    # Init
    pval_thrs_rng = 0.0:0.001:1.0
    tpr = Array{Float64}(undef, length(EFF_PERC_RNG), length(pval_thrs_rng))
    fpr = Array{Float64}(undef, length(EFF_PERC_RNG), length(pval_thrs_rng))

    # CpelNano config
    config = CpelNano.CpelNanoConfig()
    config.max_size_subreg = 350
    config.max_em_iters = 25
    config.max_em_init = 5
    config.verbose = false

    # Get vectors
    eff_size_arr, ϕ1_arr, ϕ2_arr = get_ϕs_eff_rng(eff_type, config)

    # Loop over percent of effect size
    for i = 1:length(ϕ1_arr)

        # Init
        ϕ1 = ϕ1_arr[i]
        ϕ2 = ϕ2_arr[i]
        CpelNano.print_log("Effect percentage $(eff_size_arr[i])")

        # Generate data for multiple reps
        map_out = map(i -> map_run_sim(ϕ1, ϕ2, config), 1:n_sim_reps)

        # Get pvals
        pvals_neg = []
        pvals_pos = []
        for x in map_out
            isnothing(x) && continue
            if x[1]
                push!(pvals_pos, get_pvals(x[2], eff_type))
            else
                push!(pvals_neg, get_pvals(x[2], eff_type))
            end
        end

        # Flatten vectors
        pvals_pos = vcat(pvals_pos...)
        pvals_neg = vcat(pvals_neg...)

        # Remove NaN
        is_nan = isnan.(pvals_pos)
        pvals_pos = pvals_pos[.!is_nan]
        is_nan = isnan.(pvals_neg)
        pvals_neg = pvals_neg[.!is_nan]

        # Store
        for (j, p_thrs) in enumerate(pval_thrs_rng)

            # Count true positives and true negatives
            if length(pvals_pos) > 0
                tp = sum(pvals_pos .<= p_thrs)
                fn = sum(pvals_pos .> p_thrs)
            else
                tp = NaN
                fn = NaN
            end
            if length(pvals_neg) > 0
                tn = sum(pvals_neg .> p_thrs)
                fp = sum(pvals_neg .<= p_thrs)
            else
                tn = NaN
                fp = NaN
            end

            # Store TPR & FPR
            tpr[i,j] = tp / (tp + fn + 1)
            fpr[i,j] = fp / (fp + tn + 1)

        end

    end

    # Return
    return eff_size_arr, tpr, fpr

end
function run_null_sim(n_sim_reps::Int64)

    # CpelNano config
    config = CpelNano.CpelNanoConfig()
    config.max_size_subreg = 350
    config.max_em_iters = 25
    config.max_em_init = 5
    config.verbose = false

    # Init
    ϕ = fill(0.0, 3)

    # Generate data for multiple reps
    map_out = map(i -> map_run_null_sim(ϕ, config), 1:n_sim_reps)

    # Get pvals
    pvals_tmml = []
    pvals_tnme = []
    pvals_tcmd = []
    for x in map_out
        isnothing(x) && continue
        push!(pvals_tmml, get_pvals(x, "mml"))
        push!(pvals_tnme, get_pvals(x, "nme"))
        push!(pvals_tcmd, get_pvals(x, "cmd"))
    end

    # Flatten vectors
    pvals_tmml = vcat(pvals_tmml...)
    pvals_tnme = vcat(pvals_tnme...)
    pvals_tcmd = vcat(pvals_tcmd...)

    # Remove NaN
    is_nan = isnan.(pvals_tmml)
    pvals_tmml = pvals_tmml[.!is_nan]
    is_nan = isnan.(pvals_tnme)
    pvals_tnme = pvals_tnme[.!is_nan]
    is_nan = isnan.(pvals_tcmd)
    pvals_tcmd = pvals_tcmd[.!is_nan]

    # Return pvals
    return pvals_tmml, pvals_tnme, pvals_tcmd

end
function plt_roc(eff_size_arr, tpr, fpr, ttl)

        # Add each effect size
    plt = plot(xlabel="FPR", ylabel="TPR", title=ttl, xlim=(0, 1), ylim=(0, 1), size=(700, 700), legend=:bottomright)
    for (i, eff_size) in enumerate(eff_size_arr)
        plot!(plt, fpr[i,:], tpr[i,:],  seriestype=:line, label="$(round(eff_size;digits=2))", color=blind_friend_col[i])
    end

    # Return plot
    return plt

end
function plt_null_pvals_hist(pvals_tmml, pvals_tnme, pvals_tcmd)

    # Vars
    y_lim = 0.5
    n_bins = 20
    y_label = "percentage (%)"

    # Tmml
    h = fit(Histogram, pvals_tmml, 0:1 / n_bins:1.0)
    midpts = collect(h.edges[1]) .- 2 / n_bins / n_bins
    percs = h.weights ./ sum(h.weights)
    plt_tmml = plot(midpts, percs,  seriestype=:bar, label="", xlabel="",
        ylabel=y_label, title="Tmml", xlim=(0, 1), ylim=y_lim, size=(600, 500))

    # Tnme
    h = fit(Histogram, pvals_tnme, 0:1 / n_bins:1.0)
    midpts = collect(h.edges[1]) .- 2 / n_bins
    percs = h.weights ./ sum(h.weights)
    plt_tnme = plot(midpts, percs,  seriestype=:bar, label="", xlabel="",
        ylabel=y_label, title="Tnme", xlim=(0, 1), ylim=y_lim, size=(600, 500))

    # Tcmd
    h = fit(Histogram, pvals_tcmd, 0:1 / n_bins:1.0)
    midpts = collect(h.edges[1]) .- 2 / n_bins
    percs = h.weights ./ sum(h.weights)
    plt_tcmd = plot(midpts, percs,  seriestype=:bar, label="", xlabel="p-value",
        ylabel=y_label, title="Tcmd", xlim=(0, 1), ylim=y_lim, size=(600, 500))

    # Collage
    plt = plot(plt_tmml, plt_tnme, plt_tcmd, layout=(3, 1), size=(600, 1200))

    # Return plot
    return plt

end
    function plt_null_pvals_ecdf(pvals_tmml, pvals_tnme, pvals_tcmd)

    # Vars
    y_lim = (0, 1.0)
    x_lim = (0, 1.0)
    y_label = "ecdf"
    pval_rng = 0.01:0.01:1.0

    # Tmml
    F = ecdf(pvals_tmml)
    plt_tmml = plot(xlabel=" ", ylabel=y_label, title="Tmml", xlim=x_lim, ylim=y_lim, size=(600, 500))
    plot!(plt_tmml, pval_rng, pval_rng, seriestype=:line, label="Ideal")
    plot!(pval_rng, F.(pval_rng), seriestype=:line, label="Empirical")

    # Tnme
    F = ecdf(pvals_tnme)
    plt_tnme = plot(xlabel=" ", ylabel=y_label, title="Tnme", xlim=x_lim, ylim=y_lim, size=(600, 500))
    plot!(plt_tnme, pval_rng, pval_rng, seriestype=:line, label="Ideal")
    plot!(pval_rng, F.(pval_rng), seriestype=:line, label="Empirical")

    # Tcmd
    F = ecdf(pvals_tcmd)
    plt_tcmd = plot(xlabel="p-value", ylabel=y_label, title="Tcmd", xlim=x_lim, ylim=y_lim, size=(600, 500))
    plot!(plt_tcmd, pval_rng, pval_rng, seriestype=:line, label="Ideal")
    plot!(pval_rng, F.(pval_rng), seriestype=:line, label="Empirical")

    # Collage
    plt = plot(plt_tmml, plt_tnme, plt_tcmd, layout=(3, 1), size=(600, 1200))

    # Return plot
    return plt

end


text:
module x64

using Switch

macro rotl64(var, value)
    :(($var << $value) | ($var >> $(64-value)))
end

macro fmix64(var)
    quote
      local value = $(var)
      value $= (value >> 33)
      value = UInt64(value * UInt64(0xff51afd7ed558ccd))
      value $= (value >> 33)
      value = UInt64(value * UInt64(0xc4ceb9fe1a85ec53))
      value $= (value >> 33)
      value
    end
end

function hash128(data::Array{UInt8}, seed::UInt64)

    c1 = UInt64(0x87c37b91114253d5)
    c2 = UInt64(0x4cf5ad432745937f)

    h1 = seed
    h2 = seed

    len = UInt64(length(data))
    remainder = len & 15
    blocks = div(len, UInt64(16))
    p = convert(Ptr{UInt64}, pointer(data))

    # Body
    iterations = UInt64(blocks*2)
    for next_block::UInt64 = UInt64(1):UInt64(2):iterations
        k1 = UInt64(unsafe_load(p, next_block))
        k2 = UInt64(unsafe_load(p, next_block+1))

        k1 = k1 * c1
        k1 = @rotl64(k1, 31)
        k1 = k1 * c2
        h1 $= k1

        h1 = @rotl64(h1, 27)
        h1 = h1 + h2
        h1 = h1 * UInt64(5) + UInt64(0x52dce729)

        k2 = k2 * c2
        k2 = @rotl64(k2, 33)
        k2 = k2 * c1
        h2 $= k2

        h2 = @rotl64(h2, 31)
        h2 = h2 + h1
        h2 = h2 * UInt64(5) + UInt64(0x38495ab5)
    end

    # Tail
    last = blocks*16
    k1 = UInt64(0)
    k2 = UInt64(0)

    @switch remainder begin
        @case 15
            k2 $= UInt64(data[last+15]) << 48
        @case 14
            k2 $= UInt64(data[last+14]) << 40
        @case 13
            k2 $= UInt64(data[last+13]) << 32
        @case 12
            k2 $= UInt64(data[last+12]) << 24
        @case 11
            k2 $= UInt64(data[last+11]) << 16
        @case 10
            k2 $= UInt64(data[last+10]) << 8
        @case 9
            k2 $= UInt64(data[last+9]) << 0
            k2 = k2 * c2
            k2 = @rotl64(k2, 33)
            k2 = k2 * c1
            h2 $= k2
        @case 8
            k1 $= UInt64(data[last+8]) << 56
        @case 7
            k1 $= UInt64(data[last+7]) << 48
        @case 6
            k1 $= UInt64(data[last+6]) << 40
        @case 5
            k1 $= UInt64(data[last+5]) << 32
        @case 4
            k1 $= UInt64(data[last+4]) << 24
        @case 3
            k1 $= UInt64(data[last+3]) << 16
        @case 2
            k1 $= UInt64(data[last+2]) << 8
        @case 1
            k1 $= UInt64(data[last+1])
            k1 = k1 * c1
            k1 = @rotl64(k1, 31)
            k1 = k1 * c2
            h1 $= k1
    end

    # Finalization
    h1 $= len
    h2 $= len
    h1 = h1 + h2
    h2 = h2 + h1
    h1 = @fmix64(h1)
    h2 = @fmix64(h2)
    h1 = h1 + h2
    h2 = h2 + h1

    h1 = @rotl64(h1, 32)
    h2 = @rotl64(h2, 32)

    return (UInt128(h1) << 64) | UInt128(h2)

end

hash128(data::Array{UInt8}, seed::Int) = hash128(data, UInt64(seed))
hash128(data::Array{UInt8}) = hash128(data, UInt64(0))
hash128(data::AbstractString, seed::UInt64) = hash128(convert(Array{UInt8}, data), seed)
hash128(data::AbstractString, seed::Int) = hash128(convert(Array{UInt8}, data), UInt64(seed))
hash128(data::AbstractString) = hash128(convert(Array{UInt8}, data), UInt64(0))

end


text:
using Revise
using BayesMultiFit
BAMF=BayesMultiFit
using ReversibleJumpMCMC
const RJMCMC = ReversibleJumpMCMC
using Plots
using Distributions
using MicroscopePSFs
const PSF=MicroscopePSFs
using ProfileView
using BenchmarkTools

## simulation config
n=6             # number of emitters
μ=1000          # mean photons per emitter
sz=16           # ROI size in pixels
iterations=1000 # RJMCMC iterations
burnin=1000     # RJMCMC iterations for burn-in

## PSF config
σ=1.3f0 # Gaussian PSF Sigma in Pixels
pixelsize=1.0
psf=PSF.Gauss2D(σ,pixelsize)
 psf=PSF.InterpolatedPSF(psf,(sz*2,sz*2))

## setup prior distribution on intensity
α=Float32(4)
θ=Float32(μ/α)
g=Gamma(α,θ)
len=Int32(1024)
θ_start=Float32(1)
θ_step=Float32(5.0)
pdf_x=range(θ_start,step=θ_step,length=len)
mypdf=pdf(g,pdf_x)
plt=plot(pdf_x,mypdf)
display(plt)
prior_photons=BAMF.RJPrior(len,θ_start,θ_step,mypdf)

# set emitter positions and intensity
x=sz*rand(Float32,n)
y=sz*rand(Float32,n)

photons=Float32.(rand(g,n))
bg=1f-6
datastate=BAMF.StateFlatBg(n,x,y,photons,bg)

## Create synthetic data
data=BAMF.ArrayDD(sz)

@code_warntype BAMF.genmodel!(datastate,psf,data)
@btime BAMF.genmodel!(datastate,psf,data)
@benchmark BAMF.genmodel!(datastate,psf,data)

roi=[(i,j) for i=1:sz,j=1:sz]
@benchmark MicroscopePSFs.pdf(psf,(0,0),(0,0))
@benchmark MicroscopePSFs.pdf(psf,roi,(0,0))

## Create synthetic data
data=BAMF.ArrayDD(sz)
BAMF.genmodel!(datastate,psf,data)
BAMF.poissrnd!(data.data)
# imshow(data.data)  #look at data

## create a BAMF-type RJMCMC structure
xystd=σ/10
istd=10f0
split_std=σ/2
bndpixels=0f0
myRJ=BAMF.RJStruct(sz,psf,xystd,istd,split_std,data,bndpixels,prior_photons,BAMF.ArrayDD(sz),BAMF.ArrayDD(sz))

## setup the RJMCMC.jl model
# Jumptypes are: move, bg, add, remove, split, merge
njumptypes=6
jumpprobability=[1,0,.1,.1,.1,.1] # Model with no bg
jumpprobability=jumpprobability/sum(jumpprobability)

# create an RJMCMC structure with all model info
acceptfuns=[BAMF.accept_move,BAMF.accept_bg,BAMF.accept_add,BAMF.accept_remove,BAMF.accept_split,BAMF.accept_merge] #array of functions
propfuns=[BAMF.propose_move,BAMF.propose_bg,BAMF.propose_add,BAMF.propose_remove,BAMF.propose_split,BAMF.propose_merge] #array of functions
myRJMCMC=RJMCMC.RJMCMCStruct(burnin,iterations,njumptypes,jumpprobability,propfuns,acceptfuns)
#create an intial state
state1=BAMF.calcintialstate(myRJ)

## run chain. This is the call to the main algorithm
@code_warntype RJMCMC.buildchain(myRJMCMC,myRJ,state1)
ProfileView.@profview RJMCMC.buildchain(myRJMCMC,myRJ,state1)
@benchmark mychain=RJMCMC.buildchain(myRJMCMC,myRJ,state1)


text:
SPL(p::Missing) = NaN
SPL(p::Number) = p > 0.0 ? 10*log10(p/4e-10) : NaN

"""
    octavebandlimits(fc,n)

Compute 1/n octave band limits given center frequencies fc.
"""
function octavebandlimits(fc::Vector{T},n) where T <: Real
    fl = fu = similar(fc)
    C = 10^(3. /10.)
    fl = fc*C^(-1/(2*n))
    fu = fc*C^(1/(2*n))
    return fl, fu
end

function octavebandlimits(fc::T,n) where T <: Real
    C = 10^(3. /10.)
    fl = fc*C^(-1/(2*n))
    fu = fc*C^(1/(2*n))
    return fl, fu
end

selectfreqs(f,flim) = (f.>=flim[1]) .& (f.<=flim[2])

function octavebands_nomial(n)
    if n == 1
        return [16, 31.5, 63.0, 125.0, 250.0, 500.0, 1000.0, 2000.0, 4000.0, 8000.0, 16000.0, 32000.0]
    elseif n == 3
        return [12.5, 16, 20, 25.0, 31.5, 40.0, 50.0, 63.0, 80.0, 100.0, 125.0, 160.0, 200.0,
            250.0, 315.0, 400.0, 500.0, 630.0, 800.0, 1000.0,1250.0, 1600.0, 2000.0, 2500.0, 3150.0, 4000.0,
            5000.0, 6300.0, 8000.0, 10000.0, 12500.0, 16000.0, 20000.0, 25000.0, 31500.0, 40000.0]
    else
        error("Value of `n` not supported. Try `n = 1` or `n = 3`")
    end
end

"""
    octavebands(n,flim=(16.,40000.),nomial::Bool=false)

Compute 1/n octave band center frequencies.
"""
function octavebands(n,flim=(16.,40000.);nomial::Bool=false)
    if nomial && (n==1 || n==3)
        fc = octavebands_nomial(n)
        f_inds = selectfreqs(fc,flim)
        return fc[f_inds]
    else
        fc3 = octavebands_nomial(3)
        smin,imin = findmin(abs.(fc3.-flim[1]))
        smax,imax = findmin(abs.(fc3.-flim[2]))
        (smin != 0) ? @warn("low frequency limit is not supported and has been replaced by $(fc3[imin])") : nothing
        (smax != 0) ? @warn("high frequency limit is not supported and has been replaced by $(fc3[imax])") : nothing
        flim = (fc3[imin],fc3[imax])
        G = 10^(3. /10.)
        fc = exp.(log(flim[1]):(log(G)/n):log(flim[2]))
        return fc
    end
end

"""
    coherence_weights(CSM)

Calculate microphone weights according to mean mututal-coherence.

Amaral et al. “Slat Noise from an MD30P30N Airfoil at Extreme Angles of Attack,” AIAA Journal, vol. 56, no. 3, pp. 964–978, Mar. 2018.
"""
function coherence_weights(CSM)
    M,M,Nf = size(CSM)
    weight = zeros(M,Nf)
    for ω = 1:Nf
        for m = 1:M
            for n = 1:M
                if m != n
                    weight[m,ω] += abs(CSM[m,n,ω])^2 / abs(CSM[m,m,ω]*CSM[n,n,ω])
                end
            end
            weight[m,ω] /= M-1
        end
    end
    return weight
end

"""
    enbw(fs,win)

Equivalent Noise Band Width. Use to convert power spectrum to power spectral density and reverse. PS = PSD*ENBW
"""
enbw(fs,win) = fs*sum(abs2,win)/sum(win)^2


text:
export MPS

struct MPS{T, Tensor <: AbstractArray{T, 3}, N} <: AbstractBlock
    tensors::NTuple{N, Tensor}
end

# check if bonds match
# tensors should be an iterable
function check_bonds(tensors)
    for (prev, curr) in zip(tensors[1:end-1], tensors[2:end])
        size(prev, 3) == size(curr, 1) || return false
    end

    size(tensors[1], 1) == size(tensors[end], 3) || return false
    return true
end

function MPS(tensors::AbstractArray{T, 3}...) where T
    # TODO: custom exception
    check_bonds(tensors) || throw(ErrorException("MPS bonds mismatch"))
    MPS{T, eltype(tensors), length(tensors)}(tensors)
end

function MPS(::Type{T}, bonds::NTuple{M, I}, physicals::NTuple{N, I}) where {T, M, N, I <: Integer}
    @assert M == N + 1 "Physical legs and number of tensors mismatch"

    tensors = ntuple(i->rand(T, bonds[i], physicals[i], bonds[i+1]), N)
    MPS{T, Array{T, 3}, N}(tensors)
end

MPS(bonds::NTuple{M, I}, physicals::NTuple{N, I}) where {M, I, N} =
    MPS(Float64, bonds, physicals)
MPS(::Type{T}, bonds::NTuple{N, I}) where {T, N, I <: Integer} = MPS(T, bonds, ntuple(x->2, N-1))
MPS(::Type{T}, bonds::Integer...) where T = MPS(T, bonds)
MPS(bonds::Integer...) = MPS(Float64, bonds)

# MPS should be iteratable

import Base: length, eltype, start, next, done, getindex

length(mps::MPS{T, Tensor, N}) where {T, Tensor, N} = N
eltype(mps::MPS{T, Tensor, N}) where {T, Tensor, N} = Tensor
start(mps::MPS) = start(mps.tensors)
next(mps::MPS, state) = next(mps.tensors, state)
done(mps::MPS, state) = next(mps.tensor, state)
getindex(mps::MPS, index::Integer) = getindex(mps.tensor, index)

export virtual_bonds, physical_bonds, virtual_bond, physical_bond

virtual_bonds(mps::MPS) = push!([size(each, 1) for each in mps.tensors], size(mps.tensors[end], 3))

function virtual_bond(mps::MPS{T, TT, N}, index::Integer) where {T, TT, N}
    index == N+1 && return size(mps.tensors[end], 3)
    return size(mps.tensors[index], 1)
end

physical_bonds(mps::MPS) = [size(each, 2) for each in mps.tensors]
physical_bond(mps::MPS, index::Integer) = size(mps.tensors[index], 2)

import Base: show

"""
MPS{Float64, Array{Float64, 3}, 10}:
2 -[2]- 10 -[2]- 10 -[2]- 2
"""
function show_mps(io::IO, mps::MPS{T, TT, N}) where {T, TT, N}
    mps_str = ""
    physicals = physical_bonds(mps)
    virtuals = virtual_bonds(mps)
    for i = 1:N
        mps_str *= string(virtuals[i]) * " -[$(physicals[i])]- "
    end
    mps_str *= " $(virtuals[end])"
    print(io, mps_str)
end

@static if VERSION < v"0.7-"

function show(io::IO, ::MIME"text/plain", mps::MPS{T, TT, N}) where {T, TT, N}
    print(io, summary(mps))
    println(io, ":")
    show_mps(io, mps)
end

else # for v0.7+

function show(io::IO, ::MIME"text/plain", mps::MPS{T, TT, N}) where {T, TT, N}
    summary(io, mps)
    println(io, ":")
    show_mps(io, mps)
end

end

import Compat: LinearAlgebra

# TODO: cache the contraction
function forward(mps::MPS, input::AbstractVector)
    itr = zip(mps.tensors[1:end-1], mps.tensors[2:end],
                input[1:end-1], input[2:end])
    for (prev, curr, prev_i, curr_i) in itr
        LinearAlgebra.BLAS.gemm('N', 'N', prev[:, prev_i, :], curr[:, curr_i, :])
    end
    return itr
end


text:
##################################
# Crude, i.e., KBS1 implementation
##################################

function _iscritical(u::AbstractWord, v::AbstractWord, rewriting)
    a = rewrite_from_left(u, rewriting)
    b = rewrite_from_left(v, rewriting)
    return a ≠ b, (a,b)
end

"""
    deriverule!(rws::RewritingSystem, u::Word, v::Word[, o::Ordering=ordering(rws)])
Given a critical pair `(u, v)` with respect to `rws` adds a rule to `rws`
(if necessary) that solves the pair, i.e. makes `rws` locally confluent with
respect to `(u,v)`. See [Sims, p. 69].
"""
function deriverule!(rws::RewritingSystem{W}, u::AbstractWord, v::AbstractWord,
    o::Ordering = ordering(rws)) where W

    critical, (a, b) = _iscritical(u, v, rws)
    if critical
        simplifyrule!(a, b, o)
        push!(rws, Rule{W}(a, b, o))
    end
    return rws
end

##########################
# Naive KBS implementation
##########################

# As of now: default implementation

function _iscritical(u::AbstractWord, v::AbstractWord, rewriting, work::kbWork)
    a = rewrite_from_left!(work.lhsPair, u, rewriting)
    b = rewrite_from_left!(work.rhsPair, v, rewriting)
    return a ≠ b, (a, b)
end

"""
    deriverule!(rs::RewritingSystem, stack, work::kbWork
        [, o::Ordering=ordering(rs), deleteinactive::Bool = false])
Empty `stack` of (potentially) critical pairs by deriving and adding new rules
to `rs` resolving the pairs, i.e. maintains local confluence of `rs`.

This function may deactivate rules in `rs` if they are deemed redundant (e.g.
follow from the added new rules). See [Sims, p. 76].
"""
function deriverule!(rws::RewritingSystem{W}, stack, work::kbWork,
    o::Ordering = ordering(rws)) where W
    while !isempty(stack)
        u, v = pop!(stack)
        critical, (a, b) = _iscritical(u, v, rws, work)
        if critical
            simplifyrule!(a, b, o)
            new_rule = Rule{W}(a, b, o)
            deactivate_rules!(rws, stack, work, new_rule)
            push!(rws, new_rule)
        end
    end
end

function deactivate_rules!(rws::RewritingSystem, stack, work::kbWork, new_rule::Rule)
    for rule in rules(rws)
        rule == new_rule && continue
        (lhs, rhs) = rule
        if occursin(new_rule.lhs, lhs)
            deactivate!(rule)
            push!(stack, (first(rule), last(rule)))
        elseif occursin(new_rule.lhs, rhs)
            new_rhs = rewrite_from_left!(work.rhsPair, rhs, rws)
            update_rhs!(rule, new_rhs)
        end
    end
end

########################################
# KBS using index automata for rewriting
########################################

"""
    deriverule!(rs::RewritingSystem, at::Automaton, stack
    [,work = nothing, o::Ordering=ordering(rs))])
Adds a rule to a rewriting system and deactivates others (if necessary) that
insures that the set of rules is reduced while maintaining local confluence.
See [Sims, p. 76].
"""
function deriverule!(rs::RewritingSystem{W}, stack,
        work::kbWork, at::Automaton, o::Ordering = ordering(rs)) where {W<:AbstractWord}
    if length(stack) >= 2
        @debug "Deriving rules with stack of length=$(length(stack))"
    end
    while !isempty(stack)
        u, v = pop!(stack)
        critical, (a, b) = _iscritical(u, v, at, work)

        if critical
            simplifyrule!(a, b, o)
            new_rule = Rule{W}(a, b, o)
            push!(rs, new_rule)
            # TODO: push!(at, new_rule)
            updateautomaton!(at, rs)

            for rule in rules(rs)
                rule == new_rule && break
                (lhs, rhs) = rule
                if occursin(new_rule.lhs, lhs)
                    deactivate!(rule)
                    push!(stack, (first(rule), last(rule)))
                    # TODO: delete!(at, rule)
                    updateautomaton!(at, rs)
                elseif occursin(new_rule.lhs, rhs)
                    new_rhs = rewrite_from_left!(work.rhsPair, rhs, at)
                    update_rhs!(rule, new_rhs)
                    # TODO: push!(at, rule)
                    updateautomaton!(at, rs)
                end
            end
        end
    end
end


text:
module TestWiringDiagramAlgorithms

using Test
using Catlab.Theories, Catlab.WiringDiagrams

A, B, C, D = [ Ports([sym]) for sym in [:A, :B, :C, :D] ]
I = munit(Ports)
f = singleton_diagram(Box(:f, [:A], [:B]))
g = singleton_diagram(Box(:g, [:B], [:C]))
h = singleton_diagram(Box(:h, [:C], [:D]))

# Normalization
###############

# Normalize copies.
d = compose(mcopy(A), otimes(f,f))
normalized = compose(f, mcopy(B))
@test normalize_copy!(d) == normalized

d = compose(f, mcopy(B), otimes(g,g))
normalize_copy!(d)
normalized = compose(f, g, mcopy(C))
perm = sortperm(boxes(d); by=box->box.value)
@test is_permuted_equal(d, normalized, perm)

d = compose(mcopy(A), otimes(f,f), otimes(g,g))
normalize_copy!(d)
perm = sortperm(boxes(d); by=box->box.value)
@test is_permuted_equal(d, normalized, perm)

# Normalize deletions.
@test normalize_delete!(copy(f)) == f

d = WiringDiagram(I,I)
add_box!(d, first(boxes(f)))
@test normalize_delete!(d) == WiringDiagram(I,I)

d = WiringDiagram(A, B)
fv = add_box!(d, first(boxes(f)))
gv = add_box!(d, first(boxes(g)))
hv = add_box!(d, first(boxes(h)))
add_wires!(d, [
  (input_id(d),1) => (fv,1),
  (fv,1) => (gv,1),
  (gv,1) => (hv,1),
  (fv,1) => (output_id(d),1),
])
@test normalize_delete!(d) == f

# Normalize wiring diagrams representing morphisms in a cartesian category.
d = compose(
  mcopy(A),
  otimes(id(A),mcopy(A)),
  otimes(f,f,f),
  otimes(id(B), id(B), compose(g, delete(C)))
)
normalized = compose(f, mcopy(B))
@test normalize_cartesian!(d) == normalized

# Layout
########

@test crossing_minimization_by_sort(f, box_ids(f), sources=[input_id(d)]) == [1]

d = otimes(f,g)
@test crossing_minimization_by_sort(d, box_ids(d), sources=[input_id(d)]) == [1,2]

d = WiringDiagram(I,I)
fv1, fv2 = add_box!(d, first(boxes(f))), add_box!(d, first(boxes(f)))
gv1, gv2 = add_box!(d, first(boxes(g))), add_box!(d, first(boxes(g)))
add_wires!(d, [ Wire((fv1,1) => (gv1,1)), Wire((fv2,1),(gv2,1)) ])
@test crossing_minimization_by_sort(d, [gv1,gv2], sources=[fv1,fv2]) == [1,2]
@test crossing_minimization_by_sort(d, [fv1,fv2], targets=[gv1,gv2]) == [1,2]

end


text:
using CellularInfo, JSON

graphs = open("graphs.json", "r") do handle
    graphs = Dict{Int, Vector{Tuple{Int,Int}}}()
    for (k, v) in JSON.parse(handle)
        graphs[parse(Int, k)] = map(x -> tuple(x...), v)
    end
    graphs
end

for (rule, class) in EQUIVALENT_RULES
    for partner in class
        if graphs[rule] != graphs[partner]
            println("$rule != $(partner)")
            println(graphs[rule])
            println(graphs[partner])
        end
    end
end

max_stars_repo_name: ELIFE-ASU/dglmoore

max_issues_repo_name: ELIFE-ASU/dglmoore

max_forks_repo_name: ELIFE-ASU/dglmoore


text:
export singleton, set, card, fullset

const EntropyIndex = UInt

# Set Manipulation
# function subsets(S)
#   if S == 0
#     return []
#   else
#     cum = [S]
#     for i = 0:(n-1)
#       if (S & (1 << i)) != 0
#         append!(cum, subsets(S $ (1 << i)))
#       end
#     end
#     return cum
#   end
# end
#
# function subsetsones(S)
#   ret = zeroentropy(N)
#   ret[subsets(S)] = 1
#   return ret
# end

function Base.setdiff(S::EntropyIndex, I::EntropyIndex)
    S & (~I)
end

# S ∪ T
function Base.union(S::EntropyIndex, T::EntropyIndex)
    S | T
end

# S ∩ T
function Base.intersect(S::EntropyIndex, T::EntropyIndex)
    S & T
end

# S ⊆ T
function issubset(S::EntropyIndex, T::EntropyIndex)
    Base.setdiff(S, T) == zero(S)
end

function singleton(i::Integer)
    EntropyIndex(1) << (i-1)
end

function fullset(n::Integer)
    (EntropyIndex(1) << n) - EntropyIndex(1)
end

"""
    set(i::Integer)

Return the set of digits of the integer `i`. Note that this cannot create a set
with elements other than 1, 2, 3, 4, 5, 6, 7, 8 and 9. For instance, the set
``{2, 4, 10}`` cannot be constructed with this method, use `set([2, 4, 10])`
instead.

## Examples

To create the set ``\\{2, 4\\}``, use `set(24)` or `set(42)`.
"""
function set(i::Integer)
    ret = emptyset()
    while i > 0
        ret = union(ret, singleton(i % 10))
        i = div(i, 10)
    end
    ret
end

"""
    set(I::AbstractArray{<:Integer})

Return the set of elements of `I`.

## Examples

To create the set ``\\{2, 4\\}``, use `set([2, 4])` or `set([4, 2])`.
"""
function set(I::AbstractArray{<:Integer})
    ret = emptyset()
    for i in I
        ret = union(ret, singleton(i))
    end
    ret
end

function myin(i::Signed, I::EntropyIndex)
    (singleton(i) & I) != 0
end

function card(S::EntropyIndex)
    sum = 0
    while S > 0
        sum += S & 1
        S >>= 1
    end
    Signed(sum)
end

function emptyset()
    EntropyIndex(0)
end

function setsto(J::EntropyIndex)
    EntropyIndex(1):J
end

function mymap(map::Vector, S::EntropyIndex, n)
    T = emptyset()
    for i in 1:n
        if myin(i, S)
            T = T ∪ singleton(map[i])
        end
    end
    T
end


text:
export subd_chart, chart, getelementVertices
mutable struct subd_chart
    index::Int64
    N::Int64
    RingNodes::Vector{Int64}
    vertices::Vector{SVector{3,Float64}}
end

domain(ch::subd_chart) = ReferenceSimplex{2,Float64,3}()

function chart(Smesh::subdMesh,E)
    element = Smesh.elements[E]
    Svertices = Smesh.vertices
    nodes = element.RingNodes
    N = length(nodes)
    verticecoords = Vector{SVector{3,Float64}}(N)

    for i = 1:N
        verticecoords[i] = Svertices[nodes[i]].Coords
    end
    return chart = subd_chart(E,N,nodes,verticecoords)
end

function getelementVertices(chart::subd_chart)
    verticescoords = chart.vertices
    verticecoords3 = Vector{SVector{3,Float64}}(3)
    verticecoords3[1] = verticescoords[1]
    verticecoords3[2] = verticescoords[2]
    verticecoords3[3] = verticescoords[chart.N-6]
    return verticecoords3
end


text:
# SQRT FUNCTION
# (x+1)^0.5
function sum_sqrt(x::Array{Float64,1})
    return sum(sqrt(x+1))
end

function d_sum_sqrt(x::Array{Float64,1})
    return 0.5 * (x+1).^(-0.5)
end

function dd_sum_sqrt(x::Array{Float64,1})
    return spdiagm(-0.25 * (x+1).^(-1.5))
end

function nl_sqrt(num_var::Int64)
    return class_nl_function(num_var,sum_sqrt,d_sum_sqrt,dd_sum_sqrt);
end

# LOG function
function sum_log(x::Array{Float64,1})
    return sum(log(x+1))
end

function d_sum_log(x::Array{Float64,1})
    return (x+1).^(-1.0)
end

function dd_sum_log(x::Array{Float64,1})
    return spdiagm(-(x+1).^(-2.0))
end

function nl_log(num_var::Int64)
    return class_nl_function(num_var::Int64,sum_log,d_sum_log,dd_sum_log);
end

# Quadratic
# c'*x + 0.5 * x*Q*x + b
function nl_quad(c::AbstractArray, Q::SparseMatrixCSC{Float64,Int64}, b::Float64)
    function value(x::Array{Float64,1})
        return (c' * x)[1] + 0.5 * (x' * Q * x)[1] + b
    end

    function gradient(x::Array{Float64,1})
        return (c + Q * x)[:]
    end

    function hessian(x::Array{Float64,1})
        return Q
    end

    return class_nl_function(length(c), value, gradient, hessian)
end

function sphere_function(center_point::Float64,radius::Float64)
    b = dot(center_point,center_point) - radius
    c = -2 * center_point

    return nl_quad(c, speye(length(c)), b)
end

max_forks_repo_name: ohinder/Everest


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

mutable struct DeploymentWhatIfProperties <: SwaggerModel
    template::Any # spec type: Union{ Nothing, Any } # spec name: template
    templateLink::Any # spec type: Union{ Nothing, TemplateLink } # spec name: templateLink
    parameters::Any # spec type: Union{ Nothing, Any } # spec name: parameters
    parametersLink::Any # spec type: Union{ Nothing, ParametersLink } # spec name: parametersLink
    mode::Any # spec type: Union{ Nothing, String } # spec name: mode
    debugSetting::Any # spec type: Union{ Nothing, DebugSetting } # spec name: debugSetting
    onErrorDeployment::Any # spec type: Union{ Nothing, OnErrorDeployment } # spec name: onErrorDeployment

    function DeploymentWhatIfProperties(;template=nothing, templateLink=nothing, parameters=nothing, parametersLink=nothing, mode=nothing, debugSetting=nothing, onErrorDeployment=nothing)
        o = new()
        validate_property(DeploymentWhatIfProperties, Symbol("template"), template)
        setfield!(o, Symbol("template"), template)
        validate_property(DeploymentWhatIfProperties, Symbol("templateLink"), templateLink)
        setfield!(o, Symbol("templateLink"), templateLink)
        validate_property(DeploymentWhatIfProperties, Symbol("parameters"), parameters)
        setfield!(o, Symbol("parameters"), parameters)
        validate_property(DeploymentWhatIfProperties, Symbol("parametersLink"), parametersLink)
        setfield!(o, Symbol("parametersLink"), parametersLink)
        validate_property(DeploymentWhatIfProperties, Symbol("mode"), mode)
        setfield!(o, Symbol("mode"), mode)
        validate_property(DeploymentWhatIfProperties, Symbol("debugSetting"), debugSetting)
        setfield!(o, Symbol("debugSetting"), debugSetting)
        validate_property(DeploymentWhatIfProperties, Symbol("onErrorDeployment"), onErrorDeployment)
        setfield!(o, Symbol("onErrorDeployment"), onErrorDeployment)
        o
    end
end # type DeploymentWhatIfProperties

const _property_map_DeploymentWhatIfProperties = Dict{Symbol,Symbol}(Symbol("template")=>Symbol("template"), Symbol("templateLink")=>Symbol("templateLink"), Symbol("parameters")=>Symbol("parameters"), Symbol("parametersLink")=>Symbol("parametersLink"), Symbol("mode")=>Symbol("mode"), Symbol("debugSetting")=>Symbol("debugSetting"), Symbol("onErrorDeployment")=>Symbol("onErrorDeployment"))
const _property_types_DeploymentWhatIfProperties = Dict{Symbol,String}(Symbol("template")=>"Any", Symbol("templateLink")=>"TemplateLink", Symbol("parameters")=>"Any", Symbol("parametersLink")=>"ParametersLink", Symbol("mode")=>"String", Symbol("debugSetting")=>"DebugSetting", Symbol("onErrorDeployment")=>"OnErrorDeployment")
Base.propertynames(::Type{ DeploymentWhatIfProperties }) = collect(keys(_property_map_DeploymentWhatIfProperties))
Swagger.property_type(::Type{ DeploymentWhatIfProperties }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_DeploymentWhatIfProperties[name]))}
Swagger.field_name(::Type{ DeploymentWhatIfProperties }, property_name::Symbol) =  _property_map_DeploymentWhatIfProperties[property_name]

const _allowed_DeploymentWhatIfProperties_mode = ["Incremental", "Complete"]

function check_required(o::DeploymentWhatIfProperties)
    true
end

function validate_property(::Type{ DeploymentWhatIfProperties }, name::Symbol, val)
    if name === Symbol("mode")
        Swagger.validate_param(name, "DeploymentWhatIfProperties", :enum, val, _allowed_DeploymentWhatIfProperties_mode)
    end
end


text:
# An implementation of the Beaufort cipher in Julia

function beaufort()
    #Get plaintext and key
    println("Please enter the text: ")
    text = readline()
    println("Please enter the keyword: ")
    keyword = readline()
    text = uppercase(text)
    keyword = uppercase(keyword)
    key = ""

    # The keyword must be repeated until it matches the length of the text
    for i = firstindex(text):lastindex(text)
        #pos = ((c - 'A') % length(keyword)) + 1
        pos = i
        while pos > length(keyword)
            pos -= length(keyword)
        end
        key = string(key, keyword[pos])
    end

    # Now encrypt. Works the same as Vigenere, except subtracts rather than adding distance and the key and plaintext have opposite roles
    index = 1
    for char in key
        distance = text[index] - 'A'
        newchar = char - distance
        if(newchar > 'Z')
            newchar -= 26
        end
        if(newchar < 'A')
            newchar += 26
        end
        index += 1
        print(newchar)
    end
end
beaufort()


text:
##
const ZEROS_DOC = """
    @zeros(args...)

!!! note "Advanced"
        @zeros(numbertype, args...)

Call `zeros(numbertype, args...)`, where `numbertype` is the datatype selected with [`@init_parallel_kernel`](@ref) and the function `zeros` is chosen to be compatible with the package for parallelization selected with [`@init_parallel_kernel`](@ref) (zeros for Threads and CUDA.zeros for CUDA).

!!! note "Advanced"
    The `numbertype` can be explicitly passed as argument in order to be used instead of the default `numbertype` chosen with [`@init_parallel_kernel`](@ref). If no default `numbertype` was chosen [`@init_parallel_kernel`](@ref), then the argument `numbertype` is mandatory. This needs to be used with care to ensure that no datatype conversions occur in performance critical computations.
"""
@doc ZEROS_DOC
macro zeros(args...)
    check_initialized();
    esc(_zeros(args...));
end

##
const ONES_DOC = """
    @ones(args...)

!!! note "Advanced"
        @ones(numbertype, args...)

Call `ones(numbertype, args...)`, where `numbertype` is the datatype selected with [`@init_parallel_kernel`](@ref) and the function `ones` is chosen to be compatible with the package for parallelization selected with [`@init_parallel_kernel`](@ref) (ones for Threads CUDA.ones for CUDA).

!!! note "Advanced"
    The `numbertype` can be explicitly passed as argument in order to be used instead of the default `numbertype` chosen with [`@init_parallel_kernel`](@ref). If no default `numbertype` was chosen [`@init_parallel_kernel`](@ref), then the argument `numbertype` is mandatory. This needs to be used with care to ensure that no datatype conversions occur in performance critical computations.
"""
@doc ONES_DOC
macro ones(args...)
    check_initialized();
    esc(_ones(args...));
end

##
const RAND_DOC = """
    @rand(args...)

!!! note "Advanced"
        @rand(numbertype, args...)

Call `rand(numbertype, args...)`, where `numbertype` is the datatype selected with [`@init_parallel_kernel`](@ref) and the function `rand` is chosen/implemented to be compatible with the package for parallelization selected with [`@init_parallel_kernel`](@ref).

!!! note "Advanced"
    The `numbertype` can be explicitly passed as argument in order to be used instead of the default `numbertype` chosen with [`@init_parallel_kernel`](@ref). If no default `numbertype` was chosen [`@init_parallel_kernel`](@ref), then the argument `numbertype` is mandatory. This needs to be used with care to ensure that no datatype conversions occur in performance critical computations.
"""
@doc RAND_DOC
macro rand(args...)
    check_initialized();
    esc(_rand(args...));
end

## MACROS FORCING PACKAGE, IGNORING INITIALIZATION

macro zeros_cuda(args...)    check_initialized(); esc(_zeros(args...; package=PKG_CUDA)); end
macro ones_cuda(args...)     check_initialized(); esc(_ones(args...; package=PKG_CUDA)); end
macro rand_cuda(args...)     check_initialized(); esc(_rand(args...; package=PKG_CUDA)); end
macro zeros_threads(args...) check_initialized(); esc(_zeros(args...; package=PKG_THREADS)); end
macro ones_threads(args...)  check_initialized(); esc(_ones(args...; package=PKG_THREADS)); end
macro rand_threads(args...)  check_initialized(); esc(_rand(args...; package=PKG_THREADS)); end

## ALLOCATOR FUNCTIONS

function _zeros(args...; package::Symbol=get_package())
    numbertype = get_numbertype()
    if     (package == PKG_CUDA)    return :(ParallelStencil.ParallelKernel.zeros_gpu($numbertype, $(args...)))
    elseif (package == PKG_THREADS) return :(ParallelStencil.ParallelKernel.zeros_cpu($numbertype, $(args...)))
    else                            @KeywordArgumentError("$ERRMSG_UNSUPPORTED_PACKAGE (obtained: $package).")
    end
end

function _ones(args...; package::Symbol=get_package())
    numbertype = get_numbertype()
    if     (package == PKG_CUDA)    return :(ParallelStencil.ParallelKernel.ones_gpu($numbertype, $(args...)))
    elseif (package == PKG_THREADS) return :(ParallelStencil.ParallelKernel.ones_cpu($numbertype, $(args...)))
    else                            @KeywordArgumentError("$ERRMSG_UNSUPPORTED_PACKAGE (obtained: $package).")
    end
end

function _rand(args...; package::Symbol=get_package())
    numbertype = get_numbertype()
    if     (package == PKG_CUDA)    return :(ParallelStencil.ParallelKernel.rand_gpu($numbertype, $(args...)))
    elseif (package == PKG_THREADS) return :(ParallelStencil.ParallelKernel.rand_cpu($numbertype, $(args...)))
    else                            @KeywordArgumentError("$ERRMSG_UNSUPPORTED_PACKAGE (obtained: $package).")
    end
end

## RUNTIME ALLOCATOR FUNCTIONS

zeros_cpu(T0::DataType, T::DataType, args...) = Base.zeros(T, args...)               # If the user has passed a numbertype as first argument to @zeros, then the default numbertype added automatically by the @zeros macro (T0) is not used.
 ones_cpu(T0::DataType, T::DataType, args...) = Base.ones(T, args...)                # ...
 rand_cpu(T0::DataType, T::DataType, args...) = Base.rand(T, args...)                # ...

zeros_gpu(T0::DataType, T::DataType, args...) = CUDA.zeros(T, args...)               # ...
 ones_gpu(T0::DataType, T::DataType, args...) = CUDA.ones(T, args...)                # ...
 rand_gpu(T0::DataType, T::DataType, args...) = CUDA.CuArray(Base.rand(T, args...))  # ...

zeros_cpu(T0::DataType, args...) = (check_default_numbertype(T0); Base.zeros(T0, args...))
 ones_cpu(T0::DataType, args...) = (check_default_numbertype(T0); Base.ones(T0, args...))
 rand_cpu(T0::DataType, args...) = (check_default_numbertype(T0); Base.rand(T0, args...))

zeros_gpu(T0::DataType, args...) = (check_default_numbertype(T0); CUDA.zeros(T0, args...))
 ones_gpu(T0::DataType, args...) = (check_default_numbertype(T0); CUDA.ones(T0, args...))
 rand_gpu(T0::DataType, args...) = (check_default_numbertype(T0); CUDA.CuArray(Base.rand(T0, args...)))

function check_default_numbertype(T0::DataType)
    if (T0 == NUMBERTYPE_NONE) @ArgumentError("the numbertype argument is mandatory in @zeros, @ones and @rand when no default is set.") end
end


text:
#= iterators from https://lostella.github.io/blog/2018/07/25/iterative-methods-done-right =#
#Halting
struct HaltingIterable{I, F}
    iter::I
    fun::F
end

function iterate(iter::HaltingIterable)
    next = iterate(iter.iter)
    return dispatch(iter, next)
end

function iterate(iter::HaltingIterable, (instruction, state))
    if instruction == :halt return nothing end
    next = iterate(iter.iter, state)
    return dispatch(iter, next)
end

function dispatch(iter::HaltingIterable, next)
    next == nothing && return nothing
    return next[1], (iter.fun(next[1]) ? :halt : :continue, next[2])
end

halt(iter::I, fun::F) where {I, F} = HaltingIterable{I, F}(iter, fun)

#Tee
struct TeeIterable{I, F}
    iter::I
    fun::F
end

function iterate(iter::TeeIterable, args...)
    next = iterate(iter.iter, args...)
    if next !== nothing iter.fun(next[1]) end
    return next
end

tee(iter::I, fun::F) where {I, F} = TeeIterable{I, F}(iter, fun)

#Sampling
struct SamplingIterable{I}
    iter::I
    period::UInt
end

function iterate(iter::SamplingIterable, state=iter.iter)
    current = iterate(state)
    if current === nothing return nothing end
    for i = 1:iter.period-1
        next = iterate(state, current[2])
        if next === nothing return current[1], rest(state, current[2]) end
        current = next
    end
    return current[1], rest(state, current[2])
end

sample(iter::I, period) where I = SamplingIterable{I}(iter, period)

#Timing
struct StopwatchIterable{I}
    iter::I
end

function iterate(iter::StopwatchIterable)
    t0 = time_ns()
    next = iterate(iter.iter)
    return dispatch(iter, t0, next)
end

function iterate(iter::StopwatchIterable, (t0, state))
    next = iterate(iter.iter, state)
    return dispatch(iter, t0, next)
end

function dispatch(iter::StopwatchIterable, t0, next)
    if next === nothing return nothing end
    return (time_ns()-t0, next[1]), (t0, next[2])
end

stopwatch(iter::I) where I = StopwatchIterable{I}(iter)

#loop
function loop(iter)
    x = nothing
    for y in iter x = y end
    return x
end


text:
using WCSRO
using Documenter

makedocs(;
    modules=[WCSRO],
    authors="Ulrich Pototschnig <ulrich.pototschnig@hotmail.com> and contributors",
    repo="https://github.com/ulipoto/WCSRO.jl/blob/{commit}{path}#L{line}",
    sitename="WCSRO.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://ulipoto.github.io/WCSRO.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/ulipoto/WCSRO.jl",
)

max_stars_repo_name: ulipoto/WCSRO

max_issues_repo_name: ulipoto/WCSRO

max_forks_repo_name: ulipoto/WCSRO