text:
module Backup

using ..Ahorn, Maple
using Gtk, Gtk.ShortNames
using Dates

const timestampFormat = "YYYY-mm-dd HH-MM-SS"

backupTimer = nothing

function getTimestamp()
    return Dates.format(Dates.now(), timestampFormat)
end

function getFilename(side::Side)
    return joinpath(getSideName(side), getTimestamp() * ".bin")
end

function getBackupFilename(side::Side)
    return joinpath(getSideName(side), "LATEST.bin")
end

function pruneBackups(folder::String, side::Side, keep::Integer=10)
    files = Tuple{String, DateTime}[]
    targetFolder = joinpath(folder, getSideName(side))

    for fn in readdir(targetFolder)
        try
            timestamp = DateTime(splitext(fn)[1], timestampFormat)
            push!(files, (fn, timestamp))

        catch ArgumentError
            # Ignore, file with bad format
            # Could be "LATEST" file or user manual backups
        end
    end

    sort!(files, rev=true, by=v -> v[2])
    deleting = getindex.(files[keep + 1:end], 1)

    rm.(joinpath.(Ref(targetFolder), deleting))
end

function backup(side::Side, prune::Bool=true)
    backupFolder = get(Ahorn.config, "backup_directory", joinpath(Ahorn.storageDirectory, "Backups"))
    keep = get(Ahorn.config, "backup_maximum_files", 10)

    filename = getFilename(side)
    latestFilename = getBackupFilename(side)

    outFilename = joinpath(backupFolder, filename)
    outLatestFilename = joinpath(backupFolder, latestFilename)

    mkpath(dirname(outFilename))

    encodeSide(side, outFilename)
    encodeSide(side, outLatestFilename)

    if prune
        pruneBackups(backupFolder, side, keep)
    end
end

function startAutomaticBackup()
    backupRate = get(Ahorn.config, "backup_rate", 60)

    global backupTimer = Timer(0, interval=backupRate) do timer
        side = Ahorn.loadedState.side

        @Ahorn.catchall begin
            if side !== nothing
                backup(side)
            end
        end
    end
end

function openBackupDialog()
    side = Ahorn.loadedState.side
    backupFolder = get(Ahorn.config, "backup_directory", joinpath(Ahorn.storageDirectory, "Backups"))

    if side === nothing
        Ahorn.showFileOpenDialog(MenuItem(), backupFolder)
        Ahorn.loadedState.filename = ""
        Ahorn.persistence["files_lastfile"] = ""

    else
        latestFilename = Ahorn.loadedState.filename
        latestBackupFilename = getBackupFilename(side)

        # The latest filename may also be in another folder, make sure we get the right folder
        Ahorn.showFileOpenDialog(MenuItem(), dirname(joinpath(backupFolder, latestBackupFilename)))
        Ahorn.loadedState.filename = latestFilename
        Ahorn.persistence["files_lastfile"] = latestFilename
    end
end

function initBackup(persistence::Ahorn.Config)
    active = get(Ahorn.config, "backup_enabled", true)
    showDialog = get(Ahorn.config, "backup_dialog_enabled", true)

    if active && showDialog && get(persistence, "currently_running", false)
        @async begin
            if ask_dialog("Ahorn did not shut down properly.\nDo you want to restore from automatic backup?", Ahorn.window)
                openBackupDialog()
            end

            startAutomaticBackup()
        end

    else
        startAutomaticBackup()
    end

    persistence["currently_running"] = true
    Ahorn.saveConfig(persistence, true)
end

end


text:
using BenchmarkTools, Random
using Mill, Flux

using Plots, LaTeXStrings; pgfplotsx()

Random.seed!(42)
params = Base.Iterators.product([50, 100, 300], [100, 200, 300]) |> collect |> vec
for (m, bl) in params
    bags = length2bags(rand(0:bl, 20))
    A = rand(Float32, m, sum(length.(bags)))
    @show m, bl
    o = SegmentedMeanMax(m)
    @btime $o($A, $bags)
end

xticklabels = [latexstring("$a|$b") for (a,b) in params]

t1 = [37.764, 61.558, 144.835, 65.729, 88.520, 199.147, 91.258, 134.736, 433.920]
l1 = "Mill 1.5.2"
t2 = [38.907, 63.474, 146.922, 66.963, 89.890, 203.901, 89.999, 134.687, 424.865]
l2 = "Mill 2.0: bagcount on"
t3 = [36.856, 58.565, 138.007, 61.036, 80.387, 187.378, 89.691, 139.382, 427.770]
l3 = "Mill 2.0: bagcount off"
groupedbar(xticklabels, hcat(t1, t2, t3); ylabel=L"\mu s", xlabel="instance size | length of 20 bags", labels=[l1 l2 l3], title="MeanMax Aggregation")

savefig("aggregation.png")

Random.seed!(42)
params = Base.Iterators.product([50, 100, 300], [100, 200, 300]) |> collect |> vec
for (m, bl) in params
    bags = length2bags(rand(0:bl, 20))
    A = rand(Float32, m, sum(length.(bags)))
    @show m, bl
    o = SegmentedMeanMax(m)
    @btime gradient((o, A) -> o(A, $bags) |> sum, $o, $A)
end

xticklabels = [latexstring("$a|$b") for (a,b) in params]

t1 = [217.039, 387.292, 1203, 357.977, 651.138, 1985, 518.316, 1120, 3796]
l1 = "Mill 1.5.2"
t2 = [217.585, 439.348, 1379, 371.837, 646.035, 1862, 459.842, 941.553, 3609]
l2 = "Mill 2.0: bagcount on"
t3 = [193.479, 398.941, 1144, 322.763, 590.045, 1627, 454.366, 1031, 3651]
l3 = "Mill 2.0: bagcount off"
groupedbar(xticklabels, hcat(t1, t2, t3); ylabel=L"\mu s", xlabel="instance size | length of 20 bags", labels=[l1 l2 l3], title="MeanMax Aggregation (gradients)")

savefig("aggregation_grad.png")

ext: jl


text:
using AMQPClient
const VIRTUALHOST = "/"
const HOST = "127.0.0.1"

function send()
    # 1. Create a connection to the localhost or 127.0.0.1 of virtualhost '/'
    connection(; virtualhost=VIRTUALHOST, host=HOST) do conn
        # 2. Create a channel to send messages
        channel(conn, AMQPClient.UNUSED_CHANNEL, true) do chan
            # 3. Declare Exchange
            exchange = "topic_logs"
            exchange_topic = "topic"
            exchange_declare(chan, exchange, EXCHANGE_TYPE_TOPIC)
            # 4. Get input data
            if length(Base.ARGS) >= 2
                routing_key = Base.ARGS[1]
                received = join(Base.ARGS[2:end], ' ')
            else
                routing_key = "info"
                received = "Hello World"
            end
            # 5. Prepare and send data
            data = convert(Vector{UInt8}, codeunits(received))
            msg = Message(data, content_type="text/plain", delivery_mode=PERSISTENT)
            basic_publish(chan, msg; exchange=exchange, routing_key=routing_key)
            println("Message sent: $received, routing key: $routing_key")
        end
    end
end

send()

max_forks_repo_name: tomashq/rabbitmq-tutorials


text:
using EulerCharacteristic
using XUnit

@testset "Test euler characteristics for known sets"  begin include("calculation.jl") end
@testset "Test euler characteristic tracker" begin include("tracker.jl") end


text:
## Primitive polynomials bases (PolynomialBasis)

@doc raw"""
    struct PolynomialBasis

Wrapper type for the `helfem::atomic::polynomial_basis::PolynomialBasis` class, representing
a set of polynomials ``\{p_i(x)\}`` on a domain ``x \in [-1, 1]``.

# Constructors

    PolynomialBasis(basis::Symbol, nnodes)

Constructs a particular polynomial basis with a specific number of defining nodes.
`basis` must be one of: `:hermite`, `:legendre` or `:lip`.
"""
struct PolynomialBasis
    pb :: helfem.PolynomialBasis
    primbas :: Int
    nnodes :: Int

    function PolynomialBasis(basis::Symbol, nnodes::Integer)
        primbas = (basis == :hermite) ? 2 :
                  (basis == :legendre) ? 3 :
                  (basis == :lip) ? 4 : error("Invalid primitive basis name $basis")
        new(helfem.polynomial_basis(primbas, nnodes), primbas, nnodes)
    end
end

primbas_name(primbas::Int) =
    (primbas == 2) ? "HermiteBasis" :
    (primbas == 3) ? "LegendreBasis" :
    (primbas == 4) ? "LIPBasis" : error("Invalid primbas value $primbas")

primbas_name(pb::PolynomialBasis) = primbas_name(pb.primbas)

function Base.show(io::IO, pb::PolynomialBasis)
    classname = primbas_name(pb.primbas)
    order = helfem.pb_order(pb.pb)
    nbf = length(pb)
    write(io, "PolynomialBasis($(primbas_name(pb.primbas)), order=$order) with $(nbf) basis functions")
end

"""
    length(pb::PolynomialBasis) -> Integer

Returns the number of basis functions (polynomials) represented by the given polynomial basis.
"""
Base.length(pb::PolynomialBasis) = helfem.get_nbf(pb.pb)

"""
    HelFEM.nnodes(pb::PolynomialBasis) -> Integer

Returns the number of control nodes used to define the polynomial basis in the constructor.
"""
nnodes(pb::PolynomialBasis) = pb.nnodes

function (pb::PolynomialBasis)(xs)
    arma_xs = helfem.ArmaVector(collect(xs))
    return collect(helfem.pb_eval(pb.pb, arma_xs))
end

function pb_eval(pb::PolynomialBasis, xs::AbstractVector)
    npts, nbf = length(xs), length(pb)
    fs, dfs = helfem.ArmaMatrix(npts, nbf), helfem.ArmaMatrix(npts, nbf)
    helfem.pb_eval(pb.pb, helfem.ArmaVector(xs), fs, dfs)
    (fs = collect(fs), dfs = collect(dfs))
end

function controlpoints(pb::PolynomialBasis)
    # We only know how to calculate control points for the LIP basis currently
    pb.primbas == 4 || throw(ArgumentError("Can only calculate control points for LIP basis"))
    # The LIPBasis uses Gauss-Lobatto nodes, so the following only applies to LIP
    xs, _ = lobatto(pb.nnodes) # we only need the nodes, can discard the weights
    return xs
end


text:
@enum OBJECT begin
    INVALID_OBJ = -1
    INVALID_IDX = -2
    PLAYER = 0
    SPACE = 100
    WALL_OBJ = 51
    EXIT_OBJ = 52
    AGENT_OBJ = 53
    EXPLOSION = 54
    EXPLOSION2 = 55
    EXPLOSION3 =56
    EXPLOSION4= 57
    EXPLOSION5 = 58
    TRAIL = 59
    DOOR_OBJ = 200
    KEY_OBJ = 300
end


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#    API implementation for J4 orbit propagator.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

get_epoch(orbp::OrbitPropagatorJ4) = orbp.j4d.epoch

function get_mean_elements(orbp::OrbitPropagatorJ4)
    j4d = orbp.j4d
    orb = KeplerianElements(j4d.epoch + j4d.Δt/86400,
                            j4d.al_k * j4d.j4_gc.R0,
                            j4d.e_k,
                            j4d.i_k,
                            j4d.Ω_k,
                            j4d.ω_k,
                            j4d.f_k)

    return orb
end

"""
    init_orbit_propagator(Val(:J4), epoch::Number, a_0::Number, e_0::Number, i_0::Number, Ω_0::Number, ω_0::Number, f_0::Number, dn_o2::Number = 0, ddn_o6::Number = 0; j4_gc::J4_GravCte{T} = j4_gc_egm08) where T
    init_orbit_propagator(Val(:J4), orb_0::Orbit, dn_o2::Number = 0, ddn_o6::Number = 0; j4_gc::J4_GravCte = j4_gc_egm08)

Initialize the J4 orbit propagator.

# Args

* `epoch`: Initial orbit epoch [Julian Day].
* `a_0`: Initial mean semi-major axis [m].
* `e_0`: Initial mean eccentricity.
* `i_0`: Initial mean inclination [rad].
* `Ω_0`: Initial mean right ascension of the ascending node [rad].
* `ω_0`: Initial mean argument of perigee [rad].
* `f_0`: Initial mean true anomaly [rad].
* `dn_o2`: (OPTIONAL) First time derivative of mean motion divided by 2
           \\[rad/s²] (**Default** = 0).
* `ddn_o6`: (OPTIONAL) Second time derivative of mean motion divided by 6
            \\[rad/s³] (**Default** = 0).
* `orb_0`: Instance of the structure `KeplerianElements` with the initial mean
           orbital elements [SI].

## Keywords

* `j4_gc`: (OPTIONAL) J4 orbit propagator gravitational constants
           (**Default** = `j4_gc_egm08`).

"""
function init_orbit_propagator(::Val{:J4}, epoch::Number, a_0::Number,
                               e_0::Number, i_0::Number, Ω_0::Number,
                               ω_0::Number, f_0::Number, dn_o2::Number = 0,
                               ddn_o6::Number = 0;
                               j4_gc::J4_GravCte{T} = j4_gc_egm08) where T

    # Create the new J4 propagator structure.
    j4d = j4_init(epoch, a_0, e_0, i_0, Ω_0, ω_0, f_0, dn_o2, ddn_o6,
                  j4_gc = j4_gc)

    # Create and return the orbit propagator structure.
    return OrbitPropagatorJ4(j4d)
end

function init_orbit_propagator(::Val{:J4},
                               orb_0::Orbit,
                               dn_o2::Number = 0,
                               ddn_o6::Number = 0;
                               j4_gc::J4_GravCte = j4_gc_egm08)
    # Convert the orbit representation to Keplerian elements.
    k_0 = convert(KeplerianElements, orb_0)

    return init_orbit_propagator(Val(:J4), k_0.t, k_0.a, k_0.e, k_0.i, k_0.Ω,
                                 k_0.ω, k_0.f, dn_o2, ddn_o6; j4_gc = j4_gc)
end

function propagate!(orbp::OrbitPropagatorJ4{T}, t::Number) where T
    # Auxiliary variables.
    j4d = orbp.j4d

    # Propagate the orbit.
    r_i, v_i = j4!(j4d, t)

    # Return.
    return r_i, v_i
end

function step!(orbp::OrbitPropagatorJ4, Δt::Number)
    # Auxiliary variables.
    j4d = orbp.j4d

    # Propagate the orbit.
    r_i, v_i = j4!(j4d, j4d.Δt + Δt)

    # Return the information about the step.
    return r_i, v_i
end


text:
using SQLite, DataFrames

export read_data, ImgData

struct ImgData<:Database
    filename::String
end
ImgData() = ImgData("imgdata.db")

function Base.getproperty(database::Database, name::Symbol, x)
    if length(String(name)) > 7 && String(name)[end-6:end] == "_kwargs"
        name = Symbol(String(name)[1:end-7])
        field = getfield(database, name)
        x = x == nothing ? Dict([(f, nothing) for f in fieldnames(field)]) : x
        ans = Dict()
        for (key, value) in pairs(x)
            ans[key] = getproperty(field, key, value)
        end
        return ans
    elseif parentmodule(fieldtype(typeof(database), name)) ∉ [Base, Core]
        typename = typeof(getfield(database, name))
        return "$typename"
    else
        return getfield(database, name)
    end
end

function read_data(filename::String, data_ids)
    db = SQLite.DB(filename)
    db = DBInterface.connect(SQLite.DB, filename)
    img_ids = data_ids["image"]
    N = size(img_ids)
    cmdline = "SELECT * FROM image WHERE rowid IN ("
    cmdline *= join(map(id -> string(id), img_ids), ", ") * ")"
    ds = DBInterface.execute(db, cmdline) |> DataFrame
    data = Dict()
    coords = hcat(map(crd -> reinterpret(Float64, crd), ds.coord)...)
    gradients = hcat(map(grd -> reinterpret(Float64, grd), ds.gradients)...)
    if typeof(ds.potentials) != Array{Missing, 1}
    potentials = hcat(map(ptls -> reinterpret(Float64, ptls), ds.potentials)...)
    data["potentials"] = potentials # A x N ;
    end
    data["coords"] = coords        # D x N ; confirmed
    data["gradients"] = gradients  # D x N ; confirmed
    data["potential"] = Array{Float64, 1}(reinterpret(Float64, ds.potential))
    return data
end

include("./manybody.jl")

max_stars_repo_name: schinavro/taps

max_issues_repo_name: schinavro/taps

max_forks_repo_name: schinavro/taps


text:
using Plots, Test

@testset "Subplot sclicing" begin
    pl = @test_nowarn plot(
        rand(4, 8),
        layout = 4,
        yscale = [:identity :identity :log10 :log10],
    )
    @test pl[1][:yaxis][:scale] == :identity
    @test pl[2][:yaxis][:scale] == :identity
    @test pl[3][:yaxis][:scale] == :log10
    @test pl[4][:yaxis][:scale] == :log10
end

@testset "Plot title" begin
    pl = plot(rand(4, 8), layout = 4, plot_title = "My title")
    @test pl[:plot_title] == "My title"
    @test pl[:plot_titleindex] == 5

    plot!(pl)
    @test pl[:plot_title] == "My title"
    @test pl[:plot_titleindex] == 5

    plot!(pl, plot_title = "My new title")
    @test pl[:plot_title] == "My new title"
    @test pl[:plot_titleindex] == 5
end


text:
"""
to_string
"""
function to_string(data::ModelSelection.ModelSelectionData, result::AllSubsetRegressionResult)
    datanames_index = ModelSelection.create_datanames_index(result.datanames)

    out = ""
    out *= @sprintf("\n")
    out *= @sprintf("══════════════════════════════════════════════════════════════════════════════\n")
    out *= @sprintf("                              Best model results                              \n")
    out *= @sprintf("══════════════════════════════════════════════════════════════════════════════\n")
    out *= @sprintf("                                                                              \n")
    out *= @sprintf("                                     Dependent variable: %s                   \n", data.depvar)
    out *= @sprintf("                                     ─────────────────────────────────────────\n")
    out *= @sprintf("                                                                              \n")
    out *= @sprintf(" Selected covariates                 Coef.")
    if result.ttest
        out *= @sprintf("        Std.         t-test")
    end
    out *= @sprintf("\n")
    out *= @sprintf("──────────────────────────────────────────────────────────────────────────────\n")

    cols = ModelSelection.get_selected_variables(Int64(result.bestresult_data[datanames_index[:index]]), data.expvars, data.intercept)

    for pos in cols
        varname = data.expvars[pos]
        out *= @sprintf(" %-35s", varname)
        out *= @sprintf(" %-10f", result.bestresult_data[datanames_index[Symbol(string(varname, "_b"))]])
        if result.ttest
            out *= @sprintf("   %-10f", result.bestresult_data[datanames_index[Symbol(string(varname, "_bstd"))]])
            out *= @sprintf("   %-10f", result.bestresult_data[datanames_index[Symbol(string(varname, "_t"))]])
        end
        out *= @sprintf("\n")
    end

    out *= @sprintf("──────────────────────────────────────────────────────────────────────────────\n")
    out *= @sprintf(" Observations                        %-10d\n", result.bestresult_data[datanames_index[:nobs]])
    out *= @sprintf(" F-statistic                         %-10f\n", result.bestresult_data[datanames_index[:F]])
    if :r2adj in result.datanames
        out *= @sprintf(" Adjusted R²                         %-10f\n", result.bestresult_data[datanames_index[:r2adj]])
    end
    for criteria in result.criteria
        if AVAILABLE_CRITERIA[criteria]["verbose_show"] && criteria != :r2adj
    out *= @sprintf(" %-30s      %-10f\n", AVAILABLE_CRITERIA[criteria]["verbose_title"], result.bestresult_data[datanames_index[criteria]])
        end
    end

    if !result.modelavg
        out *= @sprintf("──────────────────────────────────────────────────────────────────────────────\n")
    else
        out *= @sprintf("\n")
        out *= @sprintf("\n")
        out *= @sprintf("══════════════════════════════════════════════════════════════════════════════\n")
        out *= @sprintf("                            Model averaging results                           \n")
        out *= @sprintf("══════════════════════════════════════════════════════════════════════════════\n")
        out *= @sprintf("                                                                              \n")
        out *= @sprintf("                                     Dependent variable: %s                   \n", data.depvar)
        out *= @sprintf("                                     ─────────────────────────────────────────\n")
        out *= @sprintf("                                                                              \n")
        out *= @sprintf(" Covariates                          Coef.")
        if result.ttest
            out *= @sprintf("        Std.         t-test")
        end
        out *= @sprintf("\n")
        out *= @sprintf("──────────────────────────────────────────────────────────────────────────────\n")

        for varname in data.expvars
            out *= @sprintf(" %-35s", varname)
            out *= @sprintf(" %-10f", result.modelavg_data[datanames_index[Symbol(string(varname, "_b"))]])
            if result.ttest
                out *= @sprintf("   %-10f", result.modelavg_data[datanames_index[Symbol(string(varname, "_bstd"))]])
                out *= @sprintf("   %-10f", result.modelavg_data[datanames_index[Symbol(string(varname, "_t"))]])
            end
            out *= @sprintf("\n")
        end
        out *= @sprintf("\n")
        out *= @sprintf("──────────────────────────────────────────────────────────────────────────────\n")
        out *= @sprintf(" Observations                        %-10d\n", result.modelavg_data[datanames_index[:nobs]])
        out *= @sprintf(" Adjusted R²                         %-10f\n", result.modelavg_data[datanames_index[:r2adj]])
        out *= @sprintf(" F-statistic                         %-10f\n", result.modelavg_data[datanames_index[:F]])
        out *= @sprintf(" Combined criteria                   %-10f\n", result.modelavg_data[datanames_index[:order]])
        out *= @sprintf("──────────────────────────────────────────────────────────────────────────────\n")

    end

    return out
end


text:
"""
    viz_backprop(img, model; top_k=1, target_based=false, target_index)

Returns gradient generated with Vanilla Backpropagation for `top_k` classes. It’s the original saliency map algorithm for supervised deep learning from [Simonyan et. al. (2013)](https://arxiv.org/abs/1312.6034). These gradients can further be visualized using utilty functions.

# Arguments
 - `img`: source image for which gradients are to be calculated
 - `model`: model to be used for visualization
 - `top_k=1`: gradients are returned for `top_k` classes as per predictions by the model.
 - `target_based=false`: if set to `true`, function returns the maximizing gradient for the specified `target_index` class along with gradient for `top_k` classes.
 - `target_index`: defines the target class, for which the maximizing gradients are to be calculated. Required when `target_based` is set to `true`.

The function returns an of gradients for `top_k=1` predicted classes by default. If `target_based` is set to `true` and `target_index` is specified, the function returns a tuple: `(gradients_for_top_k_classes, gradient_for_target_class)`.
"""
function viz_backprop(img, model; top_k=1, target_based=false, target_index=-1)
    model_ = model[1:end-1] |> gpu
    set_training_false()
    preds, back = Zygote.pullback(model_, img)
    set_training_true()
    probs = softmax(preds)
    prob, inds = get_topk(probs, k=top_k)
    grads = []
    for (i, idx) in enumerate(inds)
      push!(grads, (back(one_hot_encode(preds, idx))[1], prob[i], idx))
    end
    if (target_based==true)
        if target_index == -1
            @warn("No Target Class defined!!")
        else
            return (grads, back(one_hot_encode(preds, target_index)))
        end
    end
    grads
end

# Integrated Backpropagation

"""
    stepped_images(img, steps)

Utility function for `viz_integrated_gradient`. Generate scaled images, ie. images multiplied by a factor `i` defined in the range of `0` to `1.0` over equal intervals of `1/steps`.
"""
function stepped_images(img, steps)
    step_size = 1.0/steps
    imgs_in = [Float32.(img.*i) for i in 0+step_size:step_size:1.0]
    return imgs_in
end

"""
    viz_integrated_gradient(img, model, steps; target_based=false, target_index)

Produces gradients generated with integrated gradients from the image scaled over `steps`. This functions uses the `viz_backprop` for the calculation of gradients for each of scaled image and then returns integrated gradients as the average over all the scaled images.

# Arguments
 - `img`: source image on which the gradient are to be calculated.
 - `model`: model to be used for visualization
 - `steps`: number of steps between 0 to 1.0 for which the image is to be scaled while calculating integrated gradients.
 - `target_based=false`: if set `true`, the integrated gradients are evaluated for the specified `target_index`.
 - `target_index`: index of the target class for which the maximized gradients are to be evaluated. Necessary condition when `target_based` is set as `true`.

"""
function viz_integrated_gradient(img, model, steps; target_based=false, target_index=-1)
    imgs_in = stepped_images(img, steps)
    integrated_gradient = zeros(Float32,size(img)...)
    if (target_based==true)
        for img_ in imgs_in
            integrated_gradient .+= viz_backprop(img_, model; top_k=1, target_based=true, target_index=target_index)[2][1]
        end
    else
        for img_ in imgs_in
            grad = viz_backprop(img_, model; top_k=1, target_based=false)[1][1]
            integrated_gradient .+= grad
        end
    end
    Float32.(integrated_gradient/steps)
end


text:
##  This View checks for the index to be L1 or the mirrored version (L2)
# and then replaces the value by half of the parent at L1
struct FourierSplit{T,N, AA<:AbstractArray{T, N}} <: AbstractArray{T,N}
    parent::AA # holds the data (or is another view)
    D::Int # dimension along which to apply to copy
    L1::Int # low index position to copy from (and half)
    L2::Int # high index positon to copy to (and half)

    # This version below is needed to avoid a split for the firs rft dimension but still return half the value
    # FFTs and other RFT dimension should use the version without L2
    function FourierSplit(parent::AA, D::Int,L1::Int,L2::Int) where {T,N, AA<:AbstractArray{T, N}}
        return new{T,N, AA}(parent, D, L1, L2)
    end
    function FourierSplit(parent::AA, D::Int,L1::Int) where {T,N, AA<:AbstractArray{T, N}}
        mid = fft_center(size(parent)[D])
        L2 = mid + (mid-L1)
        return FourierSplit(parent, D,L1,L2)
    end
end

Base.IndexStyle(::Type{FD}) where {FD<:FourierSplit} = IndexStyle(parenttype(FD))
parenttype(::Type{FourierSplit{T,N,AA}}) where {T,N,AA} = AA
parenttype(A::FourierSplit) = parenttype(typeof(A))
Base.parent(A::FourierSplit) = A.parent
Base.size(A::FourierSplit) = size(parent(A))

@inline function Base.getindex(A::FourierSplit{T,N, <:AbstractArray{T, N}}, i::Vararg{Int,N}) where {T,N}
    if i[A.D]==A.L2 || i[A.D]==A.L1 # index along this dimension A.D corrsponds to slice L2
        # not that "setindex" in the line below modifies only the index, not the array
        @inbounds return parent(A)[Base.setindex(i,A.L1, A.D)...] / 2
    else i[A.D]==A.L2
        @inbounds return parent(A)[i...]
        # @inbounds return parent(A)[i...]
    end
end

## This View checks for the index to be L1
# and then replaces the value by add the value at the mirrored position L2
struct FourierJoin{T,N, AA<:AbstractArray{T, N}} <: AbstractArray{T, N}
    parent::AA
    D::Int # dimension along which to apply to copy
    L1::Int # low index position to copy from (and half)
    L2::Int # high index positon to copy to (and half)

    # This version below is needed to avoid a split for the firs rft dimension but still return half the value
    # FFTs and other RFT dimension should use the version without L2
    function FourierJoin(parent::AA, D::Int, L1::Int, L2::Int) where {T, N, AA<:AbstractArray{T, N}}
        return new{T, N, AA}(parent, D, L1, L2)
    end

    function FourierJoin(parent::AA, D::Int,L1::Int) where {T, N, AA<:AbstractArray{T, N}}
        mid = fft_center(size(parent)[D])
        L2 = mid + (mid-L1)
        return FourierJoin(parent, D, L1, L2)
    end
end
Base.IndexStyle(::Type{FS}) where {FS<:FourierJoin} = IndexStyle(parenttype(FS))
parenttype(::Type{FourierJoin{T,N,AA}}) where {T,N,AA} = AA
parenttype(A::FourierJoin) = parenttype(typeof(A))
Base.parent(A::FourierJoin) = A.parent
Base.size(A::FourierJoin) = size(parent(A))

@inline function Base.getindex(A::FourierJoin{T,N, <:AbstractArray{T, N}}, i::Vararg{Int,N}) where {T,N}
    if i[A.D]==A.L1
        @inbounds return (parent(A)[i...] + parent(A)[Base.setindex(i, A.L2, A.D)...])
    else
        @inbounds return (parent(A)[i...])
    end
end


text:
# Model
include_model("quadruped")

# signed-distance function (box)
function ϕ_func(model::Quadruped, q)
	h = 0.25
	w = 0.5
	x = 1.25

	p_leg_1 = kinematics_2(model, q, body = :calf_1, mode = :ee)
	p_leg_2 = kinematics_2(model, q, body = :calf_2, mode = :ee)
	p_leg_3 = kinematics_3(model, q, body = :calf_3, mode = :ee)
	p_leg_4 = kinematics_3(model, q, body = :calf_4, mode = :ee)

	z1 = p_leg_1[2] - (((p_leg_1[1] > (x - w / 2.0)) && (p_leg_1[1] < (x + w / 2.0))) ? h : 0.0)
	z2 = p_leg_2[2] - (((p_leg_2[1] > (x - w / 2.0)) && (p_leg_2[1] < (x + w / 2.0))) ? h : 0.0)
	z3 = p_leg_3[2] - (((p_leg_3[1] > (x - w / 2.0)) && (p_leg_3[1] < (x + w / 2.0))) ? h : 0.0)
	z4 = p_leg_4[2] - (((p_leg_4[1] > (x - w / 2.0)) && (p_leg_4[1] < (x + w / 2.0))) ? h : 0.0)

	@SVector [z1, z2, z3, z4]
end

include(joinpath(pwd(), "models/visualize.jl"))
vis = Visualizer()
render(vis)
θ = pi / 3.5

q1 = initial_configuration(model, θ)
qT = copy(q1)
qT[1] = 2.5
visualize!(vis, model, [q1])

# Horizon
T = 51

# Time step
tf = 2.5
h = tf / (T - 1)

# Bounds

# control
u1 = initial_torque(model, q1, h)[model.idx_u]

# ul <= u <= uu
_uu = Inf * ones(model.m)
_uu[model.idx_u] = model.uU

_ul = zeros(model.m)
_ul[model.idx_u] = model.uL
ul, uu = control_bounds(model, T, _ul, _uu)

xl, xu = state_bounds(model, T,
    x1 = [q1; q1],
    xT = [qT; qT])

# Objective
include_objective(["velocity", "nonlinear_stage"])
q_ref = linear_interpolation(q1, qT, T)
visualize!(vis, model, q_ref)

x0 = configuration_to_state(q_ref)

# penalty on slack variable
obj_penalty = PenaltyObjective(1.0e5, model.m)

# quadratic tracking objective
# Σ (x - xref)' Q (x - x_ref) + (u - u_ref)' R (u - u_ref)
obj_control = quadratic_tracking_objective(
    [Diagonal(1.0 * ones(model.n)) for t = 1:T],
    [Diagonal([1.0e-5 * ones(model.nu)..., zeros(model.m - model.nu)...]) for t = 1:T-1],
    [x0[end] for t = 1:T],
    [[u1; zeros(model.m - model.nu)] for t = 1:T])

# quadratic velocity penalty
# Σ v' Q v
obj_velocity = velocity_objective(
    [Diagonal(1.0e-1 * ones(model.nq)) for t = 1:T-1],
    model.nq,
    h = h,
    idx_angle = collect([3, 4, 5, 6, 7, 8, 9, 10, 11]))

# torso height
q2_idx = (12:22)
t_h = kinematics_1(model, q1, body = :torso, mode = :com)[2]
function l_stage_torso_h(x, u, t)
    10000.0 * (kinematics_1(model, view(x, q2_idx), body = :torso, mode = :com)[2] - t_h)^2.0
end
l_terminal_torso_h(x) = 0.0
obj_torso_h = nonlinear_stage_objective(l_stage_torso_h, l_terminal_torso_h)

# torso lateral
function l_stage_torso_lat(x, u, t)
    1.0 * (kinematics_1(model, view(x, q2_idx), body = :torso, mode = :com)[1] - kinematics_1(model, view(x0[t], q2_idx), body = :torso, mode = :com)[1])^2.0
end
l_terminal_torso_lat(x) = 0.0
obj_torso_lat = nonlinear_stage_objective(l_stage_torso_lat, l_terminal_torso_lat)

# foot 1 height
function l_stage_fh1(x, u, t)
    10.0 * (kinematics_2(model, view(x, q2_idx), body = :calf_1, mode = :ee)[2] - 0.1)^2.0
end
l_terminal_fh1(x) = 0.0
obj_fh1 = nonlinear_stage_objective(l_stage_fh1, l_terminal_fh1)

# foot 2 height
function l_stage_fh2(x, u, t)
    10.0 * (kinematics_2(model, view(x, q2_idx), body = :calf_2, mode = :ee)[2] - 0.1)^2.0
end
l_terminal_fh2(x) = 0.0
obj_fh2 = nonlinear_stage_objective(l_stage_fh2, l_terminal_fh2)

# foot 3 height
function l_stage_fh3(x, u, t)
    10.0 * (kinematics_3(model, view(x, q2_idx), body = :calf_3, mode = :ee)[2] - 0.1)^2.0
end
l_terminal_fh3(x) = 0.0
obj_fh3 = nonlinear_stage_objective(l_stage_fh3, l_terminal_fh3)

# foot 4 height
function l_stage_fh4(x, u, t)
    10.0 * (kinematics_3(model, view(x, q2_idx), body = :calf_4, mode = :ee)[2] - 0.1)^2.0
end
l_terminal_fh4(x) = 0.0
obj_fh4 = nonlinear_stage_objective(l_stage_fh4, l_terminal_fh4)

obj = MultiObjective([obj_penalty,
                      obj_control,
                      obj_velocity,
                      obj_torso_h,
                      obj_torso_lat,
                      obj_fh1,
                      obj_fh2,
                      obj_fh3,
                      obj_fh4])

# Constraints
include_constraints(["contact"])
con_contact = contact_constraints(model, T)
con = multiple_constraints([con_contact])

# Problem
prob = trajectory_optimization_problem(model,
               obj,
               T,
               h = h,
               xl = xl,
               xu = xu,
               ul = ul,
               uu = uu,
               con = con)

# trajectory initialization
u0 = [[u1; 1.0e-4 * rand(model.m - model.nu)] for t = 1:T-1] # random controls

# Pack trajectories into vector
z0 = pack(x0, u0, prob)
z0 .+= 0.001 * randn(prob.num_var)

# Solve
include_snopt()

@time z̄ , info = solve(prob, copy(z0),
    nlp = :SNOPT7,
    tol = 1.0e-3, c_tol = 1.0e-3,
    time_limit = 60 * 20, mapl = 5)

check_slack(z̄, prob)
x̄, ū = unpack(z̄, prob)

# Visualize
render(vis)
visualize!(vis, model, state_to_configuration(x̄), Δt = h)
setobject!(vis["box"], GeometryBasics.HyperRectangle(Vec(0.0, 0.0, 0.0), Vec(0.5, 1.0, 0.25)))
settransform!(vis["box"], Translation(1.0, -0.5, 0))


text:
"""
    Fermi.Integrals

Module to compute and manage molecular integrals.
"""
module Integrals

using Fermi
using Fermi.Options
using Fermi.Orbitals
using LinearAlgebra
using GaussianBasis
using TensorOperations

import Base: getindex, setindex!, delete!, show

export IntegralHelper, delete!, mo_from_ao!
export BasisSet, BasisFunction

include("ERITypes.jl")

"""
    IntegralHelper

Manager for integrals computation and storage.
Accesss like a dictionary e.g.,
```
julia> ints = Fermi.Integrals.IntegralHelper()
julia> ints["S"] #Returns overlap matrix
```
A key is associated with each type of integral

    "S"           -> Overlap integral
    "T"           -> Electron kinetic energy integral
    "V"           -> Electron-nuclei attraction integral
    "ERI"         -> Electron repulsion integral

# Fields
    molecule                    Molecule object
    orbitals                    Orbitals used in the integral computation
    basis                       Basis set name
    cache                       Holds computed integrals
    eri_type                    Defines how electron repulsion integrals are handled
"""
struct IntegralHelper{T<:AbstractFloat,E<:AbstractERI,O<:AbstractOrbitals}
    molecule::Molecule
    orbitals::O
    basis::String
    cache::Dict{String,AbstractArray{T}}
    eri_type::E
end

function IntegralHelper(x...;k...)

    precision = Options.get("precision")
    if precision == "single"
        IntegralHelper{Float32}(x...; k...)
    elseif precision == "double"
        IntegralHelper{Float64}(x...; k...)
    else
        throw(FermiException("precision can only be `single` or `double`. Got $precision"))
    end
end

function IntegralHelper{T}(;molecule = Molecule(), orbitals = nothing,
                           basis = Options.get("basis"), eri_type=nothing) where T<:AbstractFloat

    if orbitals === nothing
        bset = BasisSet(basis, molecule.atoms)
        orbitals = AtomicOrbitals(bset)
    end

    # If the associated orbitals are AtomicOrbitals and DF is requested, JKFIT is set by default
    if Options.get("df") && orbitals isa AtomicOrbitals && eri_type === nothing
        eri_type = JKFIT(molecule)

    # If df is requested, but the orbitals are not AtomicOrbitals then RIFIT is set by default
    elseif Options.get("df") && eri_type === nothing

        eri_type = RIFIT(molecule)

    # Else, if eri_type is not an AbstractERI object, Chonky is used. This will overwrite invalid entries for eri_type
    elseif !(typeof(eri_type) <: AbstractERI)
        eri_type = orbitals isa AtomicOrbitals ? SparseERI() : Chonky()
    end

    cache = Dict{String, AbstractArray{T}}()
    IntegralHelper(molecule, orbitals, basis, cache, eri_type)
end

function IntegralHelper{T}(bset::BasisSet, eri_type=nothing) where T<:AbstractFloat
    IntegralHelper{T}(
        molecule = bset.molecule,
        orbitals = AtomicOrbitals(bset),
        basis = bset.basis_name,
        eri_type=eri_type
    )
end

"""
    getindex(I::IntegralHelper,entry::String)

Called when `helper["foo"]` syntax is used. If the requested entry already
exists, simply return the entry. If not, compute the requested entry.
"""
function getindex(I::IntegralHelper, entry::String)
    if haskey(I.cache, entry)
        return I.cache[entry]
    else
        compute!(I, entry)
        return I.cache[entry]
    end
end

function setindex!(I::IntegralHelper, A::AbstractArray, key::String)
    I.cache[key] = A
end

function delete!(I::IntegralHelper, keys...)
    for k in keys
        delete!(I.cache, k)
    end
    GC.gc()
end

function delete!(I::IntegralHelper)
    delete!(I, keys(I.cache)...)
end

# Integrals specific to orbital types
include("AtomicIntegrals.jl")
include("ROIntegrals.jl")

# Pretty printing
function string_repr(X::IntegralHelper{T,E,O}) where {T,E,O}
    eri_string = replace("$E", "Fermi.Integrals."=>"")
    orb_string = replace("$O", "Fermi.Orbitals."=>"")
    out = ""
    out = out*" ⇒ Fermi IntegralHelper\n"
    out = out*" ⋅ Data Type:                 $(T)\n"
    out = out*" ⋅ Basis:                     $(X.basis)\n"
    out = out*" ⋅ ERI:                       $(eri_string)\n"
    out = out*" ⋅ Orbitals:                  $(orb_string)\n"
    cache_str = ""
    for k in keys(X.cache)
        cache_str *= k*" "
    end
    out = out*" ⋅ Stored Integrals:          $(cache_str)"
    return out
end

function string_repr(X::AbstractDFERI)
    eri_string = replace("$(typeof(X))", "Fermi.Integrals."=>"")
    return "$(eri_string): $(X.basisset.basis_name)"
end

function show(io::IO, ::MIME"text/plain", X::T) where T<:Union{IntegralHelper,JKFIT,RIFIT}
    print(io, string_repr(X))
end

end # Module


text:
module DLPack

using CUDA  # For exporting to CuArrays
using PyCall

export DLArray, DLMatrix, DLVector

@enum DLDeviceType::Cint begin
    kDLCPU = 1
    kDLGPU = 2
    kDLCPUPinned = 3
    kDLOpenCL = 4
    kDLVulkan = 7
    kDLMetal = 8
    kDLVPI = 9
    kDLROCM = 10
    kDLExtDev = 12
end

struct DLContext
    device_type::Cint
    device_id::Cint
end

@enum DLDataTypeCode::Cuint begin
    kDLInt = 0
    kDLUInt = 1
    kDLFloat = 2
    kDLBfloat = 4
end

struct DLDataType
    code::Cuchar
    bits::Cuchar
    lanes::Cushort
end

Base.convert(::Type{T}, code::DLDataTypeCode) where {T <: Integer} = T(code)

jltypes_to_dtypes() = Dict(
    Int8 => DLDataType(kDLInt, 8, 1),
    Int16 => DLDataType(kDLInt, 16, 1),
    Int32 => DLDataType(kDLInt, 32, 1),
    Int64 => DLDataType(kDLInt, 64, 1),
    UInt8 => DLDataType(kDLUInt, 8, 1),
    UInt16 => DLDataType(kDLUInt, 16, 1),
    UInt32 => DLDataType(kDLUInt, 32, 1),
    UInt64 => DLDataType(kDLUInt, 64, 1),
    Float16 => DLDataType(kDLFloat, 16, 1),
    Float32 => DLDataType(kDLFloat, 32, 1),
    Float64 => DLDataType(kDLFloat, 64, 1)
)

struct DLTensor
    data::Ptr{Cvoid}
    ctx::DLContext
    ndim::Cint
    dtype::DLDataType
    shape::Ptr{Clonglong}
    strides::Ptr{Clonglong}
    byte_offset::Culonglong
end

# Defined as mutable since we need a finalizer that calls `deleter`
# to destroy its original enclosing context `manager_ctx`
mutable struct DLManagedTensor
    dl_tensor::DLTensor
    manager_ctx::Ptr{Cvoid}
    deleter::Ptr{Cvoid}
end

function DLManagedTensor(po::PyObject)
    if !pyisinstance(po, PyCall.@pyglobalobj(:PyCapsule_Type))
        throw(ArgumentError("PyObject must be a PyCapsule"))
    end

    # Replace the capsule destructor to prevent it from deleting the tensor.
    # We will use the `DLManagedTensor.deleter` instead
    PyCall.@pycheck ccall(
        (@pysym :PyCapsule_SetDestructor),
        Cint, (PyPtr, Ptr{Cvoid}),
        po, C_NULL
    )

    dlptr = PyCall.@pycheck ccall(
        (@pysym :PyCapsule_GetPointer),
        Ptr{DLManagedTensor}, (PyPtr, Ptr{UInt8}),
        po, ccall((@pysym :PyCapsule_GetName), Ptr{UInt8}, (PyPtr,), po)
    )
    manager = unsafe_load(dlptr)

    if manager.deleter != C_NULL
        delete = manager -> ccall(manager.deleter, Cvoid, (Ptr{Cvoid},), Ref(manager))
        finalizer(delete, manager)
    end

    return manager
end

struct DLArray{T, N}
    manager::DLManagedTensor

    function DLArray{DLDataType, N}(po::PyObject) where {T, N}
        manager = DLManagedTensor(po)

        if N != (n = manager.dl_tensor.ndim)
            throw(ArgumentError("Dimensionality mismatch, object ndims is $n"))
        end

        return new(manager)
    end

    function DLArray{T, N}(po::PyObject) where {T, N}
        manager = DLManagedTensor(po)

        if N != (n = manager.dl_tensor.ndim)
            throw(ArgumentError("Dimensionality mismatch, object ndims is $n"))
        elseif jltypes_to_dtypes()[T] != (D = manager.dl_tensor.dtype)
            throw(ArgumentError("Type mismatch, object dtype is $D"))
        end

        return new(manager)
    end
end

const DLVector{T} = DLArray{T, 1}
const DLMatrix{T} = DLArray{T, 2}

###########
#  Utils  #
###########

device_type(ctx::DLContext) = DLDeviceType(ctx.device_type)
device_type(tensor::DLTensor) = device_type(tensor.ctx)
device_type(manager::DLManagedTensor) = device_type(manager.dl_tensor)
device_type(array::DLArray) = device_type(array.manager)

Base.eltype(array::DLArray{T}) where {T} = T

Base.ndims(array::DLArray{T, N}) where {T, N} = N

_size(tensor::DLTensor) = tensor.shape
_size(manager::DLManagedTensor) = _size(manager.dl_tensor)
_size(array::DLArray) = _size(array.manager)

function Base.size(array::DLArray{T, N}) where {T, N}
    ptr = Base.unsafe_convert(Ptr{NTuple{N, Int64}}, _size(array))
    return unsafe_load(ptr)
end
#
function Base.size(array::DLArray{T, N}, d::Integer) where {T, N}
    if 1 ≤ d
        return d ≤ N ? size(array)[d] : Int64(1)
    end
    throw(ArgumentError("Dimension out of range"))
end

_strides(tensor::DLTensor) = tensor.strides
_strides(manager::DLManagedTensor) = _strides(manager.dl_tensor)
_strides(array::DLArray) = _strides(array.manager)

function Base.strides(array::DLArray{T, N}) where {T, N}
    ptr = Base.unsafe_convert(Ptr{NTuple{N, Int64}}, _strides(array))
    return unsafe_load(ptr)
end

byte_offset(tensor::DLTensor) = Int(tensor.byte_offset)
byte_offset(manager::DLManagedTensor) = byte_offset(manager.dl_tensor)
byte_offset(array::DLArray) = byte_offset(array.manager)

Base.pointer(tensor::DLTensor) = tensor.data
Base.pointer(manager::DLManagedTensor) = pointer(manager.dl_tensor)
Base.pointer(array::DLArray) = pointer(array.manager)

##############
#  Wrappers  #
##############

function Base.unsafe_wrap(::Type{Array}, array::DLArray{T}) where {T}
    if device_type(array) == kDLCPU
        addr = Int(pointer(array))
        return GC.@preserve array unsafe_wrap(Array, Ptr{T}(addr), size(array))
    end
    throw(ArgumentError("Only CPU arrays can be wrapped with Array"))
end

function Base.unsafe_wrap(::Type{CuArray}, array::DLArray{T}) where {T}
    if device_type(array) == kDLGPU
        addr = Int(pointer(array))
        return GC.@preserve array unsafe_wrap(CuArray, CuPtr{T}(addr), size(array))
    end
    throw(ArgumentError("Only GPU arrays can be wrapped with CuArray"))
end

end


text:
# the main tester for FINUFFT.jl

using FINUFFT

using Test
using LinearAlgebra
using Random

# run a test at requested tolerance tol, floating-point precision dtype...
function test_nufft(tol::Real, dtype::DataType)
    @assert dtype <: FINUFFT.finufftReal

    Random.seed!(1)

    T = dtype    # abbrev; we no longer infer dtype as type of tol
                 # (this would be confusing since tol can be any type)
    nj = 10      # sizes for this small test: # NU pts
    nk = 11      # targ NU pts for t3
    ms = 12      # modes x
    mt = 13      # modes y
    mu = 14      # modes z

    # nonuniform data, using the full allowed input domain [-3pi,3pi)
    x = Array{T}(3*pi*(2*rand(nj).-1.0))
    y = Array{T}(3*pi*(2*rand(nj).-1.0))
    z = Array{T}(3*pi*(2*rand(nj).-1.0))
    c = rand(Complex{T},nj)
    s = rand(T,nk)
    t = rand(T,nk)
    u = rand(T,nk)
    f = rand(Complex{T},nk)

    # uniform data
    F1D = rand(Complex{T}, ms)
    F2D = rand(Complex{T}, ms, mt)
    F3D = rand(Complex{T}, ms, mt, mu)

    modevec(m) = -floor(m/2):floor((m-1)/2+1)
    k1 = modevec(ms)
    k2 = modevec(mt)
    k3 = modevec(mu)

    errfac = 100       # allowed multiple of tol for errors rel to direct calc
    errdifffac = 10    # allowed multiple of tol for errors rel to 2nd NUFFT

    @testset "NUFFT ($T)" begin
        ## 1D
        @testset "1D" begin
            # 1D1 - here we include tests of opts-setting, vectorized, guru...
            @testset "1D1" begin
                out = zeros(Complex{T},ms)
                ref = zeros(Complex{T},ms)     # direct calc...
                for j=1:nj
                    for ss=1:ms
                        ref[ss] += c[j] * exp(1im*k1[ss]*x[j])
                    end
                end
                # Simple, writing into array, setting some non-default opts...
                nufft1d1!(x, c, 1, tol, out, debug=1, spread_sort=0, dtype=T)
                relerr_1d1 = norm(vec(out)-vec(ref), Inf) / norm(vec(ref), Inf)
                @test relerr_1d1 < errfac*tol

                # Different caller which returns array
                out2 = nufft1d1(x, c, 1, tol, ms, dtype=T)
                reldiff = norm(vec(out)-vec(out2), Inf) / norm(vec(out), Inf)
                @test reldiff < errdifffac*tol

                # guru1d1
                plan = finufft_makeplan(1,[ms;],1,1,tol,dtype=T)
                finufft_setpts(plan,x)
                out3 = finufft_exec(plan,c)
                finufft_destroy(plan)
                relerr_guru = norm(vec(out3)-vec(ref), Inf) / norm(vec(ref), Inf)
                @test relerr_guru < errfac*tol

                # guru1d1 vectorized ("many")
                ntrans = 3         # let's stack 3 transforms at once
                plan = finufft_makeplan(1,[ms;],1,ntrans,tol,dtype=T)
                finufft_setpts(plan,x)
                cstack = hcat(c,2*c,3*c);           # change the coeff vectors
                out4 = finufft_exec(plan,cstack)
                finufft_destroy(plan)
                refstack = hcat(ref,2*ref,3*ref);   # ditto
                relerr_guru_many = norm(vec(out4)-vec(refstack), Inf) / norm(vec(refstack), Inf)
                @test relerr_guru_many < errfac*tol

                # simple vectorized ("many")
                fstack = nufft1d1(x,cstack,+1,tol,ms,dtype=T)
                relerr_many = norm(vec(fstack)-vec(refstack), Inf) / norm(vec(refstack), Inf)
                @test relerr_many < errfac*tol
            end

            # 1D2
            @testset "1D2" begin
                out = zeros(Complex{T},nj)
                ref = zeros(Complex{T},nj)
                for j=1:nj
                    for ss=1:ms
                        ref[j] += F1D[ss] * exp(1im*k1[ss]*x[j])
                    end
                end
                nufft1d2!(x, out, 1, tol, F1D, dtype=T)
                relerr_1d2 = norm(vec(out)-vec(ref), Inf) / norm(vec(ref), Inf)
                @test relerr_1d2 < errfac*tol
                out2 = nufft1d2(x, 1, tol, F1D, dtype=T)
                reldiff = norm(vec(out)-vec(out2), Inf) / norm(vec(out), Inf)
                @test reldiff < errdifffac*tol
            end

            # 1D3
            @testset "1D3" begin
                out = zeros(Complex{T},nk)
                ref = zeros(Complex{T},nk)
                for k=1:nk
                    for j=1:nj
                        ref[k] += c[j] * exp(1im*s[k]*x[j])
                    end
                end
                nufft1d3!(x,c,1,tol,s,out, dtype=T)
                relerr_1d3 = norm(vec(out)-vec(ref), Inf) / norm(vec(ref), Inf)
                @test relerr_1d3 < errfac*tol
                out2 = nufft1d3(x,c,1,tol,s, dtype=T)
                reldiff = norm(vec(out)-vec(out2), Inf) / norm(vec(out), Inf)
                @test reldiff < errdifffac*tol
            end
        end

        ## 2D
        @testset "2D" begin
            @testset "2D1" begin
                # 2D1
                out = zeros(Complex{T},ms,mt)
                ref = zeros(Complex{T},ms,mt)
                for j=1:nj
                    for ss=1:ms
                        for tt=1:mt
                            ref[ss,tt] += c[j] * exp(1im*(k1[ss]*x[j]+k2[tt]*y[j]))
                        end
                    end
                end
                nufft2d1!(x, y, c, 1, tol, out, dtype=T)
                relerr_2d1 = norm(vec(out)-vec(ref), Inf) / norm(vec(ref), Inf)
                @test relerr_2d1 < errfac*tol
                out2 = nufft2d1(x, y, c, 1, tol, ms, mt, dtype=T)
                reldiff = norm(vec(out)-vec(out2), Inf) / norm(vec(out), Inf)
                @test reldiff < errdifffac*tol
            end

            @testset "2D2" begin
                # 2D2
                out = zeros(Complex{T},nj)
                ref = zeros(Complex{T},nj)
                for j=1:nj
                    for ss=1:ms
                        for tt=1:mt
                            ref[j] += F2D[ss, tt] * exp(1im*(k1[ss]*x[j]+k2[tt]*y[j]))
                        end
                    end
                end
                nufft2d2!(x, y, out, 1, tol, F2D, dtype=T)
                relerr_2d2 = norm(vec(out)-vec(ref), Inf) / norm(vec(ref), Inf)
                @test relerr_2d2 < errfac*tol
                out2 = nufft2d2(x, y, 1, tol, F2D, dtype=T)
                reldiff = norm(vec(out)-vec(out2), Inf) / norm(vec(out), Inf)
                @test reldiff < errdifffac*tol
            end

            @testset "3D3" begin
                # 2D3
                out = zeros(Complex{T},nk)
                ref = zeros(Complex{T},nk)
                for k=1:nk
                    for j=1:nj
                        ref[k] += c[j] * exp(1im*(s[k]*x[j]+t[k]*y[j]))
                    end
                end
                nufft2d3!(x,y,c,1,tol,s,t,out, dtype=T)
                relerr_2d3 = norm(vec(out)-vec(ref), Inf) / norm(vec(ref), Inf)
                @test relerr_2d3 < errfac*tol
                out2 = nufft2d3(x,y,c,1,tol,s,t, dtype=T)
                reldiff = norm(vec(out)-vec(out2), Inf) / norm(vec(out), Inf)
                @test reldiff < errdifffac*tol
            end
        end

        ## 3D
        @testset "3D" begin
            @testset "3D1" begin
                # 3D1
                out = zeros(Complex{T},ms,mt,mu)
                ref = zeros(Complex{T},ms,mt,mu)
                for j=1:nj
                    for ss=1:ms
                        for tt=1:mt
                            for uu=1:mu
                                ref[ss,tt,uu] += c[j] * exp(1im*(k1[ss]*x[j]+k2[tt]*y[j]+k3[uu]*z[j]))
                            end
                        end
                    end
                end
                nufft3d1!(x, y, z, c, 1, tol, out, dtype=T)
                relerr_3d1 = norm(vec(out)-vec(ref), Inf) / norm(vec(ref), Inf)
                @test relerr_3d1 < errfac*tol
                out2 = nufft3d1(x, y, z, c, 1, tol, ms, mt, mu, dtype=T)
                reldiff = norm(vec(out)-vec(out2), Inf) / norm(vec(out), Inf)
                @test reldiff < errdifffac*tol
            end

            @testset "3D2" begin
                # 3D2
                out = zeros(Complex{T},nj)
                ref = zeros(Complex{T},nj)
                for j=1:nj
                    for ss=1:ms
                        for tt=1:mt
                            for uu=1:mu
                                ref[j] += F3D[ss, tt, uu] * exp(1im*(k1[ss]*x[j]+k2[tt]*y[j]+k3[uu]*z[j]))
                            end
                        end
                    end
                end
                nufft3d2!(x, y, z, out, 1, tol, F3D, dtype=T)
                relerr_3d2 = norm(vec(out)-vec(ref), Inf) / norm(vec(ref), Inf)
                @test relerr_3d2 < errfac*tol
                out2 = nufft3d2(x, y, z, 1, tol, F3D, dtype=T)
                reldiff = norm(vec(out)-vec(out2), Inf) / norm(vec(out), Inf)
                @test reldiff < errdifffac*tol
            end

            @testset "3D3" begin
                # 3D3
                out = zeros(Complex{T},nk)
                ref = zeros(Complex{T},nk)
                for k=1:nk
                    for j=1:nj
                        ref[k] += c[j] * exp(1im*(s[k]*x[j]+t[k]*y[j]+u[k]*z[j]))
                    end
                end
                nufft3d3!(x,y,z,c,1,tol,s,t,u,out, dtype=T)
                relerr_3d3 = norm(vec(out)-vec(ref), Inf) / norm(vec(ref), Inf)
                @test relerr_3d3 < errfac*tol
                out2 = nufft3d3(x,y,z,c,1,tol,s,t,u, dtype=T)
                reldiff = norm(vec(out)-vec(out2), Inf) / norm(vec(out), Inf)
                @test reldiff < errdifffac*tol
            end
        end
    end
end

# Main: do the tests
# 1st arg is tolerance (no longer used to infer dtype), 2nd is dtype...
test_nufft(1e-14, Float64)
test_nufft(1e-4, Float32)


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

"""
    innermap(f::Base.Callable, A::AbstractArray{<:AbstractArray})

Nested `map` at depth 2. Equivalent to `map(X -> map(f, X) A)`.
"""
function innermap end
export innermap

innermap(f::Base.Callable, A::AbstractArray{<:AbstractArray{T,M},N}) where {T,M,N} =
    map(X -> map(f, X), A)

"""
    deepmap(f::Base.Callable, x::Any)
    deepmap(f::Base.Callable, A::AbstractArray{<:AbstractArray{<:...}})

Applies `map` at the deepest possible layer of nested arrays. If `A` is not
a nested array, `deepmap` behaves identical to `Base.map`.
"""
function deepmap end
export deepmap

deepmap(f::Base.Callable, x::Any) = map(f, x)

deepmap(f::Base.Callable, A::AbstractArray{<:AbstractArray}) =
    map(X -> deepmap(f, X), A)

"""
    flatview(A::AbstractArray)
    flatview(A::AbstractArray{<:AbstractArray{<:...}})

View array `A` in a suitable flattened form. The shape of the flattened form
will depend on the type of `A`. If the `A` is not a nested array, the return
value is `A` itself. When no type-specific method is available, `flatview`
will fall back to `Base.Iterators.flatten`.
"""
function flatview end
export flatview

@inline flatview(A::AbstractArray) = A

# TODO: Implement flatview on generic nested arrays via new `FlatView`, using
# deepgetindex to implement getindex, etc.
@inline flatview(A::AbstractArray{<:AbstractArray}) = Base.Iterators.flatten(A)

"""
    nestedview(A::AbstractArray{T,M+N}, M::Integer)
    nestedview(A::AbstractArray{T,2})

AbstractArray{<:AbstractArray{T,M},N}

View array `A` in as an `M`-dimensional array of `N`-dimensional arrays by
wrapping it into an [`ArrayOfSimilarArrays`](@ref).

It's also possible to use a `StaticVector` of length `S` as the type of the
inner arrays via

    nestedview(A::AbstractArray{T}, ::Type{StaticArrays.SVector{S}})
    nestedview(A::AbstractArray{T}, ::Type{StaticArrays.SVector{S,T}})
"""
function nestedview end
export nestedview

@inline nestedview(A::AbstractArray{T,L}, M::Integer) where {T,L} =
    ArrayOfSimilarArrays{T,M}(A)

@inline nestedview(A::AbstractArray{T,L}, ::Val{M}) where {T,L,M} =
    ArrayOfSimilarArrays{T,M}(A)

@inline nestedview(A::AbstractArray{T,2}) where {T} =
    VectorOfSimilarVectors{T}(A)

"""
    innersize(A:AbstractArray{<:AbstractArray}, [dim])

Returns the size of the element arrays of `A`. Fails if the element arrays
are not of equal size.
"""
function innersize end
export innersize

function innersize(A::AbstractArray{<:AbstractArray{T,M},N}) where {T,M,N}
    s = if !isempty(A)
        let sz_A = size(first(A))
            ntuple(i -> Int(sz_A[i]), Val(M))
        end
    else
        ntuple(_ -> zero(Int), Val(M))
    end

    let s = s
        if any(X -> size(X) != s, A)
            throw(DimensionMismatch("Shape of element arrays of A is not equal, can't determine common shape"))
        end
    end

    s
end

@inline innersize(A::AbstractArray{<:AbstractArray}, dim::Integer) =
    innersize(A)[dim]

"""
    deepgetindex(A::AbstractArray, idxs...)
    deepgetindex(A::AbstractArray{<:AbstractArray, N}, idxs...) where {N}

Recursive `getindex` on flat or nested arrays. If A is an array of arrays,
uses the first `N` entries in `idxs` on `A`, then the rest on the inner
array(s). If A is not a nested array, `deepgetindex` is equivalent to
`getindex`.

See also [`deepsetindex!`](@ref) and [`deepview`](@ref).
"""
function deepgetindex end
export deepgetindex

Base.@propagate_inbounds deepgetindex(A::AbstractArray{T,N}, idxs::Vararg{Any,N}) where {T,N} = getindex(A, idxs...)
Base.@propagate_inbounds deepgetindex(A::AbstractArray{<:AbstractArray,N}, idxs::Vararg{Any,N}) where {N} = getindex(A, idxs...)

Base.@propagate_inbounds function deepgetindex(A::AbstractArray{<:AbstractArray,N}, idxs...) where {N}
    idxs_outer, idxs_inner = split_tuple(idxs, Val{N}())
    _deepgetindex_impl(A, idxs_outer, idxs_inner)
end

Base.@propagate_inbounds _deepgetindex_impl(A::AbstractArray{<:AbstractArray}, idxs_outer::NTuple{N,Real}, idxs_inner::Tuple) where {N} =
    deepgetindex(getindex(A, idxs_outer...), idxs_inner...)

Base.@propagate_inbounds _deepgetindex_impl(A::AbstractArray{<:AbstractArray}, idxs_outer::NTuple{N,Any}, idxs_inner::Tuple) where {N} =
    _deepgetindex_tupled.(view(A, idxs_outer...), (idxs_inner,))

Base.@propagate_inbounds _deepgetindex_tupled(A::AbstractArray, idxs::Tuple) = deepgetindex(A, idxs...)

"""
    deepsetindex!(A::AbstractArray, x, idxs...)
    deepsetindex!(A::AbstractArray{<:AbstractArray,N}, x, idxs...) where {N}

Recursive `setindex!` on flat or nested arrays. If A is an array of arrays,
uses the first `N` entries in `idxs` on `A`, then the rest on the inner
array(s). If A is not a nested array, `deepsetindex!` is equivalent to
`setindex!`.

See also [`deepgetindex`](@ref) and [`deepview`](@ref).
"""
function deepsetindex! end
export deepsetindex!

Base.@propagate_inbounds deepsetindex!(A::AbstractArray{T,N}, x, idxs::Vararg{Any,N}) where {T,N} = setindex!(A, x, idxs...)
Base.@propagate_inbounds deepsetindex!(A::AbstractArray{<:AbstractArray,N}, x, idxs::Vararg{Any,N}) where {N} = setindex!(A, x, idxs...)

Base.@propagate_inbounds function deepsetindex!(A::AbstractArray{<:AbstractArray,N}, x, idxs...) where {N}
    idxs_outer, idxs_inner = split_tuple(idxs, Val{N}())
    _deepsetindex_impl!(A, x, idxs_outer, idxs_inner)
    A
end

Base.@propagate_inbounds function _deepsetindex_impl!(A::AbstractArray{<:AbstractArray}, x, idxs_outer::NTuple{N,Real}, idxs_inner::Tuple) where {N}
    B = getindex(A, idxs_outer...)
    deepsetindex!(B, x, idxs_inner...)
end

Base.@propagate_inbounds function _deepsetindex_impl!(A::AbstractArray{<:AbstractArray}, x, idxs_outer::NTuple{N,Any}, idxs_inner::Tuple) where {N}
    B = view(A, idxs_outer...)
    for i in eachindex(B, x)
        deepsetindex!(B[i], x[i], idxs_inner...)
    end
end

"""
    deepview(A::AbstractArray, idxs...)
    deepview(A::AbstractArray{<:AbstractArray, N}, idxs...) where {N}

Recursive `view` on flat or nested arrays. If A is an array of arrays,
uses the first `N` entries in `idxs` on `A`, then the rest on the inner
array(s). If A is not a nested array, `deepview` is equivalent to `view`.

See also [`deepgetindex`](@ref) and [`deepsetindex!`](@ref).
"""
function deepview end
export deepview

Base.@propagate_inbounds deepview(A::AbstractArray{T,N}, idxs::Vararg{Any,N}) where {T,N} = view(A, idxs...)
Base.@propagate_inbounds deepview(A::AbstractArray{<:AbstractArray,N}, idxs::Vararg{Any,N}) where {N} = view(A, idxs...)

Base.@propagate_inbounds function deepview(A::AbstractArray{<:AbstractArray,N}, idxs...) where {N}
    idxs_outer, idxs_inner = split_tuple(idxs, Val{N}())
    _deepview_impl(A, idxs_outer, idxs_inner)
end

Base.@propagate_inbounds _deepview_impl(A::AbstractArray{<:AbstractArray}, idxs_outer::NTuple{N,Real}, idxs_inner::NTuple{M,Real}) where {N,M} =
    deepview(getindex(A, idxs_outer...), idxs_inner...)

Base.@propagate_inbounds _deepview_impl(A::AbstractArray{<:AbstractArray}, idxs_outer::NTuple{N,Real}, idxs_inner::NTuple{M,Any}) where {N,M} =
    deepview(getindex(A, idxs_outer...), idxs_inner...)

Base.@propagate_inbounds _deepview_impl(A::AbstractArray{<:AbstractArray}, idxs_outer::NTuple{N,Any}, idxs_inner::NTuple{M,Real}) where {N,M} =
    throw(ArgumentError("deepview not supported yes with outer indices $idxs_outer and inner indices $idxs_inner"))

Base.@propagate_inbounds _deepview_impl(A::AbstractArray{<:AbstractArray}, idxs_outer::NTuple{N,Any}, idxs_inner::NTuple{M,Any}) where {N,M} =
    _deepview_tupled.(view(A, idxs_outer...), (idxs_inner,))

Base.@propagate_inbounds _deepview_tupled(A::AbstractArray, idxs::Tuple) = deepview(A, idxs...)

"""
    abstract_nestedarray_type(T_inner::Type, ::Val{ndims_tuple})

Return the type of nested `AbstractArray`s. `T_inner` specifies the element
type of the innermost layer of arrays, `ndims_tuple` specifies the
dimensionality of each nesting layer (outer arrays first).

If `ndims_tuple` is empty, the returns is the (typically scalar) type
`T_inner` itself.
"""
function abstract_nestedarray_type end
export abstract_nestedarray_type

Base.@pure function abstract_nestedarray_type(::Type{T_inner}, outer::Val{ndims_tuple}) where {T_inner,ndims_tuple}
    _abstract_nestedarray_type_impl(T_inner, ndims_tuple...)
end

Base.@pure _abstract_nestedarray_type_impl(::Type{T_inner}) where {T_inner} = T_inner

Base.@pure _abstract_nestedarray_type_impl(::Type{T_inner}, N) where {T_inner} = AbstractArray{T_inner, N}

Base.@pure _abstract_nestedarray_type_impl(::Type{T_inner}, N, M, ndims_tuple...) where {T_inner} =
    AbstractArray{<:_abstract_nestedarray_type_impl(T_inner, M, ndims_tuple...), N}


text:
@with_kw struct FreeBody{R,T} <: RigidBody{R}
    mass::T = 1.0
    J::Diagonal{T,SVector{3,T}} = Diagonal(@SVector ones(3))
    Jinv::Diagonal{T,SVector{3,T}} = Diagonal(@SVector ones(3))
end

(::FreeBody)(;kwargs...) = FreeBody{UnitQuaternion{Float64,CayleyMap},Float64}(;kwargs...)

function wrenches(model::FreeBody, x::SVector, u::SVector)
    F = forces(model, x, u)
    M = moments(model, x, u)
    return F,M
end
function forces(model::FreeBody, x::SVector, u::SVector)
    q = orientation(model, x)
    F = @SVector [u[1], u[2], u[3]]
    q*F  # world frame
end

function moments(model::FreeBody, x::SVector, u::SVector)
    return @SVector [u[4], u[5], u[6]]  # body frame
end

Base.size(::FreeBody{<:UnitQuaternion}) = 13,6
Base.size(::FreeBody{<:Rotation}) = 12,6

inertia(model::FreeBody, x, u) = model.J
inertia_inv(model::FreeBody, x, u) = model.Jinv
mass_matrix(model::FreeBody, x, u) = Diagonal(@SVector fill(model.mass,3))


text:
@graphs begin
	@node phase_4 = Start() -> p4_roll
	@node p4_roll = GetAvailableDice("""
									 Roll all action dice not in the hunt box. Place all Eye results in the hunt box and input the remaining dice.
									 """) -> p4_end
	@node p4_end = End("End of Phase") -> []

	@node adjust_dice = Start() -> adjust_roll
	@node adjust_roll = GetAvailableDice() -> adjust_end
	@node adjust_end = End() -> []
end

max_stars_repo_name: mvmorin/queller-bot

max_issues_repo_name: mvmorin/queller-bot

max_forks_repo_name: mvmorin/queller-bot


text:
# -------------------- motifs ---------------------------------
# concentration parameter for Dirichlet that generate Categorical Distributions

const motif_alpha_pc = 0.1;
const DNA_char_cap = ['A','C','G','T'];
const DNA_char_nocap = ['a','c','g','t'];
const mixture_weights_alpha = 50;               # α for generating mixture weight using Dirichlet
                                                # assume for now mixture weights are roughly the same
                                                # Note that modes are vector of unit ranges
const percent_complement = .5;                  # probability that a motif is represented in its reverse complement direction
const DNA_complement = Dict('A'=>'T','C'=>'G','G'=>'C','T'=>'A');
const data_pt_len = 100;                        # Length of each data string
const background = Categorical(0.25*ones(4));   # simple i.i.d flat background for now
const bg_array = [.25, .25, .25, .25];
const sample_near_center = 5;                   #= increment and decrement the avail. position
                                                    from the start and the end, resp. =#

const v_f64 = Vector{Float64};
const cat = Categorical{Float64, v_f64};
const u_conc = Union{Float64, v_f64, Nothing};

const max_num_read_fasta = 3000;
const nmbr_regex = r"\-?[0-9]+\.?[0-9]*";

# dummy reference
const dummy = Dict('A'=>Array{Float32}([1, 0, 0, 0]),
                   'C'=>Array{Float32}([0, 1, 0, 0]),
                   'G'=>Array{Float32}([0, 0, 1, 0]),
                   'T'=>Array{Float32}([0, 0, 0, 1]));

# product categorical

struct product_categorical
    pc::Vector{cat}      # product categorical
    len::Int64              # length of the product
    dim::Int64              # dimension for each categorical in the product
    alpha::u_conc           # concentration parameter;

    function product_categorical(L::Int,
                                 dim::Int,
                                 alpha::Union{Real, v_f64})
        if typeof(alpha) == v_f64
            length(alpha) != dim && error("alpha must be of length same as the dimension");
        end
        D = typeof(alpha) == v_f64 ? Dirichlet(alpha) : Dirichlet(dim, alpha);
        new([Categorical(rand(D)) for _ in 1:L], L, dim, alpha)
    end

    function product_categorical(mat::Matrix{T}) where T <: Real
        L, width = size(mat);
        @assert width == 4 "the input matrix must have width 4"
        new([Categorical(mat[i,:]) for i = 1:L], L, 4, nothing)
    end
end

view(v::product_categorical) = [round(v.pc[j].p[i], digits=2) for j = 1:v.len, i = 1:v.dim];
Base.rand(v::product_categorical) = [rand(v.pc[i]) for i = 1:v.len];

# motifs definitions

struct single_part_motif
    P::product_categorical # P stands for probability measure
    len::Int64
    function single_part_motif(L::T) where T <: Integer
        new(product_categorical(L, 4, motif_alpha_pc), L)
    end
    function single_part_motif(mat::Matrix{T}) where T <: Real
        new(product_categorical(mat), size(mat,1));
    end
end

# example ----------------------------------
# a = single_part_motif(5);
# view(a)
# ------------------------------------------

# an intermediate type for gap motifs
struct k_parts_motif
    P::Vector{single_part_motif}
    K::Int64
    """
    K: number of parts in the motif
    L_vec: vector of lengths L's for each product categorical
    """
    function k_parts_motif(length_vec::Vector{Int})
        K = length(length_vec);
        new([single_part_motif(length_vec[i]) for i = 1:K], K)
    end
end

#=
motif with k parts and gaps in between every adjacent parts
Input:
    L_vec:      Length of each of the parts of the motif
    gap_vec:    Length of each of the gap in between parts --
                gap_vecᵢ specifies the gap in between ith and
                (i+1)th part
    gap_dist:   Multinomial that generates the alphabets in an i.i.d fashion
=#

struct gapped_k_parts_motif
    P_motifs::k_parts_motif
    gap_len::Vector{Int}    # the maximal gap length that can occur in between each parts
    P_gap::Categorical{Float64, Vector{Float64}}
    K::Int
    function gapped_k_parts_motif(L_vec::Vector{T}, gap_vec::Vector{T}, gap_dist=bg_array) where T <: Integer
        K = length(L_vec);
        @assert length(gap_vec) == K-1 "Length of the gap_vec (2nd argument) must be length(L_vec)-1"
        new(k_parts_motif(L_vec), gap_vec, Categorical(gap_dist), K)
    end
end

struct mixture_gapped_k_parts_motifs
    modes::Vector{UnitRange{Int}} # inclusive
    mixture_weights::cat
    motif::gapped_k_parts_motif
    num_modes::Int
    function mixture_gapped_k_parts_motifs(
        L_vec::Vector{Int},
        gap_vec::Vector{Int},
        modes::Vector{UnitRange{Int64}},
        gap_dist=bg_array
    )
        num_modes = length(modes);
        new(modes,
            Categorical(rand(Dirichlet(num_modes, mixture_weights_alpha))),
            gapped_k_parts_motif(L_vec, gap_vec, gap_dist),
            num_modes
        )
    end
end

# a simpler mixture motif than the one above

struct mixture_k_parts_motifs
    modes::Vector{UnitRange{Int}} # inclusive
    mixture_weights::cat
    motif::single_part_motif
    num_modes::Int
    function mixture_k_parts_motifs(modes::Vector{UnitRange{Int64}})
        K = max_v(modes);
        num_modes = length(modes);
        new(modes,
            Categorical(rand(Dirichlet(num_modes, mixture_weights_alpha))),
            single_part_motif(K),
            num_modes
        )
    end
end

# motif type
const motif_type = Union{single_part_motif,
                         gapped_k_parts_motif,
                         mixture_k_parts_motifs,
                         mixture_gapped_k_parts_motifs,
                         Nothing
};

# type for a (simulated) single dna background string with a motif in it
const sim_dna_str_w_motif = NamedTuple{(:str, :motif_where, :mode, :complement),
                    Tuple{String, UnitRange{Int64}, Int64, Bool}
                    };
const public_dna_str_w_motif = NamedTuple{(:str, :motif_where, :mode),
                    Tuple{String, UnitRange{Int64}, Int64}
                    };


text:
using Weave
using Test

function convert_test(outfile, infile="documents/chunk_options.noweb")
  outfile = joinpath("documents/convert", outfile)
  convert_doc(infile, outfile)
  result =  read(outfile, String)
  ref =  read(outfile * ".ref", String)
  @test result == ref
  rm(outfile)
end

convert_test("chunk_options.jmd")
convert_test("chunk_options.jl")
convert_test("chunk_options.mdw")
convert_test("chunk_options_nb.mdw", "documents/chunk_options.ipynb")

# Separate test for notebook (output depends on julia version)
function contents(chunk::Weave.DocChunk)
  return join([strip(c.content) for c in chunk.content], "")
end

function contents(chunk::Weave.CodeChunk)
  return chunk.content
end

function contents(doc::Weave.WeaveDoc)
  return join([contents(chunk) for chunk in doc.chunks], "")
end

outfile = "documents/convert/chunk_options.ipynb"
infile = "documents/chunk_options.noweb"
convert_doc(infile, outfile)
input = contents(Weave.read_doc(infile))
output = contents(Weave.read_doc(outfile))
@test input == output
rm(outfile)

# Test script reader
@test contents(
    Weave.read_doc("documents/chunk_options.noweb")) == contents(
      Weave.read_doc("documents/chunk_options.jl"))


text:
struct Direction{T} end

const Left = Direction{:Left}()
const Right = Direction{:Right}()

const Beginning = Left
const Ending = Right

struct Endpoint{T, D, B <: Bound}
    endpoint::T

    function Endpoint{T,D,B}(ep::T) where {T, D, B <: Bounded}
        @assert D isa Direction
        new{T,D,B}(ep)
    end

    function Endpoint{T,D,B}(ep::Nothing) where {T, D, B <: Unbounded}
        @assert D isa Direction
        new{T,D,B}()
    end

    function Endpoint{T,D,B}(ep::Nothing) where {T, D, B <: Bounded}
        throw(MethodError(Endpoint{T,D,B}, (ep,)))
    end
end

Endpoint{T,D,B}(ep) where {T, D, B <: Bounded} = Endpoint{T,D,B}(convert(T, ep))

const LeftEndpoint{T,B} = Endpoint{T, Left, B} where {T,B}
const RightEndpoint{T,B} = Endpoint{T, Right, B} where {T,B}

LeftEndpoint{B}(ep::T) where {T,B} = LeftEndpoint{T,B}(ep)
RightEndpoint{B}(ep::T) where {T,B} = RightEndpoint{T,B}(ep)

LeftEndpoint(i::AbstractInterval{T,L,R}) where {T,L,R} = LeftEndpoint{T,L}(L !== Unbounded ? first(i) : nothing)
RightEndpoint(i::AbstractInterval{T,L,R}) where {T,L,R} = RightEndpoint{T,R}(R !== Unbounded ? last(i) : nothing)

endpoint(x::Endpoint) = isbounded(x) ? x.endpoint : nothing
bound_type(x::Endpoint{T,D,B}) where {T,D,B} = B

isclosed(x::Endpoint) = bound_type(x) === Closed
isunbounded(x::Endpoint) = bound_type(x) === Unbounded
isbounded(x::Endpoint) = bound_type(x) !== Unbounded

function Base.hash(x::Endpoint{T,D,B}, h::UInt) where {T,D,B}
    h = hash(:Endpoint, h)
    h = hash(D, h)
    h = hash(B, h)

    # Bounded endpoints can skip hashing `T` as this is handled by hashing the endpoint
    # value field. Unbounded endpoints however, should ignore the stored garbage value and
    # should hash `T`.`
    h = B !== Unbounded ? hash(x.endpoint, h) : hash(T, h)

    return h
end

Base.broadcastable(e::Endpoint) = Ref(e)

"""
    ==(a::Endpoint, b::Endpoint) -> Bool

Determine if two endpoints are equal. When both endpoints are left or right then the points
and inclusiveness must be the same.

Checking the equality of left-endpoint and a right-endpoint is slightly more difficult. A
left-endpoint and a right-endpoint are only equal when they use the same point and are
both included. Note that left/right endpoints which are both not included are not equal
as the left-endpoint contains values below that point while the right-endpoint only contains
values that are above that point.

Visualizing two contiguous intervals can assist in understanding this logic:

    [x..y][y..z] -> RightEndpoint == LeftEndpoint
    [x..y)[y..z] -> RightEndpoint != LeftEndpoint
    [x..y](y..z] -> RightEndpoint != LeftEndpoint
    [x..y)(y..z] -> RightEndpoint != LeftEndpoint
"""
function Base.:(==)(a::Endpoint, b::Endpoint)
    return (
        isunbounded(a) && isunbounded(b) ||
        a.endpoint == b.endpoint && bound_type(a) == bound_type(b)
    )
end

function Base.:(==)(a::LeftEndpoint, b::RightEndpoint)
    a.endpoint == b.endpoint && isclosed(a) && isclosed(b)
end

function Base.:(==)(a::RightEndpoint, b::LeftEndpoint)
    b == a
end

function Base.isequal(a::Endpoint, b::Endpoint)
    return (
        isunbounded(a) && isunbounded(b) ||
        isequal(a.endpoint, b.endpoint) && isequal(bound_type(a), bound_type(b))
    )
end

function Base.isequal(a::LeftEndpoint, b::RightEndpoint)
    isequal(a.endpoint, b.endpoint) && isclosed(a) && isclosed(b)
end

function Base.isequal(a::RightEndpoint, b::LeftEndpoint)
    isequal(b, a)
end

function Base.isless(a::LeftEndpoint, b::LeftEndpoint)
    return (
        !isunbounded(b) && (
            isunbounded(a) ||
            a.endpoint < b.endpoint ||
            a.endpoint == b.endpoint && isclosed(a) && !isclosed(b)
        )
    )
end

function Base.isless(a::RightEndpoint, b::RightEndpoint)
    return (
        !isunbounded(a) && (
            isunbounded(b) ||
            a.endpoint < b.endpoint ||
            a.endpoint == b.endpoint && !isclosed(a) && isclosed(b)
        )
    )
end

function Base.isless(a::LeftEndpoint, b::RightEndpoint)
    return (
        isunbounded(a) ||
        isunbounded(b) ||
        a.endpoint < b.endpoint
    )
end

function Base.isless(a::RightEndpoint, b::LeftEndpoint)
    return (
        !isunbounded(a) && !isunbounded(b) &&
        (
            a.endpoint < b.endpoint ||
            a.endpoint == b.endpoint && !(isclosed(a) && isclosed(b))
        )
    )
end

# Comparisons between Scalars and Endpoints
Base.:(==)(a, b::Endpoint) = a == b.endpoint && isclosed(b)
Base.:(==)(a::Endpoint, b) = b == a

function Base.isless(a, b::LeftEndpoint)
    return (
        !isunbounded(b) && (
            a < b.endpoint ||
            a == b.endpoint && !isclosed(b)
        )
    )
end

function Base.isless(a::RightEndpoint, b)
    return (
        !isunbounded(a) &&
        (
            a.endpoint < b ||
            a.endpoint == b && !isclosed(a)
        )
    )
end

Base.isless(a, b::RightEndpoint) = isunbounded(b) || a < b.endpoint
Base.isless(a::LeftEndpoint, b)  = isunbounded(a) || a.endpoint < b


text:
##################
# Define Physics #
##################

@present Flow2DQuantities <: ExtendedFreeExtCalc2D begin
  X::Space

  k::Hom(Form1(X), Form1(X))    # diffusivity (usually scalar multiplication)
  kᵨ::Hom(Form1(X), Form1(X))   # REMOVE: Diffusion of density
  R₀::Hom(Form0(X), Form0(X))    # Ideal gas constant (usually scalar multiplication)
  kᵥ::Hom(Form1(X), Form1(X))    # viscosity (usually scalar multiplication)
  e2t::Hom(Form0(X), Form0(X))    # scale between internal energy and temperature (ideal gas)
  t2e::Hom(Form0(X), Form0(X))    # inverse of e2t (ideal gas)
  e2p::Hom(Form0(X), Form0(X))    # # scale between internal energy × density and pressure (ideal gas)
  L₀::Hom(Form1(X)⊗DualForm2(X), DualForm2(X))
  L₁::Hom(Form1(X)⊗DualForm1(X), DualForm1(X))
  L₁′::Hom(Form1(X)⊗Form1(X), Form1(X))
  ∂₀::Hom(Form0(X), Form0(X)) # all_wall boundary condition
  ∂₀₊::Hom(Form0(X), Form0(X)) # top/bottom wall boundary condition
  ∂₀ₛ::Hom(Form0(X), Form0(X)) # Sticking boundary condition
  ∂₁ₛ::Hom(Form1(X), Form1(X)) # Sticking boundary condition
  ∂₁ₗ₊::Hom(Form1(X), Form1(X)) # Sticking boundary condition
  ∂₁ₑ::Hom(Form1(X), Form1(X)) # In/Out edge flow boundary condition
  ∂ₜₐ::Hom(Form0(X), Form0(X))  # Temperature advection boundary condition
  ∂ₜ::Hom(Form0(X), Form0(X))  # internal energy boundary condition
  ∂ₚ::Hom(Form0(X), Form0(X))  # Pressure boundary condition
  ∂ₚₐ::Hom(Form0(X), Form0(X))  # Pressure advection boundary condition
  ∂ᵥ::Hom(Form1(X), Form1(X))  # Velocity boundary condition
  dneg₁::Hom(DualForm1(X), DualForm1(X)) # In/Out edge flow boundary condition
  neg₁::Hom(Form1(X), Form1(X)) # In/Out edge flow boundary condition
  dneg₀::Hom(DualForm0(X), DualForm0(X)) # In/Out edge flow boundary condition
  neg₀::Hom(Form0(X), Form0(X)) # In/Out edge flow boundary condition
  mask₁ₑ::Hom(Form1(X), Form1(X)) # In/Out edge flow boundary condition
  mask₀ₗ::Hom(Form0(X), Form0(X)) # In/Out edge flow boundary condition
  mask₀ₛ::Hom(Form0(X), Form0(X)) # In/Out edge flow boundary condition
  half::Hom(Form0(X), Form0(X)) # half
  third::Hom(Form1(X), Form1(X))    # multiplication by 1/3
  i₁′::Hom(Form1(X)⊗Form1(X), Form0(X))
  mult₀::Hom(Form0(X)⊗Form0(X), Form0(X)) # elementwise multiplication
  div₀::Hom(Form0(X)⊗Form0(X), Form0(X)) # elementwise division
  div₁::Hom(Form1(X)⊗Form1(X), Form1(X)) # elementwise division
  avg₀₁::Hom(Form0(X), Form1(X))
  sub₀::Hom(Form0(X)⊗Form0(X), Form0(X))
end

Diffusion = @decapode Flow2DQuantities begin
  (T, Ṫ)::Form0{X}
  ϕ::DualForm1{X}

  # Fick's first law
  ϕ ==  ⋆₁{X}(k(d₀{X}(T)))
  # Diffusion equation
  Ṫ == ⋆₀⁻¹{X}(dual_d₁{X}(ϕ))
end

Advection = @decapode Flow2DQuantities begin
  (T, Ṫ)::Form0{X}
  V::Form1{X}
  Ṫ == neg₀(⋆₀⁻¹{X}(L₀(V, ⋆₀{X}(T))))
end

Superposition = @decapode Flow2DQuantities begin
  (Ṫ₁, Ṫ₂, Ṫ, T)::Form0{X}
  Ṫ == Ṫ₁ + Ṫ₂
  ∂ₜ{Form0{X}}(T) == Ṫ
end

compose_diff_adv = @relation (T, Ṫ, V) begin
  diffusion(T, Ṫ₁)
  advection(T, Ṫ₂, V)
  superposition(Ṫ₁, Ṫ₂, Ṫ, T)
end
#=
DiffusionAdvection = oapply(compose_diff_adv,
                  [OpenDiagram(Diffusion, [:T, :Ṫ]),
                   OpenDiagram(Advection, [:T, :Ṫ, :V]),
                   OpenDiagram(Superposition, [:Ṫ₁, :Ṫ₂, :Ṫ, :T])])

dwd_res = diag2dwd(diff_adv.functor, in_vars = [:V, :T])
draw_equations(diff_adv.functor)
draw_dwd(dwd_res)
=#
NavierStokes = @decapode Flow2DQuantities begin
  (V, V̇, G)::Form1{X}
  (T, ρ, ṗ, p)::Form0{X}
  V̇ == neg₁(L₁′(V, V)) +
        div₁(kᵥ(Δ₁{X}(V) + third(d₀{X}(δ₁{X}(V)))), avg₀₁(ρ)) +
        d₀{X}(half(i₁′(V, V))) +
        neg₁(div₁(d₀{X}(p),avg₀₁(ρ))) +
        G
  ∂ₜ{Form1{X}}(V) == V̇
  ṗ == neg₀(⋆₀⁻¹{X}(L₀(V, ⋆₀{X}(p))))# + ⋆₀⁻¹{X}(dual_d₁{X}(⋆₁{X}(kᵨ(d₀{X}(ρ)))))
  ∂ₜ{Form0{X}}(p) == ṗ
end

Energy = @decapode Flow2DQuantities begin
  (V)::Form1{X}
  (ρ, p, T, Ṫ, Ṫₐ, Ṫ₁, bc₀)::Form0{X}

  ρ == div₀(p, R₀(T))
  Ṫₐ == neg₀(⋆₀⁻¹{X}(L₀(V, ⋆₀{X}(T))))
  ∂ₜₐ(Ṫₐ) == bc₀
  Ṫ == Ṫₐ + Ṫ₁
  ∂ₜ{Form0{X}}(T) == Ṫ
end

BoundaryConditions = @decapode Flow2DQuantities begin
  (V, V̇, bc₁)::Form1{X}
  (Ṫ, ṗ, bc₀)::Form0{X}
  # no-slip edges
  ∂ᵥ(V̇) == bc₁
  # No change on left/right boundaries
  ∂ₜ(Ṫ) == bc₀
  ∂ₚ(ṗ) == bc₀
end

dwd_res = diag2dwd(NavierStokes)
draw_equations(NavierStokes)
draw_dwd(dwd_res)

compose_heat_xfer = @relation (V, ρ) begin
  flow(V, V̇, T, ρ, ṗ, p)
  energy(Ṫ, V, ρ, p, T, Ṫ₁)
  diffusion(T, Ṫ₁)
  bcs(Ṫ, ṗ, V, V̇)
end

HeatXfer = oapply(compose_heat_xfer,
                  [OpenDiagram(NavierStokes, [:V, :V̇, :T, :ρ, :ṗ, :p]),
                   OpenDiagram(Energy, [:Ṫ, :V, :ρ, :p, :T, :Ṫ₁]),
                   OpenDiagram(Diffusion, [:T, :Ṫ]),
                   OpenDiagram(BoundaryConditions, [:Ṫ, :ṗ, :V, :V̇])])

draw_equations(HeatXfer.functor)

new_dwd = diag2dwd(HeatXfer.functor, in_vars=[:T, :p, :V, :G])
draw_dwd(new_dwd)

max_stars_repo_name: AlgebraicJulia/DECAPODES-Benchmarks

max_issues_repo_name: AlgebraicJulia/DECAPODES-Benchmarks

max_forks_repo_name: AlgebraicJulia/DECAPODES-Benchmarks


text:
using MomentTensors, Test, StaticArrays

# Comparison with example at http://www.larskrieger.de/mopad/#getall
@testset "Decomposition" begin
    let m = MT([ 3 -5 10
                -5  1  4
                10  4  2])
        d = decompose(m)
        @testset "Has field $f" for f in (
                :iso, :dev, :dc, :clvd, :iso_m0, :dev_m0, :prop_iso, :prop_dev,
                :prop_dc, :prop_clvd, :m0)
            @test f in fieldnames(typeof(d))
        end
        @test d.iso == MT(2, 2, 2, 0, 0, 0)
        @test d.prop_iso ≈ 0.13 atol=0.01
        @test d.dev == MT([ 1 -5 10
                           -5 -1  4
                           10  4  0])
        @test d.prop_dev ≈ 0.87 atol=0.01
        @test d.dc ≈ MT([1.36 -2.97  7.30
                        -2.97 -1.77  1.95
                         7.30  1.95  0.41]) atol=0.01
        @test d.prop_dc ≈ 0.55 atol=0.01
        @test d.clvd ≈ MT([-0.36 -2.03  2.70
                           -2.03  0.77  2.05
                            2.70  2.05 -0.41]) atol=0.01
        @test d.prop_clvd ≈ 0.32 atol=0.01
        @test d.m0 ≈ 14.9031089939
        @test mw(d.m0) ≈ -5.25114874821
        @test eps_non_dc(m) ≈ -0.34250995536253115
    end

    @testset "_sortperm_abs_3" begin
        @test MomentTensors._sortperm_abs_3(SVector(1.0, 2.0, 3.0)) == SVector(1, 2, 3)
        @test MomentTensors._sortperm_abs_3(SVector(-3.0, -2.0, 1)) == SVector(3, 2, 1)
    end
end


text:
#using Pkg; Pkg.activate("..")

using OceanTurb, Printf, StaticArrays, BenchmarkTools

Ns = (128,)

function atomic_forward_euler_benchmark(m, Δt)
    @printf "% 24s:" "update"
    @btime OceanTurb.update!($(m.bcs), $(m.timestepper.eqn), $(m.solution), $m)

    @printf "% 24s:" "calc explicit rhs"
    @btime OceanTurb.calc_explicit_rhs!($(m.timestepper.rhs), $(m.timestepper.eqn), $(m.solution), $m)

    @printf "% 24s:" "update solution"
    @btime OceanTurb.forward_euler_step!($(m.timestepper.rhs), $(m.solution), $Δt)

    @printf "% 24s:" "tick"
    @btime OceanTurb.tick!($(m.clock), $Δt)
    return nothing
end

function atomic_backward_euler_benchmark(m, Δt)
    @printf "% 24s:" "update"
    @btime OceanTurb.update!($(m.bcs), $(m.timestepper.eqn), $(m.solution), $m)

    @printf "% 24s:" "calc implicit rhs"
    @btime OceanTurb.calc_implicit_rhs!($(m.timestepper.rhs), $(m.timestepper.eqn), $(m.solution), $m)

    @printf "% 24s:" "calc diffusive lhs"
    @btime OceanTurb.calc_diffusive_lhs!($(m.timestepper.lhs), $(m.timestepper.eqn.M), $(m.timestepper.eqn.K), $(m.solution), $Δt, $m)

    @printf "% 24s:" "update solution"
    @btime OceanTurb.backward_euler_step!($(m.timestepper.rhs), $(m.timestepper.lhs), $(m.solution), $Δt)

    @printf "% 24s:" "tick"
    @btime OceanTurb.tick!($(m.clock), $Δt)
    return nothing
end

@printf "\nAtomic timestepping benchmarking of default configuration...\n"
Δt = 1.1
for N in Ns

    @printf "N: % 4d, forward Euler\n" N
    m = KPP.Model(N=N, stepper=:ForwardEuler)
    atomic_forward_euler_benchmark(m, Δt)

    @printf "N: % 4d, backward Euler\n" N
    m = KPP.Model(N=N, stepper=:BackwardEuler)
    atomic_backward_euler_benchmark(m, Δt)
end

@printf "\nAtomic timestepping benchmarking of configuration with one constant boundary condition...\n"
Δt = 1.1
for N in Ns

    @printf "N: % 4d, forward Euler\n" N
    m = KPP.Model(N=N, stepper=:ForwardEuler)
    m.bcs.T.top = BoundaryCondition(Flux, 1.0)
    atomic_forward_euler_benchmark(m, Δt)

    @printf "N: % 4d, backward Euler\n" N
    m = KPP.Model(N=N, stepper=:BackwardEuler)
    m.bcs.T.top = BoundaryCondition(Flux, 1.0)
    atomic_backward_euler_benchmark(m, Δt)
end


text:
using LinearAlgebra

function find_extreme_points(points::Array{Vector{T}, 1},
                             zmin::Union{Missing, Vector{T}}=missing) where T <: Real
    M = length(points[1])

    @show(zmin)
    if ismissing(zmin)
        zmin = zeros(M)
    end
    @show zmin

    extremepts = [-1 for j = 1:M]
    for j = 1:M
        w = 0.00001 * ones(M)
        w[j] = 1.0
        minval = -Inf
        for (i, p) in enumerate(points)
            val = maximum((p-zmin) ./ w)
            if val > minval
                extremepts[j] = i
                minval = val
            end
        end
    end

    return extremepts
end

function find_hyperplane_params(points::Array{Vector{T}, 1}, b=1) where T <: Real
    M = length(points[1])

    interceptvals = zeros(M)

    Z = transpose(reduce(hcat, points))
    bv = b * ones(M)

    a = Z\bv
    return a
end

function find_axis_intercepts(points::Array{Vector{T}, 1}, b=1) where T <: Real
    return 1 ./ find_hyperplane_params(points, b)
end

function associate_points(
        canpoints::Array{Vector{T}, 1},
        refpoints::Array{Vector{T}, 1}) where T <: Real
    N = length(canpoints)
    R = length(refpoints)
    # for each candidate point, want to find the closest reference point,
    # measured as perpendicular distance to line through ref point and origin

    closest_ref_index = zeros(Int64, N)
    closest_ref_dist = zeros(Float64, N)
    dists = zeros(R)
    for n = 1:N
        s = canpoints[n]
        for r = 1:R
            w = refpoints[r]
            dists[r] = norm( s - (dot(w, s)/dot(w, w)) * w )
        end
        closest_ref_index[n] = argmin(dists)
        closest_ref_dist[n] = dists[closest_ref_index[n]]
    end
    return closest_ref_index, closest_ref_dist
end

max_stars_repo_name: phawthorne/echidna

max_issues_repo_name: phawthorne/echidna

max_forks_repo_name: phawthorne/echidna


text:
@testset "MultipleSequenceAlignment" begin

    @testset "Type Hierarchy" begin

        @test AbstractAlignedObject <: AbstractMatrix{Residue}
        @test AbstractMultipleSequenceAlignment <: AbstractAlignedObject
        @test AbstractAlignedSequence <: AbstractAlignedObject
        @test MultipleSequenceAlignment <: AbstractMultipleSequenceAlignment
        @test AnnotatedMultipleSequenceAlignment <: AbstractMultipleSequenceAlignment
        @test AlignedSequence <: AbstractAlignedSequence
        @test AnnotatedAlignedSequence <: AbstractAlignedSequence

        @test AnnotatedMultipleSequenceAlignment <: AnnotatedAlignedObject
        @test AnnotatedAlignedSequence <: AnnotatedAlignedObject

        @test MultipleSequenceAlignment <: UnannotatedAlignedObject
        @test AlignedSequence <: UnannotatedAlignedObject
    end

    @testset "Creation" begin

        M = reshape(reinterpret(Residue,collect(1:21)),(3,7))
        m = NamedArray(M)

        T = permutedims(M, [2,1])

        S = reshape(reinterpret(Residue,collect(1:21)),(1,21))
        s = NamedArray(S)

        @test MultipleSequenceAlignment(M) == MultipleSequenceAlignment(m)
        @test AnnotatedMultipleSequenceAlignment(M) == AnnotatedMultipleSequenceAlignment(m)

        @test AlignedSequence(S) == AlignedSequence(s)
        @test AnnotatedAlignedSequence(S) == AnnotatedAlignedSequence(s)

        @test MultipleSequenceAlignment(M) != MultipleSequenceAlignment(T)
        @test AnnotatedMultipleSequenceAlignment(M) != AnnotatedMultipleSequenceAlignment(T)

        @test_throws AssertionError AlignedSequence(M)
        @test_throws AssertionError AnnotatedAlignedSequence(M)
    end

    @testset "MSA & sequences" begin

        M = reshape(reinterpret(Residue,collect(1:21)),(3,7))
        msa = MultipleSequenceAlignment(M)
        annotated_msa = AnnotatedMultipleSequenceAlignment(M)

        S = reshape(reinterpret(Residue,collect(1:21)),(1,21))
        sequence = AlignedSequence(S)
        annotated_sequence = AnnotatedAlignedSequence(S)

        @testset "Getters" begin

            @test isa(annotations(annotated_msa), Annotations)
            @test isa(annotations(annotated_sequence), Annotations)

            for object in (msa, annotated_msa, sequence, annotated_sequence)
                @test isa(namedmatrix(object), NamedResidueMatrix{Array{Residue,2}})
            end
        end

        @testset "Dimension names" begin

            for object in (msa, annotated_msa, sequence, annotated_sequence)
                @test dimnames(namedmatrix(object)) == ["Seq","Col"]
            end
        end

        @testset "Convert" begin

            msa2annot = convert(AnnotatedMultipleSequenceAlignment, msa)
            annot2msa = convert(MultipleSequenceAlignment, annotated_msa)

            seq2annot = convert(AnnotatedAlignedSequence, sequence)
            annot2seq = convert(AlignedSequence, annotated_sequence)

            @test isa(annotations(msa2annot), Annotations)
            @test isa(annotations(seq2annot), Annotations)

            @test_throws MethodError annotations(annot2msa)
            @test_throws MethodError annotations(annot2seq)

            @test namedmatrix(msa2annot) == namedmatrix(msa)
            @test namedmatrix(annot2msa) == namedmatrix(annotated_msa)
            @test namedmatrix(seq2annot) == namedmatrix(sequence)
            @test namedmatrix(annot2seq) == namedmatrix(annotated_sequence)
        end

        @testset "AbstractArray Interface" begin

            @test size(msa) == (3,7)
            @test size(annotated_msa) == (3,7)

            @test size(sequence) == (1,21)
            @test size(annotated_sequence) == (1,21)

            for object in (msa, annotated_msa, sequence, annotated_sequence)
                @test length(object) == 21
            end
        end

        @testset "Indexing" begin

            for object in (msa, annotated_msa)
                for i in 1:3, j in 1:7
                    @test object[string(i),string(j)] == object[i,j]
                end
            end

            for object in (sequence, annotated_sequence)
                for j in 1:7
                    @test object["1",string(j)] == object[1,j]
                    @test object[j] == object[1,j]
                    # Special sequence indexing:
                    @test object[string(j)] == object[1,j]
                end
            end
        end

        @testset "Show" begin

            out = IOBuffer()

            show(out, msa)
            str = String(take!(out))
            @test startswith(str,"MultipleSequenceAlignment : ")
            @test occursin("Seq", str)
            @test occursin("Col", str)
            @test length(split(str,'\n')) == 6

            show(out, annotated_msa)
            str = String(take!(out))
            @test startswith(str,"AnnotatedMultipleSequenceAlignment with 0 annotations : ")
            @test occursin("Seq", str)
            @test occursin("Col", str)
            @test length(split(str,'\n')) == 6

            show(out, sequence)
            str = String(take!(out))
            @test startswith(str,"AlignedSequence : ")
            @test occursin("Seq", str)
            @test occursin("Col", str)
            @test length(split(str,'\n')) == 4

            show(out, annotated_sequence)
            str = String(take!(out))
            @test startswith(str,"AnnotatedAlignedSequence with 0 annotations : ")
            @test occursin("Seq", str)
            @test occursin("Col", str)
            @test length(split(str,'\n')) == 4
        end

        @testset "Transpose" begin

            @test transpose(msa) == permutedims(msa)
            @test size(transpose(msa)) == (7,3)
            @test size(transpose(annotated_msa)) == (7,3)

            @test transpose(sequence) == permutedims(sequence)
            @test size(transpose(sequence)) == (21,1)
            @test size(transpose(annotated_sequence)) == (21,1)
        end

        @testset "Get residues" begin

            @test getresidues(msa) == M
            @test getresidues(annotated_msa) == M
            @test getresidues(sequence) == S
            @test getresidues(annotated_sequence) == S

            for object in (msa, annotated_msa, sequence, annotated_sequence)
                @test isa(getresidues(object), Matrix{Residue})
            end

            for object in (msa, annotated_msa, sequence, annotated_sequence)
                @test getresiduesequences(msa) == [res"ADEIMSY",
                                                   res"RCGLFTV",
                                                   res"NQHKPW-"]
            end
        end

        @testset "Size" begin

            for object in (M, NamedArray(M), msa, annotated_msa)
                @test ncolumns(object) == 7
                @test nsequences(object) == 3
            end

            for object in (S, NamedArray(S), sequence, annotated_sequence)
                @test ncolumns(object) == 21
                @test nsequences(object) == 1
            end
        end

        @testset "Get sequences" begin

            for object in (msa, annotated_msa)
                for seq in 1:3
                    @test getsequence(object, string(seq)) == getsequence(msa, seq)
                    @test size(getsequence(msa, seq)) == (1,7)
                end
            end
        end

        @testset "Sequence names" begin

            for object in (M, NamedArray(M), msa, annotated_msa)
                @test sequencenames(object) == ["1","2","3"]
            end
        end

        @testset "Column names" begin

            for object in (M, NamedArray(M), msa, annotated_msa)
                @test columnnames(object) == ["1","2","3","4","5","6","7"]
            end
        end

        @testset "Column mapping" begin

            for object in (NamedArray(M), msa, annotated_msa)
                @test getcolumnmapping(object) == [1,2,3,4,5,6,7]
            end
        end

        @testset "Sequence as string" begin

            for object in (M, NamedArray(M), msa, annotated_msa)
                @test stringsequence(object, 1) == "ADEIMSY"
                @test stringsequence(getsequence(object,1)) == "ADEIMSY"
            end

            for object in (NamedArray(M), msa, annotated_msa)
                @test stringsequence(msa, "1") == "ADEIMSY"
            end
        end

        @testset "Copy and setindex!" begin

            copy_msa = copy(msa)
            deepcopy_msa = deepcopy(msa)

            copy_annotated_msa = copy(annotated_msa)
            deepcopy_annotated_msa = deepcopy(annotated_msa)

            for x in [copy_msa, deepcopy_msa, copy_annotated_msa, deepcopy_annotated_msa]
                x[1,:] = res"YSMIEDA"
                @test vec(x[1,:]) == res"YSMIEDA"

                x["2",:] = res"YSMIEDA"
                @test vec(x["2",:]) == res"YSMIEDA"

                x[:,1] = res"YYY"
                @test vec(x[:,1]) == res"YYY"

                x[:,"2"] = res"YYY"
                @test vec(x[:,"2"]) == res"YYY"

                x[end,end] = 'X'
                @test x[end,end] == XAA

                @test length(unique(x)) != 21
            end

            copy_seq = copy(sequence)
            deepcopy_seq = deepcopy(sequence)

            copy_annotated_seq = copy(annotated_sequence)
            deepcopy_annotated_seq = deepcopy(annotated_sequence)

            for x in [copy_seq, deepcopy_seq, copy_annotated_seq, deepcopy_annotated_seq]
                x[1] = XAA
                @test x[1] == XAA

                x["1","2"] = XAA
                @test x["1","2"] == XAA

                x[end] = 'X'
                @test x[end] == XAA

                # Special setindex! for sequences:
                x["3"] = GAP
                @test x["3"] == GAP

                @test length(unique(x)) == 19
            end

            @test length(unique(msa)) == 21
            @test length(unique(annotated_msa)) == 21
            @test length(unique(sequence)) == 21
            @test length(unique(annotated_sequence)) == 21
        end
    end
end


text:
include("dbloader.jl")

# some basic features: distance to destination, available beer styles, distance to starting point
# all features are precomputed and L2-normalized
struct Features
    pairwise_distance::Matrix{Float64}
    beer_count::Vector{Float64}
    home_distance::Vector{Float64}
    scores::Vector{Float64}
    dim::Int
    function Features(coord::Coordinates, locinfo::LocationInfo)
        pairwise_distance = mapslices(normalize, coord.pairwise; dims=1) * -1
        beer_count = normalize(size.(collect(locinfo.beers), 1))
        home_distance = Vector{Float64}(undef, coord.dim)
        scores = Vector{Float64}(undef, coord.dim)
        new(pairwise_distance, beer_count, home_distance, scores, 3)
    end
end

function update!(feat::Features, pointwise)
    copy!(feat.home_distance, pointwise)
    normalize!(feat.home_distance)
    feat.home_distance .*= -1
end

function calc_scores!(f::Features, betas, id)
    @inbounds for i ∈ eachindex(f.scores)
        f.scores[i] = f.home_distance[i]*betas[1] + f.pairwise_distance[id, i]*betas[2] + f.beer_count[i]*betas[3]
    end
end

# different implementation needed when not traveling brewery->brewery
function calc_scores!(f::Features, betas)
    @inbounds for i ∈ eachindex(f.scores)
        f.scores[i] = f.home_distance[i]*betas[1] + f.home_distance[i]*betas[2] + f.beer_count[i]*betas[3]
    end
end