text:
using RecipesBase
import RecipesBase.apply_recipe

using LazySets.Approximations: overapproximate, PolarDirections

# global values
DEFAULT_COLOR = :auto
DEFAULT_ALPHA = 0.5
DEFAULT_LABEL = ""
DEFAULT_GRID = true
DEFAULT_ASPECT_RATIO = :none
PLOT_PRECISION = 1e-3
PLOT_POLAR_DIRECTIONS = 40
DEFAULT_PLOT_LIMIT = 1000

function _extract_limits(p::RecipesBase.AbstractPlot)
    lims = Dict()
    if length(p) > 0
        subplot = p[1]
        for symbol in [:x, :y]
            lims[symbol] = subplot[Symbol(symbol,:axis)][:lims]
        end
    else
        lims[:x] = :auto
        lims[:y] = :auto
    end
    return lims
end

function _extract_extrema(p::RecipesBase.AbstractPlot)
    extrema = Dict()
    if length(p) > 0
        subplot = p[1]
        for symbol in [:x, :y]
            emin = subplot[Symbol(symbol,:axis)][:extrema].emin
            emax = subplot[Symbol(symbol,:axis)][:extrema].emax
            extrema[symbol] = (emin, emax)
        end
    else
        extrema[:x] = (-Inf, Inf)
        extrema[:y] = (-Inf, Inf)
    end
    return extrema
end

function _update_plot_limits!(lims, X::LazySet)
    box = box_approximation(X)
    box_min = low(box)
    box_max = high(box)
    for (idx,symbol) in enumerate([:x, :y])
        if lims[symbol] != :auto
            # width of the plotting window including the new set in the `symbol` direction
            width = max(lims[symbol][2], box_max[idx]) -
                    min(lims[symbol][1], box_min[idx])

            # scaling factor for beautification
            ϵ = 0.05
            offset = width*ϵ

            # extend the current plot limits if the new set (plus a small offset) falls outside
            lims[symbol] = (min(lims[symbol][1], box_min[idx] - offset),
                            max(lims[symbol][2], box_max[idx] + offset))
        end
    end
    nothing
end

function _set_auto_limits_to_extrema!(lims, extr)
    # if the limit is :auto, set the limits to the current extrema
    for symbol in [:x, :y]
        if lims[symbol] == :auto
            emin, emax = extr[symbol]
            # if the extrema are (-Inf,Inf), i.e. an empty plot,
            # set it to (-1.,1.)
            lmin = isinf(emin) ? -1.0 : emin
            lmax = isinf(emax) ? 1.0 : emax
            lims[symbol] = (lmin, lmax)
        end
    end
    # otherwise keep the old limits
    nothing
end

function _bounding_hyperrectangle(lims, N)
    low_lim = [lims[:x][1] - DEFAULT_PLOT_LIMIT, lims[:y][1] - DEFAULT_PLOT_LIMIT]
    high_lim = [lims[:x][2] + DEFAULT_PLOT_LIMIT, lims[:y][2] + DEFAULT_PLOT_LIMIT]
    return Hyperrectangle(low=convert.(N,low_lim), high=convert.(N,high_lim))
end

"""
    plot_list(list::AbstractVector{VN}, [ε]::N=N(PLOT_PRECISION),
              [Nφ]::Int=PLOT_POLAR_DIRECTIONS, [fast]::Bool=false; ...)
        where {N, VN<:LazySet{N}}

Plot a list of convex sets.

### Input

- `list` -- list of convex sets (1D or 2D)
- `ε`    -- (optional, default: `PLOT_PRECISION`) approximation error bound
- `Nφ`   -- (optional, default: `PLOT_POLAR_DIRECTIONS`) number of polar
            directions (used to plot lazy intersections)
- `fast` -- (optional, default: `false`) switch for faster plotting but without
            individual plot recipes (see notes below)

### Notes

For each set in the list we apply an individual plot recipe.

The option `fast` provides access to a faster plotting scheme where all sets in
the list are first converted to polytopes and then plotted in one single run.
This, however, is not suitable when plotting flat sets (line segments,
singletons) because then the polytope plot recipe does not deliver good results.
Hence by default we do not use this option.
For plotting a large number of (non-flat) polytopes, we highly advise activating
this option.

### Examples

```julia
julia> B1 = BallInf(zeros(2), 0.4);

julia> B2 = BallInf(ones(2), 0.4);

julia> plot([B1, B2])
```

Some of the sets in the list may not be plotted precisely but rather
overapproximated first.
The second argument `ε` controls the accuracy of this overapproximation.

```julia
julia> Bs = [BallInf(zeros(2), 0.4), Ball2(ones(2), 0.4)];

julia> plot(Bs, 1e-3)  # default accuracy value (explicitly given for clarity)

julia> plot(Bs, 1e-2)  # faster but less accurate than the previous call
```
"""
@recipe function plot_list(list::AbstractVector{VN}, ε::N=N(PLOT_PRECISION),
                           Nφ::Int=PLOT_POLAR_DIRECTIONS, fast::Bool=false
                          ) where {N, VN<:LazySet{N}}
    if fast
        label --> DEFAULT_LABEL
        grid --> DEFAULT_GRID
        if DEFAULT_ASPECT_RATIO != :none
            aspect_ratio --> DEFAULT_ASPECT_RATIO
        end
        seriesalpha --> DEFAULT_ALPHA
        seriescolor --> DEFAULT_COLOR
        seriestype --> :shape

        first = true
        x = Vector{N}()
        y = Vector{N}()
        for Xi in list
            if Xi isa Intersection
                res = plot_recipe(Xi, ε, Nφ)
            else
                # hard-code overapproximation here to avoid individual
                # compilations for mixed sets
                Pi = overapproximate(Xi, ε)
                vlist = transpose(hcat(convex_hull(vertices_list(Pi))...))
                if isempty(vlist)
                    @warn "overapproximation during plotting was empty"
                    continue
                end
                res = vlist[:, 1], vlist[:, 2]
                if length(res[1]) > 2
                    # add first vertex to "close" the polygon
                    push!(res[1], vlist[1, 1])
                    push!(res[2], vlist[1, 2])
                end
            end
            if isempty(res)
                continue
            else
                x_new, y_new = res
            end
            if first
                first = false
            else
                push!(x, N(NaN))
                push!(y, N(NaN))
            end
            append!(x, x_new)
            append!(y, y_new)
        end
        x, y
    else
        for Xi in list
            if Xi isa Intersection
                @series Xi, ε, Nφ
            else
                @series Xi, ε
            end
        end
    end
end

"""
    plot_lazyset(X::LazySet{N}, [ε]::N=N(PLOT_PRECISION); ...) where {N}

Plot a convex set.

### Input

- `X` -- convex set
- `ε` -- (optional, default: `PLOT_PRECISION`) approximation error bound

### Notes

See [`plot_recipe(::LazySet)`](@ref).

For polyhedral set types (subtypes of `AbstractPolyhedron`), the argument `ε` is
ignored.

### Examples

```julia
julia> B = Ball2(ones(2), 0.1);

julia> plot(B, 1e-3)  # default accuracy value (explicitly given for clarity)

julia> plot(B, 1e-2)  # faster but less accurate than the previous call
```
"""
@recipe function plot_lazyset(X::LazySet{N}, ε::N=N(PLOT_PRECISION)) where {N}
    if dim(X) == 1
        plot_recipe(X, ε)
    else
        label --> DEFAULT_LABEL
        grid --> DEFAULT_GRID
        if DEFAULT_ASPECT_RATIO != :none
            aspect_ratio --> DEFAULT_ASPECT_RATIO
        end
        seriesalpha --> DEFAULT_ALPHA
        seriescolor --> DEFAULT_COLOR

        # extract limits and extrema of already plotted sets
        p = plotattributes[:plot_object]
        lims = _extract_limits(p)
        extr = _extract_extrema(p)

        if !isbounded(X)
            _set_auto_limits_to_extrema!(lims, extr)
            X = intersection(X, _bounding_hyperrectangle(lims, eltype(X)))

        # if there is already a plotted set and the limits are fixed,
        # automatically adjust the axis limits (e.g. after plotting a unbounded set)
        elseif length(p) > 0
            _update_plot_limits!(lims, X)
        end

        xlims --> lims[:x]
        ylims --> lims[:y]

        res = plot_recipe(X, ε)
        if isempty(res)
            res
        else
            x, y = res
            if length(x) == 1 ||
                    (length(x) == 2 && norm([x[1], y[1]] - [x[2], y[2]]) ≈ 0)
                # single point
                seriestype := :scatter
            elseif length(x) == 2
                # flat line segment
                linecolor   --> DEFAULT_COLOR
                markercolor --> DEFAULT_COLOR
                markershape --> :circle
                seriestype := :path
            else
                seriestype := :shape
            end
            x, y
        end
    end
end

"""
    plot_singleton(S::AbstractSingleton{N}, [ε]::N=zero(N); ...) where {N}

Plot a singleton.

### Input

- `S` -- singleton
- `ε` -- (optional, default: `0`) ignored, used for dispatch

### Examples

```julia
julia> plot(Singleton([0.5, 1.0]))
```
"""
@recipe function plot_singleton(S::AbstractSingleton{N}, ε::N=zero(N)) where {N}
    label --> DEFAULT_LABEL
    grid --> DEFAULT_GRID
    if DEFAULT_ASPECT_RATIO != :none
        aspect_ratio --> DEFAULT_ASPECT_RATIO
    end
    seriesalpha --> DEFAULT_ALPHA
    seriescolor --> DEFAULT_COLOR
    seriestype := :scatter

    # update manually set plot limits if necessary
    p = plotattributes[:plot_object]
    if length(p) > 0
        lims = _extract_limits(p)
        _update_plot_limits!(lims, S)
        xlims --> lims[:x]
        ylims --> lims[:y]
    end

    plot_recipe(S, ε)
end

"""
    plot_emptyset(∅::EmptySet, [ε]::N=zero(N); ...)

Plot an empty set.

### Input

- `∅` -- empty set
- `ε` -- (optional, default: `0`) ignored, used for dispatch
"""
@recipe function plot_emptyset(∅::EmptySet{N}, ε::N=zero(N)) where {N}
    label --> DEFAULT_LABEL
    grid --> DEFAULT_GRID
    if DEFAULT_ASPECT_RATIO != :none
        aspect_ratio --> DEFAULT_ASPECT_RATIO
    end

    plot_recipe(∅)
end

"""
    plot_intersection(cap::Intersection{N}, [ε]::N=zero(N),
                      [Nφ]::Int=PLOT_POLAR_DIRECTIONS) where {N}

Plot a lazy intersection.

### Input

- `cap`  -- lazy intersection
- `ε`    -- (optional, default `0`) ignored, used for dispatch
- `Nφ`   -- (optional, default: `PLOT_POLAR_DIRECTIONS`) number of polar
            directions used in the template overapproximation

### Notes

This function is separated from the main `LazySet` plot recipe because iterative
refinement is not available for lazy intersections (since it uses the support
vector (but see
[#1187](https://github.com/JuliaReach/LazySets.jl/issues/1187))).

Also note that if the set is a *nested* intersection, you may have to manually
overapproximate this set before plotting (see
`LazySets.Approximations.overapproximate` for details).

### Examples

```julia
julia> using LazySets.Approximations

julia> X = Ball2(zeros(2), 1.) ∩ Ball2(ones(2), 1.5);  # lazy intersection

julia> plot(X)
```

You can specify the accuracy of the overapproximation of the lazy intersection
by passing a higher value for `Nφ`, which stands for the number of polar
directions used in the overapproximation.
This number can also be passed to the `plot` function directly.

```julia
julia> plot(overapproximate(X, PolarDirections(100)))

julia> plot(X, -1., 100)  # equivalent to the above line
```
"""
@recipe function plot_intersection(cap::Intersection{N},
                                   ε::N=zero(N),
                                   Nφ::Int=PLOT_POLAR_DIRECTIONS
                                  ) where {N}
    label --> DEFAULT_LABEL
    grid --> DEFAULT_GRID
    if DEFAULT_ASPECT_RATIO != :none
        aspect_ratio --> DEFAULT_ASPECT_RATIO
    end
    seriesalpha --> DEFAULT_ALPHA
    seriescolor --> DEFAULT_COLOR
    seriestype := :shape

    # extract limits and extrema of already plotted sets
    p = plotattributes[:plot_object]
    lims = _extract_limits(p)
    extr = _extract_extrema(p)

    if !isbounded(cap)
        _set_auto_limits_to_extrema!(lims, extr)
        bounding_box = _bounding_hyperrectangle(lims, eltype(cap))
        if !isbounded(cap.X)
            bounded_X = Intersection(bounding_box, cap.X)
            cap = Intersection(bounded_X, cap.Y)
        else
            cap = Intersection(bounding_box, cap.Y)
        end

    # if there is already a plotted set and the limits are fixed,
    # automatically adjust the axis limits (e.g. after plotting a unbounded set)
    elseif length(p) > 0
        _update_plot_limits!(lims, cap)
    end

    xlims --> lims[:x]
    ylims --> lims[:y]

    plot_recipe(cap, ε)
end

# non-convex sets

@recipe function plot_union(cup::UnionSet{N}, ε::N=zero(N)) where {N}
    @series [cup.X, cup.Y], ε
end

@recipe function plot_union(cup::UnionSetArray{N}, ε::N=zero(N)) where {N}
    @series array(cup), ε
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

function f28825()
    eval(:(1+1))
end

ext: jl

max_stars_repo_name: TimoLarson/julia

max_issues_repo_name: AsForMe/julia

max_forks_repo_name: AsForMe/julia


text:
module TimeseriesSurrogates

using Distributions
using StatsBase
using InplaceOps
using AbstractFFTs
using DSP
using Interpolations
using Wavelets
using Requires
using Plots

# Example systems
include("testsystems.jl")

# Periodogram interpolation
include("interpolation.jl")

# Visualization routine for time series + surrogate + periodogram/acf/histogram
include("surrogate_plot.jl")

# The different surrogate routines
include("randomshuffle.jl")
include("randomphases.jl")
include("randomamplitudes.jl")
include("aaft.jl")
include("iaaft.jl")
include("wiaaft.jl")

# Define and export plot routines for all combinations of example processes and surrogate
# types
processes = (:AR1, :NSAR2, :NLNS, :randomwalk, :SNLST)
surrogate_methods = (:randomshuffle, :randomphases, :randomamplitudes, :aaft, :iaaft)
include("plots_and_anim.jl")

include("uncertaindatasets.jl")

export NLNS, NSAR2, AR1, randomwalk, SNLST,
        randomshuffle, randomamplitudes, randomphases, aaft, iaaft, wiaaft

end # module


text:
# Reference: Chirp Z-Transform Spectral Zoom Optimization with MATLAB
# http://prod.sandia.gov/techlib/access-control.cgi/2005/057084.pdf

# For spectral zoom, set:
#   W = exp(-j*2*pi*(f2-f1)/(m*fs));
#   A = exp(j*2*pi*f1/fs);
# where:
#   f1 = start freq
#   f2 = end freq
#   m  = length of x
#   fs = sampling rate
function czt(x::Vector{Complex128}, m::Int, w::Complex128, a::Complex128)
    # TODO: add argument valdiation
    # TODO: figure out why output isn't matching FFT
    n = length(x)
    N = [0:n-1].+n
    NM = [-(n-1):(m-1)].+n
    M = [0:m-1].+n

    nfft = nextpow2(n+m-1)
    W2 = w.^(([-(n-1):max(m-1,n-1)].^2)/2)

    fg = zeros(Complex128, nfft)
    fg[1:n] = x.*(a.^-(N.-n)).*W2[N]
    fg = fft(fg)

    fw = zeros(Complex128, nfft)
    fw[1:length(NM)] = 1./W2[NM]
    fw = fft(fw)
    gg = ifft(fg.*fw)

    return gg[M].*W2[M]
end

function czt(x::Vector{Complex128}, m::Int)
    # TODO: add argument valdiation
    czt( x, m, exp(-im*2*pi/m), 1.0+0.0im)
end


text:
using FastRef
using LazySets
import FastRef: forward_network, forward_affine_map, ishull

nnet = read_nnet("nnet/86442.nnet")

delta = 0.25

solver1 = MaxSens(delta)
solver2 = HullGrid(delta)
solver3 = DimGrid(delta)
solver4 = FastGrid(delta)
solver5 = SpeGuid(delta)
solver6 = HullTree(delta)
solver7 = DimTree(delta)
solver8 = FastTree(delta)

in_hyper = Hyperrectangle(fill(1.0, 8), fill(1.0, 8))
out_hyper = Hyperrectangle(fill(0.0, 2), fill(10.0, 2))
problem = Problem(nnet, in_hyper, out_hyper)

#solver1

time1 = 0

solve(solver1, problem)
for i = 1:10
    timed_result =@timed solve(solver1, problem)
    print("MaxSens - test " * string(i) * " - Time: " * string(timed_result.time) * " s")
    print(" - Output: " * string(timed_result.value) * "\n")
    global time1 += timed_result.time
end

print("Average time: " * string(time1/10) * " s\n\n")

#solver2

time2 = 0

solve(solver2, problem)
for i = 1:10
    timed_result =@timed solve(solver2, problem)
    print("HullGrid - test " * string(i) * " - Time: " * string(timed_result.time) * " s")
    print(" - Output: " * string(timed_result.value) * "\n")
    global time2 += timed_result.time
end

print("Average time: " * string(time2/10) * " s\n\n")

#solver3

time3 = 0

solve(solver3, problem)
for i = 1:10
    timed_result =@timed solve(solver3, problem)
    print("DimGrid - test " * string(i) * " - Time: " * string(timed_result.time) * " s")
    print(" - Output: " * string(timed_result.value) * "\n")
    global time3 += timed_result.time
end

print("Average time: " * string(time3/10) * " s\n\n")

#solver4

time4 = 0

solve(solver4, problem)
for i = 1:10
    timed_result =@timed solve(solver4, problem)
    print("FastGrid - test " * string(i) * " - Time: " * string(timed_result.time) * " s")
    print(" - Output: " * string(timed_result.value) * "\n")
    global time4 += timed_result.time
end

print("Average time: " * string(time4/10) * " s\n\n")

#solver5

time5 = 0

solve(solver5, problem)
for i = 1:10
    timed_result =@timed solve(solver5, problem)
    print("SpeGuid - test " * string(i) * " - Time: " * string(timed_result.time) * " s")
    print(" - Output: " * string(timed_result.value) * "\n")
    global time5 += timed_result.time
end

print("Average time: " * string(time5/10) * " s\n\n")

#solver6

time6 = 0

solve(solver6, problem)
for i = 1:10
    timed_result =@timed solve(solver6, problem)
    print("HullTree - test " * string(i) * " - Time: " * string(timed_result.time) * " s")
    print(" - Output: " * string(timed_result.value) * "\n")
    global time6 += timed_result.time
end

print("Average time: " * string(time6/10) * " s\n\n")

#solver7

time7 = 0

solve(solver7, problem)
for i = 1:10
    timed_result =@timed solve(solver7, problem)
    print("DimTree - test " * string(i) * " - Time: " * string(timed_result.time) * " s")
    print(" - Output: " * string(timed_result.value) * "\n")
    global time7 += timed_result.time
end

print("Average time: " * string(time7/10) * " s\n\n")

#solver8

time8 = 0

solve(solver8, problem)
for i = 1:10
    timed_result =@timed solve(solver8, problem)
    print("FastTree - test " * string(i) * " - Time: " * string(timed_result.time) * " s")
    print(" - Output: " * string(timed_result.value) * "\n")
    global time8 += timed_result.time
end

print("Average time: " * string(time8/10) * " s\n\n")


text:
using Pkg

Pkg.add("CSV")
Pkg.add("DataFrames")
Pkg.add("Symbolics")
Pkg.add("Optim")
Pkg.add("LineSearches")
Pkg.add("BenchmarkTools")
Pkg.add("Statistics")


text:
# getindex

function Base.get(f::Function, t::NDSparse{T,D}, idxs::D) where {T,D<:Tuple}
    flush!(t)
    i = searchsorted(t.index, convertkey(t, idxs))
    return length(i) == 1 ? t.data[first(i)] : f()
end

function Base.get(t::NDSparse{T,D}, idxs::D, default) where {T,D<:Tuple}
    get(t, idxs) do
        default
    end
end

getindex(t::NDSparse, idxs...) = (flush!(t); _getindex(t, idxs))

_getindex(t::NDSparse{T,D}, idxs::D) where {T,D<:Tuple} = _getindex_scalar(t, idxs)
_getindex(t::NDSparse, idxs::Tuple{Vararg{Real}}) = _getindex_scalar(t, idxs)

function _getindex_scalar(t, idxs)
    get(t, idxs) do
        throw(KeyError(idxs))
    end
end

function Base.haskey(t::NDSparse{T, D}, idxs::D) where {T,D<:Tuple}
    flush!(t)
    i = searchsorted(t.index, convertkey(t, idxs))
    return length(i) == 1
end

# branch instead of diagonal dispatch to avoid ambiguities
_in(x, y) = isa(x,typeof(y)) ? isequal(x, y) : in(x, y)
_in(x, ::Colon) = true
_in(x, v::AbstractVector) = (idx=searchsortedfirst(v, x); idx<=length(v) && v[idx]==x)
_in(x, v::AbstractString) = x == v
_in(x, v::Symbol) = x === v
_in(x, v::Number) = isequal(x, v)

# test whether row r is within product(idxs...)
@inline row_in(cs, r::Integer, idxs) = _row_in(cs[1], r, idxs[1], tail(cs), tail(idxs))
@inline _row_in(c1, r, i1, rI, ri) = _in(c1[r],i1) & _row_in(rI[1], r, ri[1], tail(rI), tail(ri))
@inline _row_in(c1, r, i1, rI::Tuple{}, ri) = _in(c1[r],i1)

range_estimate(col, idx) = 1:length(col)
range_estimate(col::AbstractVector{T}, idx::T) where {T} = searchsortedfirst(col, idx):searchsortedlast(col,idx)
range_estimate(col, idx::AbstractArray) = searchsortedfirst(col,first(idx)):searchsortedlast(col,last(idx))
range_estimate(col::Columns, idx::AbstractArray) = searchsortedfirst(col,first(idx)):searchsortedlast(col,last(idx))

const _fwd = Base.Order.ForwardOrdering()

range_estimate(col, idx, lo, hi) = 1:length(col)
range_estimate(col::AbstractVector{T}, idx::T, lo, hi) where {T} =
    searchsortedfirst(col, idx, lo, hi, _fwd):searchsortedlast(col, idx, lo, hi, _fwd)
range_estimate(col, idx::AbstractArray, lo, hi) =
    searchsortedfirst(col, first(idx), lo, hi, _fwd):searchsortedlast(col, last(idx), lo, hi, _fwd)

isconstrange(col, idx) = false
isconstrange(col::AbstractVector{T}, idx::T) where {T} = true
isconstrange(col, idx::AbstractArray) = isequal(first(idx), last(idx))

function range_estimate(I::Columns, idxs)
    r = range_estimate(columns(I)[1], idxs[1])
    i = 1; n = length(idxs)
    while i < n && isconstrange(columns(I)[i], idxs[i])
        i += 1
        r = intersect(r, range_estimate(columns(I)[i], idxs[i], first(r), last(r)))
    end
    return r
end

function _getindex(t::NDSparse, idxs)
    I = t.index
    cs = astuple(columns(I))
    if fieldcount(typeof(idxs)) !== fieldcount(typeof(columns(I)))
        error("wrong number of indices")
    end
    for idx in idxs
        isa(idx, AbstractVector) && (issorted(idx) || error("indices must be sorted for ranged/vector indexing"))
    end
    out = convert(Vector{Int32}, range_estimate(I, idxs))
    filter!(i->row_in(cs, i, idxs), out)
    keepdims = filter(i->eltype(columns(t.index)[i]) != typeof(idxs[i]), 1:length(idxs))
    NDSparse(Columns(map(x->x[out], getsubfields(columns(I), keepdims))), t.data[out], presorted=true)
end

# iterators over indices - lazy getindex

"""
`where(arr::NDSparse, indices...)`

Returns an iterator over data items where the given indices match. Accepts the
same index arguments as `getindex`.
"""
function where(d::NDSparse, idxs::Vararg{Any,N}) where N
    I = d.index
    cs = astuple(columns(I))
    data = d.data
    rng = range_estimate(I, idxs)
    (data[i] for i in Iterators.Filter(r->row_in(cs, r, idxs), rng))
end

"""
`update!(f::Function, arr::NDSparse, indices...)`

Replace data values `x` with `f(x)` at each location that matches the given
indices.
"""
function update!(f::Union{Function,Type}, d::NDSparse, idxs::Vararg{Any,N}) where N
    I = d.index
    cs = astuple(columns(I))
    data = d.data
    rng = range_estimate(I, idxs)
    for r in rng
        if row_in(cs, r, idxs)
            data[r] = f(data[r])
        end
    end
    d
end

pairs(d::NDSparse) = (d.index[i]=>d.data[i] for i in 1:length(d))

"""
`pairs(arr::NDSparse, indices...)`

Similar to `where`, but returns an iterator giving `index=>value` pairs.
`index` will be a tuple.
"""
function pairs(d::NDSparse, idxs::Vararg{Any,N}) where N
    I = d.index
    cs = astuple(columns(I))
    data = d.data
    rng = range_estimate(I, idxs)
    (I[i]=>data[i] for i in Iterators.Filter(r->row_in(cs, r, idxs), rng))
end

# setindex!

setindex!(t::NDSparse, rhs, idxs...) = _setindex!(t, rhs, idxs)

# assigning to an explicit set of indices --- equivalent to merge!

setindex!(t::NDSparse, rhs, I::Columns) = setindex!(t, fill(rhs, length(I)), I) # TODO avoid `fill`

setindex!(t::NDSparse, rhs::AbstractVector, I::Columns) = merge!(t, NDSparse(I, rhs, copy=false))

# assigning a single item

_setindex!(t::NDSparse{T,D}, rhs::AbstractArray, idxs::D) where {T,D} = _setindex_scalar!(t, rhs, idxs)
_setindex!(t::NDSparse, rhs::AbstractArray, idxs::Tuple{Vararg{Real}}) = _setindex_scalar!(t, rhs, idxs)
_setindex!(t::NDSparse{T,D}, rhs, idxs::D) where {T,D} = _setindex_scalar!(t, rhs, idxs)
#_setindex!(t::NDSparse, rhs, idxs::Tuple{Vararg{Real}}) = _setindex_scalar!(t, rhs, idxs)

function _setindex_scalar!(t, rhs, idxs)
    foreach(push!, columns(t.index_buffer), idxs)
    push!(t.data_buffer, rhs)
    t
end

# vector assignment: works like a left join

_setindex!(t::NDSparse, rhs::NDSparse, idxs::Tuple{Vararg{Real}}) = _setindex!(t, rhs.data, idxs)
_setindex!(t::NDSparse, rhs::NDSparse, idxs) = _setindex!(t, rhs.data, idxs)

function _setindex!(d::NDSparse{T,D}, rhs::AbstractArray, idxs) where {T,D}
    for idx in idxs
        isa(idx, AbstractVector) && (issorted(idx) || error("indices must be sorted for ranged/vector indexing"))
    end
    flush!(d)
    I = d.index
    data = d.data
    ll = length(I)
    p = product(idxs...)
    elem = iterate(p)
    elem === nothing && return d
    R, s = elem
    i = j = 1
    L = I[i]
    while i <= ll
        c = cmp(L, R)
        if c < 0
            i += 1
            L = I[i]
        elseif c == 0
            data[i] = rhs[j]
            i += 1
            L = I[i]
            j += 1
            elem = iterate(p, s)
            elem === nothing && break
            R, s = elem
        else
            j += 1
            elem = iterate(p, s)
            elem === nothing && break
            R, s = elem
        end
    end
    return d
end

# broadcast assignment of a single value into all matching locations

function _setindex!(d::NDSparse{T,D}, rhs, idxs) where {T,D}
    for idx in idxs
        isa(idx, AbstractVector) && (issorted(idx) || error("indices must be sorted for ranged/vector indexing"))
    end
    flush!(d)
    I = d.index
    cs = astuple(columns(I))
    data = d.data
    rng = range_estimate(I, idxs)
    for r in rng
        if row_in(cs, r, idxs)
            data[r] = rhs
        end
    end
    d
end

"""
`flush!(arr::NDSparse)`

Commit queued assignment operations, by sorting and merging the internal temporary buffer.
"""
function flush!(t::NDSparse)
    if !isempty(t.data_buffer)
        # 1. form sorted array of temp values, preferring values added later (`right`)
        temp = NDSparse(t.index_buffer, t.data_buffer, copy=false, agg=right)

        if any(isshared, _cols_tuple(keys(t)))
            t.index = copy(keys(t))
        end
        if any(isshared, _cols_tuple(values(t)))
            t.data = copy(values(t))
        end

        # 2. merge in
        _merge!(t, temp, right)

        # 3. clear buffer
        empty!(t.index_buffer)
        empty!(t.data_buffer)
    end
    nothing
end


text:
#############################
## This file extends the linear functionality in src/triVolterraFullKernel.jl
## to nonlinear and integro-differential equations. This has many different variations.
## We hard-coded some standard ones here but it is straightforward to modify it for more general cases.

#############################
##   Provides a measure of distance between guess and true solution of non-linear 2nd kind Volterra where f(u) is the nonlinearity
####
function triNonLinearVolterraObjective(uguess,f,K,g::Fun,depth)
    u = Fun(Jacobi(0,1,0..1),uguess)
    fofu = pad(Fun(x->f(u(x)),Jacobi(0,1,0..1)).coefficients,depth)
    V = Conversion(Jacobi(0,0,0..1),Jacobi(0,1,0..1))[1:depth,1:depth]*reflectPabtoPba(depth)*WLoweringP01P00(depth)*triVolterraFullKernelOpP01(K,depth,true)
    return pad(uguess,depth)-V*fofu-pad(g.coefficients,depth)
end

#############################
##   Provides a measure of distance between guess and true solution of non-linear, integro-differential
##   2nd kind Volterra where f(u) is the nonlinearity, DD is the order of the derivative and
##   Init is a vector of length D containing the initial conditions.
##
##   Currently hard-coded for 1st and 2nd derivative problems only but easy to extend to general case when needed.
####
function triNonLinearIntegroDiffVolterraObjective(uguess,f,K,g,depth,DD,init,lincomb,kerneldepth)
    u = Fun(Jacobi(0,1,0..1),uguess)
    gF = Fun(x->g(x),Jacobi(DD,DD+1, 0..1))
    fofu = pad(Fun(x->lincomb(x)+f(u(x)),Jacobi(0,1,0..1)).coefficients,depth)
    Kfun = Fun((x,y)->K(x,y),JacobiTriangle())
    if length(Kfun.coefficients) == 1    # Check if kernel is a constant, in which case explicit monomial method is used
        V = triVolterraOpP01([Kfun.coefficients[1],0.,0.,0.],depth,true)
    else
        V = triVolterraFullKernelOpP01(K,depth,true,kerneldepth)
    end
    V = reflectPabtoPba(depth)*WLoweringP01P00(depth)*V
    V = Conversion(Jacobi(0,0,0..1),Jacobi(DD,1+DD,0..1))[1:depth,1:depth]*V[1:depth,1:depth]
    if DD==1
        return pad([u(0)-init[1];pad((Derivative(DD)*u).coefficients,depth)-V*fofu-pad(gF.coefficients,depth)],depth)
    end
    if DD==2
        return pad([u(0)-init[1];(Derivative(DD-1)*u)(0)-init[2];pad((Derivative(DD)*u).coefficients,depth)-V*fofu-pad(gF.coefficients,depth)],depth)
    end
end
function triNonLinearIntegroDiffVolterraObjective(uguess,f,K,g,depth,DD,init,lincomb)
    u = Fun(Jacobi(0,1,0..1),uguess)
    gF = Fun(x->g(x),Jacobi(DD,DD+1, 0..1))
    fofu = pad(Fun(x->lincomb(x)+f(u(x)),Jacobi(0,1,0..1)).coefficients,depth)
    Kfun = Fun((x,y)->K(x,y),JacobiTriangle())
    if length(Kfun.coefficients) == 1    # Check if kernel is a constant, in which case explicit monomial method is used
        V = triVolterraOpP01([Kfun.coefficients[1],0.,0.,0.],depth,true)
    else
        V = triVolterraFullKernelOpP01(K,depth,true)
    end
    V = reflectPabtoPba(depth)*WLoweringP01P00(depth)*V
    V = Conversion(Jacobi(0,0,0..1),Jacobi(DD,1+DD,0..1))[1:depth,1:depth]*V[1:depth,1:depth]
    if DD==1
        return pad([u(0)-init[1];pad((Derivative(DD)*u).coefficients,depth)-V*fofu-pad(gF.coefficients,depth)],depth)
    end
    if DD==2
        return pad([u(0)-init[1];(Derivative(DD-1)*u)(0)-init[2];pad((Derivative(DD)*u).coefficients,depth)-V*fofu-pad(gF.coefficients,depth)],depth)
    end
end
#############################
##   This is the pochhammer symbol function from FastTransforms.jl. It's reproduced here directly for testing purposes.
####
function pochhammer(x::Number,n::Integer)
    ret = one(x)
    if n≥0
        for i=0:n-1
            ret *= x+i
        end
    else
        ret /= pochhammer(x+n,-n)
    end
    ret
end

#############################
##   alternI: Diagonal matrix with entries (-1)^n
####
function alternI(N::Integer)
    altI=BandedMatrix{Float64}(undef, (N, N), (0,0))
    for n=1:N
        altI[n,n]=(-1)^(n+1)
    end
    return altI
end

#############################
##   EvalJPat1mat: Operator which evaluates a coefficient vector of
##              0..1 Jacobi polynomials in P_n^(α,β)basis at x=0 if given α.
##              Combine with alternI(N) and use β instead of α to evaluate at x=1.
##              Be mindful that in the ApproxFun version used here (though not in later versions), we have Jacobi(β,α) = P_n^(α,β)(x)
####
function EvalJPat1mat(α,N::Integer)
    poch=BandedMatrix{Float64}(undef, (N, N), (0,0))
    for n=1:N
        poch[n,n]=FastTransforms.pochhammer(1+α,n-1)/SpecialFunctions.gamma(n);
    end
    return poch
end
function EvalatLHS(β,N)
    return Ones(1,N)*alternI(N)*EvalJPat1mat(β,N)
end
function EvalatRHS(α,N)
    return Ones(1,N)*EvalJPat1mat(α,N)
end

#############################
## Some specific commonly occuring speccial cases of the above for higher efficiency computations, evaluate LHS of specific polynomial bases
## Prime in name indicates number of derivatives. Can also just use the general function but this is a bit faster if you know what you will need.
####
function DirectEvalLHSP10at0(N::Integer)
    vek=ones(1,N)
    for n=1:N
        vek[1,n]=(-1)^(n-1)
    end
    return vek
end
function DirectEvalLHSP10atPrime0(N::Integer)
    vek=ones(1,N)
    for n=1:N
        vek[1,n]=(-1)^(n-1)*n
    end
    return vek
end
function DirectEvalLHSP10atPPrime0(N::Integer)
    vek=ones(1,N)
    for n=1:N
        vek[1,n]=(-1)^(n-1)*n*(n+1)/2
    end
    return vek
end
function DirectEvalLHSP10atPPPrime0(N::Integer)
    vek=ones(1,N)
    for n=1:N
        vek[1,n]=(-1)^(n-1)*n*(n+1)*(n+2)/6
    end
    return vek
end
function DirectEvalLHSP10atPPPPrime0(N::Integer)
    vek=ones(1,N)
    for n=1:N
        vek[1,n]=(-1)^(n-1)*(n+3)*(n+2)*(n+1)*n/24
    end
    return vek
end


text:
function duftau(t,
		fft_u,  fft_v,
		fft_du, fft_dv,
		A1, A2,
                matr, conjmatr,
		sigma, llambda, Ntaumm)

    sigma = 1

    u = ifft(exp.(1im * t .* A1) .* fft_u .* matr)
    v = ifft(exp.(1im * t .* A1) .* fft_v .* matr)

    du = ifft(exp(1im * t .* A1) .* fft_du .* matr)
    dv = ifft(exp(1im * t .* A1) .* fft_dv .* matr)

    z  = (u  .+ conj(v) ) / 2
    dz = (du .+ conj(dv)) / 2

    fz1 = 2 * abs(z) .^ 2 * dz .+ z .^ 2 .* conj(dz)

    champu = -1im * llambda * A2 .* exp.(-1im * t .* A1) .* fft(conjmatr .* fz1)
    champv = -1im * llambda * A2 .* exp.(-1im * t .* A1) .* fft(conjmatr .* conj(fz1))

    return champu, champv

end


text:
using Test

include("../nodes/node_gaussian.jl")

@testset "constructor" begin

    f = NodeGaussian("x0", "x1", 1., 1., "f1")

    @test f.node_id == "f1"
    @test f.transition == 1.0
    @test f.precision == 1.0
    @test f.edge_data_id == "x0"
    @test f.edge_mean_id == "x1"

    #TODO: test invalid inputs
end

max_stars_repo_name: wmkouw/schedulefreeMP

max_issues_repo_name: wmkouw/schedulefreeMP

max_forks_repo_name: wmkouw/schedulefreeMP


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

# INFLUENCE FUNCTION

function influence(obj::MLE, w::AbstractWeights)
    s = score(obj)
    j = jacobian(obj, w)
    return - s / j'
end

function influence(obj::TwoStageModel, w::AbstractWeights)
    c = influence(first_stage(obj), w) * crossjacobian(obj, w)'
    s = score(obj)
    j = jacobian(obj, w)
    return - (s + c) / j'
end

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

# FINITE-SAMPLE ADJUSTMENT

varcorrection(corr::Homoscedastic, w::AbstractWeights) = 1

function varcorrection(corr::Heteroscedastic, w::AbstractWeights)
    sum(w) * varcorrection(w, true)
end

varcorrection(corr::Clustered, w::AbstractWeights) = corr.n_clusters / (corr.n_clusters - 1)

function varcorrection(corr::CrossCorrelated, w::AbstractWeights)
    sum(w) * varcorrection(w, true)
end

# SET COVARIANCE MATRIX

_vcov!(obj::ParModel, V::AbstractMatrix) = (switch_stage(obj).V = V)

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

# COVARIANCE MATRIX FOR HOMOSCEDASTIC MLE

function _vcov!(obj::MLE, corr::Homoscedastic, w::AbstractWeights)
    obj.V = getcorr(obj).expected ? inv(crossprod(score(obj), w)) : - inv(jacobian(obj, w))
end

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

# COVARIANCE MATRIX FOR MLE AND TWO-STAGE MODELS

for m in (MLE, TwoStageModel)
    @eval begin

        function _vcov!(obj::$m, corr::Heteroscedastic, w::UnitWeights)
            ψ = influence(obj, w)
            V = crossprod(ψ)
            corr.corrected && lmul!(varcorrection(corr, w), V)
            _vcov!(obj, V)
        end

        function _vcov!(obj::$m, corr::Heteroscedastic, w::FrequencyWeights)
            ψ = influence(obj, w)
            V = crossprod(ψ, w)
            corr.corrected && lmul!(varcorrection(corr, w), V)
            _vcov!(obj, V)
        end

        function _vcov!(obj::$m, corr::Heteroscedastic, w::AbstractWeights)
            ψ = influence(obj, w) ; lmul!(Diagonal(w), ψ)
            V = crossprod(ψ)
            corr.corrected && lmul!(varcorrection(corr, w), V)
            _vcov!(obj, V)
        end

        function _vcov!(obj::$m, corr::Clustered, w::UnitWeights)
            ψ = corr.matrix * influence(obj, w)
            V = crossprod(ψ)
            corr.corrected && lmul!(varcorrection(corr, w), V)
            _vcov!(obj, V)
        end

        function _vcov!(obj::$m, corr::Clustered, w::AbstractWeights)
            ψ = corr.matrix * influence(obj, w) ; lmul!(Diagonal(w), ψ)
            V = crossprod(ψ)
            corr.corrected && lmul!(varcorrection(corr, w), V)
            _vcov!(obj, V)
        end

        function _vcov!(obj::$m, corr::CrossCorrelated, w::UnitWeights)
            ψ = influence(obj, w)
            V = crossprod(ψ, corr.matrix)
            corr.corrected && lmul!(varcorrection(corr, w), V)
            _vcov!(obj, V)
        end

        function _vcov!(obj::$m, corr::CrossCorrelated, w::AbstractWeights)
            ψ = influence(obj, w) ; lmul!(Diagonal(w), ψ)
            V = crossprod(ψ, corr.matrix)
            corr.corrected && lmul!(varcorrection(corr, w), V)
            _vcov!(obj, V)
        end
    end
end

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

# GMM WEIGHT MATRIX

wmatrix(obj::GMM, corr::Heteroscedastic, ::UnitWeights) = crossprod(score(obj))
wmatrix(obj::GMM, corr::Clustered, ::UnitWeights) = crossprod(corr.matrix * score(obj))
wmatrix(obj::GMM, corr::CrossCorrelated, ::UnitWeights) = crossprod(score(obj), corr.matrix)
wmatrix(obj::GMM, corr::Heteroscedastic, w::FrequencyWeights) = crossprod(score(obj), w)

function wmatrix(obj::GMM, corr::Heteroscedastic, w::AbstractWeights)
    s = score(obj) ; lmul!(Diagonal(w), s)
    crossprod(s)
end

function wmatrix(obj::GMM, corr::Clustered, w::AbstractWeights)
    s = corr.matrix * score(obj) ; lmul!(Diagonal(w), s)
    crossprod(s)
end

function wmatrix(obj::GMM, corr::CrossCorrelated, w::AbstractWeights)
    s = score(obj) ; lmul!(Diagonal(w), s)
    crossprod(s, corr.matrix)
end

# COVARIANCE MATRIX FOR GMM

function _vcov!(obj::GMM, corr::CorrStructure, w::AbstractWeights)

    j = jacobian(obj, w)
    S = wmatrix(obj, corr, w)

    if obj.method == "Method of moments"
        V = j \ (S / j')
    elseif obj.method == "One-step GMM"
        ψ = j / (j' * j)
        V = crossprod(ψ, S)
    elseif obj.method == "Two-step GMM"
        ω = obj.W \ j
        ψ = ω / (j' * ω)
        V = crossprod(ψ, S)
    elseif obj.method == "Optimal GMM"
        ω = S \ j
        V = inv(j' * ω)
    end

    corr.corrected && lmul!(varcorrection(corr, w), V)
    _vcov!(obj, V)
end

max_stars_repo_name: lbittarello/Microeconometrics

max_issues_repo_name: lbittarello/Microeconometrics

max_forks_repo_name: lbittarello/Microeconometrics


text:
module Manager

include("types.jl")
include("merge-tables.jl")

export MergeTwoTablesTask,
       MergeFunctionStrategyName, NAIVE_JOIN, SIGMOID_JOIN,
       GridMergeParameters
export merge_two_tables

end  # module Manager


text:
using DelayEmbeddings
using StaticArrays
using Test

# Download some test timeseries

# Old:
repo = "https://raw.githubusercontent.com/JuliaDynamics/NonlinearDynamicsTextbook/master/exercise_data"
tsfolder = joinpath(@__DIR__, "timeseries")
todownload = ["$n.csv" for n in 1:4]

mkpath(tsfolder)
for a in todownload
    download(repo*"/"*a, joinpath(tsfolder, a))
end

#New:
todownload = ["test_time_series_lorenz_standard_N_10000_multivariate.csv", "test_time_series_roessler_N_10000_multivariate.csv"]
repo = "https://raw.githubusercontent.com/JuliaDynamics/JuliaDynamics/master/timeseries"
for a in todownload
    download(repo*"/"*a, joinpath(tsfolder, a))
end

ti = time()

diffeq = (atol = 1e-9, rtol = 1e-9, maxiters = typemax(Int))

@testset "DelayEmbeddings tests" begin
    include("dataset_tests.jl")
    include("embedding_tests.jl")
    include("traditional/delaytime_test.jl")
    include("traditional/embedding_dimension_test.jl")
    include("unified/test_pecora.jl")
    include("unified/uzal_cost_test.jl")
    include("unified/test_pecuzal_embedding.jl")

    # The following methods have been tested throughly and also published in research.
    # We know that they work, but unfortunately the tests we have written
    # about them are not good.
    # See https://github.com/JuliaDynamics/DelayEmbeddings.jl/issues/95
    # include("unified/mdop_tests.jl")
    # include("unified/test_garcia.jl")

end

ti = time() - ti
println("\nTest took total time of:")
println(round(ti, sigdigits=3), " seconds or ", round(ti/60, sigdigits=3), " minutes")


text:
using SymPy

function lagrange(order)

  is=collect(-floor(Int,order/2) : 1 : floor(Int,order/2));
  x=symbols("x");
  L=x^0;
  Ls=Array(typeof(L),(order,));

  for j = 1 : order
    ith=is[j];
    for i in is
      if i!= ith
        L=L*(x-i)/(ith-i);
      end
    end
    Ls[j]=copy(L);
    L=x^0;
  end

  return Ls;

end

function fd(order,d;t=Rational)
  Ls=lagrange(order);
  x =symbols("x");

  is=collect(-floor(Int,order/2) : 1 : floor(Int,order/2));
  D=Array(typeof(x),(order,order));
  for i = 1 : order
    for j = 1 : order
      D[i,j] =diff(Ls[j],x,d)(x=>is[i]);
    end
  end

  return Array{t}(D);

end

order=4
D=fd(order+1,2,t=Float64);

max_stars_repo_name: ReidAtcheson/FiniteDifferenceCoeffs

max_issues_repo_name: ReidAtcheson/FiniteDifferenceCoeffs

max_forks_repo_name: ReidAtcheson/FiniteDifferenceCoeffs


text:
#### OUTPUT WRITER ####

"""
`OutputWriter` contains information about requested model outputs.

# Fields
- `format::Symbol` : exported file format. `:VTK`, `:JLD2` or `:MAT`
- `path::String` : path to the folder
- `outputs::Dict{Symbol,Function}` : dictionary mapping `Symbol`s output names to functions of `r` (rheology) and `s` (state)
- `frequency::F` : save output every `frequency` seconds. To disable frequency sampling use `nothing` as field value
- `interval::I` : save output every `interval` iterations. To disable interval sampling use `nothing` as field value
- `data::Dict{Symbol,Vector{Float64}}` : preallocated output data
- `last_output_time::Float64` : last simulation time where an export occured
- `opened_path::Bool` : `True` if the folder has already been created
"""
mutable struct OutputWriter{TF,F,I}
    path::String
    outputs::Dict{Symbol,Function} # functions of r (rheology) and s (state)
    frequency::F
    interval::I
    data::Dict{Symbol,Vector{Float64}} # preallocation of output data
    last_output_time::Float64
    force_path::Bool
end

# Define aliases for differents file formats
VTKOutputWriter{F,I} = OutputWriter{:VTK,F,I}
JLD2OutputWriter{F,I} = OutputWriter{:JLD2,F,I}
MATOutputWriter{F,I} = OutputWriter{:MAT,F,I}

# Multiple output formats :
struct MixedOutputWriter{OF1,F1,I1,OF2,F2,I2}
    ow1::OutputWriter{OF1,F1,I1}
    ow2::OutputWriter{OF2,F2,I2}
end

"""
    OutputWriter(format, model, path, outputs ; frequency = nothing, interval = nothing)

`OutputWriter` constructor. One of the kwargs `frequency` or `interval` has to be given.

# Positional arguments
- `format::Symbol` : exported file format. `:VTK`, `:JLD2` or `:MAT`
- `model::Model`
- `path::String` : path to the folder
- `outputs::Dict{Symbol,Function}` : dictionary mapping `Symbol`s output names to functions of `r` (rheology) and `s` (state)

# Keyword arguments
- `frequency::F` : save output every `frequency` seconds
- `interval::I` : save output every `interval` iterations
- `force_path::Bool = true` : If false and path already exists, a new path variant is created
"""
function OutputWriter(format, model, path, outputs, frequency, interval, force_path)

    ncells = getncells(model.grid)
    data = Dict{Symbol,Vector{Float64}}()
    for (key,value) in pairs(outputs)
        data[key] = Vector{Float64}(undef,ncells)
    end
    if frequency isa Nothing
        @assert interval isa Real
    else
        @assert interval isa Nothing
    end
    TF = typeof(frequency)
    TI = typeof(interval)

    path = create_dir(path, force_path)

    # loading relevant dependencies
    (format == :JLD2) && (@info("loading FileIO package") ; @eval(Rheologies, import FileIO))
    (format == :MAT) && (@info("loading MAT package") ; @eval(Rheologies, import MAT))

    return OutputWriter{format,TF,TI}(path, outputs, frequency, interval, data, model.clock.tspan[1], force_path)
end
OutputWriter(format, model, path, outputs, ::Nothing, ::Nothing, force_path) = @error "Please provide one of the two kwargs `frequency` or `interval`"
OutputWriter(format, model, path, outputs ; frequency = nothing, interval = nothing, force_path = true) = OutputWriter(format, model, path, outputs, frequency, interval, force_path)

VTKOutputWriter(model, path, outputs ; frequency = nothing, interval = nothing, force_path = true) = OutputWriter(:VTK, model, path, outputs, frequency, interval, force_path)
JLD2OutputWriter(model, path, outputs ; frequency = nothing, interval = nothing, force_path = true) = OutputWriter(:JLD2, model, path, outputs, frequency, interval, force_path)
MATOutputWriter(model, path, outputs ; frequency = nothing, interval = nothing, force_path = true) = OutputWriter(:MAT, model, path, outputs, frequency, interval, force_path)

# NOT NEEDED ANYMORE
# function JuAFEM.reinit!(ow::OutputWriter,model)
#     ow.last_output_time = model.clock.tspan[1]
#     ow.opened_path = false
# end
#
# function JuAFEM.reinit!(mow::MixedOutputWriter,model)
#     JuAFEM.reinit!(mow.ow1,model)
#     JuAFEM.reinit!(mow.ow2,model)
# end
### write_output! ###

write_output!(model, u, ::Nothing) = nothing

function write_output!(model, u, ow::OutputWriter{TF,F,I}) where {TF,F,I}
    c = model.clock

    if c.current_time != c.tspan[2]
        if (I <: Real) & (F <: Real)
            c1 = ((model.clock.iter-1)%ow.interval != 0)
            c2 = (model.clock.current_time - ow.last_output_time < ow.frequency)
            !(c1 || c2) && return nothing
        elseif I <: Real
            ((model.clock.iter-1)%ow.interval != 0) && (return nothing)
        elseif F <: Real
            (model.clock.current_time - ow.last_output_time < ow.frequency) && (return nothing)
        end
    end
    # reinitialize data arrays
    for value in values(ow.data)
        fill!(value,0.0)
    end

    # fill data
    nqp = length(model.material_state[1])
    for (el, cell_states) in enumerate(model.material_state)
        r = model.material_properties[el]
        for state in cell_states
            s = state
            for (key,value) in pairs(ow.outputs)
                ow.data[key][el] += value(r,s)/nqp
            end
        end
    end

    # filename
    filename = "iter_"*@sprintf("%04d",model.clock.iter)*"_time_$(model.clock.current_time)"

    # export
    export_sim(filename, model, u, ow)

    # update outputwriter output time
    ow.last_output_time = model.clock.current_time

    return nothing
end

function write_output!(model, u, mixed_ow::MixedOutputWriter)
    write_output!(model, u, mixed_ow.ow1)
    write_output!(model, u, mixed_ow.ow2)
end

# function create_dir!(ow::OutputWriter)
#     path = ow.path
#     if ispath(path)
#         ow.force_path && (return nothing) # we reuse the existing folder
#         if occursin(r"\([0-9]+\)", path[end-4:end])
#             ind1 = findlast('(',path)
#             ind2 = findlast(')',path)
#             num = parse(Int,path[ind1+1:ind2-1])
#             ow.path = path[1:ind1-1]*"("*string(num+1)*")"
#             create_dir!(ow)
#         else
#             ow.path = path*"(1)"
#             create_dir!(ow)
#         end
#     else
#         splitp = splitpath(path)
#         folder = splitp[end]
#         folderdir = joinpath(splitp[1:end-1]...)
#         cd(folderdir)
#         mkdir(folder)
#         ow.opened_path = true
#         return nothing
#     end
# end

function create_dir(path, force_path)
    if ispath(path)
        if force_path
            non_jl_files = filter!(file -> (!occursin(".jl",file) | occursin(".jld",file)),cd(readdir,path))
            if !isempty(non_jl_files)
                @info "directory already contains non .jl files. Do you want to erase them (y/n) ? "
                erase_input(path,non_jl_files)
            end
            return path
        else
            if occursin(r"\([0-9]+\)", path[end-4:end])
                ind1 = findlast('(',path)
                ind2 = findlast(')',path)
                num = parse(Int,path[ind1+1:ind2-1])
                path = path[1:ind1-1]*"("*string(num+1)*")"
                create_dir(path, force_path)
            else
                path = path*"(1)"
                create_dir(path, force_path)
            end
        end
    else
        splitp = splitpath(path)
        folder = splitp[end]
        folderdir = joinpath(splitp[1:end-1]...)
        cd(folderdir)
        mkdir(folder)
        return path
    end
end

function erase_input(path,files)
    key = readline()
    if (key == "y") | (key == "Y")
        @info "All non .jl files removed from directory"
        rm.(joinpath.(Ref(path),files))
    elseif (key == "n") | (key == "N")
        @info "Existing non .jl files left in place"
    else
        @warn "input key is non recognized, please type \"y\" or \"n\""
        erase_input(path,files)
    end
    return nothing
end

function export_sim(filename, model, u, ow::VTKOutputWriter)
    vtk_grid(joinpath(ow.path, filename), model.dofhandler) do vtkfile
        vtk_point_data(vtkfile, model.dofhandler, u) # displacement field
        vtk_point_data(vtkfile, model.dirichlet_bc) # dirichlet boundary conditions
        for (key,value) in pairs(ow.data)
            vtk_cell_data(vtkfile, value, string(key))
        end
    end
    println("VTK saved")
end

function export_sim(filename, model, u, ow::JLD2OutputWriter)
    ow.data[:u] = Vector{eltype(u)}(undef,length(u))
    ow.data[:u] .= u

    # convert keys from symbols to strings
    data_dict = Dict{String,Any}(string(k)=>v  for (k,v) in pairs(ow.data))
    data_dict["t"] = model.clock.current_time # add simulation time
    data_dict["res"] = spatial_resolution(model.grid)
    data_dict["domain_size"] = get_spatial_domain_size(model.grid)

    FileIO.save(joinpath(ow.path, filename*".jld2"), data_dict)
    println("JLD2 saved")
end

function export_sim(filename, model, u, ow::MATOutputWriter)
    ow.data[:u] .= u
    # convert keys from symbols to strings
    data_dict = Dict{String,Any}(string(k)=>v  for (k,v) in pairs(ow.data))
    data_dict["t"] = model.clock.current_time # add simulation time
    data_dict["res"] = [ i for i in spatial_resolution(model.grid)]
    data_dict["domain_size"] = [ i for i in get_spatial_domain_size(model.grid)]
    #save(joinpath(ow.path, filename*".mat"), data_dict)
    println(typeof.(data_dict[key] for key in keys(data_dict)))
    MAT.matwrite(joinpath(ow.path, filename*".mat"), data_dict)
    println("MAT saved")
end

function spatial_resolution(grid::Grid{dim,Cell{2,4,4}}) where {dim}
    lx, ly = get_spatial_domain_size(grid)

    cell = grid.cells[1]
    el_lx = grid.nodes[cell.nodes[2]].x[1] - grid.nodes[cell.nodes[1]].x[1]
    el_ly = grid.nodes[cell.nodes[4]].x[2] - grid.nodes[cell.nodes[1]].x[2]

    return Int(round(lx/el_lx)), Int(round(ly/el_ly))
end

function get_spatial_domain_size(grid)
    x_nodes = [ node.x[1] for node in grid.nodes ]
    y_nodes = [ node.x[2] for node in grid.nodes ]
    lx = maximum(x_nodes) - minimum(x_nodes)
    ly = maximum(y_nodes) - minimum(y_nodes)
    return lx, ly
end


text:
export PERM_NONE,
    PERM_ALL,
    has_permission,
    permissions_in,
    reply,
    split_message,
    plaintext,
    heartbeat_ping,
    upload_file,
    set_game,
    @fetch,
    @fetchval,
    @deferred_fetch,
    @deferred_fetchval

const CRUD_FNS = :create, :retrieve, :update, :delete

"""
Regex expressions for [`split_message`](@ref) to not break Discord formatting.
"""
const STYLES = [
    r"```.+?```"s, r"`.+?`", r"~~.+?~~", r"(_|__).+?\1", r"(\*+).+?\1",
]

"""
Bitwise permission flags.
More details [here](https://discordapp.com/developers/docs/topics/permissions#permissions-bitwise-permission-flags).
"""
@enum Permission begin
    PERM_CREATE_INSTANT_INVITE=1<<0
    PERM_KICK_MEMBERS=1<<1
    PERM_BAN_MEMBERS=1<<2
    PERM_ADMINISTRATOR=1<<3
    PERM_MANAGE_CHANNELS=1<<4
    PERM_MANAGE_GUILD=1<<5
    PERM_ADD_REACTIONS=1<<6
    PERM_VIEW_AUDIT_LOG=1<<7
    PERM_VIEW_CHANNEL=1<<10
    PERM_SEND_MESSAGES=1<<11
    PERM_SEND_TTS_MESSAGES=1<<12
    PERM_MANAGE_MESSAGES=1<<13
    PERM_EMBED_LINKS=1<<14
    PERM_ATTACH_FILES=1<<15
    PERM_READ_MESSAGE_HISTORY=1<<16
    PERM_MENTION_EVERYONE=1<<17
    PERM_USE_EXTERNAL_EMOJIS=1<<18
    PERM_CONNECT=1<<20
    PERM_SPEAK=1<<21
    PERM_MUTE_MEMBERS=1<<22
    PERM_DEAFEN_MEMBERS=1<<23
    PERM_MOVE_MEMBERS=1<<24
    PERM_USE_VAD=1<<25
    PERM_PRIORITY_SPEAKER=1<<8
    PERM_STREAM=1<<9
    PERM_CHANGE_NICKNAME=1<<26
    PERM_MANAGE_NICKNAMES=1<<27
    PERM_MANAGE_ROLES=1<<28
    PERM_MANAGE_WEBHOOKS=1<<29
    PERM_MANAGE_EMOJIS=1<<30
end
@boilerplate Permission :export

const PERM_NONE = 0
const PERM_ALL = |(Int.(instances(Permission))...)

"""
    has_permission(perms::Integer, perm::Permission) -> Bool

Determine whether a bitwise OR of permissions contains one [`Permission`](@ref).

## Examples
```jldoctest; setup=:(using Discord)
julia> has_permission(0x0420, PERM_VIEW_CHANNEL)
true

julia> has_permission(0x0420, PERM_ADMINISTRATOR)
false

julia> has_permission(0x0008, PERM_MANAGE_ROLES)
true
```
"""
function has_permission(perms::Integer, perm::Permission)
    admin = perms & Int64(PERM_ADMINISTRATOR) == Int64(PERM_ADMINISTRATOR)
    has = perms & Int64(perm) == Int64(perm)
    return admin || has
end

"""
    permissions_in(m::Member, g::Guild, ch::DiscordChannel) -> Int64

Compute a [`Member`](@ref)'s [`Permission`](@ref)s in a [`DiscordChannel`](@ref).
"""
function permissions_in(m::Member, g::Guild, ch::DiscordChannel)
    !ismissing(m.user) && m.user.id == g.owner_id && return PERM_ALL

    # Get permissions for @everyone.
    idx = findfirst(r -> r.name == "@everyone", g.roles)
    everyone = idx === nothing ? nothing : g.roles[idx]
    perms = idx === nothing ? Int64(0) : everyone.permissions
    perms & Int64(PERM_ADMINISTRATOR) == Int64(PERM_ADMINISTRATOR) && return PERM_ALL

    roles = idx === nothing ? m.roles : [everyone.id; m.roles]

    # Apply role overwrites.
    for role in roles
        idx = findfirst(
            o -> o.type === OT_ROLE && o.id == role,
            coalesce(ch.permission_overwrites, Overwrite[]),
        )
        if idx !== nothing
            o = ch.permission_overwrites[idx]
            perms &= ~o.deny
            perms |= o.allow
        end
    end

    # Apply user-specific overwrite.
    if !ismissing(m.user)
        idx = findfirst(
            o -> o.type === OT_MEMBER && o.id == m.user.id,
            coalesce(ch.permission_overwrites, Overwrite[]),
        )
        if idx !== nothing
            o = ch.permission_overwrites[idx]
            perms &= ~o.deny
            perms |= o.allow
        end
    end

    return perms
end

Base.print(io::IO, c::DiscordChannel) = print(io, "<#$(c.id)>")
Base.print(io::IO, r::Role) = print(io, "<@&$(r.id)>")
Base.print(io::IO, u::User) = print(io, "<@$(u.id)>")
function Base.print(io::IO, m::Member)
    if ismissing(m.user)
        print(io, something(coalesce(m.nick, "<unknown member>"), "<unknown member>"))
    elseif ismissing(m.nick) || m.nick === nothing
        print(io, m.user)
    else
        print(io, "<@!$(m.user.id)>")
    end
end
function Base.print(io::IO, e::Emoji)
    s = if e.id === nothing
        coalesce(e.require_colons, false) ? ":$(e.name):" : e.name
    else
        coalesce(e.animated, false) ? "<a:$(e.name):$(e.id)>" : "<:$(e.name):$(e.id)>"
    end
    print(io, s)
end

"""
    reply(
        c::Client,
        m::Message,
        content::Union{AbstractString, AbstractDict, NamedTuple, Embed};
        at::Bool=false,
    ) -> Future{Response}

Reply (send a message to the same [`DiscordChannel`](@ref)) to a [`Message`](@ref).
If `at` is set, then the message is prefixed with the sender's mention.
"""
function reply(c::Client, m::Message, content::AbstractString; at::Bool=false)
    at && !ismissing(m.author) && (content = string(m.author, " ", content))
    return create_message(c, m.channel_id; content=content)
end

function reply(
    c::Client,
    m::Message,
    embed::Union{AbstractDict, NamedTuple, Embed};
    at::Bool=false,
)
    return if at && !ismissing(m.author)
        create_message(c, m.channel_id; content=string(m.author), embed=embed)
    else
        create_message(c, m.channel_id; embed=embed)
    end
end

"""
    filter_ranges(u::Vector{UnitRange{Int}})

Filter a list of ranges, discarding ranges included in other ranges from the list.

# Example
```jldoctest; setup=:(using Discord)
julia> Discord.filter_ranges([1:5, 3:8, 1:20, 2:16, 10:70, 25:60, 5:35, 50:90, 10:70])
4-element Vector{UnitRange{Int64}}:
 1:20
 5:35
 50:90
 10:70
```
"""
function filter_ranges(u::Vector{UnitRange{Int}})
    v = fill(true, length(u))
    for i in 1:length(u)
        if !all(m -> (m[1] == i) || (u[i] ⊈ m[2]),
                m for m in enumerate(u) if v[m[1]] == true)
            v[i] = false
        end
    end
    return u[v]
end

"""
    split_message(text::AbstractString; chunk_limit::UInt=2000,
                  extrastyles::Vector{Regex}=Vector{Regex}(),
                  forcesplit::Bool = true) -> Vector{String}

Split a message into chunks with at most chunk_limit length, preserving formatting.

The `chunk_limit` has as default the 2000 character limit of Discord's messages,
but can be changed to any nonnegative integer.

Formatting is specified by [`STYLES`](@ref)) and can be aggregated
with the `extrastyles` argument.

Discord limits messages to 2000, so the code forces split if format breaking
cannot be avoided. If desired, however, this behavior can be lifter by setting
`forcesplit` to false.

## Examples
```jldoctest; setup=:(using Discord)
julia> split_message("foo")
1-element Vector{String}:
 "foo"

julia> split_message(repeat('.', 1995) * "**hello, world**")[2]
"**hello, world**"

julia> split_message("**hello**, *world*", chunk_limit=10)
2-element Vector{String}:
 "**hello**,"
 "*world*"

julia> split_message("**hello**, _*beautiful* world_", chunk_limit=15)
┌ Warning: message was forced-split to fit the desired chunk length limit 15
└ @ Main REPL[66]:28
3-element Vector{String}:
 "**hello**,"
 "_*beautiful* wo"
 "rld_"

julia> split_message("**hello**, _*beautiful* world_", chunk_limit=15, forcesplit=false)
┌ Warning: message could not be split into chunks smaller than the length limit 15
└ @ Main REPL[66]:32
2-element Vector{String}:
 "**hello**,"
 "_*beautiful* world_"

julia> split_message("**hello**\n=====\n", 12)
2-element Vector{String}:
 "**hello**\n=="
 "==="

julia> split_message("**hello**\n≡≡≡≡≡\n", chunk_limit=12, extrastyles = [r"\n≡+\n"])
2-element Vector{String}:
 "**hello**"
 "≡≡≡≡≡"
"""
function split_message(text::AbstractString; chunk_limit::Int=2000,
                       extrastyles::Vector{Regex}=Vector{Regex}(),
                       forcesplit::Bool = true)
    chunks = String[]

    while !isempty(text)
        if length(text) ≤ chunk_limit
            push!(chunks, strip(text))
            return chunks
        end
        # get ranges associated with the formattings
        # mranges = vcat(findall.(union(STYLES, extrastyles),Ref(text))...) can't use findall in julia 1.0 and 1.1 ...
        mranges = [m.offset:m.offset+ncodeunits(m.match)-1 for m in vcat(collect.(eachmatch.(union(STYLES, extrastyles), text))...)]
        # filter ranges to eliminate inner formattings
        franges = filter_ranges(mranges)
        # get ranges that get split apart by the chunk limit - there should be only one, unless text is ill-formatted
        splitranges = filter(r -> (length(text[1:r[1]]) ≤ chunk_limit) & (length(text[1:r[end]]) > chunk_limit), franges)

        if length(splitranges) > 0
            stop = minimum(map(r -> prevind(text, r[1]), splitranges))
        end

        if length(splitranges) == 0
            # get highest valid unicode index if no range is split apart
            stop = maximum(filter(n -> length(text[1:n])≤chunk_limit, thisind.(Ref(text), 1:ncodeunits(text))))
        elseif (stop == 0) && (forcesplit == true)
            # get highest valid unicode if format breaking cannot be avoided and forcesplit is true
            stop = maximum(filter(n -> length(text[1:n])≤chunk_limit, thisind.(Ref(text), 1:ncodeunits(text))))
            # @warn "message was forced-split to fit the desired chunk length limit $chunk_limit"
        elseif stop == 0
            # give up at this point if current chunk cannot be split and `forcesplit` is set to false
            push!(chunks, strip(text))
            # @warn "message could not be split into chunks smaller than the length limit $chunk_limit"
            return chunks
        end

        # push chunk and select remaining text
        push!(chunks, strip(text[1:stop]))
        text = strip(text[nextind(text, stop):end])
    end

    return chunks
end

"""
    plaintext(m::Message) -> String
    plaintext(c::Client, m::Message) -> String

Get the [`Message`](@ref) contents with any [`User`](@ref) mentions replaced with their
plaintext. If a [`Client`](@ref) is provided, [`DiscordChannel`](@ref)s [`Role`](@ref) are
also replaced. However, only channels and roles stored in state are replaced; no API
requests are made.
"""
function plaintext(m::Message)
    content = m.content

    for u in coalesce(m.mentions, User[])
        name = "@$(u.username)"
        content = replace(content, "<@$(u.id)>" => name)
        content = replace(content, "<@!$(u.id)>" => name)
    end

    return content
end

function plaintext(c::Client, m::Message)
    content = m.content

    for u in coalesce(m.mentions, User[])
        member = get(c.state, Member; guild=m.guild_id, user=u.id)
        nick = if member !== nothing && member.nick isa String
            "@$(member.nick)"
        else
            "@$(u.username)"
        end
        content = replace(content, "<@$(u.id)>" => "@$(u.username)")
        content = replace(content, "<@!$(u.id)>" => "@$nick")
    end

    guild = get(c.state, Guild; guild=m.guild_id)
    if guild !== nothing
        for r in coalesce(m.mention_roles, Snowflake[])
            role = get(c.state, Role; guild=m.guild_id, role=r)
            if role !== nothing
                content = replace(content, "<@&$r>" => "@$(role.name)")
            end
        end

        for cap in unique(eachmatch(r"<#(\d+?)>", content))
            ch = get(c.state, DiscordChannel; channel=parse(Snowflake, first(cap.captures)))
            if ch !== nothing
                content = replace(content, cap.match => "#$(ch.name)")
            end
        end
    end

    return content
end

"""
    heartbeat_ping(c::Client) -> Nullable{Period}

Get the [`Client`](@ref)'s ping time to the gateway. If the client is not connected, or no
heartbeats have been sent/acknowledged, `nothing` is returned.
"""
function heartbeat_ping(c::Client)
    isopen(c) || return nothing
    zero = DateTime(0)
    return c.last_hb == zero || c.last_ack == zero ? nothing : c.last_ack - c.last_hb
end

"""
    upload_file(c::Client, ch::DiscordChannel, path::AbstractString; kwargs...) -> Message

Send a [`Message`](@ref) with a file [`Attachment`](@ref). Any keywords are passed on to
[`create_message`](@ref).
"""
function upload_file(c::Client, ch::DiscordChannel, path::AbstractString; kwargs...)
    return create_message(c, ch.id; kwargs..., file=open(path))
end

"""
    set_game(
        c::Client,
        game::AbstractString;
        type::Union{ActivityType, Int}=AT_GAME,
        since::Nullable{Int}=c.presence["since"],
        status::Union{PresenceStatus, AbstractString}=c.presence["status"],
        afk::Bool=c.presence["afk"],
        kwargs...,
    ) -> Bool

Shortcut for [`update_status`](@ref) to set the [`Client`](@ref)'s [`Activity`](@ref). Any
additional keywords are passed into the `activity` section.
"""
function set_game(
    c::Client,
    game::AbstractString;
    type::Union{ActivityType, Int}=AT_GAME,
    since::Nullable{Int}=c.presence["since"],
    status::Union{PresenceStatus, AbstractString}=c.presence["status"],
    afk::Bool=c.presence["afk"],
    kwargs...,
)
    activity = merge(Dict("name" => game, "type" => type), kwargs)
    return update_status(c, since, activity, status, afk)
end

"""
    @fetch [functions...] block

Wrap all calls to the specified CRUD functions ([`create`](@ref), [`retrieve`](@ref),
[`update`](@ref), and [`delete`](@ref)) with `fetch` inside a block. If no functions are
specified, all CRUD functions are wrapped.

## Examples
Wrapping all CRUD functions:
```julia
@fetch begin
    guild_resp = create(c, Guild; name="foo")
    guild_resp.ok || error("Request for new guild failed")
    channel_resp = retrieve(c, DiscordChannel, guild_resp.val)
end
```
Wrapping only calls to `retrieve`:
```julia
@fetch retrieve begin
    resp = retrieve(c, DiscordChannel, 123)
    future = create(c, Message, resp.val; content="foo")  # Behaves normally.
end
```
"""
macro fetch(exs...)
    validate_fetch(exs...)
    fns = length(exs) == 1 ? CRUD_FNS : exs[1:end-1]
    ex = wrapfn!(exs[end], fns, :fetch)
    quote
        $ex
    end
end

"""
    @fetchval [functions...] block

Identical to [`@fetch`](@ref), but calls are wrapped with [`fetchval`](@ref) instead.
"""
macro fetchval(exs...)
    validate_fetch(exs...)
    fns = length(exs) == 1 ? CRUD_FNS : exs[1:end-1]
    ex = wrapfn!(exs[end], fns, :fetchval)
    quote
        $ex
    end
end

"""
    @deferred_fetch [functions...] block

Identical to [`@fetch`](@ref), but `Future`s are not `fetch`ed until the **end** of the
block. This is more efficient, but only works when there are no data dependencies in the
block.

## Examples
This will work:
```julia
@deferred_fetch begin
    guild_resp = create(c, Guild; name="foo")
    channel_resp = retrieve(c, DiscordChannel, 123)
end
```
This will not, because the second call is dependent on the first value:
```julia
@deferred_fetch begin
    guild_resp = create(c, Guild; name="foo")
    channels_resp = retrieve(c, DiscordChannel, guild_resp.val)
end
```
"""
macro deferred_fetch(exs...)
    validate_fetch(exs...)
    fns = length(exs) == 1 ? CRUD_FNS : exs[1:end-1]
    ex = deferfn!(exs[end], fns, :fetch)
    quote
        $ex
    end
end

"""
    @deferred_fetchval [functions...] block

Identical to [`@deferred_fetch`](@ref), but `Future`s have [`fetchval`](@ref) called on
them instead of `fetch`.
"""
macro deferred_fetchval(exs...)
    validate_fetch(exs...)
    fns = length(exs) == 1 ? CRUD_FNS : exs[1:end-1]
    ex = deferfn!(exs[end], fns, :fetchval)
    quote
        $ex
    end
end

# Validate the arguments to CRUD macros.
function validate_fetch(exs...)
    if !(exs[end] isa Expr && exs[end].head === :block)
        throw(ArgumentError("Final argument must be a block"))
    end
    if !all(fn -> fn in CRUD_FNS, exs[1:end-1])
        throw(ArgumentError("Only CRUD functions can be wrapped"))
    end
end

# Wrap calls to certain functions in a call to another function.
wrapfn!(ex, ::Tuple, ::Symbol) = esc(ex)
function wrapfn!(ex::Expr, fns::Tuple, with::Symbol)
    if ex.head === :call && ex.args[1] in fns
        ex = :($(esc(with))($(esc(ex))))
    else
        map!(arg -> wrapfn!(arg, fns, with), ex.args, ex.args)
    end
    return ex
end

# Defer fetching a Future until the end of a block.
deferfn!(ex, ::Tuple) = (esc(ex), Pair{Symbol, Symbol}[])
function deferfn!(ex::Expr, fns::Tuple)
    renames = Pair{Symbol, Symbol}[]

    if ex.head === :(=) && ex.args[2] isa  Expr && ex.args[2].args[1] in fns
        newsym = gensym(ex.args[1])
        push!(renames, ex.args[1] =>  newsym)
        ex.args[1] = newsym
        map!(esc, ex.args, ex.args)
    else
        for i in eachindex(ex.args)
            ex.args[i], rs = deferfn!(ex.args[i], fns)
            append!(renames, rs)
        end
    end

    return ex, renames
end
function deferfn!(ex, fns::Tuple, deferred::Symbol)
    ex, renames = deferfn!(ex, fns)
    repls = map(r -> :($(esc(r[1])) = $(esc(deferred))($(esc(r[2])))), renames)
    append!(ex.args, repls)
    return ex
end


text:
mutable struct Tokenizer
    str::String
    i::Int
    offset::Int
    tokens::Array{Any,1}
    Tokenizer(str) = new(str, 1, 0, [])
end

# ==============================================================================
isDigit(l) =   '0' ≤ l ≤ '9'
isHex(l) = ('a' ≤ l ≤ 'f') || ('A' ≤ l ≤ 'F') || isDigit(l)
isUpper(l) = 'A' ≤ l ≤ 'Z'
isLower(l) = 'a' ≤ l ≤ 'z'
isAlpha(l) = isLower(l) || isUpper(l) || l == '-'
# ==============================================================================
# some convenience functions
# ==============================================================================
is(p, symbol) = (p.i < length(p.str)) &&  (p.str[p.i] == symbol)
at(p) = p.str[p.i + p.offset]
range(p) = p.str[p.i:p.i + (p.offset-1)] # Why -1 ???
getIt(p, symbol)= ( push!(p.tokens, symbol); p.i+=1)
# ==============================================================================
function parseHex(v) # Parse hex as RGB[A] values
    if length(v) == 3
        color = [parse(Int, "0x" * v[1:1]), parse(Int, "0x" * v[2:2]), parse(Int, "0x" * v[3:3])]
    elseif length(value) == 6
        color = [parse(Int, "0x" * v[1:2]), parse(Int, "0x" * v[3:4]), parse(Int, "0x" * v[5:6])]
    elseif length(value) == 8
        color = [parse(Int, "0x" * v[1:2]), parse(Int, "0x" * v[3:4]), parse(Int, "0x" * v[5:6]), parse(Int, "0x" * v[7:8])]
    else
        println("Error parsing Hex to color values")
    end
    color
end
# ==============================================================================
function getHex(p)
  p.i+=1
  p.offset = 1 # Assume 0 is true
  while isHex(at(p))
      p.i + p.offset < length(p.str) ?  p.offset+=1 : break
  end
   push!(p.tokens, Hex(parseHex(range(p))));
   p.i += p.offset
end
# ==============================================================================
function getAlpha(p)
  p.offset = 1 # Assume 0 is true
  while isAlpha(at(p)) || isDigit(at(p))
      p.i + p.offset < length(p.str) ?  p.offset+=1 : break
  end
   push!(p.tokens, Alpha(range(p)));
   p.i += p.offset
end
# ==============================================================================
function getDigit(p)
  p.offset = 1 # Assume 0 is true
  while isDigit(at(p))
    p.i + p.offset < length(p.str) ?  p.offset+=1 : break
  end
  push!(p.tokens, Digit(range(p)));
  p.i += p.offset
end
# ==============================================================================
#  Dot Digit Open Close Value End Unit Dotdot Hex At
#  Hash Dollar Reference Alpha Var Attr Tag Class
function getString(p, symbol)
    p.i+=1 # Trim off the quotes.
    p.offset = 1
    while at(p) != symbol
        p.i + p.offset < length(p.str) ? p.offset+=1 : break
    end #
    push!(p.tokens, Str(range(p)));
    p.offset+=1
    p.i += p.offset
end
# ==============================================================================
# Main function to tokenize
# ==============================================================================
function tokenize(s)
  t = Tokenizer(s)

while t.i < length(t.str) # lastindex

    (is(t, ' ') || is(t, '\t')) && getIt(t, Space(t.str[t.i:t.i]))
    (is(t, ',') || is(t, ';') || is(t, '\r') || is(t, '\n')) && getIt(t, End(t.str[t.i:t.i]))

    isDigit(t.str[t.i]) && getDigit(t)
    isAlpha(t.str[t.i]) && getAlpha(t)
    is(t, '#')          && getHex(t)
    is(t, '\'')         && getString(t, '\'')
    is(t, '\"')         && getString(t, '\"')
    is(t, '{')          && getIt(t, Open("{"))
    is(t, '}')          && getIt(t, Close("}"))
    is(t, ':')          && getIt(t, Dotdot(":"))
    is(t, '%')          && getIt(t, Unit("%"))
    is(t, '$')          && getIt(t, Dollar("\$"))
    is(t, '@')          && getIt(t, At("@"))
    is(t, '.')          && getIt(t, Dot("."))
end
  return t
end


text:
num(x::Array{<:Rational}) = map(y->y.num, x)
den(x::Array{<:Rational}) = map(y->y.den, x)

# Transform a rational matrix into an integer matrix by multiplying each row
function intmat(A::Matrix{<:Rational})
    func = z -> round.(Int, z*(lcm(den(z))//gcd(num(z))))
    mapslices(func, A; dims=2)
end

function collinear(a::Array{Int}, b::Array{Int})
    length(a) == length(b) || error("Input arrays must have the same dimenisions")
    a_n, b_n = a/gcd(a), b/gcd(b)
    if a_n==b_n || a_n==-b_n; return true; end
    return false
end

# Normalises the columns of a matrix
function norm_cols(A::AbstractMatrix{<:Integer})
    norm_A = copy(A)
    for i in 1:size(A,2)
        if (g = gcd(A[:,i])) > 1; norm_A[:, i] = div.(norm_A[:, i], g); end
    end
    return norm_A
end

# Check matrices are the same up to permuting columns
function check_columns_same(A::Matrix{T}, B::Matrix{T}) where T<:Real
    m,n  = size(A)
    if size(B) != (m,n) return false end
    sA = Set{Array{T, 1}}()
    sB = Set{Array{T, 1}}()
    for i in 1:n
        push!(sA, A[:, i])
        push!(sB, B[:, i])
    end
    return sA == sB
end

# Check columns of two matrices span same linear subspace
# It is assumed columns in each matrix are linearly independent
function check_column_span_same(A::Matrix{T}, B::Matrix{T}) where T<:Real
    m, n = size(A)
    if size(B) != (m,n) return false end
    rankA = rank(A)
    return rank(A) == rank([A B])
end

# Tests whether all elements of vector are zero
function is_zero(x::Vector{T}) where T<:Real
    z = zero(T)
    for xi in x
        if xi != z return false end
    end
    return true
end

# Verify columns of input matrix are orthogonal (no zero columns allowed)
function check_columns_orthogonal(A::Matrix{T}) where T<:Real
    m, n = size(A)
    for i in 1:n
        if is_zero(A[:,i]) return false end
        for j in 1:(i-1)
            if dot(A[:,i], A[:,j]) != zero(T) return false end
        end
    end
    return true
end

# Use a modified version of Gram-Schmitd to orthonalize columns of integer matrix
function gram_schmidt(A::Matrix{T}) where T<:Integer
    m, n = size(A)
    B = Array{T}(undef, m, n)
    for i in 1:n
        sq_norms = T[dot(B[:,j], B[:,j]) for j in 1:(i-1)]
        prod_sq_norms = prod(sq_norms)
        B[:,i] = prod_sq_norms * A[:,i]
        for j in 1:(i-1)
            B[:,i] -= div(prod_sq_norms,sq_norms[j]) * dot(A[:,i], B[:,j]) * B[:,j]
        end
        g = gcd(B[:,i])
        if g > 1
            B[:, i] = div.(B[:, i], g)
        end
    end
    return B
end

# Transforms rays so that they are orthogonal to all elements in a linear subspace
# Arguments:
# `A::Matrix{T}`: matrix of rays (each column is a ray)
# `B::Matrix{T}`: matrix whose columns form a basis of linear subspace
function normalise_rays(A::Matrix{T}, B::Matrix{T}) where T<:Integer
    m, n = size(A)
    k = size(B,2) # Size of basis
    new_rays = Array{T}(undef, m, n)
    B = gram_schmidt(B) # must have orthogonal basis for linear subspace
    sq_norms = T[dot(B[:,j], B[:,j]) for j in 1:k]
    prod_sq_norms = prod(sq_norms)
    for i in 1:n
        new_rays[:,i] = prod_sq_norms * A[:,i]
        for j in 1:k
            new_rays[:,i] -= div(prod_sq_norms,sq_norms[j]) * dot(A[:,i], B[:,j]) * B[:,j]
            g = gcd(new_rays[:,i])
            if g > 1;
                new_rays[:, i] = div.(new_rays[:, i], g)
            end
        end
    end
    return new_rays
end


text:
struct RainGauge
    name::String
    id::Int
    lat::Float64
    lon::Float64
    steps::Array{String}
    precip::Array{Float64}
end

in_path = ARGS[1]
out_path = ARGS[2]

file_out = open(out_path, "w")
data_files = readdir(in_path)

gauges = Dict()

for (index, file_name) in enumerate(data_files)
    file_path = in_path*"/"*file_name
    file = open(file_path, "r")
    file_content = read(file, String)
    data = eachmatch(r"(\d{2}/\d{2}/\d{4},\d{2}:\d{2}:\d{2}),(\d{1,4}.\d{1,4}),150", file_content)
    gauge_content = [parse(Float64, x.captures[2]) for x in data]
    timesteps = [replace(x.captures[1], "," => "-") for x in data]
    timesteps = [x[7:10]*"-"*x[4:5]*"-"*x[1:2]*" "*x[12:end] for x in timesteps]
    station_name = match(r"Station Name, (.*)\r", file_content).captures[1]
    station_num = parse(Int, match(r"Station Number, (.*)\r", file_content).captures[1])
    latitude = parse(Float64, match(r"Latitude,(-?\d{1,4}.{1,9})", file_content).captures[1])
    longitude = parse(Float64, match(r"Longitude,(-?\d{1,4}.{1,9})", file_content).captures[1])
    gauges[index] = RainGauge(station_name, station_num, latitude, longitude, timesteps, gauge_content)
end

write(file_out, "name,id,latitude,longitude,"*join(gauges[1].steps, ",")*"\n")
for gauge in gauges
    gauge = gauge[2]
    write(file_out, gauge.name*","*string(gauge.id)*","*string(gauge.lat)*","*string(gauge.lon)*","*join(gauge.precip, ",")*"\n")
end

max_stars_repo_name: smart-facility/smart-faware

max_issues_repo_name: smart-facility/smart-faware

max_forks_repo_name: smart-facility/smart-faware


text:
# This file was generated, do not modify it. # hide
@load LogisticClassifier pkg=MLJLinearModels
X2 = select(X, Not([:Year, :Today]))
clf = machine(LogisticClassifier(), X2, y)


text:
immutable mbedtls_mpi
    s::Cint
    n::Csize_t
    p::Ptr{Cuint}
end

immutable mbedtls_rsa_context
    ver::Cint
    len::Csize_t
    N::mbedtls_mpi
    E::mbedtls_mpi
    D::mbedtls_mpi
    P::mbedtls_mpi
    Q::mbedtls_mpi
    # More fields follow, but omitted here, since they
    # are not required for this wrapper
end

type RSA
    data::Ptr{mbedtls_rsa_context}

    function RSA(padding=MBEDTLS_RSA_PKCS_V21, hash_id=MD_MD5)
        ctx = new()
        ctx.data = Libc.malloc(1000)
        ccall((:mbedtls_rsa_init, MBED_CRYPTO), Void,
            (Ptr{Void}, Cint, Cint),
            ctx.data, padding, hash_id)
        finalizer(ctx, ctx->begin
            ccall((:mbedtls_rsa_free, MBED_CRYPTO), Void, (Ptr{Void},), ctx.data)
            Libc.free(ctx.data)
        end)
        ctx
    end
end

function mpi_import!(mpi::Ptr{mbedtls_mpi}, b::BigInt)
    # Export from GMP
    size = ndigits(b, 2)
    nbytes = div(size+8-1,8)
    data = Array(UInt8, nbytes)
    count = Ref{Csize_t}(0)
    ccall((:__gmpz_export,:libgmp), Ptr{Void},
            (Ptr{Void}, Ptr{Csize_t}, Cint, Csize_t, Cint, Csize_t, Ptr{BigInt}),
            data, count, 1, 1, 1, 0, &b)
    @assert count[] == nbytes
    # Import into mbedtls
    @err_check ccall((:mbedtls_mpi_read_binary, MBED_CRYPTO), Cint,
        (Ptr{mbedtls_mpi}, Ptr{UInt8}, Csize_t),
        mpi, data, nbytes)
end

function mpi_size(mpi::Ptr{mbedtls_mpi})
    ccall((:mbedtls_mpi_size, MBED_CRYPTO), Csize_t, (Ptr{mbedtls_mpi},), mpi)
end

# 0.4 had a buggy implementation of fieldoffset (off by one)
if !isdefined(Base, :fieldoffset) && isdefined(Base, :field_offset)
    fieldoffset(T, i) = Base.field_offset(T, i - 1)
end

function pubkey_from_vals!(ctx::RSA, e::BigInt, n::BigInt)
    Nptr = Ptr{mbedtls_mpi}(ctx.data+fieldoffset(mbedtls_rsa_context,3 #= :N =#))
    mpi_import!(Nptr, n)
    mpi_import!(Ptr{mbedtls_mpi}(ctx.data+fieldoffset(mbedtls_rsa_context,4 #= :E =#)), e)
    nptr_size = mpi_size(Nptr)
    unsafe_store!(Ptr{Csize_t}(ctx.data+fieldoffset(mbedtls_rsa_context,2 #=:len =#)),
        nptr_size)
    @err_check ccall((:mbedtls_rsa_check_pubkey, MBED_CRYPTO), Cint,
        (Ptr{Void},), ctx.data)
    ctx
end

function verify(ctx::RSA, hash_alg::MDKind, hash, signature, rng = nothing; using_public=true)
    (!using_public && rng == nothing) &&
        error("Private key verification requires the rng")
    # All errors, including validation errors throw
    @err_check ccall((:mbedtls_rsa_pkcs1_verify, MBED_CRYPTO), Cint,
        (Ptr{Void}, Ptr{Void}, Ptr{Void}, Cint, Cint, Csize_t, Ptr{UInt8}, Ptr{UInt8}),
        ctx.data,
        rng == nothing ? C_NULL : c_rng,
        rng == nothing ? C_NULL : pointer_from_objref(rng),
        using_public ? 0 : 1,
        hash_alg, sizeof(hash), hash, signature)
end

function gen_key!(ctx::RSA, f_rng, p_rng, nbits, exponent)
    @err_check ccall((:mbedtls_rsa_gen_key, MBED_CRYPTO), Cint,
        (Ptr{Void}, Ptr{Void}, Ptr{Void}, Cint, Cint),
        ctx.data, f_rng, p_rng, nbits, exponent)
    ctx
end

function gen_key(rng::AbstractRNG, nbits=2048, exponent=65537)
    ctx = RSA()
    gen_key!(ctx, c_rng, pointer_from_objref(rng), nbits, exponent)
    ctx
end

function public(ctx::RSA, input, output)
    @err_check ccall((:mbedtls_rsa_public, MBED_CRYPTO), Cint,
        (Ptr{Void}, Ptr{Void}, Ptr{Void}), ctx.data, input, output)
    output
end

function private(ctx::RSA, f_rng, p_rng, input, output)
    @err_check ccall((:mbedtls_rsa_private, MBED_CRYPTO), Cint,
        (Ptr{Void}, Ptr{Void}, Ptr{Void}, Ptr{Void}, Ptr{Void}),
        ctx.data, f_rng, p_rng, input, output)
    output
end

function private(ctx::RSA, rng::AbstractRNG, input, output)
    private(ctx, c_rng, pointer_from_objref(rng), input, output)
end


text:
using Mimi
using MimiDICE2010
using Distributions

"""
Case: We want to do an SCC calculation with MimiDICE2010, which consists of running both a base and modified model
(the latter being a model including an additional emissions pulse,
see the create_marginal_model function or create your own two models).
We then take the difference between the consumption level in these two models and
obtain the discounted net present value to get the SCC.
"""

## STEP 1 Defining typical variables for a simulation, number of trials and simulation definition

# Define the number of trials
N = 100

# define your simulation (defaults to Monte Carlo sampling)
sd = @defsim begin
    t2xco2 = Truncated(Gamma(6.47815626,0.547629469), 1.0, Inf) # a dummy distribution
end

## STEP 2 Payload object

"""
Simulation definitions can hold a user-defined payload object which is not used or modified by Mimi.
In this example, we will use the payload to hold an array of pre-computed
discount factors that we will use in the SCC calculation, as well as a storage array for saving the SCC values.
"""
# Choose what year to calculate the SCC for
scc_year = 2015
year_idx = findfirst(isequal(scc_year), MimiDICE2010.model_years)

# Pre-coimpute the discount factors for each discount rate
discount_rates = [0.03, 0.05, 0.07]
nyears = length(MimiDICE2010.model_years)
discount_factors = [[zeros(year_idx - 1)... [(1/(1 + r))^((t-year_idx)*10) for t in year_idx:nyears]...] for r in discount_rates]

# Create an array to store the compute SCC in each trial for each discount rate
scc_results = zeros(N, length(discount_rates))

# Set the payload object in the simulation definition
my_payload_object = (discount_factors, scc_results) # In this case, the payload object is a tuple which holds both both arrays
Mimi.set_payload!(sd, my_payload_object)

## STEP 3 Post-trial function
"""
User may ned to perform other calculations before or after each trial is run
For example, the SCC is calculated using two models, so this calculation needs to happen in a post-trial function
"""
function my_scc_calculation(sim_inst::SimulationInstance, trialnum::Int, ntimesteps::Int, tup::Nothing)
    mm = sim_inst.models[1]
    discount_factors, scc_results = Mimi.payload(sim_inst)  # Unpack the payload object

    marginal_damages = mm[:neteconomy, :C] * -1 * 10^12 * 12/44 # convert from trillion $/ton C to $/ton CO2; multiply by -1 to get positive value for damages
    for (i, df) in enumerate(discount_factors)
        scc_results[trialnum, i] = sum(df .* marginal_damages .* 10)
    end
end

# STEP 4: RUN THE SIMULATION
# Build the marginal model
mm = MimiDICE2010.get_marginal_model(year = scc_year)   # The additional emissions pulse will be added in the specified year

# Run
si = run(sd, mm, N; trials_output_filename = "ecs_sample.csv", post_trial_func = my_scc_calculation)

# View the scc_results by retrieving them from the payload object
scc_results = Mimi.payload(si)[2]   # Recall that the SCC array was the second of two arrays we stored in the payload tuple

max_stars_repo_name: tlorans/SustainableFinance

max_issues_repo_name: tlorans/SustainableFinance

max_forks_repo_name: tlorans/SustainableFinance


text:
@testset "acoustic sources" begin
    a2 = Acoustic(0.1,0.1 + 0.0im,2)
    a2_host = Acoustic(1.3,1.5 + 0.0im,2)

    # Create two point sources
    source_position = [0.0,1.0]
    amplitude = 1.0
    s1 = point_source(a2, source_position, amplitude)
    s2 = point_source(a2, 2.0*source_position, amplitude)

    # Create new souce as a linear combination of two other sources
    s3 = 2*s1 + s2;

    # Check that the field is indeed a linear conbination
    x = [1.0, 1.0]
    @test field(s3,x,1.0) == 2*field(s1,x,1.0) + field(s2,x,1.0)

    # Test regular basis expansion of the source
    ω = 0.8
    basis_order = 7
    centre =  [1.0,0.0]

    s3_expand = source_expand(s3, centre; basis_order = 7)

    xs = [centre + 0.1 .* [cos(τ),sin(τ)] for τ = 0.0:0.3:1.5]
    @test norm([field(s3,x,ω) - s3_expand(x,ω) for x in xs]) < 1e-7*norm([field(s3,x,ω) for x in xs])

    source = plane_source(a2_host; position = SVector(-10.0,0.0), direction = SVector(1.0,0.0), amplitude = 1.0)
    s_expand = source_expand(source, centre; basis_order = 4)
    @test norm([field(source,x,ω) - s_expand(x,ω) for x in xs]) < 2e-9*norm([field(source,x,ω) for x in xs])

    # Test DimensionMismatch between 3D and 2D

    direction = SVector(-10.0,1.0)
    position = SVector(rand(2)...)
    amplitude = rand() + rand() * im

    a3_host = Acoustic(3, ρ = 1.3, c = 1.5 + 0.0im)
    @test_throws(DimensionMismatch,PlaneSource(a3_host; direction = direction))
    @test_throws(DimensionMismatch,PlaneSource(a3_host; amplitude = [1.0,2.0]))

    # Test equivalence between plane-sources
    psource = PlaneSource(a2_host; direction = direction, position =  position, amplitude = amplitude)

    source = plane_source(a2_host; direction = direction, position =  position, amplitude = amplitude)

    pf = field(psource)
    sf = field(source)

    xs = [rand(2) for i = 1:10]
    ωs = [rand() for i = 1:10]

    errs = map(xs,ωs) do x,ω
        abs(pf(x,ω) - sf(x,ω))
    end

    @test maximum(errs) < 2.0 * eps(Float64)

    # test constructors for 3D acoustics checks
    a3_host = Acoustic(3, ρ = 1.3, c = 1.5 + 0.0im)

    direction = rand(3)
    position = rand(3)
    amplitude = rand() + 0.1

    # check expansion in regular spherical waves
    psource = plane_source(a3_host; direction = direction, position =  position, amplitude = amplitude)
    pointsource = point_source(a3_host, position, amplitude)

    # Check that the field converges to its regular basis expansion around centre
    k = real(ω / psource.medium.c)
    centre = position + (4.0pi / k) .* rand(3)
    x = centre + (0.1 / k) .* rand(3)

    s_expand = source_expand(psource, x; basis_order = 4)

    # test if the source is equal to its series expansion in a regular basis.
    @test norm(field(psource,x,ω) - s_expand(x,ω)) < 1e-10 * abs(field(psource,x,ω))

    # create source through a regular expansion of spherical waves

    order = 3
    len = basisorder_to_basislength(typeof(a3_host),order)

    position = rand(3)

    source = regular_spherical_source(a3_host, rand(len) ./ (1:len);
        amplitude = rand(Complex{Float64}),
        position = position
    );

    k = real(ω / source.medium.c)
    centre = position + (4.0pi / k) .* rand(3)
    x = centre + (0.1 / k) .* rand(3)

    s_expand = source_expand(source, x; basis_order = 4)

    @test norm(field(source,x,ω) - s_expand(x,ω)) < 1e-10 * abs(field(source,x,ω))

end


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct HydroDispatchCascade <: HydroCascade
        name::String
        available::Bool
        bus::PSY.Bus
        active_power::Float64
        reactive_power::Float64
        rating::Float64
        prime_mover::PSY.PrimeMovers
        active_power_limits::NamedTuple{(:min, :max), Tuple{Float64, Float64}}
        reactive_power_limits::Union{Nothing, PSY.Min_Max}
        ramp_limits::Union{Nothing, NamedTuple{(:up, :down), Tuple{Float64, Float64}}}
        time_limits::Union{Nothing, NamedTuple{(:up, :down), Tuple{Float64, Float64}}}
        base_power::Float64
        operation_cost::PSY.OperationalCost
        storage_target::Float64
        conversion_factor::Float64
        upstream::Vector{NamedTuple{(:unit, :lag, :multiplier), Tuple{PSY.HydroGen, Int64, Float64}}}
        services::Vector{PSY.Service}
        dynamic_injector::Union{Nothing, PSY.DynamicInjection}
        ext::Dict{String, Any}
        time_series_container::InfrastructureSystems.TimeSeriesContainer
        internal::IS.InfrastructureSystemsInternal
    end

# Arguments
- `name::String`
- `available::Bool`
- `bus::PSY.Bus`
- `active_power::Float64`
- `reactive_power::Float64`, validation range: `reactive_power_limits`, action if invalid: `warn`
- `rating::Float64`: Thermal limited MVA Power Output of the unit. <= Capacity, validation range: `(0, nothing)`, action if invalid: `error`
- `prime_mover::PSY.PrimeMovers`: Prime mover technology according to EIA 923
- `active_power_limits::NamedTuple{(:min, :max), Tuple{Float64, Float64}}`
- `reactive_power_limits::Union{Nothing, PSY.Min_Max}`, action if invalid: `warn`
- `ramp_limits::Union{Nothing, NamedTuple{(:up, :down), Tuple{Float64, Float64}}}`: ramp up and ramp down limits in MW (in component base per unit) per minute, validation range: `(0, nothing)`, action if invalid: `error`
- `time_limits::Union{Nothing, NamedTuple{(:up, :down), Tuple{Float64, Float64}}}`: Minimum up and Minimum down time limits in hours, validation range: `(0, nothing)`, action if invalid: `error`
- `base_power::Float64`: Base power of the unit in MVA, validation range: `(0, nothing)`, action if invalid: `warn`
- `operation_cost::PSY.OperationalCost`: Operation Cost of Generation [`OperationalCost`](@ref)
- `storage_target::Float64`: Storage target at the end of simulation as ratio of storage capacity.
- `conversion_factor::Float64`: Conversion factor from flow/volume to energy: m^3 -> p.u-hr.
- `upstream::Vector{NamedTuple{(:unit, :lag, :multiplier), Tuple{PSY.HydroGen, Int64, Float64}}}`: unit: upstream units; lag: duration in number of periods between upstream release and downstream availability; multiplier: relationship between upstream energy release and downstream energy availability
- `services::Vector{PSY.Service}`: Services that this device contributes to
- `dynamic_injector::Union{Nothing, PSY.DynamicInjection}`: corresponding dynamic injection device
- `ext::Dict{String, Any}`
- `time_series_container::InfrastructureSystems.TimeSeriesContainer`: internal time_series storage
- `internal::IS.InfrastructureSystemsInternal`: power system internal reference, do not modify
"""
mutable struct HydroDispatchCascade <: HydroCascade
    name::String
    available::Bool
    bus::PSY.Bus
    active_power::Float64
    reactive_power::Float64
    "Thermal limited MVA Power Output of the unit. <= Capacity"
    rating::Float64
    "Prime mover technology according to EIA 923"
    prime_mover::PSY.PrimeMovers
    active_power_limits::NamedTuple{(:min, :max), Tuple{Float64, Float64}}
    reactive_power_limits::Union{Nothing, PSY.Min_Max}
    "ramp up and ramp down limits in MW (in component base per unit) per minute"
    ramp_limits::Union{Nothing, NamedTuple{(:up, :down), Tuple{Float64, Float64}}}
    "Minimum up and Minimum down time limits in hours"
    time_limits::Union{Nothing, NamedTuple{(:up, :down), Tuple{Float64, Float64}}}
    "Base power of the unit in MVA"
    base_power::Float64
    "Operation Cost of Generation [`OperationalCost`](@ref)"
    operation_cost::PSY.OperationalCost
    "Storage target at the end of simulation as ratio of storage capacity."
    storage_target::Float64
    "Conversion factor from flow/volume to energy: m^3 -> p.u-hr."
    conversion_factor::Float64
    "unit: upstream units; lag: duration in number of periods between upstream release and downstream availability; multiplier: relationship between upstream energy release and downstream energy availability"
    upstream::Vector{
        NamedTuple{(:unit, :lag, :multiplier), Tuple{PSY.HydroGen, Int64, Float64}},
    }
    "Services that this device contributes to"
    services::Vector{PSY.Service}
    "corresponding dynamic injection device"
    dynamic_injector::Union{Nothing, PSY.DynamicInjection}
    ext::Dict{String, Any}
    "internal time_series storage"
    time_series_container::InfrastructureSystems.TimeSeriesContainer
    "power system internal reference, do not modify"
    internal::IS.InfrastructureSystemsInternal
end

function HydroDispatchCascade(
    name,
    available,
    bus,
    active_power,
    reactive_power,
    rating,
    prime_mover,
    active_power_limits,
    reactive_power_limits,
    ramp_limits,
    time_limits,
    base_power,
    operation_cost = PSY.TwoPartCost(0.0, 0.0),
    storage_target = 1.0,
    conversion_factor = 1.0,
    upstream = NamedTuple{(:unit, :lag, :multiplier), Tuple{PSY.HydroGen, Int64, Float64}}[],
    services = PSY.Device[],
    dynamic_injector = nothing,
    ext = Dict{String, Any}(),
    time_series_container = InfrastructureSystems.TimeSeriesContainer(),
)
    HydroDispatchCascade(
        name,
        available,
        bus,
        active_power,
        reactive_power,
        rating,
        prime_mover,
        active_power_limits,
        reactive_power_limits,
        ramp_limits,
        time_limits,
        base_power,
        operation_cost,
        storage_target,
        conversion_factor,
        upstream,
        services,
        dynamic_injector,
        ext,
        time_series_container,
        IS.InfrastructureSystemsInternal(),
    )
end

function HydroDispatchCascade(;
    name,
    available,
    bus,
    active_power,
    reactive_power,
    rating,
    prime_mover,
    active_power_limits,
    reactive_power_limits,
    ramp_limits,
    time_limits,
    base_power,
    operation_cost = PSY.TwoPartCost(0.0, 0.0),
    storage_target = 1.0,
    conversion_factor = 1.0,
    upstream = NamedTuple{(:unit, :lag, :multiplier), Tuple{PSY.HydroGen, Int64, Float64}}[],
    services = PSY.Device[],
    dynamic_injector = nothing,
    ext = Dict{String, Any}(),
    time_series_container = InfrastructureSystems.TimeSeriesContainer(),
    internal = IS.InfrastructureSystemsInternal(),
)
    HydroDispatchCascade(
        name,
        available,
        bus,
        active_power,
        reactive_power,
        rating,
        prime_mover,
        active_power_limits,
        reactive_power_limits,
        ramp_limits,
        time_limits,
        base_power,
        operation_cost,
        storage_target,
        conversion_factor,
        upstream,
        services,
        dynamic_injector,
        ext,
        time_series_container,
        internal,
    )
end

# Constructor for demo purposes; non-functional.
function HydroDispatchCascade(::Nothing)
    HydroDispatchCascade(;
        name = "init",
        available = false,
        bus = PSY.Bus(nothing),
        active_power = 0.0,
        reactive_power = 0.0,
        rating = 0.0,
        prime_mover = PSY.PrimeMovers.HY,
        active_power_limits = (min = 0.0, max = 0.0),
        reactive_power_limits = nothing,
        ramp_limits = nothing,
        time_limits = nothing,
        base_power = 0.0,
        operation_cost = PSY.TwoPartCost(nothing),
        storage_target = 0.0,
        conversion_factor = 0.0,
        upstream = NamedTuple{
            (:unit, :lag, :multiplier),
            Tuple{PSY.HydroGen, Int64, Float64},
        }[],
        services = PSY.Device[],
        dynamic_injector = nothing,
        ext = Dict{String, Any}(),
        time_series_container = InfrastructureSystems.TimeSeriesContainer(),
    )
end

InfrastructureSystems.get_name(value::HydroDispatchCascade) = value.name

PowerSystems.get_available(value::HydroDispatchCascade) = value.available

PowerSystems.get_bus(value::HydroDispatchCascade) = value.bus

PowerSystems.get_active_power(value::HydroDispatchCascade) =
    get_value(value, value.active_power)

PowerSystems.get_reactive_power(value::HydroDispatchCascade) =
    get_value(value, value.reactive_power)

PowerSystems.get_rating(value::HydroDispatchCascade) = get_value(value, value.rating)

PowerSystems.get_prime_mover(value::HydroDispatchCascade) = value.prime_mover

PowerSystems.get_active_power_limits(value::HydroDispatchCascade) =
    get_value(value, value.active_power_limits)

PowerSystems.get_reactive_power_limits(value::HydroDispatchCascade) =
    get_value(value, value.reactive_power_limits)

PowerSystems.get_ramp_limits(value::HydroDispatchCascade) =
    get_value(value, value.ramp_limits)

PowerSystems.get_time_limits(value::HydroDispatchCascade) = value.time_limits

PowerSystems.get_base_power(value::HydroDispatchCascade) = value.base_power

PowerSystems.get_operation_cost(value::HydroDispatchCascade) = value.operation_cost

PowerSystems.get_storage_target(value::HydroDispatchCascade) = value.storage_target

PowerSystems.get_conversion_factor(value::HydroDispatchCascade) = value.conversion_factor
"""Get [`HydroDispatchCascade`](@ref) `upstream`."""
get_upstream(value::HydroDispatchCascade) = value.upstream

PowerSystems.get_services(value::HydroDispatchCascade) = value.services

PowerSystems.get_dynamic_injector(value::HydroDispatchCascade) = value.dynamic_injector

PowerSystems.get_ext(value::HydroDispatchCascade) = value.ext

InfrastructureSystems.get_time_series_container(value::HydroDispatchCascade) =
    value.time_series_container

PowerSystems.get_internal(value::HydroDispatchCascade) = value.internal

InfrastructureSystems.set_name!(value::HydroDispatchCascade, val) = value.name = val

PowerSystems.set_available!(value::HydroDispatchCascade, val) = value.available = val

PowerSystems.set_bus!(value::HydroDispatchCascade, val) = value.bus = val

PowerSystems.set_active_power!(value::HydroDispatchCascade, val) = value.active_power = val

PowerSystems.set_reactive_power!(value::HydroDispatchCascade, val) =
    value.reactive_power = val

PowerSystems.set_rating!(value::HydroDispatchCascade, val) = value.rating = val

PowerSystems.set_prime_mover!(value::HydroDispatchCascade, val) = value.prime_mover = val

PowerSystems.set_active_power_limits!(value::HydroDispatchCascade, val) =
    value.active_power_limits = val

PowerSystems.set_reactive_power_limits!(value::HydroDispatchCascade, val) =
    value.reactive_power_limits = val

PowerSystems.set_ramp_limits!(value::HydroDispatchCascade, val) = value.ramp_limits = val

PowerSystems.set_time_limits!(value::HydroDispatchCascade, val) = value.time_limits = val

PowerSystems.set_base_power!(value::HydroDispatchCascade, val) = value.base_power = val

PowerSystems.set_operation_cost!(value::HydroDispatchCascade, val) =
    value.operation_cost = val

PowerSystems.set_storage_target!(value::HydroDispatchCascade, val) =
    value.storage_target = val

PowerSystems.set_conversion_factor!(value::HydroDispatchCascade, val) =
    value.conversion_factor = val
"""Set [`HydroDispatchCascade`](@ref) `upstream`."""
set_upstream!(value::HydroDispatchCascade, val) = value.upstream = val

PowerSystems.set_services!(value::HydroDispatchCascade, val) = value.services = val

PowerSystems.set_ext!(value::HydroDispatchCascade, val) = value.ext = val

InfrastructureSystems.set_time_series_container!(value::HydroDispatchCascade, val) =
    value.time_series_container = val


text:
push!(LOAD_PATH, "@stdlib")
import Pkg; Pkg.add("Conda");
using Conda
try
    run(`which nvcc`)
    ENV["GPU"] = 1
    Pkg.build("ADCME")
catch
end
using ADCME

@info "Install Boost"
CONDA = get_conda()
run(`$CONDA install libboost==1.73.0=h3ff78a5_11`)
# run(`$CONDA install boost==1.73.0`)

@info "Install AMGCL"
UNZIP = joinpath(ADCME.BINDIR, "unzip")
if !isdir("$(@__DIR__)/amgcl")
    download("https://github.com/ddemidov/amgcl/archive/master.zip", "$(@__DIR__)/amgcl.zip")
    run(`$UNZIP -o $(@__DIR__)/amgcl.zip -d $(@__DIR__)`)
    mv("$(@__DIR__)/amgcl-master","$(@__DIR__)/amgcl", force=true)
    rm("$(@__DIR__)/amgcl.zip")
end

@info "Build Custom Operators"
change_directory("CustomOps/build")
require_file("build.ninja") do
    ADCME.cmake()
end
ADCME.make()