text:
module CSVFiles

using TextParse, IteratorInterfaceExtensions, TableTraits, TableTraitsUtils,
    DataValues, FileIO, HTTP, TableShowUtils, CodecZlib
import IterableTables

export load, save, File, @format_str

struct CSVFile
    filename::String
    delim
    keywords
end

struct CSVStream
    io
    delim
    keywords
end

function Base.show(io::IO, source::CSVFile)
    TableShowUtils.printtable(io, getiterator(source), "CSV file")
end

function Base.show(io::IO, ::MIME"text/html", source::CSVFile)
    TableShowUtils.printHTMLtable(io, getiterator(source))
end

Base.showable(::MIME"text/html", source::CSVFile) = true

function Base.show(io::IO, ::MIME"application/vnd.dataresource+json", source::CSVFile)
    TableShowUtils.printdataresource(io, getiterator(source))
end

Base.showable(::MIME"application/vnd.dataresource+json", source::CSVFile) = true

function Base.show(io::IO, source::CSVStream)
    TableShowUtils.printtable(io, getiterator(source), "CSV file")
end

function Base.show(io::IO, ::MIME"text/html", source::CSVStream)
    TableShowUtils.printHTMLtable(io, getiterator(source))
end

Base.showable(::MIME"text/html", source::CSVStream) = true

function Base.show(io::IO, ::MIME"application/vnd.dataresource+json", source::CSVStream)
    TableShowUtils.printdataresource(io, getiterator(source))
end

Base.showable(::MIME"application/vnd.dataresource+json", source::CSVStream) = true

function fileio_load(f::FileIO.File{FileIO.format"CSV"}, deprecated_delim=nothing; delim=deprecated_delim===nothing ? ',' : deprecated_delim, args...)
    if deprecated_delim!==nothing
        deprecated_delim!=delim && error("deprecated_delim and delim can not both be used at the same time.")
        Base.depwarn("The positional `delim` keyword in the `load` function is deprecated. Instead use the keyword argument `delim`.", :CSVFiles)
    end

    return CSVFile(f.filename, delim, args)
end

function fileio_load(f::FileIO.File{FileIO.format"TSV"}, deprecated_delim=nothing; delim=deprecated_delim===nothing ? '\t' : deprecated_delim, args...)
    if deprecated_delim!==nothing
        deprecated_delim!=delim && error("deprecated_delim and delim can not both be used at the same time.")
        Base.depwarn("The positional `delim` keyword in the `load` function is deprecated. Instead use the keyword argument `delim`.", :CSVFiles)
    end

    return CSVFile(f.filename, delim, args)
end

function fileio_load(s::FileIO.Stream{FileIO.format"CSV"}, deprecated_delim=nothing; delim=deprecated_delim===nothing ? ',' : deprecated_delim, args...)
    if deprecated_delim!==nothing
        deprecated_delim!=delim && error("deprecated_delim and delim can not both be used at the same time.")
        Base.depwarn("The positional `delim` keyword in the `load` function is deprecated. Instead use the keyword argument `delim`.", :CSVFiles)
    end

    return CSVStream(s.io, delim, args)
end

function fileio_load(s::FileIO.Stream{FileIO.format"TSV"}, deprecated_delim=nothing; delim=deprecated_delim===nothing ? '\t' : deprecated_delim, args...)
    if deprecated_delim!==nothing
        deprecated_delim!=delim && error("deprecated_delim and delim can not both be used at the same time.")
        Base.depwarn("The positional `delim` keyword in the `load` function is deprecated. Instead use the keyword argument `delim`.", :CSVFiles)
    end

    return CSVStream(s.io, delim, args)
end

IteratorInterfaceExtensions.isiterable(x::CSVFile) = true
TableTraits.isiterabletable(x::CSVFile) = true
TableTraits.supports_get_columns_copy_using_missing(x::CSVFile) = true

IteratorInterfaceExtensions.isiterable(x::CSVStream) = true
TableTraits.isiterabletable(x::CSVStream) = true
TableTraits.supports_get_columns_copy_using_missing(x::CSVStream) = true

function _loaddata(file)
    if startswith(file.filename, "https://") || startswith(file.filename, "http://")
        response = HTTP.get(file.filename)
        data = String(response.body)
        return TextParse._csvread(data, file.delim; stringarraytype=Array, file.keywords...)
    else
        return csvread(file.filename, file.delim; stringarraytype=Array, file.keywords...)
    end
end

function IteratorInterfaceExtensions.getiterator(file::CSVFile)
    res = _loaddata(file)

    it = TableTraitsUtils.create_tableiterator([i for i in res[1]], [Symbol(i) for i in res[2]])

    return it
end

function TableTraits.get_columns_copy_using_missing(file::CSVFile)
    columns, colnames = _loaddata(file)
     return NamedTuple{(Symbol.(colnames)...,), Tuple{typeof.(columns)...}}((columns...,))
end

function IteratorInterfaceExtensions.getiterator(s::CSVStream)
    res = TextParse.csvread(s.io, s.delim; stringarraytype=Array, s.keywords...)

    it = TableTraitsUtils.create_tableiterator([i for i in res[1]], [Symbol(i) for i in res[2]])

    return it
end

function TableTraits.get_columns_copy_using_missing(s::CSVStream)
    columns, colnames = TextParse.csvread(s.io, s.delim; stringarraytype=Array, s.keywords...)
    return NamedTuple{(Symbol.(colnames)...,), Tuple{typeof.(columns)...}}((columns...,))
end

function Base.collect(x::CSVFile)
    return collect(getiterator(x))
end

function Base.collect(x::CSVStream)
    return collect(getiterator(x))
end

include("csv_writer.jl")

end # module


text:
URL_BASE = "https://www.philadelphiafed.org/-/media/research-and-data/real-time-center/survey-of-professional-forecasters/data-files/files/"

function url_individual(survey)
  return URL_BASE * "individual_" * survey * ".xlsx"
end

function url_mean(survey, type)
  if type == "level"
    what = "level"
  elseif type == "growth"
    what = "growth"
  elseif type in ["pct", "pct_chg"]
    what = "pct_chg"
  elseif type in ["avg_chg", "ave", "ave_chg"]
    what = "ave_chg"
  else
    error(type, " not implemented yet")
  end
  url = URL_BASE * "mean_" * survey * "_" * what * ".xlsx"
  return url
end

function url_median(survey, type)
  if type == "level"
    what = "level"
  elseif type == "pct"
    what = "pct_chg"
  elseif type == "ave"
    what = "ave_chg"
  else
    error(type, " not implemented yet")
  end
  url = URL_BASE * "median_" * survey * "_" * what * ".xlsx"
  return url
end

function url_dispersion(survey)
  return URL_BASE * "dispersion_" * survey * ".xlsx"
end

ext: jl


text:
# ------------------------------------
# read noumber of nodes
# ------------------------------------
function read_nodenum(skipnum)
    """
    xmax: number of nodes in xi direction
    ymax: number of nodes in eta direction
    """
    fff = []
    open("grid/nodesnum", "r") do f
        fff = read(f,String)
    end
    fff = split(fff,"\n",keepempty=false)
    num_nodes = length(fff) - skipnum

    for i in 1+skipnum:length(fff)
        fff[i] = replace(fff[i]," \r" => "")
    end

    temp = split(fff[2]," ")
    xmax = parse(Int64,temp[1])
    ymax = parse(Int64,temp[2])

    return xmax, ymax
end

# ------------------------------------
# read nodes
# ------------------------------------
function read_nodes(skipnum,xmax,ymax)
    fff = []
    open("grid/nodes", "r") do f
        fff = read(f,String)
    end
    fff = split(fff,"\n",keepempty=false)
    num_nodes = length(fff) - skipnum

    for i in 1+skipnum:length(fff)
        fff[i] = replace(fff[i]," \r" => "")
    end

    nodes = zeros(xmax,ymax,2)
    for i in 1:num_nodes
        temp = split(fff[i+skipnum]," ")

        xnum = parse(Int64,temp[1])
        ynum = parse(Int64,temp[2])
        nodes[xnum,ynum,1] = parse(Float64,temp[3])
        nodes[xnum,ynum,2] = parse(Float64,temp[4])
    end
    return nodes
end

# ------------------------------------
# read vecAx and vecAy
# ------------------------------------
function read_vecA(skipnum,xmax,ymax)
    fff = []
    open("grid/vecAx", "r") do f
        fff = read(f,String)
    end
    fff = split(fff,"\n",keepempty=false)
    num_cell = length(fff) - skipnum

    for i in 1+skipnum:length(fff)
        fff[i] = replace(fff[i]," \r" => "")
    end

    vecAx = zeros(xmax, ymax-1, 2)
    for i in 1:num_cell
        temp = split(fff[i+skipnum]," ")

        x = parse(Int64,temp[1])
        y = parse(Int64,temp[2])
        vecAx[x,y,1] = parse(Float64,temp[3])
        vecAx[x,y,2] = parse(Float64,temp[4])
    end

    fff = []
    open("grid/vecAy", "r") do f
        fff = read(f,String)
    end
    fff = split(fff,"\n",keepempty=false)
    num_cell = length(fff) - skipnum

    for i in 1+skipnum:length(fff)
        fff[i] = replace(fff[i]," \r" => "")
    end

    vecAy = zeros(xmax-1, ymax, 2)
    for i in 1:num_cell
        temp = split(fff[i+skipnum]," ")

        x = parse(Int64,temp[1])
        y = parse(Int64,temp[2])
        vecAy[x,y,1] = parse(Float64,temp[3])
        vecAy[x,y,2] = parse(Float64,temp[4])
    end
    return vecAx, vecAy
end

function read_allgrid()
    skip = 1
    xmax, ymax   = read_nodenum(skip)
    nodes        = read_nodes(skip, xmax, ymax)
    vecAx, vecAy = read_vecA(skip, xmax, ymax)
    println("fin reading grid")
    return xmax, ymax, nodes, vecAx, vecAy
end


text:
function encode(msg::Vector{UInt8})
    tuptuo = [0x00]
    lastzero = 0x01
    for x in reverse(msg)
        if iszero(x)
            push!(tuptuo, lastzero)
            lastzero = 0x00
        else
            push!(tuptuo, x)
        end
        lastzero += 0x01
    end
    push!(tuptuo, lastzero)
    return reverse(tuptuo)
end
encode(x::UInt8) = encode([x])
encode(output, msg) = write(output, encode(msg)...)

next!(x) = read(x, UInt8)

function decode!(msg, input, n)
    code = next!(input)
    iszero(code) && return decode!(msg, input, n)
    i = 1
    while true
        code - 1 > n - i + 1 && return decode!(msg, input, n)
        for _ in 1:code - 1
            tmp = next!(input)
            iszero(tmp) && return decode!(msg, input, n)
            msg[i] = tmp
            i += 1
        end
        code = next!(input)
        if iszero(code)
            if i == n + 1
                return msg
            else
                return decode!(msg, input, n)
            end
        else
            if i > n
                return decode!(msg, input, n)
            else
                msg[i] = 0x00
                i += 1
            end
        end
    end
end


text:
using Espresso: matchex
using LinearAlgebra: I
using MacroTools: @capture
using InteractiveUtils: subtypes

import Base: sort

# ========================
# @map macro
# ========================

"""
    map(ex, args)

Return an instance of the map type corresponding to the given expression.

### Input

- `ex`   -- an expression defining the map, in the form of an anonymous function
- `args` -- additional optional arguments

### Output

A map that best matches the given expression.

### Examples

Let us first create a linear map using this macro:

```jldoctest
julia> @map x -> [1 0; 0 0]*x
LinearMap{Int64,Array{Int64,2}}([1 0; 0 0])
```

We can create an affine system as well:

```jldoctest
julia> @map x -> [1 0; 0 0]*x + [2, 0]
AffineMap{Int64,Array{Int64,2},Array{Int64,1}}([1 0; 0 0], [2, 0])
```

Additional arguments can be passed to `@map` using the function-call form, i.e.
separating the arguments by commas, and using parentheses around the macro call.
For example, an identity map of dimension 5 can be defined as:

```jldoctest
julia> @map(x -> x, dim=5)
IdentityMap(5)
```
A state constraint on such map can be specified passing the additional argument
`x ∈ X`.

An identity map can alternatively be created by giving a the size of the identity
matrix as `I(n)`, for example:

```jldoctest
julia> @map x -> I(5)*x
IdentityMap(5)
```
"""
macro map(ex, args...)
    dimension = nothing
    x = (ex.args)[1]
    rhs = (ex.args)[2].args[2]
    if x == rhs
        # identity map, eg: @map(x -> x, dim=2)
        if @capture(args[1], (dim = dim_) | (dim: dim_) )
            dimension = dim
        else
            throw(ArgumentError("cannot parse dimension of identity map"))
        end

        # constrained identity map: @map(x -> x, dim=2, x ∈ X)
        if length(args) > 1 && @capture(args[2], x ∈ X_)
            return Expr(:call, ConstrainedIdentityMap, esc(:($(dimension))), esc(:($(X))))
        else
            return Expr(:call, IdentityMap, esc(:($(dimension))))
        end
    end

    throw(ArgumentError("unable to match the given expression to a known map type"))
end

macro map(ex)
    x = (ex.args)[1]
    rhs = (ex.args)[2].args[2]

    # x -> I(n)*x
    # (this rule is more specific than x -> Ax so it should come before it)
    MT = IdentityMap
    pat = Meta.parse("I(_n) * $x")
    matched = matchex(pat, rhs)
    matched != nothing &&
        return Expr(:call, :($MT), esc(:($(matched[:_n]))))

    # x -> Ax
    MT = LinearMap
    pat = Meta.parse("_A * $x")
    matched = matchex(pat, rhs)
    matched != nothing &&
        return Expr(:call, :($MT), esc(:($(matched[:_A]))))

    # x -> Ax + b
    MT = AffineMap
    pat = Meta.parse("_A * $x + _b")
    matched = matchex(pat, rhs)
    matched != nothing &&
        return Expr(:call, :($MT), esc(:($(matched[:_A]))), esc(:($(matched[:_b]))))

    throw(ArgumentError("unable to match the given expression to a known map type"))
end

# ========================
# @system macro
# ========================

const TYPES_CONTINUOUS = subtypes(AbstractContinuousSystem)
const FIELDS_CONTINUOUS = fieldnames.(TYPES_CONTINUOUS)

const TYPES_DISCRETE = subtypes(AbstractDiscreteSystem)
const FIELDS_DISCRETE = fieldnames.(TYPES_DISCRETE)

types_table(::Type{AbstractContinuousSystem}) = TYPES_CONTINUOUS
types_table(::Type{AbstractDiscreteSystem}) = TYPES_DISCRETE

fields_table(::Type{AbstractContinuousSystem}) = FIELDS_CONTINUOUS
fields_table(::Type{AbstractDiscreteSystem}) = FIELDS_DISCRETE

"""
    _corresponding_type(AT::Type{<:AbstractSystem}, fields::Tuple)

Return the system type whose field names match those in `fields`.

### Input

- `AT`     -- abstract system type, which can be either `AbstractContinuousSystem`
              or `AbstractDiscreSystem`
- `fields` -- tuple of field names

### Output

The system type (either discrete or continous, depending on `AT`) whose fields
names correspond to those in `fields`, or an error if the `fields` do not match
any known system type.

### Examples

```jldoctest
julia> using MathematicalSystems: _corresponding_type

julia> _corresponding_type(AbstractContinuousSystem, ((:A),))
LinearContinuousSystem

julia> _corresponding_type(AbstractContinuousSystem, ((:A), (:B), (:X), (:U)))
ConstrainedLinearControlContinuousSystem
```
"""
function _corresponding_type(AT::Type{<:AbstractSystem}, fields::Tuple)
    TYPES = types_table(AT)
    FIELDS = fields_table(AT)
    idx = findall(x -> issetequal(fields, x), FIELDS)
    if isempty(idx)
        throw(ArgumentError("the entry $(fields) does not match a " *
                            "`MathematicalSystems.jl` structure"))
    end
    @assert length(idx) == 1 "found more than one compatible system type"
    return TYPES[idx][1]
end

"""
    _capture_dim(expr)

Return the tuple containing the dimension(s) in `expr`.

### Input

- `expr` -- symbolic expression that can be of any of the following forms:

    - `:x` or `:(x)` -- state dimension
    - `:(x, u)`      -- state and input dimension
    - `:(x, u, w)`   -- state, input and noise dimensions

### Output

- The scalar `x` if `expr` specifies the state dimension.
- The vector `[x, u]` if `expr` specifies state and input dimension.
- The vector `[x, u, w]` if `expr` specifies state, input and noise dimensions.
"""
function _capture_dim(expr)
    # the order of the `if` clauses is important, as e.g.
    # `@capture(expr, (x_, u_, w_))` is a particular case of `@capture(expr, (x_))`.
    if @capture(expr, (x_, u_, w_))
        dims = [x, u, w]
    elseif @capture(expr, (x_, u_))
        dims = [x, u]
    elseif @capture(expr, (x_))
        dims = x
    else
        throw(ArgumentError("the dimensions in expression $expr could not be parsed; " *
            "see the documentation for valid examples"))
    end
    return dims
end

"""
    is_equation(expr)

Return `true` if the given expression `expr` corresponds to an equation `lhs = rhs`
and `false` otherwise. This function just detects the presence of the symbol `=`.

### Input

- `expr` -- expression

### Output

A `Bool` indicating whether `expr` is an equation or not.
"""
function is_equation(expr)
    return @capture(expr, lhs_ = rhs_)
end

# returns `(equation, AT, state)` if `expr` if the given expression `expr`
# corresponds to a dynamic equation, either of the form `x⁺ = rhs` or `x' = rhs`
# in the former case `AT = AbstractDiscreteSystem` and in the latter
# `AT = AbstractContinuousSystem`. `equation` is `x = rhs` where the superscript
# is stripped from the expression and the state corresponds to `x`,
# otherwise, return `(nothing, nothing, nothing)`
# Additionally, the `lhs` can be in the form `E*x'` or `E*x⁺` where `E` can be
# any variable name.
function strip_dynamic_equation(expr)
    expr_str = string(expr)
    stripped_equation = ""
    if occursin("'", expr_str)
        AT = AbstractContinuousSystem
        stripped_equation =  replace(expr_str, "'" => "")
    elseif occursin("⁺", expr_str)
        AT = AbstractDiscreteSystem
        stripped_equation =  replace(expr_str, "⁺" => "")
    end
    stripped_expr = Meta.parse(stripped_equation)
    # check if `stripped_expr` is an equation and extract `lhs` and `rhs` if so
    !@capture(stripped_expr, lhs_ = rhs_) && return (nothing, nothing, nothing)

    # extract the name of the state variable from the lhs
    if @capture(lhs, (E_*x_)) || @capture(lhs, (x_))
        state = x
        return (stripped_expr, AT, state)
    end
    return (nothing, nothing, nothing)
end

function _parse_system(expr::Expr)
    return _parse_system((expr,))
end

function _parse_system(exprs::NTuple{N, Expr}) where {N}
    # define default dynamic equation, unknown abstract system type,
    # and empty list of constraints
    dynamic_equation = nothing
    AT = nothing
    constraints = Vector{Expr}()

    # define defaults for state, noise and input symbol and default dimension
    default_input_var = :u
    default_noise_var = :w
    state_var = nothing
    input_var = nothing
    noise_var = nothing
    dimension = nothing
    initial_state = nothing # for initial-value problems

    # main loop to parse the subexpressions in exprs
    for ex in exprs

        if is_equation(ex)  # parse an equation
            (stripped, abstract_system_type, subject) = strip_dynamic_equation(ex)
            if subject != nothing
                dynamic_equation = stripped
                state_var = subject
                AT = abstract_system_type
                # if the stripped system has the structure x_ = A_*x_ + B_*u_ or
                # one of the other patterns, handle u_ as input variable
                if @capture(stripped, (x_ = A_*x_ + B_*u_) |
                                      (x_ = x_ + B_*u_) |
                                      (x_ = A_*x_ + B_*u_ + c_) |
                                      (x_ = x_ + B_*u_ + c_) |
                                      (x_ = f_(x_, u_)) )
                    if (f == :+) || (f == :-) || (f == :*)
                        # pattern x_ = f_(x_, u_) also catches the cases:
                        # x_ = x_ + u_, x_ = x_ - u_ and x_ = x_*u_
                        # where u_ doesn't necessarily need to be the input
                    else
                        input_var = u
                    end
                end

            elseif @capture(ex, (dim = (f_dims_)) | (dims = (f_dims_)))
                dimension = _capture_dim(f_dims)

            elseif @capture(ex, (input = u_) | (u_ = input))
                input_var = u

            elseif @capture(ex, (noise = w_) | (w_ = noise))
                noise_var = w

            else
                throw(ArgumentError("could not properly parse the equation $ex; " *
                                    "see the documentation for valid examples"))
            end

        elseif @capture(ex, x_(0) ∈ X0_)
            # TODO? handle equality, || @capture(ex, x_(0) = X0_)
            if x != state_var
                throw(ArgumentError("the initial state assignment, $x(0), does "*
                "not correspond to the state variable $state_var"))
            end
            initial_state = X0

        elseif @capture(ex, state_ ∈ Set_)  # parse a constraint
            push!(constraints, ex)

        elseif @capture(ex, (input: u_) | (u_: input))  # parse an input symbol
            input_var = u

        elseif @capture(ex, (noise: w_) | (w_: noise))  # parse a noise symbol
            noise_var = w

        elseif @capture(ex, (dim: (f_dims_)) | (dims: (f_dims_)))  # parse a dimension
            dimension = _capture_dim(f_dims)

        else
            throw(ArgumentError("the expression $ex could not be parsed; " *
                                "see the documentation for valid examples"))
        end
    end

    # error handling for the given equations
    dynamic_equation == nothing && throw(ArgumentError("the dynamic equation was not found"))

    # error handling for the given set constraints
    nsets = length(constraints)
    nsets > 3 && throw(ArgumentError("cannot parse $nsets set constraints"))

    # error handling for variable names
    state_var == nothing && throw(ArgumentError("the state variable was not found"))
    got_input_var = input_var != nothing
    got_noise_var = noise_var != nothing
    if got_input_var && (state_var == input_var)
         throw(ArgumentError("state and input variables have the same name `$(state_var)`"))
    elseif got_noise_var && (state_var == noise_var)
         throw(ArgumentError("state and noise variables have the same name `$(state_var)`"))
    elseif got_input_var && got_noise_var && (input_var == noise_var)
         throw(ArgumentError("input and noise variables have the same name `$(input_var)`"))
    end

    # assign default values
    if !got_input_var
        input_var = default_input_var
    end
    if !got_noise_var
        noise_var = default_noise_var
    end

    return dynamic_equation, AT, constraints,
           state_var, input_var, noise_var, dimension, initial_state
end

"""
    extract_dyn_equation_parameters(equation, state, input, noise, dim, AT)

Extract the value and field parameter from the dynamic equation `equation` according
to the variable `state`, `input` and `noise`.

For the right-hand side of the dynamic equation, this function returns a vector of tuples
containing some elements from the list
- `(:A_user, :A)`
- `(:B_user, :B)`
- `(:c_user, :c)`
- `(:D_user, :D)`
- `(:f_user, :f)`
- `(:statedim_user :statedim)`
- `(:inputdim_user :inputdim)`
- `(:noisedim_user :noisedim)`
and for the left-hand side, it returns either an empty vector `Any[]` or
`[(:E_user, :E)]` where the first argument of the tuple corresponds to the value
and the second argument of the tuple corresponds to the field parameter.

### Input

- `equation` -- dynamic equation
- `state`    -- state variable
- `input`    -- input variable
- `noise`    -- noise variable
- `dim`      -- dimensionality
- `AT`       -- abstract system type

### Output

Two arrays of tuples containing the value and field parameters for the right-hand
and left-hand side of the dynamic equation `equation`.
"""
function extract_dyn_equation_parameters(equation, state, input, noise, dim, AT)
    @capture(equation, lhs_ = rhscode_)
    lhs_params = Vector{Tuple{Any, Symbol}}()
    rhs_params = Vector{Tuple{Any, Symbol}}()
    # if a * is used on the lhs, the rhs is a code-block
    if  @capture(lhs, E_*x_)
        push!(lhs_params, (E, :E))
        rhs = rhscode.args[end]
    else
        rhs = rhscode
    end

    # if rhs is parsed as addition => affine system which is controlled, noisy or both
    if @capture(rhs, A_ + B__)
        # parse summands of rhs and add * if needed
        summands = add_asterisk.([A, B...], Ref(state), Ref(input), Ref(noise))
        push!(rhs_params, extract_sum(summands, state, input, noise)...)
    # if rhs is a function call except `*` or `-` => black-box system
    elseif @capture(rhs, f_(a__)) && f != :(*) && f != :(-)
        # the dimension argument needs to be a iterable
        (dim == nothing) && throw(ArgumentError("for a blackbox system, the dimension has to be defined"))
        dim_vec = [dim...]
        push!(rhs_params, extract_blackbox_parameter(rhs, dim_vec)...)

    # if rhs is a single term => affine systm (e.g. A*x, Ax, 2x, x or 0)
    else
        # if not specified, assume dim = 1
        dim = (dim == nothing) ? 1 : dim
        if rhs == state  # => rhs = x
            if AT == AbstractDiscreteSystem
                push!(rhs_params, (dim, :statedim))
            elseif AT == AbstractContinuousSystem
                push!(rhs_params, (I(dim), :A))
            end
        elseif rhs == :(0) && AT == AbstractContinuousSystem # x' = 0
            push!(rhs_params, (dim, :statedim))
        else
            if @capture(rhs, -var_) # => rhs = -x
                if state == var
                    push!(rhs_params, (-1.0*I(dim), :A))
                end
            else
                rhs = add_asterisk(rhs, state, input, noise)
                if @capture(rhs, array_ * var_)
                    if state == var # rhs = A_x_ or rhs= A_*x_
                        value = tryparse(Float64, string(array))
                        if value == nothing # e.g. => rhs = Ax
                            push!(rhs_params, (array, :A))
                        else # => e.g., rhs = 2x
                            push!(rhs_params, (value*I(dim), :A))
                        end
                    else
                        throw(ArgumentError("if there is only one term on the "*
                                  "right-hand side, the state needs to be contained"))
                    end
                end
            end
        end
    end
    return lhs_params, rhs_params
end

"""
    add_asterisk(summand, state::Symbol, input::Symbol, noise::Symbol)

Checks if expression `summand` contains `state`, `input` or `noise` at its end.
If so, a multiplication expression, e.g. `Expr(:call, :*, :A, :x) is created.
If not, `summand` is returned.

### Input

- `summand` -- expressions
- `state`   -- state variable
- `input`   -- input variable
- `noise`   -- noise variable

### Output

Multiplication expression or symbol.

### Example

```jldoctest
julia> using MathematicalSystems: add_asterisk

julia> add_asterisk(:(A1*x), :x, :u, :w)
:(A1 * x)

julia> add_asterisk(:(c1), :x, :u, :w)
:c1

julia> add_asterisk(:(Ax1), :x1, :u, :w)
:(A * x1)

julia> add_asterisk(:(Awb), :x1, :u, :wb)
:(A * wb)

julia> add_asterisk(:(A1u), :x, :u, :w)
:(A1 * u)

julia> add_asterisk(:(A1ub), :x, :u, :w)
:A1ub
```
"""
function add_asterisk(summand, state::Symbol, input::Symbol, noise::Symbol)
    if @capture(summand, A_ * x_)
        return summand
    end

    str = string(summand)
    statestr = string(state); lenstate = length(statestr)
    inputstr = string(input); leninput = length(inputstr)
    noisestr = string(noise); lennoise = length(noisestr)

    # if summand contains the state, input or noise variable at the end and has
    # one or more additional characters (i.e. length(str) > length(state) for the state)
    # a `*` is added in between
    if lenstate < length(str) && str[(end-lenstate+1):end] == statestr
        return Meta.parse(str[1:end-length(statestr)]*"*$state")

    elseif leninput < length(str) && str[(end-leninput+1):end] == inputstr
        return Meta.parse(str[1:end-length(inputstr)]*"*$input")

    elseif lennoise < length(str) && str[(end-lennoise+1):end] == noisestr
        return Meta.parse(str[1:end-length(noisestr)]*"*$noise")

    else # summand is returned which is either a constant term or the state, input or noise variable
        return summand
    end
end

"""
    extract_sum(summands, state::Symbol, input::Symbol, noise::Symbol)

Extract the variable name and field name for every element of `summands` which
corresponds to the elements of the right-hand side of an affine system.

### Input

- `summands` -- array of expressions
- `state`    -- state variable
- `input`    -- input variable
- `noise`    -- noise variable

### Output

Array of tuples of symbols with variable name and field name.

### Example

```jldoctest
julia> using MathematicalSystems: extract_sum

julia> extract_sum([:(A1*x)], :x, :u, :w)
1-element Array{Tuple{Any,Symbol},1}:
 (:(hcat(A1)), :A)

julia> extract_sum([:(A1*x), :(B1*u), :c], :x, :u, :w)
3-element Array{Tuple{Any,Symbol},1}:
 (:(hcat(A1)), :A)
 (:(hcat(B1)), :B)
 (:(vcat(c)), :c)

julia> extract_sum([:(A1*x7), :( B1*u7), :( B2*w7)], :x7, :u7, :w7)
3-element Array{Tuple{Any,Symbol},1}:
 (:(hcat(A1)), :A)
 (:(hcat(B1)), :B)
 (:(hcat(B2)), :D)
```

### Notes

If an element of `summands` is a multiplication expression `lhs*rhs`,
return `lhs` as variable name and `:A` as field name if `rhs==state`,
`:B` as field name if `rhs==input` and `:D` as field name if `rhs==noise`.

If an element of `summands` is a symbol, and not equal to `input` or `noise`,
the symbol is the variable name and the field name is `:c`. If it is equal to
`input`, the variable name is a `IdentityMultiple(I,state_dim)` where `state_dim`
is extracted from the state matrix (i.e. take the symbol `lhs` of `lhs*rhs` where
`rhs==state` which corresponds to the state matrix and generate the expression
`state_dim = size(lhs,1)` which is evaluated in the scope where `@system` is called)
and the field name is `:B`.

Similiarily, if the element is equal to `noise`, the variable name is
`IdentityMultiple(I, state_dim)` and the field name is `:D`.
"""
function extract_sum(summands, state::Symbol, input::Symbol, noise::Symbol)
    params = Tuple{Any,Symbol}[]
    state_dim = 1
    got_state_dim = false

    num_state_assignments = 0
    num_input_assignments = 0
    num_noise_assignments = 0
    num_const_assignments = 0

    for summand in summands
        if @capture(summand, array_ * var_)
            if state == var
                push!(params, (Expr(:call, :hcat, array), :A))
                # obtain "state_dim" for later using in IdentityMultiple
                state_dim =  Expr(:call, :size, :($array), 1)
                got_state_dim = true
                num_state_assignments += 1

            elseif input == var
                push!(params, (Expr(:call, :hcat, array), :B))
                 num_input_assignments += 1

            elseif noise == var
                push!(params, (Expr(:call, :hcat, array), :D))
                num_noise_assignments += 1

            else
                throw(ArgumentError("in the dynamic equation, the expression "*
                "$summand does not contain the state $state, the input $input "*
                "or the noise term $noise"))
            end
        elseif @capture(summand, array_)
            identity = :(I($state_dim))
            # if array == variable: field value equals identity
            if state == array
                push!(params, (identity, :A))
                num_state_assignments += 1
            elseif input == array
                push!(params, (identity, :B))
                num_input_assignments += 1
            elseif noise == array
                push!(params, (identity, :D))
                num_noise_assignments += 1
            else
                push!(params, (Expr(:call, :vcat, array), :c))
                num_const_assignments += 1
            end
        end
    end
    num_const_assignments > 1 && throw(ArgumentError("there is more than one constant term"))
    num_state_assignments > 1 && throw(ArgumentError("there is more than one state term `$state`"))
    num_input_assignments > 1 && throw(ArgumentError("there is more than one input term `$input`"))
    num_noise_assignments > 1 && throw(ArgumentError("there is more than one noise term `$noise`"))

    return params
end

# Extract the variable name of the function of the rhs of an equation if it has
# the form `f_(x_)`, `f_(x_,u_)` or `f_(x_,u_,w_)`.
# In addition to the variable name of the function and the field name `:f`,
# depending on the number of input arguments, the field names `:statedim`, `:inputdim`
# and `:noisedim` with corresponding variable names (which are the number of state,
# input and noise dimensions provided as input to the macro) is returned.
function extract_blackbox_parameter(rhs, dim::AbstractVector)
    if @capture(rhs, f_(x_))
        @assert length(dim) == 1
        return [(f, :f), (dim[1], :statedim)]
    elseif @capture(rhs, f_(x_,u_))
        @assert length(dim) == 2
        return [(f, :f), (dim[1], :statedim),
                         (dim[2], :inputdim)]
    elseif @capture(rhs, f_(x_,u_,w_))
        @assert length(dim) == 3
        return [(f, :f), (dim[1], :statedim),
                         (dim[2], :inputdim),
                         (dim[3], :noisedim)]
    end
end

# Take the vectors of tuples providing the variable and fields names for the
# lhs, the rhs and the sets, combine them, and return a vector of field names
# and a vector of variable names.
function constructor_input(lhs, rhs, set)
    rhs_fields = [tuple[2] for tuple in rhs]
    lhs_fields = [tuple[2] for tuple in lhs]
    set_fields = [tuple[2] for tuple in set]
    (length(unique(set_fields)) != length(set_fields)) &&
        throw(ArgumentError("There is some ambiguity in the set definition"))
    rhs_var_names = [tuple[1] for tuple in rhs]
    lhs_var_names = [tuple[1] for tuple in lhs]
    set_var_names = [tuple[1] for tuple in set]
    field_names = (rhs_fields..., lhs_fields..., set_fields...)
    var_names = (rhs_var_names..., lhs_var_names..., set_var_names...)
    return field_names, var_names
end

# extract the variable name and the field name for a set expression. The method
# checks wether the set belongs to the `state`, `input` or `noise` and returns a
# tuple of symbols where the field name is either `:X`, `:U` or `:W` and
# the variable name is the value parsed as Set_
function extract_set_parameter(expr, state, input, noise) # input => to check set definitions
    if @capture(expr, x_(0) ∈ Set_)
        return Set, :x0

    elseif @capture(expr, x_ ∈ Set_)
        if x == state
            return  Set, :X
        elseif x == input
            return  Set, :U
        elseif x == noise
            return  Set, :W
        else
            error("$expr is not a valid set constraint definition; it does not contain"*
                   "the state $state, the input $input or noise term $noise")
        end
    end
    throw(ArgumentError("the set entry $(expr) does not have the correct form `x_ ∈ X_`"))
end

"""
    sort(parameters::Vector, order::Tuple)

Filter and sort the vector `parameters` according to `order`.

### Input

- `parameters` -- vector of tuples
- `order`      -- tuple of symbols

### Output

A new vector of tuples corresponding to `parameters` filtered and sorted according
to `order`.

### Examples

```jldoctest
julia> parameters= [(:U1, :U), (:X1, :X), (:W1, :W)];

julia> sort(parameters, (:X, :U, :W))
3-element Array{Tuple{Any,Symbol},1}:
 (:X1, :X)
 (:U1, :U)
 (:W1, :W)

julia>  parameters = [(:const, :c), (:A, :A)];

julia> sort(parameters, (:A, :B, :c, :D))
2-element Array{Tuple{Any,Symbol},1}:
 (:A, :A)
 (:const, :c)
```

### Notes

`parameters` is a vector that contains tuples whose second element
is considered for the sorting according to `order`.

If a value of `order` is not contained in `parameters`, the corresponding entry of
`order` will be omitted.
"""
function sort(parameters::Vector{<:Tuple{Any, Symbol}}, order::NTuple{N, Symbol}) where {N}
    order_parameters = Vector{Tuple{Any, Symbol}}()
    for ordered_element in order
        for tuple in parameters
            if tuple[2] == ordered_element
                push!(order_parameters, tuple)
            end
        end
    end
    return order_parameters
end

function _get_system_type(dyn_eq, AT, constr, state, input, noise, dim)
    lhs, rhs = extract_dyn_equation_parameters(dyn_eq, state, input, noise, dim, AT)
    ordered_rhs = sort(rhs, (:A, :B, :c, :D, :f, :statedim, :inputdim, :noisedim))
    ordered_set = sort(extract_set_parameter.(constr, state, input, noise), (:X, :U, :W))
    field_names, var_names = constructor_input(lhs, ordered_rhs, ordered_set)
    sys_type = _corresponding_type(AT, field_names)
    return sys_type, var_names
end

"""
    system(expr...)

Return an instance of the system type corresponding to the given expressions.

### Input

- `expr` -- expressions separated by commas which define the dynamic equation,
            the constraint sets or the dimensionality of the system

### Output

A system that best matches the given expressions.

### Notes

**Terms.** The expression `expr` contains one or more of the following sub-expressions:

- dynamic equation, either continuous, e.g.`x' = Ax`, or discrete, e.g. `x⁺ = Ax`
- set constraints, e.g. `x ∈ X`
- input constraints, e.g. `u ∈ U`
- dimensionality, e.g. `dim: (2,1)` or `dim = 1`
- specification of the input variable, e.g. `input: u` or `input = u`
- specification of the noise variable, e,g, `noise: w` or `noise = w`

The macro call is then formed by separating the previous sub-expressions
(which we simply call *terms* hereafter), as in:

```julia
@system(dynamic eq., set constr., input constr., input specif., noise spec., dimens.)
```

The different terms that compose the system's definition do not have to appear in
any particular order. Moreover, the only mandatory term is the dynamic equation;
the other terms are optional and default values may apply depending on the system
type; this is explained next.

**Dynamic equation.** The time derivative in a continuous system is specified
by using `'`, as in `x' = A*x`. A discrete system is specified using `⁺` (which
can be written with the combination of keys `\\^+[TAB]`), as in `x⁺ = A*x`.
Moreover, the asterisk denoting matrix-vector products is optional. For instance,
both `x' = Ax` and `x' = A*x` are parsed as the linear continuous system whose
state matrix is `A`. The matrix is supposed to be defined at the call site.

**Default values.** When the dynamic equation is parsed, the variable on the
left-hand side is interpreted as the state variable. The input variable is by
default `u` and the noise variable is by default `w`. If we want to change the
default name of the input variable, this can be done by adding the term
`input: var` (or equivalently, `input=var)` where `var` corresponds to the new
name of the input variable, eg. `@system(x' = A*x + B*v, input:v)`.
Similarly, a noise variable is specified with `noise: var` or `noise=var`.

**Exceptions.** The following exceptions and particular cases apply:

- If the right-hand side has the form `A*x + B*foo`, `A*x + B*foo + c` or
  `f(x, foo)`, the equation is parsed as a controlled linear (affine) or controlled
  black-box system with input `foo`. Note that in this case, `input` variable does
  not correspond to the default value of `u`, but `foo` is parsed as being the input.

- If the left-hand side contains a multiplicative term in the form `E*x⁺` or `E*x'`,
  the equation is parsed as an algebraic system. In this case, the asterisk
  `*` operator is mandatory.

- Systems of the form `x' = α*x` where `α` is a scalar are parsed as linear
  systems. The default dimension is `1` and `α` is parsed as `Float64`;
  if the system is higher-dimensional, use `dim`, as in `x' = 2x, dim=3`.

### Examples

Let us first create a continuous linear system using this macro:

```jldoctest system_macro
julia> A = [1. 0; 0 1.];

julia> @system(x' = A*x)
LinearContinuousSystem{Float64,Array{Float64,2}}([1.0 0.0; 0.0 1.0])
```

A discrete system is defined by using `⁺`:

```jldoctest system_macro
julia> @system(x⁺ = A*x)
LinearDiscreteSystem{Float64,Array{Float64,2}}([1.0 0.0; 0.0 1.0])
```

Additionally, a set definition `x ∈ X` can be added to create a constrained system.
For example, a discrete controlled affine system with constrained states and
inputs writes as:

```jldoctest system_macro
julia> using LazySets

julia> B = Matrix([1. 0.5]');

julia> c = [1., 1.5];

julia> X = BallInf(zeros(2), 10.);

julia> U = BallInf(zeros(1), 2.);

julia> @system(x' = A*x + B*u + c, x ∈ X, u ∈ U)
ConstrainedAffineControlContinuousSystem{Float64,Array{Float64,2},Array{Float64,2},Array{Float64,1},BallInf{Float64,Array{Float64,1}},BallInf{Float64,Array{Float64,1}}}([1.0 0.0; 0.0 1.0], [1.0; 0.5], [1.0, 1.5], BallInf{Float64,Array{Float64,1}}([0.0, 0.0], 10.0), BallInf{Float64,Array{Float64,1}}([0.0], 2.0))
```
For the creation of a black-box system, the state, input and noise dimensions have
to be defined separately. For a constrained controlled black-box system, the macro
writes as

```jldoctest system_macro
julia> f(x, u) = x + u;

julia> @system(x⁺ = f(x, u), x ∈ X, u ∈ U, dim: (2,2))
ConstrainedBlackBoxControlDiscreteSystem{typeof(f),BallInf{Float64,Array{Float64,1}},BallInf{Float64,Array{Float64,1}}}(f, 2, 2, BallInf{Float64,Array{Float64,1}}([0.0, 0.0], 10.0), BallInf{Float64,Array{Float64,1}}([0.0], 2.0))
```
"""
macro system(expr...)
    try
        dyn_eq, AT, constr, state, input, noise, dim, x0 = _parse_system(expr)
        sys_type, var_names = _get_system_type(dyn_eq, AT, constr, state, input, noise, dim)
        sys = Expr(:call, :($sys_type), :($(var_names...)))
        if x0 == nothing
            return esc(sys)
        else
            ivp = Expr(:call, InitialValueProblem, :($sys), :($x0))
            return esc(ivp)
        end
    catch ex
        if  isa(ex, ArgumentError)
            return :(throw($ex))
        else
            throw(ex)
        end
    end
end

"""
    ivp(expr...)

Return an instance of the initial-value problem type corresponding to the given
expressions.

### Input

- `expr` -- expressions separated by commas which define the dynamic equation,
            the constraint sets or the dimensionality of the system, and the set
            of initial states (required)

### Output

An initial-value problem that best matches the given expressions.

### Notes

This macro behaves like the `@system` macro, the sole difference being that in
`@ivp` the constraint on the set of initial states is mandatory. For the technical
details we refer to the documentation of [`@system`](@ref).

The macro can also be called with a `system` argument of type `AbstractSystem`
in the form `@ivp(system, state(0) ∈ initial_set)`.

### Examples

```jldoctest ivp_macro
julia> p = @ivp(x' = -x, x(0) ∈ [1.0]);

julia> typeof(p)
InitialValueProblem{LinearContinuousSystem{Float64,IdentityMultiple{Float64}},Array{Float64,1}}

julia> initial_state(p)
1-element Array{Float64,1}:
 1.0

julia> sys = @system(x' = [1 0; 0 1] * x);

julia> @ivp(sys, x(0) ∈ [-1, 1])
InitialValueProblem{LinearContinuousSystem{Int64,Array{Int64,2}},Array{Int64,1}}(LinearContinuousSystem{Int64,Array{Int64,2}}([1 0; 0 1]), [-1, 1])
```
"""
macro ivp(expr...)
    try
        # check if the @ivp macro is called with an AbstractSystem argument
        if parses_as_system(expr[1])
            sys = expr[1]
            @capture(expr[2], x_(0) ∈ x0_) || throw(ArgumentError("malformed epxpression")) # TODO handle equality
            ivp = Expr(:call, InitialValueProblem, :($(expr[1])), :($x0))
            return esc(ivp)
        else
            dyn_eq, AT, constr, state, input, noise, dim, x0 = _parse_system(expr)
            sys_type, var_names = _get_system_type(dyn_eq, AT, constr, state, input, noise, dim)
            sys = Expr(:call, :($sys_type), :($(var_names...)))
            if x0 == nothing
                return throw(ArgumentError("an initial-value problem should define the " *
                            "initial states, but such expression was not found"))
            else
                ivp = Expr(:call, InitialValueProblem, :($sys), :($x0))
                return esc(ivp)
            end
        end
    catch ex
        if  isa(ex, ArgumentError)
            return :(throw($ex))
        else
            throw(ex)
        end
    end
end

function parses_as_system(expr)
    # The argument is a system if
    # 1) it is directly defined using @system in the macro call
    if @capture(expr, @system(def_))
        return true
    end
    # 2) it is a variable, i.e. it is not a dynamic equation
    if !@capture(expr, a_ = b_)
        return true
    end
    return false
end


text:
# Solves the Poisson equation using the Mixed finite element method
using Revise
using AdFem
using DelimitedFiles
using SparseArrays
using PyPlot

λ = 1.0
μ = 1.0
n = 50
mmesh = Mesh(n, n, 1/n, degree = BDM1)
a = 1/2μ
b = -λ/(2μ*(2μ+2λ))

TestCase = [
    (
        (x,y)->begin;x*y*(x - 1)*(y - 1);end,
        (x,y)->begin;x*y*(x - 1)*(y - 1);end,
        (x,y)->begin;2.0*x^2 + 8.0*x*y - 6.0*x + 6.0*y^2 - 10.0*y + 2.0;end,
        (x,y)->begin;6.0*x^2 + 8.0*x*y - 10.0*x + 2.0*y^2 - 6.0*y + 2.0;end,
    ),
    (
        (x,y)->begin;x^2*y^2*(x - 1)*(y^2 - 1);end,
        (x,y)->begin;x^2*y^2*(x - 1)*(y^2 - 1);end,
        (x,y)->begin;12.0*x^3*y^2 - 2.0*x^3 + 24.0*x^2*y^3 - 12.0*x^2*y^2 - 12.0*x^2*y + 2.0*x^2 + 18.0*x*y^4 - 16.0*x*y^3 - 18.0*x*y^2 + 8.0*x*y - 6.0*y^4 + 6.0*y^2;end,
        (x,y)->begin;36.0*x^3*y^2 - 6.0*x^3 + 24.0*x^2*y^3 - 36.0*x^2*y^2 - 12.0*x^2*y + 6.0*x^2 + 6.0*x*y^4 - 16.0*x*y^3 - 6.0*x*y^2 + 8.0*x*y - 2.0*y^4 + 2.0*y^2;end,
    ),
    (
        (x,y)->begin;x^2*y^2*(x - 1)*(y^2 - 1);end,
        (x,y)->begin;x*y*(x - 1)*(y - 1);end,
        (x,y)->begin;12.0*x^3*y^2 - 2.0*x^3 - 12.0*x^2*y^2 + 2.0*x^2 + 18.0*x*y^4 - 18.0*x*y^2 + 8.0*x*y - 4.0*x - 6.0*y^4 + 6.0*y^2 - 4.0*y + 2.0;end,
        (x,y)->begin;24.0*x^2*y^3 - 12.0*x^2*y + 6.0*x^2 - 16.0*x*y^3 + 8.0*x*y - 6.0*x + 2.0*y^2 - 2.0*y;end,
    )

]

for k = 1:length(TestCase)
    @info "Running TestCase $k..."
    ufunc, vfunc, gfunc, hfunc = TestCase[k]

    A = compute_fem_bdm_mass_matrix(a*ones(get_ngauss(mmesh)), b*ones(get_ngauss(mmesh)), mmesh)
    B = compute_fem_bdm_div_matrix(mmesh)
    C = compute_fem_bdm_skew_matrix(mmesh)

    D = [A B' C'
        B spzeros(2mmesh.nelem, 3mmesh.nelem)
        C spzeros(mmesh.nelem, 3mmesh.nelem)]

    gD = bcedge(mmesh)
    t1 = eval_f_on_gauss_pts(gfunc, mmesh)
    t2 = eval_f_on_gauss_pts(hfunc, mmesh)
    f1 = compute_fvm_source_term(t1, mmesh)
    f2 = compute_fvm_source_term(t2, mmesh)

    rhs = [zeros(2mmesh.ndof); f1; f2; zeros(mmesh.nelem)]

    sol = D\rhs
    u = sol[2mmesh.ndof+1:2mmesh.ndof+2mmesh.nelem]
    close("all")
    figure(figsize=(15, 10))
    subplot(231)
    title("Reference")
    xy = fvm_nodes(mmesh)
    x, y = xy[:,1], xy[:,2]
    uf = ufunc.(x, y)
    visualize_scalar_on_fvm_points(uf, mmesh)
    subplot(232)
    title("Numerical")
    visualize_scalar_on_fvm_points(u[1:mmesh.nelem], mmesh)
    subplot(233)
    title("Absolute Error")
    visualize_scalar_on_fvm_points( abs.(u[1:mmesh.nelem] - uf) , mmesh)

    subplot(234)
    xy = fvm_nodes(mmesh)
    x, y = xy[:,1], xy[:,2]
    uf = vfunc.(x, y)
    visualize_scalar_on_fvm_points(uf, mmesh)
    subplot(235)
    visualize_scalar_on_fvm_points(u[mmesh.nelem+1:end], mmesh)
    subplot(236)
    visualize_scalar_on_fvm_points( abs.(u[mmesh.nelem+1:end] - uf) , mmesh)
    savefig("varying_elasticity$k.png")
end


text:
function solve_lorenz(u0, tspan, p, time_interval=0.01)
    function lorenz(du, u, p, t)
            du[1] = p[1]*(u[2]-u[1])
            du[2] = u[1]*(p[2] - u[3]) - u[2]
            du[3] = u[1]*u[2] - p[3]*u[3]
    end
    problem = ODEProblem(lorenz, u0, tspan, p)
    sol = solve(problem, saveat=time_interval)
    return sol
end


text:
# implements lazy matrix products and sums
# in contrast to LazyArrays, it implements a different type of laziness
# LazyArrays is fully lazy and avoids computing temporaries, but assumes the
# entries of the component arrays can be computed easily and efficiently
# The structures below rely on the mul! functions of the constituent
# matrices / factorizations to define a potentially more efficient multiplication\
# IDEA: introduce temporaries and optimal matrix multiplication order to minimize complexity
# IDEA: enable rdiv! for LazyMatrixSum, Product
abstract type LazyFactorization{T} <: Factorization{T} end
# const default_tol = 1e-6

################################################################################
# in contrast to the AppliedMatrix in LazyArrays, this is not completely lazy
# in that it calculates intermediate results, important to make use of special structure
# observation: this is remarkably similar to the KroneckerProducts implementation
# allow AbstractMatOrFacOrUni? will need adjustment for size function
struct LazyMatrixProduct{T, AT, F} <: LazyFactorization{T}
    args::AT
    tol::F
    # temporaries
    # args is vector or tuple of constituent matrices
    function LazyMatrixProduct(args; tol::Real = default_tol)
        for i in 1:length(args)-1
            size(args[i], 2) == size(args[i+1], 1) || throw(DimensionMismatch("$i"))
        end
        T = promote_type(eltype.(args)...)
        new{T, typeof(args), typeof(tol)}(args, tol)
    end
end
function LazyMatrixProduct(A::AbstractMatOrFac...; tol::Real = default_tol)
    args_vec = vcat((lazy_product_flattened_args(Ai) for Ai in A)...)
    args = tuple(args_vec...) # casting to tuple makes it possible to infer types
    LazyMatrixProduct(args, tol = tol)
end
lazy_product_flattened_args(A) = [A]
lazy_product_flattened_args(A::LazyMatrixProduct) = [A.args...]

function Base.size(L::LazyMatrixProduct, i::Int)
    if i == 1
        size(L.args[1], 1)
    elseif i == 2
        size(L.args[end], 2)
    elseif i > 2
        1
    end
end
Base.size(L::LazyMatrixProduct) = size(L, 1), size(L, 2)
Base.adjoint(L::LazyMatrixProduct) = LazyMatrixProduct(reverse(adjoint.(L.args)))
issquare(A::AbstractMatOrFac) = size(A, 1) == size(A, 2)
Base.Matrix(L::LazyMatrixProduct) = prod(Matrix, L.args)
Base.AbstractMatrix(L::LazyMatrixProduct) = prod(AbstractMatrix, L.args)

# function LinearAlgebra.diag(L::LazyMatrixProduct)
#     middle = length(L.args) > 2 ? *(L.args[2:end-1]) : 1
#     n = minimum(size(L))
#     d = zeros(eltype(L), n)
#     for i in 1:n
#         d = @views L.args[1][i, :]' * middle * L.args[end][:, i]
#     end
#     return d
# end
# function LinearAlgebra.factorize(A::LazyMatrixProduct; k::Int = 16, sigma::Real = 1e-2)
#     return CGFact(A, tol = A.tol)
# end

Base.:*(x::AbstractMatrix, L::LazyMatrixProduct) = adjoint(L'*x')
function Base.:*(L::LazyMatrixProduct, x::AbstractVector) # IDEA: could consolidate with LazyMatrixSum
    T = promote_type(eltype(L), eltype(x))
    b = zeros(T, size(L, 1))
    mul!(b, L, x)
end
function Base.:*(L::LazyMatrixProduct, X::AbstractMatrix)
    T = promote_type(eltype(L), eltype(X))
    B = zeros(T, (size(L, 1), size(X, 2)))
    mul!(B, L, X)
end
# IDEA: this could be more memory-efficient if we knew that all matrices had the same size
function LinearAlgebra.mul!(y::AbstractVecOrMat, L::LazyMatrixProduct, x::AbstractVecOrMat, α::Real = 1, β::Real = 0)
    z = deepcopy(x)
    for A in reverse(L.args)
        z = A*z
    end
    @. y = α*z + (iszero(β) ? 0 : β*y)
    return y
end

################################################################################
# in contrast to the AppliedMatrix in LazyArrays, this is not completely lazy
# in that it calculates intermediate results
# TODO: should we have a tol field?
struct LazyMatrixSum{T, AT, F} <: LazyFactorization{T}
    args::AT
    tol::F
    function LazyMatrixSum(args; tol::Real = default_tol)
        all(==(size(args[1])), size.(args)) || throw(DimensionMismatch())
        T = promote_type(eltype.(args)...)
        new{T, typeof(args), typeof(tol)}(args, tol)
    end
end
function LazyMatrixSum(A::AbstractMatOrFac...; tol::Real = default_tol)
    args_vec = vcat((lazy_sum_flattened_args(Ai) for Ai in A)...)
    args = tuple(args_vec...)
    LazyMatrixSum(args, tol = tol)
end
lazy_sum_flattened_args(A) = [A]
lazy_sum_flattened_args(A::LazyMatrixSum) = [A.args...]

Base.size(L::LazyMatrixSum, i...) = size(L.args[1], i...)
Base.adjoint(L::LazyMatrixSum) = LazyMatrixSum(adjoint.(L.args))
Base.Matrix(L::LazyMatrixSum) = sum(Matrix, L.args)
Base.AbstractMatrix(L::LazyMatrixSum) = sum(AbstractMatrix, L.args)
Base.getindex(L::LazyMatrixSum, i...) = sum(x->getindex(x, i...), L.args)
LinearAlgebra.diag(L::LazyMatrixSum) = sum(diag, L.args)

Base.:*(x::AbstractMatrix, L::LazyMatrixSum) = adjoint(L'*x')
function Base.:*(L::LazyMatrixSum, x::AbstractVector)
    T = promote_type(eltype(L), eltype(x))
    b = zeros(T, size(L, 1))
    mul!(b, L, x)
end
function Base.:*(L::LazyMatrixSum, X::AbstractMatrix)
    T = promote_type(eltype(L), eltype(X))
    B = zeros(T, (size(L, 1), size(X, 2)))
    mul!(B, L, X)
end
function LinearAlgebra.mul!(y::AbstractVecOrMat, L::LazyMatrixSum, x::AbstractVecOrMat,
                            α::Real = 1, β::Real = 0)
    @. y = (iszero(β) ? 0 : β*y)
    for A in L.args
        mul!(y, A, x, α, 1)
    end
    return y
end

function LinearAlgebra.:\(A::LazyFactorization, b::AbstractVector)
    T = promote_type(eltype(A), eltype(b))
    x = zeros(T, size(A, 1))
    ldiv!(x, A, b)
end
# solve general BlockGramian via minimum residual solver
function LinearAlgebra.ldiv!(x::AbstractVector, A::LazyFactorization, b::AbstractVector; kwargs...)
    minres!(x, A, b; kwargs...)
end

# IDEA: more efficient temporary allocation
# function LinearAlgebra.mul!(y::AbstractVector, L::LazyMatrixProduct, x::AbstractVector, α::Real = 1, β::Real = 0)
#     if all(issquare, L.args)
#         z1, z2 = similar(x), similar(x)
#
#         # if odd and greater than 1, we need to start with z
#         iseven(length(L.args)) ? copyto!(y, x) : copyto!(z, x)
#         for A in L.args
#             mul!(z, A, y)
#             y, z = z, y
#         end
#
#     else
#         z = copy(x)
#         for A in L.args
#             z = A*z
#         end
#         copyto!(y, z)
#     end
#     return y
# end


text:
using Base.Test
using Sodium

function parse_hexstring(s::String)
    out = zeros(Uint8, int(length(s)/2))
    for i in 1:length(s)/2
        offset = i * 2
        out[i] = parseint(s[offset-1:offset], 16)
    end
    out
end

hexstring(arr::Array{Uint8,1}) = join([@sprintf("%02.x", i) for i in arr])

@test 0 == sodium_init()

println("Testing with sodium library version : " * bytestring(sodium_version_string()))

tstr = "The quick brown fox jumps over the lazy dog"

print("Testing crypto_hash_sha256 : ")

out = zeros(Uint8, 32)

tstr_sha256 = "d7a8fbb307d7809469ca9abcb0082e4f8d5651e46d3cdb762d02d0bf37c9e592"
refval = parse_hexstring(tstr_sha256)

@test 0 == crypto_hash_sha256(out, tstr, length(tstr))
@test refval == out
println("PASSED")

print("Testing crypto_hash_sha512 : ")
tstr_sha512 = "07e547d9586f6a73f73fbac0435ed76951218fb7d0c8d788a309d785436bbb642e93a252a954f23912547d1e8a3b5ed6e1bfd7097821233fa0538f3db854fee6"
out = zeros(Uint8, 64)
@test 0 == crypto_hash_sha512(out, tstr, length(tstr))
refval = parse_hexstring(tstr_sha512)
@test refval == out
println("PASSED")

print("Testing crypto_auth_hmacsha256 : ")
tstr_hmacsha256 = "f7bc83f430538424b13298e6aa6fb143ef4d59a14946175997479dbc2d1a3cd8"
out = zeros(Uint8, 32)
key = zeros(Uint8, 32)
key[1] = 'k'
key[2] = 'e'
key[3] = 'y'

@test 0 == crypto_auth_hmacsha256(out, tstr, length(tstr), key)
refval = parse_hexstring(tstr_hmacsha256)
@test refval == out
println("PASSED")


text:
using Documenter, StochasticIntegrals

makedocs(
    format = Documenter.HTML(),
    sitename = "StochasticIntegrals",
    modules = [StochasticIntegrals],
    pages = ["index.md"]
)

deploydocs(
    repo   = "github.com/s-baumann/StochasticIntegrals.jl.git",
    target = "build",
    deps   = nothing,
    make   = nothing
)


text:
"""
$(SIGNATURES)

Average number of sigma bound exceedings

Returns the average number of unbiased states / innovation values that exceed the ⨦σ bound of the given covariance
"""
function mean_num_sigma_bound_exceedings(state_over_time::Vector{Vector{T}}, covariance_over_time::Vector{Matrix{T}}) where T
    mean(map((x, P) -> abs.(x) .> sqrt.(diag(P)), state_over_time, covariance_over_time))
end

function mean_num_sigma_bound_exceedings(state_over_time::Vector{T}, covariance_over_time::Vector{T}) where T
    mean(map((x, P) ->  abs(x) > sqrt(P), state_over_time, covariance_over_time))
end

"""
$(SIGNATURES)

Innovation magnitude bound test (σ-bound-test)

Tests if approximately 68% of state values lie within the ⨦σ bound
"""
function sigma_bound_test(state_over_time, covariance_over_time)
    isapprox.(mean_num_sigma_bound_exceedings(state_over_time, covariance_over_time), .32, atol = .015)
end

"""
$(SIGNATURES)

Innovation magnitude bound test (2σ-bound-test)

Tests if approximately 95% of state values lie within the ⨦2σ bound
"""
function two_sigma_bound_test(state_over_time, covariance_over_time)
    isapprox.(mean_num_sigma_bound_exceedings(state_over_time, 4 .* covariance_over_time), .05, atol = .008)
end

"""
$(SIGNATURES)

Normalized innovation squared (NIS) Test

Double-tailed siginicance test with false alarm probability α = 0.05

Calculates confidence interval [r1 r2] and tests Prob{ ∑ NIS values)} ∈ [r1 r2] ∣ H_0 ) = 1 - α
with Hypothesis H_0: N * ∑ NIS values ∼ χ^2_{dof}
     dof (degree of freedom): N * m (N: window length, m: dimension of state vector)
"""
function nis_test(nis_over_time, dof)
    sum_of_nis = sum(nis_over_time)

    r1 = cquantile(Chisq(dof), .975)
    r2 = cquantile(Chisq(dof), .025)

    (sum_of_nis >= r1) && (sum_of_nis <= r2)
end

"""
$(SIGNATURES)

Normalized innovation squared (NIS)

Returns NIS-value for a single innovation sequence `seq` and its variance `var`
"""
function calc_nis(seq, var)
    dot(seq, var \ seq)
end

function calc_nis(mu::AbstractMeasurementUpdate)
    calc_nis(get_innovation(mu), mu.innovation_covariance)
end


text:
type ModelName

end

function ModelName()

end

function train!(model::ModelName, X::Matrix, y::Vector)

end

function predict(model::ModelName,
                 x::Matrix)
    n = size(x,1)
    res = zeros(n)
    for i = 1:n
        res[i] = predict(model, x[i,:])
    end
    return res
end

function predict(model::ModelName,
                 x::Vector)
end

function test_ModelName()
    X_train, X_test, y_train, y_test = make_cla()
    model = ModelName()
    train!(model,X_train, y_train)
    predictions = predict(model,X_test)
    print("classification accuracy", accuracy(y_test, predictions))
end


text:
# Reactor Separator Model
id = :PerezGalvan2017f

pL = [0.0; 0.0]
pU = [0.08; 0.16]
PerezGalvan2017f_x0(p::Vector{T}) where T = T[0.5; 0.0; 0.0; 0.0; p[1]; p[2]]
function PerezGalvan2017f_f(du::Vector{T}, u::Vector{T}, p::Vector{T}, t) where T
    k = 0.06; xF0 = 0.0; D = 1.0; H = 63.33
    α = 7.5; B = 1.2; L = 1.704; F = 1.0
    V = pi*(H/4)*D^2    # THIS MAY BE INCORRECT...
    xF = (F*xF0 + B*u[2])/(F + B)
    y2 = α*u[2]/(1.0 + (α-1.0)*u[2])
    y3 = α*u[3]/(1.0 + (α-1.0)*u[3])
    y4 = α*u[4]/(1.0 + (α-1.0)*u[4])
    y5 = α*u[5]/(1.0 + (α-1.0)*u[5])
    du[1] = (F+B)*(xF - u[1])/H + k*u[1]*(1.0 - u[1])
    du[2] = (L + F + B)*u[3] - B*u[2] - V*y2
    du[3] = (L + F + B)*(u[4] - u[3]) + V*(y2 - y3)
    du[4] = (F+B)*u[1] + L*u[5] - (L + F + B)*u[4] +V*(y3 - y4)
    du[5] = L*(u[6] - u[5]) + V*(y4 - y5)
    du[6] = -(L + D)*u[6] + V*y5
    return
end
tspan = (0.0, 100.0)

prob = ODERelaxProb(PerezGalvan2017d_f!, tspan, PerezGalvan2017d_x0, pL, pU)
url = "https://www.sciencedirect.com/science/article/abs/pii/S0098135417300923"
source = "PerezGalvan2017"
desc = "Reactor Separator Model"


text:
@testset "1006.clumsy-factorial.jl" begin
    @test clumsy(4) == 7
    @test clumsy(10) == 12
end


text:
avxt() = false

macro threaded(code)
    return esc(:(
        if $(@__MODULE__).avxt()
            @avxt($code)
        else
            @threads($code)
        end
    ))
end

# multithreaded kernels
function copy_allthreads(C, A; kwargs...)
    @assert length(C) == length(A)
    @threaded for i in eachindex(C, A)
        @inbounds C[i] = A[i]
    end
    return nothing
end

function scale_allthreads(B, C, s; kwargs...)
    @assert length(C) == length(B)
    @threaded for i in eachindex(C)
        @inbounds B[i] = s * C[i]
    end
    return nothing
end

function add_allthreads(C, A, B; kwargs...)
    @assert length(C) == length(B) == length(A)
    @threaded for i in eachindex(C)
        @inbounds C[i] = A[i] + B[i]
    end
    return nothing
end

function triad_allthreads(A, B, C, s; kwargs...)
    @assert length(C) == length(B) == length(A)
    @threaded for i in eachindex(C)
        @inbounds A[i] = B[i] + s * C[i]
    end
    return nothing
end


text:
module VortexHelperSwitchBlock
using ..Ahorn, Maple

const colorNames = Dict{String, Int}(
    "Blue" => 0,
    "Rose" => 1,
    "Orange" => 2,
    "Lime" => 3
)

@mapdef Entity "VortexHelper/SwitchBlock" SwitchBlock(
			x::Integer, y::Integer,
			width::Integer=8, height::Integer=8,
			index::Integer=0)

const placements = Ahorn.PlacementDict(
    "Switch Block ($index - $color) (Vortex Helper)" => Ahorn.EntityPlacement(
        SwitchBlock,
        "rectangle",
        Dict{String, Any}(
            "index" => index
        )
    ) for (color, index) in colorNames
)

Ahorn.editingOptions(entity::SwitchBlock) = Dict{String, Any}(
    "index" => colorNames
)

Ahorn.minimumSize(entity::SwitchBlock) = 16, 16
Ahorn.resizable(entity::SwitchBlock) = true, true

Ahorn.selection(entity::SwitchBlock) = Ahorn.getEntityRectangle(entity)

const colors = Dict{Int, Ahorn.colorTupleType}(
    1 => (255, 50, 101, 255) ./ 255,
	2 => (255, 149, 50, 255) ./ 255,
	3 => (156, 255, 50, 255) ./ 255
)
const defaultColor = (50, 50, 255, 255) ./ 255
const borderMultiplier = (0.9, 0.9, 0.9, 1)

const frame = "objects/VortexHelper/onoff/solid"

function getSwitchBlockRectangles(room::Maple.Room)
    entities = filter(e -> e.name == "VortexHelper/SwitchBlock", room.entities)
    rects = Dict{Int, Array{Ahorn.Rectangle, 1}}()

    for e in entities
        index = get(e.data, "index", 0)
        rectList = get!(rects, index) do
            Ahorn.Rectangle[]
        end

        push!(rectList, Ahorn.Rectangle(
            Int(get(e.data, "x", 0)),
            Int(get(e.data, "y", 0)),
            Int(get(e.data, "width", 8)),
            Int(get(e.data, "height", 8))
        ))
    end

    return rects
end

function notAdjacent(entity::SwitchBlock, ox, oy, rects)
    x, y = Ahorn.position(entity)
    rect = Ahorn.Rectangle(x + ox + 4, y + oy + 4, 1, 1)

    for r in rects
        if Ahorn.checkCollision(r, rect)
            return false
        end
    end

    return true
end

function drawSwitchBlock(ctx::Ahorn.Cairo.CairoContext, entity::SwitchBlock, room::Maple.Room)
    switchBlockRectangles = getSwitchBlockRectangles(room)

    x, y = Ahorn.position(entity)

    width = Int(get(entity.data, "width", 32))
    height = Int(get(entity.data, "height", 32))

    tileWidth = ceil(Int, width / 8)
    tileHeight = ceil(Int, height / 8)

    index = Int(get(entity.data, "index", 0))
    color = get(colors, index, defaultColor)

    rect = Ahorn.Rectangle(x, y, width, height)
    rects = get(switchBlockRectangles, index, Ahorn.Rectangle[])

    if !(rect in rects)
        push!(rects, rect)
    end

    for x in 1:tileWidth, y in 1:tileHeight
        drawX, drawY = (x - 1) * 8, (y - 1) * 8

        closedLeft = !notAdjacent(entity, drawX - 8, drawY, rects)
        closedRight = !notAdjacent(entity, drawX + 8, drawY, rects)
        closedUp = !notAdjacent(entity, drawX, drawY - 8, rects)
        closedDown = !notAdjacent(entity, drawX, drawY + 8, rects)
        completelyClosed = closedLeft && closedRight && closedUp && closedDown

        if completelyClosed
            if notAdjacent(entity, drawX + 8, drawY - 8, rects)
                Ahorn.drawImage(ctx, frame, drawX, drawY, 24, 0, 8, 8, tint=color)

            elseif notAdjacent(entity, drawX - 8, drawY - 8, rects)
                Ahorn.drawImage(ctx, frame, drawX, drawY, 24, 8, 8, 8, tint=color)

            elseif notAdjacent(entity, drawX + 8, drawY + 8, rects)
                Ahorn.drawImage(ctx, frame, drawX, drawY, 24, 16, 8, 8, tint=color)

            elseif notAdjacent(entity, drawX - 8, drawY + 8, rects)
                Ahorn.drawImage(ctx, frame, drawX, drawY, 24, 24, 8, 8, tint=color)

            else
                Ahorn.drawImage(ctx, frame, drawX, drawY, 8, 8, 8, 8, tint=color)
            end

        else
            if closedLeft && closedRight && !closedUp && closedDown
                Ahorn.drawImage(ctx, frame, drawX, drawY, 8, 0, 8, 8, tint=color)

            elseif closedLeft && closedRight && closedUp && !closedDown
                Ahorn.drawImage(ctx, frame, drawX, drawY, 8, 16, 8, 8, tint=color)

            elseif closedLeft && !closedRight && closedUp && closedDown
                Ahorn.drawImage(ctx, frame, drawX, drawY, 16, 8, 8, 8, tint=color)

            elseif !closedLeft && closedRight && closedUp && closedDown
                Ahorn.drawImage(ctx, frame, drawX, drawY, 0, 8, 8, 8, tint=color)

            elseif closedLeft && !closedRight && !closedUp && closedDown
                Ahorn.drawImage(ctx, frame, drawX, drawY, 16, 0, 8, 8, tint=color)

            elseif !closedLeft && closedRight && !closedUp && closedDown
                Ahorn.drawImage(ctx, frame, drawX, drawY, 0, 0, 8, 8, tint=color)

            elseif !closedLeft && closedRight && closedUp && !closedDown
                Ahorn.drawImage(ctx, frame, drawX, drawY, 0, 16, 8, 8, tint=color)

            elseif closedLeft && !closedRight && closedUp && !closedDown
                Ahorn.drawImage(ctx, frame, drawX, drawY, 16, 16, 8, 8, tint=color)
            end
        end
    end
end

Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::SwitchBlock, room::Maple.Room) = drawSwitchBlock(ctx, entity, room)

end

max_stars_repo_name: catapillie/VortexHelper

max_issues_repo_name: catapillie/VortexHelper

max_forks_repo_name: catapillie/VortexHelper


text:
#
# Copyright (c) 2021 Tobias Thummerer, Lars Mikelsons, Josef Kircher
# Licensed under the MIT license. See LICENSE file in the project root for details.
#

using FMI

# our simulation setup
t_start = 0.0
t_stop = 8.0

# this FMU runs under Windows/Linux
pathToFMU = joinpath(dirname(@__FILE__), "../model/OpenModelica/v1.17.0/SpringFrictionPendulum1D.fmu")

# this FMU runs only under Windows
if Sys.iswindows()
    pathToFMU = joinpath(dirname(@__FILE__), "../model/Dymola/2020x/SpringFrictionPendulum1D.fmu")
end

# load the FMU container
myFMU = fmiLoad(pathToFMU)

# print some useful FMU-information into the REPL
fmiInfo(myFMU)

# make an instance from the FMU
fmiInstantiate!(myFMU; loggingOn=true)

# setup the experiment, start time = 0.0 (optional for setup=true)
#fmiSetupExperiment(myFMU, t_start)

# enter and exit initialization (optional for setup=true)
#fmiEnterInitializationMode(myFMU)
#fmiExitInitializationMode(myFMU)

# run the FMU in mode Model-Exchange (ME) with adaptive step sizes, result values are stored in `solution`
solution, _ = fmiSimulateME(myFMU, t_start, t_stop; setup=true)

# plot the results
fmiPlot(myFMU, solution)

# unload the FMU, remove unpacked data on disc ("clean up")
fmiUnload(myFMU)


text:
module Stack

s = zeros(Int64, 30)
depth = 0

function reset()
  global depth = 0
  global s = zeros(s)
end

function push(addr)
  global depth

  if depth + 1 > 30
    throw( OverflowError() )
  end

  depth = depth + 1

  s[depth] = addr
end

function pop()
  global depth

  if depth - 1 < 0
    throw( OverflowError() )
  end

  addr = s[depth]
  depth = depth - 1

  return addr
end

end

max_stars_repo_name: kinoroy/KCPSM-Simulator

max_issues_repo_name: kinoroy/KCPSM-Simulator

max_forks_repo_name: kinoroy/KCPSM-Simulator


text:
const PPLC = Union{Path{Port},Path{Link},Path{Component}}

"""
Representation of the routing resources of a `TopLevel`.
"""
struct RoutingGraph
    """
    Adjacency information of routing resources, encode as a
    `LightGraphs.SimpleDiGraph`.
    """
    graph::SimpleDiGraph{Int64}

    """
    Translation information mapping elements on the parent `TopLevel` to indices
    in `graph`.

    Implemented as a `Dict{Path{<:Union{Port,Link,Component}}, Int64}` where the
    values in the dict are the vertex index in `graph` of the key.
    """
    map::Dict{PPLC,Int64}

    function RoutingGraph(graph::AbstractGraph, map::Dict{P,Int64}) where {P<:Path}

        #= Make sure that all of the values in the port map are valid. =#
        num_vertices = nv(graph)
        for v in values(map)
            if v > num_vertices
                error("Map value $v greater than number of vertices: $(nv(graph))")
            end
        end
        m = Dict{PPLC,Int}(map)
        return new(graph, m)
    end
end

"""
Return the `map` of a [`RoutingGraph`](@ref)

Method List
-----------
$(METHODLIST)
"""
getmap(graph::RoutingGraph) = graph.map

function translate_routes(r::RoutingGraph, graphs::Vector{<:SparseDiGraph})
    map_reversed = rev_dict(r.map)

    routes = map(graphs) do g
        route = SparseDiGraph{PPLC}()
        # Add vertices
        for v in vertices(g)
            path = map_reversed[v]
            add_vertex!(route, path)
        end
        # Add edges
        for src in vertices(g), dst in outneighbors(g, src)
            add_edge!(route, map_reversed[src], map_reversed[dst])
        end

        return route
    end

    return routes
end

# Make a vertex for each port and a single vertex in the middle.
# Connect all inputs with an edge input -> middle
#
# Connect all outputs with an edge middle -> output
function build_routing_mux(c::Component)
    # Assert that this component has not children or intra-component routing
    assert_no_children(c)
    assert_no_intrarouting(c)

    g = DiGraph(length(c.ports) + 1)
    # Create a path for node #1, the node in the middle of the mux
    m = Dict{Path,Int}(Path{Component}() => 1)

    for (i, (p, v)) in enumerate(c.ports)
        # Add 1 to i to account for the Path{Component}() we already added.
        node_index = i + 1
        m[Path{Port}(p)] = node_index

        if v.class == Output
            add_edge!(g, 1, node_index)
        elseif v.class == Input
            add_edge!(g, node_index, 1)
        else
            throw(KeyError(v.class))
        end
    end

    return RoutingGraph(g, m)
end

# Create a graph with vertices for each port and no edges.
function build_routing_blackbox(c::Component)
    # Get all of the ports for this component.
    ports = keys(c.ports)
    # Instantiate the number of vertices equal to the number of ports.
    g = DiGraph(length(ports))
    # Sequentially assign port mappings.
    if length(ports) == 0
        m = Dict{Path,Int}()
    else
        m = Dict{Path,Int}(Path{Port}(p) => i for (i, p) in enumerate(ports))
    end
    return RoutingGraph(g, m)
end

const __ROUTING_DICT = Dict("mux" => build_routing_mux)

# Constructor dispatch based on primitive type. Default to BlackBox
function routing_skeleton(c::Component)::RoutingGraph
    f = get(__ROUTING_DICT, c.primitive, build_routing_blackbox)
    return f(c)
end

function routing_skeleton(tl::TopLevel)::RoutingGraph
    # Return an empty graph
    g = DiGraph(0)
    m = Dict{Path,Int}()
    return RoutingGraph(g, m)
end

function record!(new_dict::Dict, old_dict::Dict, offset::Integer, extension)
    for (path, index) in old_dict
        # Update path
        new_path = catpath(extension, path)
        # Update index
        new_index = index + offset
        # Record the new item.
        haskey(new_dict, new_path) && error()
        new_dict[new_path] = new_index
    end
end

function add_subgraphs!(top::RoutingGraph, prefixes, subgraphs::Vector{<:RoutingGraph})
    # Iterate through the prefixes and subgraphs.
    for (prefix, subgraph) in zip(prefixes, subgraphs)
        # Get an offset calculated by the number of vertices in the current graph.
        offset = nv(top.graph)
        # Add the number of vertices in the subgraph to g.
        add_vertices!(top.graph, nv(subgraph.graph))
        # Create dictionary records for each of the new nodes.
        record!(top.map, subgraph.map, offset, prefix)
        # Record all the edges from the subgraph - account for offset.
        for edge in edges(subgraph.graph)
            add_edge!(top.graph, src(edge) + offset, dst(edge) + offset)
        end
    end
    return nothing
end

function add_links!(top, c::AbstractComponent)
    # Iterate through all links for the component.
    for (key, link) in c.links
        # Create a node for the link
        add_vertex!(top.graph)
        link_index = nv(top.graph)
        # Create an entry in the link table
        linkpath = Path{Link}(key)
        top.map[linkpath] = link_index
        # Connect the link to all its attached ports by looking through the
        # PortPaths in the link and getting the reference vertex in the graph
        # from the portmap
        for port in sources(link)
            src = top.map[port]
            add_edge!(top.graph, src, link_index)
        end
        for port in dests(link)
            dst = top.map[port]
            add_edge!(top.graph, link_index, dst)
        end
    end
    return nothing
end

function splicegraphs(component::AbstractComponent, topgraph::RoutingGraph, subgraphs)
    add_subgraphs!(topgraph, keys(component.children), subgraphs)
    add_links!(topgraph, component)
    return topgraph
end

function routing_graph(
    component::AbstractComponent, memoize = true, memo = Dict{String,RoutingGraph}()
)

    # Check if the graph for this component has been memoized and return the
    # memoized result if it has.
    if memoize && haskey(memo, component.name)
        return memo[component.name]
    end

    # If this component has not children, no need to route each of the children,
    # just build a routing skeleton for it.
    if length(component.children) == 0
        graph = routing_skeleton(component)

        # Recursively build graphs for each of the children of this component.
        # Once all the graphs for the children have been created, splice them into
        # a skeleton graph for this component.
    else
        subgraphs = [routing_graph(i, memoize, memo) for i in children(component)]
        topgraph = routing_skeleton(component)
        graph = splicegraphs(component, topgraph, subgraphs)
    end
    # Memoize this result
    memoize && memoize!(component, memo, component.name, graph)
    safety_check(component, graph)
    return graph
end

memoize!(::Component, md::Dict, key, value) = (md[key] = value)
memoize!(::TopLevel, md::Dict, key, value) = nothing

safety_check(::Component, g) = nothing
function safety_check(::TopLevel, g)
    maprev = rev_dict_safe(g.map)
    for (k, v) in maprev
        length(v) > 1 && throw(ErrorException("$k => $v"))
    end
    @debug routing_graph_info(g)
end

function routing_graph_info(g)
    return """
           Number of Vertices: $(nv(g.graph))
           Number of Edges: $(ne(g.graph))
           Entries in Map Dictionary: $(length(g.map))
           """
end


text:
using LinearAlgebra
using StaticArrays
using FastExpm
"""
Contains code to perform time evolution
"""

abstract type IntegratorType end

Base.@kwdef struct Piecewise{TS,EXPM} <: IntegratorType
    n_slices::TS = 1
    expm_method::EXPM
end

Base.@kwdef struct Continuous{OPTS} <: IntegratorType
    ode_options::OPTS
end

_get_integration_func() = pw_evolve

"""
Given a set of Hamiltonians (drift and control) compute the propagator for use with static arrays
"""
function pw_evolve(A, B, x, n_pulses, dt, n_slices, U0)
    U = U0
    for i = 1:n_slices
        # compute the propagator
        Htot = A
        for j = 1:n_pulses
            @views Htot = Htot + B[j] .* x[j, i]
        end
        U = exp((-1.0im * dt) .* Htot) * U
    end
    U
end

"""
Given a set of Hamiltonians compute the piecewise evolution saving the propagator for each time slice
"""
#TODO maybe rewrite this so its not allocation out inside but we allocate it ourselves, make it a vector of static arrays?
function pw_evolve_save(A::T, B, x, n_pulses, dt, n_slices) where {T<:StaticMatrix}
    out = Vector{T}(undef, n_slices)
    for i = 1:n_slices
        # compute the propagator
        Htot = A
        for j = 1:n_pulses
            @views Htot = Htot + B[j] * x[j, i]
        end
        @views U = exp(-1.0im * dt * Htot)# * U
        out[i] = U
    end
    out
end

"""
An almost allocation free (1 alloc in my tests) version of the Hamiltonian saver, since the update is inplace we have to be careful with the definition of out!
"""
function pw_ham_save!(A, B, x, n_pulses, n_slices, out)
    K = n_pulses
    Htot = zeros(eltype(A),size(A))

    @views @inbounds for i = 1:n_slices
        Htot .= 0.0 .* Htot
        @inbounds for j = 1:K
            @. Htot = Htot + B[j] * x[j, i]
        end
        @. out[i] = Htot + A
    end
end

"""
An almost allocation free (1 alloc in my tests) version of the Generator saving function, this saves the generators of propagators (implementation might differ for prob)
"""
function pw_gen_save!(A, B, x, n_pulses, n_slices, duration, out)
    K = n_pulses
    Htot = zeros(eltype(A),size(A))
    dt = duration / n_slices

    @views @inbounds for i = 1:n_slices
        Htot .= 0.0 .* Htot
        @inbounds for j = 1:K
            @. Htot = Htot + B[j] * x[j, i]
        end
        @. out[i] = (Htot + A) .* (-1.0im * dt)
    end
end

"""
Compute the propagator using fastExpm and save it in the array defined in out
"""
function pw_prop_save!(A, B, x, n_pulses, n_slices, dt, out)
    K = n_pulses
    # Htot = copy(A)
    Htot = zeros(eltype(A),size(A))

    @inbounds for i = 1:n_slices
        Htot .= 0.0 .* Htot
        @inbounds for j = 1:K
            Htot .= Htot .+ B[j] .* x[j, i]
        end
        out[i] .= exp((-1.0im * dt) .* (Htot .+ A))
    end
end


text:
@testset "GoogLeNet" begin
    m = googlenet()
    @test size(m(rand(Float32, 224, 224, 3, 50))) == (1000, 50)
end


text:
using RegularExpressions
using Documenter: deploydocs, makedocs

makedocs(
    sitename = "RegularExpressions.jl",
    modules = [RegularExpressions],
    doctest = false,
)
deploydocs(repo = "github.com/bramtayl/RegularExpressions.jl.git")


text:
export  subpopulationAgeReport

"""
```
subpopulationAgeReport(
    mpSim::MPsim,
    timeGrid::Vector{Float64},
    ageRes::Real,
    ageType::Symbol,
    subpopulations::Subpopulation... )
```
This function generates reports for the distribution of the ages of the personnel members in the valid subpopulations in `subpopulations` of the manpower simulation `mpSim`. The reports are generated on the grid of time points `timeGrid` and on a grid of ages with resoluion `ageRes`.

The type of age that is reported is determined by `ageType`, which can take these valid values:
* `:age` for the actual age of the personnel members,
* `:tenure` for the tenure of the personnel members, and
* `:timeInNode` for the time the personnel members are in their current base node. If the subpopulation is based on the entire population, this reports the tenure instead.

This function will issue a warning and not generate any report in the following two cases:
1. There are no positive time points in the time grid;
2. None of the entered subpopulations are valid subpopulations in the simulation, where a subpopulation is invalid if its souce nodes doesn't exist in the simulation.

This function returns a `Dict` with the valid subpopulations as the keys (`String`), and a `DataFrame` (the age report for the subpopulation) as the value. In case the function issues a warning, its return value will be an empty dictionary.
"""
function subpopulationAgeReport( mpSim::MPsim, timeGrid::Vector{Float64},
    ageRes::Real, ageType::Symbol,
    subpopulations::Subpopulation... )::Dict{String,DataFrame}

    result = Dict{String,DataFrame}()

    if now( mpSim ) == 0
        @warn "Simulation hasn't started yet, can't make report."
        return result
    end  # if now( mpSim ) == 0

    timeGrid = timeGrid[0.0 .<= timeGrid .<= now( mpSim )]
    timeGrid = unique( sort( timeGrid, rev = true ) )

    if isempty( timeGrid )
        @warn "No valid time points in time grid, cannot generate report."
        return result
    end  # if isempty( timeGrid )

    if timeGrid[end] > 0.0
        push!( timeGrid, 0.0 )
    end  # if timeGrid[end] > 0.0

    reverse!( timeGrid )

    if ageRes <= 0
        @warn "Age resolution must be > 0, cannot generate report."
        return result
    end  # if ageRes <= 0

    if ageType ∉ [:age, :tenure, :timeInNode]
        @warn "Unknown type of data requested, cannot generate report."
        return result
    end  # if ageType ∉ [:age, :tenure, :timeInNode]

    subpopulations = filter( collect( subpopulations ) ) do subpopulation
        sourceNode = subpopulation.sourceNode
        return ( sourceNode == "active" ) ||
            haskey( mpSim.baseNodeList, sourceNode ) ||
            haskey( mpSim.compoundNodeList, sourceNode )
    end  # filter( ... ) do subpopulation

    if isempty( subpopulations )
        return result
    end  # if isempty( subpopulations )

    personnelAges = map( tPoint -> getAgesAtTime( mpSim, tPoint, subpopulations,
        ageType ), timeGrid )
    # ! Vector broadcast doesn't work properly for some reason.
    personnelAges = hcat( personnelAges... )
    result = processSubpopulationAges( personnelAges, subpopulations, timeGrid,
        Float64( ageRes ) )
    return result

end  # subpopulationAgeReport( mpSim, timeGrid, ageRes, ageType,
     #   subpopulations )

"""
```
subpopulationAgeReport(
    mpSim::MPsim,
    timeRes::Real,
    ageRes::Real,
    ageType::Symbol,
    subpopulations::Subpopulation... )
```
This function generates reports for the distribution of the ages of the personnel members in the valid subpopulations in `subpopulations` of the manpower simulation `mpSim`. The reports are generated on a grid of time points with resolution `timeRes` and on a grid of ages with resoluion `ageRes`.

The type of age that is reported is determined by `ageType`, which can take these valid values:
* `:age` for the actual age of the personnel members,
* `:tenure` for the tenure of the personnel members, and
* `:timeInNode` for the time the personnel members are in their current base node. If the subpopulation is based on the entire population, this reports the tenure instead.

This function will issue a warning and not generate any report in the following two cases:
1. There are no positive time points in the time grid;
2. None of the entered subpopulations are valid subpopulations in the simulation, where a subpopulation is invalid if its souce nodes doesn't exist in the simulation.

This function returns a `Dict` with the valid subpopulations as the keys (`String`), and a `DataFrame` (the age report for the subpopulation) as the value. In case the function issues a warning, its return value will be an empty dictionary.
"""
subpopulationAgeReport( mpSim::MPsim, timeRes::Real, ageRes::Real,
    ageType::Symbol,
    subpopulations::Subpopulation... )::Dict{String,DataFrame} =
    subpopulationAgeReport( mpSim, generateTimeGrid( mpSim, timeRes ), ageRes,
    ageType, subpopulations... )

include( joinpath( repPrivPath, "subpopulationage.jl" ) )


text:
# -----------------------------------------------
# Linear Regression routines
# -----------------------------------------------

"""
Linear regression function
"""

function linregress(X::Array{T}, Y::Array{T}, regtype=OLS(); fast=false) where T<:Number
    if isdefined(regtype, :intercept)
        X = regtype.intercept ? addintercept(X) : nointercept(X)
    end
    fitted_beta, iter = regtype(X, Y)
    residuals = Y - aprod(X, fitted_beta)
    if regtype.robust == false || typeof(regtype.robust) == BasicVariance
        variance = BasicVariance(X, residuals, length(fitted_beta))
    elseif typeof(regtype.robust) == HCEVariance
        variance = HCEVariance(X, residuals)
    end
    if fast
        return (fitted_beta, variance)
    elseif iter == 0
        # The regression algorithm is not iterative:
        return linearfitresult(regtype, X, Y, fitted_beta, residuals, variance, regtype.robust)
    else
        # The regression algorithm is iterative:
        return linearfitresult(regtype, X, Y, fitted_beta, residuals, variance, regtype.robust, iter)
    end
end

function linregress(X::Array{T1}, Y::Array{T2}, regtype=OLS(); fast=false) where {T1<:Number, T2<:Number}
    X, Y = promote(X, Y)
    return linregress(X, Y, regtype; fast=fast)
end

# #########   OLS   #########

function (olsreg::OLS)(X::Array{T}, Y::Array{T}) where {T<:Number}
    # Simple OLS form:
    return (X' * X) \ (X' * Y), 0
end

# #########   GLS   #########

function (glsreg::GLS)(X::Array{T}, Y::Array{T}, omega::Array{T}) where {T<:Number}
    # Use the supplied Ω:
    return (X' * inv(omega) * X) \ (X' * inv(omega) * Y), 0
end

# #########   FGLS   #########

function (fglsreg::FGLS)(X::Array{T}, Y::Array{T}) where {T<:Number}
    # (1) Execute OLS regression to compute residuals:
    residuals = Y - X * ((X' * X) \ (X' * Y))
    # (2) Use residuals to estimate Ω:
    omega = diagm(squeeze(residuals, 2).^2)
    return (X' * inv(omega) * X) \ (X' * inv(omega) * Y), 0
end

# #########   Iterated FGLS   #########

function (iteratedfglsreg::IteratedFGLS)(X::Array{T}, Y::Array{T}) where {T<:Number}
    itfgls_beta = (X' * X) \ (X' * Y)
    residuals = Y - aprod(X, itfgls_beta)
    # (2) Use residuals to estimate Ω iteratively:
    iter::Int = 0
    for it = 1:iteratedfglsreg.iterations
        # println("it[$it] ")
        iter += 1
        omega = diagm(squeeze(residuals, 2).^2)
        try
            itfgls_beta = (X' * inv(omega) * X) \ (X' * inv(omega) * Y)
            # println("\nitgls_beta[$itfgls_beta]")
            residuals = Y - aprod(X, itfgls_beta)
        catch
            break
        end
    end
    return itfgls_beta, iter
end


text:
"""
    add_hash_cache!(graph, endpoints=[], uncacheable=[]
                    [; compression=DEFAULT_COMPRESSION, cachedir=DEFAULT_CACHE_DIR])

Optimizes a delayed execution graph `graph::DispatchGraph`
by wrapping individual nodes in load-from-disk on execute-and-store
wrappers depending on the state of the disk cache and of the graph.
The function modifies inplace the input graph and should be called on
the same unmodified `graph` after each execution and *not* on the
modified graph.

Once the original graph is modified, calling `run!` on it will,
if the cache is already present, load the top most consistent key
or alternatively re-run and store the outputs of nodes which have new state.

# Arguments
  * `exec::Executor` the `Dispatcher.jl` executor
  * `graph::DispatchGraph` input dispatch graph
  * `endpoints::AbstractVector` leaf nodes for which caching will occur;
nodes that depend on these will not be cached. The nodes can be specified
either by label of by the node object itself
  * uncacheable::AbstractVector` nodes that will never be cached and
will always be executed (these nodes are still hashed and their hashes
influence upstream node hashes as well)

# Keyword arguments
  * `compression::String` enables compression of the node outputs.
Available options are `"none"`, for no compression, `"bz2"` or `"bzip2"`
for BZIP compression and `"gz"` or `"gzip"` for GZIP compression
  * `cachedir::String` the cache directory.

Note: This function should be used with care as it modifies the input
      dispatch graph. One way to handle this is to make a function that
      generates the dispatch graph and calling `add_hash_cache!` each time on
      the distict, functionally identical graphs.
"""
function add_hash_cache!(graph::DispatchGraph,
                         endpoints::AbstractVector=[],
                         uncacheable::AbstractVector=[];
                         compression::String=DEFAULT_COMPRESSION,
                         cachedir::String=DEFAULT_CACHE_DIR)
    # Checks
    if isempty(endpoints)
        @warn "No enpoints for graph, will not process dispatch graph."
        return nothing
    end

    # Initializations
    _endpoints = map(n->get_node(graph, n), endpoints)
    _uncacheable = imap(n->get_node(graph, n), uncacheable)
    subgraph = Dispatcher.subgraph(graph, _endpoints)
    work = collect(nodes(graph))                # nodes to be traversed
    solved = Set{DispatchNode}()                # computable tasks
    node2hash = Dict{DispatchNode, String}()    # node => hash mapping
    storable = Set{String}()                    # hashes of nodes with storable output
    updates = Dict{DispatchNode, DispatchNode}()

    # Load hashchain
    hashchain = load_hashchain(cachedir, compression=compression)
    # Traverse dispatch graph
    while !isempty(work)
        node = popfirst!(work)
        deps = dependencies(node)
        if isempty(deps) || issubset(deps, solved)
            # Node is solvable
            push!(solved, node)
            _hash_node, _hash_comp = node_hash(node, node2hash)
            node2hash[node] = _hash_node
            skipcache = node in _uncacheable || !(node in subgraph.nodes)
            # Wrap nodes
            if _hash_node in keys(hashchain) && !skipcache &&
                    !(_hash_node in storable)
                # Hash match and output cacheable
                wrap_to_load!(updates, node, _hash_node,
                              cachedir=cachedir,
                              compression=compression)
            elseif _hash_node in keys(hashchain) && skipcache
                # Hash match and output *non-cachable*
                wrap_to_store!(updates, node, _hash_node,
                               cachedir=cachedir,
                               compression=compression,
                               skipcache=skipcache)
            else
                # Hash miss
                hashchain[_hash_node] = _hash_comp
                push!(storable, _hash_node)
                wrap_to_store!(updates, node, _hash_node,
                               cachedir=cachedir,
                               compression=compression,
                               skipcache=skipcache)
            end
        else
            # Non-solvable node
            push!(work, node)
        end
    end
    # Update graph
    for i in 1:length(graph.nodes)
        graph.nodes[i] = updates[graph.nodes[i]]
    end
    # Write hashchain
    store_hashchain(hashchain, cachedir, compression=compression)
    return updates
end

"""
    run!(exec, graph, endpoints, uncacheable=[]
         [;compression=DEFAULT_COMPRESSION, cachedir=DEFAULT_CACHE_DIR])

Runs the `graph::DispatchGraph` and loads or executes and stores the
outputs of the nodes in the subgraph whose leaf nodes are given by
`endpoints`. Nodes in `uncachable` are not locally cached.

# Arguments
  * `exec::Executor` the `Dispatcher.jl` executor
  * `graph::DispatchGraph` input dispatch graph
  * `endpoints::AbstractVector` leaf nodes for which caching will occur;
nodes that depend on these will not be cached. The nodes can be specified
either by label of by the node object itself
  * `uncacheable::AbstractVector` nodes that will never be cached and will
always be executed (these nodes are still hashed and their hashes influence
upstream node hashes as well)

# Keyword arguments
  * `compression::String` enables compression of the node outputs.
Available options are `"none"`, for no compression, `"bz2"` or `"bzip2"`
for BZIP compression and `"gz"` or `"gzip"` for GZIP compression
  * `cachedir::String` The cache directory.

# Examples
```julia
julia> using Dispatcher
       using DispatcherCache

       # Some functions
       foo(x) = begin sleep(1); x end
       bar(x) = begin sleep(1); x+1 end
       baz(x,y) = begin sleep(1); x-y end

       # Make a dispatch graph out of some operations
       op1 = @op foo(1)
       op2 = @op bar(2)
       op3 = @op baz(op1, op2)
       D = DispatchGraph(op3)
# DispatchGraph({3, 2} directed simple Int64 graph,
# NodeSet(DispatchNode[
# Op(DeferredFuture at (1,1,241),baz,"baz"),
# Op(DeferredFuture at (1,1,239),foo,"foo"),
# Op(DeferredFuture at (1,1,240),bar,"bar")]))

julia> # First run, writes results to disk (lasts 2 seconds)
       result_node = [op3]  # the node for which we want results
       cachedir = "./__cache__"  # directory does not exist
       @time r = run!(AsyncExecutor(), D,
                      result_node, cachedir=cachedir)
       println("result (first run) = \$(fetch(r[1].result.value))")
# [info | Dispatcher]: Executing 3 graph nodes.
# [info | Dispatcher]: Node 1 (Op<baz, Op<foo>, Op<bar>>): running.
# [info | Dispatcher]: Node 2 (Op<foo, Int64>): running.
# [info | Dispatcher]: Node 3 (Op<bar, Int64>): running.
# [info | Dispatcher]: Node 2 (Op<foo, Int64>): complete.
# [info | Dispatcher]: Node 3 (Op<bar, Int64>): complete.
# [info | Dispatcher]: Node 1 (Op<baz, Op<foo>, Op<bar>>): complete.
# [info | Dispatcher]: All 3 nodes executed.
#   2.029992 seconds (11.53 k allocations: 1.534 MiB)
# result (first run) = -2

julia> # Secod run, loads directly the result from ./__cache__
       @time r = run!(AsyncExecutor(), D,
                      [op3], cachedir=cachedir)
       println("result (second run) = \$(fetch(r[1].result.value))")
# [info | Dispatcher]: Executing 1 graph nodes.
# [info | Dispatcher]: Node 1 (Op<baz>): running.
# [info | Dispatcher]: Node 1 (Op<baz>): complete.
# [info | Dispatcher]: All 1 nodes executed.
#   0.005257 seconds (2.57 k allocations: 478.359 KiB)
# result (second run) = -2

julia> readdir(cachedir)
# 2-element Array{String,1}:
#  "cache"
#  "hashchain.json"
```
"""
function run!(exec::Executor,
              graph::DispatchGraph,
              endpoints::AbstractVector,
              uncacheable::AbstractVector=[];
              compression::String=DEFAULT_COMPRESSION,
              cachedir::String=DEFAULT_CACHE_DIR)
    # Make a copy of the input graph that will be modified
    # and mappings from the original nodes to the copies
    tmp_graph = Base.deepcopy(graph)
    node2tmpnode = Dict(graph.nodes[i] => tmp_graph.nodes[i]
                        for i in 1:length(graph.nodes))
    # Construct the endpoints and uncachable node lists
    # for the copied dispatch graph
    tmp_endpoints = [node2tmpnode[get_node(graph, node)]
                     for node in endpoints]
    tmp_uncacheable = [node2tmpnode[get_node(graph, node)]
                       for node in uncacheable]
    # Modify input graph
    updates = add_hash_cache!(tmp_graph,
                              tmp_endpoints,
                              tmp_uncacheable,
                              compression=compression,
                              cachedir=cachedir)
    # Run temporary graph
    return run!(exec, [updates[e] for e in tmp_endpoints])
end


text:
function zero_closure(sys::MomentEquations, binary_vars::Array{Int,1}=Int[])

    closure = OrderedDict()
    closure_exp = OrderedDict()

    if !isempty(binary_vars)
        sys = bernoulli_moment_eqs(sys, binary_vars)
    end

    if typeof(sys) == CentralMomentEquations
        for i in sys.iter_q
            closure[sys.M[i]] = 0
            closure_exp[sys.M[i]] = 0
        end
    else
        μ = copy(sys.μ)
        μ_symbolic = copy(sys.μ)
        raw_to_central = raw_to_central_moments(sys.N, sys.q_order)

        unique_iter_q = unique(sort(i) for i in sys.iter_q)
        sub = Dict()

        for i in unique_iter_q
            μ[i] = -(raw_to_central[i]-μ[i])
            closure[sys.μ[i]] = simplify(μ[i], expand=true)
            μ[i] = simplify(substitute(μ[i], closure_exp), expand=true)
            closure_exp[sys.μ[i]] = μ[i]

            perms = collect(multiset_permutations(i, length(i)))[2:end]

            for iter_perm in perms

                iter_perm_ind = sortperm(sortperm(iter_perm))
                for r in sys.iter_all
                    sub[sys.μ[r]] = sys.μ[r[iter_perm_ind]]
                end

                iter_perm = Tuple(iter_perm)
                closure[sys.μ[iter_perm]] = substitute(closure[sys.μ[i]], sub)
                closure_exp[sys.μ[iter_perm]] = substitute(closure_exp[sys.μ[i]], sub)
            end

        end
    end

    close_eqs(sys, closure_exp, closure, true)

end