text:
# by Lazaro Alonso
using GLMakie, FileIO
GLMakie.activate!()
# See following links for references:
# https://twitter.com/AnsonBiggs/status/1444823816031510529/photo/1
# https://github.com/JuliaAnimators/Javis.jl/blob/master/examples/gravities.jl
let

    diameters = Dict(
        "Mercury" => 4879,
        "Venus" => 12104,
        "Earth" => 12756,
        "Moon" => 3475,
        "Mars" => 6792,
        "Jupiter" => 142984,
        "Saturn" => 120536,
        "Uranus" => 51118,
        "Neptune" => 49528,
        "Pluto" => 2370
    )
    struct Body
        position::Int
        name::String
        gravity::Real
        color::String
        radius::Real
    end
    Body(pos, name, grav, color) = Body(pos, name, grav, color, log(diameters[name])*2)
    bodies = [
        Body(1, "Mercury", 3.7, "snow4")
        Body(2, "Venus", 8.9, "navajowhite")
        Body(3, "Earth", 9.8, "lightskyblue")
        Body(4, "Moon", 1.6, "gainsboro")
        Body(5, "Mars", 3.7, "orangered")
        Body(6, "Jupiter", 23.1, "olive")
        Body(7, "Saturn", 9, "burlywood")
        Body(8, "Uranus", 8.7, "cyan3")
        Body(9, "Neptune", 11, "dodgerblue")
        Body(10, "Pluto", 0.7, "rosybrown4")
    ]
    framerate = 30
    height = 1000 # meters

    frames = let
        # Get planet with slowest acceleration
        slowest = [p.gravity for p in bodies] |> minimum
        # Kinematic equation to determine seconds to fall from height
        time = 2 * height / slowest |> sqrt
        # Calculate total frames and add a few seconds at the end to display final result
        ceil(Int, time * framerate) + framerate * 5
    end

    textsize = 28
    f = Node(1)
    set_theme!(theme_black())
    fig = Figure(resolution = (2400, 1600))
    axs = [Axis(fig[i,j]) for j in 1:11, i in 1:3]
    pathBodies = "./scripts/bodies/"
    for (idx, body) in enumerate(bodies)
        # Kinematic equations to calculate how many frames each planet will be active
        t_final = 2 * height / body.gravity |> sqrt
        v_final = body.gravity * t_final
        frame_final = ceil(Int, t_final * framerate)
        # Calculate time of current frame for kinematic math
        time = @lift($f / framerate)
        # Calculate planets current position
        body_y = @lift(height - 0.5 * body.gravity * $time^2)
        sec = @lift($body_y <=0  ? round(t_final, digits=1) : round($time, digits = 1))
        velocity = @lift($body_y <=0  ? round(v_final, digits=1) : round(body.gravity * $time, digits = 1))
        body_y = @lift($body_y <=0 ? 0 : $body_y)

        # Set the planets position
        scatter!(axs[idx+1, 2], @lift(Point2f(1, $body_y)), color = body.color, markersize = 25)
        text!(axs[idx+1, 1], @lift(string($sec, "s")),
            position = Point2f(1,1),
            textsize = textsize, color = @lift($body_y <=0 ? "orange" : "grey90"),
            align = (:center, :center))
        text!(axs[idx+1, 1],
            @lift(string($velocity, "m/s")),
            position = Point2f(1,0), textsize = textsize,
            color = @lift($body_y <=0 ? "orange" : "grey90"),
            align = (:center, :center))

        ylims!(axs[idx+1,2], -15, height+15)

        bname = pathBodies*body.name * ".png"
        img = rotr90(load(bname))
        image!(axs[idx+1, 3], img)
        mx = maximum(size(img))
        mn = minimum(size(img))
        limits!(axs[idx+1,3], -50, mx, -50, mx)
        scatter!(axs[idx+1,3], Point2f(mx÷2, mx÷2 -50), color = :transparent,
            strokecolor = body.color, markersize = 170, strokewidth = @lift($body_y <=0 ? 2 : 0))
        # Set text that is static during entire planet translation
        text!(axs[idx+1, 1], body.name, position = Point2f(1, 3),
            textsize = textsize, align = (:center, :center))
        text!(axs[idx+1, 1], string(body.gravity, "m/s²"),
            position = Point2f(1, 2), textsize = textsize,
            align = (:center, :center))
        limits!(axs[idx+1, 1], 0, 3, -0.5, 3.5)
    end
    ylims!(axs[1,2], -15, height+15)
    text!(axs[1,2], "1 Km", position = Point2f(0, 1000),
        textsize = textsize, align = (:left, :center))
    text!(axs[1,2], "0 Km", position = Point2f(0, 0),
        textsize = textsize, align = (:left, :center))

    text!(axs[1,1], "Body:", position = Point2f(0, 3),
        textsize = textsize, align = (:left, :center))
    text!(axs[1,1], "Acceleration:", position = Point2f(0, 2),
        textsize = textsize, align = (:left, :center))
    text!(axs[1,1], "Time:", position = Point2f(0, 1),
        textsize = textsize, align = (:left, :center))
    text!(axs[1,1], "Velocity:", position = Point2f(0, 0),
        textsize = textsize, align = (:left, :center))
    limits!(axs[1, 1], 0, 3, -0.5, 3.5)

    [hidedecorations!(axs[idx,k]) for idx in 1:11, k in 1:3]
    [hidespines!(axs[idx,k]) for idx in 1:11, k in 1:3]

    Label(fig[0,:], "Falling Bodies in the Solar System",
        tellwidth = false, textsize = 70, color = :grey50)
    Label(fig[end+1,:], "Twitter: @LazarusAlon", color =:white,
        tellwidth = false, textsize = 30, halign = :right,
        font ="noto-sands-bold")
    Label(fig[1,:], "using Makie", color = "#f7c530",
        tellwidth = false, textsize = 60, halign = :left,
        font ="noto-sands-bold")
    rowgap!(fig.layout, 0)
    colgap!(fig.layout, 0)
    rowsize!(fig.layout, 2, Auto(0.2))
    rowsize!(fig.layout, 4, Auto(0.2))
    fig

    record(fig, joinpath(@__DIR__, "output", "gravities.mp4"), framerate = 2*framerate) do io
        for i in 1:frames
            f[] = i
            recordframe!(io)  # record a new frame
        end
    end
end

max_stars_repo_name: lazarusA/BeautifulMakie

max_issues_repo_name: lazarusA/BeautifulMakie

max_forks_repo_name: lazarusA/BeautifulMakie


text:
export SQG

"""
    SQG

Data structure to store arrays and parameters for the
Surface-Quasi-Geostrophic model

- rho : Background density
- g : Gravity
- buoyancy_freq_n : Background stratification
- odg_b : Amplitude of the buoyancy
"""
mutable struct SQG

    grid :: Grid
    rho :: Float64
    g :: Float64
    buoyancy_freq_n :: Float64
    odg_b :: Float64
    hv_order :: Int
    hv_val :: Float64

    b̂ :: Array{ComplexF64, 2}
    ψ̂ :: Array{ComplexF64, 2}

    u_x :: Array{Float64, 2}
    u_y :: Array{Float64, 2}
    û_x :: Array{ComplexF64, 2}
    û_y :: Array{ComplexF64, 2}

    a :: Array{Float64, 2}
    â :: Array{ComplexF64, 2}

    function SQG( grid, f0 )

        rho=1e3
        g = 9.81
        buoyancy_freq_n = 3 * f0
        odg_b = 1e-3
        hv_order = 8
        hv_val = 0

        nx, ny = grid.nx, grid.ny

        b̂ = zeros(ComplexF64, nx÷2+1, ny)
        ψ̂ = zeros(ComplexF64, nx÷2+1, ny)

        u_x = zeros(nx, ny)
        u_y = zeros(nx, ny)
        û_x = zeros(ComplexF64, nx÷2+1, ny)
        û_y = zeros(ComplexF64, nx÷2+1, ny)

        a = zeros(Float64, nx, ny)
        â = zeros(ComplexF64, nx÷2+1, ny)

        new( grid, rho, g, buoyancy_freq_n, odg_b, hv_order, hv_val,
             b̂, ψ̂, u_x, u_y, û_x, û_y, a, â)

    end

end

export update_velocities!

"""
    update_velocities!(u_x, u_y, model, b̂)

Compute the streamfunction and the velocity from the
buoyancy according to the Surface-Quasi-Geostrophic model.

"""
function update_velocities!( sqg :: SQG )

    nx, ny = sqg.grid.nx, sqg.grid.ny
    sqg.ψ̂ .= sqg.b̂ .* sqg.grid.on_k ./ sqg.buoyancy_freq_n

    sqg.û_x .= - 1im .* sqg.grid.ky .* sqg.ψ̂
    sqg.û_y	.=  1im .* sqg.grid.kx .* sqg.ψ̂

    sqg.u_x .= irfft(sqg.û_x, nx)
    sqg.u_y .= irfft(sqg.û_y, nx)

end

export update_advection_term!

"""
    update_advection_term(model)

Compute the Fourier transform of the partial derivative
of the advection term plus the hyperviscosity term

"""
function update_advection_term!(sqg)

    nx = sqg.grid.nx

    update_velocities!( sqg )

    # Advection term

    sqg.a .=  sqg.u_x .* irfft(-1im .* sqg.grid.kx .* sqg.b̂, nx)
    sqg.a .+= sqg.u_y .* irfft(-1im .* sqg.grid.ky .* sqg.b̂, nx)

    sqg.â .= rfft(sqg.a)

    # Summing Hyperviscosity term

    sqg.â .-= sqg.hv_val .* (sqg.grid.k.^sqg.hv_order) .* sqg.b̂

end

export TimeSolver

struct TimeSolver

    b̂ :: Array{ComplexF64, 2}
    db̂ :: Array{ComplexF64, 2}

    function TimeSolver( nx, ny)

        b̂ = zeros(ComplexF64, nx÷2+1, ny)
        db̂ = zeros(ComplexF64, nx÷2+1, ny)

        new( b̂, db̂ )

    end

end

export step!

"""
    step!(model, time_solver, dt)

Time integration of the advection equation of b
by 4th order Runge-Kutta method
"""
function step!( sqg :: SQG, rk4 :: TimeSolver, dt )

    rk4.b̂ .= sqg.b̂

    update_advection_term!(sqg)

    sqg.b̂ .= rk4.b̂ .+ sqg.â .* dt / 2
    rk4.db̂ .= sqg.â ./ 2

    update_advection_term!(sqg)

    sqg.b̂ .= rk4.b̂ .+ sqg.â .* dt / 2
    rk4.db̂ .+= sqg.â

    update_advection_term!(sqg)

    sqg.b̂ .= rk4.b̂ .+ sqg.â .* dt
    rk4.db̂ .+= sqg.â

    update_advection_term!(sqg)

    rk4.db̂ .+= sqg.â ./ 2

    sqg.b̂ .= rk4.b̂ .+ dt / 3 .* rk4.db̂

end

ext: jl


text:
"""
    Maybe{T}()
    Maybe(value::T)

A structure which can hold a value of type `T`, or represent the absence of a value.

The API is optimised for speed over memory usage, by allowing a function that may otherwise
return `Union{T, Nothing}` to instead always return `Maybe{T}`, and hence be type-stable.
"""
struct Maybe{T}
    # TODO see the changes in this merge request:
    #   https://github.com/iamed2/ResultTypes.jl/commit/24e0aa779a29376ce3320cd8f106c5416aaff4c1
    # It may be possible to improve performance by using Union{Some{T},Nothing} instead of
    # this partial initialisation approach.
    valid::Bool
    value::T

    Maybe{T}() where {T} = new{T}(false)
    Maybe(value::T) where {T} = new{T}(true, value)
end

"""
    valid(x::Maybe) -> Bool

Return true iff `x` holds a value.
"""
valid(x::Maybe) = x.valid

"""
    value(x::Maybe{T}) -> T

Returns the value stored in `x`, or throws an `ArgumentError` if `!valid(x)`.

Note that, in a tight loop, it is preferable to use a combination of calls to `valid` and
`unsafe_value`, as it will generate more optimal code.
"""
value(x::Maybe) = valid(x) ? unsafe_value(x) : throw(ArgumentError("$x has no value."))

"""
    unsafe_value(x::Maybe{T}) -> T

Returns the value stored in `x`.

It is "unsafe" when `!valid(x)`, in that the return value of this function is undefined. If
`T` is a reference type, calling this function will result in an `UndefRefError` being
thrown.
"""
unsafe_value(x::Maybe) = x.value

function Base.show(io::IO, x::Maybe{T}) where {T}
    if valid(x)
        print(io, "Maybe{$T}($(x.value))")
    else
        print(io, "Maybe{$T}()")
    end
end


text:
"""
    DBAclustResult(centers,clustids,result)

Holds results of a DBAclust run.

"""
mutable struct DBAclustResult{T}
    centers::T
    clustids::Array{Int}
    converged::Bool
    iterations::Int
    dbaresult::DBAResult
end

"""
    dbaclust(
        sequences,
        nclust::Int,
        dtwdist::DTWDistance;
        n_init::Int           = 1,
        iterations::Int       = 100,
        inner_iterations::Int = 10,
        rtol::Float64         = 1e-4,
        rtol_inner::Float64   = rtol,
        n_jobs::Int           = 1,
        show_progress::Bool   = true,
        store_trace::Bool     = true,
        i2min::AbstractVector = [],
        i2max::AbstractVector = [],
    )

# Arguments:
- `nclust`: Number of clsuters
- `n_init`: Number of initialization tries
- `inner_iterations`: Number of iterations in the inner alg.
- `i2min`: Bounds on the warping path
"""
function dbaclust(
    sequences,
    nclust::Int,
    dtwdist::DTWDistance;
    n_init::Int           = 1,
    iterations::Int       = 100,
    inner_iterations::Int = 10,
    rtol::Float64         = 1e-4,
    rtol_inner::Float64   = rtol,
    show_progress::Bool   = true,
    store_trace::Bool     = true,
    i2min::AbstractVector = [],
    i2max::AbstractVector = [],
)

    best_result = []
    best_cost = []
    for i = 1:n_init
        dbaclust_result = dbaclust_single(
            sequences,
            nclust,
            dtwdist;
            iterations = iterations,
            inner_iterations = inner_iterations,
            rtol = rtol,
            rtol_inner = rtol_inner,
            show_progress = show_progress,
            store_trace = store_trace,
            i2min = i2min,
            i2max = i2max,
        )
        if isempty(best_cost) || dbaclust_result.dbaresult.cost < best_cost
            best_result = deepcopy(dbaclust_result)
            best_cost = best_result.dbaresult.cost
        end
    end #1:n_init

    return best_result
end

"""
    avgseq, results = dbaclust_single(sequences, dist; kwargs...)

Perfoms a single DTW Barycenter Averaging (DBA) given a collection of `sequences`
and the current estimate of the average sequence.

Example usage:

    x = [1,2,2,3,3,4]
    y = [1,3,4]
    z = [1,2,2,4]
    avg,result = dba([x,y,z], DTW(3))
"""
function dbaclust_single(
    sequences::AbstractVector,
    nclust::Int,
    dtwdist::DTWDistance;
    init_centers::AbstractVector = dbaclust_initial_centers(
        sequences,
        nclust,
        dtwdist
    ),
    iterations::Int       = 100,
    inner_iterations::Int = 10,
    rtol::Float64         = 1e-4,
    rtol_inner::Float64   = rtol,
    show_progress::Bool   = true,
    store_trace::Bool     = true,
    i2min::AbstractVector = [],
    i2max::AbstractVector = [],
)

    T = floattype(eltype(sequences))
    # rename for convienence
    avgs   = init_centers
    N = length(avgs[1])

    # check initial centers have the same length
    if !all(length(a) == N for a in avgs)
        throw(ArgumentError("all initial centers should be the same length"))
    end

    # dimensions
    nseq      = length(sequences)
    maxseqlen = maximum([length(s) for s in sequences])

    # TODO switch to ntuples?
    counts    = [zeros(Int, N) for _ = 1:nclust]
    sums      = [Array{T}(undef,N) for _ = 1:nclust]

    # cluster assignments for each sequence
    clus_asgn = Array{Int}(undef,nseq)
    c         = 0

    # arrays storing path through dtw cost matrix
    i1, i2    = Int[], Int[]

    # variables storing optimization progress
    converged       = false
    iter            = 0
    inner_iter      = 0
    converged_inner = false
    last_cost       = Inf
    total_cost      = 0.0
    cost_trace      = Float64[]
    costs           = Array{Float64}(undef,nseq)

    # main loop ##
    if show_progress
        prog = ProgressMeter.ProgressThresh(rtol, 2)
        ProgressMeter.update!(
            prog,
            Inf;
            showvalues = [
                (:iteration, iter),
                (Symbol("max iteration"), iterations),
                (:cost, total_cost),
            ],
        )
    end#showprogress

    while !converged && iter < iterations

        # first, update cluster assignments based on nearest
        # centroid (measured by dtw distance). Keep track of
        # total cost (sum of all distances to centers).
        total_cost = 0.0
        for s = 1:nseq
            # process sequence s
            seq = sequences[s]

            # find cluster assignment for s
            costs[s] = Inf
            for c_ = 1:nclust
                # if one of the two is empty, use unconstrained window. If both are nonempty, but not the same lenght, distpath will throw error
                if isempty(i2min) && isempty(i2max)
                    cost, i1_, i2_ = distpath(dtwdist, avgs[c_], seq)
                else
                    cost, i1_, i2_ = distpath(dtwdist, avgs[c_], seq, i2min, i2max)
                end
                if cost < costs[s]
                    # store cluster, and match indices
                    c        = c_
                    i1       = i1_
                    i2       = i2_
                    costs[s] = cost
                end
            end

            # s was assigned to cluster c
            clus_asgn[s] = c
            cnt          = counts[c]
            sm           = sums[c]
            avg          = avgs[c]

            # update stats for barycentric average for
            # the assigned cluster
            for t = 1:length(i2)
                cnt[i1[t]] += 1
                sm[i1[t]]  += seq[i2[t]]
            end
        end

        # if any centers are unused, and reassign them to the sequences
        # with the highest cost
        unused = setdiff(1:nclust, unique(clus_asgn))
        if !isempty(unused)
            # reinitialize centers
            for c in unused
                avgs[c] = deepcopy(sequences[argmax(costs)])
                for s = 1:nseq
                    seq = sequences[s]
                    if isempty(i2min) && isempty(i2max)
                        cost, = distpath(dtwdist, avgs[c], seq)
                    else
                        cost, = distpath(dtwdist, avgs[c], seq, i2min, i2max)
                    end
                    if costs[s] > cost
                        costs[s] = cost
                    end
                end
            end
            # we need to reassign clusters, start iteration over
            continue
        end

        # store history of cost while optimizing (optional)
        total_cost = sum(costs)
        store_trace && push!(cost_trace, total_cost)

        # check convergence
        Δ = (last_cost - total_cost) / total_cost
        if Δ < rtol
            converged = true
        else
            last_cost = total_cost
        end

        # update barycenter estimates
        for (a, s, c) in zip(avgs, sums, counts)
            for t = 1:N
                c[t] == 0 && continue
                a[t] = s[t] / c[t]
            end
            # zero out sums and counts for next iteration
            s .= 0
            c .= 0
        end

        # add additional inner dba iterations

        for i = 1:nclust
            seqs            = view(sequences, clus_asgn .== i)
            inner_iter      = 0
            converged_inner = false
            oldcost         = 1.0e100
            while !converged_inner && inner_iter < inner_iterations
                newcost = dba_iteration!(
                    sums[i],
                    avgs[i],
                    counts[i],
                    seqs,
                    dtwdist;
                    i2min = i2min,
                    i2max = i2max,
                )
                copy!(avgs[i], sums[i])
                inner_iter += 1
                δ = (oldcost - newcost) / oldcost
                if δ < rtol_inner
                    converged_inner = true
                else
                    oldcost = newcost
                end
            end
        end

        # update progress bar
        iter += 1
        show_progress && ProgressMeter.update!(
            prog,
            Δ;
            showvalues = [
                (:iteration, iter),
                (Symbol("max iteration"), iterations),
                (:cost, total_cost),
            ],
        )
    end

    return DBAclustResult(
        avgs,
        clus_asgn,
        converged,
        iter,
        DBAResult(total_cost, converged_inner, inner_iter, cost_trace),
    )
end

"""
   dbaclust_initial_centers(sequences, nclust, dtwdist::DTWDistance)

Uses kmeans++ (but with dtw distance) to initialize the centers
for dba clustering.
"""
function dbaclust_initial_centers(
    sequences::AbstractVector,
    nclust::Int,
    dtwdist::DTWDistance,
)

    # number of sequences in dataset
    nseq          = length(sequences)
    # distance of each datapoint to each center
    dists         = zeros(nclust, nseq)
    # distances to closest center
    min_dists     = zeros(1, nseq)
    # choose a center uniformly at random
    center_ids    = zeros(Int, nclust)
    center_ids[1] = rand(1:nseq)

    # assign the rest of the centers
    for c = 1:(nclust-1)

        # first, compute distances for the previous center
        cent = sequences[center_ids[c]]
        for (i, seq) in enumerate(sequences)
            # this distance will be zero
            i == center_ids[c] && continue
            # else, compute dtw distance
            dists[c, i], = distpath(dtwdist, seq, cent)
        end

        # for each sequence, find distance to closest center
        minimum!(min_dists, dists[1:c, :])

        min_dists .= abs2.(min_dists)

        # sample the next center
        center_ids[c+1] = sample(1:nseq, Weights(view(min_dists, :)))
    end

    # return list of cluster centers
    return [deepcopy(sequences[c]) for c in center_ids]
end


text:
using .Snopt

export SNOPT

"""
    SNOPT(;options=Dict(), names=Snopt.Names(), warmstart=nothing)

Use Snopt as the optimizer

# Arguments
- `options::Dict`: options for Snopt.  see Snopt docs.
- `names::Snopt.Names`: custom names for function and variables.
- `warmstart::Snopt.Start`: a warmstart object (one of the outputs of Snopt.Outputs)
"""
struct SNOPT{T1,T2,T3} <: AbstractSolver
    options::T1
    names::T2
    warmstart::T3
end

SNOPT(;options=Dict(), names=Snopt.Names(), warmstart=nothing) = SNOPT(options, names, warmstart)

# wrapper for SNOPT
function optimize(snopt::SNOPT, cache, x0, lx, ux, lg, ug, rows, cols)

    function fsnopt!(g, df, dg, x, deriv)
        f = evaluate!(g, df, dg, x, cache)
        fail = false  # TODO
        # TODO: use deriv

        return f, fail
    end

    if !isnothing(snopt.warmstart)
        x0 = snopt.warmstart
    end

    return snopta(fsnopt!, x0, lx, ux, lg, ug, rows, cols,
        snopt.options, names=snopt.names)
end


text:
#_________________________________________________________________________#
#  Whale Optimization Algorithm (WOA) source codes demo 0.1               #
#                                                                         #
#                                                                         #
#  Programmer: Jesús Mejía                                                #
#  base on matlab code of S. Mirjalili                                    #
#                                                                         #
#                                                                         #
#   Main paper: S. Mirjalili, A. Lewis                                    #
#               The Whale Optimization Algorithm,                         #
#               Advances in Engineering Software , in press,              #
#               DOI: http://dx.doi.org/10.1016/j.advengsoft.2016.01.008   #
#                                                                         #
#_________________________________________________________________________#

mutable struct WOA <: AbstractParameters
  N::Int
end

"""
    WOA(;N = 30, information = Information(), options = Options())

Parameters for the Whale Optimization Algorithm. `N` is the population size (number of
whales).

# Example

```jldoctest
julia> f(x) = sum(x.^2)
f (generic function with 1 method)

julia> optimize(f, [-1 -1 -1; 1 1 1.0], WOA())
+=========== RESULT ==========+
  iteration: 500
    minimum: 3.9154600000000003e-100
  minimizer: [8.96670478694908e-52, -1.9291317455298046e-50, 4.3113080446722046e-51]
    f calls: 15000
 total time: 0.0134 s
+============================+

julia> optimize(f, [-1 -1 -1; 1 1 1.0], WOA(N = 100))
+=========== RESULT ==========+
  iteration: 500
    minimum: 1.41908e-145
  minimizer: [9.236161414012512e-74, -3.634919950380001e-73, 3.536831799149254e-74]
    f calls: 50000
 total time: 0.0588 s
+============================+

```
"""
function WOA(;N = 30, information = Information(), options = Options())
  parameters = WOA(N)

  Algorithm(
            parameters,
            information = information,
            options = options,
           )

end

function initialize!(
    status,
    parameters::WOA,
    problem::AbstractProblem,
    information::Information,
    options::Options,
    args...;
    kargs...
  )

  lb, ub= problem.bounds[1,:], problem.bounds[2,:]

  #Initialize the positions of search agents
  max_it = 500
  options.iterations = options.iterations == 0 ? max_it : options.iterations
  options.f_calls_limit = options.f_calls_limit == 0 ?
  options.iterations * parameters.N : options.f_calls_limit

  P = generate_population(parameters.N, problem)
  best_sol = deepcopy(get_best(P))
  status = State(best_sol, P)
  #=
  status.population = P
  status.best_sol = deepcopy(get_best(status.population))
  =#
  status.f_calls = parameters.N

  status

end

function update_state!(
    status::State,
    parameters::WOA,
    problem::AbstractProblem,
    information::Information,
    options::Options,
    args...;
    kargs...
  )

  # Return back the search agents that go beyond the boundaries of the search space
  Leader_score = minimum(status)
  Leader_pos = minimizer(status)

  Max_iter = options.iterations
  N = parameters.N
  D = size(problem.bounds, 2)
  t = status.iteration

  a=2-t*((2)/Max_iter) # a decreases linearly fron 2 to 0 in Eq. (2.3)

  # a2 linearly dicreases from -1 to -2 to calculate t in Eq. (3.12)
  a2=-1+t*((-1)/Max_iter)

  X_new = zeros(N,D)

  # Update the Position of search agents
  for i=1:N
    r1 = rand() # r1 is a random number in [0,1]
    r2 = rand() # r2 is a random number in [0,1]

    A = 2a*r1 - a  # Eq. (2.3) in the paper
    C = 2r2        # Eq. (2.4) in the paper

    b = 1               #  parameters in Eq. (2.5)
    l = (a2-1)*rand()+1   #  parameters in Eq. (2.5)

    p = rand()        # p in Eq. (2.6)

    x = copy(status.population[i].x)
    for j = 1:D
      if p < 0.5
        if abs(A) >= 1
          rand_leader_index = floor(Int, N* rand()+1)
          X_rand  = status.population[rand_leader_index].x
          D_X_rand= abs(C*X_rand[j] - x[j]  ) # Eq. (2.7)
          x[j] = X_rand[j]-A*D_X_rand       # Eq. (2.8)

        else
          D_Leader=abs(C*Leader_pos[j] - x[j] ) # Eq. (2.1)
          x[j] =Leader_pos[j]-A*D_Leader      # Eq. (2.2)
        end

      else

        distance2Leader=abs(Leader_pos[j]-x[j])
        # Eq. (2.5)
        x[j] =distance2Leader*exp(b.*l).*cos(l.*2*pi)+Leader_pos[j]

      end

    end # for j

    reset_to_violated_bounds!(x, problem.bounds)
    X_new[i,:] = x
  end # for i

  status.population = create_solutions(X_new, problem)
  sol = get_best(status.population)
  if is_better(sol, status.best_sol)
    status.best_sol = deepcopy(sol)
  end

end

function final_stage!(
    status::State,
    parameters::WOA,
    problem::AbstractProblem,
    information::Information,
    options::Options,
    args...;
    kargs...
  )
  status.final_time = time()
end

# The Whale Optimization Algorithm


text:
module SwitchingAnalysis

using Reexport
@reexport using DataFrames, CSV, CategoricalArrays
@reexport using Statistics, StatsBase,Plots, StatsPlots, KernelDensity, HypothesisTests
@reexport using GLM, MixedModels
@reexport using Survival

using Distributions

export mac_gdrive, linux_gdrive, files_dir, figs_dir, columns_types
export drug_colors, protocol_colors, Treatment_dict, Plotting_position
export Prew, Poutcome, Pobservations, Pprotocol, Pnext
export dropnan, dropnan!, Protocol_colors!, Drug_colors!, jump_missing, jump_NaN, skipmean,binquantile, TimeFromLeaving
export conf_ints, trim_conf_ints, trim_conf_ints!
export process_streaks
export count_cases, count_bouts!, process_bouts
export CIq, wilcoxon
export mediansurvival_analysis, survivalrate_algorythm, hazardrate_algorythm, function_analysis, KaplanMeier, survival_analysis
export KDensity, trim_dist, trim_dist!
export ODC, calculate_odc
export Likelyhood_Ratio_test
export summarize, effect_size, MVT, MVT_scatter, plot_wilcoxon, WebersLaw, plot_wilcoxon_odc, plot_odc, odc_regression, plot_deltafromrew, plot_QODC

include("constants.jl")
include("process_bouts.jl")
include("process_streaks.jl")
include("SignedRank_analysis.jl")
include("survival_analysis.jl")
include("plotting_functions.jl")
include("probabilities.jl")
include("utilities.jl")
include("KDensity.jl")
include("Omission_Duty_Cycle.jl")
include("DOFtests.jl")

end # module


text:
#fdtd.jl
using DelimitedFiles

#
# single 1d fdtd step in-place
#
function FDTDStep1D!(e::Array{T}, h::Array{T}, invϵ::Array{T}, Δx::T, Δt::T) where T
    ne = length(e)
    nh = length(h)
    dtdxFactor = Δt / Δx

    #e-field pass
    e[1] = e[end] = zero(T) #ensure boundary condition
    @inbounds for i = 2:(ne - 1) #spare first and last element
        e[i] += dtdxFactor * invϵ[i] * (h[i - 1] - h[i])
    end
    #h-field pass
    @inbounds for i = 1:nh
        h[i] += dtdxFactor * (e[i] - e[i + 1])
    end
end

#
# fdtd 1d simulation
#
function FDTD1D(steps::Int, snapshots::Array{Int}, n::Int, Δx::Float64, Δt::Float64, ϵBorder::Float64 = Inf)

    # spatial coordinates (not needed for calculation)
    x0 = (n - 1.0) / 2.0 * Δx # artificial midpoint
    x_e = collect(range(-x0 / 2.0         , length = n,     step = Δx))
    x_h = collect(range(-x0 / 2.0 + Δx / 2, length = n - 1, step = Δx))

    #initial fields: gaussian packet
    e = exp.(-0.5e2 * (n - 1) * Δx .* x_e .^ 2)
    h = exp.(-0.5e2 * (n - 2) * Δx .* x_h .^ 2)
    invϵ = map(x -> if (x > ϵBorder) 1.0 / 4.0 else 1.0 end, x_e)

    #initialize the result variables and snapshot-track-variables
    snapshotIndex = 1
    nextSnapshotIndex = snapshots[1]
    eResult = Array{Float64}(undef, length(snapshots), n)
    hResult = Array{Float64}(undef, length(snapshots), n - 1)

    #lets go
    for i = 1:steps
        FDTDStep1D!(e, h, invϵ, Δx, Δt) #magic

        if i == nextSnapshotIndex #if we reached a snapshot point:
            #copy the current fields
            for m = 1:n
                eResult[snapshotIndex, m] = e[m]
            end
            for m = 1:(n - 1)
                hResult[snapshotIndex, m] = h[m]
            end
            if snapshotIndex < length(snapshots) #and advance the snapshot-track-variables
                snapshotIndex += 1
                nextSnapshotIndex = snapshots[snapshotIndex]
            end
        end
    end

    #return out the results
    (efield = eResult, xe = x_e, hfield = hResult, xh = x_h, snaps = snapshots)
end

#
# save result in csv file
#
function SaveFDTDResult(resultFile, result)
    snapshots = result.snaps
    open(resultFile, "w") do io
        writedlm(io, [result.snaps, result.xe, result.xh], ',')
        writedlm(io, result.efield, ',')
        writedlm(io, result.hfield, ',')
    end
end

max_stars_repo_name: JulienKluge/ComputationalPhotonics

max_issues_repo_name: JulienKluge/ComputationalPhotonics

max_forks_repo_name: JulienKluge/ComputationalPhotonics


text:
#---------------------------------------------#
# Geometrical transformations for the segment #
#---------------------------------------------#

function getSegTransform(A::Array{Float64, 1}, B::Array{Float64,1})
    #= Returns (u0, u) the transform data from [-1, 1] to a 2d segment [A, B] =#
    return (A + B)/2, (B-A)/2
end

 @inline function segRefToSeg2d(xref::Float64, A::Array{Float64,1},
    B::Array{Float64, 1})
    #= Transformation from the 1d [-1, 1] reference segment to a 2-dimensional
        segment [A, B] where A and B are the nodes =#
    return A + (xref + 1)/2*(B-A)
end

@inline function segRefToSeg2d(xref::Float64, u0::Array{FLoat64, 1},
    u::Array{Float64})
    #= Affine transformation ans = u0 + xref*u where u0 and u are obtained from
        getSegTransform() =#
    return u0 + xref*u
end

#----------------------------------------------#
# Geometrical transformations for the triangle #
#----------------------------------------------#
function getTriTransform(T::Array{Float64,2})
    #= Returns the matrices required for the computation of the transformation
        from the triangle {(0,0),(1,0),(1,0)} to {A, B, C} where A, B and C are
        stored columnwise in T =#
    return T[:, 1], [T[1,2]-T[1,1] T[1,3]-T[1,1];T[2,2]-T[2,1] T[2,3]-T[2,1]]
end

@inline function triRefToAny(xref::Array{Float64,1}, T::Array{Float64, 2})
    #= Transformation from the reference triangle {(0,0),(1,0),(1,0)} to
        T={A, B, C} =#
    return T[:, 1] + [T[1,2]-T[1,1] T[1,3]-T[1,1];T[2,2]-T[2,1] T[2,3]-T[2,1]]*xref
end

@inline function triRefToAny(xref::Array{Float64}, u0::Array{Float64,1},
    u::Array{Float64})
    #= Transformation from the reference triangle {(0,0),(1,0),(1,0)} to
        T={A, B, C} using the matrices returned by getTriTransform() =#
    return u0 + u*xref
end


text:
const KeyType = Union{Symbol, Int}

const Arguments = Vector{Any}
const NamedArguments = Dictionary{Symbol, Any}
const MixedArguments = Dictionary{KeyType, Any}

# make a copy with distinct keys
set(d::AbstractDictionary, ps::Pair...) = merge(d, dictionary(ps))

function filterkeys(f, d::AbstractDictionary)
    idxs = findall(f, keys(d))
    return getindices(d, idxs)
end

function separate(f, d::AbstractDictionary)
    d1, d2 = empty(d), empty(d)
    for (k, v) in pairs(d)
        target = ifelse(f(v), d1, d2)
        insert!(target, k, v)
    end
    return d1, d2
end

valid_options(; options...) = valid_options(values(options))

function valid_options(nt)
    ks = filter(k -> nt[k] !== automatic, keys(nt))
    return NamedTuple{ks}(nt)
end


text:
# using PowerModels
using JuMP
using Suppressor
import Ipopt
#TODO: move RCall to run_qc_relax
@suppress using RCall # use RCall since hitandrun has an efficient implementation there,
## the prefered Julia lib would be MAMBA but it does not have this sampling. See:
## https://mambajl.readthedocs.io/en/latest/samplers.html?highlight=sample#sampling-functions

const TOLERANCE = 1e-4;

# ================ module begin ===============================
function create_pol(N::Integer)
"""Create N-dim unit cube."""
    rSTDO = R"""
    library(volesti)
    P <- gen_cube($N,'H')
    P$b <- cbind(matrix(0,1,$N), matrix(1,1,$N))
    """ # mondify cube from [-1 1] to [0 1]
    return rcopy(R"P$A"),rcopy(R"P$b")
end

function create_scaled_pol(N::Integer,U::AbstractArray)
"""Create polytope with increased upper bounds."""
    rSTDO = R"""
    library(volesti)
    P <- gen_cube($N,'H')
    P$b <- cbind(matrix(0,1,$N), matrix(c(unlist($U)),1,$N))
    """ # mondify cube from [-1 1] to [0 1]
    return rcopy(R"P$A"),rcopy(R"P$b")
end

function add_to_pol(Al::AbstractArray,bl::Union{AbstractArray{<:Number},Number})
    R"""
    P$A <- rbind(P$A, $Al)
    P$b <- c(P$b, $bl)"""
    nothing
end

function get_pol()
    A = rcopy(R"P$A")
    b = rcopy(R"P$b")
    return A::Array{Float64,2}, b::Union{AbstractArray{<:Number},Number}
end

function sample_pol(number_of_samples::Integer=1)
    samples = rcopy(R"sample_points(P, $number_of_samples)")
    return samples::Array{Float64,2}
end

function get_pol_volume()
    vol = rcopy(R"volume(P)")
    return vol::Number
end

function sample_polyhedron(N::Integer,A::AbstractArray=[],b::AbstractArray=[])
"""Get N number of uniformly distributed samples from the convex polytope
    defined by A and b. The polytope is most likely defined in the R session
    already. If A and b are not given just call sample_pol"""
    if ~isempty(A)
        # define polytope by A and b using R call
    end
    samples = transpose(sample_pol(N))
    return samples::AbstractArray
end
# ================ module end ===============================

function tt(v, A)
"""General purpose method that implements NxM matrix multiply by 1xM vector, so
    that each column of the matrix is multiplied by the corresponding value of
    the vector."""
    r = similar(A)
    @inbounds for j = 1:size(A,2)
        @simd for i = 1:size(A,1)
            r[i,j] = v[j] * A[i,j]
        end
    end
    r
end

function get_header_names(variables::AbstractArray)
    header = [] # used as header in x_opt results csv
    push!(header, JuMP.name.(variables))
    header = map((x) -> uppercase(replace(x, r"0_|\[|\]" => "")),header[1]) # prettify
    return [header]::AbstractArray
end

function get_slack_idx(power_model::AbstractPowerModel)
    # Slack type == 3
    bus_idx =  [k for (k,v) in power_model.data["bus"] if v["bus_type"] == 3]
    bus_idx = parse(Int64,bus_idx[1])
    gen_idx = [k for (k,v) in power_model.data["gen"] if v["gen_bus"] == bus_idx]
    return parse(Int64, gen_idx[1])
end

function run_qc_relax(pm::AbstractPowerModel, number_of_iterations::Integer, vol::Integer=0)
    """ Iteratively solve modified QC-AC-OPF
    inputs: power model and number of iterations"""
    # Initialize variables
    start = time();
    info(logger,"QC relaxation tolerance set to $TOLERANCE .")
    optimal_setpoints = []
    vars = []
    # build optimization problem
    # get variables: PG\{slack} and Vm (all generators)
    # slack: bus_type == 3
    #
    # to get all: vars = [pm.var[:nw][0][:pg].data; pm.var[:nw][0][:vm].data];
    #
    slack_gen_idx = get_slack_idx(pm)
    for i=1:length(pm.var[:nw][0][:pg].data)
        if i != slack_gen_idx
            push!(vars, JuMP.variable_by_name(pm.model, string("0_pg[",i,"]")))
        end
    end
    gen_indexes = map(x -> x["gen_bus"], values(pm.data["gen"]))
    for g in gen_indexes
        push!(vars, JuMP.variable_by_name(pm.model, string("0_vm[",g,"]")))
    end
    # extract header names
    header = get_header_names(vars) # call before normalization

    N = length(vars);
#TODO: Currently the variables of the optimization problem are scaled between
# [0;1] so we sample from a unit cube. if create_scaled_pol() is used the polytope
# polytope gets scaled and we can remove the scaling of the vars. This later
# might prove to be the proper solution as I could not make sure that scaling
# only vars is sufficient.
    nFactor = JuMP.upper_bound.(vars)
    # auxiliary variables to bridge AffExpr <-> NLconstraint
    #   See "Syntax notes" http://www.juliaopt.org/JuMP.jl/dev/nlp/#Nonlinear-Modeling-1)
    @variable(pm.model,aux[1:N])
    @constraint(pm.model,aux .== vars./nFactor)
    vars = vars./nFactor; # normalize after aux is defined

    # Create x_hat random sample
    create_pol(N) # Create initial unit cube
    x_hat = sample_pol() # Initial sample from N unit cube

    # Build optimization problem
    # Alter the QCCR problem
    @variable(pm.model, r);
    @objective(pm.model, Min, r);
    @NLparameter(pm.model, x_hat_p[i = 1:N] == x_hat[i]) # updating the parameter has performance gain vs rebuilding the model
    @NLconstraint(pm.model, con_sphere, sqrt(sum((aux[i]-x_hat_p[i])^2 for i in 1:N))<= r);

    debug(logger,"number of it: $number_of_iterations")
    for k = 1:number_of_iterations # Julia for loop has it's own scope, x_hat not accessible
        # Execute optimization
        result = optimize_model!(pm, optimizer=optimizer_with_attributes(Ipopt.Optimizer, "print_level" => 0))
        # Calculate outputs
        x_opt = JuMP.value.(vars);
        n_normT = transpose(x_opt[:] - x_hat); # This can be simply -> objective_value(pm.model) NO, in the article that is not abs() but l^1 norm
        debug(logger,string("Objective_value: " ,JuMP.value(r)))
        push!(optimal_setpoints, (x_opt[:] .* nFactor)' )
        if !(isapprox(JuMP.value(r), 0; atol=TOLERANCE))
            # Update results
            debug(logger,string("new A row: ", n_normT))
            debug(logger,string("new b val: ", n_normT*x_opt))
            add_to_pol(n_normT, n_normT*x_opt)
            debug(logger,string("Polytope values |get_pol() \n", get_pol() ))
        end
        try
            x_hat = sample_pol()
            debug(logger,"New points sampled.")
        catch e
            debug(logger,typeof(e))
            debug(logger,e)
            return A,b
        end
        if vol > 0
            if k % vol == 0 || k == 1
                # save file
                A,b = get_pol()
                if k == 1
                    writedlm(string("polytope_A_",0,".csv"), A, ',')
                    writedlm(string("polytope_b_",0,".csv"), b, ',')
                else
                    writedlm(string("polytope_A_",k,".csv"), A, ',')
                    writedlm(string("polytope_b_",k,".csv"), b, ',')
                end
            end
        end
        # print("Volume: ", get_volume(N,A,b))
        JuMP.set_value.(x_hat_p, x_hat)
        debug(logger,"New x_hat JuMP parameter values set.")
    end
    A,b = get_pol()
    debug(logger,"Return polytope.")
    debug(logger, string("Variable Names: ", header))
    return A,b,optimal_setpoints,header,nFactor
end

function mod_acopf!(pm::AbstractPowerModel)
"""Modify the objective function and add NL constraint to AC-OPF problem."""
    vars = []
    slack_gen_idx = get_slack_idx(pm)
    for i=1:length(pm.var[:nw][0][:pg].data)
        if i != slack_gen_idx
            push!(vars, JuMP.variable_by_name(pm.model, string("0_pg[",i,"]")))
        end
    end
    gen_indexes = map(x -> x["gen_bus"], values(pm.data["gen"]))
    for g in gen_indexes
        push!(vars, JuMP.variable_by_name(pm.model, string("0_vm[",g,"]")))
    end
    pm.model[:vars] = vars;
    N = length(vars);
    x_hat = zeros(1,N);

    @variable(pm.model, r);
    @objective(pm.model, Min, r);
    @NLparameter(pm.model, x_hat_p[i = 1:N] == x_hat[i]) # current row of samples
    pm.model[:x_hat_p] = x_hat_p;

    @NLconstraint(pm.model, con_sphere,
        sqrt(sum((vars[i]-x_hat_p[i])^2 for i in 1:N))<= r);

    return nothing
end

function run_mod_acopf(pm::AbstractPowerModel,sp::AbstractArray)
    """Run modified ACOPF with setpoints defined by sp"""
    power_model = pm; # do not modify original model
    JuMP.set_value.(power_model.model[:x_hat_p],sp)
    result = optimize_model!(power_model, optimizer=optimizer_with_attributes(Ipopt.Optimizer, "print_level" => 0))

    return JuMP.value.(power_model.model[:vars])
end

max_stars_repo_name: timonviola/PSDC

max_issues_repo_name: timonviola/PSDC

max_forks_repo_name: timonviola/PSDC


text:
using OptEvapCool
using StatsPlots
using Printf: @sprintf
using LaTeXStrings
using Plots.PlotMeasures

function anu_crossbeam_trap(duration = 1.97, input_dir = "")
    # Physical parameters
    Np = 3e7 # Initial atom count
    T₀ = 15e-6 # Initial temperature
    species = Rb87

    # Numerical parameters
    Nt = ceil(Int64, 1e5) # Test particles
    F = Np / Nt
    Nc = 3 # Target number of test particles per cell

    # Beam parameters
    P₁ = exponential_ramp(15, 1.5, 0.8) # Watts
    P₂ = exponential_ramp(7.5, 1.5, 0.8)
    # Set final power to 1.5W to achieve BEC

    w₀ = 130e-6 # Waist (m)
    θ = ( 22.5 * π / 180 ) / 2 # Half-angle between beams

    dir1 = [sin(θ), 0, cos(θ)] # Directions
    dir2 = [-sin(θ), 0, cos(θ)]
    λ₁ = 1064e-9 # Wavelengths
    λ₂ = 1090e-9

    beam1 = GaussianBeam([0,0,0], dir1, P₁, w₀, λ₁)
    beam2 = GaussianBeam([0,0,0], dir2, P₂, w₀, λ₂)

    # Trapping frequencies
    m = species.m
    κ = kappa(species)

    Uₜ_coeff = 2 * κ / (π * w₀^2)
    Uₜ(t) = Uₜ_coeff * (P₁(t) + P₂(t)) #Trap depth

    ωx_coeff = sqrt(4 * cos(θ)^2 / (m * w₀^2))
    ωz_coeff = sqrt(4 * sin(θ)^2 / (m * w₀^2))
    ωy_coeff = sqrt(4 / (m * w₀^2))

    ωx(t) = ωx_coeff * sqrt(Uₜ(t))
    ωz(t) = ωz_coeff * sqrt(Uₜ(t))
    ωy(t) = ωy_coeff * sqrt(Uₜ(t))

    # Cloud initialisation
    positions = harmonic_boltzmann_positions(Nt, m, T₀, ωx(0), ωy(0), ωz(0))
    velocities = boltzmann_velocities(Nt, m, T₀)

    # Adjust vertical equilibrium position
    function y0(t)
        q = -4*g^2 / (ωy(t)^4 * w₀^2)
        return - w₀/2 * sqrt(- q * exp(q))
    end

    centre(t) = [0, y0(t), 0]
    positions .+= centre(0)

    # Three-body loss & background loss
    K = 4.3e-29 * 1e-12
    τ_bg = 180

    # Maximum timestep
    max_dt(t) = 0.05 * 2π / max(ωx(t), ωy(t), ωz(t))

    # Gravitaty
    acc_grav! = acceleration(gravity)
    pot_grav! = potential(gravity)

    # GAUSSIAN BEAM SIMULATION
    acc_b1! = acceleration(beam1)
    acc_b2! = acceleration(beam2)
    pot_b1! = potential(beam1)
    pot_b2! = potential(beam2)

    function accel_gb(p, s, t, o)
        fill!(o, 0.0)
        acc_b1!(p, s, t, o)
        acc_b2!(p, s, t, o)
        acc_grav!(p, s, t, o)
        return nothing
    end

    function poten_gb(p, s, t, o)
        fill!(o, 0.0)
        pot_b1!(p, s, t, o)
        pot_b2!(p, s, t, o)
        pot_grav!(p, s, t, o)
        return nothing
    end

    sensor_gb = GlobalSensor()
    measure_gb = measurer(sensor_gb, 0.001, ωx, ωy, ωz, centre)

    T_gb(t) = (length(sensor_gb.ke) > 0) ? 2 * last(sensor_gb.ke) / (3 * kB) : T₀
    evap_gb = ellipsoid_evap(ωx, ωy, ωz, T_gb, 1e-6, centre)

    gb_conditions = SimulationConditions(species, F, positions, velocities,
        accel_gb, poten_gb, evap = evap_gb, τ_bg = τ_bg, K = K)

    if input_dir == "" # Simulate new results
        @info "Gaussian beam simulation"
        final_gb_cloud = evolve(gb_conditions, duration;
            Nc = Nc, max_dt = max_dt, measure = measure_gb)
    end

    # HARMONIC SIMULATION
    harm = HarmonicField(ωx, ωy, ωz, centre)
    acc_harm! = acceleration(harm)
    pot_harm! = potential(harm)

    function accel_harm(p, s, t, o)
        fill!(o, 0.0)
        acc_harm!(p, s, t, o)
        acc_grav!(p, s, t, o)
        return nothing
    end

    function poten_harm(p, s, t, o)
        fill!(o, 0.0)
        pot_harm!(p, s, t, o)
        pot_grav!(p, s, t, o)
        return nothing
    end

    function pot_evap_harm(p, s, t, o)
        fill!(o, 0.0)
        pot_harm!(p, s, t, o)
        return nothing
    end

    evap_harm = energy_evap(Uₜ, pot_evap_harm)

    sensor_harm = GlobalSensor()
    measure_harm = measurer(sensor_harm, 0.001, ωx, ωy, ωz, centre)

    harm_conditions = SimulationConditions(species, F, positions, velocities,
        accel_harm, poten_harm, evap = evap_harm, τ_bg = τ_bg, K = K)

    if input_dir == ""
        @info "Harmonic approximation simulation"
        final_harm_cloud = evolve(harm_conditions, duration;
            Nc = Nc, max_dt = max_dt, measure = measure_harm)
    end

    if input_dir != "" # Load results
        sensor_gb = loadsensor("$input_dir/sensor-gb-data.csv")
        sensor_harm = loadsensor("$input_dir/sensor-harm-data.csv")
    end

    # Save sensor data (before plotting, in case plotting fails)
    ft = filetime()
    dir = "./results/$ft-anu-crossbeam-trap"
    mkpath(dir)

    savecsv(sensor_gb, "$dir/sensor-gb-data.csv")
    savecsv(sensor_harm, "$dir/sensor-harm-data.csv")

    # THEORY
    # Pethick model
    peth_t, peth_N, peth_T, peth_Γ = harmonic_theory(
        species, ωx, ωy, ωz, T₀, Np, Uₜ, duration, τ_bg
    )
    # Purdue model
    pur_t, pur_N, pur_T, pur_Γ = purdue_theory(
        species, ωx, ωy, ωz, T₀, Np, Uₜ, duration, τ_bg, K
    )

    # PLOTTING
    default(fontfamily="Computer Modern",
        linewidth=2, framestyle=:box, label=nothing, grid=false)
    scalefontsizes()
    scalefontsizes(1.5)

    window_time = min(1e-2, duration/2)
    gb_time = sensor_gb.time
    harm_time = sensor_harm.time
    gb_window_size = OptEvapCool.window_time_size(gb_time, window_time)
    harm_window_size = OptEvapCool.window_time_size(harm_time, window_time)
    rolling_gb_time = rollmean(gb_time, gb_window_size)
    rolling_harm_time = rollmean(harm_time, harm_window_size)

    # Colour palette
    col1 = RGB(0.0588, 0, 0.5882)
    col2 = RGB(0.9000, 0, 0.3490)
    col3 = RGB(1, 0.6510, 0)
    col4 = RGB(0.0824, 0.8392, 0)

    # Temperature
    rolling_gb_temp = 2 / (3 * kB) * rollmean(sensor_gb.ke, gb_window_size)
    rolling_harm_temp = 2 / (3 * kB) * rollmean(sensor_harm.ke, harm_window_size)

    temp_order = -6
    temp_yformatter(y) = @sprintf("%.1f",y/(10.0^temp_order))

    temperature_plt = plot(rolling_gb_time, rolling_gb_temp,
        xlabel = "Time (s)",
        ylabel = L"\textrm{Temperature\ \ }({}\times10^{%$temp_order}\mathrm{K})",
        label = false,
        linecolor = col1,
        ls = :solid,
        ylims = (0, Inf),
        minorticks = 5,
        yformatter = temp_yformatter,
        dpi = 300)
    plot!(rolling_harm_time, rolling_harm_temp,
          label = false, linecolor = col2, ls= :dash)
    plot!(peth_t, peth_T, label = false, linecolor = col3, ls = :dot)
    plot!(pur_t, pur_T, label = false, linecolor = col4, ls = :dashdot)

    # Number
    gb_Np = sensor_gb.Nt .* sensor_gb.F
    harm_Np = sensor_harm.Nt .* sensor_harm.F

    num_order = 6
    num_yformatter(y) = @sprintf("%.1f",y/(10.0^num_order))

    number_plt = plot(gb_time, gb_Np,
        xlabel = "Time (s)",
        ylabel = L"\textrm{Number\ \ }({}\times10^{%$num_order})",
        label = false,
        linecolor = col1,
        ls = :solid,
        ylims = (0, Inf),
        minorticks = 5,
        yformatter = num_yformatter,
        dpi = 300)
    plot!(harm_time, harm_Np,
          label = false, linecolor = col2, ls= :dash)
    plot!(peth_t, peth_N, label = false, linecolor = col3, ls = :dot)
    plot!(pur_t, pur_N, label = false, linecolor = col4, ls = :dashdot)

    # Collrate
    rolling_gb_timesteps = (
        gb_time[gb_window_size:end] - gb_time[1:end - gb_window_size + 1]
    )
    rolling_harm_timesteps = (
        harm_time[harm_window_size:end] - harm_time[1:end - harm_window_size + 1]
    )
    rolling_gb_Nt = rollmean(sensor_gb.Nt, gb_window_size)
    rolling_harm_Nt = rollmean(sensor_harm.Nt, harm_window_size)
    rolling_gb_Γ = (2 * rolling(sum, sensor_gb.coll, gb_window_size) ./
        (rolling_gb_Nt .* rolling_gb_timesteps)
    )
    rolling_harm_Γ = (2 * rolling(sum, sensor_harm.coll, harm_window_size) ./
        (rolling_harm_Nt .* rolling_harm_timesteps)
    )

    collrate_plt = plot(rolling_gb_time, rolling_gb_Γ,
        xlabel = "Time (s)",
        ylabel = "Collision rate (Hz)",
        label = false,
        linecolor = col1,
        ls = :solid,
        dpi = 300)
    plot!(rolling_harm_time, rolling_harm_Γ,
          label = false, linecolor = col2, ls= :dash)
    plot!(peth_t, peth_Γ, label = false, linecolor = col3, ls = :dot)
    plot!(pur_t, pur_Γ, label = false, linecolor = col4, ls = :dashdot)

    # Phase space density
    rolling_gb_n0 = rollmean(sensor_gb.n0 .* sensor_gb.F, gb_window_size)
    rolling_harm_n0 = rollmean(sensor_harm.n0 .* sensor_harm.F, harm_window_size)

    rolling_gb_psd = rolling_gb_n0 .* (
        (2*π*h̄^2 ./ (m * kB .* rolling_gb_temp)) .^1.5
    )
    rolling_harm_psd = rolling_harm_n0 .* (
        (2*π*h̄^2 ./ (m * kB .* rolling_harm_temp)) .^1.5
    )

    psd_plt = plot(rolling_gb_time, rolling_gb_psd,
        xlabel = "Time (s)",
        ylabel = "Phase space density",
        label = false,
        linecolor = col1,
        ls = :solid,
        ylims = (0, Inf),
        minorticks = 5,
        dpi = 300)
    plot!(rolling_harm_time, rolling_harm_psd,
          label = false, linecolor = col2, ls= :dash)

    # Number vs temp.
    rolling_gb_Np = rollmean(gb_Np, gb_window_size)
    rolling_harm_Np = rollmean(harm_Np, harm_window_size)

    numtemp_plt = plot(rolling_gb_temp, rolling_gb_Np,
        xlabel = L"\textrm{Temperature\ }(\mu K)",
        ylabel = L"\textrm{Number\ \ }({}\times10^{%$num_order})",
        yformatter = num_yformatter,
        xformatter = temp_yformatter,
        label = false,
        linecolor = col1,
        ls = :solid,
        ylims = (0, Inf),
        xlims = (0, Inf),
        minorticks = 5,
        right_margin = 5mm,
        dpi = 300)
    #=
    plot!(rolling_harm_temp, rolling_harm_Np,
          label = false, linecolor = col2, ls= :dash)
    plot!(peth_T, peth_N, label = false, linecolor = col3, ls = :dot)
    plot!(pur_T, pur_N, label = false, linecolor = col4, ls = :dashdot)
    =#
    # Experimental Results
    scatter!([700e-9], [6e6], color = "black", label = false,
        markersize = 8, markershape = :utriangle)
    #vline!([2.6], line = (:black, 5))

    # Number vs PSD.
    numpsd_plt = plot(rolling_gb_psd, rolling_gb_Np,
        xlabel = "Phase space density",
        ylabel = L"\textrm{Number\ \ }({}\times10^{%$num_order})",
        yformatter = num_yformatter,
        label = false,
        linecolor = col1,
        ls = :solid,
        ylims = (0, Inf),
        xlims = (0, Inf),
        minorticks = 5,
        right_margin = 5mm,
        dpi = 300)
    # Experimental Results
    expt_T = 700e-9
    expt_N = 6e6
    expt_n0 = expt_N * ωx(1.97) * ωy(1.97) * ωz(1.97) * (m / (2π*kB*expt_T))^1.5
    expt_λdB =sqrt(2*π*h̄^2 / (m * kB * expt_T))
    expt_psd = expt_n0 * expt_λdB^3
    scatter!([expt_psd], [expt_N], color = "black", label = false,
        markersize = 8, markershape = :utriangle)
    #vline!([2.6], line = (:black, 5))

    # Save plots
    savefig(temperature_plt, "$dir/temp.png")
    savefig(collrate_plt, "$dir/collrate.png")
    savefig(number_plt, "$dir/number.png")
    savefig(psd_plt, "$dir/psd.png")
    savefig(numtemp_plt, "$dir/numtemp.png")
    savefig(numpsd_plt, "$dir/numpsd.png")

    return nothing
end


text:
#===============================================================================

    Adaptation schemes for adjusting parameter update steps. Contains
        ✔   Adaptation of uniform random walkers by changing ϵ to target the
            acceptance rate
        ✔   Haario-type adaptive schemes that adjusts the correlation of
            a Gaussian random walker to match that of the empirical covariance
        ✗   Adaptations for MALA algorithms, google search gives
            https://link.springer.com/article/10.1007/s11009-006-8550-0 (compare
            to STAN, decide accordingly)
        ✗   Adaptations for Hamiltonian Monte Carlo, google search gives
            http://proceedings.mlr.press/v28/wang13e.pdf (compare to STAN,
            decide accordingly)

===============================================================================#

#===============================================================================
                Adaptations for the Uniform Random Walker
===============================================================================#

"""
    NoAdaptation

A struct-flag for indicating that no adaptation is to be done.
"""
struct NoAdaptation <: Adaptation end

"""
    mutable struct AdaptationUnifRW{T} <: Adaptation
        proposed::Int64
        accepted::Int64
        target_accpt_rate::Float64
        adapt_every_k_steps::Int64
        scale::T
        min::T
        max::T
        offset::T
        N::Int64
    end

A struct containing information about the way in which to adapt the steps of the
uniform random walker. `proposed` and `accepted` are the internal counters that
keep track of the number of proposed and accepted samples. `target_accpt_rate`
is the acceptance rate of the Metropolis-Hastings steps that is supposed to be
targetted. `min` is used to enforce the minimal allowable range that the random
walker can sample from, `max` enforces the maximum. `offset` introduces some
delay in the start of decelerting the adaptation extent and `scale` is a scaling
parameter for adaptation speed. `N` is the number of coordinates of the random
walker.
"""
mutable struct AdaptationUnifRW{T} <: Adaptation
    proposed::Int64
    accepted::Int64
    target_accpt_rate::Float64
    adapt_every_k_steps::Int64
    scale::T
    min::T
    max::T
    offset::T
    N::Int64

    function AdaptationUnifRW{T}(
            target_accpt_rate, adapt_every_k_steps, scale::T,
            min, max, offset, N,
        ) where T
        new{T}(
            0, 0, target_accpt_rate, adapt_every_k_steps, scale,
            min, max, offset, N
        )
    end
end

"""
    AdaptationUnifRW(θ::K; kwargs...) where K

Define an adaptation scheme for a random walker that does updates on the
parameter of shape and type `θ`. The following named parameters can be specified

# Arguments
- `adapt_every_k_steps=100`: number of steps based on which adaptation happens
- `target_accpt_rate=0.234`: acceptance rate of MH step that is to be targetted
- `scale=1.0`: scaling of the adaptation
- `min=1e-12`: minimum allowable (half)-range of the uniform sampler
- `max=1e7`: maximum allowable (half)-range of the unifor sampler
- `offset=1e2`: number of adaptation steps after which the shrinking of adaptive
                steps is supposed to start.
"""
function AdaptationUnifRW(θ::K; kwargs...) where K
    non_vec_kwargs = [:adapt_every_k_steps, :target_accpt_rate]
    trgt_ar = _assure_scalar( get(kwargs, non_vec_kwargs[2], 0.234) )
    steps = _assure_scalar( get(kwargs, non_vec_kwargs[1], 100) )
    possible_vec_kwargs = filter(u->!(u.first in non_vec_kwargs), kwargs)

    ulengths = unique( map(x->length(x), values(possible_vec_kwargs)) )

    lul = length(ulengths)
    @assert lul <= 2

    cat = ( lul == 0 || maximum(ulengths) == 1 ? Val(:scalar) : Val(:nonscalar) )

    _AdaptationUnifRW(
        ismutable(K), cat, Val(length(θ)), trgt_ar, steps;
        possible_vec_kwargs...
    )
end

"""
    _adpt_rw_fill_defaults(f::Function)

Return a vector with set parameters, setting them according to the "recipe" `f`,
falling back on default values if necessary.
"""
function _adpt_rw_fill_defaults(f::Function)
    [
        f(s, v) for (s, v) in
        [(:scale, 1.0), (:min, 1e-12), (:max, 1e7), (:offset, 1e2)]
    ]
end

"""
    _adpt_rw_fill_defaults(; kwargs...)

Return a vector with set scalar parameters, falling back on default values if
necessary.
"""
function _adpt_rw_fill_defaults(; kwargs...)
    _adpt_rw_fill_defaults( (s,v)->first(get(kwargs, s, v)) )
end

"""
    _adpt_rw_fill_defaults(N; kwargs...)

Return a vector with set vector parameters, falling back on default values if
necessary.
"""
function _adpt_rw_fill_defaults(N; kwargs...)
    _adpt_rw_fill_defaults(
        (s,v)->_upgrade_to_vec( get(kwargs, s, repeat([v], N)), N )
    )
end

"""
    _adpt_rw_fill_defaults(N; kwargs...)

Return a vector with set `SVector` parameters, falling back on default values if
necessary.
"""
function _adpt_rw_fill_defaults(n::Val{N}; kwargs...) where N
    z = zero(SVector{N,Float64})
    _adpt_rw_fill_defaults(
        (s,v)->_upgrade_to_svec( get(kwargs, s, z.+v), n )
    )
end

"""
    _AdaptationUnifRW(
            ::Any, ::Val{:scalar}, ::Val{N}, trgt_ar, steps; kwargs...
        ) where N

Internal constructor of AdaptationUnifRW, populating entries with scalars.
"""
function _AdaptationUnifRW(
        ::Any, ::Val{:scalar}, ::Val{N}, trgt_ar, steps; kwargs...
    ) where N
    AdaptationUnifRW{Float64}(
        trgt_ar, steps, _adpt_rw_fill_defaults(; kwargs...)..., N
    )
end

"""
    _AdaptationUnifRW(
            ::Val{true}, ::Val{:nonscalar}, ::Val{N}, trgt_ar, steps; kwargs...
        ) where N

Internal constructor of AdaptationUnifRW, populating entries with vectors.
"""
function _AdaptationUnifRW(
        ::Val{true}, ::Val{:nonscalar}, ::Val{N}, trgt_ar, steps; kwargs...
    ) where N
    AdaptationUnifRW{Vector{Float64}}(
        trgt_ar, steps, _adpt_rw_fill_defaults(N; kwargs...)..., N
    )

end

"""
    _AdaptationUnifRW(
            ::Val{false}, ::Val{:nonscalar}, v::Val{N}, trgt_ar, steps; kwargs...
        ) where N

Internal constructor of AdaptationUnifRW, populating entries with `SVector`s.
"""
function _AdaptationUnifRW(
        ::Val{false}, ::Val{:nonscalar}, v::Val{N}, trgt_ar, steps; kwargs...
    ) where N
    AdaptationUnifRW{SVector{N,Float64}}(
        trgt_ar, steps, _adpt_rw_fill_defaults(v; kwargs...)..., N
    )
end

"""
    Base.:(==)(a::AdaptationUnifRW{T}, b::AdaptationUnifRW{S}) where{T,S}

Convenience comparison of `AdaptationUnifRW` instances. Used only for tests.
"""
function Base.:(==)(a::AdaptationUnifRW{T}, b::AdaptationUnifRW{S}) where{T,S}
    T != S && return false
    fns = fieldnames(AdaptationUnifRW)
    for fn in fns
        getfield(a, fn) != getfield(b, fn) && return false
    end
    true
end

"""
    isequal_except(
        ::AdaptationUnifRW{T}, b::AdaptationUnifRW{S}, args...
        ) where{T,S}

Convenience comparison of `AdaptationUnifRW` instances. Used only for tests.
The fields specified by `args` are excluded from comparison.
"""
function isequal_except(
        a::AdaptationUnifRW{T},
        b::AdaptationUnifRW{S},
        args...
    ) where{T,S}
    T != S && return false
    fns = fieldnames(AdaptationUnifRW)
    for fn in fns
        fn in args && (continue)
        getfield(a, fn) != getfield(b, fn) && return false
    end
    true
end

"""
    acceptance_rate(adpt::AdaptationUnifRW)

Compute current acceptance rate of the Metropolis-Hastings update step
"""
acceptance_rate(adpt::AdaptationUnifRW) = (
    (adpt.proposed == 0) ? 0.0 : adpt.accepted/adpt.proposed
)

"""
    acceptance_rate!(adpt::AdaptationUnifRW)

Destructive computation of a current acceptance rate that also resets the
number of proposals and accepted samples to zeros.
"""
function acceptance_rate!(adpt::AdaptationUnifRW)
    a_r = acceptance_rate(adpt)
    reset!(adpt)
    a_r
end

"""
    reset!(adpt::AdaptationUnifRW)

Reset the number of proposals and accepted samples to zero.
"""
function DataStructures.reset!(adpt::AdaptationUnifRW)
    adpt.proposed = 0
    adpt.accepted = 0
end

"""
    readjust!(rw::UniformRandomWalk, adpt::AdaptationUnifRW, mcmc_iter)

Adaptive readjustment of the range for sampling uniforms by the random walker.
"""
function readjust!(rw::UniformRandomWalk, adpt, mcmc_iter)
    δ = compute_δ(adpt, mcmc_iter)
    a_r = acceptance_rate!(adpt)
    ϵ = compute_ϵ(rw.ϵ, adpt, a_r, δ)
    rw.ϵ = ϵ
    ϵ
end

#=
function register!(updt, adpt::AdaptationUnifRW, accepted::AbstractArray{<:Bool}, v)
    register!(updt, adpt, first(accepted), v)
end
=#

"""
    register!(adpt::AdaptationUnifRW, accepted::Bool, ::Any)

Register the result of the acceptance decision in the Metropolis-Hastings step.
"""
function register!(updt, adpt::AdaptationUnifRW, accepted::Bool, ::Any)
    adpt.accepted += accepted
    adpt.proposed += 1
end

"""
    time_to_update(adpt::AdaptationUnifRW)

Return true if it's the time to update the ϵ parameter of Uniform random walks
"""
function time_to_update(::Val{true}, adpt::AdaptationUnifRW)
    adpt.proposed >= adpt.adapt_every_k_steps
end

"""
    compute_δ(p, mcmc_iter)

δ decreases roughly proportional to scale/sqrt(iteration)
"""
function compute_δ(p, mcmc_iter)
    p.scale/sqrt(
        max(
            1.0,
            mcmc_iter/p.adapt_every_k_steps-p.offset
        )
    )
end

"""
    compute_ϵ(ϵ_old, p, a_r, δ, flip=1.0, f=identity, finv=identity)

ϵ is moved by δ to adapt to target acceptance rate
"""
function compute_ϵ(ϵ_old, p, a_r, δ, flip=1.0, f=identity, finv=identity)
    ϵ = finv.(f.(ϵ_old) .+ flip.*(2 .*(a_r > p.target_accpt_rate)-1).*δ)
    ϵ = max.(min.(ϵ,  p.max), p.min)    # trim excessive updates
end

#===============================================================================
                        Haario-type adaptive scheme
===============================================================================#
@doc raw"""
    mutable struct HaarioTypeAdaptation{T,TF} <: Adaptation
        mean::Vector{T}
        cov::Matrix{T}
        adapt_every_k_steps::Int64
        scale::Float64
        N::Int64
        M::Int64
        fλ::TF
    end

A struct containing information about the way in which to adapt the steps of the
Gaussian random walker in a scheme:
```math
X_{i+1} = λZ+(1-λ)W,\quad\mbox{where}\quad Z ∼ N(X_i, c\texttt{Id}),
```
and
```math
W ∼ N(X_i, Σ_A),
```
with $Σ_A$ the covariance matrix that the adaptive scheme aims to learn. `mean`
and `cov` are the empirical mean and covariance of the (appropriately
log-transformed) local `state` vector. Every `adapt_every_k_steps` number of
steps an adaptation is performed. `scale` is a scaling for going from the
empirical covariance to the covariance of the Gaussian random walker. `N` is the
total number of terms based on which the `mean` and `cov` were computed, `M` is
the number of updates since the last call to adaptation and `fλ` is a function
for determining the `λ` weight.

    HaarioTypeAdaptation(
            state::Vector{T};
            adapt_every_k_steps=100,
            scale = 2.38^2,
            f::TF=((x,y,z)->x),
        ) where {T,TF}

Base constructor.
"""
mutable struct HaarioTypeAdaptation{T,TF} <: Adaptation
    mean::Vector{T}
    cov::Matrix{T}
    adapt_every_k_steps::Int64
    scale::Float64
    N::Int64
    M::Int64
    fλ::TF

    function HaarioTypeAdaptation(
            state::Vector{T};
            adapt_every_k_steps=100,
            scale = 2.38^2,
            f::TF=((x,y,z)->x),
        ) where {T,TF}
        new{T,TF}(
            zero(state),
            zeros(T, (length(state), length(state))),
            adapt_every_k_steps,
            scale,
            1,
            0,
            f
        )
    end
end

register_only_on_my_turn(::Val{false}, ::HaarioTypeAdaptation) = false

function register_only_on_my_turn(::Val{true}, updt::HaarioTypeAdaptation)
    updt.M += 1 # register an update on my turn
    false
end

function register!(updt, adpt::HaarioTypeAdaptation, ::Any, θ)
    remove_constraints!(updt.rw.gsn_A, θ)
    old_sum_sq = (adpt.N-1)/adpt.N * adpt.cov + adpt.mean * adpt.mean'
    adpt.mean .= adpt.mean .* (adpt.N/(adpt.N+1)) .+ θ ./ (adpt.N+1)
    new_sum_sq = old_sum_sq + (θ * θ')/adpt.N
    adpt.cov .= new_sum_sq - (adpt.N+1)/adpt.N*(adpt.mean * adpt.mean')
    reimpose_constraints!(updt.rw.gsn_A, θ)
    adpt.N += 1
end

function time_to_update(::Val{true}, adpt::HaarioTypeAdaptation)
    ttu = (adpt.M >= adpt.adapt_every_k_steps)
    ttu && (adpt.M = 0) # reset the counter
    ttu
end

function readjust!(rw::GaussianRandomWalkMix, adpt, mcmc_iter)
    Σ = 2.38^2/length(rw)*adpt.cov
    rw.gsn_B.Σ = Symmetric(Σ)
    rw.λ = adpt.fλ(rw.λ, adpt.N, mcmc_iter)
end


text:
# Description: miscellaneous functions

using ArrayViews

export rmfile, printbacktrace, smallmatvec!, smallmatvec, smallmatTvec!,
        smallmatvec_kernel!, smallmatTvec_kernel!,
        smallmatTvec, smallmatmat!,
        smallmatmat, smallmatmatT!, smallmatmatT, smallmatTmat!, smallmatTmat,
        smallmatvec_revv!,
        smallmatmat_kernel!, smallmatmatT_kernel!, smallmatTmat_kernel!,
        checkZeroRows,
        checkZeroColumns, checkIdenticalColumns, checkSparseColumns,
        checkSparseRows, findLarge, isSymmetric, make_symmetric!,
        getBranchName, getTimeString, isFieldDefined, get_parallel_fname,
        hasExtension, joinpath_ascii

@doc """
 ### Tools rmfile

 Removes the file if it exists.  This is only necessary because it is an error
   to remove a file that does not exist.  Does not work on directories

   Arguments:
     fname::AbstractString  :  name of file

"""->
function rmfile(fname::AbstractString; recursive=false)
  if isfile(fname)
    rm(fname, recursive=recursive)
  end
  return nothing
end

# this function is both slow and dangerous
function readdlmc(fname)
  map(x->eval(parse(x)),readcsv(fname))
end

function printbacktrace(file=STDOUT)
# print the current call stack
  bt = backtrace()
  s = sprint(io->Base.show_backtrace(io, bt))
  println(file,"backtrace: ", s)

  return s
end

#------------------------------------------------------------------------------
# Ax = b

"""
  Computes the product Ax = b

  This is a wrapper around different matrix-vector multiplication routines.
  It calls a special method for very small matrices and vectors, falling back
  to BLAS otherwise.  Performance tests indicate the check for which function
  to call is not that expensive

  **Inputs**

   * A: the matrix
   * x: the vector

  **Inputs/Outputs**

   * b: the output vector

  Aliasing restrictions: same as BLAS, b cannot alias anything
"""
function smallmatvec!(A::AbstractArray{T,2},
           x::AbstractArray{T2,1}, b::AbstractArray{T3, 1}) where {T, T2, T3}

  m, n = size(A)

  # extrapolating from the square matrix case, where my algorithm
  # is faster for m = n <= 8
  # for mixed float-int operations, the arrays need to be
  # reinterpratable
  if m*n <= 65
    smallmatvec_kernel!(A, x, b)
  else
    A_mul_B!(b, A, x)
#    LinAlg.BLAS.gemv!('N', 1.0, A, x, 0.0, b)
  end

  return b
end

# if either of the arguments that get reinerpreted are a ContiguousView,
# don't call BLAS
function smallmatvec!(A::ContiguousView{Complex128,2},
           x::AbstractArray{Float64,1}, b::ContiguousView{Complex128, 1})
# this is the case where we can't reinterpret as required by the BLAS interface

  smallmatvec_kernel!(A, x, b)

  return b
end

function smallmatvec!(A::ContiguousView{Complex128,2},
           x::AbstractArray{Float64,1}, b::AbstractArray{Complex128, 1})
# this is the case where we can't reinterpret as required by the BLAS interface

  smallmatvec_kernel!(A, x, b)

  return b
end

function smallmatvec!(A::AbstractArray{Complex128,2},
           x::AbstractArray{Float64,1}, b::ContiguousView{Complex128, 1})
# this is the case where we can't reinterpret as required by the BLAS interface

  smallmatvec_kernel!(A, x, b)

  return b
end

function smallmatvec!(A::ROContiguousView{Complex128,2, P},
           x::AbstractArray{Float64,1}, b::ContiguousView{Complex128, 1, P2}) where {P<:Array, P2<:Array}
# this is the case where we can't reinterpret as required by the BLAS interface

  smallmatvec_kernel!(A, x, b)

  return b
end

function smallmatvec!(A::ROContiguousView{Complex128,2, P},
           x::AbstractArray{Float64,1}, b::AbstractArray{Complex128, 1}) where P<:Array
# this is the case where we can't reinterpret as required by the BLAS interface

  smallmatvec_kernel!(A, x, b)

  return b
end

"""
  Performs b = alpha*A*x + beta*b.  This works well for small matrices

  **Inputs**

   * A: matrix
   * x: vector
   * alpha: scalar
   * beta: scalar

  **Inputs/Outputs**

   * b: vector
"""
function smallmatvec_kernel!(A::AbstractArray{T,2},
           x::AbstractArray{T2,1}, b::AbstractArray{T3, 1},
           alpha::Number=1, beta::Number=0) where {T, T2, T3}
  (m,n) = size(A)
  xm = length(x)
  bm = length(b)

  @assert n == xm
  @assert m == bm

  # overwrite b, first column only
  @inbounds begin
    @simd for i=1:m
      b[i] = beta*b[i]
    end

    @simd for i=1:n  # loop across columns
       ax_i = alpha*x[i]
      @simd for j=1:m  # loop down columns
        b[j] += A[j,i]*ax_i
      end
    end

  end  # end begin inbounds

  return b

end   # end of smallmatvec! function

# reverse mode to back propigate b to x
function smallmatvec_revv!(A::AbstractArray{T,2},
           x_bar::AbstractArray{T2,1},
           b_bar::AbstractArray{T3, 1}) where {T, T2, T3}

  (m,n) = size(A)
  xm = length(x_bar)
  bm = length(b_bar)

  @assert n == xm
  @assert m == bm

  @inbounds begin
    # reverse mode always updates the output, so no need to special case
    # the i = 1 to overwrite it
    for i=1:n
      @simd for j=1:m
        x_bar[i] += b_bar[j]*A[j, i]
      end
    end

  end  # end begin inbounds

  return x_bar
end

function smallmatvec(A::AbstractArray{T,2}, x::AbstractArray{T2, 1}) where {T, T2}
  (m,n) = size(A)
  T3 = promote_type(T, T2)
  b =zeros(T3, m)
  smallmatvec!(A, x, b)
end

#------------------------------------------------------------------------------
# A.'x = b

"""
  This function does A.'*x = b.  Note this uses the simple transpose of A, not
  the Hermetian transpose.

  This function call BLAS for large matrices and uses a special algorithm for
  small matrices.

  TODO: double check that the cutoff is correct.

  **Inputs**

   * A: the matrix
   * x: the vector

  **Inputs/Outputs**

   * b: vector overwritten with the result
"""
function smallmatTvec!(A::AbstractMatrix, x::AbstractVector, b::AbstractVector)
# currently this doesn't call BLAS for the complex-real case, so its ok to
# pass a ContiguousView to At_mul_B
  m, n = size(A)

  if m*n <= 65
    smallmatTvec_kernel!(A, x, b)
  else
    At_mul_B!(b, A, x)
  end

  return b
end

"""
  Performs b = alpha*A*x + beta*b.  Works well for small matrices
"""
function smallmatTvec_kernel!(A::AbstractMatrix, x::AbstractVector,
                              b::AbstractVector, alpha::Number=1, beta::Number=0)
  (m,n) = size(A)
  xm = length(x)
  bm = length(b)

  @assert m == xm
  @assert n == bm

  @inbounds begin
    for i=1:n  # loop over columns of A
      # perform action for each column of A
      # over write each entry
      b[i] = alpha*A[1, i]*x[1] + beta*b[i]

      # accumulate over rest of column
      @simd for j=2:m
        b[i] += alpha*A[j, i]*x[j]
      end
    end

  end

  return b
end

function smallmatTvec(A::AbstractArray{T, 2}, x::AbstractArray{T2, 1}) where {T, T2}

  (m,n) = size(A)
  T3 = promote_type(T, T2)
  b = zeros(T3, n)
  smallmatTvec!(A, x, b)
end

#------------------------------------------------------------------------------
# matrix-matrix multiply

"""
  Performs a matrix-matrix multiply.  Uses a special algorithm for small
  matrices, otherwise calls BLAS (when possible).

  Julia's `A_mul_B` does not support the alpha and beta arguments of BLAS, so
  its somewhat difficult to make that the interface for `smallmatmat`.
  Use `smallmatmat_kernel` instead

  **Inputs**

   * A: first matrix
   * x: second matrix

  **Inputs/Inputs**

   * x: output matrix, overwritten

   Aliasing restrictions: no aliasing allowed
"""
function smallmatmat!(A::AbstractArray{T, 2},
                      x::AbstractArray{T2, 2},
                      b::AbstractArray{T3, 2}) where {T, T2, T3}
  m, n = size(A)
  if m < 6 && n < 9
    smallmatmat_kernel!(A, x, b)
  else
    A_mul_B!(b, A, x)
  end

  return b
end

# don't call BLAS when reinterpret isn't supported

function smallmatmat!(A::ContiguousView{Complex128, 2},
                             x::AbstractArray{T2, 2},
                             b::ContiguousView{Complex128, 2}) where T2

  smallmatmat_kernel!(A, x, b)

  return b
end

function smallmatmat!(A::AbstractArray{Complex128, 2},
                             x::AbstractArray{T2, 2},
                             b::ContiguousView{Complex128, 2}) where T2

  smallmatmat_kernel!(A, x, b)

  return b
end

function smallmatmat!(A::ContiguousView{Complex128, 2},
                             x::AbstractArray{T2, 2},
                             b::AbstractArray{Complex128, 2}) where T2

  smallmatmat_kernel!(A, x, b)

  return b
end

function smallmatmat!(A::ROContiguousView{Complex128, 2, P},
                   x::AbstractArray{T2, 2},
                   b::ContiguousView{Complex128, 2}) where {T2, P<:Array}

  smallmatmat_kernel!(A, x, b)

  return b
end

function smallmatmat!(A::ROContiguousView{Complex128, 2, P},
                   x::AbstractArray{T2, 2},
                   b::AbstractArray{Complex128, 2}) where {T2, P<:Array}

  smallmatmat_kernel!(A, x, b)

  return b
end

function smallmatmat_kernel!(A::AbstractArray{T, 2},
                      x::AbstractArray{T2, 2},
                      b::AbstractArray{T3, 2},
                      alpha::Number=1, beta::Number=0) where {T, T2, T3}
# TODO: this comment needs to be a @doc
# multiplies matrix A by matrix x, writing the solution to matrix b
# both dimensions of A and the final dimension of x are used for looping
# the array sizes are not checked explicitly
# this uses the same technique as smallmatvec!, simply multiplying A by the columns
# of x repeatedly, without making any copies
# Faster than Julia's built-in mat-mat for matrices without m=n~=28

  (m,n) = size(A)
  (xn, p) = size(x)
  (bm, bn) = size(b)

  @assert n == xn
  @assert m == bm
  @assert p == bn

  @inbounds begin
    for k=1:p  # loop over the column vectors of x
      # overwrite b, first column only
      @simd for i=1:m
        b[i, k] = beta*b[i, k]
        #b[i, k] = x[1, k]*A[i, 1]
      end

      for i=1:n  # loop across columns
        xval = alpha*x[i, k]
        @simd for j=1:m  # loop down columns
          b[j, k] += A[j,i]*xval
        end
      end
    end

  end   # end begin inbounds

  return b

end     # end of smallmatmat! function

function smallmatmat(A::AbstractArray{T,2}, x::AbstractArray{T2, 2}) where {T, T2}
  (m,n) = size(A)
  (xn, p) = size(x)
  T3 = promote_type(T, T2)
  b = zeros(T3, m, p)
  smallmatmat!(A, x, b)
end

#------------------------------------------------------------------------------
# matrix-matrix transposed multiplication

"""
  Performs A*(x.') = b.  Uses a special algorithm for small matrices,
  otherwise calls BLAS

  **Inputs**

   * A: the first matrix
   * x: the matrix to be transposed

  **Inputs/Outputs**

   * b: output matrix, overwritten
"""
function smallmatmatT!(A::AbstractArray{T, 2},
                       x::AbstractArray{T2, 2},
                       b::AbstractArray{T3, 2}) where {T, T2, T3}
  m, n = size(A)
  if m < 7 && n < 10
    smallmatmatT_kernel!(A, x, b)
  else
    A_mul_Bt!(b, A, x)
  end

  return b
end

# don't call BLAS when reinterpret isn't supported

function smallmatmatT!(A::ContiguousView{Complex128, 2},
                             x::AbstractArray{T2, 2},
                             b::ContiguousView{Complex128, 2}) where T2

  smallmatmatT_kernel!(A, x, b)

  return b
end

function smallmatmatT!(A::AbstractArray{Complex128, 2},
                             x::AbstractArray{T2, 2},
                             b::ContiguousView{Complex128, 2}) where T2

  smallmatmatT_kernel!(A, x, b)

  return b
end

function smallmatmatT!(A::ContiguousView{Complex128, 2},
                             x::AbstractArray{T2, 2},
                             b::AbstractArray{Complex128, 2}) where T2

  smallmatmatT_kernel!(A, x, b)

  return b
end

function smallmatmatT!(A::ROContiguousView{Complex128, 2, P},
                   x::AbstractArray{T2, 2},
                   b::ContiguousView{Complex128, 2}) where {T2, P<:Array}

  smallmatmatT_kernel!(A, x, b)

  return b
end

function smallmatmatT!(A::ROContiguousView{Complex128, 2, P},
                   x::AbstractArray{T2, 2},
                   b::AbstractArray{Complex128, 2}) where {T2, P<:Array}

  smallmatmatT_kernel!(A, x, b)

  return b
end

function smallmatmatT_kernel!(A::AbstractArray{T, 2},
                       x::AbstractArray{T2, 2},
                       b::AbstractArray{T3, 2},
                       alpha::Number=1, beta::Number=0) where {T, T2, T3}
# TODO: this comment needs to be a @doc
# multiplies A by x.', storing result in b

  (p, xn) = size(x)
  (m,n) = size(A)
  (bm, bn) = size(b)

  @assert n == xn
  @assert m == bm
  @assert p == bn

  @inbounds begin
    # overwrite b
    for i=1:p
      @simd for j=1:m
        b[j, i] = beta*b[j, i]
      end
    end

    # add to b
    for j=1:n  # loop across remaining columns of A
      for i=1:m  # loop down a column
        a_i = alpha*A[i, j]
        # multiply this entry by row i of x.'
        @simd for k=1:p
          b[i,k] += a_i*x[k, j]
        end
      end
    end

  end   # end begin inbounds

  return b

end     # end of smallmatmatT! function

function smallmatmatT(A::AbstractArray{T,2}, x::AbstractArray{T2, 2}) where {T, T2}
  (m,n) = size(A)
  (p, xn) = size(x)
  T3 = promote_type(T, T2)
  b = zeros(T3, m, p)
  smallmatmatT!(A, x, b)
end

#------------------------------------------------------------------------------
# matrix transposed-matrix multiply

"""
  Performs A.'*x = b.  Uses a special algorithm for small matrices, otherwise
  calls BLAS

  **Inputs**

   * A: the matrix to be transposed
   * x: the second matrix

  **Inputs/Outputs**

   * b: output matrix, overwritten

  Aliasing restrictions: no aliasing
"""
function smallmatTmat!(A::AbstractMatrix{T},
                              x::AbstractMatrix{T2},
                              b::AbstractMatrix{T3}) where {T, T2, T3}
  m, n = size(A)
  if m < 10 && n < 7
    smallmatTmat_kernel!(A, x, b)
  else
    At_mul_B!(b, A, x)
  end

  return b
end

# this doesn't call BLAS in the Complex128-Float64 case, so no need to deal
# with reinterpret issues

#TODO: performance test
function smallmatTmat_kernel!(A::AbstractMatrix{T}, x::AbstractMatrix{T2},
                       b::AbstractMatrix{T3},
                       alpha::Number=1, beta::Number=0) where {T, T2, T3}

  n, m = size(A)
  xn, p = size(x)
  bm, bn = size(b)

  @assert n == xn
  @assert m == bm
  @assert p == bn

  # overwrite b
  @inbounds @simd for i=1:p
    @simd for j=1:m
      b[j, i] = beta*b[j, i]
    end
  end

  @inbounds @simd for i=1:p  # columns of x
    @simd for j=1:m  # columns of A
      # zero b here?
      @simd for k=1:n
        b[j, i] += alpha*A[k, j]*x[k, i]
      end
    end
  end

  return nothing
end

function smallmatTmat(A::AbstractMatrix{T}, x::AbstractMatrix{T2}) where {T, T2}

  n, m = size(A)
  n, p = size(x)
  T3 = promote_type(T, T2)
  b = zeros(m, p)
  smallmatTmat!(A, x, b)
  return b
end

function checkZeroRows(A::AbstractArray{T,2}, tol::AbstractFloat) where T <: Number
# checks each row of a matrix for zeros
# 2d matrices only
# returns the integer number of zero rows, and a Bool
# array telling which rows have all zeros
  (m,n) = size(A)

  zero_mask = zeros(Bool, m)  # record whether each row has all zeros
  zero_row = zeros(Bool, n)  # record results for a row
  for i=1:m
    fill!(zero_row, false)
    for j=1:n
      if abs(A[i,j]) < tol  # if zero entry
        zero_row[j] = true
      end
    end  # end loop over column

    rowsum = sum(zero_row)
    zero_mask[i] = (rowsum == n)  # true if all zeros in row
  end  # end loop over rows

  return sum(zero_mask), zero_mask

end   # end of checkZeroRows function

function checkZeroColumns(A::AbstractArray{T,2}, tol::AbstractFloat) where T <: Number
# TODO: @doc this
# checks each column of a matrix for zeros
# 2d matrices only
# returns integer number of zero rows, and a Bool
# array telling which rows have all zeros
  (m,n) = size(A)

  zero_mask = zeros(Bool, n)  # record whether each column has all zeros
  zero_row = zeros(Bool, m)  # record results for a column
  for i=1:n
    fill!(zero_row, false)
    for j=1:m
      if abs(A[j, i]) < tol  # if zero entry
        zero_row[j] = true
      end
    end  # end loop over row

    rowsum = sum(zero_row)
    zero_mask[i] = (rowsum == n)  # true if all zeros in row
  end  # end loop over columns

  return sum(zero_mask), zero_mask

end   # end checkZeroColumns function

function checkIdenticalColumns(A::AbstractArray{T,2},
                               colnum::Integer,
                               tol::AbstractFloat) where T <: Number
# TODO: @doc this
# checks which columns are identical to column number col
# returns number of column identical and an array of bools telling which ones
# does not count column colnum as being the same as itself

  (m,n) = size(A)
  cnt = 0

  col = sview(A, :, colnum)
  is_same = zeros(Bool, n)

  for i=1:n  # loop over columns

    if i == colnum  # skip the specified column
      continue
    end
    col_i = sview(A, :, i)
    diff_norm = norm(col - col_i)/m

    if diff_norm < tol
      is_same[i] = true
      cnt += 1
    end
  end     # end of for loop over columns

  return cnt, is_same

end     # end checkIdenticalColumns

function findLarge(A::AbstractArray{T,2}, tol::AbstractFloat) where T <: Number
# TODO: @doc this

  (m,n) = size(A)
  cnt = 0
  for i=1:m
    for j=1:n
      if A[i,j] > tol
        println(i, ", ", j)
        cnt += 1
      end
    end
  end

  return cnt

end     # end findLarge

function checkSparseColumns(A::AbstractArray{T,2},
sparsity_bnds::AbstractArray{T2, 2},
tol::AbstractFloat) where {T <: Number, T2 <: Integer}
# TODO: @doc this
# checks that all entries outside the range specified by sparsity_bnds
# are zero
# returns the number of columns with out of bounds entries and an array of
# bools specifying which ones

  (m,n) = size(A)
  out_of_bounds = zeros(Bool, n)

  for i=1:n  # loop over columns
    min = sparsity_bnds[1, i]
    max = sparsity_bnds[2, i]

    for j=1:(min -1)
      entry_j = A[j, i]
      if abs(entry_j) > tol
        out_of_bounds[i] = true
        println("entry ", j, ", ", i, " is non zero")
        break
      end
    end

    for j=(max+1):m
      entry_j = A[j, i]
      if abs(entry_j) > tol
        out_of_bounds[i] = true
        println("entry ", j, ", ", i, " is non zero")
        break
      end
    end

  end # end loop over columns

  cnt = sum(out_of_bounds)

  return cnt, out_of_bounds

end     # end of checkSparseColumns function

function checkSparseRows(A::AbstractArray{T,2},
sparsity_bnds::AbstractArray{T2, 2},
tol::AbstractFloat) where {T <: Number, T2 <: Integer}
# checks that all entries outside the range specified by sparsity_bnds
# are zero
# returns the number of rows with out of bounds entries and an array of
# bools specifying which ones

  (m,n) = size(A)
  out_of_bounds = zeros(Bool, m)

  for i=1:m  # loop over columns
    min = sparsity_bnds[1, i]
    max = sparsity_bnds[2, i]

    for j=1:(min -1)
      entry_j = A[i, j]
      if abs(entry_j) > tol
        out_of_bounds[i] = true
        println("entry ", i, ", ", j, " is non zero")
        break
      end
    end

    for j=(max+1):n
      entry_j = A[i, j]
      if abs(entry_j) > tol
        out_of_bounds[i] = true
        println("entry ", i, ", ", j, " is non zero")
        break
      end
    end

  end     # end loop over columns

  cnt = sum(out_of_bounds)

  return cnt, out_of_bounds

end     # end of checkSparseRows function

@doc """
### ODLCommonTools.isSymmetric

  This function checks if an array is symmetric or not, using
  the specified tolerance for comparing if two entries are
  equal.

  Inputs:
    A : an array to check for symmetry, must be possible to access all
        entries
    tol: tolerance for floating point equality

  Outputs:
    val: a Bool indicating if A is symmetric

"""->
function isSymmetric(A::AbstractArray, tol=1e-14)

  (m,n) = size(A)
  val = true
  for j=1:n
    for i=1:j
      val = val && (abs(A[i,j] - A[j,i]) < tol)
    end
  end

  return val
end

@doc """
ODLCommonTools.make_symmetric!

  This function copies the lower triangle of a matrix to the upper triangle,
  ensuring it is exactly symmetric.  It does not checking to determine if
  the matrix is close so symmetric beforehand

  Inputs/Outputs:
    A : a matrix of any type

"""->
function make_symmetric!(A::AbstractMatrix)
# make the matrix symmetrix, no questions asked
# the lower triangle of the matrix is copied to the upper

  for i=1:size(A,1)
    for j=1:(i-1)
      A[j, i] = A[i,j]
    end
  end

  return nothing
end

#----------------------------------------------------------
export FIFOQueue, front
import Base.push!, Base.pop!, Base.length, Base.isempty,
       Base.resize!, Base.empty!
mutable struct FIFOQueue{T}
  s::Array{T, 1}  # array of values
  tail::Int  # current tail of que (points to most recently inserted elements)
  head::Int  # points to least recently inserted elements (next to pop)
  fac::Float64  # factor by which to expand que when it runs out of space

  function FIFOQueue{T}(; size_hint=1, fac=1.4) where T
    # size_hint = initial size of array
    size_hint = max(size_hint, 1)
    arr = Array{T}(size_hint)
    tail = 0
    head = 1
   return new(arr, tail, head, fac)
  end

end

function push!(que::FIFOQueue{T}, val::T) where T
  # check size
  len = length(que.s)
  if (que.tail + 1) > len
    resize_que!(que)
  end

  que.tail += 1

  # do insertion
  que.s[que.tail] = val

  return nothing
end

function pop!(que::FIFOQueue{T}) where T
  # add check that head >= tail ?
  pos = que.head
  val = que.s[pos]
  que.head += 1

  return val
end

function pop!(que::FIFOQueue{T}, vals::AbstractArray{T, 1}) where T
# remove last n elements, where n = length(vals)

  n = length(vals)
  for i=1:n
    vals[i] = pop!(que)
  end

  return nothing
end

# function front retrieve element without removing
function front(que::FIFOQueue{T}) where T
  return que.s[que.head]
end

# function popN!  # remove multiple elements

function length(que::FIFOQueue)
  return que.tail - que.head + 1
end

# for general collection compatability
endof(que::FIFOQueue) = length(que)

function isempty(que::FIFOQueue)
  return que.head > que.tail
end

function empty!(que::FIFOQueue)
  n = length(que)

  for i=1:n
    pop!(que)
  end

  return nothing
end

function resize!(que::FIFOQueue, new_size)
  # change size of que, if doing so will not remove elements
  if new_size > que.tail
    resize!(que.s, new_size)
  else
    println(STDERR, "Warning: not resizing que")
  end

  return nothing
end

# for resizing the que when it runs out of space
function resize_que!(que::FIFOQueue)

  # decide if shifting the queue or resizing it is better
  if que.head > cld(length(que), 10)
    pos = 1
    for i=que.head:que.tail
      que.s[pos] = que.s[i]
      pos += 1
    end
    que.head = 1
    que.tail = pos - 1
  else  # resize
    new_size = convert(Int, ceil(length(que)*que.fac))
    resize!(que.s, new_size)
  end
end

function getBranchName(dir=pwd())
# get the name of the current branch of the git repo in the specified directory
  nm = readstring(`git rev-parse --abbrev-ref HEAD`)
  return nm[1:end-1]  # remove newline

return end

function getTimeString()
  t = now()
  y = Dates.year(t)
  m = Dates.month(t)
  d = Dates.day(t)
  h = Dates.hour(t)
  minutes = Dates.minute(t)

  return "$y-$m-$d $h:$minutes"
end
@doc """

  Wrapper and generalization of Base.isdefined.
  This function takes an object and one or more symbols and checks if
  the fields of the object with those names are defined, but throws an
  exception if the symbol is not a field of the type.  Returns true
  iff all field names are defined.
"""
function isFieldDefined(obj, req_fieldnames...)

  if length(req_fieldnames) == 0
    throw(ErrorException("must check at least one fieldname"))
  end

  fnames = fieldnames(obj)
  obj_type = typeof(obj)
  alldefined = true
  for i=1:length(req_fieldnames)
    fname_i = Symbol(req_fieldnames[i])  # convert to symbol if possible

    if !(fname_i in fnames)
      throw(ErrorException("fieldname $fname_i is not a field of $obj_type"))
    end

    alldefined = alldefined && isdefined(obj, fname_i)
  end

  return alldefined
end

"""
  Given a file name (including extension), adds an _comm_rank, where comm_rank
  is an MPI communicator rank, to the end of the name, before the extension.
  For example, foo.dat -> foo_0.dat

  Also works if the file name does not have an extension:

  foo -> foo_0.dat

  Inputs:
    fname: original file name, including extension
    comm_rank: communicator rank

  Outputs:
    fname_stub: new file name, including extension
"""
function get_parallel_fname(fname::String, comm_rank)

  # figure out where file extension starts
  len = length(fname)
  sep_loc = 0  # location of separator
  for i=len:-1:1
    if fname[i] == '.'
      sep_loc = i
    end
    if fname[i] == '/'  # linux only
      break
    end
  end

  if sep_loc == 0  # fname does not have an extension
    fname_stub = string(fname, "_", comm_rank)
  else
    fname_stub = fname[1:(sep_loc-1)]  # grab the name before the extension
    fname_stub *= "_$comm_rank"  # add the comm_rank
    fname_stub *= fname[sep_loc:end]
  end

  return fname_stub
end

"""
  Returns the file name and extension as separate strings.
  Ex. output.dat -> output, .dat

  If the file does not have an extension the second string will be empty
  **Input**

   * fname: file name to split

  **Outputs**

   * fstem: file name, without extension
   * fext: file extension (including the period)
"""
function split_fname(fname::AbstractString)

  sep_loc = 0
  for i=length(fname):-1:1
    if fname[i] == '.'
      sep_loc = i
    end
  end

  if sep_loc == 0  # no extension
    fstem = identity(fname)
    fext = ""
  else
    fstem = fname[1:(sep_loc-1)]
    fext = fname[sep_loc:end]
  end

  return fstem, fext
end

"""
  Returns true if the file name as an extesion, false otherwise

  **Inputs**

   * fname: AbstractString

  **Outputs**

   * Bool
"""
function hasExtension(fname::AbstractString)
  fstem, fext = split_fname(fname)

  return length(fext) > 0
end

"""
  Wrapper around joinpath() that always returns an String.
"""
function joinpath_ascii(str::String...)

  return ascii(joinpath(str...))
end

"""
  Prepends the given string to another string as though it is a Linux path
  (ie. using a colon delimiter)

  Inputs:

    path: the existing path, can be an empty string
    new_entry: string to add to the path (should not contain the delimiter)

  Outputs:

    new_path: the results of prepend operation
"""
function prepend_path(path::AbstractString, new_entry::AbstractString)

  if path == ""
    new_path = new_entry
  else
    new_path = string(new_entry, ":", path)
  end

  return new_path
end

"""
  Like [`prepend_path`](@ref), but adds the new entry to the end rather than
  the beginning
"""
function append_path(path::AbstractString, new_entry::AbstractString)

  if path == ""
    new_path = new_entry
  else
    new_path = string(path, ":", new_entry)
  end

  return new_path
end


text:
@testset "accumulate" begin
    @testset "basic" begin
        @system SAccumulate(Controller) begin
            a => 1 ~ track
            b(a) => a + 1 ~ accumulate
        end
        s = instance(SAccumulate)
        @test s.a' == 1 && s.b' == 0
        update!(s)
        @test s.a' == 1 && s.b' == 2
        update!(s)
        @test s.a' == 1 && s.b' == 4
        update!(s)
        @test s.a' == 1 && s.b' == 6
    end

    @testset "cross reference" begin
        @system SAccumulateXRef(Controller) begin
            a(b) => b + 1 ~ accumulate
            b(a) => a + 1 ~ accumulate
        end
        s = instance(SAccumulateXRef)
        @test s.a' == 0 && s.b' == 0
        update!(s)
        @test s.a' == 1 && s.b' == 1
        update!(s)
        @test s.a' == 3 && s.b' == 3
        update!(s)
        @test s.a' == 7 && s.b' == 7
    end

    @testset "cross reference mirror" begin
        @system SAccumulateXrefMirror1(Controller) begin
            a(b) => b + 1 ~ accumulate
            b(a) => a + 2 ~ accumulate
        end
        @system SAccumulateXrefMirror2(Controller) begin
            a(b) => b + 2 ~ accumulate
            b(a) => a + 1 ~ accumulate
        end
        s1 = instance(SAccumulateXrefMirror1); s2 = instance(SAccumulateXrefMirror2)
        @test s1.a' == s2.b' == 0 && s1.b' == s2.a' == 0
        update!(s1); update!(s2)
        @test s1.a' == s2.b' == 1 && s1.b' == s2.a' == 2
        update!(s1); update!(s2)
        @test s1.a' == s2.b' == 4 && s1.b' == s2.a' == 5
        update!(s1); update!(s2)
        @test s1.a' == s2.b' == 10 && s1.b' == s2.a' == 11
    end

    @testset "time" begin
        @system SAccumulateTime(Controller) begin
            t(x=context.clock.time) => 0.5x ~ track(u"hr")
            a => 1 ~ track
            b(a) => a + 1 ~ accumulate
            c(a) => a + 1 ~ accumulate(time=t)
        end
        s = instance(SAccumulateTime)
        @test s.a' == 1 && s.b' == 0 && s.c' == 0
        update!(s)
        @test s.a' == 1 && s.b' == 2 && s.c' == 1
        update!(s)
        @test s.a' == 1 && s.b' == 4 && s.c' == 2
        update!(s)
        @test s.a' == 1 && s.b' == 6 && s.c' == 3
    end

    @testset "unit hour" begin
        @system SAccumulateUnitHour(Controller) begin
            a => 1 ~ accumulate(u"hr")
        end
        s = instance(SAccumulateUnitHour)
        @test iszero(s.a')
        update!(s)
        @test s.a' == 1u"hr"
        update!(s)
        @test s.a' == 2u"hr"
    end

    @testset "unit day" begin
        @system SAccumulateUnitDay(Controller) begin
            a => 1 ~ accumulate(u"d")
        end
        s = instance(SAccumulateUnitDay)
        @test iszero(s.a')
        update!(s)
        @test s.a' == 1u"hr"
        update!(s)
        @test s.a' == 2u"hr"
    end

    @testset "when" begin
        @system SAccumulateWhen(Controller) begin
            t(context.clock.tick) ~ track::int
            f ~ preserve(parameter)
            w(t, f) => t < f ~ flag
            a => 1 ~ accumulate
            b => 1 ~ accumulate(when=w)
            c => 1 ~ accumulate(when=!w)
        end
        n = 5
        s = instance(SAccumulateWhen; config=:0 => :f => n)
        simulate!(s, stop=n)
        @test s.a' == n
        @test s.b' == n
        @test s.c' == 0
        simulate!(s, stop=n)
        @test s.a' == 2n
        @test s.b' == n
        @test s.c' == n
    end

    @testset "minmax" begin
        @system SAccumulateMinMax(Controller) begin
            a => 1 ~ accumulate(min=2, max=5)
        end
        r = simulate(SAccumulateMinMax, stop=4)
        @test r[!, :a] == [2, 3, 4, 5, 5]
    end

    @testset "reset" begin
        @system SAccumulateReset(Controller) begin
            t(context.clock.tick) ~ track::int
            r(t) => t % 4 == 0 ~ flag
            a => 1 ~ accumulate
            b => 1 ~ accumulate(reset=r)
            c => 1 ~ accumulate(reset=r, init=10)
        end
        r = simulate(SAccumulateReset, stop=10)
        @test r[!, :a] == [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
        @test r[!, :b] == [0, 1, 2, 3, 4, 1, 2, 3, 4, 1, 2]
        @test r[!, :c] == [10, 11, 12, 13, 14, 11, 12, 13, 14, 11, 12]
    end

    @testset "transport" begin
        @system SAccumulateTransport(Controller) begin
            a(a, b) => -max(a - b, 0) ~ accumulate(init=10)
            b(a, b, c) => max(a - b, 0) - max(b - c, 0) ~ accumulate
            c(b, c) => max(b - c, 0) ~ accumulate
        end
        s = instance(SAccumulateTransport)
        @test s.a' == 10 && s.b' == 0 && s.c' == 0
        update!(s)
        @test s.a' == 0 && s.b' == 10 && s.c' == 0
        update!(s)
        @test s.a' == 0 && s.b' == 0 && s.c' == 10
        update!(s)
        @test s.a' == 0 && s.b' == 0 && s.c' == 10
    end

    @testset "distribute" begin
        @system SAccumulateDistribute(Controller) begin
            s(x=context.clock.time) => (100u"hr^-1" * x) ~ track
            d1(s) => 0.2s ~ accumulate
            d2(s) => 0.3s ~ accumulate
            d3(s) => 0.5s ~ accumulate
        end
        s = instance(SAccumulateDistribute)
        c = s.context
        @test c.clock.time' == 0u"hr" && s.s' == 0 && s.d1' == 0 && s.d2' == 0 && s.d3' == 0
        update!(s)
        @test c.clock.time' == 1u"hr" && s.s' == 100 && s.d1' == 0 && s.d2' == 0 && s.d3' == 0
        update!(s)
        @test c.clock.time' == 2u"hr" && s.s' == 200 && s.d1' == 20 && s.d2' == 30 && s.d3' == 50
        update!(s)
        @test c.clock.time' == 3u"hr" && s.s' == 300 && s.d1' == 60 && s.d2' == 90 && s.d3' == 150
        update!(s)
        @test c.clock.time' == 4u"hr" && s.s' == 400 && s.d1' == 120 && s.d2' == 180 && s.d3' == 300
    end
end