text:
# Descriptive Statistics

#############################
#
#   Location
#
#############################

# Geometric mean
"""
    geomean(a)

Return the geometric mean of a collection.
"""
geomean(a) = exp(mean(log, a))

# Harmonic mean
"""
    harmmean(a)

Return the harmonic mean of a collection.
"""
harmmean(a) = inv(mean(inv, a))

# Generalized mean
"""
    genmean(a, p)

Return the generalized/power mean with exponent `p` of a real-valued array,
i.e. ``\\left( \\frac{1}{n} \\sum_{i=1}^n a_i^p \\right)^{\\frac{1}{p}}``, where `n = length(a)`.
It is taken to be the geometric mean when `p == 0`.
"""
function genmean(a, p::Real)
    if p == 0
        return geomean(a)
    end

    # At least one of `x` or `p` must not be an int to avoid domain errors when `p` is a negative int.
    # We choose `x` in order to exploit exponentiation by squaring when `p` is an int.
    r = mean(a) do x
        float(x)^p
    end
    return r^inv(p)
end

# compute mode, given the range of integer values
"""
    mode(a, [r])
    mode(a::AbstractArray, wv::AbstractWeights)

Return the mode (most common number) of an array, optionally
over a specified range `r` or weighted via a vector `wv`.
If several modes exist, the first one (in order of appearance) is returned.
"""
function mode(a::AbstractArray{T}, r::UnitRange{T}) where T<:Integer
    isempty(a) && throw(ArgumentError("mode is not defined for empty collections"))
    len = length(a)
    r0 = r[1]
    r1 = r[end]
    cnts = zeros(Int, length(r))
    mc = 0    # maximum count
    mv = r0   # a value corresponding to maximum count
    for i = 1:len
        @inbounds x = a[i]
        if r0 <= x <= r1
            @inbounds c = (cnts[x - r0 + 1] += 1)
            if c > mc
                mc = c
                mv = x
            end
        end
    end
    return mv
end

"""
    modes(a, [r])::Vector
    mode(a::AbstractArray, wv::AbstractWeights)::Vector

Return all modes (most common numbers) of an array, optionally over a
specified range `r` or weighted via vector `wv`.
"""
function modes(a::AbstractArray{T}, r::UnitRange{T}) where T<:Integer
    r0 = r[1]
    r1 = r[end]
    n = length(r)
    cnts = zeros(Int, n)
    # find the maximum count
    mc = 0
    for i = 1:length(a)
        @inbounds x = a[i]
        if r0 <= x <= r1
            @inbounds c = (cnts[x - r0 + 1] += 1)
            if c > mc
                mc = c
            end
        end
    end
    # find all values corresponding to maximum count
    ms = T[]
    for i = 1:n
        @inbounds if cnts[i] == mc
            push!(ms, r[i])
        end
    end
    return ms
end

# compute mode over arbitrary iterable
function mode(a)
    isempty(a) && throw(ArgumentError("mode is not defined for empty collections"))
    cnts = Dict{eltype(a),Int}()
    # first element
    mc = 1
    mv, st = iterate(a)
    cnts[mv] = 1
    # find the mode along with table construction
    y = iterate(a, st)
    while y !== nothing
        x, st = y
        if haskey(cnts, x)
            c = (cnts[x] += 1)
            if c > mc
                mc = c
                mv = x
            end
        else
            cnts[x] = 1
            # in this case: c = 1, and thus c > mc won't happen
        end
        y = iterate(a, st)
    end
    return mv
end

function modes(a)
    isempty(a) && throw(ArgumentError("mode is not defined for empty collections"))
    cnts = Dict{eltype(a),Int}()
    # first element
    mc = 1
    x, st = iterate(a)
    cnts[x] = 1
    # find the mode along with table construction
    y = iterate(a, st)
    while y !== nothing
        x, st = y
        if haskey(cnts, x)
            c = (cnts[x] += 1)
            if c > mc
                mc = c
            end
        else
            cnts[x] = 1
            # in this case: c = 1, and thus c > mc won't happen
        end
        y = iterate(a, st)
    end
    # find values corresponding to maximum counts
    return [x for (x, c) in cnts if c == mc]
end

# Weighted mode of arbitrary vectors of values
function mode(a::AbstractVector, wv::AbstractWeights{T}) where T <: Real
    isempty(a) && throw(ArgumentError("mode is not defined for empty collections"))
    length(a) == length(wv) ||
        throw(ArgumentError("data and weight vectors must be the same size, got $(length(a)) and $(length(wv))"))

    # Iterate through the data
    mv = first(a)
    mw = first(wv)
    weights = Dict{eltype(a), T}()
    for (x, w) in zip(a, wv)
        _w = get!(weights, x, zero(T)) + w
        if _w > mw
            mv = x
            mw = _w
        end
        weights[x] = _w
    end

    return mv
end

function modes(a::AbstractVector, wv::AbstractWeights{T}) where T <: Real
    isempty(a) && throw(ArgumentError("mode is not defined for empty collections"))
    length(a) == length(wv) ||
        throw(ArgumentError("data and weight vectors must be the same size, got $(length(a)) and $(length(wv))"))

    # Iterate through the data
    mw = first(wv)
    weights = Dict{eltype(a), T}()
    for (x, w) in zip(a, wv)
        _w = get!(weights, x, zero(T)) + w
        if _w > mw
            mw = _w
        end
        weights[x] = _w
    end

    # find values corresponding to maximum counts
    return [x for (x, w) in weights if w == mw]
end

#############################
#
#   quantile and friends
#
#############################

"""
    percentile(x, p)

Return the `p`th percentile of a collection `x`, i.e. `quantile(x, p / 100)`.
"""
percentile(x, p) = quantile(x, p * 0.01)

"""
    nquantile(x, n::Integer)

Return the n-quantiles of collection `x`, i.e. the values which
partition `v` into `n` subsets of nearly equal size.

Equivalent to `quantile(x, [0:n]/n)`. For example, `nquantiles(x, 5)`
returns a vector of quantiles, respectively at `[0.0, 0.2, 0.4, 0.6, 0.8, 1.0]`.
"""
nquantile(x, n::Integer) = quantile(x, (0:n)/n)

"""
    quantilerank(itr, value; method=:inc)

Compute the quantile position in the [0, 1] interval of `value` relative to collection `itr`.

Different definitions can be chosen via the `method` keyword argument.
Let `count_less` be the number of elements of `itr` that are less than `value`,
`count_equal` the number of elements of `itr` that are equal to `value`, `n` the length of `itr`,
`greatest_smaller` the highest value below `value` and `smallest_greater` the lowest value above `value`.
Then `method` supports the following definitions:

- `:inc` (default): Return a value in the range 0 to 1 inclusive.
Return `count_less / (n - 1)` if `value ∈ itr`, otherwise apply interpolation based on
definition 7 of quantile in Hyndman and Fan (1996)
(equivalent to Excel `PERCENTRANK` and `PERCENTRANK.INC`).
This definition corresponds to the lower semi-continuous inverse of
[`quantile`](@ref) with its default parameters.

- `:exc`: Return a value in the range 0 to 1 exclusive.
Return `(count_less + 1) / (n + 1)` if `value ∈ itr` otherwise apply interpolation
based on definition 6 of quantile in Hyndman and Fan (1996)
(equivalent to Excel `PERCENTRANK.EXC`).

- `:compete`: Return `count_less / (n - 1)` if `value ∈ itr`, otherwise
return `(count_less - 1) / (n - 1)`, without interpolation
(equivalent to MariaDB `PERCENT_RANK`, dplyr `percent_rank`).

- `:tied`: Return `(count_less + count_equal/2) / n`, without interpolation.
Based on the definition in Roscoe, J. T. (1975)
(equivalent to `"mean"` kind of SciPy `percentileofscore`).

- `:strict`: Return `count_less / n`, without interpolation
(equivalent to `"strict"` kind of SciPy `percentileofscore`).

- `:weak`: Return `(count_less + count_equal) / n`, without interpolation
(equivalent to `"weak"` kind of SciPy `percentileofscore`).

!!! note
    An `ArgumentError` is thrown if `itr` contains `NaN` or `missing` values
    or if `itr` contains fewer than two elements.

# References
Roscoe, J. T. (1975). [Fundamental Research Statistics for the Behavioral Sciences]
(http://www.bryanburnham.net/wp-content/uploads/2014/07/Fundamental-Statistics-for-the-Behavioral-Sciences-v2.0.pdf#page=57)",
2nd ed., New York : Holt, Rinehart and Winston.

Hyndman, R.J and Fan, Y. (1996) "[Sample Quantiles in Statistical Packages]
(https://www.amherst.edu/media/view/129116/original/Sample+Quantiles.pdf)",
*The American Statistician*, Vol. 50, No. 4, pp. 361-365.

# Examples
```julia
julia> using StatsBase

julia> v1 = [1, 1, 1, 2, 3, 4, 8, 11, 12, 13];

julia> v2 = [1, 2, 3, 5, 6, missing, 8];

julia> v3 = [1, 2, 3, 4, 4, 5, 6, 7, 8, 9];

julia> quantilerank(v1, 2)
0.3333333333333333

julia> quantilerank(v1, 2, method=:exc), quantilerank(v1, 2, method=:tied)
(0.36363636363636365, 0.35)

# use `skipmissing` for vectors with missing entries.
julia> quantilerank(skipmissing(v2), 4)
0.5

# use broadcasting with `Ref` to compute quantile rank for multiple values
julia> quantilerank.(Ref(v3), [4, 8])
2-element Vector{Float64}:
 0.3333333333333333
 0.8888888888888888
```
"""
function quantilerank(itr, value; method::Symbol=:inc)
    ((value isa Number && isnan(value)) || ismissing(value)) &&
        throw(ArgumentError("`value` cannot be NaN or missing"))
    any(x -> ismissing(x) || (x isa Number && isnan(x)), itr) &&
        throw(ArgumentError("`itr` cannot contain missing or NaN entries"))

    count_less = count_equal = n = 0
    greatest_smaller = smallest_greater = value
    for x in itr
        if x == value
            count_equal += 1
        elseif x < value
            count_less += 1
            if greatest_smaller == value || greatest_smaller < x
                greatest_smaller = x
            end
        else
            if smallest_greater == value || smallest_greater > x
                smallest_greater = x
            end
        end
        n += 1
    end

    n == 0 && throw(ArgumentError("`itr` is empty. Pass a collection with at least two elements"))
    n == 1 && throw(ArgumentError("`itr` has only 1 value. Pass a collection with at least two elements"))

    if method == :inc
        if greatest_smaller == value
            return 0.0
        elseif count_equal > 0
            return count_less / (n - 1)
        elseif smallest_greater == value
            return 1.0
        else
            lower = (count_less - 1) / (n - 1)
            upper = count_less / (n - 1)
            ratio = (value - greatest_smaller) / (smallest_greater - greatest_smaller)
            return lower + ratio * (upper - lower)
        end
    elseif method == :exc
        if count_less == 0 && count_equal == 0
            return 0.0
        elseif count_less == 0
            return 1.0 / (n + 1)
        elseif count_equal > 0
            return (count_less + 1) / (n + 1)
        elseif smallest_greater == value
            return 1.0
        else
            lower = count_less / (n + 1)
            upper = (count_less + 1) / (n + 1)
            ratio = (value - greatest_smaller) / (smallest_greater - greatest_smaller)
            return lower + ratio * (upper - lower)
        end
    elseif method == :compete
        if value > maximum(itr)
            return 1.0
        elseif value ≤ minimum(itr)
            return 0.0
        else
            value ∈ itr && (count_less += 1)
            return (count_less - 1) / (n - 1)
        end
    elseif method == :tied
        return (count_less + count_equal/2) / n
    elseif method == :strict
        return count_less / n
    elseif method == :weak
        return (count_less + count_equal) / n
    else
        throw(ArgumentError("method=:$method is not valid. Pass :inc, :exc, :compete, :tied, :strict or :weak."))
    end
end

"""
    percentilerank(itr, value; method=:inc)

Return the `q`th percentile of `value` in collection `itr`, i.e. [`quantilerank(itr, value)`](@ref) * 100.

See the [`quantilerank`](@ref) docstring for more details.
"""
percentilerank(itr, value; method::Symbol=:inc) = quantilerank(itr, value, method=method) * 100

#############################
#
#   Dispersion
#
#############################

# span, i.e. the range minimum(x):maximum(x)
"""
    span(x)

Return the span of a collection, i.e. the range `minimum(x):maximum(x)`.
The minimum and maximum of `x` are computed in one pass using `extrema`.
"""
span(x) = ((a, b) = extrema(x); a:b)

# Variation coefficient: std / mean
"""
    variation(x, m=mean(x))

Return the coefficient of variation of collection `x`, optionally specifying
a precomputed mean `m`. The coefficient of variation is the ratio of the
standard deviation to the mean.
"""
variation(x, m) = stdm(x, m) / m
variation(x) = ((m, s) = mean_and_std(x); s/m)

# Standard error of the mean: std / sqrt(len)
# Code taken from var in the Statistics stdlib module

# faster computation of real(conj(x)*y)
realXcY(x::Real, y::Real) = x*y
realXcY(x::Complex, y::Complex) = real(x)*real(y) + imag(x)*imag(y)

"""
    sem(x; mean=nothing)
    sem(x::AbstractArray[, weights::AbstractWeights]; mean=nothing)

Return the standard error of the mean for a collection `x`.
A pre-computed `mean` may be provided.

When not using weights, this is the (sample) standard deviation
divided by the sample size. If weights are used, the
variance of the sample mean is calculated as follows:

* `AnalyticWeights`: Not implemented.
* `FrequencyWeights`: ``\\frac{\\sum_{i=1}^n w_i (x_i - \\bar{x_i})^2}{(\\sum w_i) (\\sum w_i - 1)}``
* `ProbabilityWeights`: ``\\frac{n}{n-1} \\frac{\\sum_{i=1}^n w_i^2 (x_i - \\bar{x_i})^2}{\\left( \\sum w_i \\right)^2}``

The standard error is then the square root of the above quantities.

# References

Carl-Erik Särndal, Bengt Swensson, Jan Wretman (1992). Model Assisted Survey Sampling.
New York: Springer. pp. 51-53.
"""
function sem(x; mean=nothing)
    if isempty(x)
        # Return the NaN of the type that we would get for a nonempty x
        T = eltype(x)
        _mean = mean === nothing ? zero(T) / 1 : mean
        z = abs2(zero(T) - _mean)
        return oftype((z + z) / 2, NaN)
    elseif mean === nothing
        n = 0
        y = iterate(x)
        value, state = y
        # Use Welford algorithm as seen in (among other places)
        # Knuth's TAOCP, Vol 2, page 232, 3rd edition.
        _mean = value / 1
        sse = real(zero(_mean))
        while y !== nothing
            value, state = y
            y = iterate(x, state)
            n += 1
            new_mean = _mean + (value - _mean) / n
            sse += realXcY(value - _mean, value - new_mean)
            _mean = new_mean
        end
    else
        n = 1
        y = iterate(x)
        value, state = y
        sse = abs2(value - mean)
        while (y = iterate(x, state)) !== nothing
            value, state = y
            n += 1
            sse += abs2(value - mean)
        end
    end
    variance = sse / (n - 1)
    return sqrt(variance / n)
end

function sem(x::AbstractArray; mean=nothing)
    if isempty(x)
        # Return the NaN of the type that we would get for a nonempty x
        T = eltype(x)
        _mean = mean === nothing ? zero(T) / 1 : mean
        z = abs2(zero(T) - _mean)
        return oftype((z + z) / 2, NaN)
    end
    return sqrt(var(x; mean=mean, corrected=true) / length(x))
end

function sem(x::AbstractArray, weights::UnitWeights; mean=nothing)
    if length(x) ≠ length(weights)
        throw(DimensionMismatch("array and weights do not have the same length"))
    end
    return sem(x; mean=mean)
end

# Weighted methods for the above
sem(x::AbstractArray, weights::FrequencyWeights; mean=nothing) =
    sqrt(var(x, weights; mean=mean, corrected=true) / sum(weights))

function sem(x::AbstractArray, weights::ProbabilityWeights; mean=nothing)
    if isempty(x)
        # Return the NaN of the type that we would get for a nonempty x
        return var(x, weights; mean=mean, corrected=true) / 0
    else
        _mean = mean === nothing ? Statistics.mean(x, weights) : mean
        # sum of squared errors = sse
        sse = sum(Broadcast.instantiate(Broadcast.broadcasted(x, weights) do x_i, w
            return abs2(w * (x_i - _mean))
        end))
        n = count(!iszero, weights)
        return sqrt(sse * n / (n - 1)) / sum(weights)
    end
end

# Median absolute deviation
@irrational mad_constant 1.4826022185056018 BigFloat("1.482602218505601860547076529360423431326703202590312896536266275245674447622701")

"""
    mad(x; center=median(x), normalize=true)

Compute the median absolute deviation (MAD) of collection `x` around `center`
(by default, around the median).

If `normalize` is set to `true`, the MAD is multiplied by
`1 / quantile(Normal(), 3/4) ≈ 1.4826`, in order to obtain a consistent estimator
of the standard deviation under the assumption that the data is normally distributed.
"""
function mad(x; center=nothing, normalize::Union{Bool, Nothing}=nothing, constant=nothing)
    isempty(x) && throw(ArgumentError("mad is not defined for empty arrays"))
    T = eltype(x)
    # Knowing the eltype allows allocating a single array able to hold both original values
    # and differences from the center, instead of two arrays
    S = isconcretetype(T) ? promote_type(T, typeof(middle(zero(T)))) : T
    x2 = x isa AbstractArray ? copyto!(similar(x, S), x) : collect(S, x)
    c = center === nothing ? median!(x2) : center
    if isconcretetype(T)
        x2 .= abs.(x2 .- c)
    else
        x2 = abs.(x2 .- c)
    end
    m = median!(x2)
    if normalize isa Nothing
        Base.depwarn("the `normalize` keyword argument will be false by default in future releases: set it explicitly to silence this deprecation", :mad)
        normalize = true
    end
    if !isa(constant, Nothing)
        Base.depwarn("keyword argument `constant` is deprecated, use `normalize` instead or apply the multiplication directly", :mad)
        m * constant
    elseif normalize
        m * mad_constant
    else
        m
    end
end

"""
    StatsBase.mad!(x; center=median!(x), normalize=true)

Compute the median absolute deviation (MAD) of array `x` around `center`
(by default, around the median), overwriting `x` in the process.
`x` must be able to hold values of generated by calling `middle` on its elements
(for example an integer vector is not appropriate since `middle` can produce
non-integer values).

If `normalize` is set to `true`, the MAD is multiplied by
`1 / quantile(Normal(), 3/4) ≈ 1.4826`, in order to obtain a consistent estimator
of the standard deviation under the assumption that the data is normally distributed.
"""
function mad!(x::AbstractArray;
              center=median!(x),
              normalize::Union{Bool,Nothing}=true,
              constant=nothing)
    isempty(x) && throw(ArgumentError("mad is not defined for empty arrays"))
    x .= abs.(x .- center)
    m = median!(x)
    if normalize isa Nothing
        Base.depwarn("the `normalize` keyword argument will be false by default in future releases: set it explicitly to silence this deprecation", :mad)
        normalize = true
    end
    if !isa(constant, Nothing)
        Base.depwarn("keyword argument `constant` is deprecated, use `normalize` instead or apply the multiplication directly", :mad)
        m * constant
    elseif normalize
        m * mad_constant
    else
        m
    end
end

# Interquartile range
"""
    iqr(x)

Compute the interquartile range (IQR) of collection `x`, i.e. the 75th percentile
minus the 25th percentile.
"""
iqr(x) = (q = quantile(x, [.25, .75]); q[2] - q[1])

# Generalized variance
"""
    genvar(X)

Compute the generalized sample variance of `X`. If `X` is a vector, one-column matrix,
or other iterable, this is equivalent to the sample variance.
Otherwise if `X` is a matrix, this is equivalent to the determinant of the covariance
matrix of `X`.

!!! note
    The generalized sample variance will be 0 if the columns of the matrix of deviations
    are linearly dependent.
"""
genvar(X::AbstractMatrix) = size(X, 2) == 1 ? var(vec(X)) : det(cov(X))
genvar(itr) = var(itr)

# Total variation
"""
    totalvar(X)

Compute the total sample variance of `X`. If `X` is a vector, one-column matrix,
or other iterable, this is equivalent to the sample variance.
Otherwise if `X` is a matrix, this is equivalent to the sum of the diagonal elements
of the covariance matrix of `X`.
"""
totalvar(X::AbstractMatrix) = sum(var(X, dims=1))
totalvar(itr) = var(itr)

#############################
#
#   Z-scores
#
#############################

function _zscore!(Z::AbstractArray, X::AbstractArray, μ::Real, σ::Real)
    # Z and X are assumed to have the same size
    iσ = inv(σ)
    if μ == zero(μ)
        for i = 1 : length(X)
            @inbounds Z[i] = X[i] * iσ
        end
    else
        for i = 1 : length(X)
            @inbounds Z[i] = (X[i] - μ) * iσ
        end
    end
    return Z
end

@generated function _zscore!(Z::AbstractArray{S,N}, X::AbstractArray{T,N},
                             μ::AbstractArray, σ::AbstractArray) where {S,T,N}
    quote
        # Z and X are assumed to have the same size
        # μ and σ are assumed to have the same size, that is compatible with size(X)
        siz1 = size(X, 1)
        @nextract $N ud d->size(μ, d)
        if size(μ, 1) == 1 && siz1 > 1
            @nloops $N i d->(d>1 ? (1:size(X,d)) : (1:1)) d->(j_d = ud_d ==1 ? 1 : i_d) begin
                v = (@nref $N μ j)
                c = inv(@nref $N σ j)
                for i_1 = 1:siz1
                    (@nref $N Z i) = ((@nref $N X i) - v) * c
                end
            end
        else
            @nloops $N i X d->(j_d = ud_d ==1 ? 1 : i_d) begin
                (@nref $N Z i) = ((@nref $N X i) - (@nref $N μ j)) / (@nref $N σ j)
            end
        end
        return Z
    end
end

function _zscore_chksize(X::AbstractArray, μ::AbstractArray, σ::AbstractArray)
    size(μ) == size(σ) || throw(DimensionMismatch("μ and σ should have the same size."))
    for i=1:ndims(X)
        dμ_i = size(μ,i)
        (dμ_i == 1 || dμ_i == size(X,i)) || throw(DimensionMismatch("X and μ have incompatible sizes."))
    end
end

"""
    zscore!([Z], X, μ, σ)

Compute the z-scores of an array `X` with mean `μ` and standard deviation `σ`.
z-scores are the signed number of standard deviations above the mean that an
observation lies, i.e. ``(x - μ) / σ``.

If a destination array `Z` is provided, the scores are stored
in `Z` and it must have the same shape as `X`. Otherwise `X` is overwritten.
"""
function zscore!(Z::AbstractArray{ZT}, X::AbstractArray{T}, μ::Real, σ::Real) where {ZT<:AbstractFloat,T<:Real}
    size(Z) == size(X) || throw(DimensionMismatch("Z and X must have the same size."))
    _zscore!(Z, X, μ, σ)
end

function zscore!(Z::AbstractArray{<:AbstractFloat}, X::AbstractArray{<:Real},
                 μ::AbstractArray{<:Real}, σ::AbstractArray{<:Real})
    size(Z) == size(X) || throw(DimensionMismatch("Z and X must have the same size."))
    _zscore_chksize(X, μ, σ)
    _zscore!(Z, X, μ, σ)
end

zscore!(X::AbstractArray{<:AbstractFloat}, μ::Real, σ::Real) = _zscore!(X, X, μ, σ)

zscore!(X::AbstractArray{<:AbstractFloat}, μ::AbstractArray{<:Real}, σ::AbstractArray{<:Real}) =
    (_zscore_chksize(X, μ, σ); _zscore!(X, X, μ, σ))

"""
    zscore(X, [μ, σ])

Compute the z-scores of `X`, optionally specifying a precomputed mean `μ` and
standard deviation `σ`. z-scores are the signed number of standard deviations
above the mean that an observation lies, i.e. ``(x - μ) / σ``.

`μ` and `σ` should be both scalars or both arrays. The computation is broadcasting.
In particular, when `μ` and `σ` are arrays, they should have the same size, and
`size(μ, i) == 1  || size(μ, i) == size(X, i)` for each dimension.
"""
function zscore(X::AbstractArray{T}, μ::Real, σ::Real) where T<:Real
    ZT = typeof((zero(T) - zero(μ)) / one(σ))
    _zscore!(Array{ZT}(undef, size(X)), X, μ, σ)
end

function zscore(X::AbstractArray{T}, μ::AbstractArray{U}, σ::AbstractArray{S}) where {T<:Real,U<:Real,S<:Real}
    _zscore_chksize(X, μ, σ)
    ZT = typeof((zero(T) - zero(U)) / one(S))
    _zscore!(Array{ZT}(undef, size(X)), X, μ, σ)
end

zscore(X::AbstractArray{<:Real}) = ((μ, σ) = mean_and_std(X); zscore(X, μ, σ))
zscore(X::AbstractArray{<:Real}, dim::Int) = ((μ, σ) = mean_and_std(X, dim); zscore(X, μ, σ))

#############################
#
#   entropy and friends
#
#############################

"""
    entropy(p, [b])

Compute the entropy of a collection of probabilities `p`,
optionally specifying a real number `b` such that the entropy is scaled by `1/log(b)`.
Elements with probability 0 or 1 add 0 to the entropy.
"""
function entropy(p)
    if isempty(p)
        throw(ArgumentError("empty collections are not supported since they do not " *
                            "represent proper probability distributions"))
    end
    return -sum(xlogx, p)
end

entropy(p, b::Real) = entropy(p) / log(b)

"""
    renyientropy(p, α)

Compute the Rényi (generalized) entropy of order `α` of an array `p`.
"""
function renyientropy(p::AbstractArray{T}, α::Real) where T<:Real
    α < 0 && throw(ArgumentError("Order of Rényi entropy not legal, $(α) < 0."))

    s = zero(T)
    z = zero(T)
    scale = sum(p)

    if α ≈ 0
        for i = 1:length(p)
            @inbounds pi = p[i]
            if pi > z
                s += 1
            end
        end
        s = log(s / scale)
    elseif α ≈ 1
        for i = 1:length(p)
            @inbounds pi = p[i]
            if pi > z
                s -= pi * log(pi)
            end
        end
        s = s / scale
    elseif (isinf(α))
        s = -log(maximum(p))
    else # a normal Rényi entropy
        for i = 1:length(p)
            @inbounds pi = p[i]
            if pi > z
                s += pi ^ α
            end
        end
        s = log(s / scale) / (1 - α)
    end
    return s
end

"""
    crossentropy(p, q, [b])

Compute the cross entropy between `p` and `q`, optionally specifying a real
number `b` such that the result is scaled by `1/log(b)`.
"""
function crossentropy(p::AbstractArray{<:Real}, q::AbstractArray{<:Real})
    length(p) == length(q) || throw(DimensionMismatch("Inconsistent array length."))

    # handle empty collections
    if isempty(p)
        Base.depwarn(
            "support for empty collections will be removed since they do not " *
            "represent proper probability distributions",
            :crossentropy,
        )
        # return zero for empty arrays
        return xlogy(zero(eltype(p)), zero(eltype(q)))
    end

    # use pairwise summation (https://github.com/JuliaLang/julia/pull/31020)
    broadcasted = Broadcast.broadcasted(xlogy, vec(p), vec(q))
    return - sum(Broadcast.instantiate(broadcasted))
end

crossentropy(p::AbstractArray{<:Real}, q::AbstractArray{<:Real}, b::Real) =
    crossentropy(p,q) / log(b)

"""
    kldivergence(p, q, [b])

Compute the Kullback-Leibler divergence from `q` to `p`,
also called the relative entropy of `p` with respect to `q`,
that is the sum `pᵢ * log(pᵢ / qᵢ)`. Optionally a real number `b`
can be specified such that the divergence is scaled by `1/log(b)`.
"""
function kldivergence(p::AbstractArray{<:Real}, q::AbstractArray{<:Real})
    length(p) == length(q) || throw(DimensionMismatch("Inconsistent array length."))

    # handle empty collections
    if isempty(p)
        Base.depwarn(
            "support for empty collections will be removed since they do not "*
            "represent proper probability distributions",
            :kldivergence,
        )
        # return zero for empty arrays
        pzero = zero(eltype(p))
        qzero = zero(eltype(q))
        return xlogy(pzero, zero(pzero / qzero))
    end

    # use pairwise summation (https://github.com/JuliaLang/julia/pull/31020)
    broadcasted = Broadcast.broadcasted(vec(p), vec(q)) do pi, qi
        # handle pi = qi = 0, otherwise `NaN` is returned
        piqi = iszero(pi) && iszero(qi) ? zero(pi / qi) : pi / qi
        return xlogy(pi, piqi)
    end
    return sum(Broadcast.instantiate(broadcasted))
end

kldivergence(p::AbstractArray{<:Real}, q::AbstractArray{<:Real}, b::Real) =
    kldivergence(p,q) / log(b)

#############################
#
#   summary
#
#############################

struct SummaryStats{T<:Union{AbstractFloat,Missing}}
    mean::T
    min::T
    q25::T
    median::T
    q75::T
    max::T
    nobs::Int
    nmiss::Int
end

"""
    summarystats(a)

Compute summary statistics for a real-valued array `a`. Returns a
`SummaryStats` object containing the mean, minimum, 25th percentile,
median, 75th percentile, and maxmimum.
"""
function summarystats(a::AbstractArray{T}) where T<:Union{Real,Missing}
    # `mean` doesn't fail on empty input but rather returns `NaN`, so we can use the
    # return type to populate the `SummaryStats` structure.
    s = T >: Missing ? collect(skipmissing(a)) : a
    m = mean(s)
    R = typeof(m)
    n = length(a)
    ns = length(s)
    qs = if ns == 0
        R[NaN, NaN, NaN, NaN, NaN]
    elseif T >: Missing
        quantile!(s, [0.00, 0.25, 0.50, 0.75, 1.00])
    else
        quantile(s, [0.00, 0.25, 0.50, 0.75, 1.00])
    end
    SummaryStats{R}(m, qs..., n, n - ns)
end

function Base.show(io::IO, ss::SummaryStats)
    println(io, "Summary Stats:")
    @printf(io, "Length:         %i\n", ss.nobs)
    ss.nobs > 0 || return
    @printf(io, "Missing Count:  %i\n", ss.nmiss)
    @printf(io, "Mean:           %.6f\n", ss.mean)
    @printf(io, "Minimum:        %.6f\n", ss.min)
    @printf(io, "1st Quartile:   %.6f\n", ss.q25)
    @printf(io, "Median:         %.6f\n", ss.median)
    @printf(io, "3rd Quartile:   %.6f\n", ss.q75)
    @printf(io, "Maximum:        %.6f\n", ss.max)
end

"""
    describe(a)

Pretty-print the summary statistics provided by [`summarystats`](@ref):
the mean, minimum, 25th percentile, median, 75th percentile, and
maximum.
"""
DataAPI.describe(x) = describe(stdout, x)
function DataAPI.describe(io::IO, a::AbstractArray{T}) where T<:Union{Real,Missing}
    show(io, summarystats(a))
    println(io, "Type:           $(string(eltype(a)))")
end
function DataAPI.describe(io::IO, a::AbstractArray)
    println(io, "Summary Stats:")
    println(io, "Length:         $(length(a))")
    println(io, "Type:           $(string(eltype(a)))")
    println(io, "Number Unique:  $(length(unique(a)))")
    return
end


text:
using AWS.EC2
using AWS

include("config.jl")

instances = ec2_instances_by_tag("ClusterName", ec2_clustername; running_only=false)
println("Terminating EC2 nodes : ", instances)

ec2_terminate(instances)
ec2_show_status(instances)

ext: jl


text:
using Resumable
using Base.Test

# First, kick off the docker container that contains the reference Node
# Resumable.js implementation, so we can interact with it.
container_id = readchomp(`docker run --rm -p 3000:3000/tcp -d staticfloat/resumeablejs_test_server`)
info("Test server running in container $(container_id[1:10])")

# Generate a 5MB+ file that will cause us to upload 5 (and a little more) chunks
bigfile = open("bigfile", "w")
write(bigfile, rand(UInt8, 5*1024*1024 + 1337))

try
    # Upload the big file to the Node.JS server
    Resumable.Client.upload_file("bigfile", "http://localhost:3000/upload")

    # Do it again to ensure that we don't need to upload those pieces again
    Resumable.Client.upload_file("bigfile", "http://localhost:3000/upload")

    # TODO verify the big file uploaded properly
finally
    info("Showing logs for docker test container:")
    run(`docker logs $container_id`)

    info("Killing docker server container:")
    success(`docker kill $container_id`)
end


text:
@testset "SBML import and conversion" begin
    sbmlm = load_sbml_model(model_paths["ecoli_core_model.xml"])
    m = convert(CoreModel, sbmlm)

    @test size(stoichiometry(sbmlm)) == (92, 95)
    @test size(stoichiometry(m)) == (n_metabolites(sbmlm), n_reactions(sbmlm))
    @test length(m.S.nzval) == 380
    @test length.(bounds(sbmlm)) == (95, 95)
    @test length.(bounds(m)) == (95, 95)
    @test all([length(m.xl), length(m.xu), length(m.c)] .== 95)

    @test metabolites(m)[1:3] == ["M_succoa_c", "M_ac_c", "M_fru_b"]
    @test reactions(m)[1:3] == ["R_EX_fum_e", "R_ACONTb", "R_GLNS"]

    cm = convert(CoreModelCoupled, sbmlm)
    @test n_coupling_constraints(cm) == 0
end


text:
julia> """
       # Y combinator

       * `λf. (λx. f (x x)) (λx. f (x x))`
       """
       Y = f -> (x -> x(x))(y -> f((t...) -> y(y)(t...)))

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# ACS, 10/2020
# Average abs error function to be used in objective function
function cfun2(a, A)
    cf = mean(abs.(a-A))
    return cf
end
# Permutation function
function  pfun(R)
    Rs = copy(R)
    ts = sample(1:length(R), 2)
    Rs[reverse(ts)] = R[ts]
    return Rs
end
# objective function to be minimized
function efun6(a,b,c,b2,R,nlag=50,n2lag=25)
    A = autocor(R, 1:nlag)
    B = autocor(abs.(R), 1:n2lag)
    C = crosscor(R, abs.(R), -nlag:nlag)
    B2 = autocor(R.*R, 1:n2lag)
    Z = cfun2(a,A)+cfun2(b,B)+cfun2(c,C)+cfun2(b2,B2)
    return Z
end
# main opt loop: adapted from TISEN
function annelOpt(aR,bR,cR,bR2,cmin,rs0,nlag,n2lag,cgoal,maxA)
  temp = 1e-4
  maxT = 20000
  maxG = 2000
  maxS = 200
  wc   = 0.9
  afac = 0.5
  #maxA = 120e6
  rs   = copy(rs0)
  cost = copy(cmin)
  nTot  = 0
  nGood = 0
  nAll  = 0
  ra    = rs
  tini  = temp
  E     = zeros(convert(Int64,ceil(maxA/1e5)),1).+cmin
  tE    = copy(cmin)
  for tt=1:maxA
    nTot = nTot+1
    nAll = nAll+1
    if (nAll%100000) == 0
      println(nAll/1e6)
      E[convert(Int64,ceil(nAll/1e5))] = tE
    end
    cmax = cost-temp*log(rand())
    trs  = pfun(rs)
    tE   = efun6(aR,bR,cR,bR2,trs,nlag,n2lag)
    if tE<=cmax
      nGood = nGood+1
      rs    = trs
      cost  = tE
    end
    if cost<=cgoal
      ra = rs
      break
    end
  if ((nTot<maxT) && (nGood<maxG))
    if (cost<(cmin*wc))
      cmin = cost
      ra   = rs
    end
    continue
  end
  if ((temp==tini) && (nTot>1.5*nGood))
    tini = 10*temp
    temp = tini
  else
    if (nGood<=maxS)
      afac=sqrt(afac)
      maxT=maxT*sqrt(2)
      temp=tini
    else
      temp=temp*afac
    end
  end
  nTot  = 0
  nGood = 0
end
return ra,E
end

max_stars_repo_name: silvaac/SMCarticle

max_issues_repo_name: silvaac/SMCarticle

max_forks_repo_name: silvaac/SMCarticle


text:
# comparison
==(x::Posit8,y::Posit8) = ccall((:p8_eq, SoftPositPath), Bool, (Posit8,Posit8),x,y)
==(x::Posit16,y::Posit16) = ccall((:p16_eq, SoftPositPath), Bool, (Posit16,Posit16),x,y)
==(x::Posit32,y::Posit32) = ccall((:p32_eq, SoftPositPath), Bool, (Posit32,Posit32),x,y)

<(x::Posit8,y::Posit8) = ccall((:p8_lt, SoftPositPath), Bool, (Posit8,Posit8),x,y)
<(x::Posit16,y::Posit16) = ccall((:p16_lt, SoftPositPath), Bool, (Posit16,Posit16),x,y)
<(x::Posit32,y::Posit32) = ccall((:p32_lt, SoftPositPath), Bool, (Posit32,Posit32),x,y)

<=(x::Posit8,y::Posit8) = ccall((:p8_le, SoftPositPath), Bool, (Posit8,Posit8),x,y)
<=(x::Posit16,y::Posit16) = ccall((:p16_le, SoftPositPath), Bool, (Posit16,Posit16),x,y)
<=(x::Posit32,y::Posit32) = ccall((:p32_le, SoftPositPath), Bool, (Posit32,Posit32),x,y)

# for Posits with 1 exponent bit
for T in (:Posit8_1, :Posit16_1, :Posit24_1)
    @eval begin
        ==(x::$T,y::$T) = ccall((:pX1_eq, SoftPositPath), Bool, ($T,$T),x,y)
        <=(x::$T,y::$T) = ccall((:pX1_le, SoftPositPath), Bool, ($T,$T),x,y)
        <(x::$T,y::$T) = ccall((:pX1_lt, SoftPositPath), Bool, ($T,$T),x,y)
    end
end

# ==(x::T,y::T) where {T<:PositX1} = ccall((:pX1_eq, SoftPositPath), Bool, (T,T),x,y)
# <=(x::T,y::T) where {T<:PositX1} = ccall((:pX1_le, SoftPositPath), Bool, (T,T),x,y)
# <(x::T,y::T) where {T<:PositX1} = ccall((:pX1_lt, SoftPositPath), Bool, (T,T),x,y)

# for Posits with 2 exponent bits
for T in (:Posit8_2, :Posit16_2, :Posit24_2)
    @eval begin
        ==(x::$T,y::$T) = ccall((:pX2_eq, SoftPositPath), Bool, ($T,$T),x,y)
        <=(x::$T,y::$T) = ccall((:pX2_le, SoftPositPath), Bool, ($T,$T),x,y)
        <(x::$T,y::$T) = ccall((:pX2_lt, SoftPositPath), Bool, ($T,$T),x,y)
    end
end


text:
include("DislocationForces.jl")
include("DislocationMobility.jl")
include("DislocationTopology.jl")

max_stars_repo_name: dcelisgarza/DDD

max_forks_repo_name: dcelisgarza/DDD


text:
#Reference:
#https://syl1.gitbook.io/julia-language-a-concise-tutorial/language-core/custom-types

mutable struct Person
  myname::String
  age::Int64
end

struct Shoes
   shoesType::String
   colour::String
end

mutable struct Student
   s::Person
   school::String
   shoes::Shoes
   grades::Array{Float64,1}
end

function printMyActivity(self::Student)
   println("$(self.s.myname) studies at $(self.school) school")
end

struct Employee
   s::Person
   annualIncome::Float64
   company::String
   shoes::Shoes
   monthlyIncomes::Array{Float64,1}
end

function printMyActivity(self::Employee)
  println("I work at $(self.company) company")
end

gymShoes = Shoes("gym","white")
proShoes = Shoes("classical","brown")

Marc = Student(Person("Marc",15),"Nicholas School",gymShoes, [100, 90])
MrBrown = Employee(Person("Brown",45),1200.0,"ABC Corporation Inc.", proShoes, [250, 180])

John_person = Person("John",55)

John = Student(John_person,"Divinity School", proShoes, [80, 50, 95])

printMyActivity(Marc)
printMyActivity(MrBrown)
printMyActivity(John)

##
#Returning a Struct
function createStudent(name::String, age::Int64, school::String, shoes::Shoes, grades::Array{Float64,1})
   myStudent = Student(Person(name,age), school, gymShoes, grades)
   return myStudent
end

Mary = createStudent("Mary", 30, "Yale", gymShoes, [98.5, 85, 100])
printMyActivity(Mary)
print("Mary's grades are: $(Mary.grades)\n")

#Adding grades
push!(Mary.grades, 97.5)
print("Mary's grades are: $(Mary.grades)\n")

#Adding grades
push!(Mary.grades, 85.5, 88.8)
print("Mary's grades are: $(Mary.grades)\n")

Mary.school = "Nicholas School"

#Trying to modify a nonmutable struct
print("MrBrown's mobthly incomes are: $(MrBrown.monthlyIncomes)\n")
push!(MrBrown.monthlyIncomes, 850)
print("MrBrown's mobthly incomes are: $(MrBrown.monthlyIncomes)\n")

MrBrown.company = "GE"
FUCR.genoffon

C:\Users\rapiduser\Ali_Julia

max_stars_repo_name: power-systems-op/grace-arpae

max_issues_repo_name: power-systems-op/grace-arpae

max_forks_repo_name: power-systems-op/grace-arpae


text:
function link_f(net::FastNet,c,rp)
    @inbounds linkpos=net.kcstart[c]+rp
    @inbounds net.kid[linkpos]
end

function link_f(net::FastNet,rp)
    @inbounds net.kid[rp]
end

function linkstate_f(net::FastNet,kid)
    @inbounds net.kstate[kid]
end

function makelink_f!(net::FastNet,src,dst)
    @inbounds begin
        kid=link_f(net,net.L,1)
        net.ksrc[kid]=src
        net.kdst[kid]=dst
        cls=net.ttable[net.nstate[src],net.nstate[dst]]
        linkstate_f!(net,kid,cls)
        net.knexts[kid]=net.nout[src]
        net.knextd[kid]=net.nin[dst]
        net.nout[src]=kid
        net.nin[dst]=kid
    end
    kid
end

function destroylink_f!(net::FastNet,kid)
    lastsrc=0
    lastdst=0
    @inbounds begin
        srcn=net.ksrc[kid]
        dstn=net.kdst[kid]
        cur=net.nout[srcn]
        while cur!=kid
            lastsrc=cur
            cur=net.knexts[cur]
        end
        cur=net.nin[dstn]
        while cur!=kid
            lastdst=cur
            cur=net.knextd[cur]
        end
        if lastsrc!=0
            net.knexts[lastsrc]=net.knexts[kid]
        else
            net.nout[srcn]=net.knexts[kid]
        end
        if lastdst!=0
            net.knextd[lastdst]=net.knextd[kid]
        else
            net.nin[dstn]=net.knextd[kid]
        end
        linkstate_f!(net,kid,net.L)
    end
    nothing
end

function linkexists_f(net::FastNet,kid)
    @inbounds net.kstate[kid]!=net.L
end

function linksrc_f(net::FastNet,kid)
    @inbounds net.ksrc[kid]
end

function linkdst_f(net::FastNet,kid)
    @inbounds net.kdst[kid]
end

function nextlinkout_f(net::FastNet,kid)
    @inbounds net.knexts[kid]
end

function nextlinkin_f(net::FastNet,kid)
    @inbounds net.knextd[kid]
end

function countlinks_f(net::FastNet)
    @inbounds net.K-net.kclen[net.L]
end

function countlinks_f(net::FastNet,cls::Integer)
    @inbounds net.kclen[cls]
end

function linkcounts_f(net::FastNet)
    n=net.L-2
    @inbounds net.kclen[1:n]
end

function countlinks_f(net::FastNet,cls::Union{Vector,Tuple})
    ret=0
    @inbounds begin
        for c in cls
            ret+=net.kclen[c]
        end
    end
    ret
end

function randomlink_f(net::FastNet)
    link_f(net,rand(net.rng,1:countlinks_f(net)))
end

function randomlink_f(net::FastNet,cls::Integer)
    link_f(net,cls,rand(net.rng,1:countlinks_f(net,cls)))
end

function randomlink_f(net::FastNet,cls::Union{Array,Tuple})
    tot=countlinks_f(net,cls)
    r=rand(net.rng,1:tot)
    old=r
    i=1
    cr=old-countlinks_f(net,cls[i])
    while cr>0
        i+=1
        old=cr
        cr=old-countlinks_f(net,cls[i])
    end
    link_f(net,cls[i],old)
end


text:
# CSS style constructors

# ---------------- add <style></style>
"""
    tag_style( CssStylesStr::String ; type::String = "text/css" )

add a pair of <style> tags for your CSS style string.
"""
tag_style( CssStylesStr::String ; type::String = "text/css" ) = begin
    return PairedHtmlTag( "style", nothing, nothing, nothing,
    singleAttributes = String[],
    pairedAttributes = Dict{String,String}( "type" => add_doublequotations(type) ),
    content = [BlankHtmlTag(CssStylesStr),]  )
end # tag_style

# ---------------- link to external CSS file
"""
    link_css( href::String )

linking an external css file. returns a <link /> tag string.
"""
link_css( href::String ) = tag_link( "stylesheet", "text/css", href )

# ------------------------ JS script string
"""
    CSS_str(script)
"""
macro CSS_str(script)
    return quote
        tag_style( "\n" * string($script) * "\n" )
    end
end # CSS_str

#

max_stars_repo_name: Clpr/HtmlConstructor

max_issues_repo_name: Clpr/HtmlConstructor

max_forks_repo_name: Clpr/HtmlConstructor


text:
module DocsHelper
export @doc_signature, isnodocumentationfound, mygetdoc
using Markdown

"""
checks whether doc string is the default "No documentation found."
"""
function isnodocumentationfound(doc)
  doc isa Markdown.MD &&
  length(doc.content) >= 1 &&
  doc.content[1] isa Markdown.Paragraph &&
  length(doc.content[1].content) == 1 &&
  doc.content[1].content[1] == "No documentation found."
end

"""
    mygetdoc(binding, sig)

Like Docs.doc but returns nothing if no information is found
"""
function mygetdoc(binding::Base.Docs.Binding, sig::Type = Union{})
    if Base.Docs.defined(binding)
        result = Base.Docs.getdoc(Base.Docs.resolve(binding), sig)
        result === nothing || return result
    end
    results, groups = Base.Docs.DocStr[], Base.Docs.MultiDoc[]
    # Lookup `binding` and `sig` for matches in all modules of the docsystem.
    for mod in Base.Docs.modules
        dict = Base.Docs.meta(mod)
        if haskey(dict, binding)
            multidoc = dict[binding]
            push!(groups, multidoc)
            for msig in multidoc.order
                sig <: msig && push!(results, multidoc.docs[msig])
            end
        end
    end
    if isempty(groups)
        # When no `MultiDoc`s are found that match `binding` then we check whether `binding`
        # is an alias of some other `Binding`. When it is we then re-run `doc` with that
        # `Binding`, otherwise if it's not an alias then we generate a summary for the
        # `binding` and display that to the user instead.
        alias = Base.Docs.aliasof(binding)
        alias == binding ? nothing : mygetdoc(alias, sig)
    elseif isempty(results)
        # There was at least one match for `binding` while searching. If there weren't any
        # matches for `sig` then we don't return anything.
        nothing
    else
        # Get parsed docs and concatenate them.
        md = Base.Docs.catdoc(map(Base.Docs.parsedoc, results)...)
        # Save metadata in the generated markdown.
        if isa(md, Markdown.MD)
            md.meta[:results] = results
            md.meta[:binding] = binding
            md.meta[:typesig] = sig
        end
        return md
    end
end

# Some additional convenience `doc` methods that take objects rather than `Binding`s.
mygetdoc(obj::UnionAll) = mygetdoc(Base.unwrap_unionall(obj))
mygetdoc(object, sig::Type = Union{}) = mygetdoc(Base.Docs.aliasof(object, typeof(object)), sig)
mygetdoc(object, sig...)              = mygetdoc(object, Tuple{sig...})

"""
Base.Docs.@doc cannot parse more signatures directly, but only if a true function is defined

Hence we do our own version by using the code which is normally used to extract the signature out of
function definitions.
"""
macro doc_signature(doc, expr)
  Docs.objectdoc(__source__, __module__, doc, :(nothing), expr, signature(expr))
end

# TODO open a ticket for this change that tv are not added separately as signatures:
function signature!(tv, expr::Expr)
    is_macrocall = Docs.isexpr(expr, :macrocall)
    if is_macrocall || Docs.isexpr(expr, :call)
        sig = :(Union{Tuple{}})
        first_arg = is_macrocall ? 3 : 2 # skip function arguments
        for arg in expr.args[first_arg:end]
            Docs.isexpr(arg, :parameters) && continue
            if Docs.isexpr(arg, :kw) # optional arg
                push!(sig.args, :(Tuple{$(sig.args[end].args[2:end]...)}))
            end
            push!(sig.args[end].args, Docs.argtype(arg))
        end
        if Docs.isexpr(expr.args[1], :curly) && isempty(tv)
            append!(tv, tvar.(expr.args[1].args[2:end]))
        end
        # for i = length(tv):-1:1
        #     push!(sig.args, :(Tuple{$(tv[i].args[1])}))
        # end
        for i = length(tv):-1:1
            sig = Expr(:where, sig, tv[i])
        end
        return sig
    elseif Docs.isexpr(expr, :where)
        append!(tv, Docs.tvar.(expr.args[2:end]))
        return signature!(tv, expr.args[1])
    else
        return signature!(tv, expr.args[1])
    end
end
signature!(tv, @nospecialize(other)) = :(Union{})
signature(expr::Expr) = signature!([], expr)
signature(@nospecialize other) = signature!([], other)

end # module


text:
using Statistics

const MYEPS = 20*eps()

@testset "Bodies" begin

  bn = nothing
  @test numpts(bn) == 0

  p = Plate(1,101)
  dx, dy = diff(p)
  @test maximum(dx) ≈ minimum(dx) ≈ 0.01
  @test maximum(dy) ≈ minimum(dy) ≈ 0.0

  dxc, dyc = centraldiff(p)
  @test dxc[1] ≈ 0.005 && dxc[101] ≈ 0.005 && maximum(dxc[2:100]) ≈ 0.01
  @test maximum(dyc) ≈ minimum(dyc) ≈ 0.0
  @test sum(dlength(p)) ≈ sum(dlengthmid(p))

  c = Rectangle(1,2,101)
  dx, dy = diff(c)
  @test length(dx) == 600
  @test sum(dlength(c)) ≈ 12.0

  c = Square(1,0.01)
  @test isapprox(mean(dlength(c)),0.01,atol=1e-4)

  c = Ellipse(1,2,0.01)
  @test isapprox(mean(dlength(c)),0.01,atol=1e-4)

  c = Ellipse(1,2,100)
  nx, ny = normalmid(c)
  @test nx[1] == ny[26] == -nx[51] == -ny[76] == 1.0
  @test abs(sum(nx)) < 1000.0*eps(1.0)
  @test abs(sum(ny)) < 1000.0*eps(1.0)

end

@testset "Shifted rectangle" begin

  b = Rectangle(0.5,1.0,11,shifted=true)
  nx, ny = normalmid(b)
  nx, ny = normalmid(b)
  @test all(isapprox.(nx[1:10],0.0,atol=MYEPS)) &&
        all(isapprox.(nx[11:30],1.0,atol=MYEPS)) &&
        all(isapprox.(nx[31:40],0.0,atol=MYEPS)) &&
        all(isapprox.(nx[41:60],-1.0,atol=MYEPS))
  @test all(isapprox.(ny[1:10],-1.0,atol=MYEPS)) &&
        all(isapprox.(ny[11:30],0.0,atol=MYEPS)) &&
        all(isapprox.(ny[31:40],1.0,atol=MYEPS)) &&
        all(isapprox.(ny[41:60],0.0,atol=MYEPS))

  ds = dlengthmid(b)
  @test all(ds .≈ 0.1)

  T = RigidTransform((1.0,-1.0),π/2)
  T(b)

  nx, ny = normalmid(b)
  @test all(isapprox.(nx[1:10],1.0,atol=MYEPS)) &&
        all(isapprox.(nx[11:30],0.0,atol=MYEPS)) &&
        all(isapprox.(nx[31:40],-1.0,atol=MYEPS)) &&
        all(isapprox.(nx[41:60],0.0,atol=MYEPS))
  @test all(isapprox.(ny[1:10],0.0,atol=MYEPS)) &&
        all(isapprox.(ny[11:30],1.0,atol=MYEPS)) &&
        all(isapprox.(ny[31:40],0.0,atol=MYEPS)) &&
        all(isapprox.(ny[41:60],-1.0,atol=MYEPS))

  ds = dlengthmid(b)
  @test all(ds .≈ 0.1)

end

@testset "Lists" begin

  n1 = 101
  n2 = 201
  p = Plate(1,n1)
  c = Rectangle(1,2,n2)

  bl = BodyList([p,c])
  @test length(bl) == 2

  @test eltype(bl) == Body

  m1 = RigidBodyMotion(complex(0.0),0.0)
  m2 = RigidBodyMotion(RigidBodyTools.PitchHeave(1.0, 11.0, 0.2, 0.0, 0.0, 0.5, 1.0, 0.0))

  ml = RigidMotionList([m1,m2])
  @test length(ml) == 2
  @test eltype(ml) == RigidBodyMotion

  t1 = RigidTransform((0.0,0.0),0.0)
  t2 = RigidTransform((1.0,0.0),π/2)

  tl = RigidTransformList([t1,t2])
  push!(tl,t1)
  @test length(tl) == 3
  @test eltype(tl) == RigidTransform

  v = rand(numpts(bl))
  @test numpts(bl) == length(p)+length(c)
  @test v[1:length(p)] == view(v,bl,1)
  @test v[(length(p)+1):(length(p)+length(c))] == view(v,bl,2)

  tl = RigidTransformList([t1,t2])
  tl(bl)
  @test bl[1].cent == tl[1].trans
  @test bl[2].cent == tl[2].trans
  @test bl[1].α == tl[1].α
  @test bl[2].α == tl[2].α

  b1 = Rectangle(0.5,1.0,11,shifted=true)
  b2 = Rectangle(0.5,1.0,11,shifted=false)
  bl = BodyList()
  push!(bl,b1)
  push!(bl,b2)
  nx, ny = normalmid(bl)
  nxb1, nyb1 = normalmid(b1)
  nxb2, nyb2 = normalmid(b2)
  @test nx[1:60] == nxb1 && ny[1:60] == nyb1 && nx[61:120] == nxb2 && ny[61:120] == nyb2

end

@testset "Assign velocity" begin
    b = Circle(1.0,100)
    T = RigidTransform((rand(),rand()),0.0)
    T(b)

    ċ = rand(ComplexF64)
    m = RigidBodyMotion(ċ,1.0)
    u,v = assign_velocity(b,m,0.0)
    @test u[26] ≈ -1.0+real(ċ) atol = 1e-14
    @test v[51] ≈ -1.0+imag(ċ) atol = 1e-14

    u2, v2 = m(0.0,b)
    @test u == u2 && v == v2

    b2 = Circle(1.0,100)
    T2 = RigidTransform((rand(),rand()),0.0)
    T2(b2)
    ċ2 = rand(ComplexF64)
    m2 = RigidBodyMotion(ċ2,1.0)

    bl = BodyList([b,b2])
    ml = RigidMotionList([m,m2])
    u, v = assign_velocity(bl,ml,0.0)
    @test u[26] ≈ -1.0+real(ċ) atol = 1e-14
    @test u[126] ≈ -1.0+real(ċ2) atol = 1e-14
    @test v[151] ≈ -1.0+imag(ċ2) atol = 1e-14

    u2, v2 = ml(0.0,bl)
    @test u2 == u && v2 == v

    vel = rigidbodyvelocity(ml,0.0)
    @test vel[3] == 1.0
    @test vel[1]+im*vel[2] ≈ ċ atol = 1e-14
    @test vel[4]+im*vel[5] ≈ ċ2 atol = 1e-14
    @test vel[6] == 1.0

    x = rand(15)
    tl = RigidTransformList(x)
    @test tl[4].trans[1] == x[10]
    @test tl[4].trans[2] == x[11]
    @test tl[4].α == x[12]

end


text:
using BenchmarkTools, Plots, Plots.PlotMeasures, LaTeXStrings
using BenchmarkTools: minimum, median

SUITE = BenchmarkGroup()
model = "POWERTRAIN"
cases = ["DTN01", "DTN02", "DTN03", "DTN04", "DTN05", "DTN06"]

SUITE[model] = BenchmarkGroup()

include("Powertrain.jl")
validation = []

# ----------------------------------------
#  DTN01 (dense time)
# ----------------------------------------
prob_DTN01 = powertrain_homog(θ=1)
@time sol_DTN01 = _solve_powertrain(prob_DTN01) # warm-up run
property = req(sol_DTN01)
push!(validation, Int(property))
SUITE[model][cases[1]] = @benchmarkable _solve_powertrain($prob_DTN01)

# ----------------------------------------
#  DTN02 (dense time)
# ----------------------------------------
prob_DTN02 = powertrain_homog(θ=2)
@time sol_DTN02 = _solve_powertrain(prob_DTN02) # warm-up run
property = req(sol_DTN02)
push!(validation, Int(property))
SUITE[model][cases[2]] = @benchmarkable _solve_powertrain($prob_DTN02)

# ----------------------------------------
#  DTN03 (dense time)
# ----------------------------------------
prob_DTN03 = powertrain_homog(θ=3)
@time sol_DTN03 = _solve_powertrain(prob_DTN03) # warm-up run
property = req(sol_DTN03)
push!(validation, Int(property))
SUITE[model][cases[3]] = @benchmarkable _solve_powertrain($prob_DTN03)

# ----------------------------------------
#  DTN04 (dense time)
# ----------------------------------------
prob_DTN04 = powertrain_homog(θ=4)
@time sol_DTN04 = _solve_powertrain(prob_DTN04) # warm-up run
property = req(sol_DTN04)
push!(validation, Int(property))
SUITE[model][cases[4]] = @benchmarkable _solve_powertrain($prob_DTN04)

# ----------------------------------------
#  DTN05 (dense time)
# ----------------------------------------
prob_DTN05 = powertrain_homog(θ=5)
@time sol_DTN05 = _solve_powertrain(prob_DTN05) # warm-up run
property = req(sol_DTN05)
push!(validation, Int(property))
SUITE[model][cases[5]] = @benchmarkable _solve_powertrain($prob_DTN05)

# ----------------------------------------
#  DTN06 (dense time)
# ----------------------------------------
prob_DTN06 = powertrain_homog(θ=6)
@time sol_DTN06 = _solve_powertrain(prob_DTN06) # warm-up run
property = req(sol_DTN06)
push!(validation, Int(property))
SUITE[model][cases[6]] = @benchmarkable _solve_powertrain($prob_DTN06)

# ==============================================================================
# Execute benchmarks and save benchmark results
# ==============================================================================

# tune parameters
tune!(SUITE)

# run the benchmarks
results = run(SUITE, verbose=true)

# return the sample with the smallest time value in each test
println("minimum time for each benchmark:\n", minimum(results))

# return the median for each test
println("median time for each benchmark:\n", median(results))

# export runtimes
runtimes = Dict()
for (i, c) in enumerate(cases)
   t = median(results[model][c]).time * 1e-9
   runtimes[c] = round(t, digits=4)
end

for (i, c) in enumerate(cases)
   print(io, "JuliaReach, $model, $c, $(validation[i]), $(runtimes[c])\n")
end

# ==============================================================================
# Plot
# ==============================================================================

fig = Plots.plot()
Plots.plot!(fig, sol_DTN03, vars=(1, 3),
           color=:blue, alpha=0.5, lw=1.0, linecolor=:black,
           tickfont=font(30, "Times"), guidefontsize=45,
           xlab=L"x_1",
           ylab=L"x_3",
           xtick=([-0.1, 0.0, 0.1, 0.2],
                  [L"-0.1", L"0.0", L"-0.1", L"0.2"]),
           ytick=([20, 40, 60, 80],
                  [L"20", L"40", L"60", L"80"]),
           xlims=(-0.1, 0.2), ylims=(-5.0, 85),
           bottom_margin=6mm, left_margin=5mm, right_margin=12mm, top_margin=3mm,
           size=(1000, 1000))
savefig(fig, joinpath(@__DIR__, "ARCH-COMP21-JuliaReach-Powertrain_DTN03"))

fig = Plots.plot()
Plots.plot!(fig, sol_DTN05, vars=(1, 3),
           color=:blue, alpha=0.5, lw=1.0, linecolor=:black,
           tickfont=font(30, "Times"), guidefontsize=45,
           xlab=L"x_1",
           ylab=L"x_3",
           xtick=([-0.1, 0.0, 0.1, 0.2],
                  [L"-0.1", L"0.0", L"-0.1", L"0.2"]),
           ytick=([20, 40, 60, 80],
                  [L"20", L"40", L"60", L"80"]),
           xlims=(-0.1, 0.2), ylims=(-5.0, 85),
           bottom_margin=6mm, left_margin=5mm, right_margin=12mm, top_margin=3mm,
           size=(1000, 1000))
savefig(fig, joinpath(@__DIR__, "ARCH-COMP21-JuliaReach-Powertrain_DTN05"))

sol_DTN01 = nothing
sol_DTN02 = nothing
sol_DTN03 = nothing
sol_DTN04 = nothing
sol_DTN05 = nothing
sol_DTN06 = nothing

GC.gc()


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    CylinderSurface(radius, bottom, top)

A circular cylinder surface embedded in R³ with given `radius`,
delimited by `bottom` and `top` planes.

    CylinderSurface(radius, segment)

Alternatively, construct a right circular cylinder with given `radius`
and `segment` between origin of `bottom` and `top` planes.

    CylinderSurface(radius)

Finally, construct a right vertical circular cylinder with given `radius`.

See https://en.wikipedia.org/wiki/Cylinder.
"""
struct CylinderSurface{T} <: Primitive{3,T}
  radius::T
  bot::Plane{T}
  top::Plane{T}
end

function CylinderSurface(radius::T, segment::Segment{3,T}) where {T}
  a, b = extrema(segment)
  v    = b - a
  bot  = Plane(a, v)
  top  = Plane(b, v)
  CylinderSurface(radius, bot, top)
end

function CylinderSurface(radius::T) where {T}
  _0 = (T(0), T(0), T(0))
  _1 = (T(0), T(0), T(1))
  segment = Segment(_0, _1)
  CylinderSurface(radius, segment)
end

paramdim(::Type{<:CylinderSurface}) = 2

isconvex(::Type{<:CylinderSurface}) = true

radius(c::CylinderSurface) = c.radius

axis(c::CylinderSurface) = Line(origin(c.bot), origin(c.top))

planes(c::CylinderSurface) = (c.bot, c.top)

function isright(c::CylinderSurface{T}) where {T}
  # cylinder is right if axis
  # is aligned with plane normals
  a = axis(c)
  d = a(T(1)) - a(T(0))
  v = normal(c.bot)
  w = normal(c.top)
  isparallelv = isapprox(norm(d × v), zero(T), atol=atol(T))
  isparallelw = isapprox(norm(d × w), zero(T), atol=atol(T))
  isparallelv && isparallelw
end

boundary(::CylinderSurface) = nothing


text:
mutable struct SquareLattice
    L::Int
    point_list
    coordinates_next

    # initialize an arbitrary SquareLattice with a given square length
    SquareLattice(L::Int) = (
        lat=new();
        lat.L = L;
        lat.point_list = Array{Any, 2}(undef, L, L);
        lat.coordinates_next = [mod(i+1, 1:L) for i in 1:L];
        lat
    )
    # initialize a specific SquareLattice with given data of specific form
    # here we only consider the S = 1/2, S_z in {1, -1}
    function SquareLattice(twoD_list :: Array)
        size_of_list = size(twoD_list)
        @assert size_of_list[1] == size_of_list[2]
        lat=new()
        L = size_of_list[1]
        lat.L = L
        lat.point_list = twoD_list
        lat.coordinates_next = [mod(i+1, 1:L) for i in 1:L]
        # print("ALL spins one half? ",all_spin_one_half(lat.point_list))
        @assert all_spin_one_half(lat.point_list)
        lat
    end
end

function SumOver()
    throw("NotImplemented")
end

function print_lat(lat:: SquareLattice)
    print("Square Lattice with length ", lat.L, "\nPoint List: ", lat.point_list, "\n")
end

function set_spin!(lat::SquareLattice, position::Tuple{Int, Int}, value::Int)
    @assert (value==1 || value==-1)
    lat.point_list[position...] = value
end

function get_spin(lat::SquareLattice, position::Tuple{Int, Int})
    lat.point_list[position...]
end

function get_spins(lat::SquareLattice, pos_list::Array{Tuple{Int, Int}, 1})
    spin_list = Array{Int, 1}(undef, length(pos_list))
    for i in 1:length(spin_list)
        spin_list[i] = get_spin(lat, pos_list[i])
    end
    spin_list
end

function flip!(lat::SquareLattice, position::Tuple{Int, Int})
    @assert length(position)==2
    lat.point_list[position...] *= -1
end


text:
#
# Grid data access
# ================
#
# You may want to use the this package to obtain grid data for your application.
# The [ExtendableGrids](https://github.com/j-fu/ExtendableGrids.jl) package provides a lightweight grid data structure which
# in its core is a `Dict{Type,Any}` and the possibility to [dispatch the return type on the key](https://j-fu.github.io/ExtendableGrids.jl/stable/tdict/).
#
function extract_2d()
    builder=SimplexGridBuilder(Generator=Triangulate)

    p1=point!(builder,0,0)
    p2=point!(builder,1,0)
    p3=point!(builder,1,1)
    p4=point!(builder,0,1)

    facetregion!(builder,1)
    facet!(builder,p1,p2)
    facetregion!(builder,2)
    facet!(builder,p2,p3)
    facetregion!(builder,3)
    facet!(builder,p3,p4)
    facetregion!(builder,4)
    facet!(builder,p4,p1)

    grid=simplexgrid(builder,maxvolume=0.25)

    @show grid[Coordinates]
    @show grid[CellNodes]
    @show grid[CellRegions]
    @show grid[BFaceNodes]
    @show grid[BFaceRegions]
    grid
end
# The output of this call is this:
# ```julia
# grid[Coordinates] = [0.0 1.0 1.0 0.0 0.5; 0.0 0.0 1.0 1.0 0.5]
# grid[CellNodes] = Int32[2 4 5 1; 3 1 3 2; 5 5 4 5]
# grid[CellRegions] = Int32[1, 1, 1, 1]
# grid[BFaceNodes] = Int32[2 3 4 1; 1 2 3 4]
# grid[BFaceRegions] = Int32[1, 2, 3, 4]
# ```
# Thus, the grid is described by five arrays:
# - `2 x npoints` coordinates,
# - `3 x ntriangles` connectivity,
# - `ntriangles` region markers,
# - `2 x nfacets` boundary faces,
# - `nfacets` boundary face markers.


text:
struct LoLa{A<:AbstractMatrix,T<:Tuple,Fs<:Tuple}
    w_E::A
    w_ds::T
    w_d_reducers::Fs
end

function LoLa(n::Int, w_d_reducers::Tuple; init=Flux.glorot_uniform)
    return LoLa(init(n, n), ntuple(_ -> init(n, n), length(w_d_reducers)), w_d_reducers)
end

function Flux.update!(opt, l::LoLa, dl)
    if haskey(dl, :w_E)
        Flux.update!(opt, l.w_E, dl[:w_E])
    end
    if haskey(dl, :w_ds)
        for i in eachindex(dl[:w_ds])
            Flux.update!(opt, l.w_ds[i], dl[:w_ds][i])
        end
    end
    return l
end

Flux.functor(l::LoLa) = (l.w_E, l.w_ds), x -> LoLa(x[1], x[2], l.w_d_reducers)
Flux.params!(p::Zygote.Params, l::LoLa, seen=IdSet()) = push!(p, l)

m²(k) = k[4]^2 - k[1]^2 - k[2]^2 - k[3]^2
p_T(k) = hypot(k[1], k[2])
E(k) = k[4]

using Compat

include("lola_kernel.jl")

E!(_E::AbstractVector, w, k::AbstractVector) = @tullio _E[i] = w[i, j] * E(k[j])
E!(_E::AbstractMatrix, w, k::AbstractMatrix) = @tullio _E[i, l] = w[i, j] * E(k[j, l])

slice(res, i) = view(res, i, Base.tail(axes(res))...)

function _lola3(l, k)
    T = eltype(eltype(k))
    res = similar(k, T, 3 + length(l.w_ds), axes(k)...)
    map!(m², slice(res, 1), k)
    map!(p_T, slice(res, 2), k)
    E!(slice(res, 3), l.w_E, k)
    _k = reinterpret(reshape, T, k)
    return res, _k
end

function (l::LoLa)(k)
    res, _k = _lola3(l, k)
    for i in 1:length(l.w_ds)
        wd!(slice(res, 3 + i), l.w_ds[i], _k, l.w_d_reducers[i])
    end
    return res
end

function dw_E(Δ, k)
    if ndims(k) == 1
        @tullio dw_E[i, j] := Δ[3, i] * E(k[j])
    else
        @tullio dw_E[i, l] := Δ[3, i, j] * E(k[l, j])
    end
end

function _dE(w_E, Δ)
    if ndims(Δ) == 2
        @tullio dE[i] := w_E[j, i] * Δ[3, j]
    else
        @tullio dE[i, l] := w_E[j, i] * Δ[3, j, l]
    end
end

function dk(l, k, Δ, Ω, pullbacks_wd)
    dE = _dE(l.w_E, Δ)
    dk = similar(k)
    map!(dk, k, slice(Δ, 1), slice(Δ, 2), slice(Ω, 2), dE) do k_i, Δ1, Δ2, Ω2, dE
        2Δ1 * SA[-k_i[1], -k_i[2], -k_i[3], k_i[4]] +
            SA[Δ2 / Ω2 * k_i[1], Δ2 / Ω2 * k_i[2], 0, dE]
    end
    _dk = reinterpret(reshape, eltype(Δ), dk)
    for i in 1:length(l.w_ds)
        pullbacks_wd[i][2](_dk, slice(Δ, 3 + i))
    end
    return dk
end
function dw_ds(l, Δ, pullbacks_wd)
    ntuple(length(l.w_ds)) do i
        dw_d = reinterpret(reshape, eltype(Δ), zero(l.w_ds[i]))
        pullbacks_wd[i][1](dw_d, slice(Δ, 3 + i))
        return dw_d
    end
end

function ChainRulesCore.rrule(l::LoLa, k)
    Ω, _k = _lola3(l, k)
    pullbacks_wd = ntuple(length(l.w_ds)) do i
        _, pb_w, pb_k = wd_adjoint!(slice(Ω, 3 + i), l.w_ds[i], _k, l.w_d_reducers[i])
        return pb_w, pb_k
    end

    function lola_pullback(Δ)
        return (
            Composite{typeof(l)}(; w_E=dw_E(Δ, k), w_ds=dw_ds(l, Δ, pullbacks_wd)),
            dk(l, k, Δ, Ω, pullbacks_wd),
        )
    end
    return Ω, lola_pullback
end


text:
# Mu' = f
type ODEProblem{uType,tType,isinplace,F,C,MM} <: AbstractODEProblem{uType,tType,isinplace}
  f::F
  u0::uType
  tspan::Tuple{tType,tType}
  callback::C
  mass_matrix::MM
end

function ODEProblem(f,u0,tspan; iip = isinplace(f,3),callback=nothing,mass_matrix=I)
  ODEProblem{typeof(u0),eltype(tspan),iip,typeof(f),typeof(callback),typeof(mass_matrix)}(f,u0,tspan,callback,mass_matrix)
end

# Mu' = f[1] + f[2] + ...
type SplitODEProblem{uType,tType,isinplace,F,C,MM} <: AbstractODEProblem{uType,tType,isinplace}
  f::F
  u0::uType
  tspan::Tuple{tType,tType}
  callback::C
  mass_matrix::MM
end

function SplitODEProblem(f,u0,tspan; iip = isinplace(f[2],3),callback=nothing,mass_matrix=I)
  @assert typeof(f) <: Tuple
  SplitODEProblem{typeof(u0),eltype(tspan),iip,typeof(f),typeof(callback),typeof(mass_matrix)}(f,u0,tspan,callback,mass_matrix)
end

# M[i]*u[i]' = f[i]
type PartitionedODEProblem{uType,tType,isinplace,F,C,MM} <: AbstractODEProblem{uType,tType,isinplace}
  f::F
  u0::uType
  tspan::Tuple{tType,tType}
  callback::C
  mass_matrix::MM
end

function PartitionedODEProblem(f,u0,tspan; iip = isinplace(f[1],3),callback=nothing,mass_matrix=nothing)
  @assert typeof(f) <: Tuple
  @assert typeof(u0) <: Tuple
  if mass_matrix == nothing
    _mm = ((I for i in 1:length(f))...)
  end
  PartitionedODEProblem{typeof(u0),eltype(tspan),iip,typeof(f),typeof(callback),typeof(_mm)}(f,u0,tspan,callback,_mm)
end

# u'' = f(t,u,du,ddu)
function SecondOrderODEProblem(f,u0,du0,tspan; iip = isinplace(f,4),callback=nothing,mass_matrix=I)
  if iip
    f1 = function (t,u,v,du)
      du .= v
    end
  else
    f1 = function (t,u,v)
      v
    end
  end
  _f = (f1,f)
  _u0 = (u0,du0)
  _mass_matrix = (mass_matrix,I)
  PartitionedODEProblem{typeof(_u0),eltype(tspan),iip,typeof(_f),typeof(callback),typeof(_mass_matrix)}(_f,_u0,tspan,callback,_mass_matrix)
end


text:
"""
    GeoMeantoRelEntrBridge{T}

The `geometric mean cone` is representable with a relative entropy constraint and a
nonnegative auxiliary variable, since ``u \\le \\prod_{i=1}^n w_i^{1/n}`` is equivalent
to ``y \\ge 0`` and ``0 \\le u + y \\le \\prod_{i=1}^n w_i^{1/n}``, and the latter
inequality is equivalent to ``1 \\le \\prod_{i=1}^n (\\frac{w_i}{u + y})^{1/n}``, which
is equivalent to ``0 \\le \\sum_{i=1}^n \\log (\\frac{w_i}{u + y})^{1/n}``, which is
equivalent to ``0 \\ge \\sum_{i=1}^n (u + y) \\log (\\frac{u + y}{w_i})``. Thus
``(u, w) \\in GeometricMeanCone(1 + n)`` is representable as ``y \\ge 0``,
``(0, w, (u + y) e) \\in RelativeEntropyCone(1 + 2n)``, where ``e`` is a vector of ones.
"""
struct GeoMeantoRelEntrBridge{T,F,G,H} <: AbstractBridge
    y::MOI.VariableIndex
    nn_index::CI{F,MOI.Nonnegatives} # for y >= 0
    relentr_index::CI{G,MOI.RelativeEntropyCone}
end
function bridge_constraint(
    ::Type{GeoMeantoRelEntrBridge{T,F,G,H}},
    model::MOI.ModelLike,
    f::H,
    s::MOI.GeometricMeanCone,
) where {T,F,G,H}
    f_scalars = MOIU.eachscalar(f)
    (y, nn_index) = MOI.add_constrained_variables(model, MOI.Nonnegatives(1))
    w_func = MOIU.vectorize(fill(
        MOIU.operate(+, T, f_scalars[1], MOI.SingleVariable(y[1])),
        MOI.dimension(s) - 1,
    ))
    relentr_func = MOIU.operate(
        vcat,
        T,
        zero(MOI.ScalarAffineFunction{Float64}),
        f_scalars[2:end],
        w_func,
    )
    relentr_index = MOI.add_constraint(
        model,
        relentr_func,
        MOI.RelativeEntropyCone(MOI.output_dimension(relentr_func)),
    )
    return GeoMeantoRelEntrBridge{T,F,G,H}(y[1], nn_index, relentr_index)
end

function MOI.supports_constraint(
    ::Type{<:GeoMeantoRelEntrBridge{T}},
    ::Type{<:MOI.AbstractVectorFunction},
    ::Type{MOI.GeometricMeanCone},
) where {T}
    return true
end
function MOIB.added_constrained_variable_types(::Type{<:GeoMeantoRelEntrBridge})
    return [(MOI.Nonnegatives,)]
end
function MOIB.added_constraint_types(
    ::Type{<:GeoMeantoRelEntrBridge{T,F,G}},
) where {T,F,G}
    return [(G, MOI.RelativeEntropyCone)]
end
function concrete_bridge_type(
    ::Type{<:GeoMeantoRelEntrBridge{T}},
    H::Type{<:MOI.AbstractVectorFunction},
    ::Type{MOI.GeometricMeanCone},
) where {T}
    F = MOI.VectorOfVariables
    S = MOIU.scalar_type(H)
    G = MOIU.promote_operation(
        vcat,
        T,
        T,
        S,
        MOIU.promote_operation(+, T, S, MOI.SingleVariable),
    )
    return GeoMeantoRelEntrBridge{T,F,G,H}
end

# Attributes, Bridge acting as a model
MOI.get(bridge::GeoMeantoRelEntrBridge, ::MOI.NumberOfVariables) = 1
function MOI.get(bridge::GeoMeantoRelEntrBridge, ::MOI.ListOfVariableIndices)
    return [bridge.y]
end
function MOI.get(
    bridge::GeoMeantoRelEntrBridge{T,F},
    ::MOI.NumberOfConstraints{F,MOI.Nonnegatives},
) where {T,F}
    return 1
end
function MOI.get(
    bridge::GeoMeantoRelEntrBridge{T,F,G},
    ::MOI.NumberOfConstraints{G,MOI.RelativeEntropyCone},
) where {T,F,G}
    return 1
end
function MOI.get(
    bridge::GeoMeantoRelEntrBridge{T,F},
    ::MOI.ListOfConstraintIndices{F,MOI.Nonnegatives},
) where {T,F}
    return [bridge.nn_index]
end
function MOI.get(
    bridge::GeoMeantoRelEntrBridge{T,F,G},
    ::MOI.ListOfConstraintIndices{G,MOI.RelativeEntropyCone},
) where {T,F,G}
    return [bridge.relentr_index]
end

# References
function MOI.delete(model::MOI.ModelLike, bridge::GeoMeantoRelEntrBridge)
    MOI.delete(model, bridge.relentr_index)
    MOI.delete(model, bridge.nn_index)
    return MOI.delete(model, bridge.y)
end

# Attributes, Bridge acting as a constraint
function MOI.get(
    model::MOI.ModelLike,
    ::MOI.ConstraintFunction,
    bridge::GeoMeantoRelEntrBridge{T,F,G,H},
) where {T,F,G,H}
    relentr_func = MOIU.eachscalar(MOI.get(
        model,
        MOI.ConstraintFunction(),
        bridge.relentr_index,
    ))
    d = div(length(relentr_func) - 1, 2)
    u_func = MOIU.remove_variable(
        MOIU.operate(-, T, relentr_func[end], MOI.SingleVariable(bridge.y)),
        bridge.y,
    )
    w_func = relentr_func[2:(1+d)]
    return MOIU.convert_approx(H, MOIU.operate(vcat, T, u_func, w_func))
end
function MOI.get(
    model::MOI.ModelLike,
    ::MOI.ConstraintSet,
    bridge::GeoMeantoRelEntrBridge,
)
    return MOI.GeometricMeanCone(
        1 + div(
            MOI.dimension(MOI.get(
                model,
                MOI.ConstraintSet(),
                bridge.relentr_index,
            )) - 1,
            2,
        ),
    )
end
function MOI.supports(
    ::MOI.ModelLike,
    ::Union{MOI.ConstraintPrimalStart,MOI.ConstraintDualStart},
    ::Type{<:GeoMeantoRelEntrBridge},
)
    return true
end
function MOI.get(
    model::MOI.ModelLike,
    attr::Union{MOI.ConstraintPrimal,MOI.ConstraintPrimalStart},
    bridge::GeoMeantoRelEntrBridge,
)
    relentr_primal = MOI.get(model, attr, bridge.relentr_index)
    d = div(length(relentr_primal) - 1, 2)
    y_val = MOI.get(model, attr, bridge.nn_index)[1]
    u_val = sum(relentr_primal[(2+d):end]) / d - y_val
    return vcat(u_val, relentr_primal[2:(1+d)])
end
function MOI.set(
    model::MOI.ModelLike,
    attr::MOI.ConstraintPrimalStart,
    bridge::GeoMeantoRelEntrBridge,
    value,
)
    u_val = value[1]
    u_pos = -min(zero(u_val), u_val)
    d = length(value) - 1
    MOI.set(model, MOI.VariablePrimalStart(), bridge.y, u_pos)
    MOI.set(model, attr, bridge.nn_index, [u_pos])
    MOI.set(
        model,
        attr,
        bridge.relentr_index,
        vcat(0, value[2:end], fill(u_pos + u_val, d)),
    )
    return
end
# Given a is dual on y >= 0 and (b, c, d) is dual on RelativeEntropyCone constraint,
# dual on (u, w) in GeometricMeanCone is (-a, c). Note that sum(d) = -a, so we could
# instead use (sum(d), c).
function MOI.get(
    model::MOI.ModelLike,
    attr::Union{MOI.ConstraintDual,MOI.ConstraintDualStart},
    bridge::GeoMeantoRelEntrBridge,
)
    u_dual = -MOI.get(model, attr, bridge.nn_index)[1]
    relentr_dual = MOI.get(model, attr, bridge.relentr_index)
    d = div(length(relentr_dual) - 1, 2)
    w_dual = relentr_dual[2:(d+1)]
    return vcat(u_dual, w_dual)
end
# Given GeometricMeanCone constraint dual start of (u, w), constraint dual on y >= 0 is -u
# and on RelativeEntropyCone constraint is (-u/n, w, u/n * (log.(w/geomean(w)) .+ 1)).
# Note log.(w/geomean(w)) = log.(w) .- sum(log, w) / n.
function MOI.set(
    model::MOI.ModelLike,
    ::MOI.ConstraintDualStart,
    bridge::GeoMeantoRelEntrBridge,
    value,
)
    d = length(value) - 1
    u = value[1]
    w = value[2:end]
    MOI.set(model, MOI.ConstraintDualStart(), bridge.nn_index, [-u])
    relentr_dual =
        vcat(-u / d, value[2:end], u / d * (log.(w) .+ (1 - sum(log, w) / d)))
    MOI.set(
        model,
        MOI.ConstraintDualStart(),
        bridge.relentr_index,
        relentr_dual,
    )
    return
end


text:
Move = Pair{CartesianIndex{2}, CartesianIndex{2}}

"""
    initialize_board()

Return an empty board.
"""
function initialize_board()
    board = zeros(Int, 5, 5)
    return board
end

"""
    _robot_possible_moves(board::Matrix{Int}, position::CartesianIndex{2})

Return all the possible moves that can be done on `board` from a robot starting
at `position`.
"""
function _robot_possible_moves(board::Matrix{Int}, position::CartesianIndex{2})
    @assert board[position] in [1,2,3,4,5,6]
    moves = Vector{Pair{CartesianIndex{2}, CartesianIndex{2}}}()
    for (i,j) in [(-1,0),(0,1),(1,0),(0,-1)]
        ending_position = CartesianIndex(position[1]+i, position[2]+j)
        while (ending_position[1] >= 1 && ending_position[1] <= 5) && (ending_position[2] >= 1 && ending_position[2] <= 5) && board[ending_position] == 0
            ending_position = CartesianIndex(ending_position[1]+i, ending_position[2]+j)
        end
        if ending_position[1] != 0 && ending_position[1] != 6 && ending_position[2] != 0 && ending_position[2] != 6 && board[ending_position] != 0
            ending_position = CartesianIndex(ending_position[1]-i, ending_position[2]-j)
            if position != ending_position
                push!(moves, position => ending_position)
            end
        end
    end
    return moves
end

"""
    get_red_position(board::Matrix{Int})

Return the position of the Red robot.
"""
function get_red_position(board::Matrix{Int})
    return findfirst(x->x==1, board)
end

"""
    get_robot_positions(board::Matrix{Int})

Return a vector of positions with the positions of all the robots.
"""
function get_robot_positions(board::Matrix{Int})
    return findall(x->x!=0, board)
end

"""
    get_possible_moves(board::Matrix{Int})

Return a vector of all possible moves that can be made on the `board`.
"""
function get_possible_moves(board::Matrix{Int})
    moves = Vector{Pair{CartesianIndex{2}, CartesianIndex{2}}}()
    for position in get_robot_positions(board)
        append!(moves, _robot_possible_moves(board, position))
    end
    return moves
end

"""
    is_finished(board::Matrix{Int})

Check if Red is in the center of the `board`.
"""
function is_finished(board::Matrix{Int})
    if 1 in board[3,3]
        return true
    end
    return false
end

"""
    move_robot!(board::Matrix{Int}, starting_position::CartesianIndex{2}, ending_position::CartesianIndex{2})
    move_robot!(board::Matrix{Int}, move::Pair{CartesianIndex{2}, CartesianIndex{2}})

Move the robot in `starting_position` or `move[1]` to the `ending_position` or
`move[2]`.
"""
function move_robot!(board::Matrix{Int}, starting_position::CartesianIndex{2}, ending_position::CartesianIndex{2})
    @assert board[starting_position] != 0
    @assert board[ending_position] == 0
    board[ending_position] = board[starting_position]
    board[starting_position] = 0
    return board
end
move_robot!(board::Matrix{Int}, move::Pair{CartesianIndex{2}, CartesianIndex{2}}) =
    move_robot!(board, move[1], move[2])

"""
    print_current_board(board::Matrix{Int})

Print the current `board` in a simplified format where `X` corresponds to
robots, `O` to the Red one and `+` to the center of the board.
"""
function print_current_board(board::Matrix{Int})
    for i in 1:5
        for j in 1:5
            if i == 3 && j == 3 && board[i,j] == 0
                print("+")
            elseif board[i,j] == 0
                print(" ")
            elseif board[i,j] == 1
                print("O")
            else
                print("X")
            end
        end
        print("\n")
    end
    return nothing
end

max_stars_repo_name: remi-garcia/LunarLockout

max_issues_repo_name: remi-garcia/LunarLockout

max_forks_repo_name: remi-garcia/LunarLockout


text:
using LightGraphs
using Test

include("../../src/utils/maximalcliques.jl")

@testset "utils/maximalcliques" begin
    @testset "refine/refine_set" begin
        settest = [1, 2, 3, 4]
        @test refine(settest, 1) == [[1], [2, 3, 4]]
        @test refine(settest, 4) == [[4], [1, 2, 3]]
        @test refine(settest, 5) == [settest]
        @test refine(settest, settest) == [settest]
        settest = [[1, 2, 3], [2, 3, 4], [5, 6]]
        @test refineset(settest, 2) == [[2], [1, 3], [2], [3, 4], [5, 6]]
        @test refineset(settest, 7) == settest
    end
    @testset "lbfs/perfect elimination ordering" begin
        graph = path_graph(10)
        add_edge!(graph, 1, 3)
        add_edge!(graph, 1, 9)
        add_edge!(graph, 3, 5)
        add_edge!(graph, 5, 7)
        add_edge!(graph, 7, 9)
        add_edge!(graph, 3, 9)
        add_edge!(graph, 3, 7)
        @test compute_lbfs(graph) == [1, 2, 3, 9, 7, 5, 4, 8, 10, 6]
        @test compute_perfect_elimination_ordering(graph) == [6, 10, 8, 4, 5, 7, 9, 3, 2, 1]
    end
    @testset "isasubset" begin
        sets = [[1, 2, 3], [4, 5, 6]]
        @test isasubset([1, 2], sets)
        @test !isasubset([3, 4], sets)
    end
    @testset "get prior/later/closest neighbors" begin
        graph = path_graph(10)
        add_edge!(graph, 1, 3)
        add_edge!(graph, 1, 9)
        add_edge!(graph, 3, 5)
        add_edge!(graph, 5, 7)
        add_edge!(graph, 7, 9)
        add_edge!(graph, 3, 9)
        add_edge!(graph, 3, 7)
        ordering = [1, 2, 3, 9, 7, 5, 4, 8, 10, 6]
        @test get_prior_neighbors(graph, ordering, 4) == [1, 3]
        @test get_prior_neighbors(graph, ordering, 4, [3]) == [1] # Here ordering[3] = 3 and ordering[1] = 1, the exclude parameter takes ***indexes***.
        @test get_prior_neighbors(graph, ordering, 4, [3, 1]) == []
        @test get_prior_neighbors(graph, ordering, 1) == []
        @test get_later_neighbors(graph, ordering, 4) == [7, 8, 10]
        @test get_later_neighbors(graph, ordering, 4, [5]) == [8, 10]
        @test get_later_neighbors(graph, ordering, 4, [5, 8, 9]) == []
        @test get_later_neighbors(graph, ordering, length(ordering)) == []
        @test get_closest_neighbor(graph, ordering, 1) == nothing
        @test get_closest_neighbor(graph, ordering, 4) == 3
    end
    @testset "is_complete" begin
        graph = path_graph(4)
        add_edge!(graph, 1, 4)
        add_edge!(graph, 1, 3)
        @test !is_complete(graph)
        @test is_complete(graph, 1:3)
        add_edge!(graph, 4, 2)
        @test is_complete(graph)
    end
    @testset "is_maximal_clique" begin
        graph = path_graph(4)
        add_edge!(graph, 1, 4)
        add_edge!(graph, 1, 3)
        @test is_maximal_clique(graph, [1, 2, 3])
        @test !is_maximal_clique(graph, [1, 2, 3], [[1, 2, 3, 4]]) # Not a real scenario on this graph
        @test !is_maximal_clique(graph, [1, 2, 3, 4])
    end
    @testset "get_maximal_clique" begin
        graph = path_graph(4)
        add_edge!(graph, 1, 4)
        add_edge!(graph, 1, 3)
        cliques = get_maximal_cliques(graph)
        @test isasubset([1, 2, 3], cliques)
        @test isasubset([1, 3, 4], cliques)
    end
    @testset "get_cliquetree" begin
        graph = path_graph(4)
        add_edge!(graph, 1, 4)
        add_edge!(graph, 1, 3)
        cliques = get_maximal_cliques(graph)
        cliquetree, nb_lc = get_cliquetree(cliques)
        @test cliquetree[1] == [1, 2]
    end
    @testset "is_chordal" begin
        graph = path_graph(4)
        add_edge!(graph, 1, 4)
        @test !is_chordal(graph)
        add_edge!(graph, 1, 3)
        @test is_chordal(graph)
    end
end

max_stars_repo_name: charlyalizadeh/OPFdecgen

max_issues_repo_name: charlyalizadeh/OPFdecgen

max_forks_repo_name: charlyalizadeh/OPFdecgen


text:
using SparseArrays, Distributions
import LambdaTAME: spatial_network, duplication_perturbation_noise_model, netalignmr, knearest_sparsification, successive_netalignmr, successive_netalignmr_profiled
import LambdaTAME: local_search_profiled, local_search
@testset "Post Processing" begin

    m=10
    steps=10
    p=.5

    A = spatial_network(m,2;degreedist=LogNormal(log(5),1))
    n = m + steps
    B,_ = duplication_perturbation_noise_model(A,steps, p)

    #perm = shuffle(1:size(B,1))#
    #B = B[perm,perm]
    d = 10
    U = rand(size(A,1),d)
    V = rand(size(B,1),d)
    matching = collect(enumerate(1:size(A,1)))
    L = sparse(1:min(n,m),1:min(n,m),1,n,m)
    L += sprand(size(L)...,.1) # add some fill in

    @testset "Klau Algorithm" begin
        @suppress_out begin
            @inferred netalignmr(B,A,V,U,5,KlauAlgo())
            @inferred netalignmr(B,A,V,U,matching,5,KlauAlgo())
            @inferred netalignmr(B,A,V,U,KlauAlgo())
            @inferred netalignmr(B,A,V,U,matching,KlauAlgo())
            @inferred netalignmr(B,A,L,KlauAlgo())
        end
    end

    @testset "Successive Klau Algorithm" begin
        @suppress_out begin
            @inferred successive_netalignmr(B,A,V,U,matching,SuccessiveKlauAlgo(k=5,iterDelta=10,successive_iter=3))
            @inferred successive_netalignmr_profiled(B,A,V,U,matching,SuccessiveKlauAlgo(k=5,iterDelta=10,successive_iter=3))
        end
    end

    @testset "k nearest sparsification" begin
        @inferred knearest_sparsification(V,U,matching,5)
    end

    @testset "Local Search" begin

        A_ten = DistributedTensorConstruction.tensor_from_graph(A,3,Clique())
        B_ten = DistributedTensorConstruction.tensor_from_graph(B,3,Clique())

        k = 15
        @suppress_out begin
            @inferred local_search(A,B,A_ten,B_ten,U,V,Dict(matching),d)
            @inferred local_search_profiled(A,B,A_ten,B_ten,U,V,Dict(matching),d)
            @inferred local_search(A,B,A_ten,B_ten,U,V,Dict(matching),LocalSearch())
            @inferred local_search_profiled(A,B,A_ten,B_ten,U,V,Dict(matching),LocalSearch())
        end

    end

end

max_stars_repo_name: chuckcol/LambdaTAME

max_issues_repo_name: chuckcol/LambdaTAME

max_forks_repo_name: chuckcol/LambdaTAME


text:
module rigid_box_examples
using FinEtools
using FinEtoolsAcoustics
import Arpack: eigs

function rigid_box_Q4_example()
    println("""

    Example from Boundary element acoustics: Fundamentals and computer codes, TW Wu, page 123.
    Internal resonance problem. Reference frequencies: 90.7895, 181.579, 215.625, 233.959, 272.368, 281.895
    Quadrilateral mesh.
    """)

    t0 = time()

    rho=1.21*1e-9;# mass density
    c =345.0*1000;# millimeters per second
    bulk= c^2*rho;
    Lx=1900.0;# length of the box, millimeters
    Ly=800.0; # length of the box, millimeters
    n=14;#
    neigvs=18;
    OmegaShift=10.0;

    fens,fes = Q4block(Lx,Ly,n,n); # Mesh

    geom = NodalField(fens.xyz)
    P = NodalField(zeros(size(fens.xyz,1),1))

    numberdofs!(P)

    femm = FEMMAcoust(IntegDomain(fes, GaussRule(2, 2)),
    MatAcoustFluid(bulk,rho))

    S = acousticstiffness(femm, geom, P);
    C = acousticmass(femm, geom, P);

    d,v,nev,nconv =eigs(C+OmegaShift*S, S; nev=neigvs, which=:SM)
    d = d .- OmegaShift;
    fs=real(sqrt.(complex(d)))./(2*pi)
    println("Eigenvalues: $fs [Hz]")

    println("Total time elapsed = ",time() - t0,"s")

    File =  "rigid_box.vtk"
    scalarllist = Any[]
    for n  = 1:15
        push!(scalarllist, ("Pressure_mode_$n", v[:,n]));
    end
    vtkexportmesh(File, connasarray(fes), geom.values, FinEtools.MeshExportModule.VTK.Q4;
    scalars=scalarllist)
    @async run(`"paraview.exe" $File`)

    true

end # rigid_box_Q4_example

function allrun()
    println("#####################################################")
    println("# rigid_box_Q4_example ")
    rigid_box_Q4_example()
end # function allrun

end # module rigid_box_examples


text:
using TextUnidecode, BenchmarkTools
const SUITE = BenchmarkGroup()
SUITE["南无阿弥陀佛"] = @benchmarkable unidecode("南无阿弥陀佛")
SUITE["really only ascii all fine!"] = @benchmarkable unidecode("really only ascii all fine!")


text:
using Test
import Thermodynamics
const TD = Thermodynamics
const TP = TD.Parameters
const TDTP = TD.TemperatureProfiles
using ForwardDiff

import CLIMAParameters
const CP = CLIMAParameters

function get_parameter_set(::Type{FT}) where {FT}
    toml_dict = CP.create_toml_dict(FT; dict_type = "alias")
    aliases = string.(fieldnames(TP.ThermodynamicsParameters))
    param_pairs = CP.get_parameter_values!(toml_dict, aliases, "Thermodynamics")
    param_set = TP.ThermodynamicsParameters{FT}(; param_pairs...)
    logfilepath = joinpath(@__DIR__, "logfilepath_$FT.toml")
    CP.log_parameter_information(toml_dict, logfilepath)
    return param_set
end

const param_set_Float64 = get_parameter_set(Float64)
const param_set_Float32 = get_parameter_set(Float32)
parameter_set(::Type{Float64}) = param_set_Float64
parameter_set(::Type{Float32}) = param_set_Float32

@testset "TemperatureProfiles - DecayingTemperatureProfile" begin
    for FT in [Float32, Float64]
        param_set = parameter_set(FT)
        _grav = FT(TP.grav(param_set))
        _R_d = FT(TP.R_d(param_set))
        _MSLP = FT(TP.MSLP(param_set))

        z = collect(range(FT(0), stop = FT(25e3), length = 100))
        n = 7
        _T_virt_surf_range = collect(range(FT(260), stop = FT(320), length = n))
        _T_min_ref_range = collect(range(FT(170), stop = FT(230), length = n))
        _H_t_range = collect(range(FT(6e3), stop = FT(13e3), length = n))

        profiles = []
        for (_T_virt_surf, _T_min_ref, _H_t) in
            zip(_T_virt_surf_range, _T_min_ref_range, _H_t_range)
            profiles = [
                TDTP.DecayingTemperatureProfile{FT}(
                    param_set,
                    _T_virt_surf,
                    _T_min_ref,
                    _H_t,
                ),
                TDTP.DryAdiabaticProfile{FT}(
                    param_set,
                    _T_virt_surf,
                    _T_min_ref,
                ),
                TDTP.IsothermalProfile(param_set, _T_virt_surf),
                TDTP.IsothermalProfile(param_set, FT),
            ]

            for profile in profiles
                args = profile.(Ref(param_set), z)
                T_virt = first.(args)
                p = last.(args)

                # Test that surface pressure is equal to the
                # specified boundary condition (MSLP)
                mask_z_0 = z .≈ 0
                @test all(p[mask_z_0] .≈ _MSLP)

                function log_p_over_MSLP(_z)
                    _, _p = profile(param_set, _z)
                    return log(_p / _MSLP)
                end
                ∇log_p_over_MSLP =
                    _z -> ForwardDiff.derivative(log_p_over_MSLP, _z)
                # Uses density computed from pressure derivative
                # and ideal gas law to reconstruct virtual temperature
                T_virt_rec = -_grav ./ (_R_d .* ∇log_p_over_MSLP.(z))
                @test all(T_virt_rec .≈ T_virt)
            end
        end

    end
end

rm(joinpath(@__DIR__, "logfilepath_Float32.toml"); force = true)
rm(joinpath(@__DIR__, "logfilepath_Float64.toml"); force = true)


text:
export iid

iid(n) = dist -> (dist ^ n)
iid(n::Int...) = dist -> (dist ^ n)


text:
#using Test
using EquationEditor
#import EquationEditor: LogLevel, Debug, Info, Wslog, with_logger
#import EquationEditor.WebSocketLogger
const EE = EquationEditor
#with_logger(WebSocketLogger(stderr, Debug)) do
#    @info "Info"
#    @debug "Debug"
#end
const server, task = serveEE()
EE.open_browser("chrome", url = "localhost:8000/interact.html")

# Non recursive
Revise.includet(joinpath(@__DIR__, "track_this.jl"))
# Recursive:
#Revise.track

codelogger = Revise.debug_logger(; min_level = Base.CoreLogging.Debug)

Revise.FileModules
Revise.MethodSummary

lastlog = logger.logs[end]

fieldnames(typeof(lastlog))
println(lastlog.file)

lastlog.file

lastkwargs = Dict(lastlog.kwargs)

lastmod, lastrelocatableexpr = get(lastkwargs, :deltainfo, "")

# Create a page for lastmod if OK, store it as a file, serve it, include metadata for changeable expression positions, also update over websocket.

Revise.actions(logger, line = true)

# Add these to watched files, make a proprietary file format similar to MathML, MathCad or Mathjax?
Revise.watched_files

max_stars_repo_name: hustf/EquationEditor

max_issues_repo_name: hustf/EquationEditor

max_forks_repo_name: hustf/EquationEditor


text:
module Legacy

using Compat; import Compat.String

import ..Utilities: isexpr
import ..Collector
import ..Cache
import ..Formats

include("atdoc.jl")
include("macros.jl")
include("types.jl")

end


text:
###### type definition ##########

import Base: convert, length, show, getindex, start, next, done, isempty, endof,
             size

abstract type AbstractTimeSeries end

struct TimeArray{T, N, D <: TimeType, A <: AbstractArray{T, N}} <: AbstractTimeSeries

    timestamp::Vector{D}
    values::A
    colnames::Vector{String}
    meta::Any

    function TimeArray{T, N, D, A}(
            timestamp::AbstractVector{D},
            values::A,
            colnames::Vector{String},
            meta::Any) where {T, N, D <: TimeType, A <: AbstractArray{T, N}}
        nrow, ncol = size(values, 1, 2)

        if nrow != length(timestamp)
            throw(DimensionMismatch("values must match length of timestamp"))
        elseif ncol != length(colnames)
            throw(DimensionMismatch("column names must match width of array"))
        elseif !allunique(timestamp)
            throw(ArgumentError("there are duplicate dates"))
        elseif !(issorted(timestamp) || issorted(timestamp, rev=true))
            throw(ArgumentError("dates are mangled"))
        elseif issorted(timestamp, rev=true)
            new(reverse(timestamp), flipdim(values, 1),
                replace_dupes(colnames), meta)
        else
            new(timestamp, values, replace_dupes(colnames), meta)
        end
    end
end

###### outer constructor ########

TimeArray(d::AbstractVector{D}, v::AbstractArray{T, N}, c::Vector{S},
          m::Any) where {T, N, D <: TimeType, S <: AbstractString} =
    TimeArray{T, N, D, typeof(v)}(d, v, map(String, c), m)
TimeArray(d::D, v::AbstractArray{T, N}, c::Vector{S},
          m::Any) where {T, N, D <: TimeType, S <: AbstractString} =
    TimeArray{T, N, D, typeof(v)}([d], v, map(String, c), m)

# when no column names are provided - meta is forced to nothing
TimeArray(d::AbstractVector{D}, v::AbstractArray) where {D <: TimeType} =
    TimeArray(d, v, fill("", size(v, 2)), nothing)
TimeArray(d::D, v::AbstractArray) where {D <: TimeType} =
    TimeArray([d], v, fill("", size(v, 2)), nothing)

# when no meta is provided
TimeArray(d::AbstractVector{D}, v::AbstractArray, c) where {D <: TimeType} =
    TimeArray(d, v, c, nothing)
TimeArray(d::D, v::AbstractArray, c) where {D <: TimeType} =
    TimeArray([d], v, c, nothing)

###### conversion ###############

convert(::Type{TimeArray{Float64, 1}}, x::TimeArray{Bool, 1}) =
    TimeArray(x.timestamp, map(Float64, x.values), x.colnames, x.meta)
convert(::Type{TimeArray{Float64, 2}}, x::TimeArray{Bool, 2}) =
    TimeArray(x.timestamp, map(Float64, x.values), x.colnames, x.meta)

convert(x::TimeArray{Bool, 1}) = convert(TimeArray{Float64, 1}, x::TimeArray{Bool, 1})
convert(x::TimeArray{Bool, 2}) = convert(TimeArray{Float64, 2}, x::TimeArray{Bool, 2})

###### length ###################

length(ata::AbstractTimeSeries) = length(ata.timestamp)

###### size #####################

size(ta::TimeArray, dim...) = size(ta.values, dim...)

###### iterator protocol #########

start(ta::TimeArray)   = 1
next(ta::TimeArray, i) = ((ta.timestamp[i], ta.values[i, :]), i + 1)
done(ta::TimeArray, i) = (i > length(ta))
isempty(ta::TimeArray) = (length(ta) == 0)

###### show #####################

function show(io::IO, ta::TimeArray{T}) where {T}

    # variables
    nrow          = size(ta.values, 1)
    ncol          = size(ta.values, 2)
    intcatcher    = falses(ncol)
    for c in 1:ncol
        rowcheck = @. trunc(ta.values[:, c]) - ta.values[:, c] == 0
        if sum(rowcheck) == length(rowcheck)
            intcatcher[c] = true
        end
    end
    spacetime     = nrow > 0 ? strwidth(string(ta.timestamp[1])) + 3 : 3
    firstcolwidth = strwidth(ta.colnames[1])
    colwidth      = Int[]
    for m in 1:ncol
        (T == Bool) || (nrow == 0) ?
        push!(colwidth, max(strwidth(ta.colnames[m]), 5)) :
        push!(colwidth, max(strwidth(ta.colnames[m]), strwidth(@sprintf("%.2f", maximum(ta.values[:,m]))) + DECIMALS - 2))
    end

    # summary line
    @printf(io, "%dx%d %s", nrow, ncol, typeof(ta))
    nrow > 0 && @printf(io, " %s to %s", string(ta.timestamp[1]), string(ta.timestamp[end]))
    println(io, "")
    println(io, "")

   # row label line

    print(io, ^(" ", spacetime), ta.colnames[1], ^(" ", colwidth[1] + 2 -firstcolwidth))

   for p in 2:length(colwidth)
       print(io, ta.colnames[p], ^(" ", colwidth[p] - strwidth(ta.colnames[p]) + 2))
   end
   println(io, "")

  # timestamp and values line
    if nrow > 7
        for i in 1:4
            print(io, ta.timestamp[i], " | ")
        for j in 1:ncol
            T == Bool ?
            print(io, rpad(ta.values[i,j], colwidth[j] + 2, " ")) :
            intcatcher[j] & SHOWINT ?
            print(io, rpad(round(Integer, ta.values[i,j]), colwidth[j] + 2, " ")) :
            isnan(ta.values[i,j]) ?
            print(io, rpad(MISSING, colwidth[j] + 2, " ")) :
            print(io, rpad(round(ta.values[i,j], DECIMALS), colwidth[j] + 2, " "))
        end
        println(io, "")
        end

        println(io, '\u22EE')

        for i in nrow-3:nrow
            print(io, ta.timestamp[i], " | ")
        for j in 1:ncol
            T == Bool ?
            print(io, rpad(ta.values[i,j], colwidth[j] + 2, " ")) :
            intcatcher[j] & SHOWINT ?
            print(io, rpad(round(Integer, ta.values[i,j]), colwidth[j] + 2, " ")) :
            isnan(ta.values[i,j]) ?
            print(io, rpad(MISSING, colwidth[j] + 2, " ")) :
            print(io, rpad(round(ta.values[i,j], DECIMALS), colwidth[j] + 2, " "))
        end
        println(io, "")
        end
    elseif nrow > 0
        for i in 1:nrow
            print(io, ta.timestamp[i], " | ")
        for j in 1:ncol
            T == Bool ?
            print(io, rpad(ta.values[i,j], colwidth[j] + 2, " ")) :
            intcatcher[j] & SHOWINT ?
            print(io, rpad(round(Integer, ta.values[i,j]), colwidth[j] + 2, " ")) :
            isnan(ta.values[i,j]) ?
            print(io, rpad(MISSING, colwidth[j] + 2, " ")) :
            print(io, rpad(round(ta.values[i,j], DECIMALS), colwidth[j] + 2, " "))
        end
        println(io, "")
        end
    end
end

###### getindex #################

# single row
function getindex(ta::TimeArray, n::Int)
    # avoid conversion to column vector
    TimeArray(ta.timestamp[n], ta.values[n:n, :], ta.colnames, ta.meta)
end

# single row 1d
function getindex(ta::TimeArray{T, 1}, n::Int) where {T}
    TimeArray(ta.timestamp[n], ta.values[[n]], ta.colnames, ta.meta)
end

# range of rows
function getindex(ta::TimeArray, r::UnitRange{Int})
    TimeArray(ta.timestamp[r], ta.values[r, :], ta.colnames, ta.meta)
end

# range of 1d rows
function getindex(ta::TimeArray{T, 1}, r::UnitRange{Int}) where {T}
    TimeArray(ta.timestamp[r], ta.values[r], ta.colnames, ta.meta)
end

# array of rows
function getindex(ta::TimeArray, a::AbstractVector{S}) where {S <: Integer}
    TimeArray(ta.timestamp[a], ta.values[a, :], ta.colnames, ta.meta)
end

# array of 1d rows
function getindex(ta::TimeArray{T, 1}, a::AbstractVector{S}) where {T, S <: Integer}
    TimeArray(ta.timestamp[a], ta.values[a], ta.colnames, ta.meta)
end

# single column by name
function getindex(ta::TimeArray, s::AbstractString)
    n = findfirst(ta.colnames, s)
    TimeArray(ta.timestamp, ta.values[:, n], String[s], ta.meta)
end

# array of columns by name
function getindex(ta::TimeArray, args::AbstractString...)
    ns = [findfirst(ta.colnames, a) for a in args]
    TimeArray(ta.timestamp, ta.values[:, ns], String[a for a in args], ta.meta)
end

# single date
function getindex(ta::TimeArray{T, N, D}, d::D) where {T, N, D}
    idxs = searchsorted(ta.timestamp, d)
    length(idxs) == 1 ? ta[idxs[1]] : nothing
end

# multiple dates
function getindex(ta::TimeArray{T, N, D}, dates::Vector{D}) where {T, N, D}
    dates = sort(dates)
    idxs, _ = overlaps(ta.timestamp, dates)
    ta[idxs]
end

# StepRange{Date,...}
getindex(ta::TimeArray{T, N, D}, r::StepRange{D}) where {T, N, D} =
    ta[collect(r)]

getindex(ta::TimeArray, k::TimeArray{Bool, 1}) = ta[findwhen(k)]

# day of week
# getindex{T,N}(ta::TimeArray{T,N}, d::DAYOFWEEK) = ta[dayofweek(ta.timestamp) .== d]

# Define end keyword
endof(ta::TimeArray) = length(ta.timestamp)

# helper methods for inner constructor
function find_dupes_index(cnames)
    idx = Int[]
    for c in 1:length(cnames)
        if contains(string(cnames[1:c-1]), cnames[c])
            push!(idx, c)
        end
    end
    idx
end

function replace_dupes(cnames)
    n=1
    while length(unique(cnames)) != length(cnames)
        ds = find_dupes_index(cnames)
        for d in ds
            if n == 1
                cnames[d] = string(cnames[d], "_$n")
            else
                cnames[d] = string(cnames[d][1:length(cnames[d])-length(string(n))-1], "_$n")
            end
        end
    n +=1
    end
    cnames
end