text:
abstract type AbstractCallback end

struct Default <: AbstractCallback end
struct Info <: AbstractCallback end
struct CFL <: AbstractCallback end
struct StateCheck{T} <: AbstractCallback
    number_of_checks::T
end

Base.@kwdef struct JLD2State{T, V, B} <: AbstractCallback
    iteration::T
    filepath::V
    overwrite::B = true
end

Base.@kwdef struct VTKState{T, V, C, B} <: AbstractCallback
    iteration::T = 1
    filepath::V = "."
    counter::C = [0]
    overwrite::B = true
end

function create_callbacks(simulation::Simulation, odesolver)
    callbacks = simulation.callbacks

    if isempty(callbacks)
        return ()
    else
        cbvector = [create_callback(callback, simulation, odesolver) for callback in callbacks]
        return tuple(cbvector...)
    end
end

function create_callback(::Default, simulation::Simulation, odesolver)
    cb_info = create_callback(Info(), simulation, odesolver)
    cb_state_check = create_callback(StateCheck(10), simulation, odesolver)

    return (cb_info, cb_state_check)
end

function create_callback(::Info, simulation::Simulation, odesolver)
    Q = simulation.state
    timeend = simulation.time.finish
    mpicomm = MPI.COMM_WORLD

    starttime = Ref(now())
    cbinfo = ClimateMachine.GenericCallbacks.EveryXWallTimeSeconds(60, mpicomm) do (s = false)
        if s
            starttime[] = now()
        else
            energy = norm(Q)
            @info @sprintf(
                """Update
                simtime = %8.4f / %8.4f
                runtime = %s
                norm(Q) = %.16e""",
                ClimateMachine.ODESolvers.gettime(odesolver),
                timeend,
                Dates.format(convert(Dates.DateTime, Dates.now() - starttime[]), Dates.dateformat"HH:MM:SS"),
                energy
            )

            if isnan(energy)
                error("NaNs")
            end
        end
    end

    return cbinfo
end

function create_callback(::CFL, simulation::Simulation, odesolver)
    Q = simulation.state
    # timeend = simulation.time.finish
    # mpicomm = MPI.COMM_WORLD
    # starttime = Ref(now())
    cbcfl = EveryXSimulationSteps(100) do
        simtime = gettime(odesolver)

        @views begin
            ρ = Array(Q.data[:, 1, :])
            ρu = Array(Q.data[:, 2, :])
            ρv = Array(Q.data[:, 3, :])
            ρw = Array(Q.data[:, 4, :])
        end

        u = ρu ./ ρ
        v = ρv ./ ρ
        w = ρw ./ ρ

        # TODO! transform onto sphere

        ue = extrema(u)
        ve = extrema(v)
        we = extrema(w)

        @info @sprintf """CFL
                simtime = %.16e
                u = (%.4e, %.4e)
                v = (%.4e, %.4e)
                w = (%.4e, %.4e)
                """ simtime ue... ve... we...
    end

    return cbcfl
end

function create_callback(callback::StateCheck, simulation::Simulation, _...)
    sim_length = simulation.time.finish - simulation.time.start
    timestep = simulation.timestepper.timestep
    nChecks = callback.number_of_checks

    nt_freq = floor(Int, sim_length / timestep / nChecks)

    cbcs_dg = ClimateMachine.StateCheck.sccreate([(simulation.state, "state")], nt_freq)

    return cbcs_dg
end

function create_callback(output::JLD2State, simulation::Simulation, odesolver)
    # Initialize output
    output.overwrite && isfile(output.filepath) && rm(output.filepath; force = output.overwrite)

    Q = simulation.state
    mpicomm = MPI.COMM_WORLD
    iteration = output.iteration

    steps = ClimateMachine.ODESolvers.getsteps(odesolver)
    time = ClimateMachine.ODESolvers.gettime(odesolver)

    file = jldopen(output.filepath, "a+")
    JLD2.Group(file, "state")
    JLD2.Group(file, "time")
    file["state"][string(steps)] = Array(Q)
    file["time"][string(steps)] = time
    close(file)

    jldcallback = ClimateMachine.GenericCallbacks.EveryXSimulationSteps(iteration) do (s = false)
        steps = ClimateMachine.ODESolvers.getsteps(odesolver)
        time = ClimateMachine.ODESolvers.gettime(odesolver)
        @info steps, time
        file = jldopen(output.filepath, "a+")
        file["state"][string(steps)] = Array(Q)
        file["time"][string(steps)] = time
        close(file)
        return nothing
    end

    return jldcallback
end

function create_callback(output::VTKState, simulation::Simulation, odesolver)
    # Initialize output
    output.overwrite && isfile(output.filepath) && rm(output.filepath; force = output.overwrite)
    mkpath(output.filepath)

    state = simulation.state
    model = (simulation.rhs isa Tuple) ? simulation.rhs[1] : simulation.rhs

    function do_output(counter, model, state)
        mpicomm = MPI.COMM_WORLD
        balance_law = model.balance_law
        aux_state = model.state_auxiliary

        outprefix = @sprintf("%s/mpirank%04d_step%04d", output.filepath, MPI.Comm_rank(mpicomm), counter[1],)

        @info "doing VTK output" outprefix

        state_names = flattenednames(vars_state(balance_law, Prognostic(), eltype(state)))
        aux_names = flattenednames(vars_state(balance_law, Auxiliary(), eltype(state)))

        writevtk(outprefix, state, model, state_names, aux_state, aux_names)

        counter[1] += 1

        return nothing
    end

    do_output(output.counter, model, state)
    cbvtk = ClimateMachine.GenericCallbacks.EveryXSimulationSteps(output.iteration) do (init = false)
        do_output(output.counter, model, state)
        return nothing
    end

    return cbvtk
end

max_stars_repo_name: CliMA/CouplerMachine

max_issues_repo_name: CliMA/CouplerMachine


text:
module VersionsJSONUtil

using HTTP, JSON, Pkg.BinaryPlatforms, WebCacheUtilities, SHA, Lazy
import Pkg.BinaryPlatforms: triplet, arch

"Wrapper type to define two jlext methods for portable and installer Windows"
struct PortableWindows
    windows::Windows
end
PortableWindows(arch::Symbol) = PortableWindows(Windows(arch))
@forward PortableWindows.windows (up_os, tar_os, triplet, arch)

up_os(p::Windows) = "winnt"
up_os(p::MacOS) = "mac"
up_os(p::Linux) = libc(p) == :glibc ? "linux" : "musl"
up_os(p::FreeBSD) = "freebsd"
up_os(p) = error("Unknown OS for $(p)")

up_arch(p) = up_arch(arch(p))
function up_arch(arch::Symbol)
    if arch == :x86_64
        return "x64"
    elseif arch == :i686
        return "x86"
    elseif arch == :powerpc64le
        return "ppc64le"
    else
        return string(arch)
    end
end

tar_os(p::Windows) = "win$(wordsize(p))"
tar_os(p::FreeBSD) = "freebsd-$(arch(p))"
function tar_os(p::MacOS)
    if arch(p) == :aarch64
        return "macaarch$(wordsize(p))"
    else
        return "mac$(wordsize(p))"
    end
end
function tar_os(p::Linux)
    if arch(p) == :powerpc64le
        return "$(up_os(p))-ppc64le"
    else
        return "$(up_os(p))-$(arch(p))"
    end
end

jlext(p::Windows) = "exe"
jlext(p::PortableWindows) = "zip"
jlext(p::MacOS) = "dmg"
jlext(p) = "tar.gz"

# OS to use in the metadata
# The OS in the download URL for Linux with musl is "musl"
# But the OS in the metadata should be "linux"
meta_os(p) = up_os(p)
meta_os(p::Linux) = "linux"

function download_url(version::VersionNumber, platform)
    return string(
        "https://julialang-s3.julialang.org/bin/",
        up_os(platform), "/",
        up_arch(platform), "/",
        version.major, ".", version.minor, "/",
        "julia-", version, "-", tar_os(platform), ".", jlext(platform),
    )
end

# We're going to collect the combinatorial explosion of version/os-arch possible downloads.
# We don't have a nice, neat list of what is or is not available, and so we're just going to
# try and download each file, and if it exists, yay.  Otherwise, bleh.
julia_platforms = [
    Linux(:x86_64),
    Linux(:i686),
    Linux(:aarch64),
    Linux(:armv7l),
    Linux(:powerpc64le),
    Linux(:x86_64, libc = :musl),
    MacOS(:x86_64),
    MacOS(:aarch64),
    Windows(:x86_64),
    Windows(:i686),
    PortableWindows(:x86_64),
    PortableWindows(:i686),
    FreeBSD(:x86_64),
]

function vnum_maybe(x::AbstractString)
    try
        return VersionNumber(x)
    catch
        return nothing
    end
end

function is_stable(v::VersionNumber)
    return v.prerelease == () && v.build == ()
end

# Get list of tags from the Julia repo
function get_tags()
    @info("Probing for tag list...")
    tags_json_path = WebCacheUtilities.download_to_cache(
        "julia_tags.json",
        "https://api.github.com/repos/JuliaLang/julia/git/refs/tags",
    )
    JSON.parse(String(read(tags_json_path)))
end

function main(out_path)
    tags = get_tags()
    tag_versions = filter(x -> x !== nothing, [vnum_maybe(basename(t["ref"])) for t in tags])

    meta = Dict()
    for version in tag_versions
        for platform in julia_platforms
            url = download_url(version, platform)
            filename = basename(url)

            # Download this URL to a local file
            local filepath
            try
                @info("Downloading $(filename)...")
                filepath = WebCacheUtilities.download_to_cache(filename, url)
            catch e
                if isa(e, InterruptException)
                    rethrow(e)
                end
                continue
            end

            tarball_hash_path = hit_file_cache("$(filename).sha256") do tarball_hash_path
                open(filepath, "r") do io
                    open(tarball_hash_path, "w") do hash_io
                        write(hash_io, bytes2hex(sha256(io)))
                    end
                end
            end
            tarball_hash = String(read(tarball_hash_path))

            # Initialize overall version key, if needed
            if !haskey(meta, version)
                meta[version] = Dict(
                    "stable" => is_stable(version),
                    "files" => Vector{Dict}(),
                )
            end

            # Test to see if there is an asc signature:
            asc_signature = nothing
            if !isa(platform, MacOS) && !isa(platform, Windows)
                try
                    asc_url = string(url, ".asc")
                    @info("Downloading $(basename(asc_url))")
                    asc_filepath = WebCacheUtilities.download_to_cache(basename(asc_url), asc_url)
                    asc_signature = String(read(asc_filepath))
                catch e
                    if isa(e, InterruptException)
                        rethrow(e)
                    end
                end
            end

            # Build up metadata about this file
            kind = "archive"
            if endswith(filename, ".exe")
                kind = "installer"
            end
            file_dict = Dict(
                "triplet" => triplet(platform),
                "os" => meta_os(platform),
                "arch" => string(arch(platform)),
                "version" => string(version),
                "sha256" => tarball_hash,
                "size" => filesize(filepath),
                "kind" => kind,
                "url" => url,
            )
            # Add in `.asc` signature content, if applicable
            if asc_signature !== nothing
                file_dict["asc"] = asc_signature
            end

            # Right now, all we have are archives, but let's be forward-thinking
            # and make this an array of dictionaries that is easy to extensibly match
            push!(meta[version]["files"], file_dict)

            # Write out new versions of our versions.json as we go
            open(out_path, "w") do io
                JSON.print(io, meta, 2)
            end

            # Delete downloaded file
            rm(filepath)
        end
    end
end

end # module

ext: jl


text:
mutable struct HashDictionary{I,T} <: AbstractDictionary{I, T}
    indices::HashIndices{I}
    values::Vector{T}

    HashDictionary{I, T}(indices::HashIndices{I}, values::Vector{T}, ::Nothing) where {I, T} = new(indices, values)
end

"""
    HashDictionary{I, T}()

Construct an empty `HashDictionary` with index type `I` and element type `T`. This type of
dictionary uses hashes for fast lookup and insertion, and is both mutable and insertable.
(See `issettable` and `isinsertable`).
"""
function HashDictionary{I, T}(; sizehint::Int = 16) where {I, T}
    indices = HashIndices{I}(; sizehint=sizehint)
    HashDictionary{I, T}(indices, Vector{T}(undef, length(indices.slots)), nothing)
end
HashDictionary{I}() where {I} = HashDictionary{I, Any}()
HashDictionary() = HashDictionary{Any}()

"""
    HashDictionary{I, T}(indices, undef::UndefInitializer)

Construct a `HashDictionary` with index type `I` and element type `T`. The container is
initialized with `keys` that match the values of `indices`, but the values are unintialized.
"""
function HashDictionary{I, T}(indices, ::UndefInitializer) where {I, T}
    return HashDictionary{I, T}(HashIndices{I}(indices), undef)
end

function HashDictionary{I, T}(h::HashIndices{I}, ::UndefInitializer) where {I, T}
    return HashDictionary{I, T}(h, Vector{T}(undef, length(h.slots)), nothing)
end

function HashDictionary{I, T}(indices::HashIndices{I}, values) where {I, T}
    vals = Vector{T}(undef, length(indices.slots))
    d = HashDictionary{I, T}(indices, vals, nothing)

    @inbounds for (i, v) in zip(tokens(indices), values)
        vals[i] = v
    end

    return d
end

"""
    HashDictionary(indices, values)
    HashDictionary{I}(indices, values)
    HashDictionary{I, T}(indices, values)

Construct a `HashDictionary` with indices from `indices` and values from `values`, matched
in iteration order.
"""
function HashDictionary{I, T}(indices, values) where {I, T}
    iter_size = Base.IteratorSize(indices)
    if iter_size isa Union{Base.HasLength, Base.HasShape}
        d = HashDictionary{I, T}(; sizehint = length(indices)*2)
    else
        d = HashDictionary{I, T}()
    end

    for (i, v) in zip(indices, values)
        insert!(d, i, v)
    end

    return d
end
function HashDictionary{I}(indices, values) where {I}
    if Base.IteratorEltype(values) === Base.EltypeUnknown()
        # TODO: implement automatic widening from iterators of Base.EltypeUnkown
        values = collect(values)
    end

    return HashDictionary{I, eltype(values)}(indices, values)
end

function HashDictionary(indices, values)
    if Base.IteratorEltype(indices) === Base.EltypeUnknown()
        # TODO: implement automatic widening from iterators of Base.EltypeUnkown
        indices = collect(indices)
    end

    return HashDictionary{eltype(indices)}(indices, values)
end

"""
    HashDictionary(dict::AbstractDictionary)
    HashDictionary{I}(dict::AbstractDictionary)
    HashDictionary{I, T}(dict::AbstractDictionary)

Construct a copy of `dict` with the same keys and values.

(For copying an `AbstractDict` or other iterable of `Pair`s, see `dictionary`).
"""
HashDictionary(dict::AbstractDictionary) = HashDictionary(keys(dict), dict)
HashDictionary{I}(dict::AbstractDictionary) where {I} = HashDictionary{I}(keys(dict), dict)
HashDictionary{I, T}(dict::AbstractDictionary) where {I, T} = HashDictionary{I, T}(keys(dict), dict)

"""
    dictionary(iter)

Construct a new `AbstractDictionary` from an iterable `iter` of key-value `Pair`s. The
default container type is `HashDictionary`.
"""
function dictionary(iter)
    if Base.IteratorEltype(iter) === Base.EltypeUnknown()
        # TODO: implement automatic widening from iterators of Base.EltypeUnkown
        iter = collect(iter)
    end
    _dictionary(eltype(iter), iter)
end

dictionary(p1::Pair, p2::Pair...) = dictionary((p1, p2...))

function _dictionary(::Type{Pair{I, T}}, iter) where {I, T}
    iter_size = Base.IteratorSize(iter)
    if iter_size isa Union{Base.HasLength, Base.HasShape}
        d = HashDictionary{I, T}(; sizehint = length(iter)*2)
    else
        d = HashDictionary{I, T}()
    end

    for (i, v) in iter
        insert!(d, i, v)
    end

    return d
end

## Implementation

Base.keys(d::HashDictionary) = d.indices
isinsertable(d::HashDictionary) = true
issettable(d::HashDictionary) = true

@propagate_inbounds function gettoken(d::HashDictionary{I}, i::I) where {I}
    return gettoken(keys(d), i)
end

@inline function gettokenvalue(d::HashDictionary, token)
    return @inbounds d.values[token]
end

function istokenassigned(d::HashDictionary, token)
    return isassigned(d.values, token)
end

@inline function settokenvalue!(d::HashDictionary{I, T}, token, value::T) where {I, T}
    @inbounds d.values[token] = value
    return d
end

function gettoken!(d::HashDictionary{T}, key::T) where {T}
    indices = keys(d)
    (token, values) = _gettoken!(indices, d.values, key)
    if token < 0
        (token, values) = _insert!(indices, values, key, -token)
        d.values = values
        return (false, token)
    else
        d.values = values
        return (true, token)
    end
end

function Base.copy(d::HashDictionary{I, T}) where {I, T}
    return HashDictionary{I, T}(d.indices, copy(d.values), nothing)
end

tokenized(d::HashDictionary) = d.values

function Base.empty!(d::HashDictionary)
    empty!(d.indices)
    empty!(d.values)
    resize!(d.values, length(keys(d).slots))
    return d
end

function deletetoken!(d::HashDictionary{I, T}, token) where {I, T}
    deletetoken!(keys(d), token)
    isbitstype(T) || ccall(:jl_arrayunset, Cvoid, (Any, UInt), d.values, token-1)
    return d
end

function Base.sizehint!(d::HashDictionary, sz::Int)
    d.values = _sizehint!(d.indices, d.values, sz)
    return d
end

function Base.rehash!(d::HashDictionary, newsz::Int = length(d.inds))
    _rehash!(d.indices, d.values, newsz)
    return d
end

Base.filter!(pred, d::HashDictionary) = Base.unsafe_filter!(pred, d)

# `HashDictionary` is the default mutable AbstractDictionary
# If given some other index type, create a new `HashIndices` copy of the indices
function Base.similar(indices::AbstractIndices{I}, ::Type{T}) where {I, T}
    return similar(HashIndices{I}(indices), T)
end

# For `HashIndices` we don't copy the indices, we allow the `keys` to remain identical (`===`)
function Base.similar(indices::HashIndices{I}, ::Type{T}) where {I, T}
    return HashDictionary{I, T}(indices, undef)
end

# `HashDictionary` is the default insertable AbstractDictionary
function Base.empty(::AbstractDictionary, ::Type{I}, ::Type{T}) where {I, T}
    return HashDictionary{I, T}()
end


text:
setupfolders() =
    mkpath.([tablefolder, plotfolder])

mutable struct PlotSaver
    name
    button
end

PlotSaver(btn = "Save figure as"; label = "name") =
    PlotSaver(textbox(label), button(btn))

observe(p::PlotSaver) = observe(p.button)

name(p::PlotSaver) = observe(p.name).val
filename(p::PlotSaver) = joinpath(plotfolder, name(p))

mutable struct TableSaver
    name
    button
    checkbox
end

TableSaver(btn = "Save table as"; label = "name", checkbox = "select data") =
    TableSaver(textbox(label), button(btn), InteractBase.checkbox(true, label = checkbox))

observe(p::TableSaver) = observe(p.button)

name(p::TableSaver) = observe(p.name).val
filename(p::TableSaver) = joinpath(tablefolder, name(p))
isselected(t::TableSaver) = observe(t.checkbox).val

function _save(data, style, fn, sel)
    t1 = sel ? build_table(data, style) : data[]
    JuliaDB.save(t1, fn)
end


text:
julia> i = 0
0

julia> while true
           println(i)
           i += 1
           i % 6 == 0 && break
       end
0
1
2
3
4
5

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# coding: utf-8

# In[1]:

#a simple day trading game
#day trader is only allowed to make at maximum two trades
#the strategy is long only
#lets find out the maximum profit

#more details can be found in the following link
# https://www.geeksforgeeks.org/maximum-profit-by-buying-and-selling-a-share-at-most-twice/

#an alternative version in dynamic programming exists
#its done by using a different approach
#strongly recommend you to take a look
# https://github.com/je-suis-tm/recursion-and-dynamic-programming/blob/master/stock%20trading%20dynamic%20programming.jl

# In[2]:

#compute the maximum profit for long only strategy
function compute_profit(prices)

    #initialize maximum price with the close price
    max_price=prices[end]

    #initialize minimum price with the open price
    min_price=prices[1]

    #initialize indices
    left=1
    right=length(prices)
    minind=1
    maxind=length(prices)

    #we have two indices moving at the same time
    #one from left to find the minimum value
    #the other from right to find the maximum value
    #we are not looking for global minimum
    #we only want local minimum before the global maximum
    while left<maxind && right>minind

        #when a larger value is found
        #update accordingly
        if prices[right]>max_price

            max_price=prices[right]
            maxind=copy(right)

        end

        #the same applies to a smaller value
        if prices[left]<min_price

            min_price=prices[left]
            minind=copy(left)

        end

        left+=1
        right-=1

    end

    #maximum profit
    profit=max_price-min_price

    #when we lose money
    #abort the trade
    if profit<0

        profit=0

    end

    return profit

end

# In[3]:

#there are two scenarios to maximize the profit
#one trade or two trades
#since we can execute two transactions
#we split the array at an iterating index i into half
#we find out the maximum profit we can obtain from both halves
function stock_trading(prices,ind)

    #prevent index error
    if ind+1==length(prices)

        return 0

    end

    #split
    upper=prices[1:ind]
    lower=prices[ind+1:end]

    #compute profit
    profit=compute_profit(lower)+compute_profit(upper)

    #compare recursively
    return max(profit,stock_trading(prices,ind+1))

end

# In[4]:

stock_trading([10,22,5,75,65,80],1)

# In[5]:

stock_trading([2,30,15,10,8,25,80],1)

# In[6]:

stock_trading([100,30,15,10,8,25,80],1)

# In[7]:

stock_trading([90,70,35,11,5],1)

max_stars_repo_path: stock trading recursion.jl

max_issues_repo_path: stock trading recursion.jl

max_forks_repo_path: stock trading recursion.jl


text:
@testset "GCN" begin

include("definitions/gcn.jl")

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

@testset "IR" begin

@testset "kernel calling convention" begin
    kernel() = return

    ir = sprint(io->gcn_code_llvm(io, kernel, Tuple{}; dump_module=true))
    @test !occursin("amdgpu_kernel", ir)

    ir = sprint(io->gcn_code_llvm(io, kernel, Tuple{};
                                         dump_module=true, kernel=true))
    @test occursin("amdgpu_kernel", ir)
end

@testset "alloca addrspace" begin
    function kernel(i)
        sink(i) # sink provides an alloca in addrspace 0
        return
    end

    ir = sprint(io->gcn_code_llvm(io, kernel, Tuple{Int64}; dump_module=true))
    @test occursin(r"alloca i64, (align 8, )?addrspace\(5\)$"m, ir)
    @test !occursin(r"alloca i64(, align \d)?$"m, ir)
end

end

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

@testset "assembly" begin

@testset "skip scalar trap" begin
    workitem_idx_x() = ccall("llvm.amdgcn.workitem.id.x", llvmcall, Int32, ())
    trap() = ccall("llvm.trap", llvmcall, Nothing, ())
    function kernel()
        if workitem_idx_x() > 1
            trap()
        end
        return
    end

    asm = sprint(io->gcn_code_native(io, kernel, Tuple{}))
    @test occursin("s_trap 2", asm)
    @test_skip occursin("s_cbranch_execz", asm)
    if Base.libllvm_version < v"9"
        @test_broken occursin("v_readfirstlane", asm)
    end
end

@testset "child functions" begin
    # we often test using @noinline child functions, so test whether these survive
    # (despite not having side-effects)
    @noinline child(i) = sink(i)
    function parent(i)
        child(i)
        return
    end

    asm = sprint(io->gcn_code_native(io, parent, Tuple{Int64}; dump_module=true))
    @test occursin(r"s_add_u32.*julia_child_.*@rel32@lo\+4", asm)
    @test occursin(r"s_addc_u32.*julia_child_.*@rel32@hi\+4", asm)
end

@testset "kernel functions" begin
    @noinline nonentry(i) = sink(i)
    function entry(i)
        nonentry(i)
        return
    end

    asm = sprint(io->gcn_code_native(io, entry, Tuple{Int64}; dump_module=true, kernel=true))
    @test occursin(r"\.amdhsa_kernel _Z\d*julia_entry", asm)
    @test !occursin(r"\.amdhsa_kernel _Z\d*julia_nonentry", asm)
    @test occursin(r"\.type.*julia_nonentry_\d*,@function", asm)
end

@testset "child function reuse" begin
    # bug: depending on a child function from multiple parents resulted in
    #      the child only being present once

    @noinline child(i) = sink(i)
    function parent1(i)
        child(i)
        return
    end

    asm = sprint(io->gcn_code_native(io, parent1, Tuple{Int}; dump_module=true))
    @test occursin(r"\.type.*julia_[[:alnum:]_.]*child_\d*,@function", asm)

    function parent2(i)
        child(i+1)
        return
    end

    asm = sprint(io->gcn_code_native(io, parent2, Tuple{Int}; dump_module=true))
    @test occursin(r"\.type.*julia_[[:alnum:]_.]*child_\d*,@function", asm)
end

@testset "child function reuse bis" begin
    # bug: similar, but slightly different issue as above
    #      in the case of two child functions
    @noinline child1(i) = sink(i)
    @noinline child2(i) = sink(i+1)
    function parent1(i)
        child1(i) + child2(i)
        return
    end
    gcn_code_native(devnull, parent1, Tuple{Int})

    function parent2(i)
        child1(i+1) + child2(i+1)
        return
    end
    gcn_code_native(devnull, parent2, Tuple{Int})
end

@testset "indirect sysimg function use" begin
    # issue #9: re-using sysimg functions should force recompilation
    #           (host fldmod1->mod1 throws, so the GCN code shouldn't contain a throw)

    # NOTE: Int32 to test for #49

    function kernel(out)
        wid, lane = fldmod1(unsafe_load(out), Int32(32))
        unsafe_store!(out, wid)
        return
    end

    asm = sprint(io->gcn_code_native(io, kernel, Tuple{Ptr{Int32}}))
    @test !occursin("jl_throw", asm)
    @test !occursin("jl_invoke", asm)   # forced recompilation should still not invoke
end

@testset "LLVM intrinsics" begin
    # issue #13 (a): cannot select trunc
    function kernel(x)
        unsafe_trunc(Int, x)
        return
    end
    gcn_code_native(devnull, kernel, Tuple{Float64})
end

@test_broken "exception arguments"
#= FIXME: _ZNK4llvm14TargetLowering20scalarizeVectorStoreEPNS_11StoreSDNodeERNS_12SelectionDAGE
@testset "exception arguments" begin
    function kernel(a)
        unsafe_store!(a, trunc(Int, unsafe_load(a)))
        return
    end

    gcn_code_native(devnull, kernel, Tuple{Ptr{Float64}})
end
=#

@test_broken "GC and TLS lowering"
#= FIXME: in function julia_inner_18528 void (%jl_value_t addrspace(10)*): invalid addrspacecast
@testset "GC and TLS lowering" begin
    @eval mutable struct PleaseAllocate
        y::Csize_t
    end

    # common pattern in Julia 0.7: outlined throw to avoid a GC frame in the calling code
    @noinline function inner(x)
        sink(x.y)
        nothing
    end

    function kernel(i)
        inner(PleaseAllocate(Csize_t(42)))
        nothing
    end

    asm = sprint(io->gcn_code_native(io, kernel, Tuple{Int}))
    @test occursin("gpu_gc_pool_alloc", asm)

    # make sure that we can still ellide allocations
    function ref_kernel(ptr, i)
        data = Ref{Int64}()
        data[] = 0
        if i > 1
            data[] = 1
        else
            data[] = 2
        end
        unsafe_store!(ptr, data[], i)
        return nothing
    end

    asm = sprint(io->gcn_code_native(io, ref_kernel, Tuple{Ptr{Int64}, Int}))

    @test !occursin("gpu_gc_pool_alloc", asm)
end
=#

@testset "float boxes" begin
    function kernel(a,b)
        c = Int32(a)
        # the conversion to Int32 may fail, in which case the input Float32 is boxed in order to
        # pass it to the @nospecialize exception constructor. we should really avoid that (eg.
        # by avoiding @nospecialize, or optimize the unused arguments away), but for now the box
        # should just work.
        unsafe_store!(b, c)
        return
    end

    ir = sprint(io->gcn_code_llvm(io, kernel, Tuple{Float32,Ptr{Float32}}))
    @test occursin("jl_box_float32", ir)
    gcn_code_native(devnull, kernel, Tuple{Float32,Ptr{Float32}})
end

end

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

end


text:
include("gitutils.jl")
include("github.jl")
include("dash.jl")
include("components.jl")
include("deploy.jl")

const HELP = """
Usage: generate.jl [--deploy] [--help]

Options:

--deploy   Deploy resources to repo specified in `Sources.toml`
--help     Print out this message.
"""

function generate(ARGS, sources, build_dir, artifact_file)
    if "--help" in ARGS
        println(HELP)
        return nothing
    end

    is_deploy =  "--deploy" in ARGS

    deploy_dir = joinpath(build_dir, "deploy")

    rm(build_dir, force = true, recursive = true)
    mkdir(build_dir)

    cd(build_dir)
    @info "cloning dash..."
    v = install_dash(sources["dash"]["url"], sources["dash"]["tag"])
    dash_version = VersionNumber(v)
    @info "dash succefully installed" dash_version

    deploy_repo = sources["deploy"]["repo"]

    init_deploy_repo(deploy_repo, deploy_dir)

    old_info = old_build_info(deploy_dir)
    new_info = make_new_info(dash_version, old_info)
    @info "filling resources data..."
    fill_deploy_resources(sources, deploy_dir, new_info)
    @info "resources data filled"

    artifact_hash = create_artifact() do dir
        res_dir = resources_dir(deploy_dir)
        for name in readdir(res_dir)
            cp(
                joinpath(res_dir, name),
                joinpath(dir, name)
            )
        end
    end
    @info "artifact created" artifact_hash
    if is_deploy
        @info "pushing repo to https://github.com/$(deploy_repo)..."
        push_repo(deploy_repo, deploy_dir, deploy_tagname(new_info))
        @info "repo pushed"

        @info "creating tagball..."
        tarball_hash = archive_artifact(artifact_hash, joinpath(build_dir, tarball_name(new_info)))
        @info "tagball created" tarball_hash

        @info "binding artifact..."
        download_info = (
            "https://github.com/$(deploy_repo)/releases/download/$(deploy_tagname(new_info))/$(tarball_name(new_info))",
            tarball_hash
        )
        bind_artifact!(artifact_file, "dash_resources",  artifact_hash, force = true, download_info = [download_info])
        @info "artifact binding succefully to" artifact_file

        cd(deploy_dir) do
            upload_to_releases(deploy_repo,
                deploy_tagname(new_info),
                realpath(joinpath(build_dir, tarball_name(new_info)))
            )
        end
    else
        unbind_artifact!(artifact_file, "dash_resources")
        bind_artifact!(artifact_file, "dash_resources",  artifact_hash, force = true, download_info = nothing)
        @info "artifact binding succefully to" artifact_file
    end
    @info "resource generation done!"

end


text:
function run_clust_dbaclust_centroid(
      region::String,
      opt_problems::Array{String},
      norm_op::String,
      norm_scope::String,
      n_clust_ar::Array,
      n_init::Int,
      iterations::Int;
      rad_sc_min::Int=0,
      rad_sc_max::Int=3,
      inner_iterations::Int=30
    )

    n_clust_min = minimum(n_clust_ar)
    n_clust_max = maximum(n_clust_ar)

    # read in original data
    seq = load_pricedata(region)

    # initial points
    n_dbaclust= n_init # number of dbaclust runs (each with n_init below)

    n_init = 1 # number of initial guesses for each dbaclust run / should be set to 1 for each experiment. dbaclust is deterministic once an initial guess has been formed. We sample the initial guesses on the outer run.

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

    # create directory where data is saved
    try
      mkdir("outfiles")
    catch
     #
    end

    # save settings in txt file
    df = DataFrame()
    df[:n_clust_min]=n_clust_min
    df[:n_clust_max]=n_clust_max
    df[:n_init]=n_init
    df[:n_dbaclust]=n_dbaclust
    df[:rad_sc_min]=rad_sc_min
    df[:rad_sc_max]=rad_sc_max
    df[:iterations]=iterations
    df[:inner_iterations]=inner_iterations
    df[:region]=region

    writetable(joinpath("outfiles",string("parameters.txt")),df)

     # iterate through settings
    for n_clust=n_clust_min:n_clust_max
      for rad_sc=rad_sc_min:rad_sc_max
        for i = 1:n_dbaclust

          rmin,rmax = sakoe_chiba_band(rad_sc,24)

           ##########################
          # normalized clustering hourly
          seq_norm, hourly_mean, hourly_sdv = z_normalize(seq,scope="sequence")
          tic()
          results = dbaclust(seq_norm,n_clust,n_init,ClassicDTW();iterations=iterations,inner_iterations=inner_iterations,rtol=1e-5,show_progress=false,store_trace=false,i2min=rmin,i2max=rmax)
          el_time = toq()
          println("Elapsed time: ",el_time ," ; n_clust=",n_clust," rad_sc=",rad_sc," i=",i)
          flush(STDOUT)

          centers_norm = results.centers
          clustids = results.clustids
          centers = undo_z_normalize(seq_to_array(centers_norm),hourly_mean,hourly_sdv;idx=clustids)

           # save results to txt

          writetable(joinpath("outfiles",string("dbaclust_k_",n_clust,"_scband_",rad_sc,"_ninit_",n_init,"_it_",iterations,"_innerit_",inner_iterations,"_",i,"_cluster.txt")),DataFrame(centers'),separator='\t',header=false)
          writetable(joinpath("outfiles",string("dbaclust_k_",n_clust,"_scband_",rad_sc,"_ninit_",n_init,"_it_",iterations,"_innerit_",inner_iterations,"_",i,"_clustids.txt")),DataFrame(id=clustids),separator='\t',header=false)
          writetable(joinpath("outfiles",string("dbaclust_k_",n_clust,"_scband_",rad_sc,"_ninit_",n_init,"_it_",iterations,"_innerit_",inner_iterations,"_",i,"_cost.txt")),DataFrame(cost=results.dbaresult.cost),separator='\t',header=false)
          writetable(joinpath("outfiles",string("dbaclust_k_",n_clust,"_scband_",rad_sc,"_ninit_",n_init,"_it_",iterations,"_innerit_",inner_iterations,"_",i,"_it.txt")),DataFrame(iterations=results.iterations),separator='\t',header=false)
          writetable(joinpath("outfiles",string("dbaclust_k_",n_clust,"_scband_",rad_sc,"_ninit_",n_init,"_it_",iterations,"_innerit_",inner_iterations,"_",i,"_innerit.txt")),DataFrame(inner_iterations=results.dbaresult.iterations),separator='\t',header=false)

        end
      end
    end

    # TODO  - import dbaclust_res_to_jld2.jl and automatically generate jld2 files.
    return 0
end # function


text:
module TestProgressMeter

export insertProgress, removeProgress

@doc raw"
    insertProgress(f = 'runtests.jl'; toplevel = true, ntests = nothing, loops = true, pmargs = ())
recursively inserts `next!` statements before each test that's included in `f`
and initializes the `ProgressMeter` at the beginning of the file according to the
following rules:
- if `ntests` - the number of tests - is provided, a progressbar will be shown with equal weight to each of the `ntests` tests
- if `ntests` is not provided but `loops=false`, all `@test` will be counted and a progressbar will be shown with equal weight to each test
- if neither `ntests` is provided nor `loops=true`, a simple counter will be shown of how many tests have been run

via `pmargs`, keyword-arguments to the `ProgressMeter`-object can be provided such as e.g.
```julia
julia> insertProgress(\"example/runtest.jl\", pmargs = (desc = \"Text next to errorbar\", dt = 0.1), ntests = 10)
```
to get a progressbar that updates up to ever `0.1` seconds and has the text in `desc` next to it.
"
function insertProgress(f = "runtests.jl"; toplevel = true,
                                            ntests = nothing,
                                            loops = true,
                                            pmargs = ())
    counter = 0
    isfile(f) || throw(ArgumentError("$f is not a file"))

    lines = readlines(f)

    # if ProgressMeter is already used, assume its from TestProgressMeter and strip
    # all lines that were inserted
    if toplevel && any(x -> occursin("using ProgressMeter", x), lines)
        removeProgress(f)
        insertProgress(f, toplevel = toplevel, ntests = ntests, loops = loops, pmargs = pmargs)
        return nothing
    end

    for i in 1:length(lines)
        m = match(r"^([ \t]*)@test[^s]", lines[i])
        if m != nothing
            lines[i] = string(m.captures[],"@isdefined(pmobj) && next!(pmobj)\n", lines[i])
            counter += 1
        end
        m = match(r"^[^\#]*(include\(\")([^\")]*)(\"\))", lines[i])
        if m != nothing
            nf  = joinpath(dirname(f), m.captures[2])
            counter += insertProgress(nf; toplevel = false)
        end
    end
    if toplevel
        if length(pmargs) == 1
            pmargsstring = chop(string(pmargs), head=1, tail=2)
        else
            pmargsstring = chop(string(pmargs), head=1, tail=1)
        end
        if ntests isa Integer
            counter = ntests
            pushfirst!(lines,
                "pmobj = Progress($counter, $(pmargsstring))")
        elseif loops == false
            pushfirst!(lines,
             "pmobj = Progress($counter, $(pmargsstring))")
        else
            pushfirst!(lines,
             "pmobj = ProgressUnknown($(pmargsstring))")
            lastend = findlast(x -> occursin("end", x), lines)
            insert!(lines, lastend, "@isdefined(pmobj) && finish!(pmobj)")
        end
        pushfirst!(lines, "using ProgressMeter")
        push!(lines, "println()")
    end

    open(f, "w") do io
        write(io, join(lines,'\n'))
    end
    return counter
end

@doc raw"
    removeProgress(f = 'runtests.jl')
assuming `f` is a `runtest` folder, `removeProgress`
removes all statements inserted by `insertProgress`.
"
function removeProgress(f = "runtests.jl")
    isfile(f) || throw(ArgumentError("$f is not a file"))

    lines = readlines(f)
    toremove = falses(length(lines))
    for i in 1:length(lines)
        toremove[i] =   occursin("using ProgressMeter", lines[i]) ||
                        occursin("pmobj = Progress", lines[i]) ||
                        occursin("next!(pmobj)", lines[i]) ||
                        occursin("finish!(pmobj)", lines[i])

        m = match(r"(include\(\")([^\")]*)(\"\))", lines[i])
        if m != nothing
            dir = dirname(f)
            nf  = joinpath(dir, m.captures[2])
            removeProgress(nf)
        end
    end
    lines[end] == "println()" && (toremove[end] = true)

    deleteat!(lines, toremove)
    open(f, "w") do io
        write(io, join(lines,'\n'))
    end
    return nothing
end

end # module


text:
# This file was generated, do not modify it. # hide
G2 = SearchGraph(; dist, db=MatrixDatabase(X))
callbacks2 = SearchGraphCallbacks(hyperparameters=OptimizeParameters(kind=ParetoRadius()))  # it uses distances instead of recall, it will be faster but lower quality
@elapsed index!(G2; callbacks=callbacks2)

max_stars_repo_name: sadit/SimilaritySearchDemos

max_issues_repo_name: sadit/SimilaritySearchDemos

max_forks_repo_name: sadit/SimilaritySearchDemos


text:
module NeuralNetDiffEq

using Reexport, Statistics
@reexport using DiffEqBase

using Flux, Zygote, DiffEqSensitivity, ForwardDiff, Random, Distributions
using DiffEqFlux, Adapt, CuArrays, StochasticDiffEq
import Tracker, Optim

abstract type NeuralNetDiffEqAlgorithm <: DiffEqBase.AbstractODEAlgorithm end

struct TerminalPDEProblem{G,F,Mu,Sigma,X,T,P,A,UD,K} <: DiffEqBase.DEProblem
    g::G
    f::F
    μ::Mu
    σ::Sigma
    X0::X
    tspan::Tuple{T,T}
    p::P
    A::A
    u_domain::UD
    kwargs::K
    TerminalPDEProblem(g,f,μ,σ,X0,tspan,p=nothing;A=nothing,u_domain=nothing,kwargs...) = new{typeof(g),typeof(f),
                                                         typeof(μ),typeof(σ),
                                                         typeof(X0),eltype(tspan),
                                                         typeof(p),typeof(A),typeof(u_domain),typeof(kwargs)}(
                                                         g,f,μ,σ,X0,tspan,p,A,u_domain,kwargs)
end

Base.summary(prob::TerminalPDEProblem) = string(nameof(typeof(prob)))

function Base.show(io::IO, A::TerminalPDEProblem)
  println(io,summary(A))
  print(io,"timespan: ")
  show(io,A.tspan)
end

struct KolmogorovPDEProblem{ F, G, Phi, X , T , D ,P,U0, ND} <: DiffEqBase.DEProblem
    f::F
    g::G
    phi::Phi
    xspan::Tuple{X,X}
    tspan::Tuple{T,T}
    d::D
    p::P
    u0::U0
    noise_rate_prototype::ND
    KolmogorovPDEProblem( f, g, phi , xspan , tspan , d, p=nothing, u0=0 , noise_rate_prototype= nothing) = new{typeof(f),typeof(g),typeof(phi),eltype(tspan),eltype(xspan),typeof(d),typeof(p),typeof(u0),typeof(noise_rate_prototype)}(f,g,phi,xspan,tspan,d,p,u0,noise_rate_prototype)
end

Base.summary(prob::KolmogorovPDEProblem) = string(nameof(typeof(prob)))
function Base.show(io::IO, A::KolmogorovPDEProblem)
  println(io,summary(A))
  print(io,"timespan: ")
  show(io,A.tspan)
  print(io,"xspan: ")
  show(io,A.xspan)
  println(io , "μ")
  show(io , A.f)
  println(io,"Sigma")
  show(io , A.g)
end

include("ode_solve.jl")
include("pde_solve.jl")
include("pde_solve_ns.jl")
include("kolmogorov_solve.jl")
include("stopping_solve.jl")

export NNODE, TerminalPDEProblem, NNPDEHan, NNPDENS, KolmogorovPDEProblem, NNKolmogorov, NNStopping

end # module


text:
using Test
using TimeseriesPrediction
using Statistics, LinearAlgebra
import MultivariateStats
const MS = MultivariateStats

println("Reconstruction Tests")
@testset "cubic_shell_embedding" begin
    for γ=[0,4], τ ∈ [1,20], B=[1,3], k=[1,3], Φ=1:4
        BC=ConstantBoundary(10.)
        @testset "γ=$γ, τ=$τ, B=$B, k=$k, Φ=$Φ" begin
                #Ugly way of creating Φ dim array
                s = [rand(Float64,([10 for i=1:Φ]...,)) for i=1:10]
                emb = cubic_shell_embedding(s,γ,τ,B,k,BC)
                @test emb ==  SpatioTemporalEmbedding(s,(γ=γ,τ=τ,B=B,k=k,bc=BC))
                #Check Embedding Dimension X
                X = (γ+1)*(2B+1)^Φ
                @test typeof(emb) <: SpatioTemporalEmbedding{Φ,ConstantBoundary{Float64},X}
                @test length(emb.τ) == X
                @test length(emb.β) == X

                #check spatial extent
                @test emb.inner == TimeseriesPrediction.Region(
                                ([B*k+1 for i=1:Φ]...,),
                                ([10-B*k for i=1:Φ]...,))

                @test emb.whole == TimeseriesPrediction.Region((ones(Int,Φ)...,),(10ones(Int,Φ)...,) )
        end
    end

    #Checking a few individual Reconstructions
    @testset "Order of rec. points" begin
        γ=1; τ=1; B=1; k=1; c=10; BC=ConstantBoundary(c); Φ=2
        data = [reshape(1+i:9+i, 3,3) for i∈[0,9]]
        emb = cubic_shell_embedding(data, γ,τ,B,k,BC)
        t = 1
        α = CartesianIndex(2,2)
        rvec = zeros(18);
        emb(rvec, data, t, α)
        @test 1:18 == rvec
        α = CartesianIndex(1,1)
        rec = [c,c,c,c,1,2,c,4,5,c,c,c,c,10,11,c,13,14]
        emb(rvec, data, t, α)
        @test rec == rvec
    end
end

@testset "light_cone_embedding 1D" begin
    Φ = 1
    bc=ConstantBoundary(10.)
    s = [rand(Float64,([10 for i=1:Φ]...,)) for i=1:10]
    CI = CartesianIndex

    em = light_cone_embedding(s, 0, 1, 0, 1, bc)
    @test em.β == [CI(0)]

    em = light_cone_embedding(s, 0, 1, 1, 1, bc)
    @test em.β == CI.(-1:1)

    em = light_cone_embedding(s, 0, 1, 2, 1, bc)
    @test em.β == CI.(-2:2)

    em = light_cone_embedding(s, 0, 1, 2, 2, bc)
    @test em.β == CI.(-2:2)

    em = light_cone_embedding(s, 0, 1, 0, 1, bc)
    @test em.β == [CI(0)]

    em = light_cone_embedding(s, 1, 1, 2, 1, bc)
    @test em.β == vcat(CI.(-3:3), CI.(-2:2))

    em = light_cone_embedding(s, 1, 2, 2, 1, bc)
    @test em.β == vcat(CI.(-4:4), CI.(-2:2))

    em = light_cone_embedding(s, 1, 1, 2, 2, bc)
    @test em.β == vcat(CI.(-4:4), CI.(-2:2))

    em = light_cone_embedding(s, 2, 2, 2, 1, bc)
    @test em.β == vcat(CI.(-6:6), CI.(-4:4), CI.(-2:2))

    @test light_cone_embedding(s, 2, 2, 2, 1, bc) == STE(s,(γ=2, τ=2, r=2, c=1, bc=bc))
end

@testset "light_cone_embedding 2D" begin
    Φ = 2
    bc=ConstantBoundary(10.)
    s = [rand(Float64,([10 for i=1:Φ]...,)) for i=1:10]

    CI = CartesianIndex
    em = light_cone_embedding(s, 0, 1, 0, 0, bc)
    @test em.β == [CI(0,0)]

    em = light_cone_embedding(s, 0, 1, 0, 1, bc)
    @test em.β == [CI(0,0)]

    em = light_cone_embedding(s, 0, 20, 1, 10, bc)
    @test CI(0,1) ∈ em.β
    @test CI(0,0) ∈ em.β
    @test CI(0,-1) ∈ em.β
    @test CI(1,0) ∈ em.β
    @test CI(-1,0) ∈ em.β

    em = light_cone_embedding(s, 0, 5, 1.5, 17, bc)
    @test CI(0,1) ∈ em.β
    @test CI(0,0) ∈ em.β
    @test CI(0,-1) ∈ em.β
    @test CI(1,0) ∈ em.β
    @test CI(-1,0) ∈ em.β
    @test CI(-1,-1) ∈ em.β
    @test CI(-1,1) ∈ em.β
    @test CI(1,1) ∈ em.β
    @test CI(1,-1) ∈ em.β
end

include("system_defs.jl")

println("Testing PCA Functions")
@testset "PCAEmbedding" begin
    U,V = barkley(400;tskip=200)
    let γ=5, τ=1, B=1,k=1
        BC = PeriodicBoundary()
        em = cubic_shell_embedding(U,γ,τ,B,k,BC)
        RR = reconstruct(U, em)
        meanv = mean(mean.(U))
        R = reshape(reinterpret(Float64, RR.data), (size(RR)[2], size(RR)[1]))
        R .-= meanv
        covmat = Statistics.cov(R,dims =2)
        pcaem = PCAEmbedding(U,em)
        #Compare cov matrices
        @test norm(pcaem.covmat-covmat) < 1e-4
        #Compare Outdim
        drmodel = MS.fit(MultivariateStats.PCA,R)
        @test MS.outdim(drmodel) == MS.outdim(pcaem.drmodel)
        #Compare singular values
        @test MS.principalvars(drmodel) ≈ MS.principalvars(pcaem.drmodel) atol=0.5

    end
end


text:
using Test, YaoBlocks, YaoBase
using LuxurySparse

@testset "baisc" begin
    @test repeat(5, X, 2, 3) isa RepeatedBlock
    @test repeat(5, X, [2, 3]) isa RepeatedBlock
    @test repeat(X, (2,3))(5) == repeat(5, X, (2,3))
    @test repeat(X)(5) == repeat(5, X, 1:5)
    reg = rand_state(5)
    @test apply(reg, repeat(5, X, ())) == reg
    rp = RepeatedBlock(5, X, (1, 2, 3))
    @test iscommute(RepeatedBlock(5, X, (1, 2, 3)), repeat(5, Val(im)*X, (1, 2, 3)))
    @test_throws QubitMismatchError iscommute(RepeatedBlock(5, X, (1, 2, 3)), repeat(3, Val(im)*X, (1, 2, 3)))
    @test iscommute(RepeatedBlock(5, X, (1, 2, 3)), repeat(5, Val(im)*X, (1, 2)))
    @test isreflexive(rp)
    @test ishermitian(rp)
    @test isunitary(rp)
    @test (chsubblocks(rp, [Z]) |> subblocks .== [Z]) |> all
    @test occupied_locs(rp) == (1, 2, 3)
    @test rp |> copy == rp
    @test YaoBlocks.PropertyTrait(rp) == YaoBlocks.PreserveAll()

    @test repeat(10, H, 1:10) == repeat(10, H, Tuple(1:10))
    @test_throws ArgumentError repeat(10, SWAP, (1, 3))

    @test_throws QubitMismatchError apply!(rand_state(2), repeat(10, X, (1, 2, 3)))
    @test_throws QubitMismatchError apply!(rand_state(2), repeat(10, put(1, 1 => X), ()))

    @test mat(repeat(10, X, ())) == IMatrix{1<<10}()
end


text:
	@enum  CartesianDirection XDirection = 1 YDirection ZDirection
	@enum  CartesianAxis XAxis = 1 YAxis ZAxis
	@enum  FieldComponent XComponent = 1 YComponent ZComponent


text:
using PyGen
using Base.Test

# write your own tests here
@pygen function λ()
    yield(10)
end

# Test for \n in first line issue
@test isdefined(:λ) == true

# Test that the poodle function doesn't error out
@pygen function 🐩()
    yield("🐩")
end

@test isdefined(:🐩) == true

# Test that one of these generators actually works
@pygen function squares(n)
    i = 0
    while i <= n
        yield(i^2)
        i += 1
    end
end

@test sum(squares(10)) == sum(i^2 for i in 1:10)


text:
using Test
using PolynomialBasis
using ImplicitDomainQuadrature
# using Revise
using CutCell
include("useful_routines.jl")

function bulk_modulus(l, m)
    return l + 2m / 3
end

function lame_lambda(k, m)
    return k - 2m / 3
end

function analytical_coefficient_matrix(inradius, outradius, ls, ms, lc, mc)
    a = zeros(3, 3)
    a[1, 1] = inradius
    a[1, 2] = -inradius
    a[1, 3] = -1.0 / inradius
    a[2, 1] = 2 * (lc + mc)
    a[2, 2] = -2 * (ls + ms)
    a[2, 3] = 2ms / inradius^2
    a[3, 2] = 2(ls + ms)
    a[3, 3] = -2ms / outradius^2
    return a
end

function analytical_coefficient_rhs(ls, ms, theta0)
    r = zeros(3)
    Ks = bulk_modulus(ls, ms)
    r[2] = -Ks * theta0
    r[3] = Ks * theta0
    return r
end

struct AnalyticalSolution
    inradius::Any
    outradius::Any
    center::Any
    A1c::Any
    A1s::Any
    A2s::Any
    ls::Any
    ms::Any
    lc::Any
    mc::Any
    theta0::Any
    function AnalyticalSolution(inradius, outradius, center, ls, ms, lc, mc, theta0)
        a = analytical_coefficient_matrix(inradius, outradius, ls, ms, lc, mc)
        r = analytical_coefficient_rhs(ls, ms, theta0)
        coeffs = a \ r
        new(
            inradius,
            outradius,
            center,
            coeffs[1],
            coeffs[2],
            coeffs[3],
            ls,
            ms,
            lc,
            mc,
            theta0,
        )
    end
end

function radial_displacement(A::AnalyticalSolution, r)
    if r <= A.inradius
        return A.A1c * r
    else
        return A.A1s * r + A.A2s / r
    end
end

function (A::AnalyticalSolution)(x)
    relpos = x - A.center
    r = sqrt(relpos' * relpos)
    ur = radial_displacement(A, r)
    if ur ≈ 0.0
        [0.0, 0.0]
    else
        costheta = (x[1] - A.center[1]) / r
        sintheta = (x[2] - A.center[2]) / r
        u1 = ur * costheta
        u2 = ur * sintheta
        return [u1, u2]
    end
end

function shell_radial_stress(ls, ms, theta0, A1, A2, r)
    return (ls + 2ms) * (A1 - A2 / r^2) + ls * (A1 + A2 / r^2) - (ls + 2ms / 3) * theta0
end

function shell_circumferential_stress(ls, ms, theta0, A1, A2, r)
    return ls * (A1 - A2 / r^2) + (ls + 2ms) * (A1 + A2 / r^2) - (ls + 2ms / 3) * theta0
end

function shell_out_of_plane_stress(ls, ms, A1, theta0)
    return 2 * ls * A1 - (ls + 2ms / 3) * theta0
end

function core_in_plane_stress(lc, mc, A1)
    return (lc + 2mc) * A1 + lc * A1
end

function core_out_of_plane_stress(lc, A1)
    return 2 * lc * A1
end

function rotation_matrix(x, r)
    costheta = x[1] / r
    sintheta = x[2] / r
    Q = [
        costheta -sintheta
        sintheta costheta
    ]
    return Q
end

function shell_stress(A::AnalyticalSolution, x)
    relpos = x - A.center
    r = sqrt(relpos' * relpos)
    Q = rotation_matrix(relpos, r)

    srr = shell_radial_stress(A.ls, A.ms, A.theta0, A.A1s, A.A2s, r)
    stt = shell_circumferential_stress(A.ls, A.ms, A.theta0, A.A1s, A.A2s, r)

    cylstress = [
        srr 0.0
        0.0 stt
    ]

    cartstress = Q * cylstress * Q'
    s11 = cartstress[1, 1]
    s22 = cartstress[2, 2]
    s12 = cartstress[1, 2]
    s33 = shell_out_of_plane_stress(A.ls, A.ms, A.A1s, A.theta0)

    return [s11, s22, s12, s33]
end

function core_stress(A::AnalyticalSolution)
    s11 = core_in_plane_stress(A.lc, A.mc, A.A1c)
    s33 = core_out_of_plane_stress(A.lc, A.A1c)
    return [s11, s11, 0.0, s33]
end

function exact_stress(A::AnalyticalSolution, x)
    relpos = x - A.center
    r = sqrt(relpos' * relpos)
    if r < A.inradius
        return core_stress(A)
    else
        return shell_stress(A, x)
    end
end

function onleftboundary(x, L, W)
    return x[1] ≈ 0.0
end

function onbottomboundary(x, L, W)
    return x[2] ≈ 0.0
end

function onrightboundary(x, L, W)
    return x[1] ≈ L
end

function ontopboundary(x, L, W)
    return x[2] ≈ W
end

function onboundary(x, L, W)
    return onbottomboundary(x, L, W) ||
           onrightboundary(x, L, W) ||
           ontopboundary(x, L, W) ||
           onleftboundary(x, L, W)
end

function solve_and_compute_displacement_error(
    width,
    center,
    inradius,
    outradius,
    stiffness,
    theta0,
    nelmts,
    polyorder,
    numqp,
    penaltyfactor,
)

    lambda1, mu1 = CutCell.lame_coefficients(stiffness, +1)
    lambda2, mu2 = CutCell.lame_coefficients(stiffness, -1)
    transfstress = CutCell.plane_strain_transformation_stress(lambda1, mu1, theta0)

    analyticalsolution =
        AnalyticalSolution(inradius, outradius, center, lambda1, mu1, lambda2, mu2, theta0)

    dx = width / nelmts
    meanmoduli = 0.5 * (lambda1 + lambda2 + mu1 + mu2)
    penalty = penaltyfactor / dx * meanmoduli

    basis = TensorProductBasis(2, polyorder)
    mesh = CutCell.Mesh([0.0, 0.0], [width, width], [nelmts, nelmts], basis)
    levelset = InterpolatingPolynomial(1, basis)
    levelsetcoeffs = CutCell.levelset_coefficients(
        x -> -circle_distance_function(x, center, inradius),
        mesh,
    )

    cutmesh = CutCell.CutMesh(levelset, levelsetcoeffs, mesh)
    cellquads = CutCell.CellQuadratures(levelset, levelsetcoeffs, cutmesh, numqp)
    interfacequads = CutCell.InterfaceQuadratures(levelset, levelsetcoeffs, cutmesh, numqp)
    facequads = CutCell.FaceQuadratures(levelset, levelsetcoeffs, cutmesh, numqp)

    bilinearforms = CutCell.BilinearForms(basis, cellquads, stiffness, cutmesh)
    interfacecondition =
        CutCell.incoherent_interface_condition(basis, interfacequads, stiffness, cutmesh, penalty)

    leftdisplacementbc = CutCell.DisplacementComponentCondition(
        x->analyticalsolution(x)[1],
        basis,
        facequads,
        stiffness,
        cutmesh,
        x -> onleftboundary(x, width, width),
        [1.0, 0.0],
        penalty,
    )
    bottomdisplacementbc = CutCell.DisplacementCondition(
        analyticalsolution,
        basis,
        facequads,
        stiffness,
        cutmesh,
        x -> onbottomboundary(x, width, width),
        penalty,
    )
    topdisplacementbc = CutCell.DisplacementCondition(
        analyticalsolution,
        basis,
        facequads,
        stiffness,
        cutmesh,
        x -> ontopboundary(x, width, width),
        penalty,
    )

    sysmatrix = CutCell.SystemMatrix()
    sysrhs = CutCell.SystemRHS()

    CutCell.assemble_bilinear_form!(sysmatrix, bilinearforms, cutmesh)
    CutCell.assemble_interface_condition!(sysmatrix, interfacecondition, cutmesh)
    CutCell.assemble_bulk_transformation_linear_form!(
        sysrhs,
        transfstress,
        basis,
        cellquads,
        cutmesh,
    )
    CutCell.assemble_incoherent_interface_transformation_rhs!(
        sysrhs,
        transfstress,
        basis,
        interfacequads,
        cutmesh,
    )
    CutCell.assemble_traction_force_linear_form!(
        sysrhs,
        x -> exact_stress(analyticalsolution, x)[[1, 3]],
        basis,
        facequads,
        cutmesh,
        x -> onrightboundary(x, width, width),
    )
    CutCell.assemble_traction_force_component_linear_form!(
        sysrhs,
        x -> -exact_stress(analyticalsolution, x)[3],
        basis,
        facequads,
        cutmesh,
        x -> onleftboundary(x, width, width),
        [0.0, 1.0],
    )
    CutCell.assemble_penalty_displacement_bc!(
        sysmatrix,
        sysrhs,
        leftdisplacementbc,
        cutmesh,
    )
    CutCell.assemble_penalty_displacement_bc!(
        sysmatrix,
        sysrhs,
        bottomdisplacementbc,
        cutmesh,
    )
    CutCell.assemble_penalty_displacement_bc!(sysmatrix, sysrhs, topdisplacementbc, cutmesh)
    CutCell.assemble_penalty_displacement_component_transformation_rhs!(
        sysrhs,
        transfstress,
        basis,
        facequads,
        cutmesh,
        x -> onleftboundary(x, width, width),
        [1.0, 0.0],
    )
    CutCell.assemble_penalty_displacement_transformation_rhs!(
        sysrhs,
        transfstress,
        basis,
        facequads,
        cutmesh,
        x -> onbottomboundary(x, width, width),
    )
    CutCell.assemble_penalty_displacement_transformation_rhs!(
        sysrhs,
        transfstress,
        basis,
        facequads,
        cutmesh,
        x -> ontopboundary(x, width, width),
    )

    matrix = CutCell.make_sparse(sysmatrix, cutmesh)
    rhs = CutCell.rhs(sysrhs, cutmesh)

    sol = matrix \ rhs
    disp = reshape(sol, 2, :)

    err = mesh_L2_error(disp, analyticalsolution, basis, cellquads, cutmesh)
    den = integral_norm_on_cut_mesh(analyticalsolution, cellquads, cutmesh, 2)

    return err ./ den
end

lambda1, mu1 = 100.0, 80.0
lambda2, mu2 = 80.0, 60.0
theta0 = -0.067
stiffness = CutCell.HookeStiffness(lambda1, mu1, lambda2, mu2)

width = 1.0
penaltyfactor = 1e2

polyorder = 2
numqp = required_quadrature_order(polyorder) + 2

center = [width / 2, width / 2]
inradius = width / 4
outradius = width

powers = [3,4,5]
nelmts = [2^p + 1 for p in powers]

err = [
    solve_and_compute_displacement_error(
        width,
        center,
        inradius,
        outradius,
        stiffness,
        theta0,
        ne,
        polyorder,
        numqp,
        penaltyfactor,
    ) for ne in nelmts
]

dx = 1.0 ./ nelmts
u1err = [er[1] for er in err]
u2err = [er[2] for er in err]

u1rate = convergence_rate(dx,u1err)
u2rate = convergence_rate(dx,u2err)

@test allapprox(u1rate,repeat([3.0],length(u1rate)),0.15)
@test allapprox(u2rate,repeat([3.0],length(u2rate)),0.15)


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

using JuliaFEM
using JuliaFEM.Testing

function get_testproblems(u, la)
    nodes = Dict{Int64, Node}(
        1 => [0.0, 0.0],
        2 => [2.0, 0.0],
        3 => [1.0, 2.0],
        4 => [0.0, 2.0],
        5 => [2.5, 0.0],
        6 => [4.5, 0.0],
        7 => [4.5, 1.0],
        8 => [2.5, 1.0])
    displacement = Dict{Int64, Vector{Float64}}()
    reaction_force = Dict{Int64, Vector{Float64}}()
    for i=1:8
        displacement[i] = u[:,i]
        reaction_force[i] = la[:,i]
    end
    bc5 = Seg2([3, 2])
    bc6 = Seg2([8, 5])
    update!([bc5, bc6], "geometry", nodes)
    update!([bc5, bc6], "displacement", displacement)
    update!([bc5, bc6], "reaction force", reaction_force)
    bc5["master elements"] = [bc6]
    contact1 = Problem(Mortar, "contact between bodies", 2, "displacement")
    contact2 = Problem(Mortar, "contact between bodies", 2, "displacement")
    contact2.properties.formulation = :forwarddiff
    contact2.assembly.u = vec(u)
    contact2.assembly.la = vec(la)
    push!(contact1, bc5, bc6)
    push!(contact2, bc5, bc6)
    return contact1, contact2
end

#= TODO: Fix test
@testset "test linearization of contact force in undeformed state" begin
    u = zeros(2, 8)
    la = zeros(2, 8)
    contact1, contact2 = get_testproblems(u, la)
    assemble!(contact1, 0.0)
    assemble!(contact2, 0.0)
    @test isapprox(full(contact1.assembly.C1), full(contact2.assembly.C1))
    @test isapprox(full(contact1.assembly.K), full(contact2.assembly.K))
end
=#


text:
# deprecated functions

using Base: @deprecate

@deprecate add!(a::Matrix, b::AbstractPDMat) pdadd!(a, b)

@deprecate add_scal!(a::Matrix, b::AbstractPDMat, c::Real) pdadd!(a, b, c)

@deprecate add_scal(a::Matrix, b::AbstractPDMat, c::Real) pdadd(a, b, c)

@deprecate full(x::AbstractPDMat) Matrix(x)


text:
#__precompile__()
module GraphVisualize
# Credits for original code: Simon Danisch
# Contributors: Carlo Lucibello

export plot

using GLVisualize
using LightGraphs
using GraphLayout
using GeometryTypes, GLAbstraction, Reactive, GLWindow, GLFW
import Base: push!

type GraphPlot
    glast::Graph
    g::Graph
    resolution::Tuple{Int,Int}
    point_size::Float32

    #SIGNALS
    s_g::Signal{Graph}
    s_vertex_pos::Signal
end

function GraphPlot(g::Graph; resolution = (900, 900), point_size = 15f0)
    s_g = Signal(g)
    s_vertex_pos = Signal(layout(g, resolution, point_size))
    plt = GraphPlot(Graph(), g, resolution, point_size, s_g, s_vertex_pos)
    return plt
end

function set_observation_fps!(plt::GraphPlot, f)
    preserve(map(delta -> begin
            if has_changed(plt.g, plt.glast)
                update!(plt)
            end
        end, f))
end

function update!(plt)
    g = plt.g; glast = plt.glast
    res = plt.resolution; psize = plt.point_size

    @assert has_changed(g, glast)
    push!(plt.s_g, g)

    if vertex_added(g, glast)
        push!(plt.s_vertex_pos, push!(value(plt.s_vertex_pos), new_pos(res, psize)))
    elseif vertex_removed(g, glast)
        # Assume that the removed vertex was exchanged with last vertex
        # This is an hack, and it is slow and buggy.
        # In order to have a graphically consistent vertex removal,
        # we need either unique vertex identifiers or some external signaling.

        lastneigs = neighbors(glast, nv(glast))
        ilast = 0
        i = 1
        for neigs in fadj(g)
            if neigs == filter(x->x!=i,lastneigs)
                neigs != neighbors(glast,i) && break
                ilast = i
            end
            i += 1
        end
        i > nv(g) && (i = ilast)
        pos = value(plt.s_vertex_pos)
        pos[i] = pos[end]
        resize!(pos, length(pos)-1)
        push!(plt.s_vertex_pos, pos)
    elseif nv(g) == nv(glast)
        # nothing
    else #nv(g) != nv(glast)
        push!(plt.s_vertex_pos, layout(g, res, psize))
    end

    plt.glast = deepcopy(g)
    sleep(1e-3) # BUG: removing this results in error messages
end

push!(plt::GraphPlot, g::Graph) = push!(plt.s_g, g)

edge_removed(g::Graph, glast::Graph) = ne(glast) == ne(g) + 1
edge_added(g::Graph, glast::Graph) = ne(glast) == ne(g) - 1
vertex_removed(g::Graph, glast::Graph) = nv(glast) == nv(g) + 1
vertex_added(g::Graph, glast::Graph) = nv(glast) == nv(g) - 1

has_changed(g::Graph, glast::Graph) = !(nv(glast) == nv(g) &&  ne(glast) == ne(g))
has_really_changed(g::Graph, glast::Graph) = !(glast == g)

function edgelist(g::Graph)
    indices = Vector{Int32}(2ne(g))
    i = 0
    for (u,v) in edges(g)
        indices[i+=1] = u
        indices[i+=1] = v
    end
    indices
end

function new_pos(resolution, psize)
    wsize = minimum(resolution)
    x, y = (2rand()-1, 2rand()-1)
    return Point2f0(2psize,2psize) + (Point2f0(x, y) + Point2f0(1.,1.))* (wsize-4psize)/2
end

function layout(g::Graph, resolution, psize)
    wsize = minimum(resolution)
    x, y = layout_spring_adj(full(adjacency_matrix(g)))
    return [Point2f0(2psize,2psize) + (Point2f0(x[i], y[i]) + Point2f0(1.,1.))* (wsize-4psize)/2 for i=1:nv(g)]
end

"""
    plot(g::Graph; observe=true, resolution = (900, 900), point_size = 15f0)

Creates an OpenGL window and draws `g` in it. If `observe==true` future modifications of
`g` will be reflected in the plot. Returns a `GraphPlot` object.
"""
function plot(g::Graph; observe= true, kw...)
    plt = GraphPlot(g, kw...)
    observe && set_observation_fps!(plt, fps(10.))  # WEIRDNESS: if I inserte this line in the GraphPlot constructor
                                                    #           it doesn't update
    plot(plt)
    return plt
end

function plot(plt::GraphPlot)

    window = glscreen(resolution=plt.resolution)

    const record_interactive = true
    points = visualize((Circle(Point2f0(0), plt.point_size), plt.s_vertex_pos))

    const point_robj = points.children[]
    const gpu_position = point_robj[:position]
    elist = map(edgelist, plt.s_g)
    lines  = visualize(gpu_position, :linesegment, indices=elist)

    # current tuple of renderobject id and index into the gpu array
    const m2id = GLWindow.mouse2id(window)
    isoverpoint = const_lift(is_same_id, m2id, point_robj)

    # inputs are a dict, materialize gets the keys out of it (equivalent to mouseposition = window.inputs[:mouseposition])
    @materialize mouse_buttons_pressed, mouseposition = window.inputs

    # single left mousekey pressed (while no other mouse key is pressed)
    key_pressed = const_lift(GLAbstraction.singlepressed, mouse_buttons_pressed, GLFW.MOUSE_BUTTON_LEFT)
    # dragg while key_pressed. Drag only starts if isoverpoint is true
    mousedragg  = GLAbstraction.dragged(mouseposition, key_pressed, isoverpoint)

    # use mousedrag and mouseid + index to actually change the gpu array with the positions
    preserve(foldp((value(m2id)..., Point2f0(0)), mousedragg) do v0, dragg
        if dragg == Vec2f0(0) # if drag just started. Not the best way, maybe dragged should return a tuple of (draggvalue, started)
            id, index = value(m2id)
            if id==point_robj.id && length(gpu_position) >= index
                p0 = gpu_position[index]
            else
                p0 = v0[3]
            end
        else
            id, index, p0 = v0
            if id==point_robj.id && length(gpu_position) >= index
                gpu_position[index] = Point2f0(p0) + Point2f0(dragg)
            end

        end
        return id, index, p0
    end)

    view(lines, window, camera=:fixed_pixel)
    view(points, window, camera=:fixed_pixel)
    @spawn renderloop(window)
end

end # module

max_stars_repo_name: JuliaGraphs/GraphVisualization

max_issues_repo_name: JuliaGraphs/GraphVisualization

max_forks_repo_name: JuliaGraphs/GraphVisualization


text:
"""
    AbstractClustering

Supertype for clusterings.
"""
abstract type AbstractClustering end

"""
    PartitionalClustering{Tc<:Integer,Tw<:Real,Ty<:Real,Tm<:Real} <: AbstractClustering

Partitional clustering.
"""
struct PartitionalClustering{Tc<:Integer,Tw<:Real,Ty<:Real} <: AbstractClustering
    r::Vector{Int}
    c::Vector{Int}
    C::Matrix{Tc}
    W::Matrix{Tw}
    Y::Matrix{Ty}
    p::NamedTuple

    function PartitionalClustering{Tc,Tw,Ty}(r::Vector{Int}, c::Vector{Int},
                                            C::Matrix{Tc}, W::Matrix{Tw},
                                            Y::Matrix{Ty}, p::NamedTuple) where {Tc<:Integer,Tw<:Real,Ty<:Real}
        size(C, 2) == size(W, 2) || throw(DimensionMismatch("dimensions of constraints and weights must match"))
        return new(r, c, C, W, Y, p)
    end
end
PartitionalClustering(r::Vector{Int}, c::Vector{Int}, C::Matrix{Tc},
                    W::Matrix{Tw}, Y::Matrix{Ty}, p::NamedTuple) where {Tc,Tw,Ty} =
    PartitionalClustering{Tc,Tw,Ty}(r, c, C, W, Y, p)

function PartitionalClustering(clust::KmeansResult)
    r = Int[]
    c = Int[]
    C = Matrix{Int}(undef, 0, 0)
    W = Matrix{Float64}(undef, 0, 0)
    k = size(clust.centers, 2)
    n = length(clust.assignments)
    Y = zeros(Int, k, n)
    for i = 1:n
        Y[clust.assignments[i], i] = 1
    end
    p = (centers=clust.centers, costs=clust.costs, counts=clust.counts,
            wcounts=clust.wcounts, totalcost=clust.totalcost,
            iterations=clust.iterations, converged=clust.converged)
    return PartitionalClustering(r, c, C, W, Y, p)
end
function PartitionalClustering(clust::KmedoidsResult)
    r = Int[]
    c = Int[]
    C = Matrix{Int}(undef, 0, 0)
    W = Matrix{Float64}(undef, 0, 0)
    k = length(clust.medoids)
    n = length(clust.assignments)
    Y = zeros(Int, k, n)
    for i = 1:n
        Y[clust.assignments[i], i] = 1
    end
    p = (medoids=clust.medoids, costs=clust.costs, counts=clust.counts,
            totalcost=clust.totalcost, iterations=clust.iterations,
            converged=clust.converged)
    return PartitionalClustering(r, c, C, W, Y, p)
end
function PartitionalClustering(clust::FuzzyCMeansResult)
    r = Int[]
    c = Int[]
    C = Matrix{Int}(undef, 0, 0)
    W = Matrix{Float64}(undef, 0, 0)
    Y = permutedims(clust.weights)
    p = (centers=clust.centers, iterations=clust.iterations,
            converged=clust.converged)
    return PartitionalClustering(r, c, C, W, Y, p)
end

# Partitional clustering equality operator
Base.:(==)(a::PartitionalClustering, b::PartitionalClustering) =
    (a.r == b.r && a.c == b.c && a.C == b.C && a.W == b.W && a.Y == b.Y
            && a.p == b.p)

# Compute hash code
Base.hash(a::PartitionalClustering, h::UInt) =
    hash(a.r, hash(a.c, hash(a.C, hash(a.W, hash(a.Y, hash(a.p,
        hash(:PartitionalClustering, h)))))))

"""
    HierarchicalClustering{Tc<:Integer,Tw<:Real} <: AbstractClustering

Hierarchical clustering.
"""
struct HierarchicalClustering{Tc<:Integer,Tw<:Real} <: AbstractClustering
    r::Vector{Int}
    c::Vector{Int}
    C::Array{Tc,3}
    W::Array{Tw,3}
    p::NamedTuple

    function HierarchicalClustering{Tc,Tw}(r::Vector{Int}, c::Vector{Int},
                                        C::Array{Tc,3}, W::Array{Tw,3},
                                        p::NamedTuple) where {Tc<:Integer,Tw<:Real}
        size(C, 2) == size(W, 2) || throw(DimensionMismatch("dimensions of constraints and weights must match"))
        return new(r, c, C, W, p)
    end
end
HierarchicalClustering(r::Vector{Int}, c::Vector{Int}, C::Array{Tc,3},
                    W::Array{Tw,3}, p::NamedTuple) where {Tc,Tw} =
    HierarchicalClustering{Tc,Tw}(r, c, C, W, p)

"""
    axes(a::AbstractClustering[, d])

Access the feature and instance identifiers. Optionally, specify dimension `d`
to get the identifiers of that dimension only.
"""
Base.axes(a::AbstractClustering) = a.r, a.c
Base.axes(a::AbstractClustering, d) = d::Integer <= 2 ? axes(a)[d] : Int[]

"""
    features(a::AbstractClustering)

Access the feature identifiers.
"""
features(a::AbstractClustering) = axes(a, 1)

"""
    instances(a::AbstractClustering)

Access the instance identifiers.
"""
Base.instances(a::AbstractClustering) = axes(a, 2)

"""
    constraints(a::AbstractClustering)

Access the contraints.
"""
constraints(a::AbstractClustering) = a.C

"""
    weights(a::AbstractClustering)

Access the weights.
"""
weights(a::AbstractClustering) = a.W

"""
    parameters(c::PartitionalClustering)

Access the parameters.
"""
parameters(a::AbstractClustering) = a.p

"""
    assignments(a::PartitionalClustering)

Access the assignments.
"""
assignments(a::PartitionalClustering) = a.Y


text:
"""
    ErrorModel(mismatch, insertion, deletion, codon_insertion, codon_deletion)

Error model for sequencing.

Each field contains the relative rate of of that kind of error. For
instance, this model breaks the error rate into 80% mismatches, 10%
codon insertions, and 10% codon deletions:
`ErrorModel(8, 0, 0, 1, 1)`.

# Fields:
- `mismatch::Real`
- `insertion::Real`
- `deletion::Real`
- `codon_insertion::Real`
- `codon_deletion::Real`

"""
struct ErrorModel
    mismatch::Real
    insertion::Real
    deletion::Real
    codon_insertion::Real
    codon_deletion::Real
end

"""error model without codon errors"""
function ErrorModel(mismatch, insertion, deletion)
    return ErrorModel(mismatch, insertion, deletion, 0, 0)
end

"""turn error rates into probabilities"""
function normalize(errors::ErrorModel)
    args = [errors.mismatch,
            errors.insertion,
            errors.deletion,
            errors.codon_insertion,
            errors.codon_deletion]
    m, i, d, ci, cd = normalize(args)
    return ErrorModel(m, i, d, ci, cd)
end

struct Scores
    mismatch::Score
    insertion::Score
    deletion::Score
    codon_insertion::Score
    codon_deletion::Score
end

"""
    Scores(errors; mismatch, insertion, deletion)

Derive alignment scores from an error model.

Takes extra penalties to add to the mismatch, insertion, and deletion
scores.

# Arguments:
- `errors::ErrorModel`:
- `mismatch::Real`: substitution
- `insertion::Real`: insertion
- `deletion::Real`: deletion

"""
function Scores(errors::ErrorModel;
                mismatch=0.0,
                insertion=0.0,
                deletion=0.0)
    args = [errors.mismatch,
            errors.insertion,
            errors.deletion,
            errors.codon_insertion,
            errors.codon_deletion]
    m, i, d, ci, cd = log10.(normalize(args))
    return Scores(Score(m + mismatch),
                  Score(i + insertion),
                  Score(d + deletion),
                  Score(ci + 3 * insertion),
                  Score(cd + 3 * deletion))
end


text:
"""
    showpolicy([io], [mime], m::MDP, p::Policy)
    showpolicy([io], [mime], statelist::AbstractVector, p::Policy)
    showpolicy(...; pre=" ")

Print the states in `m` or `statelist` and the actions from policy `p` corresponding to those states.

For the MDP version, if `io[:limit]` is `true`, will only print enough states to fill the display.
"""
function showpolicy(io::IO, mime::MIME"text/plain", m::MDP, p::Policy; pre=" ", kwargs...)
    slist = nothing
    truncated = false
    limited = get(io, :limit, false)
    rows = first(get(io, :displaysize, displaysize(io)))
    rows -= 3 # Yuck! This magic number is also in Base.print_matrix
    try
        if limited && length(states(m)) > rows
            slist = collect(take(states(m), rows-1))
            truncated = true
        else
            slist = collect(states(m))
        end
    catch ex
        @info("""Unable to pretty-print policy:
              $(sprint(showerror, ex))
              """)
        show(io, mime, m)
        return show(io, mime, p)
    end
    showpolicy(io, mime, slist, p; pre=pre, kwargs...)
    if truncated
        print(io, '\n', pre, "…")
    end
end

function showpolicy(io::IO, mime::MIME"text/plain", slist::AbstractVector, p::Policy; pre::AbstractString=" ")
    S = eltype(slist)
    sa_con = IOContext(io, :compact => true)

    if !isempty(slist)
        # print first element without a newline
        print(io, pre)
        print_sa(sa_con, first(slist), p, S)

        # print all other elements
        for s in slist[2:end]
            print(io, '\n', pre)
            print_sa(sa_con, s, p, S)
        end
    end
end

showpolicy(io::IO, m::Union{MDP,AbstractVector}, p::Policy; kwargs...) = showpolicy(io, MIME("text/plain"), m, p; kwargs...)
showpolicy(m::Union{MDP,AbstractVector}, p::Policy; kwargs...) = showpolicy(stdout, m, p; kwargs...)

function print_sa(io::IO, s, p::Policy, S::Type)
    show(IOContext(io, :typeinfo => S), s)
    print(io, " -> ")
    try
        show(io, action(p, s))
    catch ex
        showerror(IOContext(io, :limit=>true), ex)
    end
end

max_stars_repo_name: JuliaPOMDP/POMDPPolicies

max_issues_repo_name: JuliaPOMDP/POMDPPolicies

max_forks_repo_name: JuliaPOMDP/POMDPPolicies


text:
using LinearAlgebra
#using ArgParse
using Base.Threads
using SIMD
using TimerOutputs
#using Arpack
#using SparseArrays

"""
    main(tJ,S,L,nums,PBC)

-tJ:      coefficient in the Hamiltonian
-S:       spin
-L:       length of chain
-nums:    # of states of interest
PBC:      oeriodic boundary condition
"""
function main(tJ,S,L,nums,PBC;
              Slist=[0],is_show=false,retbasis=false)
    to = TimerOutput()
    ret = [0,0]
    ### prep. for calc.
    @timeit to "prep." begin
        basisvec = makeV(S,L)
        print("L=$L Full Dim.=",length(basisvec),"\t")
        vec = basisvec[1]
        vecsSz = [ [copy(vec)]]; deleteat!(vecsSz[1],1)
        bits = Int64[]
        for tvec in basisvec
            bitarr_to_int!(tvec,ret)
            bit = ret[1]
            lup = count_ones(bit)
            ldown = length(tvec) - lup
            tSz = Int64(0.5*(lup-ldown))
            for (i,Sz) in enumerate(Slist)
                if tSz != Sz; continue;end
                push!(vecsSz[1],tvec)
                push!(bits,bit)
            end
        end
        basisvec = nothing
    end

    i = 1
    tvec = vecsSz[1]; tDim = length(tvec)
    println("Dim.(Sz=0) $tDim")
    solve_1dHeisenberg(tJ,S,L,nums,PBC,tDim,tvec,to)
    if is_show;show(to, allocations = true,compact =false);println("");end

    if retbasis; return vecsSz[1];end
    return nothing
end

function bitarr_to_int!(arr::Array{Bool,1},ret)
    ret[1] = 0
    ret[2] = 2^(length(arr)-1)
    for i in eachindex(arr)
        ret[1] += ret[2]*arr[i]
        ret[2] >>= 1
    end
    return nothing
end

function solve_1dHeisenberg(tJ,S,L,nums,PBC,tDim,tvec,to)
    ## potential term
    @timeit to "potential term" begin
        Ufacs=zeros(Float64,tDim)
        @inbounds for ith=1:tDim
            @simd for j = 1:L-1
                v1 = ifelse(tvec[ith][j]  ,1.0,-1.0)
                v2 = ifelse(tvec[ith][j+1],1.0,-1.0)
                Ufacs[ith] += 0.25 * (v1 .*v2)
            end
            if PBC
                v1 = ifelse(tvec[ith][L],1.0,-1.0)
                v2 = ifelse(tvec[ith][1],1.0,-1.0)
                Ufacs[ith] += 0.25* (v1 .* v2)
            end
        end
    end
    ## hopping term
    @timeit to "hopping term" begin
        non0s=[ Int64[] for ith=1:tDim]
        facs=[Float64[] for ith=1:tDim]
        @inbounds @threads for ith=1:tDim
            TF = [false];ret=[0,0]
            tv1 = tvec[ith]
            t_non0s = non0s[ith]
            t_facs  = facs[ith]
            @inbounds for jth=ith+1:tDim
                tv2=tvec[jth]
                connectable(TF,tv1,tv2)
                if TF[1]==false;continue;end
                neighbor(TF,ret,tv1,tv2,L,PBC)
                if TF[1]
                    v1  = 0.5 .* ifelse(tv2[ret[1]],1.0,-1.0)
                    v2  = 0.5 .* ifelse(tv2[ret[2]],1.0,-1.0)
                    fac = calcfac(S,v1,v2)
                    push!(t_non0s,jth)
                    push!(t_facs,fac)
                end
            end
        end
    end
    @timeit to "Lanczos" Lanczos(tJ,L,non0s,facs,Ufacs,tDim,nums)
    return nothing
end

"""
    makeV(S,L)

preparing possible spin configurations
"""
function makeV(S,L)
    #sarr = collect(-Int64(2*S):2:Int64(2*S))
    sarr = [false,true]
    t = prodS(sarr,L)
end

"""
    prodS(sarr,N)

generating product of spin states to make configs.
"""
function prodS(sarr::Array{Bool,1},N::Int64)
    ln = length(sarr)
    tbasisvecs = copy(sarr)
    for i = 1:N-1
        nlist = []
        for tmp in tbasisvecs
            for val in sarr
                if typeof(tmp) == Bool
                    ctmp = [tmp]
                else
                    ctmp = copy(tmp)
                end
                push!(nlist,append!(ctmp,[val]))
            end
        end
        tbasisvecs = nlist
    end
    Dim = ln^N
    if length(tbasisvecs) != Dim;println("err prodS");exit();end
    return tbasisvecs
end

function calcfac(S::F,Sz_cre::F,Sz_ani::F) where{F<:Float64}
    sqrt( (S*(S+1.0) - Sz_cre*(Sz_cre+1.0)) * (S*(S+1.0) - Sz_ani*(Sz_ani-1.0)))
end

function neighbor(TF,ret,ta1,ta2,N::Int64,PBC=false) where{IA<:Array{Int64,1}}
    ij = Int64[ ]
    for i=1:length(ta1)
        dif = ta1[i]-ta2[i]
        if dif != 0
            push!(ij,i)
        end
    end
    if PBC ==false
        if abs(ij[1]-ij[2])== 1
            if ta1[ij[1]]-ta2[ij[1]] == 1
                TF[1]= true; ret[1]=ij[1]; ret[2]=ij[2]
            else
                TF[1]= true; ret[2]=ij[1]; ret[1]=ij[2]
            end
        else
            TF[1]=false ; ret[1]=0; ret[2]=0
        end
    else ### not checked
        if abs(ij[1]-ij[2])== 1 || abs(ij[1]-ij[2])==N-1
            if dif[ij[1]] == 1
                TF[1]= true; ret[1]=ij[1]; ret[2]=ij[2]
            else
                TF[1]= true; ret[2]=ij[1]; ret[1]=ij[2]
            end
        else
            TF[1]=false ; ret[1]=0; ret[2]=0
        end
    end
    nothing
end

function connectable(tf,arr1,arr2)
    s1 = 0; s2 = 0
    for (i,t1) in enumerate(arr1)
       t2 = arr2[i]
       dif = Int(t1)-Int(t2)
       s1 += abs(dif)
       s2 += dif
    end
    tf[1] = ifelse( s1 != 2 || s2 != 0, false,true)
    return nothing
end

function operateH(v::Array{Float64,1},t::Float64,
                  non0s::Array{Array{Int64,1},1},
                  facs::Array{Array{Float64,1},1},
                  Ufacs::Array{Float64,1},tDim::Int64)
    ### diagonal part
    w = Ufacs .* v
    ### non-diagonal part
    tvec = 0.0 * v
    for ith = 1:tDim-1
        coef = facs[ith]
        for (ind, jth) in enumerate(non0s[ith])
            fac = 0.5 * t * coef[ind]
            tvec[jth] += fac * v[ith]
            tvec[ith] += fac * v[jth]
        end
    end
    w += tvec
    return w
end

function Lanczos(t::Float64,N::Int64,non0s::AAI,
                 facs::AAF,Ufacs::FA,tDim::Int64,
                 numstate::Int64,lm=300) where{FA<:Array{Float64,1},
                                                     AAI<:Array{Array{Int64,1},1},
                                               AAF<:Array{Array{Float64,1},1}}
    vkm1 = zeros(Float64,tDim)
    vk   = zeros(Float64,tDim); vk[1]=1.0
    betakm1 = 0.0; betak=0.0; alphas = Float64[]; betas = Float64[]
    en = [1.e+10]; TF = [false]
    for k = 1:lm
        w = operateH(vk,t,non0s,facs,Ufacs,tDim)  - betak .* vkm1
        alphak = dot(vk,w)
        push!(alphas,alphak)
        en = convcheck(k,numstate,tDim,alphas,betas,en,TF)
        if TF[1];println("En=",en[1:minimum([numstate,tDim])], " Eg.s/N ", en[1]/N, " Lanczos it=$k"); break; end
        if k==lm
            try;print_vec("En=",en[1:minimum([tDim,numstate])]);
            catch;print_vec("En",en);end
        end
        w .-= alphak .* vk
        betakm1 = betak; betak = sqrt(dot(w,w)); push!(betas,betak)
        if betak == 0.0; println("betak=0.0");break;end
        vkm1 .= vk; vk .= w .* (1.0/betak)
    end
end

function convcheck(k::Int64,numstate::Int64,tDim::Int64,
                   alphas::FA,betas::FA,en::FA,
                   TF,
                   tol=1.e-14) where{FA<:Array{Float64,1}}
    if k < numstate
        return en
    else
        T = zeros(Float64,k,k)
        for i = 1:k
            T[i,i] =alphas[i]
        end
        for j = 2:k
            i = j-1
            T[i,j] = betas[i];T[j,i] = betas[i]
        end
        iend = minimum([tDim,numstate])
        nen = eigvals(T)[1:iend]
        if k == numstate; return nen; end
        dif = maximum( abs.(en[1:iend] .- nen[1:iend]))
        if dif < tol
            TF[1] = true
            return nen
        else
            return nen
        end
    end
end

function samplerun()
    tJ = 1.0;S = 0.5;nums = 1;PBC = false;L = 4
    for L=4:2:16
        @time main(tJ,S,L,nums,PBC)
    end
end

function dmrg()
    tJ = 1.0;S = 0.5;nums = 1;PBC = false;L = 4
    println("\n\n")
    @time v_sys = main(tJ,S,L,nums,PBC;retbasis=true)
    println("\n\n")
end

samplerun()
dmrg()


text:
################################################################################

"""
    setup(v::Array{T, 2}) where <: Number

Initiate hidden Markov model object

See also: [`process!`](@ref), [`HMM`](@ref), [`HMMParams`](@ref)
"""
function setup!(v::Array{T, 2}) where T <: Number
  noIter = 1
  hmm = HMM([zeros(size(v, 2)) for i in 1:noIter], [zeros(size(v, 1) + 1) for i in 1:noIter])
  reset!(hmm)
  return hmm
end

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

max_stars_repo_name: grabherr/HiddenMarkovModelReaders

max_issues_repo_name: grabherr/HiddenMarkovModelReaders

max_forks_repo_name: grabherr/HiddenMarkovModelReaders


text:
function plot_ks_values_ahmc_vs_mh(ks_res_ahmc::Array{Any},ks_res_mh::Array{Any},n_dim;name="results/ks_test_values_ahmc_mh")

	n_diff_dims = length(n_dim)
	n_functions = length(create_testfunction_for_dim(2,hmc_safe=true))*2

	x = Array{Array{Float64}}(undef,n_diff_dims*n_functions)
	y = Array{Array{Float64}}(undef,n_diff_dims*n_functions)
	x_dim_ticks_t1=Array{Float64}(undef,n_diff_dims*n_functions)
	x_dim_ticks_t2=Array{String}(undef,n_diff_dims*n_functions)

	xy_func_annot=Array{Array{Float64}}(undef,n_functions)
	func_annot = ["HMC","MH"]

	res_ks = [ks_res_ahmc,ks_res_mh]

	for i_function in 1:n_functions
		results = res_ks[i_function]
		for i_dim in 1:n_diff_dims
			index = i_dim+(i_function-1)*n_diff_dims
			y[index] = results[i_dim,2]
			x[index] = ones(length(y[index]))*(i_dim+(i_function-1)*(n_diff_dims+2))
			x_dim_ticks_t1[index]= (i_dim+(i_function-1)*(n_diff_dims+2))
			x_dim_ticks_t2[index]= string(n_dim[i_dim])
			if(i_dim == Int(floor((n_diff_dims/2)))+1 )
				xy_func_annot[i_function] = [(i_dim+(i_function-1)*(n_diff_dims+2))-0.4,-0.11]
			end
		end
		end

	scatter(x,y,label="",margin=5Plots.mm)
	for i in 1:length(xy_func_annot)
		annotate!(xy_func_annot[i]...,text(func_annot[i],10))
	end
	xticks!((x_dim_ticks_t1,x_dim_ticks_t2))
	savefig(name)
	savefig(string(name,"_ND.pdf"))
end

function plot_ks_values(results::Array{Any},n_dim;name="results/ks_test_values")
    n_diff_dims = length(n_dim)
    n_functions = length(create_testfunction_for_dim(2))

    x = Array{Array{Float64}}(undef,n_diff_dims*n_functions)
    y = Array{Array{Float64}}(undef,n_diff_dims*n_functions)
    x_dim_ticks_t1=Array{Float64}(undef,n_diff_dims*n_functions)
    x_dim_ticks_t2=Array{String}(undef,n_diff_dims*n_functions)

    xy_func_annot=Array{Array{Float64}}(undef,n_functions)
    func_annot = collect(keys(create_testfunction_for_dim(2)))

    for i_function in 1:n_functions
        for i_dim in 1:n_diff_dims
            index = i_dim+(i_function-1)*n_diff_dims
            y[index] = results[i_dim,1+i_function]
            x[index] = ones(length(y[index]))*(i_dim+(i_function-1)*(n_diff_dims+2))
            x_dim_ticks_t1[index]= (i_dim+(i_function-1)*(n_diff_dims+2))
            x_dim_ticks_t2[index]= string(n_dim[i_dim])
            if(i_dim == Int(floor((n_diff_dims/2)))+1 )
                xy_func_annot[i_function] = [(i_dim+(i_function-1)*(n_diff_dims+2)),-0.11]
            end
        end
    end

    scatter(x,y,label="",margin=5Plots.mm)
    for i in 1:length(xy_func_annot)
        annotate!(xy_func_annot[i]...,text(func_annot[i],10))
    end
    xticks!((x_dim_ticks_t1,x_dim_ticks_t2))
    savefig(name)
    savefig(string(name,"_ND.pdf"))
end

function plot_ahmi_values(results::Array{Any},n_dim;name="results/ahmi_values")
    n_diff_dims = length(n_dim)
    n_functions = length(create_testfunction_for_dim(2))

    x = Array{Float64}(undef,n_diff_dims*n_functions)
    y = Array{Float64}(undef,n_diff_dims*n_functions)
    yerr = Array{Float64}(undef,n_diff_dims*n_functions)

    x_ticks_t1=Array{Float64}(undef,n_diff_dims)
    x_ticks_t2=Array{String}(undef,n_diff_dims)
    y_ticks_t1=Array{Float64}(undef,n_functions*3)
    y_ticks_t2=Array{String}(undef,n_functions*3)
    xy_func_annot=Array{Array{Float64}}(undef,n_functions)
    func_annot = collect(keys(create_testfunction_for_dim(2)))

    int_plot_range = 0.2

    for i_function in 1:n_functions
     for i_dim in 1:n_diff_dims
         index = i_dim+(i_function-1)*n_diff_dims
         x[index] = i_dim
    	 y[index] = ((i_function-1)*(1))+(1-results[i_dim,1+i_function][1])
    	 yerr[index] = abs(results[i_dim,1+i_function][2])
         x_ticks_t1[i_dim] = i_dim
         x_ticks_t2[i_dim] = string(n_dim[i_dim])
     end

     y_ticks_t1[(i_function-1)*3+1]= ((i_function-1)*(1))-int_plot_range
     y_ticks_t1[(i_function-1)*3+2]= ((i_function-1)*(1))
     y_ticks_t1[(i_function-1)*3+3]= ((i_function-1)*(1))+int_plot_range
     y_ticks_t2[(i_function-1)*3+1]= string(1.0-int_plot_range)
     y_ticks_t2[(i_function-1)*3+2]= string(1.0)
     y_ticks_t2[(i_function-1)*3+3]= string(1.0+int_plot_range)
     xy_func_annot[i_function] = [-0.30,((i_function-1)*(1))]
    end

    for i in 1:n_functions
    	i_start = n_diff_dims*(i-1)+1
    	i_end = n_diff_dims*i
        if i == 1
            plot(x[i_start:i_end],y[i_start:i_end],ribbon=yerr[i_start:i_end],label="",marker='x',margin=11Plots.mm)
        else
            plot!(x[i_start:i_end],y[i_start:i_end],ribbon=yerr[i_start:i_end],label="",marker='x',margin=11Plots.mm)
        end
    end

    for i in 1:length(xy_func_annot)
        annotate!(xy_func_annot[i]...,text(func_annot[i],9),annotation_clip=false)
    end
    #print(xy_func_annot)
    ylims!(-0.5,maximum(y)+0.3)
    xticks!((x_ticks_t1,x_ticks_t2))
    yticks!((y_ticks_t1,y_ticks_t2))
    xlabel!("Number of dimension")
    savefig(name)
    savefig(string(name,"_ND.pdf"))
end

function plot_time_values(results::Array{Any},n_dim;name="results/time_values")
    n_diff_dims = length(n_dim)
    n_functions = length(create_testfunction_for_dim(2))

    x = Array{Float64}(undef,n_diff_dims*n_functions)
    y = Array{Float64}(undef,n_diff_dims*n_functions)
    yerr = Array{Float64}(undef,n_diff_dims*n_functions)

    x_ticks_t1=Array{Float64}(undef,n_diff_dims)
    x_ticks_t2=Array{String}(undef,n_diff_dims)

    func_annot = collect(keys(create_testfunction_for_dim(2)))

    int_plot_range = 0.2

    for i_function in 1:n_functions
         for i_dim in 1:n_diff_dims
             index = i_dim+(i_function-1)*n_diff_dims
             x[index] = i_dim
        	 y[index] = results[i_dim,1+i_function][1]
        	 yerr[index] = abs(results[i_dim,1+i_function][2])
             x_ticks_t1[i_dim] = i_dim
             x_ticks_t2[i_dim] = string(n_dim[i_dim])
         end
    end

    for i in 1:n_functions
    	i_start = n_diff_dims*(i-1)+1
    	i_end = n_diff_dims*i
        if i == 1
            plot(x[i_start:i_end],y[i_start:i_end],ribbon=yerr[i_start:i_end],label=func_annot[i],marker='x',margin=11Plots.mm)
        else
            plot!(x[i_start:i_end],y[i_start:i_end],ribbon=yerr[i_start:i_end],label=func_annot[i],marker='x',margin=11Plots.mm)
        end
    end

    xticks!((x_ticks_t1,x_ticks_t2))
    xlabel!("Number of dimension")
    ylabel!("time in s")
    savefig(name)
    savefig(string(name,"_ND.pdf"))
end

function create_testfunction_for_dim(i_dim::Integer,maxdim=10;hmc_safe=false)
    #sig = Matrix{Float64}([1.5^2 1.5*2.5*0.4 ; 1.5*2.5*0.4 2.5^2])
    sig = Matrix{Float64}(undef,i_dim,i_dim)
    means = Vector{Float64}(undef,i_dim)

    for j in 1:1:i_dim
        means[j] = 5+5*j
        sig[j,j] = (4*(1+j/maxdim))^2 #
    end
    for j in 1:1:i_dim
        for k in 1:1:j
            sig[j,k] = 0.2*(1+i_dim/10)*sqrt(sig[j,j])*sqrt(sig[k,k])
            sig[k,j] = sig[j,k]
        end
    end
    normal = MvNormal(means,sig)
    cauchy = BAT.MultimodalCauchy(µ=5. * i_dim, σ=0.2 * i_dim,n=i_dim)
    funnel = BAT.FunnelDistribution(a=0.5, b=0.5, n=i_dim)
	if hmc_safe
		return Dict(
	        "funnel" => funnel
	    )
	end
    return Dict(
        "normal" => normal,
        "cauchy" => cauchy,
        "funnel" => funnel
    )
end

function run_ND_benchmark(;
    n_dim = 2:2:20,
    algorithm = MetropolisHastings(),
    n_chains = 4,
    n_steps = 4*10^5,
    time_benchmark = true,
    ks_test_benchmark = true,
    ahmi_benchmark = true,
	hmc_benchmark=false)

    ahmi = Matrix{Any}(undef,length(n_dim),length(create_testfunction_for_dim(2,hmc_safe=hmc_benchmark))+1)
    ks_test = Matrix{Any}(undef,length(n_dim),length(create_testfunction_for_dim(2,hmc_safe=hmc_benchmark))+1)
    times = Matrix{Any}(undef,length(n_dim),length(create_testfunction_for_dim(2,hmc_safe=hmc_benchmark))+1)

    convergence = BrooksGelmanConvergence(
        threshold = 1.6, #Change up
        corrected = false
    )
    init = MCMCChainPoolInit(
        init_tries_per_chain = 64..1024,
        nsteps_init = 2500,
    )
    burnin = MCMCMultiCycleBurnin(
        nsteps_per_cycle = 20000,
        max_ncycles = 300
    )

    for i in 1:length(n_dim)
		print(i)
        i_dim = n_dim[i]
        ahmi[i,1] = i_dim
        ks_test[i,1] = i_dim
        times[i,1] = i_dim

        testfunctions = create_testfunction_for_dim(i_dim,n_dim[end],hmc_safe=hmc_benchmark)

        for j in 1:length(testfunctions)

            dis = testfunctions[collect(keys(testfunctions))[j]]
            iid_sample = bat_sample(dis, MCMCSampling(sampler = MetropolisHastings(), trafo = NoDensityTransform(), nchains = nchains, nsteps = nsteps)).result

            mcmc_sample = nothing
            tbf = time()
            if isa(algorithm,BAT.MetropolisHastings)
                mcmc_sample = bat_sample(
                    dis,
                    MCMCSampling(
                        sampler = algorithm,
                        trafo = NoDensityTransform(),
                        nchains = n_chains,
                        nsteps = n_steps,
                        init = init,
                        burnin = burnin,
                        convergence = convergence,
                        strict = true,
                        nonzero_weights = true
                    )
                ).result
            elseif isa(algorithm,BAT.HamiltonianMC)
                mcmc_sample = bat_sample(
                    dis,
                    MCMCSampling(sampler = algorithm, trafo = NoDensityTransform(), nchains = n_chains, nsteps = n_steps)
                ).result
            end
            taf = time()

            if time_benchmark
                time_per_run = Array{Float64}(undef,5)
                for i in 1:5
                    tbf = time()
                    if isa(algorithm,BAT.MetropolisHastings)
                        bat_sample(
                            dis,
                            MCMCSampling(
                                sampler = algorithm,
                                trafo = NoDensityTransform(),
                                nchains = n_chains,
                                nsteps = n_steps,
                                init = init,
                                burnin = burnin,
                                convergence = convergence,
                                strict = true,
                                nonzero_weights = true
                            )
                        ).result
                    elseif isa(algorithm,BAT.HamiltonianMC)
                        bat_sample(
                            dis,
                            MCMCSampling(sampler = algorithm, trafo = NoDensityTransform(), nchains = n_chains, nsteps = n_steps)
                        ).result
                    end
                    taf = time()
                    time_per_run[i] = taf-tbf
                end
                times[i,j+1] = [round(mean(time_per_run)),round(std(time_per_run))]
            else
                times[i,j+1] = round(taf-tbf)
            end

            if ahmi_benchmark
                integral = try
                    bat_integrate(mcmc_sample,AHMIntegration()).result
                catch ex
                      NaN
                end
                if !isnan(integral)
                    ahmi[i,j+1] =  [round(integral.val,digits=3),round(integral.err,digits=3)]
                else
                    ahmi[i,j+1] =  [NaN,NaN]
                end
            end

            if ks_test_benchmark
                ks_test[i,j+1] = round.(bat_compare(iid_sample, mcmc_sample).result.ks_p_values,digits=3)
            end

        end
    end

    header = ["n_dims",collect(keys(create_testfunction_for_dim(2,hmc_safe=hmc_benchmark)))...]
    tf = open("results/times_ND.txt","w")
    pretty_table(tf,times,header)
    close(tf)
    tfl = open("results/times_ND.tex","w")
    pretty_table(tfl,times,header,backend=:latex)
    close(tfl)

    if(ks_test_benchmark && !(hmc_benchmark))
        header = ["n_dims",collect(keys(create_testfunction_for_dim(2,hmc_safe=hmc_benchmark)))...]
        fks = open("results/ks_test_ND.txt","w")
        pretty_table(fks,ks_test,header)
        close(fks)
        fksl = open("results/ks_test_ND.tex","w")
        pretty_table(fksl,ks_test,header,backend=:latex)
        close(fksl)
        plot_ks_values(ks_test,n_dim)
     end
    if(time_benchmark)
        plot_time_values(times,n_dim)
    end
    if(ahmi_benchmark)
        ahmif = open("results/ahmi_ND.txt","w")
        pretty_table(ahmif,ahmi,header)
        close(ahmif)
        ahmifl = open("results/ahmi_ND.tex","w")
        pretty_table(ahmifl,ahmi,header,backend=:latex)
        close(ahmifl)
        plot_ahmi_values(ahmi,n_dim)
    end

    return [ks_test,ahmi,times]
end

function run_ks_ahmc_vs_mh(;n_dim=20:5:35,n_steps=2*10^5, n_chains=4)
    ks_res_ahmc = run_ND_benchmark(n_dim=n_dim,algorithm=HamiltonianMC(), n_steps=n_steps, n_chains=n_chains, time_benchmark=false,ahmi_benchmark=false,hmc_benchmark=true)[1]
    ks_res_mh = run_ND_benchmark(n_dim=n_dim,algorithm=MetropolisHastings(), n_steps=n_steps, n_chains=n_chains, time_benchmark=false,ahmi_benchmark=false,hmc_benchmark=true)[1]
    plot_ks_values_ahmc_vs_mh(ks_res_ahmc,ks_res_mh,n_dim)
end


text:
"""
T is for dispatch.
"""

struct Token{T}
    lineno :: Int64
    colno  :: Int32
    offset :: Int64
    str    :: String
    span   :: Int32
end


text:
using Test, IonSim
using Suppressor

@testset "vibrational_modes -- VibrationalMode" begin
    # setup system
    vm = VibrationalMode(1, [1, 1], δν=1)

    # make sure δν can appropriately handle constant or function
    t = 0:1:100
    ones = [1 for _ in t]
    @test vm.δν.(t) == ones
    @test vm._cnst_δν
    vm.δν = sin
    @test vm.δν.(t) == sin.(t)
    @test !vm._cnst_δν
    vm = VibrationalMode(1, [1, 1], δν=t->t)
    @test !vm._cnst_δν

    # test ==
    vm1 = VibrationalMode(1, [1, 1], δν=1)
    @test vm == vm1

    # test that updating cutoff also appropriately updates shape field
    vm.N = 20
    @test vm.shape == [21]

    # make sure user can't directly update shape field
    vm.shape = 100
    @test vm.shape == [21]
    @test vm.N == 20

    # make sure user can't set field N to a float or negative value
    @test_throws AssertionError vm.N = 1.0
    @test_throws AssertionError vm.N = -1

    # just make sure no errors on print
    @suppress print(vm)
    @suppress show(vm)

    # make sure user can't directly alter :mode_structure or :axis after initialization
    ms = vm.mode_structure
    vm.mode_structure = "something"
    @test vm.mode_structure == ms
    axis = vm.axis
    vm.axis = "something"
    @test vm.axis == axis

    # test indexing of VibrationalMode, which should return a fockstate
    vm.N = 1
    @test vm[0].data == ComplexF64[1; 0]
    @test vm[1].data == ComplexF64[0; 1]
end


text:
export dΦ
export dΦξ

const ϵ = 1e-7

dΦ(Φ, p) = Φ(p .+ ϵ) - Φ(p)
dΦξ(Φ, v, p) = dot(dΦ(Φ, p), v)


text:
using Zygote: @adjoint
using Zygote
using ForwardDiff

val(x::ForwardDiff.Dual) = ForwardDiff.value(x)
val(x::Float64) = x
val(x::Float32) = x

function dconv(x, a, b)
    size=Base.size
    n = length(size(a))
    x = reshape(vec(x), (size(a) .+ size(b) .- ones(Int, n))...)
    r = (
        DSP.conv(x, reverse(b))[(
            # DSP.xcorr(x, b, padmode = :none)[(
            i:j for (i, j) in zip(size(b), size(x))
        )...],
        DSP.conv(x, reverse(a))[(
            # DSP.xcorr(x, a, padmode = :none)[(
            i:j for (i, j) in zip(size(a), size(x))
        )...],
    )
    return r
end
# @adjoint DSP.conv(a, b) = DSP.conv(a, b), x -> dconv(x, a, b)
@adjoint DSP.conv(a, b) = DSP.conv(val.(a), val.(b)),
x -> dconv(val.(x), val.(a), val.(b))


text:
struct EnsembleProblem{T,T2,T3,T4,T5} <: AbstractEnsembleProblem
  prob::T
  prob_func::T2
  output_func::T3
  reduction::T4
  u_init::T5
end

EnsembleProblem(prob;
    output_func = (sol,i)-> (sol,false),
    prob_func= (prob,i,repeat)->prob,
    reduction = (u,data,I)->(append!(u,data),false),
    u_init = []) =
    EnsembleProblem(prob,prob_func,output_func,reduction,u_init)

EnsembleProblem(;prob,
    output_func = (sol,i)-> (sol,false),
    prob_func= (prob,i,repeat)->prob,
    reduction = (u,data,I)->(append!(u,data),false),
    u_init = [], p = nothing) =
    EnsembleProblem(prob,prob_func,output_func,reduction,u_init)