text:
module Javis

using Animations
import Cairo: CairoImageSurface, image
using FFMPEG
using Gtk
using GtkReactive
using Hungarian
using Images
using LaTeXStrings
using LightXML
import Luxor
import Luxor: Point, @layer, translate, rotate
using ProgressMeter
using Random
using Statistics
using VideoIO

const FRAMES_SYMBOL = [:same]

abstract type AbstractAction end
abstract type AbstractObject end
abstract type AbstractTransition end

include("structs/Video.jl")
include("structs/Easing.jl")
include("structs/RFrames.jl")
include("structs/GFrames.jl")
include("structs/Frames.jl")
include("structs/Scale.jl")

"""
    Transformation

Defines a transformation which can be returned by an object to be accessible later.
This is further explained in the Javis tutorials.

# Fields
- `point::Point`: the translation part of the transformation
- `angle::Float64`: the angle component of the transformation (in radians)
- `scale::Tuple{Float64, Float64}`: the scaling component of the transformation
"""
struct Transformation
    point::Point
    angle::Float64
    scale::Scale
end

Transformation(p, a) = Transformation(p, a, 1.0)
Transformation(p, a, s::Float64) = Transformation(p, a, (s, s))
Transformation(p, a, s::Tuple{Float64,Float64}) = Transformation(p, a, Scale(s...))

include("structs/ObjectSetting.jl")
include("structs/Object.jl")
include("structs/Transitions.jl")
include("structs/Action.jl")

"""
    Line

A type to define a line by two points. Can be used i.e. in [`projection`](@ref)
We mean the mathematic definition of a continuous line and not a segment of a line.

# Fields
- `p1::Point`: start point
- `p2::Point`: second point to define the line
"""
struct Line
    p1::Point
    p2::Point
end

"""
    Base.:*(m::Array{Float64,2}, transformation::Transformation)

Convert the transformation to a matrix and multiplies m*trans_matrix.
Return a new Transformation
"""
function Base.:*(m::Array{Float64,2}, transformation::Transformation)
    p = transformation.point
    θ = transformation.angle
    s = transformation.scale
    trans_matrix = [
        s.x*cos(θ) -sin(θ) p.x
        sin(θ) s.y*cos(θ) p.y
        0 0 1
    ]
    res = m * trans_matrix
    return Transformation(Point(gettranslation(res)...), getrotation(res), getscale(res))
end

include("util.jl")
include("luxor_overrides.jl")
include("backgrounds.jl")
include("svg2luxor.jl")
include("morphs.jl")
include("action_animations.jl")
include("javis_viewer.jl")
include("latex.jl")
include("object_values.jl")

"""
    projection(p::Point, l::Line)

Return the projection of a point to a line.
"""
function projection(p::Point, l::Line)
    # move line to origin and describe it as a vector
    o = l.p1
    v = l.p2 - o
    # point also moved to origin
    x = p - o

    # scalar product <x,v>/<v,v>
    c = (x.x * v.x + x.y * v.y) / (v.x^2 + v.y^2)
    return c * v + o
end

"""
    preprocess_frames!(objects::Vector{<:AbstractObject})

Computes the frames for each object and action based on the user defined frames that the
user can provide like [`RFrames`](@ref), [`GFrames`](@ref) and `:same`.

This function needs to be called before calling [`get_javis_frame`](@ref) as it computes
the actual frames for objects and actions.

# Returns
- `frames::Array{Int}` - list of all frames normally 1:...

# Warning
Shows a warning if some frames don't have a background.
"""
function preprocess_frames!(objects::Vector{<:AbstractObject})
    compute_frames!(objects)

    for (i, object) in enumerate(objects)
        compute_frames!(object.actions; parent = object, parent_counter = i)
    end

    # get all frames
    frames = Int[]
    for object in objects
        append!(frames, collect(get_frames(object)))
    end
    frames = unique(frames)
    if !(frames ⊆ CURRENT_VIDEO[1].background_frames)
        @warn(
            "Some of the frames don't have a background. In this case: $(setdiff(frames, CURRENT_VIDEO[1].background_frames)))"
        )
    end

    if isempty(CURRENT_OBJECT)
        push!(CURRENT_OBJECT, objects[1])
    else
        CURRENT_OBJECT[1] = objects[1]
    end
    return frames
end

"""
    render(
        video::Video;
        framerate=30,
        pathname="javis_GIBBERISH.gif",
        tempdirectory="",
        liveview=false,
        ffmpeg_loglevel="panic"
    )

Renders all previously defined [`Object`](@ref) drawings to the user-defined `Video` as a gif or mp4.

# Arguments
- `video::Video`: The video which defines the dimensions of the output

# Keywords
- `framerate::Int`: The frame rate of the video
- `pathname::String`: The path for the rendered gif or mp4 (i.e `output.gif` or `output.mp4`)
    - **Default:** The animation is rendered as a gif with the `javis_` prefix and some gibberish afterwards
- `tempdirectory::String`: The folder where each frame is stored
    Defaults to a temporary directory when not set
- `liveview::Bool`: Causes a live image viewer to appear to assist with animation development
- `ffmpeg_loglevel::String`:
    - Can be used if there are errors with ffmpeg. Defaults to panic:
    All other options are described here: https://ffmpeg.org/ffmpeg.html
"""
function render(
    video::Video;
    framerate = 30,
    pathname = "javis_$(randstring(7)).gif",
    liveview = false,
    tempdirectory = "",
    ffmpeg_loglevel = "panic",
)
    objects = video.objects
    frames = preprocess_frames!(objects)

    if liveview == true
        _javis_viewer(video, length(frames), objects)
        return "Live preview started."
    end

    path, ext = "", ""
    if !isempty(pathname)
        path, ext = splitext(pathname)
    end
    render_mp4 = ext == ".mp4"
    codec_props = [:priv_data => ("crf" => "22", "preset" => "medium")]
    if render_mp4
        video_io = Base.open("temp.stream", "w")
    end
    video_encoder = nothing
    # if we render a gif and the user hasn't set a tempdirectory
    if !render_mp4 && isempty(tempdirectory)
        tempdirectory = mktempdir()
    end

    filecounter = 1
    @showprogress 1 "Rendering frames..." for frame in frames
        frame_image = convert.(RGB, get_javis_frame(video, objects, frame))
        if !isempty(tempdirectory)
            Images.save("$(tempdirectory)/$(lpad(filecounter, 10, "0")).png", frame_image)
        end
        if render_mp4
            if frame == first(frames)
                video_encoder = prepareencoder(
                    frame_image,
                    framerate = framerate,
                    AVCodecContextProperties = codec_props,
                )
            end
            appendencode!(video_encoder, video_io, frame_image, filecounter)
        end
        filecounter += 1
    end

    isempty(pathname) && return
    if ext == ".gif"
        # generate a colorpalette first so ffmpeg does not have to guess it
        ffmpeg_exe(`-loglevel $(ffmpeg_loglevel) -i $(tempdirectory)/%10d.png -vf
                    palettegen $(tempdirectory)/palette.png`)
        # then apply the palette to get better results
        ffmpeg_exe(
            `-loglevel $(ffmpeg_loglevel) -framerate $framerate -i $(tempdirectory)/%10d.png -i
               $(tempdirectory)/palette.png -lavfi paletteuse -y $pathname`,
        )
    elseif ext == ".mp4"
        finishencode!(video_encoder, video_io)
        close(video_io)
        mux("temp.stream", pathname, framerate; silent = true)
    else
        @error "Currently, only gif and mp4 creation is supported. Not a $ext."
    end
    return pathname
end

"""
    get_javis_frame(video, objects, frame)

Get a frame from an animation given a video object, its objects, and frame.

If one wants to use this without calling [`render`](@ref), [`preprocess_frames!`](@ref)
needs to be called before. That way each object and action has the correct frames it should
be applied to.

# Arguments
- `video::Video`: The video which defines the dimensions of the output
- `objects::Vector{Object}`: All objects that are performed
- `frame::Int`: Specific frame to be returned

# Returns
- `Array{ARGB32, 2}` - request frame as a matrix
"""
function get_javis_frame(video, objects, frame)
    background_settings = ObjectSetting()
    Drawing(video.width, video.height, :image)
    origin()
    origin_matrix = cairotojuliamatrix(getmatrix())
    # this frame needs doing, see if each of the scenes defines it
    for object in objects
        # if object is not in global layer this sets the background_settings
        # from the parent background object
        update_object_settings!(object, background_settings)
        CURRENT_OBJECT[1] = object
        if frame in get_frames(object)
            # check if the object should be part of the global layer (i.e Background)
            # or in its own layer (default)
            in_global_layer = get(object.opts, :in_global_layer, false)::Bool
            if !in_global_layer
                @layer begin
                    draw_object(object, video, frame, origin_matrix)
                end
            else
                draw_object(object, video, frame, origin_matrix)
                # update origin_matrix as it's inside the global layer
                origin_matrix = cairotojuliamatrix(getmatrix())
            end
        end
        # if object is in global layer this changes the background settings
        update_background_settings!(background_settings, object)
    end
    img = image_as_matrix()
    finish()
    return img
end

"""
    draw_object(object, video, frame, origin_matrix)

Is called inside the [`render`](@ref) and does everything handled for an `AbstractObject`.
It is a 4-step process:
- translate to the start position
- call the relevant actions
- call the object function
- save the result of the object if wanted inside `video.defs`
"""
function draw_object(object, video, frame, origin_matrix)
    # translate the object to it's starting position.
    # It's better to draw the object always at the origin and use `star_pos` to shift it
    translate(get_position(object.start_pos))

    # reset change keywords
    empty!(object.change_keywords)

    # first compute and perform the global transformations of this object
    # relative frame number for actions
    rel_frame = frame - first(get_frames(object)) + 1
    # call currently active actions and their transformations
    for action in object.actions
        if rel_frame in get_frames(action)
            action.func(video, object, action, rel_frame)
        elseif rel_frame > last(get_frames(action)) && action.keep
            # call the action on the last frame i.e. disappeared things stay disappeared
            action.func(video, object, action, last(get_frames(action)))
        end
    end

    # set the defaults for the frame like setline() and setopacity()
    # which can depend on the actions
    set_object_defaults!(object)

    # if the scale would be 0.0 `show_object` is set to false => don't show the object
    # (it wasn't actually scaled to 0 because it would break Cairo :D)
    cs = get_current_setting()
    !cs.show_object && return

    res = object.func(video, object, frame; collect(object.change_keywords)...)
    current_global_matrix = cairotojuliamatrix(getmatrix())
    # obtain current matrix without the initial matrix part
    current_matrix = inv(origin_matrix) * current_global_matrix

    # if a transformation let's save the global coordinates
    if res isa Point
        trans = current_matrix * Transformation(res, 0.0, 1.0)
        object.result[1] = trans
    elseif res isa Transformation
        trans = current_matrix * res
        object.result[1] = trans
    else # just save the result such that it can be used as one wishes
        object.result[1] = res
    end
end

"""
    set_object_defaults!(object)

Set the default object values
- line_width and calls `Luxor.setline`.
- opacity and calls `Luxor.opacity`.
- scale and calls `Luxor.scale`.
"""
function set_object_defaults!(object)
    cs = object.current_setting
    current_line_width = cs.line_width * cs.mul_line_width
    Luxor.setline(current_line_width)
    current_opacity = cs.opacity * cs.mul_opacity
    Luxor.setopacity(current_opacity)

    desired_scale = cs.desired_scale * cs.mul_scale
    scaleto(desired_scale)
end

const LUXOR_DONT_EXPORT = [
    :boundingbox,
    :Boxmaptile,
    :Sequence,
    :setline,
    :setopacity,
    :fontsize,
    :get_fontsize,
    :scale,
    :text,
]

# Export each function from Luxor
for func in names(Luxor; imported = true)
    if !(func in LUXOR_DONT_EXPORT)
        eval(Meta.parse("import Luxor." * string(func)))
        eval(Expr(:export, func))
    end
end

export render, latex
export Video, Object, Background, Action, RFrames, GFrames
export Line, Transformation
export val, pos, ang, scl, get_value, get_position, get_angle, get_scale
export projection, morph_to
export appear, disappear, rotate_around, follow_path, change
export rev
export scaleto
export act!
export anim_translate, anim_rotate, anim_rotate_around, anim_scale

# custom override of luxor extensions
export setline, setopacity, fontsize, get_fontsize, scale, text

end


text:
using StatsBase
using PyPlot

function meshgrid(x, y)
    X = [i for i in x, j in 1:length(y)]
    Y = [j for i in 1:length(x), j in y]
    return X, Y
end

x = -5 : 0.13 : 5
y = -5 : 0.13 : 5
(X,Y) = meshgrid(x, y)
Z = X + Y*im

Mu1 = log.(Z)./(2*pi)
Mu2 = log.(-sqrt.(Z))./(2*pi)
Mu3 = log.(sqrt.(sqrt.(Z)))./(2*pi)

mean( abs.(Mu2) ./ abs.(Mu1))

theta = -pi : 0.01 : pi
r = 0.5
z = r .* exp.(im*theta)
z2 = sqrt.(z)
z3 = -z2

figure()
plot(theta,  real.(z), imag.(z) )
plot(theta,  real.(z2), imag.(z2), "--" )
axis("equal")

figure()
plot(theta, abs.(log.(z)) )
plot(theta,  abs.(log.(z2)), "--"  )
plot(theta,  abs.(log.(z3)), ":"  )

ext: jl


text:
export face_finder

"""
`face_finder(G::SimpleGraph)` will attempt to find a face of the graph.
Presently, it is *extremely* limited in what it can do.

*If* the graph `G` is planar and 3-connected, and *if* a face is
returned, it probably is an actual face in some planar embedding of `G`.

See the `README` for this module for more detail.
"""
function face_finder(G::SimpleGraph{T}) where T
    if is_acyclic(G)
        return T[]
    end
    G = trim(G,1)  # remove twigs

    comps = collect(components(G))

    for A in comps
        GG = induce(G,A)
        C = _finder_connected(GG)
        if length(C) > 0
            return C
        end
    end
    return T[]

end

function _finder_connected(G::SimpleGraph{T}) where T
    C = girth_cycle(G)
    if all(deg(G).==2)
        return C
    end

    H = deepcopy(G)
    for v in C
        delete!(H,v)
    end
    if is_connected(H)
        return C
    end
    return face_finder(H)
end


text:
#=
push!(LOAD_PATH, "<path to directory WarpFlow>") # <--- use this if Julia can't find WarpFlow
include("scripts/example1.jl")
=#
using WarpFlow
using PyPlot

# set the data and the target
dim = 1
tmpx = shuffle([rand(45); randn(175)/10 .+ .8])
nPhi = length(tmpx)
X = Array{Float64,1}[[(tmpx[i]- minimum(tmpx))/maximum(tmpx)] for i=1:nPhi]

# set the target function by overwriting the default...for v0.5 just pass a function.
target(x) = WarpFlow.targetUnif1d(x; targSig = 0.1, limit = 0.35, center = 0.5) # this sets an alias to a function giving in the file targets.jl

# generate Flow object
kappa  = Array{Float64,1}[X[Int(i)]  for i in 1:round(nPhi/4)]
nKap   = length(kappa)
y0     = Flow(kappa, WarpFlow.array1(dim, nKap), X, WarpFlow.array2eye(dim, nPhi), dim)

# function to save images
function plotrun(;fignum = 1, save = false)
	x_grd = linspace(-0.2, 1.2, 300)
	phix_grd_0  = Array{Float64,1}[[x_grd[i]] for i=1:length(x_grd)]
	Dphix_grd_0 = Array{Float64,2}[ones(1,1)  for i=1:length(x_grd)]
	yplt0 = Flow(y0.kappa, y0.eta_coeff, phix_grd_0, Dphix_grd_0, dim)
	dydt(t,y)= WarpFlow.d_forward_dt(y, sigma)
	(t1,yplt1)=WarpFlow.ode23_abstract_end(dydt,[0,1], yplt0) # Flow y0 forward to time 1
	det_grd = Float64[abs(yplt1.Dphix[i][1]) for i=1:length(x_grd)]
	den,    = target(yplt1.phix)
	est_den = det_grd.*den

	fig = figure()
	plt[:hist]([pnt[1] for pnt in X], 50, normed=1, histtype="stepfilled")
	plot(x_grd, est_den)
	if save
		plt[:savefig]("simulations/example1_v$fignum.pdf",dpi=180)
		plt[:close](fig)
	end
	return nothing
end

#  gradient ascent on eta_coeff
lambda, sigma = 1.0, 0.1
for counter = 1:35
	@time z0 = WarpFlow.get_grad(y0, target, lambda, sigma)
	y0 = y0 + 0.002 * z0
end

plotrun()

#
# #  gradient ascent on eta_coeff
# lambda, sigma = 0.1, 0.05
# for counter = 1:35
# 	@time z0 = get_grad(y0, lambda, sigma)
# 	y0 = y0 + 0.002 * z0
# end
#
# plotrun()

max_stars_repo_name: EthanAnderes/MleDensityWarp

max_issues_repo_name: EthanAnderes/MleDensityWarp

max_forks_repo_name: EthanAnderes/MleDensityWarp


text:
rand_seis_unit() = rand() > 0.5 ? "m" : "m/s"
rand_seis_cc() = rand(rand() > 0.2 ? zne : nvc)

const cDict = Dict{Char, Function}(
  'H' => i->rand_seis_cc(),
  'L' => i->rand_seis_cc(),
  'N' => i->rand_seis_cc(),
  'A' => i->rand() > 0.5 ? 'N' : 'E',
  'B' => i->'_',
  'D' => i->rand(oidfhu),
  'F' => i->rand(zne),
  'G' => i->rand(nvc),
  'I' => i->rand(oid),
  'J' => i->rand(nvc),
  'K' => i->rand(oid),
  'M' => i->rand(nvc),
  'O' => i->'_',
  'P' => i->rand(zne),
  'Q' => i->'_',
  'R' => i->'_',
  'S' => i->rand(zne),
  'T' => i->'Z',
  'U' => i->'_',
  'V' => i->'_',
  'W' => i->rand() > 0.5 ? 'S' : 'D',
  'Z' => i->rand(icfo),
  )

const uDict = Dict{Char, Function}(
    'H' => i->rand_seis_unit(),
    'N' => i->"m/s2",
    'L' => i->rand_seis_unit(),
    'A' => i->"rad",
    'B' => i->"m",
    'D' => i->"Pa",
    'F' => i->"T",
    'G' => i->"m/s2",
    'I' => i->"%",
    'J' => i->rand(junits),
    'K' => i->rand() > 0.5 ? "Cel" : "K",
    'M' => i->"m",
    'O' => i->"m/s",
    'P' => i->rand_seis_unit(),
    'Q' => i->"V",
    'R' => i->rand_seis_unit(),
    'S' => i->"m/m",
    'T' => i->"m",
    'U' => i->"%{cloud_cover}",
    'V' => i->"m3/m3",
    'W' => i->i=='S' ? "m/s" : "{direction_vector}",
    'Z' => i->rand_seis_unit(),
  )

  """
     (cha, u) = iccodes_and_units(b::Char, s::Bool)

Using band code `b`, generate quasi-sane random instrument code `i` and channel code `c`, returning channel string `cha` = `b`*`i`*`c` and unit string `u`. If `s=true`, use only seismic intrument codes.
"""
function iccodes_and_units(b::Char, s::Bool)
  i = rand(s ? hln : iclist)
  c = cDict[i](b)
  u = uDict[i](c)
  return string(b, i, c), u
end


text:
abstract type AbstractModel end
abstract type AbstractPolarizedModel <: AbstractModel end

"""
    $(TYPEDEF)
This trait specifies whether the model is a *primitive*

# Notes
This will likely turn into a trait in the future so people
can inject their models into Comrade more easily.
"""
abstract type PrimitiveTrait end

"""
    $(TYPEDEF)
Trait for primitive model
"""
struct IsPrimitive end
"""
    $(TYPEDEF)
Trait for not-primitive model
"""
struct NotPrimitive end

"""
    isprimitive(::Type)
Dispatch function that specifies whether a type is a primitive Comrade model.
This function is used for dispatch purposes when composing models.

# Notes
If a user is specifying their own model primitive model outside of Comrade they need
to specify if it is primitive

```julia
struct MyPrimitiveModel end
ComradeBase.isprimitive(::Type{MyModel}) = ComradeBase.IsPrimitive()
```

"""
function isprimitive end

@inline isprimitive(::Type{<:AbstractModel}) = NotPrimitive()

"""
    DensityAnalytic
Internal type for specifying the nature of the model functions.
Whether they can be easily evaluated pointwise analytic. This
is an internal type that may change.
"""
abstract type DensityAnalytic end

"""
    $(TYPEDEF)
Defines a trait that a states that a model is analytic.
This is usually used with an abstract model where we use
it to specify whether a model has a analytic fourier transform
and/or image.
"""
struct IsAnalytic <: DensityAnalytic end

"""
    $(TYPEDEF)
Defines a trait that a states that a model is analytic.
This is usually used with an abstract model where we use
it to specify whether a model has does not have a easy analytic
fourier transform and/or intensity function.
"""
struct NotAnalytic <: DensityAnalytic end

"""
    visanalytic(::Type{<:AbstractModel})
Determines whether the model is pointwise analytic in Fourier domain, i.e. we can evaluate
its fourier transform at an arbritrary point.

If `IsAnalytic()` then it will try to call `visibility_point` to calculate the complex visibilities.
Otherwise it fallback to using the FFT that works for all models that can compute an image.

"""
@inline visanalytic(::Type{<:AbstractModel}) = NotAnalytic()

"""
    imanalytic(::Type{<:AbstractModel})
Determines whether the model is pointwise analytic in the image domain, i.e. we can evaluate
its intensity at an arbritrary point.

If `IsAnalytic()` then it will try to call `intensity_point` to calculate the intensity.
"""
@inline imanalytic(::Type{<:AbstractModel}) = IsAnalytic()

#=
    This is internal function definitions for how to
    compose whether a model is analytic. We need this
    for composite models.
=#
@inline Base.:*(::IsAnalytic, ::IsAnalytic) = IsAnalytic()
@inline Base.:*(::IsAnalytic, ::NotAnalytic) = NotAnalytic()
@inline Base.:*(::NotAnalytic, ::IsAnalytic) = NotAnalytic()
@inline Base.:*(::NotAnalytic, ::NotAnalytic) = NotAnalytic()

"""
    $(SIGNATURES)
Function that computes the pointwise visibility if the model has the trait
in the fourier domain `IsAnalytic()`. Otherwise it will use the FFTW fallback.
"""
function visibility_point end

"""
    $(SIGNATURES)
Function that computes the pointwise intensity if the model has the trait
in the image domain `IsAnalytic()`. Otherwise it will use construct the image in visibility
space and invert it.
"""
function intensity_point end

"""
    $(SIGNATURES)
Computes the intensity map of model by modifying the input IntensityMap object
"""
function intensitymap! end

"""
    $(SIGNATURES)
Computes the intensity map of model. This version requires additional information to
construct the grid.

# Example

```julia
m = Gaussian()
# field of view
fovx, fovy = 5.0
fovy = 5.0
# number of pixels
nx, ny = 128

img = intensitymap(m, fovx, fovy, nx, ny; pulse=DeltaPulse())
```
"""
function intensitymap end

# """
#     ($SIGNATURES)
# A threaded version of the intensitymap function.

# # Notes
# If using autodiff this won't play well with zygote given the mutation. This will be fixed in
# future versions.
# """
# function tintensitymap end

# """
#     ($SIGNATURES)
# A threaded version of the intensitymap! function.

# # Notes
# If using autodiff this won't play well with zygote given the mutation. This will be fixed in
# future versions.
# """
# function tintensitymap! end


text:
#TODO: higher order gaps

################################################################################
### DISPERSION MANAGERS
################################################################################
#@code_warntype checked
"""
    build_dispersions!(sim; num_wg_bands=15, num_bloch=17, num_free_bands=2)
"""
function build_dispersions!(sim::Simulation;
    interpolation=:cubic, β_start=0, β_stop=20, nβ=100,
    num_wg_modes=maximum([sim.sct.channels[i].quantum_number for i ∈ eachindex(sim.sct.channels)])+1,
    free_band_zone=:half, parallel=nprocs()>1, num_bloch=17, num_free_bands=2, num_wg_bands_multiplier=num_free_bands)

    num_wg = length(sim.sys.waveguides)
    if isempty(sim.sct.channels)
        throw(ErrorException("no channels specified"))
    end
    for i ∈ eachindex(sim.sys.waveguides)
        ch_inds = findall([sim.sct.channels[j].waveguide==sim.sys.waveguides[i] for j ∈ eachindex(sim.sct.channels)])
        if any(isempty.([sim.sct.channels[j].dispersion for j ∈ ch_inds]))
            println("Building dispersion for waveguide $i")

            wg_bands, ks, bands, gaps = waveguide_dispersion(sim, sim.sys.waveguides[i];
            interpolation=interpolation, β_start=β_start, β_stop=β_stop, nβ=nβ,
            num_wg_modes=num_wg_modes, free_band_zone=free_band_zone, parallel=parallel,
            num_free_bands=num_free_bands, num_bloch=num_bloch, num_wg_bands_multiplier=num_wg_bands_multiplier)

            if isempty(wg_bands)
                throw(ErrorException("no guided modes found. Either there are none in the first bandgap, or increase num_wg_modes from $num_wg_modes or num_wg_bands_multiplier from $num_wg_bands_multiplier"))
            end
            for j ∈ ch_inds
                if isempty(sim.sct.channels[j].dispersion)
                    if length(wg_bands) < sim.sct.channels[j].quantum_number
                        throw(ErrorException("channel $j (waveguide $i) has higher mode number than computed. ensure the higher mode exists in the regime you are working."))
                    end
                    push!(sim.sct.channels[j].dispersion, wg_bands[sim.sct.channels[j].quantum_number])
                    if !isempty(gaps[1])
                        for l ∈ eachindex(gaps[1])
                            push!(sim.sct.channels[j].gaps, [gaps[1][l], gaps[2][l]])
                        end
                    else
                        nothing
                    end
                end
            end
        end
    end
    return nothing
end

#@code_warntype checked
"""
    remove_dispersions!(sim)

removes dispersion curves from simulation
"""
function remove_dispersions!(sim::Simulation)
    for i ∈ 1:length(sim.sct.channels)
        for j ∈ 1:length(sim.sct.channels[i].dispersion)
            pop!(sim.sct.channels[i].dispersion)
        end
    end
    return nothing
end

"""
    waveguide_dispersion(sim, waveguide)
"""
function waveguide_dispersion(sim::Simulation, waveguide::Int;
    interpolation=:cubic, β_start=0, β_stop=20, nβ=31,
    num_wg_modes=maximum([sim.sct.channels[i].quantum_number for i ∈ eachindex(sim.sct.channels)])+1,
    free_band_zone=:half, parallel=nprocs()>1, num_bloch=17, num_free_bands=2,
    num_wg_bands_multiplier=2, disp_opt=true)

    if isempty(sim.sys.waveguides)
        nothing
    elseif ishalfspace(sim.sys.domains[ScalarFDFD.get_waveguide_domains(sim, waveguide)][1])
        wg_dispersion, ks, bands, gaps = ScalarFDFD.halfspace_waveguide_dispersion(sim, waveguide)
    elseif any(ispc.(sim.sys.domains[get_waveguide_domains(sim,waveguide)]))
        wg_dispersion, ks, bands, gaps = ScalarFDFD.pc_waveguide_dispersion(sim, waveguide;
            num_wg_bands_multiplier=num_wg_bands_multiplier, num_bloch=num_bloch,
            num_free_bands=num_free_bands, parallel=parallel, free_band_zone=free_band_zone,
            interpolation=interpolation, disp_opt=disp_opt)
    elseif ScalarFDFD.isplanar(sim.sys.domains[ScalarFDFD.get_waveguide_domains(sim, waveguide)[1]])
        wg_sim = extract_waveguide_simulation(sim, waveguide)
        wg_dispersion, ks, bands, gaps = ScalarFDFD.planar_waveguide_dispersion(wg_sim;
            nβ=nβ, num_wg_modes=num_wg_modes, β_start=β_start, β_stop=β_stop, interpolation=interpolation)
    else
        throw(ArgumentError("unrecognized waveguide type"))
    end

    return wg_dispersion, ks, bands, gaps
end

################################################################################
### WAVEGUIDE DISPERSIONS
################################################################################
"""
    wg_dispersion, bands, gaps, ks = trivial_waveguide_dispersion(sim, waveguide, num_bloch)
"""
function halfspace_waveguide_dispersion(sim::Simulation, waveguide)

    wg_dispersion = Array{AbstractInterpolation}(undef,1)

    ks = 0 .+ 0:1
    wg_inds = ScalarFDFD.get_waveguide_domains(sim, waveguide)

    bands = Array{Float64}(undef, 2, 1)
    bands[:,1] = sim.sys.domains[wg_inds[1]].n₁*ks

    wg_dispersion[1] = extrapolate(interpolate(bands[:,1], BSpline(Linear())), Line())
    wg_dispersion[1] = scale(wg_dispersion[1], ks)

    return wg_dispersion, ks, Array{Float64}(undef,0,0), (Float64[], Float64[])
end

"""
    wg_bands = pc_waveguide_dispersion(sim::Simulation, waveguide, gaps; num_free_bands=2, num_wg_bands=15, num_bloch=17, interpolation=:cubic)
"""
function pc_waveguide_dispersion(sim::Simulation, waveguide::Int; num_wg_bands_multiplier,
    num_bloch, num_free_bands, interpolation, parallel, free_band_zone, disp_opt=true)

    wg_sim = extract_waveguide_simulation(sim, waveguide)
    num_wg_bands=round(Int,num_wg_bands_multiplier*max(wg_sim.lat.a/wg_sim.lat.b,wg_sim.lat.b/wg_sim.lat.a))
    pc_sim = Simulation(wg_sim.sys.domains[2]; bnd=Boundary(bc=:p), dis=wg_sim.dis)

    if wg_sim.lat.a < wg_sim.lat.b
        ka_bloch = -π/wg_sim.lat.a:2π/wg_sim.lat.a/(num_bloch-1):π/wg_sim.lat.a
        kb_bloch = fill(0,length(ka_bloch))
        k_bloch = ka_bloch
    else
        kb_bloch = -π/wg_sim.lat.b:2π/wg_sim.lat.b/(num_bloch-1):π/wg_sim.lat.b
        ka_bloch = fill(0,length(kb_bloch))
        k_bloch = kb_bloch
    end

    _, gaps = band_structure(pc_sim, num_bloch; num_bands=num_free_bands, parallel=parallel, interpolation=interpolation, zone=free_band_zone, calc_type="waveguide band structure ", disp_opt=disp_opt)
    free_bands, _, _ = band_structure(pc_sim, (ka_bloch, kb_bloch); num_bands=num_free_bands, parallel=parallel, interpolation=interpolation, pg=Progress(Int(1e5), 1e5), disp_opt=disp_opt)
    _, _, bands = band_structure(wg_sim, (ka_bloch, kb_bloch); num_bands=num_wg_bands, parallel=parallel, calc_type="waveguide dispersion ", disp_opt=disp_opt)

    mode_bool = fill(false, num_wg_bands)
    for i ∈ eachindex(bands)
        # this is a crude way to determine guided modes TODO
        gap_size = gaps[2][1]-gaps[1][1]
        mode_bool[i] = any(gaps[1][1]+.1*gap_size .< bands[i].(1:.01:num_bloch) .< gaps[2][1]-.1*gap_size)
    end
    mode_inds = findall(mode_bool)

    wg_dispersion = Array{AbstractInterpolation}(undef,length(mode_inds))
    for j in 1:length(mode_inds)
        wg_dispersion[j] = scale(bands[mode_inds[j]], k_bloch)
    end
    return wg_dispersion, k_bloch, free_bands, gaps
end

#@code_warntype check
"""
    wg_dispersion, ks, bands, gaps = planar_waveguide_dispersion(sim, nβ, num_wg_modes, k_start, k_stop, interpolation)
"""
function planar_waveguide_dispersion(sim::Simulation; nβ, num_wg_modes, β_start, β_stop, interpolation)

    K = Array{ComplexF64}(undef, nβ, num_wg_modes)
    ks = β_start .+ (0:nβ-1)*(β_stop-β_start)/(nβ-1)
    for i ∈ eachindex(ks)
        K[i,:] = ScalarFDFD.planar_kl(sim, ks[i], .001, num_wg_modes)[1]
    end

    wg_dispersion = Array{AbstractInterpolation}(undef, num_wg_modes)
    if interpolation == :linear
        interpolator = BSpline(Linear())
    elseif interpolation == :quadratic
        interpolator = BSpline(Quadratic(Line(OnGrid())))
    elseif interpolation == :cubic
        interpolator = BSpline(Cubic(Line(OnGrid())))
    else
        throw(ArgumentError("invalid order $(interpolation)"))
    end
    for i ∈ eachindex(wg_dispersion)
        wg_dispersion[i] = extrapolate(interpolate(real(K[:,i]), interpolator), Line())
        wg_dispersion[i] = scale(wg_dispersion[i], ks)
    end

    return wg_dispersion, ks, Array{Float64}(undef,0,0), (Float64[], Float64[])
end


text:
using LaserTag
using Random
using Test
using POMDPModels
using POMDPModelTools
using POMDPPolicies
using POMDPSimulators
using ParticleFilters
using POMDPs

@time p = gen_lasertag()

pol = RandomPolicy(p, rng=MersenneTwister(1))

sim = HistoryRecorder(max_steps=10, rng=MersenneTwister(2))

filter = SIRParticleFilter(p, 10000)

hist = simulate(sim, p, pol, filter)

tikz_pic(LaserTagVis(p))

# discrete
p = gen_lasertag()

# check observation model consistency
rng = MersenneTwister(12)
N = 1_000_000
s = initialstate(p, rng)
od = observation(p, s)
obs = [rand(rng, od) for i in 1:N]
for dir in 1:8
    counts = Dict{Int,Int}()
    for i in 1:N
        o = obs[i]
        if haskey(counts, o[dir])
            counts[o[dir]] += 1
        else
            counts[o[dir]] = 1
        end
    end
    total = 0
    for (r, count) in counts
        total += count
        if r == 0
            prob = LaserTag.cdf(od.model.cdf, dir, n_clear_cells(od.distances, dir), 0)
        else
            prob = (LaserTag.cdf(od.model.cdf, dir, n_clear_cells(od.distances, dir), r) -
                 LaserTag.cdf(od.model.cdf, dir, n_clear_cells(od.distances, dir), r-1))
        end
        try
            @test isapprox(prob*N, count, atol=10, rtol=0.1)
        catch ex
            @show prob*N
            @show count
            rethrow(ex)
        end
    end
    @test total == N
end

pol = RandomPolicy(p, rng=MersenneTwister(1))

sim = HistoryRecorder(max_steps=10, rng=MersenneTwister(2))

filter = SIRParticleFilter(p, 10000)

hist = simulate(sim, p, pol, filter)

tikz_pic(LaserTagVis(p))
render(p, first(eachstep(hist)))
io = IOBuffer()
show(io, MIME("image/png"), render(p, first(eachstep(hist))))

s = initialstate(p, MersenneTwister(4))
@inferred generate_sor(p, s, 1, MersenneTwister(4))

sp, o, r = generate_sor(p, s, 1, MersenneTwister(4))
@inferred observation(p, s, 1, sp)

show(stdout, MIME("text/plain"), LaserTagVis(cpp_emu_lasertag(4)))

include("emulate_cpp_lasertag.jl")


text:
using FermiCG
using Printf
using Test
using LinearAlgebra

#@testset "ClusteredTerms" begin
    atoms = []
    push!(atoms,Atom(1,"H",[0,0,0]))
    push!(atoms,Atom(2,"H",[0,1,0]))
    push!(atoms,Atom(3,"H",[0,0,2]))
    push!(atoms,Atom(4,"H",[0,0,3]))
    push!(atoms,Atom(5,"H",[0,0.4,4]))
    push!(atoms,Atom(6,"H",[0,0,5]))
    push!(atoms,Atom(7,"H",[0,0,6]))
    push!(atoms,Atom(8,"H",[0,0,7]))
    push!(atoms,Atom(9,"H",[0,0,8]))
    push!(atoms,Atom(10,"H",[0,0,9]))
    push!(atoms,Atom(11,"H",[0,0,10]))
    push!(atoms,Atom(12,"H",[0,0,11]))
    basis = "sto-3g"
    mol     = Molecule(0,1,atoms,basis)

    mf = FermiCG.pyscf_do_scf(mol)
    nbas = size(mf.mo_coeff)[1]
    ints = FermiCG.pyscf_build_ints(mol,mf.mo_coeff, zeros(nbas,nbas));
    #e_fci, d1_fci, d2_fci = FermiCG.pyscf_fci(ints,6,6)
    #@printf(" FCI Energy: %12.8f\n", e_fci)

    C = mf.mo_coeff
    Cl = FermiCG.localize(mf.mo_coeff,"lowdin",mf)
    FermiCG.pyscf_write_molden(mol,Cl,filename="lowdin.molden")
    S = FermiCG.get_ovlp(mf)
    U =  C' * S * Cl
    println(" Build Integrals")
    flush(stdout)
    ints = FermiCG.orbital_rotation(ints,U)
    println(" done.")
    flush(stdout)

    #clusters    = reverse([(1:4),(5:8),(9:10),(11:12)])
    #init_fspace = reverse([(2,2),(2,2),(1,1),(1,1)])
    clusters    = [(1:4),(5:8),(9:10),(11:12)]

    max_roots = 3

    clusters = [Cluster(i,collect(clusters[i])) for i = 1:length(clusters)]

    cluster_bases = FermiCG.compute_cluster_eigenbasis(ints, clusters, verbose=0, max_roots=max_roots)

    display.(clusters)
    terms = FermiCG.extract_ClusteredTerms(ints, clusters)

#    for t in keys(terms)
#        #FermiCG.print_fock_sectors(collect(t))
#        print(t)
#        for tt in terms[t]
#            display(tt)
#            display(typeof(tt))
#        end
#    end

    cluster_ops = FermiCG.compute_cluster_ops(cluster_bases, ints);

    ci_vector = FermiCG.ClusteredState(clusters)
    FermiCG.add_fockconfig!(ci_vector,[(2,2),(2,2),(1,1),(1,1)])
    FermiCG.add_fockconfig!(ci_vector,[(1,2),(3,2),(2,1),(1,1)])
    FermiCG.add_fockconfig!(ci_vector,[(2,2),(2,2),(2,1),(1,1)])
    FermiCG.add_fockconfig!(ci_vector,[(2,1),(2,2),(1,2),(1,1)])
    FermiCG.add_fockconfig!(ci_vector,[(2,1),(2,2),(1,2),(1,1)])

    #FermiCG.add_fockconfig!(ci_vector,[(2,2),(0,2),(1,1),(1,1)])
    #FermiCG.add_fockconfig!(ci_vector,[(0,2),(2,2),(1,1),(1,1)])
    #FermiCG.add_fockconfig!(ci_vector,[(2,2),(0,1),(1,2),(1,1)])
    #FermiCG.add_fockconfig!(ci_vector,[(0,2),(2,1),(1,2),(1,1)])
    #FermiCG.add_fockconfig!(ci_vector,[(3,1),(1,3),(1,1),(1,1)])
    #FermiCG.add_fockconfig!(ci_vector,[(1,3),(3,1),(1,1),(1,1)])
    #FermiCG.add_fockconfig!(ci_vector,[(3,1),(2,2),(0,2),(1,1)])
    #FermiCG.add_fockconfig!(ci_vector,[(1,3),(2,2),(2,0),(1,1)])

    #FermiCG.add_fockconfig!(ci_vector,[(1,3),(2,1),(3,1),(1,1)])
    #FermiCG.add_fockconfig!(ci_vector,[(0,2),(4,2),(1,1),(1,1)])
    #FermiCG.add_fockconfig!(ci_vector,[(3,2),(1,2),(0,1),(1,0)])
    #FermiCG.add_fockconfig!(ci_vector,[(3,2),(2,2),(0,1),(0,0)])
    #FermiCG.add_fockconfig!(ci_vector,reverse([(2,2),(2,2),(1,1),(0,0)]))
    #FermiCG.add_fockconfig!(ci_vector,reverse([(3,2),(1,2),(1,1),(0,0)]))

    #FermiCG.expand_each_fock_space!(ci_vector, cluster_bases)
    FermiCG.expand_to_full_space!(ci_vector, cluster_bases)

    display(ci_vector)
    #display(cluster_bases[2][(2,2)])

    function build(ci_vector, cluster_ops, terms)
        dim = length(ci_vector)
        H = zeros(dim, dim)

        bra_idx = 0
        for (fock_bra, configs_bra) in ci_vector.data
            for (config_bra, coeff_bra) in configs_bra
                bra_idx += 1
                ket_idx = 0
                for (fock_ket, configs_ket) in ci_vector.data
                    fock_trans = fock_bra - fock_ket

                    # check if transition is connected by H
                    if haskey(terms, fock_trans) == false
                        ket_idx += length(configs_ket)
                        continue
                    end

                    for (config_ket, coeff_ket) in configs_ket
                        ket_idx += 1
                        #ket_idx <= bra_idx || continue

                        for term in terms[fock_trans]
                            H[bra_idx, ket_idx] += FermiCG.contract_matrix_element(term, cluster_ops, fock_bra, config_bra, fock_ket, config_ket)
                        end

                        #H[ket_idx, bra_idx] = H[bra_idx, ket_idx]

                    end
                end
            end
        end
        return H
    end

    if false
        display(keys(terms))
        tmp = FermiCG.TransferConfig([(2,0),(-2,0),(0,0),(0,0)])
        display(haskey(terms,tmp))
        for i in terms[tmp]
            display(i)
        end

        throw(Exception)
    end
    @time H = build(ci_vector, cluster_ops, terms)

    maximum(abs.(H-H')) < 1e-14 || error("Hamiltonian not symmetric: ",maximum(abs.(H-H')))
    F = eigen(H)
    for (idx,Fi) in enumerate(F.values[1:min(10,length(F.values))])
        @printf(" %4i %18.13f\n", idx, Fi)
    end
    #display(round.(H[1:16,17:32]; digits=8))
    #display(round.(H[1:16,1:16]; digits=8))
    #    for i in 1:size(H,1)
#        #@printf("%18.12f\n", H[i,1])
#    end
#    fock = FermiCG.FockConfig([(2,2),(2,2),(1,1),(1,1)])
##    fock_bra = FermiCG.FockConfig([(2,2),(2,2),(1,1),(0,0)])
##    fock_ket = FermiCG.FockConfig([(3,2),(1,2),(1,1),(0,0)])
##    fock_trans = fock_bra - fock_ket
#    fock_trans = fock - fock
#    config_bra = [2,2,1,1]
#    config_ket = [1,1,1,1]
#    test = 0
#    for term in  terms[fock_trans]
#        test1 = FermiCG.contract_matrix_element(term, cluster_ops, fock, config_bra, fock, config_ket)
#        println(test1/2, term.ops)
#        global test += test1
#    end
#    println(test)
#    FermiCG.print_configs(ci_vector)
#end

max_stars_repo_name: nmayhall-vt/FermiCG

max_issues_repo_name: nmayhall-vt/FermiCG

max_forks_repo_name: nmayhall-vt/FermiCG


text:
# Convert a dims argument to a Pattern
_pattern(dims::Pattern) = dims
_pattern(dims::Tuple) = Pattern(dims)
_pattern(::Colon) = Pattern(:__)
_pattern(dims) = Pattern(:__, dims)


text:
#Generate an ouput ImagineSignal for commanding the piezo to repeat the sample sequence `one_cyc` `ncycles` times.
#Appends those samples to `pos` and also returns dummy camera and laser signals so that Imagine doesn't complain
#`ncycles` default to 20 about seconds-worth of cycles based on the sample rate of `pos`
function pos_commands(pos::ImagineSignal, one_cyc, ncycles=ceil(Int, 20.0 / ustrip(inv(samprate(pos))*length(one_cyc))))
    @assert isempty(pos)
    rig_sigs = rigtemplate(rig_name(pos); sample_rate = samprate(pos))
    add_sequence!(pos, "bidi_cycle", one_cyc)
    cycle_nsamps = length(one_cyc)
    cam = first(getcameras(rig_sigs)) #not used
    las = first(getlasers(rig_sigs)) #not used
    add_sequence!(cam, "quiet", falses(cycle_nsamps)) #just to pass validation
    for i = 1:ncycles
        append!(pos, "bidi_cycle")
        append!(cam, "quiet")
        append!(las, "quiet")
    end
    pos_mon = getname(rig_sigs, monitor_name(pos))
    return [pos;cam;las;pos_mon]
end

function pos_commands(rig::AbstractString, pos_name::AbstractString, pstart, pstop, stack_rate, ncycles=ceil(Int, 20.0/ustrip(inv(stack_rate)));
                      sample_rate = 100000Hz, lp_cutoff=3.5*stack_rate)
    mod_cyc = vcat(gen_bidi_pos(pstart, pstop, 1/stack_rate, sample_rate; lp_cutoff=lp_cutoff)...)
    return pos_commands(rig, pos_name, mod_cyc, ncycles; sample_rate = sample_rate)
end

function pos_commands(rig::AbstractString, pos_name::AbstractString, mod_cyc, ncycles=ceil(Int, 20.0/ustrip(inv(stack_rate))); sample_rate = 100000Hz)
    rig_sigs = rigtemplate(rig; sample_rate = sample_rate)
    pos = getname(rig_sigs, pos_name)
    return pos_commands(pos, mod_cyc, ncycles)
end

function append_warmup!(pos, cam, las, high_las, pos_cycle, ncycs_warmup)
    nsamps_cyc = Base.length(pos_cycle) #should always be even
    add_sequence!(cam, "cam_warmup", falses(nsamps_cyc))
    add_sequence!(las, "las_warmup", falses(nsamps_cyc))
    add_sequence!(pos, "pos_warmup", pos_cycle)
    if rig_name(pos) == "ocpi-2"
        add_sequence!(high_las, "high_las_warmup", trues(nsamps_cyc))
    end
    for i = 1:ncycs_warmup
        append!(pos, "pos_warmup")
        append!(cam, "cam_warmup")
        append!(las, "las_warmup")
        append!(high_las, "high_las_warmup")
    end
	return pos, cam, las, high_las
end

function get_trial_intervals(ctr_i, lag, sample_rate, nsamps_exp, nsamps_flash)
	nsamps_from_ctr = div(nsamps_flash,2)
	nsamps_lag = ImagineInterface.calc_num_samps(lag, sample_rate)
    flash_start = ctr_i-nsamps_lag-nsamps_from_ctr
    flash_stop = ctr_i-nsamps_lag+nsamps_from_ctr
    exp_stop = flash_stop
    exp_start = flash_stop-nsamps_exp+1
    flash_itv = ClosedInterval(flash_start, flash_stop)
    exp_itv = ClosedInterval(exp_start, exp_stop)
    return exp_itv, flash_itv
end

function check_overlap(itvs::Vector{T}) where {T<:ClosedInterval}
    for i = 2:length(itvs)
        if !isempty(intersect(itvs[i-1], itvs[i]))
            error("Pulse intervals overlap")
        end
    end
    return itvs
end

#given a measured piezo cycle, a set of desired slice locations `slice_zs`, and a set of forward and reverse lags,
#return a tuple (cam_samps, las_samps, shift) with sample vectors for controlling camera exposure and laser pulsing
function get_cyc_pulses(mon_cyc, slice_zs, fwd_lags, back_lags, exp_time, flash_time, sample_rate::HasInverseTimeUnits)
	nsamps_exp = ImagineInterface.calc_num_samps(exp_time, sample_rate)
	nsamps_flash = ImagineInterface.calc_num_samps(flash_time, sample_rate)
    @assert iseven(length(mon_cyc))
    @assert issorted(slice_zs) || issorted(slice_zs; rev=true)
    nsamps_sweep = div(length(mon_cyc),2)
    cams_fwd = ClosedInterval{Int}[]
    cams_back = ClosedInterval{Int}[]
    lasers_fwd = ClosedInterval{Int}[]
    lasers_back = ClosedInterval{Int}[]
    nslices = length(slice_zs)
    #get_trial_intervals can return intervals lying outside the cycle bounds (negative first indices and last indices greater than 2*nsamps_sweep)
    #check extrema of intervals (smin and smax below)
    #if both of these extrema lie outside of 2*nsamps_sweep, throw an error
    #if only one of these does, it's okay.  return cam_samps, las_samps, and the shift value
    for s = 1:nslices
        ctr_fwd, ctr_back = find_bidi_locations(mon_cyc, sample_rate, slice_zs[s])
        cam_fwd, las_fwd = get_trial_intervals(ctr_fwd, fwd_lags[s], sample_rate, nsamps_exp, nsamps_flash)
        cam_back, las_back = get_trial_intervals(ctr_back, back_lags[s], sample_rate, nsamps_exp, nsamps_flash)
        #cam_back+=nsamps_sweep
        #las_back+=nsamps_sweep
        push!(cams_fwd, cam_fwd)
        push!(cams_back, cam_back)
        push!(lasers_fwd, las_fwd)
        push!(lasers_back, las_back)
    end
    las_itvs = vcat(lasers_fwd, lasers_back)
    cam_itvs = vcat(cams_fwd, cams_back)
    shift = 0
    all_itvs = vcat(las_itvs, cam_itvs)
    smin = minimum(minimum.(all_itvs))
    smax = maximum(maximum.(all_itvs))
    if smin < 1 && smax > 2*nsamps_sweep
        error("Cannot create pulse sequence without overlapping pulses across cycles")
    elseif !(smin >0 && (smax <= 2*nsamps_sweep))
        if abs(smin-1) > abs(smax - 2*nsamps_sweep)
            @show shift = smin-1
        else
            @show shift = smax - 2*nsamps_sweep
        end
    end
    las_samps = ImagineInterface.gen_pulses(2*nsamps_sweep, check_overlap(las_itvs.-shift))
    cam_samps = ImagineInterface.gen_pulses(2*nsamps_sweep, check_overlap(cam_itvs.-shift))
    return cam_samps, las_samps, shift #shift is the shift that should be applied to sample vectors to get them into proper alignment
end

function append_template!(pos, cam, las, high_las, slice_zs, exp_time, flash_time; freq=0.1Hz) #appends a stack sequence to each signal input. Images only on the forward sweep.
    #pad by 5μm if we have room
    z_min = minimum(slice_zs)
    z_max = maximum(slice_zs)
    z_pad = minimum([abs(z_min - 0.0μm), abs(800μm - z_max), 5.0μm])
    z_min-=z_pad
    z_max+=z_pad
    mod_cyc = vcat(gen_bidi_pos(z_min, z_max, 1/freq, samprate(pos))...)
    z_range = z_max - z_min #should also match the command (not necessarily the monitored signal)
    pos_speed = freq * z_range #in microns per second
    pad_nsecs = z_pad / pos_speed #nsecs to cover the z_pad distance
    pad_nsamps = ImagineInterface.calc_num_samps(pad_nsecs, samprate(pos))
    flash_cyc, cam_cyc, monshft = flash_cam_cycs(mod_cyc, mod_cyc, slice_zs, pad_nsamps, flash_time, exp_time, samprate(pos))
    #only take pictures on the forward sweep
    midpoint = div(length(flash_cyc),2) + 1
    flash_cyc[midpoint:end] .= false
    cam_cyc[midpoint:end] .= false
    append!(pos, "pos_template", mod_cyc)
    append!(cam, "cam_template", cam_cyc)
    append!(las, "las_template", flash_cyc)
    append!(high_las, "highlas_template", trues(length(mod_cyc)))
end

function multistack_bidi(cam_name, las_name, pos_name, pos_cyc, mon_cyc, nstacks_fwd, slice_zs, fwd_toffsets, back_toffsets, exp_time, flash_time, sample_rate, ncycs_warmup)
	ocpi2 = rigtemplate("ocpi-2"; sample_rate = sample_rate)
	pos_new = getname(ocpi2, pos_name)
	cam_new = getname(ocpi2, cam_name)
	las_new = getname(ocpi2, "all lasers")
	high_las_new = getname(ocpi2, las_name)

	FastScanning.append_warmup!(pos_new, cam_new, las_new, high_las_new, pos_cyc, ncycs_warmup)

    @show nsamps_offset = ImagineAnalyses.mon_delay(pos_cyc, mon_cyc)
	#temporal shift to make calculations easier (will shift things back later)
    mon_cyc = circshift(mon_cyc, -nsamps_offset)
	if nsamps_offset < 0
		error("can't handle negative offset")
	elseif nsamps_offset > 0
		append!(las_new, "delay", fill(false, nsamps_offset))
		append!(cam_new, "delay")
		append!(high_las_new, "delayhigh", fill(true, nsamps_offset))
	end
	cam_samps, las_samps, shift = get_cyc_pulses(mon_cyc, slice_zs, fwd_toffsets, back_toffsets, exp_time, flash_time, sample_rate)
	@assert shift == 0 #TODO: handle this

	nsamps_cyc = length(pos_cyc)
	pos_cyc_name = "pos_cyc"
	add_sequence!(pos_new, pos_cyc_name, pos_cyc)
	add_sequence!(cam_new, "cam_cyc", cam_samps)
	add_sequence!(las_new, "las_cyc", las_samps)
	add_sequence!(high_las_new, "high_cyc", trues(nsamps_cyc))

	for i = 1:nstacks_fwd
		append!(pos_new, "pos_cyc")
		append!(cam_new, "cam_cyc")
		append!(las_new, "las_cyc")
		append!(high_las_new, "high_cyc")
	end
    if nsamps_offset > 0 #add another positioner cycle because we have too many cam and las samples
        append!(pos_new, pos_cyc_name)
        append!(cam_new, "extra_cyc", fill(false, nsamps_cyc-nsamps_offset))
        append!(las_new, "extra_cyc")
        append!(high_las_new, "extra_cychigh", fill(true, nsamps_cyc-nsamps_offset))
    end
	return [pos_new; cam_new; las_new; high_las_new]
end

function slicetiming_commands!(pos::T1, mod_cyc, mon_cyc, las, high_las, cam, lags, slice_zs) where {T1<:ImagineSignal}
    nsamps_cyc = Base.length(mod_cyc) #should always be even
    nsamps_sweep = div(nsamps_cyc,2)
    ncycs_ignore = ceil(Int, (ustrip(samprate(pos))*5) / nsamps_cyc) #5 seconds of warmup
    #@show nsamps_offset = mon_lag_nsamps(mod_cyc, mon_cyc)
    @show nsamps_offset = ImagineAnalyses.mon_delay(mod_cyc, mon_cyc)
	#temporal shift to make calculations easier (will shift things back later)
    mon_cyc = circshift(mon_cyc, -nsamps_offset)
    ntrials=1 #number of trials per lag value per z slice
    #We want to deliver the laser pulse during the global shutter time (tglobal in PCO camera manual)
    #If the user images with full chip then the global exposure segment doesn't begin until _10ms_ after the exposure starts.
    #Therefore with a 1ms laser pulse we need an 11ms exposure time, and the pulse is placed in the last 1ms
    exp_time = 0.002s
    flash_time = 0.001s
    empty!(pos)
	if nsamps_offset < 0
		error("Cannot handle negative offset (but could with a bit more coding work)")
	elseif nsamps_offset > 0
		append!(las, "delay", fill(false, nsamps_offset))
		append!(cam, "delay")
		append!(high_las, "delayhigh", fill(true, nsamps_offset))
	end
    nsamps_flash = ImagineInterface.calc_num_samps(flash_time, samprate(pos))
    if iseven(nsamps_flash)
	    nsamps_flash -= 1 #force odd so that fwd and reverse flashes are easier to align
    end
	nsamps_exp = ImagineInterface.calc_num_samps(exp_time, samprate(pos))

	#Sequence from here:
	# one slow image stack
	# several fast "warmup" piezo cycles with no camera / laser activity (ncycs_ignore)
	# ntrials * nlags piezo cycles in which we take a pair of slices for each cycle
	# repeat above step for each slice, so a total of ntrials * nlags * nslices cycles (and twice that many images)
    ncycs_ignore = ceil(Int, (ustrip(samprate(pos))*5) / nsamps_cyc) #5 seconds of warmup
	append_template!(pos, cam, las, high_las, slice_zs, exp_time, flash_time) #appends slow stack sequence to each signal input
	append_warmup!(pos, cam, las, high_las, mod_cyc, ncycs_ignore) #note: the warmup cycles are named differently but sequence is same as other fast pos cycles

    pos_cyc_name = "fast_pos_cyc"
	add_sequence!(pos, pos_cyc_name, mod_cyc)
	add_sequence!(high_las, "laser_high", trues(nsamps_cyc)) #will append one of these for every cycle as we go
	#Now iterate through slices, lags, and trials
	for s = 1:length(slice_zs)
		fwd_ctr, back_ctr = find_bidi_locations(mon_cyc, samprate(pos), slice_zs[s])
		for l = 1:Base.length(lags)
            cyc_exp_name = "cyc_exp_$(l)_slice_$s"
            cyc_flash_name = "cyc_flash_$(l)_slice_$s"
            has_partial = false #do we need to add a partial cycle before and after each cycle due to exposure pulse timing?
            cam_samps, las_samps, timeshift = get_cyc_pulses(mon_cyc, [slice_zs[s]], [lags[l];], [lags[l];], exp_time, flash_time, samprate(pos))
			add_sequence!(cam, cyc_exp_name, cam_samps)
			add_sequence!(las, cyc_flash_name, las_samps)
            #if timeshift is nonzero then the imaging pulses lie partially outside the cycle and we need an extra cycle
            if timeshift < 0
                warn("Adding an extra piezo cycle in slice $s with lag $l for exposures or flashes that begin before the cycle")
                has_partial = true
                add_sequence!(cam, "cam_partial_before_$(s)_$(l)", falses(nsamps_cyc+timeshift))
                add_sequence!(las, "las_partial_before_$(s)_$(l)", falses(nsamps_cyc+timeshift))
                add_sequence!(cam, "cam_partial_after_$(s)_$(l)",  falses(-timeshift))
                add_sequence!(las, "las_partial_after_$(s)_$(l)",  falses(-timeshift))
                @show s
                @show l
            elseif timeshift > 0
                has_partial = true
                warn("Adding an extra piezo cycle in slice $s with lag $l for exposures or flashes that end after the cycle")
                add_sequence!(cam, "cam_partial_before_$(s)_$(l)", falses(timeshift))
                add_sequence!(las, "las_partial_before_$(s)_$(l)", falses(timeshift))
                add_sequence!(cam, "cam_partial_after_$(s)_$(l)", falses(nsamps_cyc-timeshift))
                add_sequence!(las, "las_partial_after_$(s)_$(l)", falses(nsamps_cyc-timeshift))
                @show s
                @show l
            end
            if has_partial
                append!(pos, pos_cyc_name)
                append!(cam, "cam_partial_before_$(s)_$(l)")
                append!(las, "las_partial_before_$(s)_$(l)")
                append!(high_las, "laser_high")
            end
			for t = 1:ntrials
				append!(pos, pos_cyc_name)
				append!(cam, cyc_exp_name)
				append!(las, cyc_flash_name)
				append!(high_las, "laser_high")
			end
            if has_partial
                append!(cam, "cam_partial_after_$(s)_$(l)")
                append!(las, "las_partial_after_$(s)_$(l)")
            end
		end
	end
    if nsamps_offset > 0 #add another positioner cycle because we have too many cam and las samples
        append!(pos, pos_cyc_name)
        append!(cam, "extra_cyc", fill(false, nsamps_cyc-nsamps_offset))
        append!(las, "extra_cyc")
        append!(high_las, "extra_cychigh", fill(true, nsamps_cyc-nsamps_offset))
    end
    @assert length(pos) == length(cam) == length(las) == length(high_las)
    return [pos; cam; las; high_las]
end

function slicetiming_commands(pos::T1, mon_cyc, las_name::AbstractString, cam_name::AbstractString, lags, slice_zs) where {T1<:ImagineSignal}
    mod_cyc = get_samples(pos, first(sequence_names(pos))) #relies on the cycle being the first (and only) subsequence
    @assert all(sequence_names(pos).==first(sequence_names(pos)))
    rig_sigs = rigtemplate(rig_name(pos); sample_rate = samprate(pos))
    cam = getname(rig_sigs, cam_name)
    las = high_las = all_las = -1
    @assert rig_name(pos) == "ocpi-2"
    if las_name != "all lasers"
        las = getname(rig_sigs, "all lasers")
        high_las = getname(rig_sigs, las_name)
        all_las = [las;high_las]
    else
        error("Please specify a specific laser, not all lasers")
    end
    if !ImagineInterface.ispos(pos) error("Only positioner signals are currently supported.") end
    if !ImagineInterface.iscam(cam) error("$cam_name is not a camera trigger signal.") end
    if !ImagineInterface.islas(las) error("$las_name is not a laser trigger signal.") end

    slicetiming_commands!(deepcopy(pos), mod_cyc, mon_cyc, las, high_las, cam, lags, slice_zs)
end


text:
function calculate_flux_variabilty(data_dictionary)

    # --- Phase 0: GET SOME STUFF ---------------------------------------------------#
    # Get the stoichiometric_matrix from data_dictionary -
    stoichiometric_matrix = data_dictionary["stoichiometric_matrix"]
    (number_of_species,number_of_fluxes) = size(stoichiometric_matrix);

    # initialize the output array -
    calculated_flux_array = zeros(number_of_fluxes,3)

    # set some constants -
    TIME_RESTART_LIM = 60
    EXIT_FLAG_SUCCESS = 0
    TOL = 1e-6
    # ------------------------------------------------------------------------------- #

    # --- Phase 1: SOLVE THE ORIGINAL PROBLEM --------------------------------------- #

    # Setup the GLPK problem -
    lp_problem = GLPK.Prob();
    GLPK.set_prob_name(lp_problem, "sample");
    GLPK.set_obj_name(lp_problem, "objective")

    # min -or- max?
    min_flag = data_dictionary["is_minimum_flag"];
    if min_flag == true
    	GLPK.set_obj_dir(lp_problem, GLPK.MIN);
    else
    	GLPK.set_obj_dir(lp_problem, GLPK.MAX);
    end

    # Set the number of constraints and fluxes -
    GLPK.add_rows(lp_problem, number_of_species);
    GLPK.add_cols(lp_problem, number_of_fluxes);

    # Set solver parameters
    solver_parameters = GLPK.SimplexParam();
    solver_parameters.msg_lev = GLPK.MSG_ERR;
    solver_parameters.presolve = GLPK.ON;

    # Setup flux bounds -
    default_bounds_array = data_dictionary["default_flux_bounds_array"]
    (number_of_fluxes,number_of_bounds) = size(default_bounds_array)
    for flux_index = 1:number_of_fluxes

    	flux_lower_bound = default_bounds_array[flux_index,1]
    	flux_upper_bound = default_bounds_array[flux_index,2]

    	# Check bounds type ... default is DB -
    	if (flux_upper_bound == flux_lower_bound)
    		flux_constraint_type = GLPK.FX
    	else
    		flux_constraint_type = GLPK.DB
    	end

    	# flux symbol? (later use name - for now, fake it)
    	flux_symbol = "R_"*string(flux_index)

    	# Set the bounds in GLPK -
    	GLPK.set_col_name(lp_problem, flux_index, flux_symbol);
    	GLPK.set_col_bnds(lp_problem, flux_index, flux_constraint_type, flux_lower_bound, flux_upper_bound);
    end

    # Setup objective function -
    objective_coefficient_array = data_dictionary["objective_coefficient_array"]
    for (flux_index,obj_coeff) in enumerate(objective_coefficient_array)
    	# Set the objective function value in GLPK -
    	GLPK.set_obj_coef(lp_problem, flux_index, obj_coeff);
    end

    # Setup problem constraints for the metabolites -
    species_bounds_array = data_dictionary["species_bounds_array"]
    for species_index = 1:number_of_species
    	species_lower_bound = species_bounds_array[species_index,1]
    	species_upper_bound = species_bounds_array[species_index,2]
    	# defualt
    	species_constraint_type = GLPK.FX
    	if (species_lower_bound != species_upper_bound)
    		species_constraint_type = GLPK.DB
    	end
    	# set the symbol -
    	species_symbol = "x_"*string(species_index)
    	# Set the species bounds in GLPK -
    	GLPK.set_row_name(lp_problem, species_index, species_symbol);
    	GLPK.set_row_bnds(lp_problem, species_index, species_constraint_type, species_lower_bound, species_upper_bound);
    end

    # Setup the stoichiometric array -
    counter = 1;
    row_index_array = zeros(Int,number_of_species*number_of_fluxes);
    col_index_array = zeros(Int,number_of_species*number_of_fluxes);
    species_index_vector = collect(1:number_of_species);
    flux_index_vector = collect(1:number_of_fluxes);
    flat_stoichiometric_array = zeros(Float64,number_of_species*number_of_fluxes);
    for species_index in species_index_vector
    	for flux_index in flux_index_vector
    		row_index_array[counter] = species_index;
    		col_index_array[counter] = flux_index;
    		flat_stoichiometric_array[counter] = stoichiometric_matrix[species_index,flux_index];
    		counter+=1;
    	end
    end

    # Passed in the matrix
    GLPK.load_matrix(lp_problem, number_of_species*number_of_fluxes, row_index_array, col_index_array, flat_stoichiometric_array);

    # Call the solver -
    exit_flag = GLPK.simplex(lp_problem, solver_parameters);

    # Get the objective function value -
    objective_value = GLPK.get_obj_val(lp_problem);

    # Get the calculated flux values for the optimal distrubution from GLPK -
    for flux_index in flux_index_vector
    	calculated_flux_array[flux_index,1] = GLPK.get_col_prim(lp_problem, flux_index);
    end

    # Get the dual values -
    dual_value_array = zeros(Float64,number_of_fluxes);
    for flux_index in flux_index_vector
    	dual_value_array[flux_index] = GLPK.get_col_dual(lp_problem, flux_index);
    end
    # ------------------------------------------------------------------------------- #

    # --- Phase 2: SOLVE THE FVA PROBLEM -------------------------------------------- #
    # retrieve the optimal value
    target_value = 0.0
    if (GLPK.get_obj_dir(lp_problem) == GLPK.MAX)
        target_value = floor(objective_value/TOL)*(TOL);
    else
        target_value = ceil(objective_value/TOL)*(TOL);
    end

    # add the constraint from phase I (the object value) -
    m = GLPK.add_rows(lp_problem,1)
    if (GLPK.get_obj_dir(lp_problem) == GLPK.MAX)
        GLPK.set_row_bnds(lp_problem, m, GLPK.LO,target_value,0.0);
    else
        GLPK.set_row_bnds(lp_problem, m, GLPK.UP,0.0,target_value);
    end
    index_array = Int[]
    value_array = Float64[]
    for flux_index = 1:number_of_fluxes
        push!(index_array,flux_index)
        push!(value_array,GLPK.get_obj_coef(lp_problem, flux_index))
    end
    GLPK.set_mat_row(lp_problem,m,number_of_fluxes,index_array,value_array)

    # zero out the objective coefficients -
    for flux_index = 1:number_of_fluxes
        GLPK.set_obj_coef(lp_problem, flux_index, 0.0);
    end

    # set the parameters
    solver_parameters.presolve = GLPK.OFF;
    solver_parameters.msg_lev = GLPK.MSG_OFF;
    solver_parameters.tm_lim = 1000*TIME_RESTART_LIM;

    # main FVA loop -
    for round_index = 1:2
        # min or maximizing the flux?
        if round_index == 1
            GLPK.set_obj_dir(lp_problem, GLPK.MIN);
        else
            GLPK.set_obj_dir(lp_problem, GLPK.MAX);
        end
        # run for each flux
        for flux_index = 1:number_of_fluxes
            # run the calc -
            GLPK.set_obj_coef(lp_problem, flux_index, 1.0);
            # ??? so by default, it would be fastFVA?
            exit_flag = GLPK.simplex(lp_problem, solver_parameters);
            # try to retrieve results
            if (exit_flag != EXIT_FLAG_SUCCESS)
                solver_parameters.tm_lim = 10000*TIME_RESTART_LIM;
                solver_parameters.presolve = GLPK.ON
                GLPK.adv_basis(lp_problem, 0);  # not fastFVA
                exit_flag = GLPK.simplex(lp_problem, solver_parameters);
                if (exit_flag != EXIT_FLAG_SUCCESS)
                    # ok, something muy malo is happening. throw an error ...
                    error("ERROR: Unable to solve modified problem. Exiting ... ")
                end
                # reset the problem options -
                solver_parameters.presolve = GLPK.OFF;
                solver_parameters.tm_lim = 1000*TIME_RESTART_LIM;
            end
            # reset the objective coefficient back to zero -
            GLPK.set_obj_coef(lp_problem,flux_index,0.0)
            # store the results
            if (GLPK.get_obj_dir(lp_problem) == GLPK.MIN)
                calculated_flux_array[flux_index,2] = GLPK.get_obj_val(lp_problem)
            else
                calculated_flux_array[flux_index,3] = GLPK.get_obj_val(lp_problem)
            end
        end
    end
    # ------------------------------------------------------------------------------- #

    # return the array -
    return (calculated_flux_array,dual_value_array)
end


text:
# Give error for use of `*` as infix operator between tensors

# Remove `*` as infix operator between tensors
function Base.:*(S1::AbstractTensor, S2::AbstractTensor)
    error("Don't use `*` for multiplication between tensors. Use `⋅` (`\\cdot`) for single contraction and `⊡` (`\\boxdot`) for double contraction.")
end

# Remove `'*` as infix operator between tensors
function Base.Ac_mul_B(S1::AbstractTensor, S2::AbstractTensor)
    error("Don't use `A'*B`, use `tdot(A,B)` (or `A'⋅B`) instead.")
end

# Remove `.'*` as infix operator between tensors
function Base.At_mul_B(S1::AbstractTensor, S2::AbstractTensor)
    error("Don't use `A.'*B`, use `tdot(A,B)` (or `A.'⋅B`) instead.")
end

# Remove `\` as infix operator between tensors
function Base.:\(S1::AbstractTensor, S2::AbstractTensor)
    error("Don't use `A\\B`, use `inv(A) ⋅ B` instead.")
end

# Remove + and - between number and Tensor (issue #75)
Base.:+(n::Number, T::AbstractTensor) = throw(MethodError(+, (n, T)))
Base.:+(T::AbstractTensor, n::Number) = throw(MethodError(+, (T, n)))
Base.:-(n::Number, T::AbstractTensor) = throw(MethodError(-, (n, T)))
Base.:-(T::AbstractTensor, n::Number) = throw(MethodError(-, (T, n)))


text:
module DimArrays

export DimArray, DimVector, DimMatrix, dictvector, name!, nest, squeeze

const NameEl = Union{Symbol,Nothing}
const FuncEl = Union{Function,Dict,Nothing}

mutable struct DimArray{T,N,TT} <: AbstractArray{T,N}
    array::TT
    dnames::NTuple{N,NameEl}    # dimension names: symbol or nothing
    ifuncs::NTuple{N,FuncEl}    # index functions: dictionary if bounded
    cname::NameEl               # content name, if any
end

using LinearAlgebra

const DimVector{T} = DimArray{T,1}
const DimMatrix{T} = DimArray{T,2}
const DimVecOrMat{T} = Union{DimVector{T},DimMatrix{T}}
# const DimRowVector{T} = DimArray{T,2,<:RowVector{T}}
const DimAdjoint{T} = DimArray{T,2,<:Adjoint{T}}

### Getters ###
###############

naked(a::DimArray) = a.array
naked(a::AbstractArray) = a

dname(a::DimArray, d::Int) = a.dnames[d]==nothing ? defname(d) : a.dnames[d]
dname(a::AbstractArray, d::Int) = defname(d)
dnames(a::AbstractArray) = [ dname(a,d) for d=1:ndims(a) ]
dname(a::AbstractArray, dims::Tuple) = join([dname(a,d) for d in dims], ", ") # for @info verbose

function defname(d::Int)
    d==1 && return :row
    d==2 && return :col
    d==3 && return :page
    return Symbol(:dim,d)
end

function ifunc(a::DimArray, d::Int)
    ifun = a.ifuncs[d]
    ifun==nothing && return identity
    ifun isa Dict && return n -> haskey(ifun,n) ? ifun[n] : n
    return ifun
end
ifunc(a::AbstractArray, d::Int) = identity
ifuncs(a::AbstractArray) = [ ifunc(a,d) for d=1:ndims(a) ]

label(a::DimArray) = a.cname==nothing ? :content : a.cname
haslabel(a::DimArray) = a.cname != nothing

using Lazy: @forward
@forward DimArray.array  Base.size, Base.length, Base.ndims, Base.getindex, Base.setindex!, Base.append!
@forward DimArray.array  Base.first, Base.iterate, Base.lastindex

### Constructors ###
####################

ensurefuncordict(f::Function) = f
ensurefuncordict(f::Dict) = f # ifunc handles this
ensurefuncordict(f::Number) = f==1 ? identity : (i -> i*f)

using Base.Iterators: repeated

const FuncOrNumberOrDict = Union{Function, Number, Dict}
const SymbolOrString = Union{Symbol, String}

"""
    DimArray(x, names...)
    DimArray(x, names..., maps...)
Interprets symbols & strings as dimension names (in order), and numbers, functions & dictionaries as index maps (for each dimension).
Numbers are converted to functions `i -> i*n` to indicate e.g. that we have every n-th data point.
The `ndims(x)+1`-th name (if given) is equivalent to keyword `label=name`, labelling what the elements of `x` mean.
Giving too few names/maps is no problem, too many will give a warning.
"""
function DimArray(x::TT, tup::Vararg; label = nothing) where TT<:AbstractArray{T,N} where {T,N}
    names = Vector{NameEl}(repeated(nothing, N) |> collect)
    funcs = Vector{FuncEl}(repeated(nothing, N) |> collect)
    n=0; f=0;
    for z in tup
        if isa(z, FuncOrNumberOrDict)
            if f<N
                funcs[f += 1] = ensurefuncordict(z)
            else
                @warn "DimArray is ignoring function $z"
            end
        elseif isa(z, SymbolOrString)
            if n<N
                if z!=""
                    names[n += 1] = Symbol(z)
                else
                    n += 1
                end
            elseif n==N && z!=""
                label=Symbol(z)
                n += 1
            elseif z!=""
                @warn "DimArray is ignoring name $z"
            end
        elseif isa(z, Vector)
            @warn "DimArray is ignoring vector $z, this shouldn't happen!"
        elseif isa(z, Tuple)
            @warn "DimArray is ignoring tuple $z, this really shouldn't happen!"
        elseif z==nothing
            # ignore, no problem!
        end
    end
    DimArray{T,N,TT}(x, NTuple{N,NameEl}(names), NTuple{N,FuncEl}(funcs), label)
end

"""
    DimArray(x, [n1, n2, ...], [f1, f2, ...], n9)
    DimArray(x, [n1, n2, ...], [f1]; label=:n9)
Vectors of names/maps are simply splatted, then digested as above.
Content label `n9` here is kept separate, equivalent to keyword `label=:n9`.
"""
DimArray(x::AbstractArray, list::Vector; kw...) = DimArray(x, list...; kw...);
DimArray(x::AbstractArray, list::Vector, label::Symbol; kw...) = DimArray(x, list...; label=label, kw...);

DimArray(x::AbstractArray, a::Vector, b::Vector; kw...) = DimArray(x, a..., b...; kw...)
DimArray(x::AbstractArray, a::Vector, b::Vector, label::Symbol; kw...) = DimArray(x, a..., b...; label=label, kw...)

DimVector(x::AbstractVector, rest...; kw...) = DimArray(x, rest...; kw...)
DimMatrix(x::AbstractVector, rest...; kw...) = DimArray(reshape(x,:,1), rest...; kw...)
@doc @doc(DimArray)
DimMatrix(x::AbstractMatrix, rest...; kw...) = DimArray(x, rest...; kw...)

"""
    DimArray(T) = DimVector(T[])
    DimArray(sym) = DimVector(Any[], sym)
Empty `DimVector` into to which you can later `push!` things, optionally with a dimension name.
"""
DimArray(T::Type=Any, rest...) = DimArray(T[], rest...)
@doc @doc(DimArray)
DimVector(T::Type=Any, rest...) = DimArray(T[], rest...)

DimArray(s::SymbolOrString, rest...) = DimArray([], s, rest...)
DimVector(s::SymbolOrString, rest...) = DimArray([], s, rest...)

# avoid the above splatting for perfectly formed case, for functions below -- will fail on wrong-length tuples
DimArray(x::TT, dnames::Tuple, ifuncs::Tuple, cname::NameEl) where TT <: AbstractArray{T,N} where {T,N} =
    DimArray{T,N,TT}(x, dnames, ifuncs, cname);

"""
    dictvector(vec, [:first, "second", ...]) = DimVector(vec, Dict(1 => :first, 2 => "second", ...))
    dictvector(vec, [:first, ...], :axis, :content)
Convenient way to define a (short) `DimVector` whose index function is a dictionary labelling the entries.
A name for the dimension / axis, and a label for its content, can be supplied too.
"""
function dictvector(vec::AbstractVector, isym::AbstractVector, rest...)
    ifun = Dict(i => s for (i,s) in enumerate(isym))
    DimVector(vec, ifun, rest...)
end
dictvector(vec::AbstractVector, s::SymbolOrString, isym::AbstractVector, rest...) =
    dictvector(vec, isym, s, rest...)

### One-array operations ###
############################

"""
    name!(x, sym1, sym2)
Name the dimensions of `x::DimArray` to be `sym1, sym2, ...`, with the `ndims(x)+1`-th one interpreted as content label.
"""
function name!(x::DimArray, ss::Vararg{SymbolOrString} ; label = nothing)
    x.dnames = ntuple(d -> d<=length(ss) ? ifelse(ss[d]=="", nothing, Symbol(ss[d])) : x.dnames[d] , ndims(x))
    if label!=nothing
        x.cname = label
        length(ss)>ndims(x) && @warn "name! is ignoring name $(ss[ndims(x)+1])"
    elseif length(ss)>ndims(x)
        x.cname = Symbol(ss[ndims(x)+1])
    end
    length(ss)>ndims(x)+1 && @warn "name! is ignoring name $(ss[ndims(x)+2])"
    x
end

import Base: push!, map, map!

"""
    push!(x, val, sym)
For `x::DimVector` of `length(x)==n`, the result has `x[n+1]=val` with index label `n+1 => sym` added to the dictionary.
Only works if x's index function is `nothing` or a `Dict`.
See also `dictvector([val, val2, val3, ...], [sym, sym2, ...])` for directly creating such things.
"""
push!(a::DimVector, x...) = begin push!(a.array, x...); a end

function push!(a::DimVector, x, s::SymbolOrString)
    push!(a.array, x)
    if a.ifuncs[1]==nothing
        a.ifuncs = (Dict{Int,String}(),)
    end
    if a.ifuncs[1] isa Dict
        a.ifuncs[1][length(a)] = string(s)
    end
    if a.ifuncs[1] isa Function
        @warn "push! can't append index label $s as dimension index is not a Dict" maxlog=2
    end
    a
end

map!(f::Function, a::DimVector) = begin map!(f, a.array); a end
map(f::Function, a::DimVector) = DimArray(map(f, a.array), a.dnames, a.ifuncs, a.cname)

import Base: transpose, adjoint, permutedims

for op in (:transpose, :adjoint)
    @eval begin
        ($op)(a::DimMatrix) = DimArray( ($op)(a.array), rev2(a.dnames), rev2(a.ifuncs), a.cname)
        ($op)(a::DimVector) = DimArray( ($op)(a.array), (:transpose, a.dnames[1]), (nothing, a.ifuncs[1]), a.cname)
        # ($op)(a::DimRowVector) = DimArray( ($op)(a.array), (a.dnames[2],), (a.ifuncs[2],), a.cname)
        ($op)(a::DimAdjoint) = DimArray( ($op)(a.array), (a.dnames[2],), (a.ifuncs[2],), a.cname)
    end # @eval
end

rev2(tup::Tuple) = (tup[2], tup[1])

const TupOrVec = Union{Tuple, AbstractVector}

permutedims(a::DimMatrix, tup::Tuple) = permutedims(a, [tup...])

function permutedims(a::DimArray, p::AbstractVector{Int})
    @assert ndims(a)==length(p) "no valid permutation of dimensions"
    DimArray(permutedims(a.array, p), a.dnames[p], a.ifuncs[p], a.cname)
end

permutedims(a::DimArray, p::AbstractVector{<:SymbolOrString}) = permutedims(a, ensuredim.((a,), p))

permutedims(a::DimMatrix) = permutedims(a, [2,1])

import Base: size, dropdims, selectdim

function size(a::DimArray, d::Union{Int,Symbol})
    d = ensuredim(a,d)
    size(a.array, d)
end

function dropdims(a::DimArray; dims::Union{Int,Symbol,TupOrVec}, verbose=false, zerodim=false)
    d = ensuredim(a,dims)
    d==1 && ndims(a)==1 && !zerodim && return a[1]
    out = DimArray(Base.dropdims(a.array; dims=d), dropd(a.dnames,d), dropd(a.ifuncs,d), a.cname)
    verbose && @info "squeezed along $(dname(a,d)), leaving directions $(join(dnames(out),", ")) size $(size(out))"
    return out
end

squeeze(a::AbstractArray; kw...) = dropdims(a; dims=Tuple(findall(size(a) .== 1)), kw...)

function selectdim(a::DimArray, d::Union{Symbol, Int}, i::Int)
    d = ensuredim(a,d)
    DimArray(selectdim(a.array,d,i), dropd(a.dnames,d), dropd(a.ifuncs,d), a.cname)
end

ensuredim(a::DimArray, d::Int) = d
function ensuredim(a::DimArray{T,N}, s::Symbol) where {T,N}
    d = something(findfirst(isequal(s), dnames(a)))
    1<=d<=N || error("""can't use direction :$s; valid options are $(dnames(a))""")
    return d
end
ensuredim(a::DimArray, dims::TupOrVec) = ensuredim.((a,), dims)

dropd(vec::Vector, d::Int) = d>length(vec) ? vec : append!(vec[1:d-1], vec[d+1:end])
dropd(tup::Tuple, d::Int) = d>length(tup) ? tup : (tup[1:d-1]..., tup[d+1:end]...)
dropd(z, dims::TupOrVec) = length(dims)==1 ? dropd(z, dims[1]) : dropd(dropd(z, dims[end]), dims[1:end-1])

import Base: sum, maximum, minimum
import Statistics: mean, std

for op in (:sum, :mean, :std, :maximum, :minimum )
    @eval begin

        function ($op)(a::DimArray; dims::Union{Int,Symbol,TupOrVec}, squeeze=isa(dims,Symbol), verbose=false, kw...)

            d = ensuredim(a, dims)

            data = ($op)(a.array; dims=d, kw...)
            if squeeze
                out = DimArray(Base.dropdims(data; dims=d), dropd(a.dnames,d), dropd(a.ifuncs,d), a.cname)
            else
                out = DimArray(data, a.dnames, a.ifuncs, a.cname)
            end

            verbose && @info "$(($op))-ed along $(dname(a,dims)), leaving directions $(join(dnames(out),", ")) size $(size(out))"
            return out
        end

    end # @eval
end

### Scalar operations ###
#########################

import Base: +, -, *, /, \

-(a::DimArray) = DimArray(-a.array, a.dnames, a.ifuncs, a.cname)

for op in (:+, :-, :*)
    @eval begin
        ($op)(a::DimArray{T1}, s::T2) where {T1<:Number,T2<:Number} = DimArray(($op)(a.array, s), a.dnames, a.ifuncs, a.cname)
        ($op)(s::T1, a::DimArray{T2}) where {T1<:Number,T2<:Number} = DimArray(($op)(s, a.array), a.dnames, a.ifuncs, a.cname)
    end
end
/(a::DimArray{T1}, s::T2) where {T1<:Number,T2<:Number} = DimArray(a.array / s, a.dnames, a.ifuncs, a.cname)
\(s::T1, a::DimArray{T2}) where {T1<:Number,T2<:Number} = DimArray(s \ a.array, a.dnames, a.ifuncs, a.cname)

### Generators ###
##################

Base.collect(itr::Base.Generator{<:DimArray}) =
    DimArray(collect(Base.Generator(itr.f, itr.iter.array)) , itr.iter.dnames, itr.iter.ifuncs, itr.iter.cname)

### Broadcasting ###
####################

# Base.Broadcast._containertype(::Type{<:DimArray}) = DimArray
#
# Base.Broadcast.promote_containertype(::Type{DimArray}, _) = DimArray
# Base.Broadcast.promote_containertype(_, ::Type{DimArray}) = DimArray
# Base.Broadcast.promote_containertype(::Type{DimArray}, ::Type{Array}) = DimArray
# Base.Broadcast.promote_containertype(::Type{Array}, ::Type{DimArray}) = DimArray
# Base.Broadcast.promote_containertype(::Type{DimArray}, ::Type{DimArray}) = DimArray
#
# array(a::DimArray) = a.array
# array(x) = x
#
# function Base.Broadcast.broadcast_c(fun, t::Type{DimArray}, list...)
#     res = broadcast(fun, array.(list)...)
#     T = eltype(res); N = ndims(res); TT = typeof(res)
#     ## is there a DimArray of the final shape, and not made by v'?
#     for a in list
#         if isa(a, DimArray) && size(a) == size(res) && findfirst(dnames(a), :transpose)==0
#         return DimArray{T,N,TT}(res, a.dnames, a.ifuncs, a.cname) end
#     end
#     ## if not, you have to work harder:
#     rest = reconcile(list, N)
#     DimArray{T,N,TT}(res, rest...)
# end

function reconcile(list, N)
    names = Vector{NameEl}(repeated(nothing, N) |> collect)
    funcs = Vector{FuncEl}(repeated(nothing, N) |> collect)
    label = nothing
    n=0; f=0;
    for a in reverse(list) # overwrite with earliest names
        if isa(a, DimArray)
            for d=1:ndims(a)
                if ifunc(a,d) !== identity # then adopt this one
                    funcs[d] = a.ifuncs[d]
                end
                if a.dnames[d] !== nothing && a.dnames[d] !== :transpose
                    names[d] = a.dnames[d]
                end
            end
            a.cname !== nothing && (label = a.cname)
        end
    end
    for d=2:N
        if names[d] != nothing && findfirst(isequal(names[d]), names[1:d-1])!=nothing # if symbol already used, add a prime...
            for dd=N:-1:d # ... to all later occurances too
                if names[dd]==names[d]
                    names[dd] = Symbol(names[d], "′")
                end
            end
        end
    end
    NTuple{N,NameEl}(names), NTuple{N,FuncEl}(funcs), label
end

### Two-array operations ###
############################

import Base: append!, hcat, vcat, cat

for op in (:+, :-, :append!)
    @eval begin
        # adding a boring array
        ($op)(a::DimArray{T1,N}, b::AbstractArray{T2,N}) where {T1<:Number,T2<:Number,N} =
            DimArray(($op)(a.array, b), a.dnames, a.ifuncs, a.cname)

        ($op)(b::AbstractArray{T1,N}, a::DimArray{T2,N}) where {T1<:Number,T2<:Number,N} =
            DimArray(($op)(b, a.array), a.dnames, a.ifuncs, a.cname)

        # adding two DimArrays
        ($op)(a::DimArray{T1}, b::DimArray{T2}) where {T1<:Number, T2<:Number} =
            DimArray(($op)(a.array, b.array), reconcile((a,b),ndims(a))...)
    end # @eval
end

function hcat(a::DimArray, bb::Vararg{DimArray})
    out = hcat(array(a), array.(bb)...)
    DimArray(out, reconcile((a,bb...), ndims(out))...)
end

function vcat(a::DimArray, bb::Vararg{DimArray})
    out = vcat(array(a), array.(bb)...)
    DimArray(out, reconcile((a,bb...), ndims(out))...)
end

import Base: *, kron

function *(a::DimArray, b::DimArray)
    a.dnames[end]!=b.dnames[1] && a.dnames[end]!=nothing && b.dnames[1]!=nothing &&
        @warn "multiplying along dimensions with mismatched names" maxlog=2
    names = (a.dnames[1:end-1]..., b.dnames[2:end]...)
    funcs = (a.ifuncs[1:end-1]..., b.ifuncs[2:end]...)
    DimArray(a.array * b.array, names, funcs, a.cname )
end

*(a::DimArray, b::AbstractMatrix) = *(a, DimArray(b))
*(a::DimArray, b::AbstractVector) = *(a, DimArray(b))
*(a::AbstractMatrix, b::DimArray) = *(DimArray(a), b)

function kron(a::DimMatrix, b::DimMatrix)
    names = prodsomething.(a.dnames, b.dnames)
    funcs = firstsomething.(a.ifuncs, b.ifuncs)
    DimArray(kron(a.array, b.array), names, funcs, a.cname)
end

function prodsomething(list...)
    out = ""
    for s in list
        s != nothing && (out *= "_"*string(s))
    end
    out == "" && return nothing
    return Symbol(out[2:end])
end

function firstsomething(list...)
    for f in list
        f != nothing && return f
    end
    return nothing
end

### Nested Arrays ###
#####################

"""
    nest([A1, A2, ...])
    nest([A1, A2, ...], sym)
Creates an array with `ndims(A1)+1` dimensions. All the constituent arrays should be similar!
Optional 2nd argument names the new dimension, giving `DimArray(nest(...),sym)`;
this also happens if `A1 isa DimArray`.

    [A1, A2, ...] |> nest(s)
Gives a function which acts on vectors as above.
"""
nest(sym::Symbol) = vec -> nest(vec, sym)

nest(vec::AbstractVector, sym::Symbol) = nest(DimArray(vec,sym))

function nest(vec::AbstractVector) # Allow Vector{Any} but assume contents is all the same Array
    x1 = first(vec)
    isa(x1,AbstractArray) || return vec # nothing to do

    out = similar(x1, size(x1)..., length(vec), )
    for (i, x) in enumerate(vec)
        out[batchindex(out,i)...] = x
    end

    if isa(vec, DimArray) && isa(x1, DimArray)
        return DimArray(out, (x1.dnames..., vec.dnames[1]), (x1.ifuncs..., vec.ifuncs[1]), x1.cname)

    elseif isa(vec, DimArray) # and "x1" is not -- use nothings
        none = repeated(nothing, ndims(x1)) |> collect
        return DimArray(out, (none..., vec.dnames[1]), (none..., vec.ifuncs[1]), nothing)

    elseif isa(x1, DimArray) # and "vec" is not -- use "nest"
        return DimArray(out, (x1.dnames..., :nest), (x1.ifuncs..., nothing), x1.cname)

    else # no DimArrays at all
        return out
    end
end

# that's based on Flux.batch
batchindex(xs, i) = (reverse(Base.tail(reverse(axes(xs))))..., i)

function nest(arr::AbstractArray{D}) where {D}
    one = [nest(slicedim(arr, ndims(arr), c)) for c=1:size(arr,ndims(arr))]
    nest(one)
end

### Display ###
###############

function Base.show(io::IO, a::DimArray)
    print(io, "DimArray(")
    show(io, a.array)
    for d=1:ndims(a)
        a.dnames[d]!= nothing && print(io, ", :",dname(a,d))
    end
    haslabel(a) && print(io, "; label = :", label(a))
    print(io, ")")
end

function Base.summary(io::IO, a::DimArray)
    lab = haslabel(a) ? ", label = "*string(label(a))*"," : ""
    if typeof(a.array) <: Array
        ndims(a)==0 && println(io, "DimArray{",eltype(a),"}$lab of zero dimensions:")
        ndims(a)==1 && println(io, length(a),"-element DimVector{",eltype(a),"}$lab with dimension:")
        ndims(a)==2 && println(io, size(a,1),"×",size(a,2)," DimMatrix{",eltype(a),"}$lab with dimensions:")
        ndims(a)>=3 && println(io, join(string.(size(a)),"×"), " DimArray{",eltype(a),",",ndims(a),"}$lab with dimensions:")
    # elseif typeof(a.array) <: RowVector
    #     println(io, length(a),"-element DimRowVector{",eltype(a),"} with dimensions:")
    elseif typeof(a.array) <: Adjoint
        println(io, length(a),"-element DimAdjoint{",eltype(a),"} with dimensions:")
    else
        println(io, summary(a.array), " wrapped in a DimArray with:")
    end
    for d=1:ndims(a)
        if a.dnames[d]==nothing #|| a.dnames[d]==:transpose
                    print(io, "   ⭒ ") # different symbol for default names
        elseif d==1 print(io, "   ⬙ ") # ⇁⇂
        elseif d==2 print(io, "   ⬗ ")
        else        print(io, "   ◇ ")
        end
        print(io, rpad(string(dname(a,d)),4))
        if length(a)>0
            print(io, " = ", stringfew(size(a,d)), " ")
        end
        if ifunc(a,d)!=identity
            print(io, " ⟹   ", stringfew(size(a,d), ifunc(a,d)), " ")
        end
        if d<ndims(a)
            println(io, "")
        end
    end
end

function stringfew(n, f=identity)
    n<=5 && return join(stringone.(collect(1:n),f), ", ")
    return join(stringone.([1,2],f), ", ") * ", … " * join(stringone.([n-1,n],f), ", ")
end

function stringone(n, f=identity)
    f(n) isa Integer && return string(f(n))
    f(n) isa Number && return string(round(f(n), digits=3))
    return string(f(n))
end

### Plotting ###
################

using RecipesBase

@recipe function ff(a::DimVector)
    label --> ""

    # first dim is x axis
    if a.dnames[1] != nothing
        xaxis --> string(dname(a,1))
    end
    xlist = ifunc(a,1).(1:length(a))

    xlist, a.array
end

@recipe function ff(a::DimMatrix)

    # first dim is x axis
    xaxis --> string(dname(a,1))
    xlist = ifunc(a,1).(1:size(a,1))

    # second dim labels series
    ylist = string(dname(a,2))*" = " .* string.(ifunc(a,2).(1:size(a,2)))
    label --> reshape(ylist, 1,:) # because series labels must be cols of a matrix, not el of vector

    if haslabel(a)
        yaxis --> string(a.cname)
    end

    xlist, a.array
end

### Conversions ###
###################

using Requires

function __init__()

    @require NamedArrays = "86f7a689-2022-50b4-a561-43c23ac3c673" begin
        using NamedArrays: NamedArray

        axestuple(a::DimArray) = Tuple([ifunc(a,d).(collect(1:size(a,d))) for d=1:ndims(a)])
        nametuple(a::DimArray) = Tuple(dnames(a))
        Base.convert(::Type{NamedArray}, a::DimArray) = NamedArray(a.array, axestuple(a), nametuple(a))
        NamedArray(a::DimArray) = convert(NamedArray, a)

    end

    @require AxisArrays = "39de3d68-74b9-583c-8d2d-e117c070f3a9" begin
        using AxisArrays: AxisArray, Axis

        axislist(a::DimArray) = [ Axis{dname(a,d)}( ifunc(a,d).(1:size(a,d) ) ) for d=1:ndims(a) ]
        Base.convert(::Type{AxisArray}, a::DimArray) = AxisArray(a.array, axislist(a)...)
        AxisArray(a::DimArray) = convert(AxisArray, a)

    end

end

### End ###
###########

end # module


text:
using Lehmann
using DelimitedFiles
using Printf
# using Gaston  # for ploting

β = 10.0 # inverse temperature
Euv = 4.0 # UV energy cutoff for the spectral density
eps = 1.0e-8 # accuracy for DLR basis

data = DelimitedFiles.readdlm("yuanyao.txt")
τ = data[:, 1] * β
τ[1] += 1.0e-8
τ[end] -= 1.0e-8
y = data[:, 2]
err = data[:, 3]
w = 1.0 ./ err.^2

# construct DLR basis
dlr = DLR.DLRGrid(:fermi, Euv, β, eps)

kernel = Spectral.kernelT(:fermi, τ, dlr.ω, β)
# y' = kernel \cdot coeff
# use linear fit to find the optimal coeff so that chi2 = sum_i w_i * (y_i-y'_i)^2 is minimized, where w_i = 1/errorbar_i^2
coeff = kernel \ y # linear fit without reweight, this gives an estimate of the initial values for the coefficients

######################### G_dlr(τ) ################################
yp = kernel * coeff
# print comparsion for selected τ point
printstyled("Compare G(τ) <-->  G_dlr(τ)\n", color=:green)
printstyled("G(0⁺)+G(β⁻): ", y[1] + y[end], " vs ", yp[1] + yp[end], "\n", color=:red)
printstyled("maximum dlr fitting error: ", maximum(abs.(y - yp)), "\n", color=:red)
printstyled(" τi       G(τ)        G_dlr(τ)      δG\n", color=:yellow)

for i in 1:25:length(y)
	@printf("%3i    %10.8f    %10.8f   %10.8f\n", i, y[i], yp[i], abs(y[i] - yp[i]))
end

######################### G_dlr(iωn) ################################
n = collect(0:100)
kernelΩ = Spectral.kernelΩ(:fermi, n, dlr.ω, β)
Gn = kernelΩ * coeff

printstyled("Matsubara frequency G_dlr(ωn)\n", color=:green)
printstyled(" n          Re           Im\n", color=:yellow)
for i in 1:length(n)
	@printf("%3i    %10.8f   %10.8f\n", n[i], real(Gn[i]), -imag(Gn[i]))
end

######################### Σ_dlr(iωn) ################################
g0inv = @. -1im * (2 * n + 1) * π / β + 0.4
Ginv = 1.0 ./ Gn
Σ = g0inv - Ginv

printstyled("Matsubara frequency Σ_dlr(ωn)\n", color=:green)
printstyled(" n          Re           Im\n", color=:yellow)
for i in 1:length(n)
	@printf("%3i    %10.8f   %10.8f\n", n[i], real(Σ[i]), -imag(Σ[i]))
end

# set(term="qt")
# p = plot(n, -imag.(Σ))
# display(p)
# readline()


text:
function fib(n::Int)
    if n < 2
        return n
    else
        return fib(n-2) + fib(n-1)
    end
end

#  @time fib(40)


text:
# dummy allocator that passes through any requests, calling into the GC if that fails.

Base.@kwdef struct NoPool <: AbstractPool
    stream_ordered::Bool
end

function alloc(pool::NoPool, sz; stream::CuStream)
    block = nothing
    for phase in 1:3
        if phase == 2
            @pool_timeit "$phase.0 gc (incremental)" GC.gc(false)
        elseif phase == 3
            @pool_timeit "$phase.0 gc (full)" GC.gc(true)
        end

        @pool_timeit "$phase.1 alloc" begin
            block = actual_alloc(sz, phase==3; pool.stream_ordered, stream)
        end
        block === nothing || break
    end

    return block
end

function free(pool::NoPool, block; stream::CuStream)
    actual_free(block; pool.stream_ordered, stream)
    return
end

reclaim(pool::NoPool, target_bytes::Int=typemax(Int)) = return 0

cached_memory(pool::NoPool) = 0


text:
# W32API.jl
# deps WCharUTF8

VERSION >= v"0.4.0-dev+6521" && __precompile__()
module W32API

import Base
import WCharUTF8
export MBox, MBoxA, MBoxW

function MBox(Msg::AbstractString, Title::AbstractString) # expect UTF8String
  return MBoxW(WCharUTF8.UTF8toWCS(Msg), WCharUTF8.UTF8toWCS(Title))
end

function MBoxA(Msg::AbstractString, Title::AbstractString) # expect ASCIIString
  user32 = Base.Libdl.dlopen_e("user32.dll")
  MessageBoxA = Base.Libdl.dlsym(user32, :MessageBoxA)
  err = ccall(MessageBoxA, stdcall,
    UInt, (Ptr{Void}, Ptr{UInt8}, Ptr{UInt8}, UInt),
    C_NULL, [Msg.data; [0x00]], [Title.data; [0x00]], 0)
  Base.Libdl.dlclose(user32)
  return err
end

function MBoxW(Msg::AbstractString, Title::AbstractString) # expect UTF8String
  return MBoxW(WCharUTF8.UTF8toWCS(Msg), WCharUTF8.UTF8toWCS(Title))
end

function MBoxW(Msg::Array{UInt16, 1}, Title::Array{UInt16, 1}) # expect wchar_t
  user32 = Base.Libdl.dlopen_e("user32.dll")
  MessageBoxW = Base.Libdl.dlsym(user32, :MessageBoxW)
  err = ccall(MessageBoxW, stdcall,
    UInt, (Ptr{Void}, Ptr{Cwchar_t}, Ptr{Cwchar_t}, UInt),
    C_NULL, [Msg; [0x0000]], [Title; [0x0000]], 0)
  Base.Libdl.dlclose(user32)
  return err
end

end