text:
"""
A `GridGenerator` represents a family of grids.

The `generate` function maps a parameter value to a grid.
"""
abstract type GridGenerator{T} end

parameter(grid::GridGenerator) = 1

struct ModelGridGenerator{T} <: GridGenerator{T}
    modelgrid   ::  AbstractGrid{T}
end

parameter(gen::ModelGridGenerator) = dimensions(gen.modelgrid)

generate(gen::ModelGridGenerator, param) = resize(gen.modelgrid, param)

struct SubGridGenerator{T} <: GridGenerator{T}
    domain      ::  Domain{T}
    generator   ::  GridGenerator{T}
end

parameter(gen::SubGridGenerator) = parameter(gen.generator)

generate(gen::SubGridGenerator, param) = subgrid(generate(gen.generator, param), gen.domain)

struct EquispacedGrids{T} <: GridGenerator{T}
    domain  ::  Domain{T}
end

equispaced_grid(box::AbstractInterval{T}, n) where {T} =
    EquispacedGrid{T}(n, infimum(box), supremum(box))

equispaced_grid(box::ProductDomain, n) =
    ProductGrid(equispaced_grid.(elements(box), n)...)

generate(gen::EquispacedGrids, param) = equispaced_grid(gen.domain, param)

struct HexagonalGrids{T} <: GridGenerator{SVector{2,T}}
    bounding_box   ::  Vector{T}
end

function hexagonal_grid(N, bounding_box)
    Rad3Over2 = sqrt(3) / 2
    x = 0:1:round(sqrt(N)-1)
    X = [c for d in x, c in x]
    Y = [d for d in x, c in x]
    n = length(x)
    X = Rad3Over2 * X
    X /= maximum(X)
    Y[:,2:2:end] .+= 0.5
    Y /= maximum(Y)
    Lx = bounding_box[2]-bounding_box[1]
    Ly = bounding_box[4]-bounding_box[3]
    X = X[:] * Lx .+ bounding_box[1]
    Y = Y[:] * Ly .+ bounding_box[3]
    G = map(SVector, X, Y)
end

generate(gen::HexagonalGrids, param::Int) = hexagonal_grid(param, gen.bounding_box)


text:
# date
s = "2012-12-12"
d = Date(s)
r = RObject(d)
@test isa(r,RObject{RealSxp})
@test rcopy(getclass(r)) == "Date"
@test length(r) == 1
@test size(r) == (1,)
@test rcopy(Date, r) === d
@test rcopy(Date, R"as.Date($s)") == d
@test rcopy(R"identical(as.Date($s), $d)")

s = ["2001-01-01", "1111-11-11", "2012-12-12"]
d = Date.(s)
r = RObject(d)
@test isa(r,RObject{RealSxp})
@test rcopy(getclass(r)) == "Date"
@test length(r) == length(d)
@test size(r) == size(d)
@test rcopy(r) == d
@test rcopy(Array, r) == d
@test rcopy(R"as.Date($s)") == d
@test rcopy(R"identical(as.Date($s), $d)")

# SubArray
d = view(d, 1:2)
r = RObject(d)
@test isa(r,RObject{RealSxp})
@test rcopy(getclass(r)) == "Date"
@test length(r) == length(d)
@test size(r) == size(d)
@test rcopy(r) == d
@test rcopy(Array, r) == d
@test rcopy(R"as.Date($s[1:2])") == d
@test rcopy(R"identical(as.Date($s[1:2]), $d)")

d = Date[]
r = RObject(d)
@test isa(r,RObject{RealSxp})
@test rcopy(getclass(r)) == "Date"
@test isa(rcopy(Array, r), Array{Date})
@test length(r) == length(d)
@test size(r) == size(d)
@test rcopy(r) == d
@test rcopy(R"as.Date(character(0))") == Date[]

# Missing array date
d = [Date("2001-01-01"), missing, Date("2012-12-12")]
r = RObject(d)
@test isa(r,RObject{RealSxp})
@test rcopy(getclass(r)) == "Date"
@test length(r) == length(d)
@test size(r) == size(d)
@test rcopy(r)[[1,3]] == d[[1,3]]
@test ismissing(rcopy(r)[2])
@test rcopy(Array, r)[[1,3]] == d[[1,3]]
@test ismissing(rcopy(Array, r)[2])

# dateTime
s = "2012-12-12T12:12:12"
d = DateTime(s)
r = RObject(d)
@test isa(r,RObject{RealSxp})
@test rcopy(getclass(r)) == ["POSIXct", "POSIXt"]
@test rcopy(getattrib(r, "tzone")) == "UTC"
@test length(r) == 1
@test size(r) == (1,)
@test rcopy(r) === d
@test rcopy(DateTime, r) === d
@test rcopy(R"as.POSIXct($s, 'UTC', '%Y-%m-%dT%H:%M:%S')") == d
@test rcopy(R"identical(as.character($d, '%Y-%m-%dT%H:%M:%S'), $s)")

s = ["2001-01-01T01:01:01", "1111-11-11T11:11:00", "2012-12-12T12:12:12"]
d = DateTime.(s)
r = RObject(d)
@test isa(r,RObject{RealSxp})
@test rcopy(getclass(r)) == ["POSIXct", "POSIXt"]
@test rcopy(getattrib(r, "tzone")) == "UTC"
@test isa(rcopy(Array, r), Array{DateTime})
@test length(r) == length(d)
@test size(r) == size(d)
@test rcopy(r) == d
@test rcopy(Array{DateTime}, r) == d
@test rcopy(R"as.POSIXct($s, 'UTC', '%Y-%m-%dT%H:%M:%S')") == d
@test rcopy(R"identical(as.character($d, '%Y-%m-%dT%H:%M:%S'), $s)")

# SubArray
d = view(d, 1:2)
r = RObject(d)
@test isa(r,RObject{RealSxp})
@test rcopy(getclass(r)) == ["POSIXct", "POSIXt"]
@test rcopy(getattrib(r, "tzone")) == "UTC"
@test isa(rcopy(Array, r), Array{DateTime})
@test length(r) == length(d)
@test size(r) == size(d)
@test rcopy(r) == d
@test rcopy(Array{DateTime}, r) == d
@test rcopy(R"as.POSIXct($s[1:2], 'UTC', '%Y-%m-%dT%H:%M:%S')") == d
@test rcopy(R"identical(as.character($d, '%Y-%m-%dT%H:%M:%S'), $s[1:2])")

d = DateTime[]
r = RObject(d)
@test isa(r,RObject{RealSxp})
@test rcopy(getclass(r)) == ["POSIXct", "POSIXt"]
@test rcopy(getattrib(r, "tzone")) == "UTC"
@test isa(rcopy(Array, r), Array{DateTime})
@test length(r) == length(d)
@test size(r) == size(d)
@test rcopy(r) == d
@test rcopy(R"as.POSIXct(character(0))") == Date[]

# Missing array dateTime
d = [DateTime("2001-01-01T01:01:01"), missing, DateTime("2012-12-12T12:12:12")]
r = RObject(d)
@test isa(r,RObject{RealSxp})
@test rcopy(getclass(r)) == ["POSIXct", "POSIXt"]
@test length(r) == length(d)
@test size(r) == size(d)
@test rcopy(r)[[1,3]] == d[[1,3]]
@test ismissing(rcopy(r)[2])
@test rcopy(Array, r)[[1,3]] == d[[1,3]]
@test ismissing(rcopy(Array, r)[2])

ext: jl


text:
# -------------------------- #
# Standard Julia API methods #
# -------------------------- #
const ROW_COL_TYPE = Union{String,Int,Vector{Int},Colon}

_check_row_col_arg(_l::Layout, x::String, name, dim::Int=1) = x == "all" || @error "Unknown specifier $x for $name, must be an integer or `\"all\"`"
_check_row_col_arg(l::Layout, x::Int, name, dim::Int=1) = x <= size(l.subplots.grid_ref, dim) || error("$name index $x is too large for layout")
_check_row_col_arg(l::Layout, x::Vector{Int}, name, dim::Int=1) = foreach(_val -> _check_row_col_arg(l, _val, name, dim), x)
_check_row_col_arg(_l::Layout, x::Colon, name, dim::Int=1) = nothing

prep_kwarg(pair::Union{Pair,Tuple}) =
    (Symbol(replace(string(pair[1]), "_" => ".")), pair[2])
prep_kwargs(pairs::AbstractVector) = Dict(map(prep_kwarg, pairs))
prep_kwargs(pairs::AbstractDict) = Dict(prep_kwarg((k, v)) for (k, v) in pairs)

"""
    size(::PlotlyBase.Plot)

Return the size of the plot in pixels. Obtained from the `layout.width` and
`layout.height` fields.
"""
Base.size(p::Plot) = (get(p.layout.fields, :width, 800),
                      get(p.layout.fields, :height, 450))

const _TRACE_TYPES = [
    :bar, :barpolar, :box, :candlestick, :carpet, :choropleth,
    :choroplethmapbox, :cone, :contour, :contourcarpet, :densitymapbox,
    :funnel, :funnelarea, :heatmap, :heatmapgl, :histogram, :histogram2d,
    :histogram2dcontour, :icicle, :image, :indicator, :isosurface, :mesh3d, :ohlc,
    :parcats, :parcoords, :pie, :pointcloud, :sankey, :scatter, :scatter3d,
    :scattercarpet, :scattergeo, :scattergl, :scattermapbox, :scatterpolar,
    :scatterpolargl, :scatterternary, :splom, :streamtube, :sunburst,
    :surface, :table, :treemap, :violin, :volume, :waterfall
]

for t in _TRACE_TYPES
    str_t = string(t)
    code = quote
        $t(;kwargs...) = GenericTrace($str_t; kwargs...)
        $t(d::AbstractDict; kwargs...) = GenericTrace($str_t, _symbol_dict(d); kwargs...)
    end
    eval(code)
    eval(Expr(:export, t))
end

Base.copy(hf::HF) where {HF <: HasFields} = HF(deepcopy(hf.fields))
Base.copy(p::Plot) = Plot(AbstractTrace[copy(t) for t in p.data], copy(p.layout))
fork(p::Plot) = Plot(deepcopy(p.data), copy(p.layout))

# -------------- #
# Javascript API #
# -------------- #

#=

this function is internal and allows us to match plotly.js semantics in
`resytle!`. The reason is that if you try to set an attribute on a trace with
and array, Plotly.restyle expects an array of arrays.

This means that to set the :x field with [1,2,3], the json should look
something like `[[1,2,3]]`. In Julia we can get this with a one row matrix `[1
2 3]'` or a tuple of arrays `([1, 2, 3], )`. This function applies that logic
and extracts the first element from an array or a tuple before calling
setindex.

All other field types are let through directly

NOTE that the argument `i` here is _not_ the same as the argument `ind` below.
`i` tracks which index we should use when extracting an element from
`v::Union{AbstractArray,Tuple}` whereas `ind` below specifies which trace to
apply the update to. =#
function _apply_restyle_setindex!(hf::Union{AbstractDict,HasFields}, k::Symbol,
                                  v::Union{AbstractArray,Tuple}, i::Int)
    setindex!(hf, v[i], k)
end

_apply_restyle_setindex!(hf::Union{AbstractDict,HasFields}, k::Symbol, v, i::Int) =
    setindex!(hf, v, k)

#=
Wrap the vector so it repeats to be at least length N

This means

```julia
_prep_restyle_vec_setindex([1, 2], 2) --> [1, 2]
_prep_restyle_vec_setindex([1, 2], 3) --> [1, 2, 1]
_prep_restyle_vec_setindex([1, 2], 4) --> [1, 2, 1, 2]

_prep_restyle_vec_setindex((1, [42, 4]), 2) --> [1, [42, 4]]
_prep_restyle_vec_setindex((1, [42, 4]), 3) --> [1, [42, 4], 1]
_prep_restyle_vec_setindex((1, [42, 4]), 4) --> [1, [42, 4], 1, [42, 4]]
``` =#
_prep_restyle_vec_setindex(v::AbstractVector, N::Int) =
    repeat(v, outer=[ceil(Int, N / length(v))])[1:N]

# treat tuples like vectors, just like JSON.json does
_prep_restyle_vec_setindex(v::Tuple, N::Int) =
    _prep_restyle_vec_setindex(Any[i for i in v], N)

# everything else just goes through
_prep_restyle_vec_setindex(v, N::Int) = v

function _update_fields(hf::GenericTrace, i::Int, update::Dict=Dict(); kwargs...)
    # apply updates in the dict w/out `_` processing
    for (k, v) in update
        _apply_restyle_setindex!(hf.fields, k, v, i)
    end
    for (k, v) in kwargs
        _apply_restyle_setindex!(hf, k, v, i)
    end
    hf
end

"""
    relayout!(l::Layout, update::AbstractDict=Dict(); kwargs...)

Update `l` using update dict and/or kwargs
"""
function relayout!(l::Layout, update::AbstractDict=Dict(); kwargs...)
    merge!(l.fields, update)  # apply updates in the dict w/out `_` processing
    merge!(l, Layout(;kwargs...))
    l
end

function relayout!(dest::Layout, src::Layout; kwargs...)
    merge!(dest, src)
    relayout!(dest; kwargs...)
end

"""
    relayout!(p::Plot, update::AbstractDict=Dict(); kwargs...)

Update `p.layout` on using update dict and/or kwargs
"""
relayout!(p::Plot, args...; kwargs...) =
    relayout!(p.layout, args...; kwargs...)

"""
    restyle!(gt::GenericTrace, i::Int=1, update::AbstractDict=Dict(); kwargs...)

Update trace `gt` using dict/kwargs, assuming it was the `i`th ind in a call
to `restyle!(::Plot, ...)`
"""
restyle!(gt::GenericTrace, i::Int=1, update::AbstractDict=Dict(); kwargs...) =
    _update_fields(gt, i, update; kwargs...)

"""
    restyle!(p::Plot, ind::Int=1, update::AbstractDict=Dict(); kwargs...)

Update `p.data[ind]` using update dict and/or kwargs
"""
restyle!(p::Plot, ind::Int, update::AbstractDict=Dict(); kwargs...) =
    restyle!(p.data[ind], 1, update; kwargs...)

"""
    restyle!(::Plot, ::AbstractVector{Int}, ::AbstractDict=Dict(); kwargs...)

Update specific traces at `p.data[inds]` using update dict and/or kwargs
"""
function restyle!(p::Plot, inds::AbstractVector{Int},
                  update::AbstractDict=Dict(); kwargs...)
    N = length(inds)
    kw = Dict{Symbol,Any}(kwargs)

    # prepare update and kw dicts for vectorized application
    for d in (kw, update)
        for (k, v) in d
            d[k] = _prep_restyle_vec_setindex(v, N)
        end
    end

    map((ind, i) -> restyle!(p.data[ind], i, update; kw...), inds, 1:N)
end

"""
    restyle!(p::Plot, update::AbstractDict=Dict(); kwargs...)

Update all traces using update dict and/or kwargs
"""
restyle!(p::Plot, update::AbstractDict=Dict(); kwargs...) =
    restyle!(p, 1:length(p.data), update; kwargs...)

"""
The `restyle!` method follows the semantics of the `Plotly.restyle` function in
plotly.js. Specifically the following rules are applied when trying to set
an attribute `k` to a value `v` on trace `ind`, which happens to be the `i`th
trace listed in the vector of `ind`s (if `ind` is a scalar then `i` is always
equal to 1)

- if `v` is an array or a tuple (both translated to javascript arrays when
`json(v)` is called) then `p.data[ind][k]` will be set to `v[i]`. See examples
below
- if `v` is any other type (any scalar type), then `k` is set directly to `v`.

**Examples**

```julia
# set marker color on first two traces to be red
restyle!(p, [1, 2], marker_color="red")

# set marker color on trace 1 to be green and trace 2 to be red
restyle!(p, [2, 1], marker_color=["red", "green"])

# set marker color on trace 1 to be red. green is not used
restyle!(p, 1, marker_color=["red", "green"])

# set the first marker on trace 1 to red, the second marker on trace 1 to green
restyle!(p, 1, marker_color=(["red", "green"],))

# suppose p has 3 traces.
# sets marker color on trace 1 to ["red", "green"]
# sets marker color on trace 2 to "blue"
# sets marker color on trace 3 to ["red", "green"]
restyle!(p, 1:3, marker_color=(["red", "green"], "blue"))
```
"""
restyle!

function update!(
        p::Plot, ind::Union{AbstractVector{Int},Int},
        update::AbstractDict=Dict(); layout::Layout=Layout(),
        kwargs...
    )
    relayout!(p; layout.fields...)
    restyle!(p, ind, update; kwargs...)
    p
end

function update!(p::Plot, update=Dict(); layout::Layout=Layout(), kwargs...)
    update!(p, 1:length(p.data), update; layout=layout, kwargs...)
end

"""
Apply both `restyle!` and `relayout!` to the plot. Layout arguments are
specified by passing an instance of `Layout` to the `layout` keyword argument.

The `update` Dict (optional) and all keyword arguments will be passed to
restyle

## Example

```jlcon
julia> p = Plot([scatter(y=[1, 2, 3])], Layout(yaxis_title="this is y"));

julia> print(json(p, 2))
{
  "layout": {
    "margin": {
      "l": 50,
      "b": 50,
      "r": 50,
      "t": 60
    },
    "yaxis": {
      "title": "this is y"
    }
  },
  "data": [
    {
      "y": [
        1,
        2,
        3
      ],
      "type": "scatter"
    }
  ]
}

julia> update!(p, Dict(:marker => Dict(:color => "red")), layout=Layout(title="this is a title"), marker_symbol="star");

julia> print(json(p, 2))
{
  "layout": {
    "margin": {
      "l": 50,
      "b": 50,
      "r": 50,
      "t": 60
    },
    "yaxis": {
      "title": "this is y"
    },
    "title": "this is a title"
  },
  "data": [
    {
      "y": [
        1,
        2,
        3
      ],
      "type": "scatter",
      "marker": {
        "color": "red",
        "symbol": "star"
      }
    }
  ]
}
```
"""
update!

"""
    addtraces!(p::Plot, traces::AbstractTrace...)

Add trace(s) to the end of the Plot's array of data
"""
addtraces!(p::Plot, traces::AbstractTrace...) = push!(p.data, traces...)

"""
    addtraces!(p::Plot, i::Int, traces::AbstractTrace...)

Add trace(s) at a specified location in the Plot's array of data.

The new traces will start at index `p.data[i]`
"""
function addtraces!(p::Plot, i::Int, traces::AbstractTrace...)
    new_data = vcat(p.data[1:i - 1], traces..., p.data[i:end])
    p.data = new_data
end

"""
    deletetraces!(p::Plot, inds::Int...) =

Remove the traces at the specified indices
"""
deletetraces!(p::Plot, inds::Int...) =
    (p.data = p.data[setdiff(1:length(p.data), inds)])

"""
    movetraces!(p::Plot, to_end::Int...)

Move one or more traces to the end of the data array"
"""
movetraces!(p::Plot, to_end::Int...) =
    (p.data = p.data[vcat(setdiff(1:length(p.data), to_end), to_end...)])

function _move_one!(x::AbstractArray, from::Int, to::Int)
    el = splice!(x, from)  # extract the element
    splice!(x, to:to - 1, (el,))  # put it back in the new position
    x
end

"""
    movetraces!(p::Plot, src::AbstractVector{Int}, dest::AbstractVector{Int})

Move traces from indices `src` to indices `dest`.

Both `src` and `dest` must be `Vector{Int}`
"""
movetraces!(p::Plot, src::AbstractVector{Int}, dest::AbstractVector{Int}) =
    (map((i, j) -> _move_one!(p.data, i, j), src, dest); p)

function purge!(p::Plot)
    empty!(p.data)
    p.layout = Layout()
    nothing
end

function react!(p::Plot, data::AbstractVector{<:AbstractTrace}, layout::Layout)
    p.data = data
    p.layout = layout
    nothing
end

# no-op here
redraw!(p::Plot) = nothing
to_image(p::Plot; kwargs...) = nothing
download_image(p::Plot; kwargs...) = nothing

_tovec(v) = _tovec([v])
_tovec(v::Vector) = eltype(v) <: Vector ? v : Vector[v]

"""
    extendtraces!(::Plot, ::Dict{Union{Symbol,AbstractString},AbstractVector{Vector{Any}}}), indices, maxpoints)

Extend one or more traces with more data. A few notes about the structure of the
update dict are important to remember:

- The keys of the dict should be of type `Symbol` or `AbstractString` specifying
  the trace attribute to be updated. These attributes must already exist in the
  trace
- The values of the dict _must be_ a `Vector` of `Vector` of data. The outer index
  tells Plotly which trace to update, whereas the `Vector` at that index contains
  the value to be appended to the trace attribute.

These concepts are best understood by example:

```julia
# adds the values [1, 3] to the end of the first trace's y attribute and doesn't
# remove any points
extendtraces!(p, Dict(:y=>Vector[[1, 3]]), [1], -1)
extendtraces!(p, Dict(:y=>Vector[[1, 3]]))  # equivalent to above
```

```julia
# adds the values [1, 3] to the end of the third trace's marker.size attribute
# and [5,5,6] to the end of the 5th traces marker.size -- leaving at most 10
# points per marker.size attribute
extendtraces!(p, Dict("marker.size"=>Vector[[1, 3], [5, 5, 6]]), [3, 5], 10)
```

"""
function extendtraces!(p::Plot, update::AbstractDict, indices::AbstractVector{Int}=[1],
                       maxpoints=-1)
    # TODO: maxpoints not handled here
    for (ix, p_ix) in enumerate(indices)
        tr = p.data[p_ix]
        for k in keys(update)
            v = update[k][ix]
            tr[k] = push!(tr[k], v...)
        end
    end
end

"""
    prependtraces!(p::Plot, update::AbstractDict, indices::AbstractVector{Int}=[1],
                    maxpoints=-1)

The API for `prependtraces` is equivalent to that for `extendtraces` except that
the data is added to the front of the traces attributes instead of the end. See
Those docstrings for more information
"""
function prependtraces!(p::Plot, update::AbstractDict, indices::AbstractVector{Int}=[1],
                        maxpoints=-1)
    # TODO: maxpoints not handled here
    for (ix, p_ix) in enumerate(indices)
        tr = p.data[p_ix]
        for k in keys(update)
            v = update[k][ix]
            tr[k] = vcat(v, tr[k])
        end
    end
end

for f in (:extendtraces!, :prependtraces!)
    @eval begin
        $(f)(p::Plot, inds::Vector{Int}=[1], maxpoints=-1; update...) =
            ($f)(p, Dict(map(x -> (x[1], _tovec(x[2])), update)), inds, maxpoints)

        $(f)(p::Plot, ind::Int, maxpoints=-1; update...) =
            ($f)(p, [ind], maxpoints; update...)

        $(f)(p::Plot, update::AbstractDict, ind::Int, maxpoints=-1) =
            ($f)(p, update, [ind], maxpoints)
    end
end

for f in [:restyle, :relayout, :update, :addtraces, :deletetraces,
          :movetraces, :redraw, :extendtraces, :prependtraces, :purge, :react]
    f! = Symbol(f, "!")
    @eval function $(f)(p::Plot, args...; kwargs...)
        out = fork(p)
        $(f!)(out, args...; kwargs...)
        out
    end
end

function _add_trace!(p::Plot, trace::GenericTrace, row::Int, col::Int, secondary_y::Bool)
    refs = p.layout.subplots.grid_ref[row, col]
    if secondary_y && length(refs) == 1
        msg = "To use secondary_y, you must have created the Subplot with seconary_y=true"
        error(msg)
    end

    ref = refs[secondary_y ? 2 : 1]
    ref_trace = deepcopy(trace)
    merge!(ref_trace, ref.trace_kwargs)
    push!(p.data, ref_trace)
    p
end

function _add_trace!(p::Plot, trace::GenericTrace, row::ROW_COL_TYPE, col::ROW_COL_TYPE, secondary_y::Bool)
    for refs in p.layout.subplots.grid_ref[row, col]
        if secondary_y && length(refs) == 1
            msg = "To use secondary_y, you must have created the Subplot with seconary_y=true"
            error(msg)
        end
        ref = refs[secondary_y ? 2 : 1]
        ref_trace = deepcopy(trace)
        merge!(ref_trace, ref.trace_kwargs)
        push!(p.data, ref_trace)
    end
    p
end

function add_trace!(p::Plot, trace::GenericTrace; row::ROW_COL_TYPE=1, col::ROW_COL_TYPE=1, secondary_y::Bool=false)
    # NOTE: This causes the domain to not be set for the first plot (row=1, col=1)
    # if row == 1 && col == 1
    #     push!(p.data, trace)
    #     return p
    # end

    _check_row_col_arg(p.layout, row, "row", 1)
    _check_row_col_arg(p.layout, col, "col", 2)

    gridref_row_ix = row isa String ? Colon() : row
    gridref_col_ix = col isa String ? Colon() : col

    _add_trace!(p, trace, gridref_row_ix, gridref_col_ix, secondary_y)
    p
end

## internal helpers
_get_colorway(p::Plot) = _get_colorway(p.layout)
function _get_colorway(l::Layout)
    D3_colorway = [
        "#1F77B4",
        "#FF7F0E",
        "#2CA02C",
        "#D62728",
        "#9467BD",
        "#8C564B",
        "#E377C2",
        "#7F7F7F",
        "#BCBD22",
        "#17BECF",
    ]
    Cycler(get(l.template.layout, :colorway, D3_colorway))
end

_get_seq_from_template_data(p::Plot, args...; kwargs...) = _get_seq_from_template_data(p.layout, args...; kwargs...)
function _get_seq_from_template_data(
        l::Layout,
        default::Cycler,
        property_sub_path::Symbol,
        root_template_path::Symbol=:scatter
    )
    template_specs = get(l.template.data, root_template_path, Dict())
    if !isempty(template_specs)
        seq = []
        default_ix = 0
        has_any_specs = false
        for spec in template_specs
            want = get(spec, property_sub_path, Dict())
            if ismissing(want) || isempty(want)
                push!(seq, default[default_ix += 1])
            else
                has_any_specs = true
                push!(seq, want)
            end
        end

        return has_any_specs ? Cycler(seq) : default
    end
    return default
end

function _get_line_dash_seq(l::Union{Layout,Plot})
    default_line_dash = Cycler([
        "solid",
        "dot",
        "dash",
        "longdash",
        "dashdot",
        "longdashdot",
    ])
    return _get_seq_from_template_data(l, default_line_dash, :line_dash, :scatter)
end

function _get_marker_symbol_seq(l::Union{Layout,Plot})
    default = Cycler([
        "circle"
        "diamond"
        "cross"
        "triangle"
        "square"
        "x"
        "pentagon"
        "hexagon"
        "hexagon2"
        "octagon"
        "star"
        "hexagram"
        "hourglass"
        "bowtie"
        "asterisk"
        "hash"
        "y"
        "line"
    ])
    return _get_seq_from_template_data(l, default, :marker_symbol, :scatter)
end

function _get_default_seq(l::Layout, attribute::Symbol)
    _getter_funcs = Dict(
        :line_dash => _get_line_dash_seq,
        :symbol => _get_marker_symbol_seq,
        :color => _get_colorway
    )
    if attribute in keys(_getter_funcs)
        return _getter_funcs[attribute](l)
    end
    error("Don't know how to get defaults for $(attribute)")
end

function _update_all_layout_type!(l::Layout, layout_type::Symbol, with::PlotlyAttribute)
    l[layout_type] = with

    for (k, v) in pairs(l)
        if startswith(string(k), string(layout_type))
            l[k] = with
        end
    end
    l
end

const _layout_obj_updaters = [:update_xaxes! => :xaxis, :update_yaxes! => :yaxis, :update_geos! => :geo, :update_mapboxes! => :mapbox, :update_polars! => :polar, :update_scenes! => :scene, :update_ternaries! => :ternary]

for (k1, k2) in _layout_obj_updaters
    @eval $(k1)(l::Layout, with::PlotlyAttribute=attr(); kwargs...) = _update_all_layout_type!(l, $(Meta.quot(k2)), merge(with, attr(;kwargs...)))
    @eval $(k1)(p::Plot, with::PlotlyAttribute=attr(); kwargs...) = $(k1)(p.layout, merge(with, attr(;kwargs...)))
    @eval export $k1
end

const _layout_vector_updaters = [:update_annotations! => :(layout.annotations), :update_shapes! => :(layout.shapes), :update_images! => :(layout.images)]

for (k1, k2) in _layout_vector_updaters
    @eval function $(k1)(layout::Layout, with::PlotlyAttribute; kwargs...)
        final_with = merge(with, attr(;kwargs...))
        for val in $(k2)
            merge!(val, final_with)
        end
    end
    @eval $(k1)(p::Plot, with::PlotlyAttribute=attr(); kwargs...) = $(k1)(p.layout, merge(with, attr(;kwargs...)))
    @eval export $(k1)
end

function _get_subplot_grid_refs(l::Layout, row::ROW_COL_TYPE, col::ROW_COL_TYPE)
    gridref_row_ix = row isa String ? Colon() : row
    gridref_col_ix = col isa String ? Colon() : col

    out = l.subplots.grid_ref[gridref_row_ix, gridref_col_ix]
    (out[1] isa Array) ? out : [out]
end

function _add_many_layout_objects!(
        l::Layout, base_val::HasFields, direction::Char, row::ROW_COL_TYPE, col::ROW_COL_TYPE, layout_key::Symbol
    )
    _check_row_col_arg(l, row, "row", 1)
    _check_row_col_arg(l, col, "col", 2)

    known_directions = ['h', 'v', 'X', 'I']
    if !(direction in known_directions)
        error("direction must be one of $(known_directions)")
    end

    layout_vec = get(l, layout_key, [])
    for refs in _get_subplot_grid_refs(l, row, col)
        # refs is always a vector...
        ref = refs[1]

        # find xref and yref
        xax, yax = string.(ref.layout_keys)
        xid = string("x", xax[6:end])
        yid = string("y", yax[6:end])

        new_val = deepcopy(base_val)
        if direction == 'h'
            new_val.xref = "$xid domain"
            new_val.yref = yid
        end
        if direction == 'v'
            new_val.xref = xid
            new_val.yref = "$yid domain"
        end
        if direction == 'X'
            new_val.xref = xid
            new_val.yref = yid
        end
        if direction == 'I'
            new_val.xref = "$xid domain"
            new_val.yref = "$yid domain"
        end
        push!(layout_vec, new_val)
    end
    l[layout_key] = layout_vec
end

_add_many_layout_objects!(p::Plot, args...) = _add_many_layout_objects!(p.layout, args...)
_add_many_shapes!(l::Layout, args...) = _add_many_layout_objects!(l, args..., :shapes)
_add_many_shapes!(p::Plot, args...) = _add_many_shapes!(p.layout, args...)

function add_hrect!(l::Union{Plot,Layout}, y0, y1; row::ROW_COL_TYPE="all", col::ROW_COL_TYPE="all", kw...)
    base_shape = rect(0, 1, y0, y1; kw...)
    _add_many_shapes!(l, base_shape, 'h', row, col)
end

function add_hline!(l::Union{Plot,Layout}, y; row::ROW_COL_TYPE="all", col::ROW_COL_TYPE="all", kw...)
    base_shape = hline(y; kw...)
    _add_many_shapes!(l, base_shape, 'h', row, col)
end

function add_vrect!(l::Union{Plot,Layout}, x0, x1; row::ROW_COL_TYPE="all", col::ROW_COL_TYPE="all", kw...)
    base_shape = rect(x0, x1, 0, 1; kw...)
    _add_many_shapes!(l, base_shape, 'v', row, col)
end

function add_vline!(l::Union{Plot,Layout}, x; row::ROW_COL_TYPE="all", col::ROW_COL_TYPE="all", kw...)
    base_shape = vline(x; kw...)
    _add_many_shapes!(l, base_shape, 'v', row, col)
end

function add_shape!(l::Union{Plot,Layout}, base_shape; row::ROW_COL_TYPE="all", col::ROW_COL_TYPE="all")
    _add_many_shapes!(l, base_shape, 'X', row, col)
end

function add_layout_image!(l::Union{Plot,Layout}, base_img; row::ROW_COL_TYPE="all", col::ROW_COL_TYPE="all")
    _add_many_layout_objects!(l, base_img, 'I', row, col, :images)
end

export add_hrect!, add_hline!, add_vrect!, add_vline!, add_shape!, add_layout_image!


text:
using Gnuplot, LinearAlgebra
using ParameterizedFunctions, DifferentialEquations
g = @ode_def begin
  dx = σ*(y - x)
  dy = x*(ρ - z) - y
  dz = x*y - β*z
end σ ρ β
u0 = [1.0; 0.0; 0.0]
tspan = (0.0, 80.0)
p = [10.0,28.0,8/3]
prob = ODEProblem(g, u0, tspan, p)
sol = solve(prob, Tsit5(), saveat = 0.01)
x, y, z = sol[1,:], sol[2,:], sol[3,:]
tempo = sol.t
# plot
@gsp "set xyplane at -3" "set auto fix" "set grid" :-
@gsp :- x y z tempo "w l notit lc palette" #palette(:plasma)
@gsp :- "set title 'Lorenz attractor'" "set cblabel 'time'"
@gsp :- xlab = "x" ylab = "y" zlab = "z"
save(term="pngcairo size 1200,800", output="lorenz.png") # hide

max_stars_repo_name: lazarusA/gnuplot-examples

max_issues_repo_name: lazarusA/gnuplot-examples

max_forks_repo_name: lazarusA/gnuplot-examples


text:
function make_seq(panel_file::AbstractString, profile_file::AbstractString, output_folder::AbstractString, depth)
    profile = JSON.parsefile(profile_file)
    config = profile["config"]
    fill_default_config!(config)

    # load assembly
    assembly = load_assembly(config["assembly"])
    panel, panel_size = load_bed(panel_file, assembly)

    if !isdir(output_folder)
        mkpath(output_folder)
    end

    if depth<=0
        depth = config["depth"]
    end

    io = nothing
    filename = joinpath(output_folder, "SEQMAKER-$SEQMAKER_VERSION")
    filename *= "-" * basename(panel_file)
    filename *= "-" * basename(profile_file)

    if config["pair-end"]
        io = fastq_open_pair("$filename.R1.fq", "$filename.R2.fq", "w")
    else
        io = fastq_open("$filename.R1.fq", "w")
    end

    temp_max = config["template_len"]["max"]
    temp_min = config["template_len"]["min"]
    temp_len_mean = ( temp_min + temp_max )/2
    read_num = (depth * panel_size) / (temp_len_mean * config["duplication_rate"])
    read_num = Int(floor(read_num))
    rand_pos = rand(read_num)
    rand_temp_len = rand(read_num)
    for i in 1:length(rand_pos)
        panel_pos = Int(floor(rand_pos[i] * panel_size))
        # linear interpolation
        temp_len = Int(floor(rand_temp_len[i] * temp_min + (1.0-rand_temp_len[i]) * temp_max))
        chr, pos_in_chr = seek_in_panel(panel, panel_pos)
        cycles = 1
        if haskey(profile, "cnv")
            copy_num = simulate_copy_number(chr, pos_in_chr, profile["cnv"])
            cycles = Int(round((rand() * copy_num)/0.5))
        end
        for c in 1:cycles
            # vibration
            this_pos_in_chr = pos_in_chr + Int(round((rand() - 0.5) * 0.4 * temp_len))
            this_temp_len = temp_len + Int(round((rand() - 0.5) * 0.1 * temp_len))
            # sampling
            original_seq, start = sample(assembly, chr, this_pos_in_chr , this_temp_len)
            seq = simulate_snv(original_seq, chr, start, profile["snv"])
            seq = simulate_fusion(seq, chr, start, profile["fusion"], assembly)
            # simulate watson/crick strand
            if rand()> 0.5
                seq = ~seq
            end
            if config["pair-end"]
                reads = pair_end_seq(seq, config)
                for pair in reads
                    fastq_write_pair(io, pair)
                end
            end
        end
    end
    close(io)
end

function seek_in_panel(panel, panel_pos)
    # binary search to find the contig
    left = 1
    right = length(panel)
    contig = 0
    while left < right
        center = div(left+right, 2)
        rec = panel[center]
        if rec["start_in_panel"]<=panel_pos && rec["end_in_panel"]>=panel_pos
            contig = center
            break
        elseif rec["start_in_panel"]>panel_pos
            right = center - 1
        elseif rec["end_in_panel"]<panel_pos
            left = center + 1
        end
    end
    if contig == 0
        contig = left
    end
    rec = panel[contig]
    offset = panel_pos - rec["start_in_panel"]
    pos_in_chr = rec["from"] + offset
    return rec["chr"], pos_in_chr
end

function fill_default_config!(config)
    if !haskey(config, "duplication_rate")
        config["duplication_rate"] = 1.0
    end
    if !haskey(config, "random_index1")
        config["random_index1"] = false
    end
    if !haskey(config, "random_index2")
        config["random_index2"] = false
    end
end


text:
# Eliminate Combinatorics dependency. These are based on Julia4 source.

import Base: eltype, length, iterate

struct Combinas{T}
    a::T
    t::Int
end

eltype(::Type{Combinas{T}}) where {T} = Vector{eltype(T)}

length(c::Combinas) = binomial(length(c.a),c.t)

function combinas(a, t::Integer)
    if t < 0
        # generate 0 combinations for negative argument
        t = length(a)+1
    end
    Combinas(a, t)
end

function iterate(c::Combinas, s=[1:c.t;])
    if !isempty(s) && s[1] > length(c.a)-c.t+1; return nothing; end
    comb = [c.a[si] for si in s]
    if c.t == 0
        # special case to generate 1 result for t==0
        return (comb,[length(c.a)+2])
    end
    s = copy(s)
    for i = length(s):-1:1
        s[i] += 1
        if s[i] > (length(c.a) - (length(s)-i))
            continue
        end
        for j = i+1:lastindex(s)
            s[j] = s[j-1]+1
        end
        break
    end
    (comb,s)
end

struct Permutas{T}
    a::T
end

eltype(::Type{Permutas{T}}) where {T} = Vector{eltype(T)}

length(p::Permutas) = factorial(length(p.a))

permutas(a) = Permutas(a)

function iterate(p::Permutas, s=[1:length(p.a);])
    if !isempty(s) && s[1] > length(p.a); return nothing; end
    perm = [p.a[si] for si in s]
    if length(p.a) == 0
        # special case to generate 1 result for len==0
        return (perm,[1])
    end
    s = copy(s)
    k = length(s)-1
    while k > 0 && s[k] > s[k+1];  k -= 1;  end
    if k == 0
        s[1] = length(s)+1   # done
    else
        l = length(s)
        while s[k] >= s[l];  l -= 1;  end
        s[k],s[l] = s[l],s[k]
        reverse!(s,k+1)
    end
    (perm,s)
end


text:
# formatting expression

### Argument specification

struct ArgSpec
    argidx::Int
    hasfilter::Bool
    filter::Function

    function ArgSpec(idx::Int, hasfil::Bool, filter::Function)
        idx != 0 || error("Argument index cannot be zero.")
        new(idx, hasfil, filter)
    end
end

getarg(args, sp::ArgSpec) =
    (a = args[sp.argidx]; sp.hasfilter ? sp.filter(a) : a)

# pos > 0: must not have iarg in expression (use pos+1), return (entry, pos + 1)
# pos < 0: must have iarg in expression, return (entry, -1)
# pos = 0: no positional argument before, can be either, return (entry, 1) or (entry, -1)
function make_argspec(s::AbstractString, pos::Int)
    # for argument position
    iarg::Int = -1
    hasfil::Bool = false
    ff::Function = Base.identity

    if !isempty(s)
        filrange = Compat.findfirst("|>", s)
        if filrange === nothing
            iarg = parse(Int, s)
        else
            ifil = first(filrange)
            iarg = ifil > 1 ? parse(Int, s[1:prevind(s, ifil)]) : -1
            hasfil = true
            ff = m_eval(Symbol(s[ifil+2:end]))
        end
    end

    if pos > 0
        iarg < 0 || error("entry with and without argument index must not coexist.")
        iarg = (pos += 1)
    elseif pos < 0
        iarg > 0 || error("entry with and without argument index must not coexist.")
    else # pos == 0
        if iarg < 0
            iarg = pos = 1
        else
            pos = -1
        end
    end

    return (ArgSpec(iarg, hasfil, ff), pos)
end

### Format entry

struct FormatEntry
    argspec::ArgSpec
    spec::FormatSpec
end

function make_formatentry(s::AbstractString, pos::Int)
    @assert s[1] == '{' && s[end] == '}'
    sc = s[2:prevind(s, lastindex(s))]
    icolon = Compat.findfirst(isequal(':'), sc)
    if icolon === nothing  # no colon
        (argspec, pos) = make_argspec(sc, pos)
        spec = FormatSpec('s')
    else
        (argspec, pos) = make_argspec(sc[1:prevind(sc, icolon)], pos)
        spec = FormatSpec(sc[nextind(sc, icolon):end])
    end
    return (FormatEntry(argspec, spec), pos)
end

### Format expression

mutable struct FormatExpr
    prefix::UTF8Str
    suffix::UTF8Str
    entries::Vector{FormatEntry}
    inter::Vector{UTF8Str}
end

_raise_unmatched_lbrace() = error("Unmatched { in format expression.")

function find_next_entry_open(s::AbstractString, si::Int)
    slen = lastindex(s)
    p = Compat.findnext(isequal('{'), s, si)
    (p === nothing || p < slen) || _raise_unmatched_lbrace()
    while p !== nothing && s[p+1] == '{'  # escape `{{`
        p = Compat.findnext(isequal('{'), s, p+2)
        (p === nothing || p < slen) || _raise_unmatched_lbrace()
    end
    pre = p !== nothing ? s[si:prevind(s, p)] : s[si:end]
    if !isempty(pre)
        pre = replace(pre, "{{" => '{')
        pre = replace(pre, "}}" => '}')
    end
    return (p, convert(UTF8Str, pre))
end

function find_next_entry_close(s::AbstractString, si::Int)
    p = Compat.findnext(isequal('}'), s, si)
    p !== nothing || _raise_unmatched_lbrace()
    return p
end

function FormatExpr(s::AbstractString)
    slen = length(s)

    # init
    prefix = UTF8Str("")
    suffix = UTF8Str("")
    entries = FormatEntry[]
    inter = UTF8Str[]

    # scan
    (p, prefix) = find_next_entry_open(s, 1)
    if p !== nothing
        q = find_next_entry_close(s, p+1)
        (e, pos) = make_formatentry(s[p:q], 0)
        push!(entries, e)
        (p, pre) = find_next_entry_open(s, q+1)
        while p !== nothing
            push!(inter, pre)
            q = find_next_entry_close(s, p+1)
            (e, pos) = make_formatentry(s[p:q], pos)
            push!(entries, e)
            (p, pre) = find_next_entry_open(s, q+1)
        end
        suffix = pre
    end
    FormatExpr(prefix, suffix, entries, inter)
end

function printfmt(io::IO, fe::FormatExpr, args...)
    isempty(fe.prefix) || print(io, fe.prefix)
    ents = fe.entries
    ne = length(ents)
    if ne > 0
        e = ents[1]
        printfmt(io, e.spec, getarg(args, e.argspec))
        for i = 2:ne
            print(io, fe.inter[i-1])
            e = ents[i]
            printfmt(io, e.spec, getarg(args, e.argspec))
        end
    end
    isempty(fe.suffix) || print(io, fe.suffix)
end

const StringOrFE = Union{AbstractString, FormatExpr}
printfmt(io::IO, fe::AbstractString, args...) = printfmt(io, FormatExpr(fe), args...)
printfmt(fe::StringOrFE, args...) = printfmt(_stdout(), fe, args...)

printfmtln(io::IO, fe::StringOrFE, args...) = (printfmt(io, fe, args...); println(io))
printfmtln(fe::StringOrFE, args...) = printfmtln(_stdout(), fe, args...)

format(fe::StringOrFE, args...) = sprint(printfmt, fe, args...)


text:
#########################################################
## Randomized singular value and eigen- decompositions ##
########################################################

# This file provides a rudimentary implementation of the randomized singular
# value decomposition and spectral (eigen-) decomposition as described in
# [^Halko2011].

import LinearAlgebra: Eigen, SVD

export rsvd, reigen

"""
    rsvd(A, n, p=0)

Compute partial singular value decomposition of `A` using a randomized
algorithm.

# Arguments

`A`: input matrix.

`n::Int`: number of singular value/vector pairs to find.

`p::Int=0`: number of extra vectors to include in computation.

# Output

`::SVD`: singular value decomposition.

!!! warning "Accuracy"
    This variant of the randomized singular value decomposition is the most
    commonly found implementation but is not recommended for accurate
    computations, as it often has trouble finding the `n` largest singular pairs,
    but rather finds `n` large singular pairs which may not necessarily be the
    largest.

# Implementation note

This function calls `rrange`, which uses naive randomized rangefinding to
compute a basis for a subspace of dimension `n` (Algorithm 4.1 of
[^Halko2011]), followed by `svd_restricted()`, which computes the
exact SVD factorization on the restriction of `A` to this randomly selected
subspace (Algorithm 5.1 of [^Halko2011]).

Alternatively, you can mix and match your own randomized algorithm using
any of the randomized range finding algorithms to find a suitable subspace
and feeding the result to one of the routines that computes the `SVD`
restricted to that subspace.
"""
function rsvd(A, n::Int, p::Int=0)
    Q = rrange(A, n+p)
    svd_restricted(A, Q, n)
end

"""
    rsvdvals(A, n, p=0)

Compute partial singular value decomposition of `A` using a randomized
algorithm.

# Arguments

- `A`: input matrix;
- `n::Int`: number of singular value/vector pairs to find;
- `p::Int=0`: number of extra vectors to include in computation.

# Output

- `::Vector`: estimated singular values of `A`.

!!! warning "Accuracy"
    This variant of the randomized singular value decomposition is the most
    commonly found implementation but is not recommended for accurate
    computations, as it often has trouble finding the `n` largest singular pairs,
    but rather finds `n` large singular pairs which may not necessarily be the
    largest.

# Implementation note

This function calls `rrange`, which uses naive randomized rangefinding to
compute a basis for a subspace of dimension `n` (Algorithm 4.1 of
[^Halko2011]), followed by `svd_restricted()`, which computes the
exact SVD factorization on the restriction of `A` to this randomly selected
subspace (Algorithm 5.1 of [^Halko2011]).

Alternatively, you can mix and match your own randomized algorithm using
any of the randomized range finding algorithms to find a suitable subspace
and feeding the result to one of the routines that computes the `SVD`
restricted to that subspace.
"""
function rsvdvals(A, n::Int, p::Int=0)
    Q = rrange(A, n+p)
    svdvals_restricted(A, Q, n)
end

"""
    rrange(A, r=0)

Compute an orthonormal basis for a subspace of `A` of dimension `r` using
naive randomized rangefinding.

# Arguments

- `A`: Input matrix. Must support `size(A)` and premultiply.
- `r::Int = 0`: The number of basis vectors to compute.

# Output

- `::Matrix`: matrix of dimension `size(A,1) x r` containing the basis
vectors of the computed subspace of `A`.

!!! warning "Warning"
    The Reference explicitly discourages using this algorithm.

# Implementation note

Whereas [^Halko2011] recommends classical Gram-Schmidt with double
reorthogonalization, we instead compute the basis with `qr()`, which
for dense `A` computes the QR factorization using Householder reflectors.
"""
function rrange(A, l::Int=0)
    m, n = size(A)
    if l > m
	    throw(ArgumentError("Cannot find $l linearly independent vectors of $m x $n matrix"))
    end
    Ω = randn(n, l)
    Y = A*Ω
    Q = Matrix(qr!(Y).Q)
    @assert m==size(Q, 1)
    @assert l==size(Q, 2)
    return Q
end

"""
    rrange_adaptive(A, r, ϵ=eps(); maxiter=10)

Compute an orthogonal basis for a subspace of `A` of dimension `l` using adaptive
randomized rangefinding.

Similar to `rrange`, but determines the oversampling parameter adaptively given
a threshold `ϵ`.

# Arguments

- `A`: input matrix. Must support `size(A)` and premultiply.
- `r::Integer`: number of basis vectors to compute.
- `ϵ::Real = eps()`: threshold to determine adaptive fitting.

## Keywords

- `maxiter::Int = 10`: maximum number of iterations to run.

# Output

`::Matrix`: matrix of dimension `size(A,1)` x l containing the basis
vectors of the computed subspace of `A`.

# References

Algorithm 4.2 of [^Halko2011]
"""
function rrange_adaptive(A, r::Integer, ϵ::Real=eps(); maxiter::Int=10)
    m, n = size(A)

    Ω = randn(n,r)
    Y = A*Ω
    Q = zeros(m,0)

    tol = ϵ / (10*√(2/π))
    for j=1:maxiter
        Tol = maximum([norm(Y[:,i]) for i=j:(j+r-1)])
        Tol > tol || break

        y = view(Y,:,j)
        y = Y[:,j] = y - Q*(Q'y)
        q = y/norm(y)
        Q = [Q q]

        ω = randn(n)
        y = A*ω
        y = y - Q*(Q'y)
        Y = [Y y]
        Yb = view(Y, :, (j+1):(j+r-1))
        Yb = Y[:, (j+1):(j+r-1)] = Yb - q * q'Yb

        j==maxiter && warn("Maximum number of iterations reached with norm $Tol > $tol")
    end
    Q
end

basis(V) = Matrix(qr(V).Q)

"""
    rrange_si(A, l; At=A', q=0)

Compute an orthonormal basis for a subspace of `A` of dimension `l` using
randomized rangefinding by subspace iteration.

*Note:* Running with `q=0` is functionally equivalent to `rrange`.

# Arguments

- `A`: input matrix. Must support `size(A)` and premultiply.
- `l::Int`: number of basis vectors to compute.

## Keywords

- `At = A'`: transpose of `A`.
- `q::Int = 0`: number of subspace iterations.
- `p` : oversampling parameter. The number of extra basis vectors to use
  in the computation, which get discarded at the end.

# Output

- `::Matrix`: A dense matrix of dimension `size(A,1) x l` containing the basis
  vectors of the computed subspace of `A`.

# Implementation note

Whereas the Reference recommends classical Gram-Schmidt with double
reorthogonalization, we instead compute the basis with `qr()`, which
for dense A computes the QR factorization using Householder reflectors.

# References

Algorithm 4.4 of [^Halko2011]
"""
function rrange_si(A, l::Int; At=A', q::Int=0)
    n = size(A, 2)
    Ω = randn(n,l+p)
    Y = A*Ω
    Q = basis(Y)
    for iter=1:q #Do some power iterations
        Ỹ=At*Q
        Q̃=basis(Ỹ)
        Y=A*Q̃
        Q=basis(Y)
    end
    Q
end

"""
    rrange_f(A, l)

Compute an orthonormal basis for a subspace of `A` of dimension `l` using
naive randomized rangefinding using stochastic randomized Fourier transforms.

*Note:* similar to `rrange`, but does not use gaussian random matrices.

# Arguments

- `A` : input matrix. Must support `size(A)` and premultiply.
- `l::Int` : number of basis vectors to compute.
- `p::Int = 0` : oversampling parameter. The number of extra basis vectors to use
  in the computation, which get discarded at the end.

# Output

- `::Matrix`: matrix of dimension `size(A,1)` x `l` containing the basis
  vectors of the computed subspace of `A`.

# Implementation note

Whereas the Reference recommends classical Gram-Schmidt with double
reorthogonalization, we instead compute the basis with `qr()`, which
for dense `A` computes the QR factorization using Householder reflectors.

# References

Algorithm 4.5 of [^Halko2011]
"""
function rrange_f(A, l::Int)
    n = size(A, 2)
    Ω = srft(l+p)
    Y = A*Ω
    Q = Matrix(qr!(Y).Q)
end

"""
    svd_restricted(A, Q, n)

Compute the SVD factorization of `A` restricted to the subspace spanned by `Q`
using exact projection.

# Arguments

- `A`: input matrix. Must support postmultiply.
- `Q`: matrix containing basis vectors of the subspace whose restriction to is
  desired.

# Output

`::SVD`: singular value decomposition.

# References

Algorithm 5.1 of [^Halko2011]
"""
function svd_restricted(A, Q, n::Int)
    B=Q'A
    S=svd!(B)
    SVD((Q*S.U)[:, 1:n], S.S[1:n], S.Vt[1:n, :])
end

"""
    svdvals_restricted(A, Q, n)

Compute the singular values of `A` restricted to the subspace spanned by `Q`
using exact projection.

# Arguments

- `A`: input matrix. Must support postmultiply.
- `Q`: matrix containing basis vectors of the subspace whose restriction to is
  desired.

# Output

- `::Vector`: estimated singular values of `A`.

# References

Algorithm 5.1 of [^Halko201]
"""
function svdvals_restricted(A, Q, n::Int)
    B=Q'A
    S=svdvals!(B)[1:n]
end

"""
    svd_re(A, Q)

Compute the SVD factorization of `A` restricted to the subspace spanned by `Q`
using row extraction.

*Note:* Remark 5.2 of [^Halko2011] recommends input of `Q` of the form `Q=A*Ω`
where `Ω` is a sample computed by `randn(n,l)` or even `srft(l)`.

# Arguments

- `A`: input matrix. Must support postmultiply
- `Q`: matrix containing basis vectors of the subspace whose restriction to is
  desired. Need not be orthogonal or normalized.

# Output

- `::SVD`: singular value decomposition.

# See also

A faster but less accurate variant of [`svd_restricted`](@ref) which uses the
interpolative decomposition `id`.

# References

Algorithm 5.2 of [^Halko2011]
"""
function svd_re(A, Q)
    F = id(Q)
    X, J = F[:B], F[:P]
    R′, W′= qr(A[J, :])
    Z = X*R′
    S=svd(Z)
    SVD(S.U, S.S, S.Vt * W′)
end

"""
    eigen_restricted(A, Q)

Compute the spectral (`Eigen`) factorization of `A` restricted to the subspace
spanned by `Q` using row extraction.

# Arguments

- `A::Hermitian`: input matrix. Must be `Hermitian` and support pre- and post-multiply.
- `Q`: orthonormal matrix containing basis vectors of the subspace whose
  restriction to is desired.

# Output

- `::LinearAlgebra.Eigen`: eigen factorization.

# References

Algorithm 5.3 of [^Halko2011]
"""
function eigen_restricted(A::Hermitian, Q)
    B = Q'A*Q
    E = eigen!(B)
    Eigen(E.values, Q*E.vectors)
end

"""
    eigen_re(A, Q)

Compute the spectral (`Eigen`) factorization of `A` restricted to the subspace
spanned by `Q` using row extraction.

*Note:* Remark 5.2 of [^Halko2011] recommends input of `Q` of the form `Q=A*Ω`
where `Ω` is a sample computed by `randn(n,l)` or even `srft(l)`.

# Arguments

- `A::Hermitian`: input matrix. Must be `Hermitian` and support pre- and post-multiply.
- `Q`: matrix containing basis vectors of the subspace whose restriction to is
  desired. Need not be orthogonal or normalized.

# Output

- `::LinearAlgebra.Eigen`: eigen factorization.

# See also

A faster but less accurate variant of `eigen_restricted()` which uses the
interpolative decomposition `id()`.

# References

Algorithm 5.4 of [^Halko2011]
"""
function eigen_re(A::Hermitian, Q)
    X, J = id(Q)
    F = qr!(X)
    V, R = F.Q, F.R
    Z=R*A[J, J]*R'
    E=eigen(Z)
    Eigen(E.values, V*E.vectors)
end

"""
    eigen_nystrom(A, Q)

Compute the spectral (`Eigen`) factorization of `A` restricted to the subspace
spanned by `Q` using the Nyström method.

# Arguments

- `A`: input matrix. Must be positive semidefinite.
- `Q`: orthonormal matrix containing basis vectors of the subspace whose
  restriction to is desired.

# Output

- `::LinearAlgebra.Eigen`: eigen factorization.

# See also

More accurate than [`eigen_restricted`](@ref) but is restricted to matrices
that can be Cholesky decomposed.

# References

Algorithm 5.5 of [^Halko2011]
"""
function eigen_nystrom(A, Q)
    B₁=A*Q
    B₂=Q'*B₁
    C=cholesky!(Hermitian(B₂))
    F=B₁/C
    S=svd!(F)
    Eigen(S.S.^2, S.U)
end

"""
    eigfact_onepass(A, Ω)

Compute the spectral (`Eigen`) factorization of `A` using only one matrix
product involving `A`.

# Arguments

- `A::Hermitian`: input matrix.
- `Ω`: sample matrix for the column space, e.g. `randn(n, l)` or `srft(l)`.
- `Ω̃;`: sample matrix for the row space. Not neeeded for `Hermitian` matrices.
- `At = A'`: computes transpose of input matrix.

# Output

- `::Base.LingAlg.Eigen`: eigen factorization.

# References

Algorithm 5.6 of [^Halko2011]
"""
function eigen_onepass(A::Hermitian, Ω)
    Y=A*Ω; Q = Matrix(qr!(Y).Q)
    B=(Q'Y)\(Q'Ω)
    E=eigen!(B)
    Eigen(E.values, Q*E.vectors)
end
function eigen_onepass(A, Ω, Ω̃; At=A')
    Y=A *Ω; Q = Matrix(qr!(Y).Q)
    Ỹ=At*Ω; Q̃ = Matrix(qr!(Ỹ).Q)
    #Want least-squares solution to (5.14 and 5.15)
    B=(Q'Y)\(Q̃'Ω)
    B̃=(Q̃'Ỹ)\(Q'Ω̃)
    #Here is a very very very hacky way to solve the problem
    B=0.5(B + B̃')
    E=eigen!(B)
    Eigen(E.values, Q*E.vectors)
end

"""
    reigen(A, l)

Compute the spectral (`Eigen`) decomposition of `A` using a randomized
algorithm.

# Arguments

- `A`: input matrix.
- `l::Int`: number of eigenpairs to find.

# Output

- `::LinearAlgebra.Eigen`: eigen decomposition.

# Implementation note

This is a wrapper around `eigen_onepass()` which uses the randomized
samples found using `srft(l)`.
"""
reigen(A::Hermitian, l::Int) = eigen_onepass(A, srft(l))
reigen(A, l::Int) = eigen_onepass(A, srft(l), srft(l))


text:
"""
    DotTest(m_rand,d_rand,operators,parameters)
Dot product test for a vector of linear operators
See also: [`ConjugateGradients`](@ref)

"""
function DotTest(m_rand,d_rand,operators,parameters)

	m_adj = LinearOperator(d_rand,operators,parameters,adj=true)
	d_fwd = LinearOperator(m_rand,operators,parameters,adj=false)
	inner1 = InnerProduct(d_rand,d_fwd)
	inner2 = InnerProduct(m_rand,m_adj)
	println("<d_rand,d_fwd> = ",inner1)
	println("<m_rand,m_adj> = ",inner2)
	return inner1,inner2

end


text:
begin
    inK = 2.0
    γ = 5 / 3
    wL = [1.13, 0.13, 0.8]
    wR = [0.3, -0.2, 1.5]
    primL = KitBase.conserve_prim(wL, γ)
    primR = KitBase.conserve_prim(wR, γ)

    u = collect(-5.0:0.2:5.0)
    ω = ones(length(u)) / length(u)

    hL = KitBase.maxwellian(u, primL)
    bL = hL .* inK ./ (2.0 * primL[end])
    hR = KitBase.maxwellian(u, primR)
    bR = hR .* inK ./ (2.0 * primR[end])

    dt = 1e-3
    dx = 1e-2

    fw = similar(wL)
    fh = similar(hL)
    fb = similar(bL)
end

#--- fluid flux ---#
KitBase.flux_upwind(rand(), rand(), rand(2), rand(2), dt)
KitBase.flux_lax!(fw, wL, wR, γ, dt, dx)
KitBase.flux_hll!(fw, wL, wR, γ, dt)
KitBase.flux_roe!(fw, wL, wR, γ, dt)
KitBase.flux_roe!(zeros(4), [1.0, 0.3, 0.0, 1.0], [0.3, -0.1, 0.0, 2.0], γ, dt)

#--- gks flux ---#
KitBase.flux_gks(0.3, 1e-3, dt)
KitBase.flux_gks(0.3, 1e-3, dt, 1e-1, 1.0)
KitBase.flux_gks(1.0, 0.125, 1e-3, dt, 1e-2, 1e-2)

KitBase.flux_gks!(zeros(3), [1.0, 0.0, 2.0], inK, γ, 1e-3, 0.72)
KitBase.flux_gks!(zeros(4), [1.0, 0.0, 0.0, 2.0], inK, γ, 1e-3, 0.72)
KitBase.flux_gks!(zeros(5), [1.0, 0.0, 0.0, 0.0, 2.0], inK, γ, 1e-3, 0.72)

KitBase.flux_gks!(fw, wL, wR, inK, γ, 1e-3, 0.72, dt, dx, dx)
KitBase.flux_gks!(
    zeros(4),
    [1.0, 0.3, 0.0, 1.0],
    [0.3, -0.1, 0.0, 2.0],
    inK,
    γ,
    1e-3,
    0.72,
    dt,
    dx,
    dx,
    dx,
) # 2D
KitBase.flux_gks!(
    zeros(4, 2),
    hcat([1.0, 0.3, 0.0, 1.0], [1.0, 0.3, 0.0, 1.0]),
    hcat([0.3, -0.1, 0.0, 2.0], [0.3, -0.1, 0.0, 2.0]),
    inK,
    γ,
    1.0,
    1.0,
    0.5,
    1.0,
    1e-2,
    dt,
    dx,
    dx,
    dx,
    zeros(4, 2),
    zeros(4, 2),
) # mixture
KitBase.flux_gks!(
    zeros(4),
    [1.0, 0.3, 0.0, 1.0],
    [0.3, -0.1, 0.0, 2.0],
    inK,
    γ,
    1e-3,
    0.72,
    dt,
) # FR

# discrete
KitBase.flux_gks!(fw, fh, wL, wR, u, inK, γ, 1e-3, 0.72, dt, dx, dx, zeros(3), zeros(3))
KitBase.flux_gks!(fw, fh, fb, wL, wR, u, inK, γ, 1e-3, 0.72, dt, dx, dx, zeros(3), zeros(3))
KitBase.flux_gks!(
    zeros(4),
    zeros(16, 16),
    [1.0, 0.3, 0.0, 1.0],
    [0.3, -0.1, 0.0, 2.0],
    rand(16, 16),
    rand(16, 16),
    γ,
    inK,
    1e-3,
    0.72,
    dt,
    dx,
    dx,
    dx,
    zeros(4),
    zeros(4),
) # 2D
KitBase.flux_gks!(
    zeros(4),
    zeros(16, 16),
    zeros(16, 16),
    [1.0, 0.3, 0.0, 1.0],
    [0.3, -0.1, 0.0, 2.0],
    rand(16, 16),
    rand(16, 16),
    γ,
    inK,
    1e-3,
    0.72,
    dt,
    dx,
    dx,
    dx,
    zeros(4),
    zeros(4),
) # 2D

KitBase.flux_ugks!(
    fw,
    fh,
    fb,
    wL,
    hL,
    bL,
    wR,
    hR,
    bR,
    u,
    ω,
    inK,
    γ,
    1e-3,
    0.81,
    0.72,
    dt,
    dx,
    dx,
)
KitBase.flux_ugks!(
    zeros(4),
    zeros(28, 28),
    zeros(28, 28),
    [1.0, 0.3, 0.0, 1.0],
    rand(28, 28),
    rand(28, 28),
    [0.3, -0.1, 0.0, 2.0],
    rand(28, 28),
    rand(28, 28),
    rand(28, 28),
    rand(28, 28),
    ones(28, 28),
    inK,
    γ,
    1e-3,
    0.81,
    0.72,
    dt,
    dx,
    dx,
    dx,
)
KitBase.flux_ugks!(
    zeros(5, 2),
    zeros(16, 16, 2),
    zeros(16, 16, 2),
    zeros(16, 16, 2),
    hcat([1.0, 0.3, 0.0, 0.0, 1.0], [1.0, 0.3, 0.0, 0.0, 1.0]),
    rand(16, 16, 2),
    rand(16, 16, 2),
    rand(16, 16, 2),
    hcat([0.3, -0.1, 0.0, 0.0, 2.0], [0.3, -0.1, 0.0, 0.0, 2.0]),
    rand(16, 16, 2),
    rand(16, 16, 2),
    rand(16, 16, 2),
    rand(16, 16, 2),
    rand(16, 16, 2),
    ones(16, 16, 2),
    inK,
    γ,
    1.0,
    1.0,
    0.5,
    1.0,
    0.1,
    dt,
    dx,
    dx,
    dx,
)

#--- kfvs flux ---#
KitBase.flux_kfvs!(fh, hL, hR, u, dt)
KitBase.flux_kfvs!(fw, fh, hL, hR, u, ω, dt)
KitBase.flux_kfvs!(
    zeros(3, 2),
    zeros(16, 2),
    rand(16, 2),
    rand(16, 2),
    randn(16, 2),
    ones(16, 2),
    dt,
)
KitBase.flux_kfvs!(fw, fh, fb, hL, bL, hR, bR, u, ω, dt)
KitBase.flux_kfvs!(
    zeros(3, 2),
    zeros(16, 2),
    zeros(16, 2),
    rand(16, 2),
    rand(16, 2),
    rand(16, 2),
    rand(16, 2),
    randn(16, 2),
    ones(16, 2),
    dt,
)
KitBase.flux_kfvs!(zeros(4), fh, fb, zero(fh), hL, bL, bL, hR, bR, bR, u, ω, dt) # Rykov
KitBase.flux_kfvs!(
    zeros(5),
    zeros(16, 16, 16),
    rand(16, 16, 16),
    rand(16, 16, 16),
    randn(16, 16, 16),
    randn(16, 16, 16),
    randn(16, 16, 16),
    ones(16, 16, 16),
    dt,
)
KitBase.flux_kfvs!(
    zeros(5),
    zeros(16, 16, 16),
    rand(16, 16, 16),
    rand(16, 16, 16),
    randn(16, 16, 16),
    randn(16, 16, 16),
    randn(16, 16, 16),
    ones(16, 16, 16),
    dt,
    1.0,
)
KitBase.flux_kfvs!(
    zeros(5),
    zeros(16),
    zeros(16),
    zeros(16),
    zeros(16),
    rand(16),
    rand(16),
    rand(16),
    rand(16),
    rand(16),
    rand(16),
    rand(16),
    rand(16),
    randn(16),
    ones(16),
    dt,
)
KitBase.flux_kfvs!(
    zeros(5, 2),
    zeros(16, 2),
    zeros(16, 2),
    zeros(16, 2),
    zeros(16, 2),
    rand(16, 2),
    rand(16, 2),
    rand(16, 2),
    rand(16, 2),
    rand(16, 2),
    rand(16, 2),
    rand(16, 2),
    rand(16, 2),
    randn(16, 2),
    ones(16, 2),
    dt,
)
KitBase.flux_kfvs!(
    zeros(4),
    zeros(16, 16),
    rand(16, 16),
    rand(16, 16),
    randn(16, 16),
    randn(16, 16),
    ones(16, 16),
    dt,
    dx,
)
KitBase.flux_kfvs!(
    zeros(4),
    zeros(16, 16),
    zeros(16, 16),
    rand(16, 16),
    rand(16, 16),
    rand(16, 16),
    rand(16, 16),
    randn(16, 16),
    randn(16, 16),
    ones(16, 16),
    dt,
    dx,
)
# 3F2V
KitBase.flux_kfvs!(
    zeros(5),
    zeros(16, 16),
    zeros(16, 16),
    zeros(16, 16),
    rand(16, 16),
    rand(16, 16),
    rand(16, 16),
    rand(16, 16),
    rand(16, 16),
    rand(16, 16),
    randn(16, 16),
    randn(16, 16),
    ones(16, 16),
    dt,
    dx,
)
KitBase.flux_kfvs!(
    zeros(5, 2),
    zeros(16, 16, 2),
    zeros(16, 16, 2),
    zeros(16, 16, 2),
    rand(16, 16, 2),
    rand(16, 16, 2),
    rand(16, 16, 2),
    rand(16, 16, 2),
    rand(16, 16, 2),
    rand(16, 16, 2),
    randn(16, 16, 2),
    randn(16, 16, 2),
    ones(16, 16, 2),
    dt,
    dx,
)

#--- kinetic central-upwind flux ---#
# 1F1V
KitBase.flux_kcu!(fw, fh, wL, hL, wR, hR, u, ω, inK, γ, 1e-3, 0.81, 0.72, dt)
KitBase.flux_kcu!(
    hcat(fw, fw),
    hcat(fh, fh),
    hcat(wL, wL),
    hcat(hL, hL),
    hcat(wR, wR),
    hcat(hR, hR),
    hcat(u, u),
    hcat(ω, ω),
    inK,
    γ,
    1.0,
    1.0,
    0.5,
    1.0,
    0.01,
    dt,
)
# 2F1V
KitBase.flux_kcu!(fw, fh, fb, wL, hL, bL, wR, hR, bR, u, ω, inK, γ, 1e-3, 0.81, 0.72, dt)
KitBase.flux_kcu!(
    zeros(3, 2),
    zeros(16, 2),
    zeros(16, 2),
    zeros(3, 2),
    rand(16, 2),
    rand(16, 2),
    zeros(3, 2),
    rand(16, 2),
    rand(16, 2),
    rand(16, 2),
    ones(16, 2),
    inK,
    γ,
    1.0,
    1.0,
    0.5,
    1.0,
    0.01,
    dt,
)
# 4F1V
KitBase.flux_kcu!(
    zeros(5),
    zeros(16),
    zeros(16),
    zeros(16),
    zeros(16),
    zeros(5),
    rand(16),
    rand(16),
    rand(16),
    rand(16),
    zeros(5),
    rand(16),
    rand(16),
    rand(16),
    rand(16),
    rand(16),
    ones(16),
    inK,
    γ,
    1e-3,
    0.81,
    0.72,
    dt,
)
KitBase.flux_kcu!(
    zeros(5, 2),
    zeros(16, 2),
    zeros(16, 2),
    zeros(16, 2),
    zeros(16, 2),
    zeros(5, 2),
    rand(16, 2),
    rand(16, 2),
    rand(16, 2),
    rand(16, 2),
    zeros(5, 2),
    rand(16, 2),
    rand(16, 2),
    rand(16, 2),
    rand(16, 2),
    rand(16, 2),
    ones(16, 2),
    inK,
    γ,
    1.0,
    1.0,
    0.5,
    1.0,
    0.01,
    dt,
)
# 1F2V
KitBase.flux_kcu!(
    zeros(4),
    zeros(16, 16),
    zeros(4),
    rand(16, 16),
    zeros(4),
    rand(16, 16),
    rand(16, 16),
    rand(16, 16),
    ones(16, 16),
    inK,
    γ,
    1e-3,
    0.81,
    0.72,
    dt,
    dx,
)
# 2F2V
KitBase.flux_kcu!(
    zeros(4),
    zeros(16, 16),
    zeros(16, 16),
    zeros(4),
    rand(16, 16),
    rand(16, 16),
    zeros(4),
    rand(16, 16),
    rand(16, 16),
    rand(16, 16),
    rand(16, 16),
    ones(16, 16),
    inK,
    γ,
    1e-3,
    0.81,
    0.72,
    dt,
    dx,
)
# 3F2V
KitBase.flux_kcu!(
    zeros(5),
    zeros(16, 16),
    zeros(16, 16),
    zeros(16, 16),
    zeros(5),
    rand(16, 16),
    rand(16, 16),
    rand(16, 16),
    zeros(5),
    rand(16, 16),
    rand(16, 16),
    rand(16, 16),
    rand(16, 16),
    rand(16, 16),
    ones(16, 16),
    inK,
    γ,
    1e-3,
    0.81,
    0.72,
    dt,
    dx,
)
KitBase.flux_kcu!(
    zeros(5, 2),
    zeros(16, 16, 2),
    zeros(16, 16, 2),
    zeros(16, 16, 2),
    zeros(5, 2),
    rand(16, 16, 2),
    rand(16, 16, 2),
    rand(16, 16, 2),
    zeros(5, 2),
    rand(16, 16, 2),
    rand(16, 16, 2),
    rand(16, 16, 2),
    rand(16, 16, 2),
    rand(16, 16, 2),
    ones(16, 16, 2),
    inK,
    γ,
    1.0,
    1.0,
    0.5,
    1.0,
    0.01,
    dt,
    dx,
)

#--- boundary flux ---#
KitBase.flux_boundary_maxwell!(zeros(3), rand(3), [1.0, 0.0, 1.0], 2.0, 5 / 3, 1e-3, 1)
KitBase.flux_boundary_maxwell!(
    zeros(4),
    rand(4),
    [1.0, 0.0, 0.0, 1.0],
    1.0,
    5 / 3,
    1e-3,
    1e-2,
    1,
)
KitBase.flux_boundary_maxwell!(
    zeros(3),
    zeros(16),
    zeros(16),
    [1.0, 0.3, 0.89],
    rand(16),
    rand(16),
    randn(16),
    ones(16),
    inK,
    dt,
    1,
)
KitBase.flux_boundary_maxwell!(
    zeros(4),
    zeros(16, 16),
    zeros(16, 16),
    [1.0, 0.3, -0.2, 0.89],
    rand(16, 16),
    rand(16, 16),
    randn(16, 16),
    randn(16, 16),
    ones(16, 16),
    inK,
    dt,
    dx,
    1,
)
KitBase.flux_boundary_maxwell!(
    zeros(5),
    zeros(16, 16, 16),
    [1.0, 0.3, 0.1, 0.1, 0.89],
    rand(16, 16, 16),
    randn(16, 16, 16),
    randn(16, 16, 16),
    randn(16, 16, 16),
    ones(16, 16, 16),
    dt,
    1,
)

KitBase.flux_boundary_specular!(zeros(3), zeros(16), rand(16), randn(16), ones(16), dt)
KitBase.flux_boundary_specular!(
    zeros(3),
    zeros(16),
    zeros(16),
    rand(16),
    rand(16),
    randn(16),
    ones(16),
    dt,
)

#--- pure equilibrium flux ---#
KitBase.flux_equilibrium!(fw, wL, wR, inK, γ, 1e-3, 0.81, 0.72, dt, dx, dx)
KitBase.flux_equilibrium!(
    zeros(4),
    zeros(4),
    zeros(4),
    inK,
    γ,
    1e-3,
    0.81,
    0.72,
    dt,
    dx,
    dx,
    dx,
)

#--- electromagnetic flux ---#
sol = 100.0
χ = 1.0
ν = 1.0
# A^+
Ap = Array{Float64}(undef, 8, 8)
Ap[1, 1] = (sol * χ) / 2.0
Ap[7, 1] = χ / 2.0
Ap[2, 2] = sol / 2.0
Ap[6, 2] = 0.5
Ap[3, 3] = sol / 2.0
Ap[5, 3] = -1.0 / 2.0
Ap[4, 4] = (sol * ν) / 2.0
Ap[8, 4] = (sol^2 * ν) / 2.0
Ap[3, 5] = -sol^2 / 2.0
Ap[5, 5] = sol / 2.0
Ap[2, 6] = sol^2 / 2.0
Ap[6, 6] = sol / 2.0
Ap[1, 7] = (sol^2 * χ) / 2.0
Ap[7, 7] = (sol * χ) / 2.0
Ap[4, 8] = ν / 2.0
Ap[8, 8] = (sol * ν) / 2.0
# A^-
An = Array{Float64}(undef, 8, 8)
An[1, 1] = -(sol * χ) / 2.0
An[7, 1] = χ / 2.0
An[2, 2] = -sol / 2.0
An[6, 2] = 1.0 / 2.0
An[3, 3] = -sol / 2.0
An[5, 3] = -1.0 / 2.0
An[4, 4] = -(sol * ν) / 2.0
An[8, 4] = (sol^2 * ν) / 2.0
An[3, 5] = -sol^2 / 2.0
An[5, 5] = -sol / 2.0
An[2, 6] = sol^2 / 2.0
An[6, 6] = -sol / 2.0
An[1, 7] = (sol^2 * χ) / 2.0
An[7, 7] = -(sol * χ) / 2.0
An[4, 8] = ν / 2.0
An[8, 8] = -(sol * ν) / 2.0
# eigenvalues
D = [sol, sol, sol * χ, sol * ν, -sol, -sol, -sol * χ, -sol * ν]

KitBase.flux_em!(
    zeros(8),
    zeros(8),
    randn(3),
    randn(3),
    randn(3),
    randn(3),
    randn(3),
    randn(3),
    randn(3),
    randn(3),
    0.0,
    0.0,
    0.0,
    0.0,
    dx,
    dx,
    Ap,
    An,
    D,
    100.0,
    1.0,
    1.0,
    dt,
)
KitBase.flux_emx!(
    zeros(8),
    zeros(8),
    zeros(8),
    zeros(8),
    zeros(8),
    zeros(8),
    randn(3),
    randn(3),
    randn(3),
    randn(3),
    randn(3),
    randn(3),
    randn(3),
    randn(3),
    0.0,
    0.0,
    0.0,
    0.0,
    dx,
    dx,
    Ap,
    An,
    Ap,
    An,
    D,
    100.0,
    1.0,
    1.0,
    dt,
)
KitBase.flux_emy!(
    zeros(8),
    zeros(8),
    zeros(8),
    zeros(8),
    zeros(8),
    zeros(8),
    randn(3),
    randn(3),
    randn(3),
    randn(3),
    randn(3),
    randn(3),
    randn(3),
    randn(3),
    0.0,
    0.0,
    0.0,
    0.0,
    dx,
    dx,
    Ap,
    An,
    Ap,
    An,
    D,
    100.0,
    1.0,
    1.0,
    dt,
)


text:
# MAT_v4.jl
# Tools for reading MATLAB v4 files in Julia
#
# Copyright (C) 2012   Simon Kornblith
# Copyright (C) 2019   Victor Saase (modified from MAT_v5.jl)
#
# Permission is hereby granted, free of charge, to any person obtaining
# a copy of this software and associated documentation files (the
# "Software"), to deal in the Software without restriction, including
# without limitation the rights to use, copy, modify, merge, publish,
# distribute, sublicense, and/or sell copies of the Software, and to
# permit persons to whom the Software is furnished to do so, subject to
# the following conditions:
#
# The above copyright notice and this permission notice shall be
# included in all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
# EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
# MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
# NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
# LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
# OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
# WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

# MATLAB's file format documentation can be found at
# http://www.mathworks.com/help/pdf_doc/matlab/matfile_format.pdf

module MAT_v4
using BufferedStreams, HDF5, SparseArrays
import Base: read, write, close
import HDF5: names, exists

round_uint8(data) = round.(UInt8, data)
complex_array(a, b) = complex.(a, b)

mutable struct Matlabv4File <: HDF5.DataFile
    ios::IOStream
    swap_bytes::Bool
    varnames::Dict{String, Int64}

    Matlabv4File(ios, swap_bytes) = new(ios, swap_bytes)
end

const mLITTLE_ENDIAN = 0
const mBIG_ENDIAN = 1
const mVAX_DFLOAT = 2
const mVAX_GFLOAT = 3
const mGRAY = 4

const pTYPE = Dict(
    0 => Float64,
    1 => Float32,
    2 => Int32,
    3 => Int16,
    4 => UInt16,
    5 => UInt8
)

const tNUMERIC = 0
const tTEXT = 1
const tSPARSE = 2

const imagfREAL = 0
const imagfCOMPLEX = 1

read_bswap(f::IO, swap_bytes::Bool, ::Type{T}) where T =
    swap_bytes ? bswap(read(f, T)) : read(f, T)
function read_bswap(f::IO, swap_bytes::Bool, ::Type{T}, dim::Union{Int, Tuple{Vararg{Int}}}) where T
    d = read!(f, Array{T}(undef, dim))
    if swap_bytes
        for i = 1:length(d)
            @inbounds d[i] = bswap(d[i])
        end
    end
    d
end

function read_bswap(f::IO, swap_bytes::Bool, d::AbstractArray{T}) where T
    readbytes!(f, reinterpret(UInt8, d))
    if swap_bytes
        for i = 1:length(d)
            @inbounds d[i] = bswap(d[i])
        end
    end
    d
end

function checkv4(f::IO)
    M, O, P, T, mrows, ncols, imagf, namlen = MAT_v4.read_header(f, false)
    if 0<=M<=4 && O == 0 && 0<=P<=5 && 0<=T<=2 && mrows>=0 && ncols>=0 && 0<=imagf<=1 && namlen>0
        swap_bytes = false
        return (true, swap_bytes)
    else
        seek(f, 0)
        M, O, P, T, mrows, ncols, imagf, namlen = MAT_v4.read_header(f, true)
        if 0<=M<=4 && O == 0 && 0<=P<=5 && 0<=T<=2 && mrows>=0 && ncols>=0 && 0<=imagf<=1 && namlen>0
            swap_bytes = true
            return (true, swap_bytes)
        end
    end
    return (false, false)
end

# Read data type and number of bytes at the start of a data element
function read_header(f::IO, swap_bytes::Bool)
    dtype = read_bswap(f, swap_bytes, Int32)

    M = div(rem(dtype, 10000), 1000)
    O = div(rem(dtype, 1000), 100)
    P = div(rem(dtype, 100), 10)
    T = div(rem(dtype, 10), 1)

    mrows = read_bswap(f, swap_bytes, Int32)
    ncols = read_bswap(f, swap_bytes, Int32)
    imagf = read_bswap(f, swap_bytes, Int32)
    namlen = read_bswap(f, swap_bytes, Int32)

    M, O, P, T, mrows, ncols, imagf, namlen
end

# Read matrix data
function read_matrix(f::IO, swap_bytes::Bool)
    M, O, P, T, mrows, ncols, imagf, namlen = read_header(f, swap_bytes)
    if ncols == 0 || mrows == 0
        # If one creates a cell array using
        #     y = cell(m, n)
        # then MATLAB will save the empty cells as zero-byte matrices. If one creates a
        # empty cells using
        #     a = {[], [], []}
        # then MATLAB does not save the empty cells as zero-byte matrices. To avoid
        # surprises, we produce an empty array in both cases.
        return ("", Matrix{Union{}}(undef, 0, 0))
    end
    name = String(read_bswap(f, M==mBIG_ENDIAN, Vector{UInt8}(undef, namlen))[1:end-1])
    if T == tNUMERIC || T == tSPARSE
        real_data = read_bswap(f, M==mBIG_ENDIAN, Vector{pTYPE[P]}(undef, ncols*mrows))
        if T == tNUMERIC && imagf == imagfCOMPLEX
            imag_data = read_bswap(f, M==mBIG_ENDIAN, Vector{pTYPE[P]}(undef, ncols*mrows))
            data = complex.(real_data, imag_data)
        else
            data = real_data
        end
        datamat = reshape(data, Int(mrows), Int(ncols))
        if T == tNUMERIC
            return (name, datamat)
        elseif T == tSPARSE
            if size(datamat,2) == 3
                return (name, sparse(datamat[:,1], datamat[:,2], datamat[:,3]))
            else
                return (name, sparse(datamat[:,1], datamat[:,2], complex.(datamat[:,3], datamat[:,4])))
            end
        end
    elseif T == tTEXT
        if mrows > 1
            charvec = UInt8.(read_bswap(f, M==mBIG_ENDIAN, Vector{pTYPE[P]}(undef, ncols*mrows)))
            charmat = reshape(charvec, Int(mrows), Int(ncols))
            data = [String(charmat[i,:]) for i in 1:mrows]
        else
            data = String(UInt8.(read_bswap(f, M==mBIG_ENDIAN, Vector{pTYPE[P]}(undef, ncols))))
        end
        return (name, data)
    end
end

# Open MAT file for reading
matopen(ios::IOStream, endian_indicator::Bool) =
    Matlabv4File(ios, endian_indicator)

# Read whole MAT file
function read(matfile::Matlabv4File)
    seek(matfile.ios, 0)
    vars = Dict{String, Any}()
    while !eof(matfile.ios)
        (name, data) = read_matrix(matfile.ios, matfile.swap_bytes)
        vars[name] = data
    end
    vars
end

# Read only variable names
function getvarnames(matfile::Matlabv4File)
    if !isdefined(matfile, :varnames)
        seek(matfile.ios, 0)
        matfile.varnames = varnames = Dict{String, Int64}()
        while !eof(matfile.ios)
            offset = position(matfile.ios)
            M, O, P, T, mrows, ncols, imagf, namlen = read_header(matfile.ios, matfile.swap_bytes)
            f = matfile.ios

            name = String(read_bswap(f, M==mBIG_ENDIAN, Vector{UInt8}(undef, namlen))[1:end-1])
            varnames[name] = offset
            imag_offset = 0
            skip(f, mrows*ncols*sizeof(pTYPE[P]))
            if imagf == imagfCOMPLEX
                skip(f, mrows*ncols*sizeof(pTYPE[P]))
            end
        end
    end
    matfile.varnames
end

exists(matfile::Matlabv4File, varname::String) =
    haskey(getvarnames(matfile), varname)
names(matfile::Matlabv4File) =
    keys(getvarnames(matfile))

# Read a variable from a MAT file
function read(matfile::Matlabv4File, varname::String)
    varnames = getvarnames(matfile)
    if !haskey(varnames, varname)
        error("no variable $varname in file")
    end
    seek(matfile.ios, varnames[varname])
    (name, data) = read_matrix(matfile.ios, matfile.swap_bytes)
    data
end

function colvals(A::AbstractSparseMatrix)
    rows = rowvals(A)
    cols = similar(rows)
    m,n = size(A)
    for i=1:n
        for j in nzrange(A,i)
            cols[j] = i
        end
    end
    cols
end

function write(parent::Matlabv4File, name::String, s)
    M = Int(parent.swap_bytes)
    O = 0
    P = 0
    for p=keys(pTYPE)
        if eltype(s) == pTYPE[p] || eltype(s) == Complex{pTYPE[p]}
            P = p
        end
    end
    if pTYPE[P] != eltype(s) && Complex{pTYPE[P]} != eltype(s) &&
        !(s isa AbstractString && pTYPE[P] == Float64) &&
        !(s isa Vector{String} && pTYPE[P] == Float64)
        error("invalid value type when writing v4 file")
    end
    if s isa AbstractSparseMatrix
        T = tSPARSE
    elseif s isa AbstractString || s isa Vector{String}
        T = tTEXT
    else
        T = tNUMERIC
    end
    write(parent.ios, Int32(1000*M + 100*O + 10*P + T))

    mrows = 1
    ncols = 1
    if s isa AbstractVector && !(s isa Vector{String})
        ncols = length(s)
    elseif s isa AbstractMatrix
        if s isa AbstractSparseMatrix
            mrows = nnz(s)
            ncols = 3
            if eltype(s) <: Complex
                ncols = 4
            end
        else
            mrows, ncols = size(s)
        end
    elseif s isa Vector{String}
        ncols = length(s[1])
        mrows = length(s)
    elseif s isa AbstractString
        ncols = length(s)
    end
    write(parent.ios, Int32(mrows))
    write(parent.ios, Int32(ncols))

    imagf = 0
    if eltype(s) <: Complex && T == tNUMERIC
        imagf = 1
    end
    write(parent.ios, Int32(imagf))

    namlen = length(name) + 1
    write(parent.ios, Int32(namlen))

    write(parent.ios, Vector{UInt8}(name))
    write(parent.ios, UInt8(0))

    if s isa AbstractArray && T == tNUMERIC
        write(parent.ios, reshape(real(s), length(s)))
        if imagf == 1
            write(parent.ios, reshape(imag(s), length(s)))
        end
    elseif s isa Number
        write(parent.ios, real(s))
        if imagf == 1
            write(parent.ios, imag(s))
        end
    elseif s isa AbstractString
        floatarray = Float64.(Vector{UInt8}(s))
        write(parent.ios, floatarray)
    elseif s isa Vector{String}
        floatarray = Matrix{Float64}(undef, mrows, ncols)
        for (i,strel) = enumerate(s)
            floatarray[i,:] = Float64.(Vector{UInt8}(strel))
        end
        write(parent.ios, floatarray)
    elseif T == tSPARSE
        rows = rowvals(s)
        cols = colvals(s)
        vals = nonzeros(s)
        write(parent.ios, pTYPE[P].(rows))
        write(parent.ios, pTYPE[P].(cols))
        write(parent.ios, pTYPE[P].(real(vals)))
        if eltype(s) <: Complex
            write(parent.ios, pTYPE[P].(imag(vals)))
        end
    end
end

# Close MAT file
close(matfile::Matlabv4File) = close(matfile.ios)

end


text:
"""
    crop(x; to)
    crop(xs...; to)

Crop one or multiple [`AbstractRaster`](@ref) or [`AbstractRasterStack`](@ref) `x`
to match the size of the object `to`, or smallest of any dimensions that are shared.

Otherwise crop to the size of the keyword argument `to`. This can be a
`Tuple` of `Dimension` or any object that will return one from `dims(to)`.

`crop` is lazy, using a `view` into the object rather than alocating new memory.

# Keywords

- `to`: the object to crop to. If `to` keyword is passed, the smallest shared
    area of all `xs` is used.

As `crop` is lazy, `filename` and `suffix` keywords don't apply.

# Example

Cropt to another raster:

```jldoctest
using Rasters, Plots
evenness = Raster(EarthEnv{HabitatHeterogeneity}, :evenness)
rnge = Raster(EarthEnv{HabitatHeterogeneity}, :range)

# Roughly cut out New Zealand from the evenness raster
nz_bounds = X(Between(165, 180)), Y(Between(-32, -50))
nz_evenness = evenness[nz_bounds...]

# Crop range to match evenness
nz_range = crop(rnge; to=nz_evenness)
plot(nz_range)

savefig("build/crop_example.png")
# output
```

![new zealand evennes cropped]/nz_crop_example.png)

Crop to a polygon:

```julia
using Rasters, Plots, Dates, Shapefile, Downloads
using Rasters.LookupArrays

# Download a borders shapefile
shapefile_url = "https://github.com/nvkelso/natural-earth-vector/raw/master/10m_cultural/ne_10m_admin_0_countries.shp"
shapefile_name = "boundary.shp"
isfile(shapefile_name) || Downloads.download(shapefile_url, shapefile_name)
shp = Shapefile.Handle(shapefile_name).shapes[6]
argentina_range = crop(evenness; to=shp)
plot(argentina_range)
savefig("build/argentina_crop_example.png")
```

![argentina evenness cropped]/argentina_crop_example.png)

$EXPERIMENTAL
"""
function crop end
function crop(l1::RasterStackOrArray, l2::RasterStackOrArray, ls::RasterStackOrArray...; kw...)
    crop((l1, l2, ls...); kw...)
end
function crop(xs; to=nothing, kw...)
    if isnothing(to)
        to = _subsetbounds((max, min), xs)
        map(l -> _crop_to_bounds(l, to), xs)
    else
        map(l -> crop(l; to, kw...), xs)
    end
end
crop(x::RasterStackOrArray; to, kw...) = _crop_to(x, to; kw...)

# crop `A` to values of dims of `to`
_crop_to(A, to::RasterStackOrArray; kw...) = _crop_to(A, dims(to); kw...)
function _crop_to(x, to::DimTuple; atol=maybe_eps(to))
    # We can only crop to sampled dims (e.g. not categorical dims like Band)
    sampled = reduce(to; init=()) do acc, d
        lookup(d) isa AbstractSampled ? (acc..., d) : acc
    end
    wrapped_bounds = map(d ->  rebuild(d, bounds(d)), sampled)
    return _crop_to_bounds(x, wrapped_bounds)
end
function _crop_to(x, to; order=nothing)
    order, wrapped_bounds = _order_and_bounds(x, to, order)
    return _crop_to_bounds(x, wrapped_bounds)
end

function _crop_to_bounds(ser::RasterSeries, wrapped_bounds::DimTuple)
    map(x -> _crop_to_bounds(x, wrapped_bounds), ser)
end
function _crop_to_bounds(x::RasterStackOrArray, wrapped_bounds::DimTuple)
    # Take a view over the bounds
    _without_mapped_crs(x) do x1
        dimranges = map(wrapped_bounds) do d
            x_d = dims(x1, d)
            range = DD.selectindices(x_d, LA.ClosedInterval(parent(d)...))
            rebuild(x_d, range)
        end
        view(x1, dimranges...)
    end
end

"""
    extend(xs...; [to])
    extend(xs; [to])
    extend(x::Union{AbstractRaster,AbstractRasterStack}; to, kw...)

Extend one or multiple [`AbstractRaster`](@ref) to match the area
covered by all `xs`, or by the keyword argument `to`.

# Keywords

- `to`: the Raster or dims to extend to. If no `to` keyword is passed, the largest
    shared area of all `xs` is used.
- `atol`: the absolute tolerance value to use when comparing the index of x and `to`.
- `filename`: a filename to write to directly, useful for large files.
- `suffix`: a string or value to append to the filename.
    A tuple of `suffix` will be applied to stack layers. `keys(st)` are the default.

```jldoctest
using Rasters, Plots
evenness = Raster(EarthEnv{HabitatHeterogeneity}, :evenness)
rnge = Raster(EarthEnv{HabitatHeterogeneity}, :range)

# Roughly cut out South America
sa_bounds = X(Between(-88, -32)), Y(Between(-57, 13))
sa_evenness = evenness[sa_bounds...]

# Extend range to match the whole-world raster
sa_range = extend(sa_evenness; to=rnge)
plot(sa_range)

savefig("build/extend_example.png")
# output
```

![extend](extend_example.png)

$EXPERIMENTAL
"""
function extend end
function extend(l1::RasterStackOrArray, l2::RasterStackOrArray, ls::RasterStackOrArray...; kw...)
    extend((l1, l2, ls...); kw...)
end
function extend(xs; to=nothing)
    if isnothing(to)
        to = _subsetbounds((min, max), xs)
        map(l -> _extend_to_bounds(l, to), xs)
    else
        map(l -> extend(l; to), xs)
    end
end
extend(x::RasterStackOrArray; to=dims(x), kw...) = _extend_to(x, to; kw...)

_extend_to(x::RasterStackOrArray, to::RasterStackOrArray; kw...) = _extend_to(x, dims(to); kw...)
function _extend_to(x, to; order=nothing, kw...)
    order, wrapped_bounds = _order_and_bounds(x, to, order)
    all(map(s -> s isa Regular, span(x, order))) || throw(ArgumentError("All dims must have `Regular` span to be extended with a polygon"))
    return _extend_to_bounds(x, wrapped_bounds; kw...)
end
function _extend_to(A::AbstractRaster, to::DimTuple;
    filename=nothing, suffix=nothing
)
    # Calculate the range of the old array in the extended array
    ranges = _without_mapped_crs(A) do A
        _without_mapped_crs(to) do to
            map(dims(A), to) do d, t
                range = DD.selectindices(t, LA.ClosedInterval(bounds(d)...))
                rebuild(d, range)
            end
        end
    end
    # Create a new extended array
    newA = create(filename, eltype(A), to;
        suffix, parent=parent(A), missingval=missingval(A),
        name=name(A), metadata=metadata(A)
    )
    # The missingval may have changed for disk-based arrays
    if !isequal(missingval(A), missingval(newA))
        A = replace_missing(A, missingval(newA))
    end
    open(newA; write=true) do O
        # Fill it with missing/nodata values
        O .= missingval(O)
        # Copy the original data to the new array
        # Somehow this is slow from disk?
        broadcast_dims!(identity, view(O, ranges...), A)
    end
    return newA
end
function _extend_to(st::AbstractRasterStack, to::Tuple; suffix=keys(st), kw...)
    mapargs((A, s) -> _extend_to(A, to; suffix=s, kw...), st, suffix)
end

function _extend_to_bounds(ser::RasterSeries, to::DimTuple; kw...)
    map(x -> _extend_to_bounds(x, wrapped_bounds; kw...), ser)
end
function _extend_to_bounds(x, wrapped_bounds; kw...)
    newdims = map(wrapped_bounds) do wb
        d = dims(x, wb)
        b = parent(wb)
        l = lookup(d)
        # Use ranges for math because they have TwicePrecision magic
        # Define a range down to the lowest value,
        # but anchored at the existing value
        fl = LA.ordered_first(l); ll = LA.ordered_last(l)
        fb = first(b); lb = last(b)
        s = step(l)
        lowerrange = fb < first(bounds(l)) ? (fl:-s:fb) : (fl:-s:fl)
        upperrange = lb > last(bounds(l))  ? (ll: s:lb) : (ll: s:ll)
        if DD.order(l) isa ForwardOrdered
            newrange = last(lowerrange):s:last(upperrange)
        elseif order(d) isa ReverseOrdered
            newrange = last(upperrange):s:last(lowerrange)
        end
        newlookup = rebuild(l; data=newrange)
        rebuild(d, newlookup)
    end
    return _extend_to(x, newdims; kw...)
end

# Shared utils

# Get the largest or smallest dimensions in a tuple of AbstractRaster
function _subsetbounds(fs, layers)
    # Combine the dimensions of all layers
    dims = DD.combinedims(layers...; check=false)
    # Search through all the dimensions choosing the shortest
    alldims = map(DD.dims, layers)
    reduce(dims; init=()) do acc, d
        all(map(l -> hasdim(l, d), layers)) || return acc
        matchingdims = map(ds -> DD.dims(ds, (d,)), alldims)
        bounds = reduce(matchingdims) do a, b
            _choosebounds(fs, a, b)
        end
        return (acc..., rebuild(d, bounds))
    end
end

_order_and_bounds(ser::AbstractRasterSeries, to, order) = _order_and_bounds(first(ser), to, order)
function _order_and_bounds(x, to, order)
    if Tables.istable(to)
        order = _table_point_order(dims(x), to, order)
        wrapped_bounds = _wrapped_table_bounds(dims(x), to, order)
    else
        order = isnothing(order) ? DEFAULT_POINT_ORDER : order
        wrapped_bounds = _wrapped_geom_bounds(dims(x), to, order)
    end
    return order, wrapped_bounds
end

# Choose bounds from either missing dimension
# (empty Tuple) or a comparison between two 1-Tuples
_choosebounds((f1, f2), (a,)::Tuple{<:Dimension}, ::Tuple{}) = bounds(a)
_choosebounds((f1, f2), (a,)::Tuple{<:Dimension}, b::Tuple{<:Dimension}) = _choosebounds((f1, f2), bounds(a), b)
_choosebounds((f1, f2), ::Tuple{}, ::Tuple{}) = ()
_choosebounds((f1, f2), ::Tuple{}, (b,)::DimTuple) = bounds(b)
_choosebounds((f1, f2), a::Tuple, ::Tuple{}) = a
function _choosebounds((f1, f2), (a1, a2)::Tuple{<:Any,<:Any}, (b,)::Tuple{<:Dimension})
    b1, b2 = bounds(b)
    return f1(a1, b1), f2(a2, b2)
end


text:
function test_ksp()

  @testset "\n   ---testing KSP solvers---" begin
  # create vectors and matrices

    b = make_vec()

    low, high = VecGetOwnershipRange(b)
    b_global_indices = collect(low:PetscInt(high - 1))

    x = make_vec()
    low, high = VecGetOwnershipRange(x)
    x_global_indices = collect(low:PetscInt(high - 1))

    println("making matrix for ksp solve")
    A = make_mat()
    MatTranspose(A, inplace=true)

    low, high = MatGetOwnershipRange(A)
    mat_global_indices = collect(low:PetscInt(high - 1))

    # perform solve
    ksp = KSP(comm)
    SetOperators(ksp, A, A)
    SetFromOptions(ksp)
    SetUp(ksp)
    KSPSolve(ksp, b, x)
    reason = GetConvergedReason(ksp)

    @test  reason  > 0# convergence

    # copy solution back to Julia
    x_copy = zeros(PetscScalar, sys_size_local)
    VecGetValues(x, sys_size_local, x_global_indices, x_copy)
    for i=1:sys_size_local
        @test isapprox( x_copy[i], x_julia[i]) atol=1e-14
    end

    PetscDestroy(ksp)

    println("   \n--- Testing LGMRES --- ")
    # test using non default KSP method

    # perform solve
    ksp = KSP(comm)
    SetOperators(ksp, A, A)

    # test PC
    println("   \n---Testing PC Options ---")
    pc = KSPGetPC(ksp)

    PCSetType(pc, PETSc2.PCBJACOBI)
    pctype = PCGetType(pc)
    @test ( pctype )== PETSc2.PCBJACOBI
    #PCFactorSetUseInPlace(pc, PetscBool(true))

    ### Do some KSP setup

    SetFromOptions(ksp)
    SetTolerances(ksp, rtol, abstol, dtol, maxits)
    SetInitialGuessNonzero(ksp, PetscBool(true))
    SetType(ksp, PETSc2.KSPLGMRES)
    SetUp(ksp)

    ### More PC testing
    n_local, first_local, ksp_arr = PCBJacobiGetSubKSP(pc)
    ksp2 = ksp_arr[1]

    pc2 = KSPGetPC(ksp2)

    PCSetReusePreconditioner(pc2, PetscBool(true))
    @test ( PCGetReusePreconditioner(pc2) )== true

    PCFactorSetAllowDiagonalFill(pc2, PetscBool(true))
    @test ( PCFactorGetAllowDiagonalFill(pc2) )== true

    PCFactorSetLevels(pc2, PetscInt(1))
    @test ( PCFactorGetLevels(pc2) )== 1  # should be pc2

    PCSetReusePreconditioner(pc2, PetscBool(true))
    @test ( PCGetReusePreconditioner(pc2) )== true

    PCFactorSetFill(pc, PetscReal(7.0))

    #=
    PCJacobiSetType(pc, PETSc2.PC_JACOBI_ROWMAX)
    @test ( PCJacobiGetType(pc) )== PETSc2.PC_JACOBI_ROWMAX
    =#

    #=
    tmp = PCFactorGetUseInPlace(pc)
    println("tmp = ", tmp)
    @test ( PCFactorGetUseInPlace(pc) )== true
    =#

    KSPSolve(ksp, b, x)
    reason = GetConvergedReason(ksp)
    ksptype = GetType(ksp)

    @test ( ksptype )== PETSc2.KSPLGMRES

    rtol_ret, abstol_ret, dtol_ret, maxits_ret = GetTolerances(ksp)

    @test isapprox( rtol_ret, rtol)
    @test isapprox( abstol_ret, abstol)
    @test isapprox( dtol_ret, dtol)
    @test ( maxits_ret )== maxits
    @test ( GetInitialGuessNonzero(ksp) )== true
    @test  reason  > 0# convergence

    rnorm = GetResidualNorm(ksp)
    @test  rnorm  < abstol

    # copy solution back to Julia
    x_copy = zeros(PetscScalar, sys_size_local)

    VecGetValues(x, sys_size_local, x_global_indices, x_copy)
    for i=1:sys_size_local
        @test isapprox( x_copy[i], x_julia[i]) atol=1e-14
    end

    PetscDestroy(x)
    PetscDestroy(b)
    PetscDestroy(A)
    PetscDestroy(ksp)
  end

return nothing
end

test_ksp()


text:
#=
characteristic scatterplots for the paper
=#
using DrWatson
quickactivate(@__DIR__, "VelocityMTD")
using MusicManipulations, PyPlot, StatsBase, Statistics, Parameters
include(srcdir("extract_data.jl"))
include(srcdir("style.jl"))

fig, axs = subplots(1,4,figsize = (2figx, figx/2))

datasets = ["tapping", "uwe", "pgmusic", "playalongs", "playalongs"]
datanames = ["tapping",  "Uwe Meile (piano)", "PG music (piano)", "drums playalongs", "drums playalongs (lead)"]
subsets = ["all", "melody", "melody", "all", "lead"]
datano = [1,1,1,8,8]

for i in 1:5
    # i == 3 && continue
    ax = axs[min(i, 4)]
    dataset = datasets[i]
    subset = subsets[i]
    c = @dict dataset subset
    mtds, vels = get_dataset(c)
    m, v = mtds[datano[i]], vels[datano[i]]
    ax.plot(m, v, ls = "None", color = "C$i", marker = "o", alpha = 0.5)
    if i ≤ 4
        ax.text(0.5, 1.05, datanames[i], va="center", ha="center",
        transform =ax.transAxes, color = "C$i", size = 32)
    else
        ax.text(0.2, 0.9, subsets[i], va="center", ha="center",
        transform =ax.transAxes, color = "C$i", size = 32)
    end
end

axs[1].set_ylabel("velocity", size = 32)
axs[4].set_ylim(0, 140)
# i ∈ (2, 4) && ax.set_xlabel("MTD")
add_identifiers!(fig)
fig.text(0.5, 0.05, "MTD", transform = fig.transFigure, size = 32)
fig.tight_layout()
fig.subplots_adjust(bottom=0.2,wspace = 0.22, hspace = 0.2, top = 0.9)
# fig.savefig(papersdir("figs", "scatterplots"))

max_stars_repo_name: Datseris/VelocityMTD

max_issues_repo_name: Datseris/VelocityMTD

max_forks_repo_name: Datseris/VelocityMTD


text:
"""
    nondimensional_exchange_functions(
        Δw,
        Δb,
        εδ_model_vars,
    )

Returns the nondimensional entrainment and detrainment
functions following Cohen et al. (JAMES, 2020), given:
 - `Δw`, updraft - environment vertical velocity differnce
 - `Δb`, updraft - environment buoynacy differnce
 - `εδ_model_vars`, entrainment detrainment model type
"""
function nondimensional_exchange_functions(param_set, Δw, εδ_model_vars)

    μ_0 = CPEDMF.μ_0(param_set)
    μ = ICP.entrainment_sigma(param_set)
    β = CPEDMF.β(param_set)
    χ = CPEDMF.χ(param_set)

    Δb = εδ_model_vars.b_up - εδ_model_vars.b_en
    μ_ij = (χ - εδ_model_vars.a_up / (εδ_model_vars.a_up + εδ_model_vars.a_en)) * Δb / Δw
    exp_arg = μ / μ_0 * μ_ij
    D_ε = 1.0 / (1.0 + exp(-exp_arg))
    D_δ = 1.0 / (1.0 + exp(exp_arg))

    M_δ = (max((εδ_model_vars.RH_up)^β - (εδ_model_vars.RH_en)^β, 0.0))^(1.0 / β)
    M_ε = (max((εδ_model_vars.RH_en)^β - (εδ_model_vars.RH_up)^β, 0.0))^(1.0 / β)

    return D_ε, D_δ, M_δ, M_ε
end;


text:
function limit(x :: Float64)
    if abs(x) < 10000.
        x
    elseif x < 0
        -10000.
    else
        10000.
    end
end

function bp_transformation!(transformation, inputs, d, new_d)
    new_d .= deriv(transformation).(inputs)
    for i in 1:length(d)
        new_d[i] *= d[i]
        new_d[i] = limit(new_d[i])
    end
end

function bp_matrix!(inputs, matrix, d, step, new_d)
    ## step = 0.01
    ## dmatrix = reshape(d, (length(d), 1)) * reshape(inputs, (1, length(inputs)))
    ## matrix[:, :] += step * reshape(d, (length(d), 1)) * reshape(inputs, (1, length(inputs)))
    for j in 1:size(matrix, 2)
        new_d[j] = 0.0
        for i in 1:size(matrix, 1)
            matrix[i, j] += step * limit(d[i] * inputs[j])
            matrix[i, j] = limit(matrix[i, j])
            new_d[j] += matrix[i, j] * d[i]
            new_d[j] = limit(new_d[j])
        end
    end
end

function bp_function!(op :: DFunction, value :: Dict)
    bp_transformation!(op.f_transformation, value[:pre_transform], value[:d], value[:pre_transform_d])
    bp_matrix!(value[:pre_matrix], op.f_matrix, value[:pre_transform_d], op.learning_rate, value[:pre_matrix_d])
end

function bp_tree!(tree :: Tree)
    if typeof(tree.op) == DFunction
        bp_function!(tree.op, tree.value)
        for i in 1:length(tree.subtrees)
            tree.subtrees[i].value[:d] .= tree.value[:input_ds][i]
            bp_tree!(tree.subtrees[i])
        end
        ## ds = lens_split(d, [c.feature_len for c in tree.op.in_classes])
        ## for i in 1:length(tree.subtrees)
        ##     bp_tree!(tree.subtrees[i], tree.value[:ds][i])
        ## end
        ## dds = map(bp_class!, tree.op.in_classes, [t.value for t in tree.subtrees], tree.value[:ds])
        ## foreach(bp_tree!, tree.subtrees, dds)
    end
end

max_stars_repo_name: Non-Contradiction/NextGeneration

max_issues_repo_name: Non-Contradiction/NextGeneration

max_forks_repo_name: Non-Contradiction/NextGeneration


text:
function pow(x, n)
  r = 1
  while n > 0
    n -= 1
    r *= x
  end
  return r
end

# Invoking the compiler outside of the genfuncs appears
# to make the specialised versions visible inside them,
# leading to a significant first-compile speedup.
@code_adjoint pow(2, 3)

@code_adjoint ((x, y) -> sum(x.*y))([1,2,3],[4,5,6])

gradient(pow, 2, 3)
gradient((x, y) -> sum(x.*y), [1, 2, 3], [4, 5, 6])

gradient(x -> sum(Float32.(x)), [1.0])