text:
#Is this usable data?
function isconcrete{T<:Number}(x::T)
    !isna(x) && isfinite(x)
end

function isconcrete(x::(@compat Irrational))
    return true
end

function isconcrete(x)
    !isna(x)
end

function hasna(xs)
    return false
end

function hasna(xs::AbstractDataArray)
    for x in xs
        if isna(x)
            return true
        end
    end
    return false
end

function isallconcrete(xs)
    ans = true
    for x in xs
        if !isconcrete(x)
            ans = false
            break
        end
    end
    return ans
end

function concretize(xss::AbstractVector...)
    if all(map(isallconcrete, xss))
        return xss
    end

    count = 0
    for j in 1:length(xss[1])
        for xs in xss
            if !isconcrete(xs[j])
                @goto next_j1
            end
        end

        count += 1
        @label next_j1
    end

    yss = Array(AbstractVector, length(xss))
    for (i, xs) in enumerate(xss)
        yss[i] = Array(eltype(xs), count)
    end

    k = 1
    for j in 1:length(xss[1])
        for xs in xss
            if !isconcrete(xs[j])
                @goto next_j2
            end
        end

        for (i, xs) in enumerate(xss)
            yss[i][k] = xs[j]
        end
        k += 1

        @label next_j2
    end

    return tuple(yss...)
end

# How many concrete elements in an iterable
function concrete_length(xs)
    n = 0
    for x in xs
        if isconcrete(x)
            n += 1
        end
    end
    n
end

function concrete_length{T}(xs::DataArray{T})
    n = 0
    for i = 1:length(xs)
        if !xs.na[i] && isconcrete(xs.data[i]::T)
            n += 1
        end
    end
    n
end

function concrete_length(xs::Iterators.Chain)
    n = 0
    for obj in xs.xss
        n += concrete_length(obj)
    end
    n
end

function concrete_minimum(xs)
    if done(xs, start(xs))
        error("argument must not be empty")
    end

    x_min = first(xs)
    for x in xs
        if Gadfly.isconcrete(x) && isfinite(x)
            x_min = x
            break
        end
    end

    for x in xs
        if Gadfly.isconcrete(x) && isfinite(x) && x < x_min
            x_min = x
        end
    end
    return x_min
end

function concrete_maximum(xs)
    if done(xs, start(xs))
        error("argument must not be empty")
    end

    x_max = first(xs)
    for x in xs
        if Gadfly.isconcrete(x) && isfinite(x)
            x_max = x
            break
        end
    end

    for x in xs
        if Gadfly.isconcrete(x) && isfinite(x) && x > x_max
            x_max = x
        end
    end
    return x_max
end

function concrete_minmax{T<:Real}(xs, xmin::T, xmax::T)
    if eltype(xs) <: Base.Callable
        return xmin, xmax
    end

    for x in xs
        if isconcrete(x)
            xT = convert(T, x)
            if isnan(xmin) || xT < xmin
                xmin = xT
            end
            if isnan(xmax) || xT > xmax
                xmax = xT
            end
        end
    end
    xmin, xmax
end

function concrete_minmax{T<:Real, TA}(xs::DataArray{TA}, xmin::T, xmax::T)
    for i = 1:length(xs)
        if !xs.na[i]
            x = xs.data[i]::TA
            xT = convert(T, x)
            if isnan(xmin) || xT < xmin
                xmin = xT
            end
            if isnan(xmax) || xT > xmax
                xmax = xT
            end
        end
    end
    xmin, xmax
end

function concrete_minmax{T}(xs, xmin::T, xmax::T)
    for x in xs
        if isconcrete(x)
            xT = convert(T, x)
            if xT < xmin
                xmin = xT
            end
            if xT > xmax
                xmax = xT
            end
        end
    end
    xmin, xmax
end

function concrete_minmax{T, TA}(xs::DataArray{TA}, xmin::T, xmax::T)
    for i = 1:length(xs)
        if !xs.na[i]
            x = xs.data[i]::TA
            xT = convert(T, x)
            if xT < xmin
                xmin = xT
            end
            if xT > xmax
                xmax = xT
            end
        end
    end
    xmin, xmax
end

function concrete_minmax{T<:Real}(xs::Iterators.Chain, xmin::T, xmax::T)
    for obj in xs.xss
        xmin, xmax = concrete_minmax(obj, xmin, xmax)
    end
    xmin, xmax
end

function nonzero_length(xs)
    n = 0
    for x in xs
        if x != 0
            n += 1
        end
    end
    n
end

# Create a new object of type T from a with missing values (i.e., those set to
# nothing) inherited from b.
function inherit{T}(a::T, b::T)
    c = copy(a)
    inherit!(c, b)
    c
end

function inherit!{T}(a::T, b::T)
    for field in fieldnames(T)
        aval = getfield(a, field)
        bval = getfield(b, field)
        # TODO: this is a hack to let non-default labelers overide the defaults
        if aval === nothing || aval === string || aval == showoff
            setfield!(a, field, bval)
        elseif typeof(aval) <: Dict && typeof(bval) <: Dict
            merge!(aval, getfield(b, field))
        end
    end
    nothing
end

isnothing(u) = is(u, nothing)
issomething(u) = !isnothing(u)

negate(f) = x -> !f(x)

function has{T,N}(xs::AbstractArray{T,N}, y::T)
    for x in xs
        if x == y
            return true
        end
    end
    return false
end

Maybe(T) = @compat(Union{T, (@compat Void)})

function lerp(x::Float64, a, b)
    a + (b - a) * max(min(x, 1.0), 0.0)
end

# Generate a unique id, primarily for assigning IDs to SVG elements.
let next_unique_svg_id_num = 0
    global unique_svg_id
    function unique_svg_id()
        uid = @sprintf("id%d", next_unique_svg_id_num)
        next_unique_svg_id_num += 1
        uid
    end
end

# Remove any markup or whitespace from a string.
function escape_id(s::AbstractString)
    s = replace(s, r"<[^>]*>", "")
    s = replace(s, r"\s", "_")
    s
end

# Can a numerical value be treated as an integer
is_int_compatable(::Integer) = true
is_int_compatable{T <: AbstractFloat}(x::T) = abs(x) < maxintfloat(T) && float(int(x)) == x
is_int_compatable(::Any) = false

## Return a DataFrame with x, y column suitable for plotting a function.
#
# Args:
#  f: Function/Expression to be evaluated.
#  a: Lower bound.
#  b: Upper bound.
#  n: Number of points to evaluate the function at.
#
# Returns:
#  A data frame with "x" and "f(x)" columns.
#
function evalfunc(f::Function, a, b, n)
    @assert n > 1

    step = (b - a) / (n - 1)
    xs = Array(typeof(a + step), n)
    for i in 1:n
        xs[i] = a + (i-1) * step
    end

    df = DataFrame(xs, map(f, xs))
    # NOTE: 'colnames!' is the older deprecated name. 'names!' was also defined
    # but threw an error.
    try
        names!(df, [:x, :f_x])
    catch
        colnames!(df, ["x", "f_x"])
    end
    df
end

# Construct a jscall to store arbitrary data in the element
function jsdata(key::AbstractString, value::AbstractString, arg::Vector{Measure}=Measure[])
    return jscall(
        """
        data("$(escape_string(key))", $(value))
        """, arg)
end

# Construct jscall properties to store arbitrary data in plotroot elements.
function jsplotdata(key::AbstractString, value::AbstractString, arg::Vector{Measure}=Measure[])
    return jscall(
        """
        plotroot().data("$(escape_string(key))", $(value))
        """, arg)
end

function svg_color_class_from_label(label::AbstractString)
    return @sprintf("color_%s", escape_id(label))
end

"""
A faster map function for PooledDataVector
"""
function pooled_map(T::Type, f::Function, xs::PooledDataVector)
    newpool = T[f(x) for x in xs.pool]
    return T[newpool[i] for i in xs.refs]
end

using Base.Dates

# Arbitrarily order colors
function color_isless(a::Color, b::Color)
    return color_isless(convert(RGB{Float32}, a), convert(RGB{Float32}, b))
end

function color_isless(a::TransparentColor, b::TransparentColor)
    return color_isless(convert(RGBA{Float32}, a), convert(RGBA{Float32}, b))
end

function color_isless(a::RGB{Float32}, b::RGB{Float32})
    if a.r < b.r
        return true
    elseif a.r == b.r
        if a.g < b.g
            return true
        elseif a.g == b.g
            return a.b < b.b
        else
            return false
        end
    else
        return false
    end
end

function color_isless(a::RGBA{Float32}, b::RGBA{Float32})
    if color_isless(color(a), color(b))
        return true
    elseif color(a) == color(b)
        return a.alpha < b.alpha
    else
        return false
    end
end

function group_color_isless{S, T <: Colorant}(a::(@compat Tuple{S, T}),
                                              b::(@compat Tuple{S, T}))
    if a[1] < b[1]
        return true
    elseif a[1] == b[1]
        return color_isless(a[2], b[2])
    else
        return false
    end
end

"""
trim longer arrays to the size of the smallest one
and zip the arrays.
"""
function trim_zip(xs...)
    mx = max(map(length, xs)...)
    mn = min(map(length, xs)...)
    if mx == mn
        zip(xs...)
    else
        zip([length(x) == mn ? x : x[1:mn] for x in xs]...)
    end
end

ext: jl


text:
"""
    abstract type StructOfConstraints <: MOI.ModelLike end

A struct storing a subfields other structs storing constraints of different
types.

See [`Utilities.@struct_of_constraints_by_function_types`](@ref)
and [`Utilities.@struct_of_constraints_by_set_types`](@ref).
"""
abstract type StructOfConstraints <: MOI.ModelLike end

function _add_variable(model::StructOfConstraints)
    model.num_variables += 1
    return broadcastcall(_add_variable, model)
end
function _add_variables(model::StructOfConstraints, n)
    model.num_variables += n
    return broadcastcall(Base.Fix2(_add_variables, n), model)
end

function final_touch(::Nothing, index_map) end
function final_touch(model::StructOfConstraints, index_map)
    broadcastcall(Base.Fix2(final_touch, index_map), model)
    return
end

function _throw_if_cannot_delete(model::StructOfConstraints, vis, fast_in_vis)
    broadcastcall(model) do constrs
        if constrs !== nothing
            _throw_if_cannot_delete(constrs, vis, fast_in_vis)
        end
        return
    end
    return
end

function _deleted_constraints(
    callback::Function,
    model::StructOfConstraints,
    vi,
)
    broadcastcall(model) do constrs
        if constrs !== nothing
            _deleted_constraints(callback, constrs, vi)
        end
        return
    end
    return
end

function MOI.add_constraint(
    model::StructOfConstraints,
    func::F,
    set::S,
) where {F<:MOI.AbstractFunction,S<:MOI.AbstractSet}
    if !MOI.supports_constraint(model, F, S)
        throw(MOI.UnsupportedConstraint{F,S}())
    end
    return MOI.add_constraint(constraints(model, F, S), func, set)
end

function constraints(
    model::StructOfConstraints,
    ci::MOI.ConstraintIndex{F,S},
) where {F,S}
    if !MOI.supports_constraint(model, F, S)
        throw(MOI.InvalidIndex(ci))
    end
    return constraints(model, F, S)
end

function MOI.get(
    model::StructOfConstraints,
    attr::Union{MOI.ConstraintFunction,MOI.ConstraintSet},
    ci::MOI.ConstraintIndex,
)
    return MOI.get(constraints(model, ci), attr, ci)
end

function MOI.delete(model::StructOfConstraints, ci::MOI.ConstraintIndex)
    MOI.delete(constraints(model, ci), ci)
    return
end

function MOI.is_valid(
    model::StructOfConstraints,
    ci::MOI.ConstraintIndex{F,S},
) where {F,S}
    if !MOI.supports_constraint(model, F, S)
        return false
    end
    return MOI.is_valid(constraints(model, ci), ci)
end

function MOI.modify(
    model::StructOfConstraints,
    ci::MOI.ConstraintIndex,
    change::MOI.AbstractFunctionModification,
)
    MOI.modify(constraints(model, ci), ci, change)
    return
end

function MOI.set(
    model::StructOfConstraints,
    attr::Union{MOI.ConstraintFunction,MOI.ConstraintSet},
    ci::MOI.ConstraintIndex,
    func_or_set,
)
    MOI.set(constraints(model, ci), attr, ci, func_or_set)
    return
end

function MOI.get(
    model::StructOfConstraints,
    attr::MOI.ListOfConstraintTypesPresent,
)
    return broadcastvcat(model) do constrs
        if constrs === nothing
            return Tuple{Type,Type}[]
        end
        return MOI.get(constrs, attr)
    end
end

function MOI.get(
    model::StructOfConstraints,
    attr::MOI.NumberOfConstraints{F,S},
)::Int64 where {F,S}
    if !MOI.supports_constraint(model, F, S)
        return 0
    end
    return MOI.get(constraints(model, F, S), attr)
end

function MOI.get(
    model::StructOfConstraints,
    attr::MOI.ListOfConstraintIndices{F,S},
) where {F,S}
    if !MOI.supports_constraint(model, F, S)
        return MOI.ConstraintIndex{F,S}[]
    end
    return MOI.get(constraints(model, F, S), attr)
end

function MOI.is_empty(model::StructOfConstraints)
    return mapreduce_constraints(&, model, true) do constrs
        return constrs === nothing || MOI.is_empty(constrs)
    end
end

function MOI.empty!(model::StructOfConstraints)
    model.num_variables = 0
    broadcastcall(model) do constrs
        if constrs !== nothing
            MOI.empty!(constrs)
        end
        return
    end
    return
end

# Can be used to access constraints of a model
"""
    broadcastcall(f::Function, model::StructOfConstraints)

Calls `f(contrs)` for every field in `model`.
"""
function broadcastcall end

"""
    broadcastvcat(f::Function, model::StructOfConstraints)

Calls `f(contrs)` for every field in `model` and `vcat`s the results.
"""
function broadcastvcat end

"""
    mapreduce_constraints(
        f::Function,
        op::Function,
        model::StructOfConstraints,
        init,
    )

Call `mapreduce` on every field of `model` given an initial value `init`. Each
element in the map is computed as `f(x)` and the elements are reduced using
`op`.
"""
function mapreduce_constraints end

# Macro code

abstract type SymbolFS end

struct SymbolFun <: SymbolFS
    s::Union{Symbol,Expr}
    typed::Bool
end

struct SymbolSet <: SymbolFS
    s::Union{Symbol,Expr}
    typed::Bool
end

_typed(s::SymbolFS) = s.typed ? Expr(:curly, esc(s.s), esc(:T)) : esc(s.s)

# Base.lowercase is moved to Unicode.lowercase in Julia v0.7
import Unicode

function _field(s::SymbolFS)
    return Symbol(replace(Unicode.lowercase(string(s.s)), "." => "_"))
end

# Represents a union of function or set types
struct _UnionSymbolFS{S<:SymbolFS}
    s::Vector{S}
end

function _typed(s::_UnionSymbolFS)
    tt = _typed.(s.s)
    return Expr(:curly, :Union, tt...)
end

_field(s::_UnionSymbolFS) = _field(s.s[1])

function _mapreduce_field(s::Union{SymbolFS,_UnionSymbolFS})
    return :(cur = op(cur, f(model.$(_field(s)))))
end

_callfield(f, s::Union{SymbolFS,_UnionSymbolFS}) = :($f(model.$(_field(s))))

function _parse_expr(::Type{S}, expr::Symbol) where {S<:SymbolFS}
    return S(expr, false)
end

function _parse_expr(::Type{S}, expr::Expr) where {S<:SymbolFS}
    if Meta.isexpr(expr, :curly)
        @assert length(expr.args) >= 1
        if expr.args[1] == :Union
            # `Union{:A, :B}` parses as
            # `Expr(:curly, :Union, :A, :B)`
            @assert length(expr.args) >= 3
            return _UnionSymbolFS{S}(_parse_expr.(S, expr.args[2:end]))
        else
            # Typed set, e.g. `MOI.EqualTo{T}` parses as:
            # `Expr(:curly, :(MOI.EqualTo), :T)`
            @assert length(expr.args) == 2
            @assert expr.args[2] == :T
            return S(expr.args[1], true)
        end
    else
        return S(expr, false)
    end
end

"""
    struct_of_constraint_code(struct_name, types, field_types = nothing)

Given a vector of `n` `Union{SymbolFun,_UnionSymbolFS{SymbolFun}}` or
`Union{SymbolSet,_UnionSymbolFS{SymbolSet}}` in `types`, defines a subtype of
`StructOfConstraints` of name `name` and which type parameters
`{T, F1, F2, ..., Fn}` if `field_types` is `nothing` and
a `{T}` otherwise.
It contains `n` field where the `i`th field has type `Ci` if `field_types` is
`nothing` and type `field_types[i]` otherwise.
If `types` is vector of `Union{SymbolFun,_UnionSymbolFS{SymbolFun}}` (resp.
`Union{SymbolSet,_UnionSymbolFS{SymbolSet}}`) then the constraints of that
function (resp. set) type are stored in the corresponding field.

This function is used by the macros [`@model`](@ref),
[`@struct_of_constraints_by_function_types`](@ref) and
[`@struct_of_constraints_by_set_types`](@ref).
"""
function struct_of_constraint_code(struct_name, types, field_types = nothing)
    T = esc(:T)
    typed_struct = :($(struct_name){$T})
    type_parametrized = field_types === nothing
    if type_parametrized
        field_types = [Symbol("C$i") for i in eachindex(types)]
        append!(typed_struct.args, field_types)
    end
    code = quote
        mutable struct $typed_struct <: StructOfConstraints
            num_variables::Int64
        end

        function $MOIU.broadcastcall(f::Function, model::$struct_name)
            $(Expr(:block, _callfield.(Ref(:f), types)...))
            return
        end

        function $MOIU.broadcastvcat(f::Function, model::$struct_name)
            return vcat($(_callfield.(Ref(:f), types)...))
        end

        function $MOIU.mapreduce_constraints(
            f::Function,
            op::Function,
            model::$struct_name,
            cur,
        )
            return $(Expr(:block, _mapreduce_field.(types)...))
        end
    end

    for (t, field_type) in zip(types, field_types)
        field = _field(t)
        push!(code.args[2].args[3].args, :($field::Union{Nothing,$field_type}))
        fun = t isa SymbolFun ? _typed(t) : :(MOI.AbstractFunction)
        set = t isa SymbolFun ? :(MOI.AbstractSet) : _typed(t)
        constraints_code = :(
            function $MOIU.constraints(
                model::$typed_struct,
                ::Type{<:$fun},
                ::Type{<:$set},
            )::$(field_type) where {$T}
                if model.$field === nothing
                    model.$field = $(field_type)()
                    $MOI.Utilities._add_variables(
                        model.$field,
                        model.num_variables,
                    )
                end
                return model.$field
            end
        )
        if type_parametrized
            append!(constraints_code.args[1].args, field_types)
        end
        push!(code.args, constraints_code)
    end
    if !isempty(types)
        is_func = any(types) do t
            return t isa Union{SymbolFun,_UnionSymbolFS{SymbolFun}}
        end
        is_set = any(types) do t
            return t isa Union{SymbolSet,_UnionSymbolFS{SymbolSet}}
        end
        @assert xor(is_func, is_set)
        if is_func
            SuperF = :(Union{$(_typed.(types)...)})
        else
            SuperF = :(MOI.AbstractFunction)
        end
        if is_set
            SuperS = :(Union{$(_typed.(types)...)})
        else
            SuperS = :(MOI.AbstractSet)
        end
        push!(
            code.args,
            :(
                function $MOI.supports_constraint(
                    model::$struct_name{$T},
                    ::Type{F},
                    ::Type{S},
                ) where {$T,F<:$SuperF,S<:$SuperS}
                    return $MOI.supports_constraint(
                        constraints(model, F, S),
                        F,
                        S,
                    )
                end
            ),
        )
    end
    # Here we make an internal constructor for our type. It needs to be an
    # internal constructor to avoid someone trying to provide
    # `typed_struct(args...)` which will throw an error about `T not defined`.
    constructor_code = :(function $typed_struct() where {$T}
        return new{$T}(0, $([:(nothing) for _ in field_types]...))
    end)
    if type_parametrized
        # A bit confusing: we need to add the extra type parameters to `{T}`.
        # Here's the `where {T}`
        append!(constructor_code.args[1].args, field_types)
        # And here's the `new{T}`
        append!(
            constructor_code.args[2].args[end].args[1].args[1].args,
            field_types,
        )
    end
    # Here's we we put the inner constructor.
    push!(code.args[2].args[3].args, constructor_code)
    return code
end

"""
    Utilities.@struct_of_constraints_by_function_types(name, func_types...)

Given a vector of `n` function types `(F1, F2,..., Fn)` in `func_types`, defines
a subtype of `StructOfConstraints` of name `name` and which type parameters
`{T, C1, C2, ..., Cn}`.
It contains `n` field where the `i`th field has type `Ci` and stores the
constraints of function type `Fi`.

The expression `Fi` can also be a union in which case any constraint for which
the function type is in the union is stored in the field with type `Ci`.
"""
macro struct_of_constraints_by_function_types(name, func_types...)
    funcs = _parse_expr.(SymbolFun, func_types)
    return struct_of_constraint_code(esc(name), funcs)
end

"""
    Utilities.@struct_of_constraints_by_set_types(name, func_types...)

Given a vector of `n` set types `(S1, S2,..., Sn)` in `func_types`, defines
a subtype of `StructOfConstraints` of name `name` and which type parameters
`{T, C1, C2, ..., Cn}`.
It contains `n` field where the `i`th field has type `Ci` and stores the
constraints of set type `Si`.
The expression `Si` can also be a union in which case any constraint for which
the set type is in the union is stored in the field with type `Ci`.
This can be useful if `Ci` is a [`MatrixOfConstraints`](@ref) in order to
concatenate the coefficients of constraints of several different set types in
the same matrix.
"""
macro struct_of_constraints_by_set_types(name, set_types...)
    sets = _parse_expr.(SymbolSet, set_types)
    return struct_of_constraint_code(esc(name), sets)
end


text:
function draw_axes()
    axvline(x=0, color="0.8", lw=2, zorder=1)
    axhline(y=0, color="0.8", lw=2, zorder=1)
end

@doc raw"""
    viz_response(data,
                 plot_title="", plot_xlabel="time, t",
                 plot_ylabel="output, y(t)",
                 savename=nothing)

plot `data[:, :output]` vs. `data[:, :t]` to visualize the response of a system to an input.
typically the data frame, `data`, is returned from [`simulate`](@ref).

note that PyPlot.jl (matplotlib) commands can be invoked after `viz_response` to make further changes to the figure panel.
e.g. `xlim([0, 1])` can be applied after `viz_response`.

# Arguments
* `data::DataFrame`: data frame of time series data, containing a `:t` column for times and `:output` column for the outputs.
* `plot_title::Union{String, LaTeXString}`: title of plot
* `plot_xlabel::Union{String, LaTeXString}`: x-label
* `plot_ylabel::Union{String, LaTeXString}`: y-label
* `savename::Union{Nothing, String}`: filename to save as a figure in .png format (dpi 250).

# Example
```
g = 4 / (4 * s ^ 2 + 0.8 * s + 1)
u = 1 / s
data = simulate(g * u, (0.0, 50.0))
viz_response(data)
```
"""
function viz_response(data::DataFrame;
        plot_title::Union{String, LaTeXString}="",
        plot_xlabel::Union{String, LaTeXString}=L"time, $t$",
        plot_ylabel::Union{String, LaTeXString}=L"output, $y(t)$",
        savename::Union{Nothing, String}=nothing
    )

    figure()
    plot(data[:, :t], data[:, :output], zorder=100)
    xlabel(plot_xlabel)
    ylabel(plot_ylabel)
    title(plot_title)
    draw_axes()
    if ! isnothing(savename)
        tight_layout()
        if ! occursin(".png", savename)
            savename *= ".png"
        end
        savefig(savename, dpi=250, format="png")
    end
    return nothing
end

"""
    viz_poles_and_zeros(tf)

plot the zeros and poles of the transfer function `tf` in the complex plane.
"""
function viz_poles_and_zeros(tf::TransferFunction)
    z, p, k = zeros_poles_gain(tf)

    figure()
    scatter(real.(z), imag.(z), marker="o", label="zeros", color="C1", zorder=100, s=50)
    scatter(real.(p), imag.(p), marker="x", label="poles", color="C2", s=50, zorder=100)
    legend()
    draw_axes()
    xlabel("Re")
    ylabel("Im")
    title("poles and zeros")
    tight_layout()
    return nothing
end

"""
    nyquist_diagram(tf, nb_pts=500, ω_max=10.0)

plot the Nyquist diagram for a transfer function `tf` to visualize its frequency response. `s=-1` is plotted as a red `+`. `nb_pts` changes the resolution. `ω_max` gives maximum frequency considered.
"""
function nyquist_diagram(tf::TransferFunction; nb_pts::Int=500, ω_max::Float64=10.0)
    ω_neg = range(-ω_max, 0.0, length=nb_pts)
    ω_pos = range(0.0, ω_max, length=nb_pts)

    g_iω_neg = [evaluate(tf, ω_i * im) for ω_i in ω_neg]
    g_iω_pos = [evaluate(tf, ω_i * im) for ω_i in ω_pos]

    figure()
    plot(real(g_iω_neg), imag(g_iω_neg), zorder=100)
    plot(real(g_iω_pos), imag(g_iω_pos), zorder=100)
    draw_axes()
    # plot -1
    plot([-1], [0], marker="+", color="r", zorder=1000, markersize=15)
    xlabel("Re[G(iω)]")
    ylabel("Im[G(iω)]")
    title("Nyquist diagram")
    tight_layout()
end

"""
    root_locus(g_ol, max_mag_Kc=10.0, nb_pts=500)

visualize the root locus plot of an open-loop transfer function `g_ol`.

# Arguments
* `g_ol::TransferFunction`: the open-loop transfer function of the closed loop system
* `max_mag_Kc::Float64=10.0`: the maximum magnitude by which the gain of `g_ol` is
    scaled in order to see the roots traversing the plane
* `nb_pts::Int=500`: the number of gains to explore. increase for higher resolution.
"""
function root_locus(g_ol::TransferFunction;
        max_mag_Kc::Float64=10.0, nb_pts::Int=500)
    # compute zeros, poles, and gain of open loop transfer function
    z, p, k = zeros_poles_k(g_ol)

    # decide set of controller gains to explore
    #   make sure Kc is the same sign as gain of G_OL
    #   so G_OL has positive gain
    Kcs = k * 10.0 .^ range(-6, log(max_mag_Kc), length=nb_pts)
    pushfirst!(Kcs, 0.0)

    # store roots of characteristic eqn. here as Kc varies
    rloc = zeros(Complex, length(Kcs), length(p))
    fill!(rloc, NaN)

    # roots of 1 + G_OL when Kc = 0 are poles of G_OL
    rloc[1, :] = p

    # loop thru controller gains, compute roots of 1 + G_OL(s) = 0
    for i = 2:length(Kcs)
        # compute characteristic polynomial 1 + Kc * G_OL(s)
        c_poly = characteristic_polynomial(Kcs[i] * g_ol)
        # compute roots of characteristic polynomial
        roots_c_poly = roots(c_poly)
        # store roots
        entry_filled = [false for _ = 1:length(p)]
        for i_p = 1:length(roots_c_poly)
            ### find closest previous root
            # first, compute distances from each
            distance_to_previous_roots = abs.(rloc[i-1, :] .- roots_c_poly[i_p])
            distance_to_previous_roots[entry_filled] .= Inf
            # get closest
            id_closest_root = argmin(distance_to_previous_roots)
            entry_filled[id_closest_root] = true
            # store root in this entry
            rloc[i, id_closest_root] = roots_c_poly[i_p]
        end
    end

    figure()
    # plot poles; corresponds to Kc = 0
    scatter(real.(p), imag.(p), marker="x", label="poles",
            color="k", s=50, zorder=100)
    # plot zeros; corresponds to |Kc| → ∞
    if length(z) > 0
        scatter(real.(z), imag.(z), marker="o", label="zeros",
                color="k", s=50, zorder=100, facecolor="None")
    end
    # plot roots traversing plane
    for i = 1:length(p)
        plot(real.(rloc[:, i]), imag.(rloc[:, i]),
            zorder=100, color="C$(i-1)")
    end
    xlabel("Re")
    ylabel("Im")
    draw_axes()
    title("root locus")
    legend()
    tight_layout()
end

"""
    axs = bode_plot(tf, log10_ω_min=-4.0, log10_ω_max=4.0, nb_pts=300)

draw the Bode plot of a transfer function `tf` to visualize its frequency response.
returns the two axes of the plot for further tuning via `matplotlib` commands.

adjust the range of frequencies that the Bode plot presents with `log10_ω_min` and `log10_ω_max`.
increase the resolution of the Bode plot with `nb_pts`.
"""
function bode_plot(g::TransferFunction; log10_ω_min::Float64=-3.0, log10_ω_max::Float64=3.0, nb_pts::Int=300)
    ω = 10.0 .^ range(log10_ω_min, log10_ω_max, length=nb_pts)
    g_iω = [evaluate(g, im * ω_i) for ω_i in ω]
    ∠g_iω = zeros(length(g_iω))

    circle_counter = 0
    ∠g_iω[1] = angle(g_iω[1])
    for i = 2:length(g_iω)
        ∠g_iω[i] = angle(g_iω[i]) - circle_counter * 2 * π
        if ∠g_iω[i] - ∠g_iω[i-1] > π
            ∠g_iω[i] -= 2 * π
            circle_counter += 1
        end
    end

    fig, axs = subplots(2, 1, sharex=true, figsize=(8, 7))
    axs[1].plot(ω, abs.(g_iω), color="C0")
    axs[1].set_ylabel(L"$|g(i\omega)|$")
    axs[1].set_title("Bode plot")
    axs[1].set_xscale("log")
    axs[1].set_yscale("log")
    axs[2].set_xscale("log")
    axs[2].plot(ω, ∠g_iω / π, color="C1")
    axs[2].yaxis.set_major_formatter(PyPlot.matplotlib.ticker.FormatStrFormatter(L"%g$\pi$"))
    axs[2].set_ylabel(L"$\angle g(i\omega)$")
    for ax in axs
        ax.minorticks_on()
        ax.grid(b=true, which="minor", alpha=0.25)
    end
    xlabel(L"frequency, $\omega$")
    tight_layout()
    return axs
end

@doc raw"""
    mk_gif(data, plot_title="", plot_xlabel="time, t",
                 plot_ylabel="output, y(t)",
                 savename="response")

make a .gif of the process response.
`data` is a data frame with two columns, `:t` and `:output`, likely returned from [`simulate`](@ref).
accepts same arguments as [`viz_response`](@ref).
ImageMagick must be installed to create the .gif.
the .gif is saved as a file `savename`.

# Arguments
* `data::DataFrame`: data frame of time series data, containing a `:t` column for times and `:output` column for the outputs.
* `plot_title::Union{String, LaTeXString}`: title of plot
* `plot_xlabel::Union{String, LaTeXString}`: x-label
* `plot_ylabel::Union{String, LaTeXString}`: y-label
* `savename::String`: filename to save as a .gif. .gif extension automatically appended if not provided.
"""
function mk_gif(t::Array{Float64}, y::Array{Float64};
        plot_title::Union{String, LaTeXString}="",
        plot_xlabel::Union{String, LaTeXString}=L"time, $t$",
        plot_ylabel::Union{String, LaTeXString}=L"output, $y(t)$",
        savename::String="response.gif"
    )
    if length(t) > 999
        error("too many points and thus images; reduce the number of points")
    end

    # let matplotlib determine x, y lims:
    plot(t, y)
    xmin, xmax, ymin, ymax = axis()
    close()

    # name to save image i as.
    step_to_image(i::Int) = @sprintf("__y_of_t_snapshot_%03d.png", i)

    # save series of images
    for i = 2:length(t)
        viz_response(t[1:i], y[1:i], plot_title=plot_title, plot_xlabel=plot_xlabel,
            plot_ylabel=plot_ylabel)
        xlim([xmin, xmax])
        ylim([ymin, ymax])
        tight_layout()
        savefig(step_to_image(i), format="png", dpi=100)
        close()
    end

    if ! occursin(".gif", savename)
        savename *= ".gif"
    end
    try
        run(`convert -delay 20 -loop 0 __y_of_t_snapshot_\*.png $savename`)
        @info "see " * savename
    catch
        @warn "You must install ImageMagick for `convert` to run. Exiting..."
    end

    # clean up
    for i = 2:length(t)
        rm(step_to_image(i))
    end
end


text:
cd(@__DIR__)
using Pkg; Pkg.activate("."); Pkg.instantiate()

using DiffEqFlux, Flux
using DifferentialEquations
using Plots, DelimitedFiles, Statistics
using Sundials

function FENEP!(out,du,u,p,t,γd)
  # DAE definition for FENE-p
  θ₁₁,θ₂₂,θ₁₂, τ₁₁,τ₂₂,τ₁₂ = u
  λ,η,L = p
  a = L^2 /(L^2 -3)
  fλ = (L^2 + (λ/η/a)*(τ₂₂+τ₁₁))/(L^2 - 3)
  out[1] =  τ₁₁ + du[1] - 2*λ*γd(t)*τ₁₂/fλ
  out[2] =  τ₂₂ + du[2]
  out[3] =  τ₁₂ + du[3] - λ*γd(t)*τ₂₂/fλ - η/fλ * γd(t)

  out[4] = θ₁₁ - λ*τ₁₁/fλ
  out[5] = θ₂₂ - λ*τ₂₂/fλ
  out[6] = θ₁₂ - λ*τ₁₂/fλ
end

function find_σ_exact(tsave,γd)
  # finds the exact solution of the FENE-p equations given a strain rate fn and
  # time points to save at
  λ = 2.0
  L = 2.0
  η = 4.0
  p = [λ,η,L]
  u₀ = zeros(6)
  du₀ = [0.0, 0.0, η*γd(0.0)*(L^2-3)/L^2, 0.0,0.0,0.0]
  tspan = (Float64(tsave[1]),Float64(tsave[end]))
  differential_vars = [true,true,true,false,false,false]
  prob = DAEProblem((out,du,u,p,t) -> FENEP!(out,du,u,p,t,γd),du₀,u₀,tspan,
                    p=p,differential_vars=differential_vars)
  sol = solve(prob,IDA(),saveat=Float64.(tsave))
  return [Float32(σ[6]) for σ in sol.u]
end

function mode_loss(f0,f1,γd)
    # Given the NN f0 and f1 and strain rate γd, this returns a function
    # which when evaluated gives the loss.
    tspan = (0.0f0,6.2831f0)
    tsave = range(tspan[1],tspan[2],length=100)
    σ_exact = find_σ_exact(tsave,γd)
    p = []
    σ0 =0.0f0
    u0 = Tracker.param([σ0])

    dudt_(u::TrackedArray,p,t) = f1(Flux.Tracker.collect(vcat(u, [γd(t)])))
    dudt_(u::AbstractArray,p,t) = Flux.data(f1(vcat(u,[γd(t)])))

    prob = ODEProblem(dudt_,u0,tspan,p)
    loss_rd() =   begin
        P_RD = vcat(Flux.Tracker.collect(diffeq_rd(p,prob,Tsit5(),u0=u0,saveat=tsave))
                      ,γd.(tsave)')
        σ_out = [f0(P_RD[:,i])[1] for i = 1:size(P_RD,2)]
        return sum( (σ_out .- σ_exact).^2 )
    end
    return loss_rd
end

function test_NN(γd,f0,f1)
    # Test the trained NN against a new strain rate
    tspan = (0.0f0,10.0f0)
    tsave = range(tspan[1],tspan[2],length=100)
    σ_exact = find_σ_exact(tsave,γd)
    dudt_(u::AbstractArray,p,t) = Flux.data(f1(vcat(u,[γd(t)])))
    prob = ODEProblem(dudt_,[0.0f0],tspan,p)
    ode_solve = solve(prob,u0=[0.0f0],Tsit5())
    f0_ = Flux.data(f0)
    σ_approx = [f0_([ode_solve(t)[1],γd(t)]).data[1] for t in tsave]

    return tsave, σ_approx, σ_exact
end
function test_err(γd,f0,f1)
    tsave, σ_approx, σ_exact = test_NN(γd,f0,f1)
    return sum( (σ_approx .- σ_exact).^2 )
end

## --- Define NN and linear model ---
f0_n = Chain(Dense(2,4,tanh), Dense(4,1))
f1_n = Chain(Dense(2,4,tanh), Dense(4,1))
f0_l = Chain(Dense(2,1))
f1_l = Chain(Dense(2,1))
# --- Define the total loss across various strain rates ----
t_loss(f0,f1) =  sum(mode_loss(f0,f1,t -> 12.0f0*cos.(ω.*t))() for ω in 1.0f0:0.2f0:2.0f0)
data = Iterators.repeated((), 5000)
opt = ADAM(0.015)
p = []
# --- Define test strain rate here, so call back can see how all errors improve
#     with time ----
γd_test = t -> 12f0*cos.(1.5f0.*t) #+ 7f0*cos.(2.5f0.*t)
err_l = []
err_n = []
cb(err,f0,f1) = push!(err,[t_loss(f0,f1).data,test_err(γd_test,f0,f1)])
Flux.train!(() -> t_loss(f0_n,f1_n), params(f0_n,f1_n,p), data, opt,cb= () -> cb(err_n,f0_n,f1_n))
Flux.train!(() -> t_loss(f0_l,f1_l), params(f0_l,f1_l,p), data, opt,cb= () -> cb(err_l,f0_l,f1_l))

# --- Make plot of how error goes with time ------
Er = (i,er) -> [e[i] for e in er]
p1 = plot(Er(1,err_n),xscale=:log10,yscale=:log10,ylabel="Error",xlabel="Training steps",
        label="Training error, Neural net",legend=:bottomleft)
plot!(p1,Er(2,err_n),label="Testing error, Neural net")
plot!(p1,Er(1,err_l), label="Training error, linear model")
plot!(p1,Er(2,err_l),label="Testing error, linear model")

# --- Plot approximations vs true sol ----
tsave, σ_approx_n, σ_exact = test_NN(γd_test,f0_n,f1_n)
tsave, σ_approx_l, σ_exact = test_NN(γd_test,f0_l,f1_l)
p2 = plot(tsave,σ_approx_n,m=:circle,label="NN solution",ylabel="stress",xlabel="time")
plot!(p2,tsave,σ_approx_l,m=:hexagon,label="Linear model")
plot!(p2,tsave,σ_exact,label="True solution",leg=:bottomright,size=(700,400))

# program takes time to run, so save data for tweaking figures
open("er_data.txt", "w") do io # or other path
           writedlm(io, [Er(1,err_n) Er(2,err_n) Er(1,err_l) Er(2,err_l)], ',')
end;

open("plt_data.txt", "w") do io # or other path
           writedlm(io, [[t for t in tsave]  σ_approx_n σ_approx_l σ_exact], ',')
end;

# Timing
using StaticArrays, BenchmarkTools
W1 = convert(SMatrix{4,2},Tracker.data(f1_n[1].W))
b1 = convert(SVector{4},Tracker.data(f1_n[1].b))
W2 = convert(SMatrix{1,4},Tracker.data(f1_n[2].W))
b2 = convert(SVector{1},Tracker.data(f1_n[2].b))

tspan = (0.0f0,10.0f0)
tsave = range(tspan[1],tspan[2],length=100)
σ_exact = find_σ_exact(tsave,γd_test)
dudt_(u,p,t) = first(W2*broadcast(tanh,W1*@SVector([u[1],γd_test(t)])+b1)+b2)
prob = ODEProblem(dudt_,0.0f0,tspan,p)
@btime ode_solve = solve(prob,Tsit5(),saveat=tsave)
@btime find_σ_exact(tsave,γd_test)


text:
function Base.show(io::IO, r::DataFrameRow)
    labelwidth = mapreduce(n -> length(string(n)), max, _names(r)) + 2
    @printf(io, "DataFrameRow (row %d)", row(r))
    for (label, value) in pairs(r)
        print(io, '\n', rpad(label, labelwidth, ' '), something(value, ""))
    end
end


text:
using BinaryBuilder

version = v"1.1.33"

sources = [
    "https://github.com/GNOME/libxslt/archive/v$(version).tar.gz" =>
        "cb4bca0adbabf218948550bdf1f1b20120e7ae8322c7d599be4b58908946aad8",
]

script = raw"""
cd ${WORKSPACE}/srcdir/libxslt-*
./autogen.sh
./configure --prefix=${prefix} --host=${target} --without-python
make -j${nproc} install
"""

products(prefix) = [
    LibraryProduct(prefix, "libxslt", :libxslt)
    LibraryProduct(prefix, "libexslt", :libexslt)
]

platforms = supported_platforms()

dependencies = [
    "https://github.com/bicycle1885/ZlibBuilder/releases/download/v1.0.4/build_Zlib.v1.2.11.jl",
    "https://github.com/bicycle1885/XML2Builder/releases/download/v1.0.2/build_XML2Builder.v2.9.9.jl"
]

build_tarballs(ARGS, "XSLTBuilder", version, sources, script, platforms, products, dependencies)


text:
## Test models ##
module Test_model

using MultipleTesting
using Test
using Distributions

@testset "Models" begin

    @testset "BUM model" begin

        bum = BetaUniformMixtureModel(0.7)
        bum0 = BetaUniformMixtureModel(1.0) ## uniform, all H₀
        bum1 = BetaUniformMixtureModel(0.0) ## beta, all H₁

        x = collect(0:0.01:1)

        @test cdf.(Ref(bum0), x) ≈ cdf.(Ref(Uniform()), x)
        @test cdf.(Ref(bum1), x) ≈ cdf.(Ref(Beta(0.5, 1.0)), x)

        @test pdf.(Ref(bum0), x) ≈ pdf.(Ref(Uniform()), x)
        @test pdf.(Ref(bum1), x) ≈ pdf.(Ref(Beta(0.5, 1.0)), x)

        @test isin(rand(bum, 100), 0.0, 1.0)
        @test isin(rand(bum0, 100), 0.0, 1.0)
        @test isin(rand(bum1, 100), 0.0, 1.0)

        @test_throws MethodError BetaUniformMixtureModel()
        @test_throws DomainError BetaUniformMixtureModel(-0.2)
        @test_throws DomainError BetaUniformMixtureModel(1.2, 0.2, 2.0)

    end

end

end


text:
using Gomah
using Gomah: chainercv, reversedims, np, @pywith
using Flux
using Statistics
using Test

const PyResNet = chainercv.links.model.resnet.ResNet
const resnet = PyResNet(50, pretrained_model = "imagenet")

@testset "conv1.conv" begin
    dummyX = 128 * rand(Float32, 1, 3, 224, 224)
    pyconv1 = resnet.conv1.conv
    chret = reversedims(pyconv1(dummyX).array)
    flconv1 = ch2conv(pyconv1)
    Flux.testmode!(flconv1, true)
    flret = flconv1(reversedims(dummyX))

    @show pyconv1.ksize
    @show pyconv1.pad
    @show pyconv1.stride

    @test size(flret) == size(chret)
    @show mean(abs.(flret .- chret))
end

@testset "conv1.bn" begin
    dummyX = rand(Float32, (1, 64, 224, 224))
    pybn = resnet.conv1.bn
    flbn = ch2bn(pybn)
    Flux.testmode!(flbn, true)
    flret = flbn(reversedims(dummyX))
    @pywith chainer.using_config("train", false) begin
        chret = reversedims(pybn(dummyX).array)
        @show size(flret), size(chret)
        @test size(flret) == size(chret)
        @test all(isapprox.(flret, chret, atol = 1e-6))
    end
end

@testset "Conv2DBNActiv" begin
    dummyX = rand(Float32, (1, 3, 224, 224))
    pyconv1 = resnet.conv1
    flconv1 = Conv2DBNActiv(pyconv1)
    Flux.testmode!(flconv1, true)
    flret = flconv1(reversedims(dummyX))
    @pywith chainer.using_config("train", false) begin
        chret = reversedims(pyconv1(dummyX).array)
        @test size(flret) == size(chret)
        @show maximum(abs.(flret - chret))
        @test all(isapprox.(flret, chret, atol = 1e-6))
    end
end

@testset "Pool1" begin
    dummyX = rand(Float32, (1, 64, 112, 112))
    pypool1 = resnet.pool1
    flpool1 = MaxPool((3, 3), pad = (0, 1, 0, 1), stride = (2, 2))
    chret = reversedims(pypool1(dummyX).array)
    flret = flpool1(reversedims(dummyX))
    @test size(flret) == size(chret)
    @show maximum(abs.(flret - chret))
    @test all(isapprox.(flret, chret, atol = 1e-4))
end

@testset "BottleNeckA" begin
    dummyX = rand(Float32, (1, 64, 56, 56))
    pyres2A = resnet.res2.a
    flres2A = BottleNeckA(pyres2A)
    Flux.testmode!(flres2A, true)
    @pywith chainer.using_config("train", false) begin
        chret = reversedims(pyres2A(dummyX).array)
        flret = flres2A(reversedims(dummyX))
        @test size(flret) == size(chret)
        @show maximum(abs.(flret - chret))
        @test all(isapprox.(flret, chret, atol = 1e-4))
    end
end

@testset "ResBlock2" begin
    dummyX = rand(Float32, (1, 64, 56, 56))
    pyres2 = resnet.res2
    flres2 = ResBlock(pyres2)
    Flux.testmode!(flres2, true)
    @pywith chainer.using_config("train", false) begin
        chret = reversedims(pyres2(dummyX).array)
        flret = flres2(reversedims(dummyX))
        @test size(flret) == size(chret)
        @show maximum(abs.(flret - chret))
        @test all(isapprox.(flret, chret, atol = 1e-4))
    end
end

@testset "ResBlock3" begin
    dummyX = rand(Float32, (1, 256, 56, 56))
    pyres3 = resnet.res3
    flres3 = ResBlock(pyres3)
    Flux.testmode!(flres3, true)
    @pywith chainer.using_config("train", false) begin
        chret = reversedims(pyres3(dummyX).array)
        flret = flres3(reversedims(dummyX))
        @test size(flret) == size(chret)
        @show maximum(abs.(flret - chret))
        @test all(isapprox.(flret, chret, atol = 1e-4))
    end
end


text:
function instantiate(pat::PatVar, sub::Sub, rule::Rule)
    if haskey(sub, pat)
        (eclass, lit) = sub[pat]
        return (lit !== nothing ? lit : eclass)
    else
        error("unbound pattern variable $pat in rule $rule")
    end
end

function instantiate(pat::PatLiteral{T}, sub::Sub, rule::Rule) where T
    pat.val
end

function instantiate(pat::PatTypeAssertion, sub::Sub, rule::Rule)
    instantiate(pat.var, sub, rule)
end

function instantiate(pat::PatTerm, sub::Sub, rule::Rule)
    # TODO support custom types here!
    # similarterm ? ask Shashi
    meta = pat.metadata
    if meta !== nothing && meta.iscall
        Expr(:call, pat.head, map(x -> instantiate(x, sub, rule), pat.args)...)
    else
        Expr(pat.head, map(x -> instantiate(x, sub, rule), pat.args)...)
    end
end

function apply_rule!(g::EGraph, rule::UnequalRule,
        match::Match, matches::MatchesBuf, rep::Report, mod::Module)
    (_, pat, sub, id) = match
    lc = id
    rinst = instantiate(pat, sub, rule)
    rc = addexpr!(g, rinst).id
    # delete!(matches, match)
    if find(g, lc) == find(g, rc)
        @log "Contradiction!" rule
        rep.reason = Contradiction()
        return (false, rep)
    end
    return (true, nothing)
end

function apply_rule!(g::EGraph, rule::SymbolicRule,
        match::Match, matches::MatchesBuf, rep::Report,  mod::Module)
    (_, pat, sub, id) = match
    rinst = instantiate(pat, sub, rule)
    rc = addexpr!(g, rinst)
    if canprune(typeof(rule)) && rule.prune
        prune!(g, id, rc)
    else
        merge!(g, id, rc.id)
    end
    return (true, nothing)
end

function apply_rule!(g::EGraph, rule::DynamicRule,
        match::Match, matches::MatchesBuf, rep::Report,  mod::Module)
    (_, pat, sub, id) = match
    f = Rules.getrhsfun(rule, mod)
    actual_params = map(rule.patvars) do x
        (eclass, literal) = sub[x]
        literal !== nothing ? literal : eclass
    end
    r = f(geteclass(g, id), g, actual_params...)
    rc = addexpr!(g, r)

    if canprune(typeof(rule)) && rule.prune
        prune!(g, id, rc)
    else
        merge!(g, id, rc.id)
    end
    return (true, nothing)
end

using .ReportReasons
function eqsat_apply!(g::EGraph, matches::MatchesBuf,
        scheduler::AbstractScheduler, rep::Report,
        mod::Module, params::SaturationParams)::Report
    i = 0
    for match ∈ matches
        i += 1

        # if i % 300 == 0
        if params.eclasslimit > 0 && g.numclasses > params.eclasslimit
            @log "E-GRAPH SIZEOUT"
            rep.reason = EClassLimit(params.eclasslimit)
            return rep
        end

        #     if params.stopwhen()
        #         @log "Halting requirement satisfied"
        #         rep.reason = ConditionSatisfied()
        #         return rep
        #     end
        # end

        rule=match[1]
        writestep!(scheduler, rule)
        if find(g, match[4]) ∈ g.pruned
            return rep
        end
        (ok, nrep) = apply_rule!(g, rule, match, matches, rep, mod)
        if !ok
            return nrep
        end
        # println(rule)
        # println(sub)
        # println(l); println(r)
        # display(egraph.classes); println()
    end
    return rep
end


text:
module TestDoctest

import EventTracker
import EventTrackerBase
using Documenter: doctest
using Test

function test_eventtracker()
    doctest(EventTracker; manual = false)
end

function test_eventtrackerbase()
    doctest(EventTrackerBase; manual = false)
end

end  # module


text:
using Test, YaoBlocks, YaoArrayRegister
using YaoBlocks.ConstGate
using StatsBase: mean
using LinearAlgebra: tr

@testset "blockfilter expect" begin
    ghz = (ArrayReg(bit"0000") + ArrayReg(bit"1111")) |> normalize!
    obs1 = kron(4, 2 => Z)
    obs2 = kron(4, 2 => X)
    obs3 = repeat(4, X)
    @test expect(obs1, ghz) ≈ 0
    @test expect(obs2, ghz) ≈ 0
    @test expect(obs3, ghz) ≈ 1

    @test expect(obs3, zero_state(4) => kron(H, H, H, H)) ≈
          expect(obs3, zero_state(4) |> kron(H, H, H, H))
    @test expect(obs1 + obs2 + obs3, ghz) ≈ 1
    @test expect(obs1 + obs2 + obs3, repeat(ghz, 3)) ≈ [1, 1, 1]
    @test expect(2 * obs3, ghz) ≈ 2
    @test expect(2 * obs3, repeat(ghz, 3)) ≈ [2, 2, 2]

    @test blockfilter(
        ishermitian,
        chain(2, kron(2, X, P0), repeat(2, Rx(0), (1, 2)), kron(2, 2 => Rz(0.3))),
    ) == Any[X, P0, kron(2, X, P0), Rx(0), repeat(2, Rx(0), (1, 2))]
    @test blockfilter(
        b -> ishermitian(b) && b isa PrimitiveBlock,
        chain(2, kron(2, X, P0), repeat(2, Rx(0), (1, 2)), kron(2, 2 => Rz(0.3))),
    ) == [X, P0, Rx(0)]
end

@testset "density matrix" begin
    reg = rand_state(4)
    dm = reg |> density_matrix
    op = put(4, 3 => X)
    @test expect(op, dm) ≈ expect(op, reg)

    reg = rand_state(4; nbatch = 3)
    dm = reg |> density_matrix
    op = put(4, 3 => X)
    @test expect(op, dm) ≈ expect(op, reg)
end

@testset "insert_qubits!" begin
    reg = rand_state(5; nbatch = 10)
    insert_qubits!(reg, 3; nqubits = 2)
    @test reg |> nqubits == 7
    @test expect(put(7, 3 => Z), reg) .|> tr |> mean ≈ 1
    @test expect(put(7, 4 => Z), reg) .|> tr |> mean ≈ 1
end

@testset "expect" begin
    for reg in [rand_state(3, nbatch = 10), rand_state(3, nbatch = 10, no_transpose_storage = true)]
        e1 = expect(put(2, 2 => X), reg |> copy |> focus!(1, 2) |> ρ)
        e2 = expect(put(2, 2 => X), reg |> copy |> focus!(1, 2))
        e3 = expect(put(3, 2 => X), reg |> ρ)
        e4 = expect(put(3, 2 => X), reg)
        e5 = expect(put(3, 2 => -X), reg)
        @test e1 ≈ e2
        @test e1 ≈ e3
        @test e1 ≈ e4
        @test e5 ≈ -e4
    end
end

@testset "viewbatch apply" begin
    reg = zero_state(1, nbatch = 2)
    viewbatch(reg, 1) |> X
    @test reg.state[:, 1] == [0, 1]
    @test reg.state[:, 2] == [1, 0]
end

@testset "operator fidelity" begin
    @test operator_fidelity(H, H) ≈ 1
    @test operator_fidelity(im * H, X) ≈ sqrt(0.5)
    @test operator_fidelity(im * H, Z) ≈ sqrt(0.5)
    @test operator_fidelity(H, Y) ≈ 0
end

@testset "gate count, time" begin
    cphase(i, j) = control(i, j => shift(2π / (2^(i - j + 1))))
    hcphases(n, i) = chain(n, i == j ? put(i => H) : cphase(j, i) for j in i:n)
    qft_circuit(n::Int) = chain(n, hcphases(n, i) for i in 1:n)
    qc = qft_circuit(3)
    @test qc |> gatecount |> length == 2
    @test qc |> gatecount |> values |> sum == 6
    res = gatecount(repeat(5, X, (2, 3)))
    @test res |> values |> sum == 2
end


text:
# Copies figure 12 to figures directory:

run(`cp ../static/Figure_MR_v8.png  ..`)


text:
# Julia wrapper for header: /usr/include/libavcodec/dirac.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

export
    av_dirac_parse_sequence_header

function av_dirac_parse_sequence_header(dsh, buf, buf_size::Csize_t, log_ctx)
    ccall((:av_dirac_parse_sequence_header, libavcodec), Cint, (Ptr{Ptr{AVDiracSeqHeader}}, Ptr{UInt8}, Csize_t, Ptr{Cvoid}), dsh, buf, buf_size, log_ctx)
end


text:
using TaylorModels

t = TaylorModel1(3, interval(0), -0.5..0.5)

@show texp = exp(t)

@show tsin = sin(t)

s = texp + tsin
p = texp * tsin

@show TaylorModels.bound_taylor1(t)
@show TaylorModels.bound_taylor1(p)
@show TaylorModels.bound_taylor1(s)

#= Output
texp = exp(t) =  [1, 1] + [1, 1] t + [0.5, 0.5] t² + [0.166666, 0.166667] t³ + [-0, 0.00288794]
tsin = sin(t) =  [1, 1] t + [-0.166667, -0.166666] t³ + [-0.00124851, 0.00124851]
TaylorModels.bound_taylor1(t) = [-0.5, 0.5]
TaylorModels.bound_taylor1(p) = [-0.291667, 0.791667]
TaylorModels.bound_taylor1(s) = [0.124999, 2.12501]
[0.124999, 2.12501]
=#


text:
"""
GAP

"""
mutable struct GAP <: FOSAlgorithm
    α::Float64
    α1::Float64
    α2::Float64
    direct::Bool
    options
end
GAP(α=0.8, α1=1.8, α2=1.8; direct=false, kwargs...) = GAP(α, α1, α2, direct, kwargs)

struct GAPData{T1,T2} <: FOSSolverData
    tmp1::Array{Float64,1}
    tmp2::Array{Float64,1}
    constinit::Base.RefValue{Bool} # If constant term has been calculated (linesearch feature)
    S1::T1
    S2::T2
end

function init_algorithm!(alg::GAP, model::AbstractFOSModel)
    S1, S2, n, status_generator = get_sets_and_status(alg, model)
    data = GAPData(Array{Float64,1}(undef, n), Array{Float64,1}(undef, n),
            Ref(false), S1, S2)
    return data, status_generator
end

function S1constant!(y, alg::GAP, data::GAPData, mem, i)
    if !data.constinit.x
        y .= 0.0
        #TODO high accuracy?
        prox!(mem, S1, y)
        mem .= α1.*mem# .+ (1-α1).*0.0
        data.constinit.x = true
    end
    y .= mem
    return
end

function S1!(y, alg::GAP, data::GAPData, x, i, status::AbstractStatus, longstep=nothing)
    α1, S1 =  alg.α1, data.S1
    prox!(y, S1, x)
    addprojeq(longstep, y, x)
    y .= α1.*y .+ (1-α1).*x
end

function S2!(y, alg::GAP, data::GAPData, x, i, status::AbstractStatus, longstep=nothing)
    α2, S2 =  alg.α2, data.S2
    prox!(y, S2, x)
    checkstatus(status, y)
    addprojineq(longstep, y, x)
    y .= α2.*y .+ (1-α2).*x
end

function Base.step(alg::GAP, data::GAPData, x, i, status::AbstractStatus, longstep=nothing)
    α = alg.α
    tmp1, tmp2 = data.tmp1, data.tmp2
    # Relaxed projection onto S1
    S1!(tmp1, alg, data, x, i, status, longstep)
    # α1, S1 =  alg.α1, data.S1
    # prox!(tmp1, S1, x)
    # addprojeq(longstep, tmp1, x)
    # tmp1 .= α1.*tmp1 .+ (1-α1).*x
    #Relaxed projection onto S2
    S2!(tmp2, alg, data, tmp1, i, status, longstep)
    # α2, S2 =  alg.α2, data.S2
    # prox!(tmp2, S2, tmp1)
    # checkstatus(status, tmp2)
    # addprojineq(longstep, tmp2, tmp1)
    # tmp2 .= α2.*tmp2 .+ (1-α2).*tmp1
    # Relaxation
    x .= α.*tmp2 .+ (1-α).*x
    return
end

function getsol(::GAP, data::GAPData, x)
    tmp1,tmp2,S1,S2 = data.tmp1,data.tmp2,data.S1,data.S2
    prox!(tmp1, S1, x)
    prox!(tmp2, S2, tmp1)
    return tmp2
end

support_linesearch(::GAP) = Val{:Fast}

support_longstep(::GAP) = true
projections_per_step(::GAP) = (1,1)


text:
using PkgTemplates

t = Template(user="danielw2904", # my github username
             dir = pwd(), # directory of package files
             julia=v"1.5.0", # minimum julia version
             plugins=[Documenter{GitHubActions}(), # for documentation
                        GitHubActions(), # for documentation building and testing
                        Git(ssh = true, gpgsign = true) # use ssh in the git repo and sign commits with gpg key
                        ]
)

t("CoolPkg.jl")

using DocumenterTools

# This is necessary for the Documentation to deploy to github pages
DocumenterTools.genkeys(user="danielw2904", repo="CoolPkg.jl")

# Follow instructions:
# 1. Copy first SSH key to link generated
# URL looks like:
# https://github.com/danielw2904/CoolPkg.jl/settings/keys
# Key looks like:
# ssh-rsa ................. Documenter
## -> check "Allow write access"
## Add a title e.g. "Documenter"
# 2. Copy secure environment variable to link generated
# For github it looks like:
# https://github.com/danielw2904/CoolPkg.jl/settings/secrets
# Key are random characters
# Give it the name "DOCUMENTER_KEY"

# Once Actions have finished you will have "dev" docs. "stable" docs are added once you tag a version (either by registering in the JuliaRegistries or manually).


text:
module BaselineUtils

using Curiosity
using GVFHordes

import ..GVFSRHordes
import ..Learner

using Distributions

struct FollowDemon <: Learner
    demon::AbstractGVF
    action_set::UnitRange
    update::Any
    function FollowDemon(demon, action_set)
        new(demon, action_set,1)
    end
end

Curiosity.update!(learner::FollowDemon, args...) = 1

function get_action_probs(l::FollowDemon, state, obs, action)
    #NOTE: This is misleading... state_t for demons is observations
    return get(l.demon.policy, state_t = obs, action_t = action)
end

(l::FollowDemon)(obs) = [get_action_probs(l, nothing, obs, a) for a ∈ l.action_set]

struct RandomDemons <: Learner
    demons::GVFHordes.AbstractHorde
    action_set::UnitRange
    update::Any
    function RandomDemons(demons,action_set)
        new(demons,action_set,1)
    end
end

Curiosity.update!(learner::RandomDemons, args...) = 1

function get_action_probs(l::RandomDemons, state, obs, action)
    #NOTE: This is misleading... state_t for demons is observations
    π = map(gvf -> get(gvf.policy; state_t = obs, action_t = action), l.demons.gvfs)
    return sum(π) / length(l.demons)
end

(l::RandomDemons)(obs) = [get_action_probs(l, nothing, obs, a) for a ∈ l.action_set]

end

max_stars_repo_name: MatthewMcLeod/curiosity

max_issues_repo_name: MatthewMcLeod/curiosity

max_forks_repo_name: MatthewMcLeod/curiosity


text:
ext = ""

deps = joinpath(dirname(@__FILE__),"..","..","deps")

@static if is_unix()
    @static if is_linux() ext = "so" end
    @static if is_apple() ext = "dylib" end

    run(`gcc -fPIC -I$deps/csfml/include -c Window/event.c Network/Network.c Graphics/shader.c`)
    run(`gcc -L$deps -lcsfml-system -lcsfml-graphics -lcsfml-window -lcsfml-audio -lcsfml-network -shared -o $deps/libjuliasfml.$ext event.o Network.o shader.o`)
    run(`rm event.o shader.o Network.o`)
end

@static if is_windows()
    if isdir(Pkg.dir("WinRPM"))
        RPMbindir = Pkg.dir("WinRPM","deps","usr","$(Sys.ARCH)-w64-mingw32","sys-root","mingw","bin")
        ENV["PATH"]=ENV["PATH"]*";"*RPMbindir

        mingw_include = Pkg.dir("WinRPM", "deps", "usr", "$(Sys.ARCH)-w64-mingw32", "sys-root", "mingw", "include")
        deps_dir = joinpath(dirname(@__FILE__),"..","..","deps")
        run(`gcc -I$mingw_include -I$deps_dir\\sfml-binaries\\csfml\\include -c Window\\event.c Network\\Network.c Graphics\\shader.c`)
        run(`gcc -L$deps_dir\\csfml\\bin -lcsfml-network-2 -lcsfml-window-2 -lcsfml-graphics-2 -shared -o $deps_dir\\libjuliasfml.dll event.o shader.o Network.o`)
        rm("event.o"); rm("shader.o"); rm("Network.o")
    else
        println("Please install WinRPM and gcc to build this package")
    end
end


text:
import Pkg

TRAVIS_BUILD_DIR = ENV["TRAVIS_BUILD_DIR"]

pushfirst!(Base.LOAD_PATH,joinpath(TRAVIS_BUILD_DIR, "ci", "RegistryTestingTools",),);

import RegistryTestingTools

RegistryTestingTools.compare_external_registry(
    TRAVIS_BUILD_DIR,
    Pkg.RegistrySpec(
        name = "General",
        uuid = "23338594-aafe-5451-b93e-139f81909106",
        url = "https://github.com/JuliaRegistries/General.git",
        ),
    )


text:
# ---
# title: 63. Unique Paths II
# id: problem63
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Array, Dynamic Programming
# link: <https://leetcode.com/problems/unique-paths-ii/description/>
# hidden: true
# ---
#
# A robot is located at the top-left corner of a `m x n` grid (marked 'Start' in
# the diagram below).
#
# The robot can only move either down or right at any point in time. The robot
# is trying to reach the bottom-right corner of the grid (marked 'Finish' in the
# diagram below).
#
# Now consider if some obstacles are added to the grids. How many unique paths
# would there be?
#
# An obstacle and space is marked as `1` and `0` respectively in the grid.
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2020/11/04/robot1.jpg)
#
#
#
#     Input: obstacleGrid = [[0,0,0],[0,1,0],[0,0,0]]
#     Output: 2
#     Explanation: There is one obstacle in the middle of the 3x3 grid above.
#     There are two ways to reach the bottom-right corner:
#     1. Right -> Right -> Down -> Down
#     2. Down -> Down -> Right -> Right
#
#
# **Example 2:**
#
# ![](https://assets.leetcode.com/uploads/2020/11/04/robot2.jpg)
#
#
#
#     Input: obstacleGrid = [[0,1],[0,0]]
#     Output: 1
#
#
#
#
# **Constraints:**
#
#   * `m == obstacleGrid.length`
#   * `n == obstacleGrid[i].length`
#   * `1 <= m, n <= 100`
#   * `obstacleGrid[i][j]` is `0` or `1`.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
using Test
using LinearAlgebra

include("../src/NoiseOperators.jl")
@testset "Trace for vectorized matrix" begin
    N = 4
    ρ = rand(N, N) + 1im * rand(N,N)
    @test tr(ρ) ≈ trace(ρ[:])
end

@testset "Spre, spost" begin
N = 10
# ρ = rand(N, N) + 1im * rand(N,N)
# A = rand(N, N) + 1im * rand(N,N)

# @test ρ * A' ≈ reshape(sup_post(A) * ρ[:], (N,N))
# @test A * ρ ≈ reshape(sup_pre(A) * ρ[:], (N,N))

ρ = sprand(N, N, .1) + 1im *  sprand(N, N, .1)
A = sprand(N, N, .1) + 1im * sprand(N, N, .1)

@test ρ * A' ≈ reshape(sup_post(A) * ρ[:], (N,N))
@test A * ρ ≈ reshape(sup_pre(A) * ρ[:], (N,N))

@test typeof(sup_post(A)) <: typeof(A)
@test typeof(sup_pre(A)) <: typeof(A)

ρ = rand(N, N) + 1im * rand(N,N)
A = sprand(N, N, .3) + 1im * sprand(N,N, .3)
B = sprand(N, N, .3) + 1im * sprand(N,N, .3)

@test A * ρ * B' ≈ reshape(sup_pre(A) * sup_post(B) * ρ[:], (N,N))
@test A * ρ * B' ≈ reshape(sup_post(B) * sup_pre(A) * ρ[:], (N,N))
@test A * ρ * B' ≈ reshape(sup_pre_post(A, B) * ρ[:], (N, N))
@test A * ρ * A' ≈ reshape(sup_pre_post(A) * ρ[:], (N, N))

end


text:
module Z4Group0s

export Z4Group0, z4group0
import Random

using ..Z4Groups

"""
    struct Z4Group0

Type that can represent 0, +1, -1, +im, -im
and supports sufficient algebra to compute and
represent tensor products of Pauli matrices.
In particular, addition is only supported if
at least one operand is zero. This type is closed
under multiplication and this restricted addition.
"""
struct Z4Group0
    z4::Z4Group
    zero::Bool
end

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

"""
    z4group0(n::Integer)::Z4Group0

Perform one-based index `n` into `(1, -1, i, -i, 0)`.
"""
function z4group0(n::Integer)
    if n == 5
        return zero(Z4Group0)
    end
    return Z4Group0(z4group(n), true)
end

function Random.rand(rng::Random.AbstractRNG, ::Random.SamplerType{Z4Group0})
    return z4group0(rand(rng, 1:5))
end

####
#### Conversion
####

Base.promote_type(::Type{T}, ::Type{Z4Group0}) where {T <: Real} = Complex{T}
Base.promote_type(::Type{Complex{T}}, ::Type{Z4Group0}) where {T <: Real} = Complex{T}

Z4Group0(x::ConvertableNumber) = convert(Z4Group0, x)

function Base.convert(::Type{Z4Group0}, x::ConvertableNumber)
    if iszero(x)
        return zero(Z4Group0)
    end
    if isreal(x)
        if isone(x)
            z = one(Z4Group)
        elseif isone(-x)
            z = -one(Z4Group)
        else
            error("Can't convert $x to {0, 1, -1, i, -i}.")
        end
    elseif iszero(x.re)
        if  isone(x.im)
            z = Z4Group(true, false) # can probably use constant prop here.
        elseif  isone(-x.im)
            z = Z4Group(true, true)
        else
            error("Can't convert $x to {0, 1, -1, i, -i}.")
        end
    else
        error("Can't convert $x to {0, 1, -1, i, -i}.")
    end
    return Z4Group0(z, true)
end

function Base.convert(::Type{T}, x::Z4Group0) where {T<:Number}
    if iszero(x)
        return zero(T)
    end
    return convert(T, x.z4)
end

####
#### IO
####

function Base.show(io::IO, ip::Z4Group0)
    if ! ip.zero
        print(io, "0")
        return nothing
    end
    show(io, ip.z4)
    return nothing
end

####
#### Container interface
####

Base.length(::Z4Group0) = 1
Base.iterate(n::Z4Group0) = (n, nothing)
Base.iterate(n::Z4Group0, ::Nothing) = nothing

####
#### Compare / predicates
####

Base.iszero(x::Z4Group0) = ~(x.zero)

Base.zero(::Type{Z4Group0}) = Z4Group0(one(Z4Group), false)

Base.one(::Type{Z4Group0}) = Z4Group0(one(Z4Group), true)
# Consider making Z4Group0 a Number so that the following need not be defined

Base.:(==)(x::Z4Group0, y::Number) = convert(Complex{Int}, x) == y

Base.isreal(x::Z4Group0) = x.zero && isreal(x.z4)

Base.real(x::Z4Group0) = ~x.zero ? 0 : real(x.z4)

Base.imag(x::Z4Group0) = ~x.zero ? 0 : imag(x.z4)

####
#### Algebra / mathematical operations
####

function Base.:*(xz::Z4Group0, yz::Z4Group0)
    zerox = xz.zero
    zeroy = yz.zero
    zero0 = zerox & zeroy
    if ! zero0
        return zero(Z4Group0)
    end
    x = xz.z4
    y = yz.z4
    sign = zero0 & (~(~(x.minus ⊻ y.minus) ⊻ (x.imag & y.imag)))
    imag = zero0 & (x.imag ⊻ y.imag)
    return Z4Group0(Z4Group(imag, sign), zero0)
end

Base.:*(xz::Z4Group, yz::Z4Group0) = yz * xz

function Base.:*(x::Z4Group0, y::Z4Group)
    if ! x.zero
        return x
    end
    return Z4Group0(x.z4 * y, true)
end

Base.:*(x::Z4Group0, y::Number) = *(promote(y, x)...)

Base.:*(y::Number, x::Z4Group0) = x * y

function Base.:+(x::Z4Group0, y::Z4Group0)
    if iszero(x)
        return y
    elseif iszero(y)
        return x
    else
        throw(DomainError("cant do this."))
    end
end

function Base.:-(x::Z4Group0)
    z4 = Z4Group(x.z4.imag, ~x.z4.minus)
    return Z4Group0(z4, x.zero)
end

end # module Z4Group0s


text:
# julia1_1.4.jl Julia Resources
# GLM Generealized Linear Model
#Pkg.add("Gadfly")
#using Gadfly
#Pkg.add("DataFrames")
#Pkg.add("GLM")

using Gadfly
using DataFrames, GLM

data=DataFrame(x=[0.9,2.0,3.0],y=[2.0,4.0,6.5])
ols=glm(@formula(y ~ x), data, Normal(), IdentityLink())
c=coef(ols)
#  y=c[1]+c[2]*x
#  y=-0.0370091 + 2.13746*x
plot(data,x=:x,y=:y,Geom.point,intercept=[c[1]],slope=[c[2]],Geom.abline(color="red",style=:dash))

using GLM, DataFrames
data=DataFrame(x1=[0.9,2.0,3.0,3.9],x2=[0.291,0.210,0.105,0.021],y=[2.0,4.0,6.5,7.9])
ols=glm(@formula(y ~ x1 + x2), data, Normal(), IdentityLink())

coef(ols)
stderror(ols)
predict(ols)


text:
#Script for building the population pyramid into data
using DataFrames,Plots,JLD2,MAT,LinearAlgebra,Distances,StatsPlots,CSV
using LinearAlgebra:normalize,normalize!
using TransformVariables, DynamicHMC, DynamicHMC.Diagnostics, Distributions
age_cats = ["0-4","5-9","10-14","15-19","20-24","25-29","30-34","35-39","40-44","45-49","50-54","55-59","60-64","65-69","70-74","75-79","80+"]
# age_cats_in
# Age specific rates of hospitalisation if becoming a case
#Estimate from CDC COVID-19 Response Team. Severe Outcomes Among Patients with Coronavirus Disease 2019 (COVID-19) - United States, February 12-March 16, 2020. MMWR Morb. Mortal. Wkly. Rep. 69, 343–346 (2020).
hosp_rate_CDC = [mean([1.6,2.5]),#0-19 year olds
                 mean([14.3,20.8]),#20-44 yos
                 mean([21.2,28.3]),#45-54
                 mean([20.5,30.1]),#55-64
                 mean([28.6,43.5]),#65-74
                 mean([30.5,58.7]),#75-84
                 mean([31.3,70.3])]./100 #85+

ICU_rate_CDC = [0.,#0-19 year olds
              mean([2.,4.2]),#20-44 yos
              mean([5.4,10.4]),#45-54
              mean([4.7,11.2]),#55-64
              mean([8.1,18.8]),#65-74
              mean([10.5,31.0]),#75-84
              mean([6.3,29.])]./100 #85+

death_rate_CDC = [0.,#0-19 year olds
            mean([0.1,0.2]),#20-44 yos
            mean([0.5,0.8]),#45-54
            mean([1.4,2.6]),#55-64
            mean([2.7,4.9]),#65-74
            mean([4.3,10.5]),#75-84
            mean([10.4,27.3])]./100 #85+

hosp_rate_by_age = [hosp_rate_CDC[1],hosp_rate_CDC[1],hosp_rate_CDC[1],hosp_rate_CDC[1],
                    hosp_rate_CDC[2],hosp_rate_CDC[2],hosp_rate_CDC[2],hosp_rate_CDC[2],hosp_rate_CDC[2],
                    hosp_rate_CDC[3],hosp_rate_CDC[3],
                    hosp_rate_CDC[4],hosp_rate_CDC[4],
                    hosp_rate_CDC[5],hosp_rate_CDC[5],
                    hosp_rate_CDC[6],hosp_rate_CDC[7]]

ICU_rate_by_age = [ICU_rate_CDC[1],ICU_rate_CDC[1],ICU_rate_CDC[1],ICU_rate_CDC[1],
                    ICU_rate_CDC[2],ICU_rate_CDC[2],ICU_rate_CDC[2],ICU_rate_CDC[2],ICU_rate_CDC[2],
                    ICU_rate_CDC[3],ICU_rate_CDC[3],
                    ICU_rate_CDC[4],ICU_rate_CDC[4],
                    ICU_rate_CDC[5],ICU_rate_CDC[5],
                    ICU_rate_CDC[6],ICU_rate_CDC[7]]
death_rate_by_age = [death_rate_CDC[1],death_rate_CDC[1],death_rate_CDC[1],death_rate_CDC[1],
                    death_rate_CDC[2],death_rate_CDC[2],death_rate_CDC[2],death_rate_CDC[2],death_rate_CDC[2],
                    death_rate_CDC[3],death_rate_CDC[3],
                    death_rate_CDC[4],death_rate_CDC[4],
                    death_rate_CDC[5],death_rate_CDC[5],
                    death_rate_CDC[6],death_rate_CDC[7]]

ICU_rate_by_age_cond_hosp = ICU_rate_by_age./hosp_rate_by_age
death_rate_by_age_cond_ICU = death_rate_by_age./hosp_rate_by_age
death_rate_by_age_cond_ICU[isnan.(death_rate_by_age_cond_ICU)] .= 0.
# death_rate_by_age_cond_ICU
"""
Sub county population data
"""
# subcounty_population = readtable("data/distribution-of-population-by-age-sex-county-and-sub-county-kenya-2019-census-volume-iii.csv")
# subcounties = unique(subcounty_population.sub_county)
# unique(subcounty_population.Age)
# subC_to_C = Dict{String,String}()
# # for
"""
Load the MCMC outputs for the age-specific relative symptomatic rate
"""
@load "data/MCMC_results_fit_d.jld2" results_d_tau_0 results_d_tau_01 results_d_tau_025 results_d_tau_05 results_d_tau_1
trans = as((θ = as(Array, asℝ₊, 17),))#Transforms the relative values into definitely positive numbers 16th value is fixed as 1
gr()
function MCMCchain2estimates(results)
    posterior = first.(transform.(trans, results.chain))
    posterior_array = zeros(length(posterior),17)
    for i = 1:length(posterior),j=1:17
        posterior_array[i,j] = posterior[i][j]
    end
    for i = 1:length(posterior),j=1:17
        posterior_array[i,j] = posterior_array[i,j]/posterior_array[i,17]
    end
    θ̂ = [mean(posterior_array[:,i]) for i = 1:17]
    lb = [quantile(posterior_array[:,i],0.025) for i = 1:17]
    ub = [quantile(posterior_array[:,i],0.975) for i = 1:17]
    return θ̂,θ̂ .- lb,ub.-θ̂
end

d_0,lb_0,ub_0 = MCMCchain2estimates(results_d_tau_0)
d_01,lb_01,ub_01 = MCMCchain2estimates(results_d_tau_01)
d_025,lb_025,ub_025 = MCMCchain2estimates(results_d_tau_025)
d_05,lb_05,ub_05 = MCMCchain2estimates(results_d_tau_05)
d_1,lb_1,ub_1 = MCMCchain2estimates(results_d_tau_1)

# fig_d = groupedbar(hcat(θ̂_0,θ̂_01,θ̂_025,θ̂_05,θ̂_1),fillalpha = 0.7,yscale = :log10)
# fig_d = scatter(θ̂_0,fillalpha = 0.7,yscale = :log10)

fig_d = scatter(d_0,yerr = (lb_0,ub_0),fillalpha = 0.3,lab="rel. infectiousness = 0%",legend = :bottomright,
        xticks = (1:17,age_cats),
        ylabel = "Relative symptomatic rate", title = "Symptomatic rates for subclinical rel. infectiousness scenarios",
        yscale = :log10,
        ms = 10)
scatter!(fig_d,d_01,yerr = (lb_01,ub_01),fillalpha = 0.3,lab="rel. infectiousness = 10%",ms=10)
scatter!(fig_d,d_025,yerr = (lb_025,ub_025),fillalpha = 0.3,lab="rel. infectiousness = 25%",ms=10)
scatter!(fig_d,d_05,yerr = (lb_05,ub_05),fillalpha = 0.3,lab="rel. infectiousness = 50%",ms=10)
scatter!(fig_d,d_1,yerr = (lb_1,ub_1),fillalpha = 0.3,lab="rel. infectiousness = 100%",ms=10)
scatter!(size=(700,400))
xlabel!("Age group (years)")
savefig("plotting/symptomatic_rates.pdf")

"""
Load age-specific susceptibilities (for use in that scenario) and age-specific relative detection rates
for different relative infectiousness values τ = 0,0.1,0.25,0.5,1
"""
@load "data/susceptibility_rates.jld2" σ
@load "data/detection_rates_for_different_taus.jld2" d_0 d_01 d_025 d_05 d_1
rel_detection_rates = hcat(d_0,d_01,d_025,d_05,d_1)

"""
Load population sizes for each risk region
"""

popsize_risk_region = readtable("data/2019_census_age_pyramids_riskregions.csv")
N_region_age = zeros(Int64,20,17)
for i = 1:20,j = 1:17
    N_region_age[i,j] = round(Int64,popsize_risk_region[i,j+2])
end
heatmap(N_region_age)

"""
Load age mixing matrix for Kenya and convert into JLD2 array, do the same for China for baseline comparisons
"""
prop_75_79_amongst_75_plus_china = 0.7927
prop_75_79_amongst_75_plus_Kenya = sum(N_region_age[:,16])/(sum(N_region_age[:,16:17]))

function extend_and_convert_Prem_matrix(M,proportion)
    M_from_to = zeros(17,17)
    for i = 1:16,j=1:15
        M_from_to[i,j] = M[i,j]
    end
    for i = 1:16
        M_from_to[i,16] = M[i,16]*proportion
        M_from_to[i,17] = M[i,16]*(1-proportion)
    end
    for j = 1:17
        M_from_to[17,j] = M_from_to[16,j]
    end
    return Matrix(M_from_to')
end

# agemixingmatrix_table_china = readtable("data/china_baseline_age_mixing.csv")
# agemixingmatrix_china = zeros(16,16)
# for i = 1:16,j=1:16
#     agemixingmatrix_china[i,j] = agemixingmatrix_table_china[i,j]
# end
# heatmap(agemixingmatrix_china,clims = (0.,2.))
# M_China = extend_and_convert_Prem_matrix(agemixingmatrix_china,prop_75_79_amongst_75_plus_china)
# heatmap(M_China,clims = (0.,2.5))
@load "data/agemixingmatrix_china.jld2" M_China

agemixingmatrix_table = readtable("data/Kenya_age_matrix.csv")
agemixingmatrix = zeros(16,16)
for i = 1:16,j=1:16
    agemixingmatrix[i,j] = agemixingmatrix_table[i,j]
end
M_Kenya = extend_and_convert_Prem_matrix(agemixingmatrix,prop_75_79_amongst_75_plus_Kenya)
heatmap(M_Kenya,clims = (0.,2.))

@save "data/agemixingmatrix_Kenya_norestrictions.jld2" M_Kenya

agemixingmatrix_table = readtable("data/Kenya_age_matrix_home_only.csv")
agemixingmatrix = zeros(16,16)
for i = 1:16,j=1:16
    agemixingmatrix[i,j] = agemixingmatrix_table[i,j]
end
M_Kenya_ho = extend_and_convert_Prem_matrix(agemixingmatrix,prop_75_79_amongst_75_plus_Kenya)
heatmap(M_Kenya_ho,clims = (0.,2.))
@save "data/agemixingmatrix_Kenya_homeonly.jld2" M_Kenya_ho

agemixingmatrix_table = readtable("data/Kenya_age_matrix_other.csv")
agemixingmatrix = zeros(16,16)
for i = 1:16,j=1:16
    agemixingmatrix[i,j] = agemixingmatrix_table[i,j]
end
M_Kenya_other = extend_and_convert_Prem_matrix(agemixingmatrix,prop_75_79_amongst_75_plus_Kenya)

agemixingmatrix_table = readtable("data/Kenya_age_matrix_school_only.csv")
agemixingmatrix = zeros(16,16)
for i = 1:16,j=1:16
    agemixingmatrix[i,j] = agemixingmatrix_table[i,j]
end
M_Kenya_school = extend_and_convert_Prem_matrix(agemixingmatrix,prop_75_79_amongst_75_plus_Kenya)

agemixingmatrix_table = readtable("data/Kenya_age_matrix_work_only.csv")
agemixingmatrix = zeros(16,16)
for i = 1:16,j=1:16
    agemixingmatrix[i,j] = agemixingmatrix_table[i,j]
end
M_Kenya_work = extend_and_convert_Prem_matrix(agemixingmatrix,prop_75_79_amongst_75_plus_Kenya)

@save "data/agemixingmatrix_Kenya_all_types.jld2" M_Kenya M_Kenya_ho M_Kenya_other M_Kenya_school M_Kenya_work
"""
Some heatmap plots
"""
@load "data/agemixingmatrix_Kenya_all_types.jld2" M_Kenya M_Kenya_ho M_Kenya_other M_Kenya_school M_Kenya_work
heatmap(M_Kenya,clims = (0.,1))
plt_w = bar(M_Kenya_work[:,17],xticks = (1:17,age_cats),size = (800,300),title = "work contacts",ylims = (0,0.6))
plt_o = bar(M_Kenya_other[:,17],xticks = (1:17,age_cats),size = (800,300),title = "other contacts",ylims = (0,0.6))
plt_h = bar(M_Kenya_ho[:,17],xticks = (1:17,age_cats),size = (800,300),title = "home contacts",ylims = (0,0.6))
plt_s = bar(M_Kenya_school[:,17],xticks = (1:17,age_cats),size = (800,300),title = "school contacts",ylims = (0,0.6))

plt_before = heatmap(M_Kenya,xticks = (1:17,age_cats),yticks = (1:17,age_cats),size = (800,600),title = "Estimated contacts before social distancing",clims = (0,3),
                    xlabel = "Contact from individual in age group",ylabel = "Contact received by someone in age group")
plt_after = heatmap(1.5*M_Kenya_ho .+ 0.25*M_Kenya_other .+ 0.75*M_Kenya_work,xticks = (1:17,age_cats),yticks = (1:17,age_cats),size = (800,600),title = "Estimated contacts after social distancing",clims = (0,3),
            xlabel = "Contact from individual in age group",ylabel = "Contact received by someone in age group")

plt_before_and_after = groupedbar(hcat(M_Kenya[:,17], 1.5*M_Kenya_ho[:,17] .+ 0.25*M_Kenya_other[:,17] .+ 0.75*M_Kenya_work[:,17]),
                                    xticks = (1:17,age_cats),
                                    size = (800,500),
                                    title = "Estimated contact rate for a 75+ year old before and after social distancing",
                                    lab = ["Before social distancing" "After social distancing"],
                                    xlabel = "Age group 75+ year old is contacting",
                                    ylabel = "Daily number of contacts")
savefig(plt_before,"plotting/social_contacts_before_SD.pdf")
savefig(plt_after,"plotting/social_contacts_after_SD.pdf")
savefig(plt_before_and_after,"plotting/contacts_made_by_oldpeople_before_and_after_SD.pdf")
layout = @layout [a b]
plot(plt_before,plt_after,layout=layout)
plot!(size = (1200,800))

"""
Load movement matrix - then derive the P,ρ and T values
"""
@load "data/mv_matrix.jld2" mv_matrix
movements_per_person = zeros(20,20)
for i = 1:20,j=1:20
    movements_per_person[i,j] = mv_matrix[j,i]/1000#change to from col to row format for ease of multiplication with row vectors
end
heatmap(movements_per_person)
ρ = [sum(movements_per_person,dims = 1)[j]*5/30 for j = 1:20]
P_dest = zeros(20,20)#probability distribution of destination
for j = 1:20
    P_dest[:,j] = normalize(movements_per_person[:,j],1)
end
T = zeros(20,20) # combined location density matrix
for i = 1:20,j=1:20
    if i == j
        T[i,j] = 1-ρ[i]
    else
        T[i,j] = ρ[j]*P_dest[i,j]
    end
end
heatmap(T,clims = (0.,0.1))

"""
Save all the data
"""
@save "data/data_for_age_structuredmodel.jld2" N_region_age M_Kenya movements_per_person P_dest ρ T σ rel_detection_rates M_Kenya_ho hosp_rate_by_age ICU_rate_by_age_cond_hosp


text:
# This file is sometimes included as a standalone script, but we need some of
# the values from bundlepaths.jl
if !(@isdefined BUNDLES_PATH)
    include("./bundlepaths.jl")
end

# Avoid namespace pollution with let
let
    package_dir = dirname(@__DIR__)

    # NodeJS isn't a hard requirement of WebIO, but is needed to build packages,
    # so we need to install it in CI.
    if isci()
        @info "CI detected, installing NodeJS..."

        using Pkg
        Pkg.add("NodeJS")
    end

    # Don't build packages outside of a dev environment (or CI).
    if !isdev()
        @warn(
            "Can't build WebIO JS when not checked out for development. "
            * "Run `Pkg.dev(\"WebIO\")` if you want to build JS."
        )
        return
    end

    # Build the dang packages!
    using NodeJS
    package_dir = normpath(joinpath(@__DIR__, "..", "packages"))
    npm = `$(NodeJS.npm_cmd()) -C $(package_dir)`

    install_cmd = `$npm install --scripts-prepend-node-path=auto --unsafe-perm`
    @info "Installing NPM dependencies..." cmd=install_cmd
    run(install_cmd)

    args = (
        haskey(ENV, "WEBIO_WEBPACK_ARGS")
        ? [ENV["WEBIO_WEBPACK_ARGS"]]
        : []
    )
    build_cmd = `$npm run build-prod --scripts-prepend-node-path=auto --unsafe-perm -- $args`
    @info "Building packages..." cmd=build_cmd
    run(build_cmd)

    if isdir(BUNDLES_PATH)
        rm(BUNDLES_PATH, recursive=true)
    end
    mkdir(BUNDLES_PATH)

    # Copy important things to the right place
    core_bundle_out = joinpath(package_dir, "webio", "dist", "webio.bundle.js")
    if !isfile(core_bundle_out)
        @error "Cannot find WebIO core bundle: $core_bundle_out"
        error("WebIO core bundle was not built properly!")
    end
    @info "Copying $(core_bundle_out) to $(CORE_BUNDLE_PATH)..."
    cp(core_bundle_out, CORE_BUNDLE_PATH; force=true)

    generic_http_bundle_out = joinpath(package_dir, "generic-http-provider", "dist", "generic-http.bundle.js")
    @info "Copying $(generic_http_bundle_out) to $(GENERIC_HTTP_BUNDLE_PATH)..."
    cp(generic_http_bundle_out, GENERIC_HTTP_BUNDLE_PATH; force=true)

    nbextension_bundle_out = joinpath(package_dir, "jupyter-notebook-provider", "dist", "jupyter-notebook.bundle.js")
    @info "Copying $(nbextension_bundle_out) to $(JUPYTER_NBEXTENSION_PATH)..."
    cp(nbextension_bundle_out, JUPYTER_NBEXTENSION_PATH; force=true)

    mux_bundle_out = joinpath(package_dir, "mux-provider", "dist", "mux.bundle.js")
    @info "Copying $(mux_bundle_out) to $(MUX_BUNDLE_PATH)..."
    cp(mux_bundle_out, MUX_BUNDLE_PATH; force=true)

    blink_bundle_out = joinpath(package_dir, "blink-provider", "dist", "blink.bundle.js")
    @info "Copying $(blink_bundle_out) to $(BLINK_BUNDLE_PATH)..."
    cp(blink_bundle_out, BLINK_BUNDLE_PATH; force=true)
end


text:
function PointVectorValuePair(arg0::Vector{jdouble}, arg1::Vector{jdouble})
    return PointVectorValuePair((Vector{jdouble}, Vector{jdouble}), arg0, arg1)
end

function PointVectorValuePair(arg0::Vector{jdouble}, arg1::Vector{jdouble}, arg2::jboolean)
    return PointVectorValuePair((Vector{jdouble}, Vector{jdouble}, jboolean), arg0, arg1, arg2)
end

function get_point(obj::PointVectorValuePair)
    return jcall(obj, "getPoint", Vector{jdouble}, ())
end

function get_point_ref(obj::PointVectorValuePair)
    return jcall(obj, "getPointRef", Vector{jdouble}, ())
end

function get_value(obj::PointVectorValuePair)
    return jcall(obj, "getValue", Vector{jdouble}, ())
end

function get_value_ref(obj::PointVectorValuePair)
    return jcall(obj, "getValueRef", Vector{jdouble}, ())
end


text:
# logic for default rcopy

"""
`rcopy` copies the contents of an R object into a corresponding canonical Julia type.
"""
rcopy(s::SymSxpPtr) = rcopy(Symbol,s)
rcopy(s::CharSxpPtr) = rcopy(Compat.String,s)

function rcopy(s::StrSxpPtr)
    if anyna(s)
        rcopy(DataArray,s)
    elseif length(s) == 1
        rcopy(Compat.String,s)
    else
        rcopy(Array{Compat.String},s)
    end
end
function rcopy(s::RealSxpPtr)
    if anyna(s)
        rcopy(DataArray{Float64},s)
    elseif length(s) == 1
        rcopy(Float64,s)
    else
        rcopy(Array{Float64},s)
    end
end
function rcopy(s::CplxSxpPtr)
    if anyna(s)
        rcopy(DataArray{Complex128},s)
    elseif length(s) == 1
        rcopy(Complex128,s)
    else
        rcopy(Array{Complex128},s)
    end
end
function rcopy(s::LglSxpPtr)
    if anyna(s)
        rcopy(DataArray{Bool},s)
    elseif length(s) == 1
        rcopy(Bool,s)
    else
        rcopy(BitArray,s)
    end
end
function rcopy(s::IntSxpPtr)
    if isFactor(s)
        rcopy(PooledDataArray,s)
    elseif anyna(s)
        rcopy(DataArray{Int},s)
    elseif length(s) == 1
        rcopy(Cint,s)
    else
        rcopy(Array{Cint},s)
    end
end

function rcopy(s::VecSxpPtr)
    if isFrame(s)
        rcopy(DataFrame,s)
    elseif isnull(getnames(s))
        rcopy(Array{Any},s)
    else
        rcopy(Dict{Symbol,Any},s)
    end
end

rcopy(s::FunctionSxpPtr) = rcopy(Function,s)

# TODO
rcopy(l::LangSxpPtr) = l
rcopy(r::RObject{LangSxp}) = r


text:
blockbandwidths(A::AdjOrTrans) = reverse(blockbandwidths(parent(A)))
subblockbandwidths(A::AdjOrTrans) = reverse(subblockbandwidths(parent(A)))


text:
__precompile__()

module MessyTimeSeries

    # Libraries
    using Dates, Distributed, Logging;
    using LinearAlgebra, Distributions, StableRNGs, Statistics;

    # Custom dependencies
    local_path = dirname(@__FILE__);
    include("$(local_path)/types.jl");
    include("$(local_path)/methods.jl");
    include("$(local_path)/kalman.jl");
    include("$(local_path)/subsampling.jl");

    # Export types
    export IntVector, IntMatrix, IntArray, FloatVector, FloatMatrix, FloatArray, JVector, JMatrix, JArray, DiagMatrix, SymMatrix,
           KalmanSettings, KalmanStatus, OnlineKalmanStatus, DynamicKalmanStatus, SizedKalmanStatus, ARIMASettings, VARIMASettings;

    # Export methods
    export check_bounds, nan_to_missing!, error_info, verb_message, interpolate_series, forward_backwards_rw_interpolation, centred_moving_average,
            soft_thresholding, solve_discrete_lyapunov, isconverged, trimmed_mean, sum_skipmissing, mean_skipmissing, std_skipmissing, is_vector_in_matrix,
            demean, standardise, lag, companion_form, no_combinations, rand_without_replacement, get_bounded_log, get_unbounded_log, get_bounded_logit, get_unbounded_logit;

    # Export functions
    export kfilter!, kfilter_full_sample, kfilter_full_sample!, kforecast, ksmoother, fmin_uc_models, arima, varima, forecast,
           block_jackknife, optimal_d, artificial_jackknife, moving_block_bootstrap, stationary_block_bootstrap;
end


text:
function convert_arguments(LS::Type{<:LineSegments}, mesh::PolytopeVertexMesh)
    return convert_arguments(LS, materialize_edges(mesh))
end

function convert_arguments(M::Type{<:GLMakieMesh}, mesh::PolytopeVertexMesh)
    return convert_arguments(M, materialize_faces(mesh))
end

max_stars_repo_name: khurrumsaleem/MOCNeutronTransport

max_issues_repo_name: khurrumsaleem/MOCNeutronTransport

max_forks_repo_name: khurrumsaleem/MOCNeutronTransport


text:
"""
    @JShape(body, args...)

Macro that creates a user-defined shape based on `Luxor` instructions within its `begin...end` block

# Example

```julia
Object(@JShape begin
    sethue(color)
    poly(points, action, close = true)
end action = :stroke color = "red" radius = 8)
```

In this example, the arguments after the `end` (i.e. `action = :stroke color = "red" radius = 8`)
can be used inside the `begin...end` block and animated using the [`change`](@ref) action.
"""
macro JShape(body::Expr, expr...)
    leftside = Expr(
        :tuple,
        Expr(:parameters, [Expr(:kw, e.args[1], e.args[2]) for e in expr]...),
        Expr(:(...), :args),
    )
    out = Expr(:(->), leftside, body)
    esc(out)
end


text:
"""
    isPPT(x, sys :: Int, dims :: Vector) :: Bool

This function returns true if the input state x is PPT
with respect to the (sys)th system. False otherwise.
dims is a vector of the sizes of the subsystems.
"""
function isPPT(x,sys::Int,dims::Vector) :: Bool
    #We don't make sure the inputs are good because the partialtranspose function
    #will do that for us

    #Get PPT of x
    xPPT = partialtranspose(x,sys,dims)

    #Check if it is PSD
    if eigmin(xPPT) >= 0 || isapprox(0,eigmin(xPPT))
        return true
    else
        return false
    end
end

"""
    swapOperator(dim :: Int) :: Matrix{Float64}

This function is the swap operator ``\\mathbb{F}`` which is defined by the action
```math
\\mathbb{F}(u \\otimes v) = v \\otimes u \\hspace{5mm} u,v \\in \\mathcal{H}_{A} .
```
The function uses that ``\\mathbb{F} = \\sum_{a,b \\in \\Sigma} E_{a,b} \\otimes E_{b,a}``
where ``E_{a,b}`` is a square matrix of dimension ``\\Sigma`` with a one in the ``(a,b)``
entry and a ``0`` everywhere else.
"""
function swapOperator(dim :: Int) :: Matrix{Float64}
    swap_operator = zeros(dim^2,dim^2)

    for col_id in 1:dim^2
        # factoring col_id into Alice and Bob subsystem ids
        a_id = floor(Int64, (col_id-1)/dim)
        b_id = (col_id-1) % dim

        # swapping Alice and Bob subsystem ids to construct row_id
        row_id = b_id*dim + a_id + 1
        swap_operator[row_id, col_id] = 1
    end

    return swap_operator
end

"""
    permuteSubsystems(
        ρ:: Vector,
        perm::Vector{Int64},
        dims::Vector{Int64}
    ) :: Vector

This function returns the vector with the subsystems permuted. For example, given three
subspaces ``A,B,C``, and the permutation ``\\pi`` defined by ``(A,B,C) \\xrightarrow[]{\\pi} (C,A,B),``
the function implements the process:
```math
    |e_{i}\\rangle_{A} |e_j \\rangle_{B} |e_k \\rangle_{C} \\xrightarrow[]{\\pi}
    |e_{k} \\rangle_{C} |e_{i}\\rangle_{A}  |e_{j} \\rangle_{B}  ,
```
by re-indexing the vector, permuting the indices appropriately, and converting it
back into a vector.
"""
function permuteSubsystems(ρ:: Vector, perm::Vector{Int64},dims::Vector{Int64}) :: Vector
    #This is almost identical to Tony Cubitt's implementation of this function https://www.dr-qubit.org/matlab.html
    #This is largely because Julia does reshape column-wise as Matlab does
    orig_shape = size(ρ)
    num_subsys = length(perm)
    #Note certain things get reversed. This is because Julia does reshape column-wise
    permTup = Tuple((num_subsys+1) .- reverse(perm)) #Reshape requires tuples
    dimTup = Tuple(reverse(dims))
    result = reshape(permutedims(reshape(ρ,dimTup),permTup),orig_shape)
end
"""
    permuteSubsystems(
        ρ:: Matrix,
        perm::Vector{Int64},
        dims::Vector{Int64}
    ) :: Matrix

This function returns the matrix with the subsystems permuted. It is a generalization of the vector code.
For example, given three subspaces ``A,B,C``, and the permutation ``\\pi`` defined by ``(A,B,C) \\xrightarrow[]{\\pi} (C,A,B),``
the function implements the process:
```math
|e_{i}\\rangle\\langle e_{i}|_{A}\\otimes |e_j \\rangle\\langle e_j|_{B}\\otimes |e_k \\rangle\\langle e_k|_{C} \\xrightarrow[]{\\pi}
|e_{k} \\rangle\\langle e_{k}|_{C}\\otimes |e_{i}\\rangle\\langle e_{i}|_{A} \\otimes  |e_{j} \\rangle\\langle e_{j}|_{B} ,
```
by re-indexing the matrix, permuting the indices, and reconstructing the matrix.
Both bra and ket indices receive the same permutation.
"""
function permuteSubsystems(ρ:: Matrix,perm::Vector{Int64},dims::Vector{Int64}) :: Matrix
    #This is almost identical to Tony Cubitt's implementation of this function https://www.dr-qubit.org/matlab.html
    #This is largely because Julia does reshape column-wise as Matlab does
    if !isequal(size(ρ)...)
            throw(DomainError(ρ, "the input ρ is not a square matrix or a pure state"))
    end
    orig_shape = size(ρ)
    num_subsys = length(perm)
    #For the matrix version we do the same thing we just have twice as many indices
    #because we keep track of bras and kets
    permPrime = (num_subsys+1) .- reverse(perm)
    permTup = Tuple([permPrime num_subsys .+ permPrime])
    dimTup = Tuple([reverse(dims) reverse(dims)])
    result = reshape(permutedims(reshape(ρ,dimTup),permTup),orig_shape)

    return result
end
"""
    shiftOperator(d::Int64) :: Matrix

This function returns the operator that shifts the computational basis mod d
for complex Euclidean space of dimension d. That is, it returns the operator ``S``
defined by the action
```math
    S|k\\rangle = |k+1 \\mod d \\rangle
```
"""
function shiftOperator(d::Int64) :: Matrix
    id = Matrix(I,d,d)
    shift = id[[end,1:(end-1)...],:]

    return shift
end

"""
    discreteWeylOperator(m :: Int64, n :: Int64, d :: Int64)

This function returns the (m,n)^th unitary for generating the generalized
Bell basis. They are defined by their action on the computational basis:
```math
    U_{m,n}|e_{k}\\rangle = e^{2 \\pi mk i / d} |e_{k+n}\\rangle
```
These are the discrete Weyl Operator basis. See
[Section 3 of this paper](https://arxiv.org/abs/1004.1655) for further details.
!!! warning
    Different works define the discrete Weyl operators differently such that
    the phase or ordering may differ. Please check your reference.
"""
function discreteWeylOperator(m :: Int64, n :: Int64, d :: Int64) :: Matrix
    if m < 0 || n < 0
        throw(DomainError((n,m), "Make sure m,n ∈ [0,1,...,d-1]."))
    elseif m >= d || n >= d
        throw(DomainError((n,m), "Make sure m,n < d."))
    end
    λ = exp(2*pi*1im / d)
    U = zeros(ComplexF64,d,d)
    for k = 0 : d-1
        U[((k+n) % d) + 1,k+1] = λ^(k*m)
    end
    return U
end

max_forks_repo_name: ChitambarLab/cv-channel


text:
const PathLike = Union{AbstractString, AbstractPath}

######################### utilities ##############################

# Apply the combine operation on the value of two Nodes
# sometimes it maybe lazy. If one of the inputs is lazy,
# a new Thunk is set as the value of the resulting node
function apply_combine(f, x, y)
    (!hasvalue(x) && !hasvalue(y)) && return y
    try
        setvalue(x, maybe_lazy(f)(x[], y[]))
    catch err
        !(err isa UndefMergeError) && rethrow(err)
        error("$(name(x)) clashed with an existing file name at path $(path(x)).\n" *
              "Pass `combine=f` to define how to combine them.")
    end
end

# attach a node file at `path`; use `combine` to overwrite.
function attach(t, path::PathLike, t′; combine=_merge_error)
    spath = Path(path).segments
    t1 = foldl((x, acc) -> acc => [x], [t′; reverse(spath)...]) |> maketree
    merge(t, maketree(name(t)=>[t1]); combine=combine)
end

# rewrite `tree` according by performing string search and replace on every
# path use `combine` to overwrite multiple files which map to the same path
function regex_rewrite_tree(tree, from_path, to_path, combine)
    newtree = maketree(name(tree)=>[])
    # Exclude the root from the matching because
    #   1) the public API where the pattern is relative to the root and
    #   2) attach wants a path relative to the root
    rootoffset = length(canonical_path(path(tree)))+2
    for x in Leaves(tree)
        newname = replace(canonical_path(path(x))[rootoffset:end], from_path => to_path)
        dir = dirname(newname)
        dir = isempty(dir) ? "." : dir
        newtree = attach(newtree,
                         dir,
                         rename(x, basename(newname));
                         combine=combine)
    end
    newtree
end

# getindex but return the rooted tree
function _getsubtree(x, path::PathLike)
    _getsubtree(x, GlobMatch(string(path)))
end
_getsubtree(x, path) = x[path]

"""
A wrapper around function `f` which causes overwriter
functions (`combine` option) to get the nodes themselves
rather than the value stored in them.

Note that if the nodes have lazy values, the values will be Thunks.
"""
struct OnNodes
    f
end

apply_combine(f::OnNodes, x, y) = f(x, y)

function _combine(cs, combine)
    seen = Dict()
    out = []
    for c in cs
        prev = get(seen, name(c), nothing)
        if prev !== nothing
            out[end] = apply_combine(combine, c, out[end])
        else
            push!(out, c)
            seen[name(c)] = c
        end
    end
    map(identity, out)
end

# flatten folders named  "." into the parent directory.
function normdots(x::FileTree; combine=_merge_error)
    c2 = map(children(x)) do y
        z=normdots(y; combine=combine)
        name(z) == "." ? children(z) : [z]
    end |> Iterators.flatten |> collect
    FileTree(x; children=_combine(c2, combine)) |> setparent
end

normdots(x::File; kw...) = x

################################ rm ################################

"""
    diff(t1::FileTree, t2::FileTree)

For each node in `t2` remove a node in `t1` at the same path if it exists.
Returns the difference tree.
"""
function Base.diff(t1::FileTree, t2::FileTree)
    if name(t1) == name(t2)
        t2_names = name.(children(t2))
        cs = []
        for x in children(t1)
            idx = findfirst(==(name(x)), t2_names)
            if !isnothing(idx)
                if t2[idx] isa File
                    @assert x isa File
                elseif x isa FileTree && t2[idx] isa FileTree
                    d = diff(x, t2[idx])
                    if !isempty(d)
                        push!(cs, d)
                    end
                end
            else
                push!(cs, x)
            end
        end
        FileTree(t1; children=cs) |> setparent
    else
        t1
    end
end

"""
    rm(t::FileTree, pattern::Union{Glob, String, AbstractPath, Regex})

remove nodes which match `pattern` from the file tree.
"""
rm(t::FileTree, path) = diff(t, _getsubtree(t, path))

################################ cp ################################

"""
    cp(t::FileTree,
       from_path::Regex,
       to_path::SubstitutionString; combine)

copy nodes in the file tree whose path matches the `from_tree` regular expression pattern
by renaming it to `to_path` pattern. Any sub-pattern in `from_path` which is surrounded by
paranthesis will be read as a matched substring which can be accessed in the to_path
substitution pattern using \\1, \\2 etc. positional matches.

If a file overwrites an existing node after copy, `combine` will be called to combine
them together. By default `combine` will error.

## Example:

```julia
julia> t = maketree("dir" => [string(j) => [string(i)=>["data.csv"] for i = 1:2] for j=1:2])
dir/
├─ 1/
│  ├─ 1/
│  │  └─ data.csv
│  └─ 2/
│     └─ data.csv
└─ 2/
   ├─ 1/
   │  └─ data.csv
   └─ 2/
      └─ data.csv

julia> cp(t, r"^([^/]*)/([^/]*)/data.csv\$", s"\1/\2.csv")
dir/
├─ 1/
│  ├─ 1/
│  │  └─ data.csv
│  ├─ 1.csv
│  ├─ 2/
│  │  └─ data.csv
│  └─ 2.csv
└─ 2/
   ├─ 1/
   │  └─ data.csv
   ├─ 1.csv
   ├─ 2/
   │  └─ data.csv
   └─ 2.csv
```
"""
function cp(t::FileTree, from_path::Regex, to_path::SubstitutionString; combine=_merge_error)
    matches = t[from_path]
    isempty(matches) && return t

    newtree = regex_rewrite_tree(matches, from_path, to_path, combine)
    merge(t, newtree; combine=combine)
end

"""
    merge(t1::FileTree, t2::FileTree; combine)

Merge two FileTrees. If files at the same path contain values, the `combine` callback will
be called with their values to result in a new value.

If one of the dirs does not have a value, its corresponding argument will be `NoValue()`
If any of the values is lazy, the output value is lazy as well.
"""
function Base.merge(t1::FileTree, t2::FileTree; combine=_merge_error, dotnorm=true)
    bigt = if name(t1) == name(t2)
        t2_names = name.(children(t2))
        t2_merged = zeros(Bool, length(t2_names))
        cs = []
        for x in children(t1)
            idx = findfirst(==(name(x)), t2_names)
            if !isnothing(idx)
                y = t2[idx]
                if y isa FileTree
                    push!(cs, merge(x, y; combine=combine, dotnorm=false))
                else
                    push!(cs, apply_combine(combine, x, y))
                end
                t2_merged[idx] = true
            else
                push!(cs, x)
            end
        end
        FileTree(t1; children=vcat(cs, children(t2)[map(!, t2_merged)]))
    else
        FileTree(nothing, ".", [t1, t2], NoValue())
    end |> setparent
    dotnorm ? normdots(bigt; combine=combine) : bigt
end

################################ mv ################################

"""
    mv(t::FileTree,
       from_path::Regex,
       to_path::SubstitutionString; combine)

move nodes in the file tree whose path matches the `from_tree` regular expression pattern
by renaming it to `to_path` pattern. Any sub-pattern in `from_path` which is surrounded by
paranthesis will be read as a matched substring which can be accessed in the to_path
substitution pattern using \\1, \\2 etc. positional matches.

If a file overwrites an existing node after copy, `combine` will be called to combine
them together. By default `combine` will error.

## Example:

```julia
julia> t = maketree("dir" => [string(j) =>
                                [string(i)=>["data.csv"]
                                    for i = 1:2] for j=1:2])
dir/
├─ 1/
│  ├─ 1/
│  │  └─ data.csv
│  └─ 2/
│     └─ data.csv
└─ 2/
   ├─ 1/
   │  └─ data.csv
   └─ 2/
      └─ data.csv

julia> mv(t, r"^([^/]*)/([^/]*)/data.csv\$", s"\\1/\\2.csv")
dir/
├─ 1/
│  ├─ 1.csv
│  └─ 2.csv
└─ 2/
   ├─ 1.csv
   └─ 2.csv
```
"""
function mv(t::FileTree, from_path::Regex, to_path::SubstitutionString; combine=_merge_error)
    matches = t[from_path]
    isempty(matches) && return t

    newtree = regex_rewrite_tree(matches, from_path, to_path, combine)
    merge(diff(t, matches), newtree; combine=combine)
end
function Base.merge(x::Node, y::Node; combine=_merge_error)
    name(x) == name(y) ? apply_combine(combine, x, y) : FileTree(nothing, ".", [x,y], NoValue())
end

################################# touch #############################

function _mknode(T, t, path::PathLike, value)
    spath = Path(path).segments
    subdir = if T == File
        T(nothing, spath[end], value)
    else
        T(nothing, spath[end], [], value)
    end

    if length(spath) == 1
        merge(t, maketree(name(t) => [subdir]))
    else
        p = joinpath(spath[1:end-1]...)
        attach(t, p, subdir)
    end
end

"""
    touch(t::FileTree, path::AbstractString; value)

Create an file node at `path` in the tree. Does not contain any value by default.
"""
function touch(t::FileTree, path::PathLike; value=NoValue())
    _mknode(File, t, path, value)
end

################################# mkpath #############################

"""
    mkpath(t::FileTree, path::AbstractString; value)

Create a directory node at `path` in the tree. Does not contain any value by default.
"""
function mkpath(t::FileTree, path::PathLike; value=NoValue())
    _mknode(FileTree, t, path, value)
end

struct UndefMergeError end
_merge_error(x, y) = throw(UndefMergeError())

################################# clip #############################

"""
    clip(t, n; combine)

Remove `n` top-level directories. `combine` will be called
to merge any nodes with equal names found at any level being clipped.
"""
function clip(t, n; combine=_merge_error)
    n==0 && return t
    cs = map(children(t)) do x
        y = clip(x, n-1)
    end
    reduce((x,y)->merge(x,y,combine=combine), cs) |> setparent
end