text:
module AcousticFeatures

using AlphaStableDistributions
using AxisArrays
using DSP
using FFTW
using ImageFiltering: BorderArray, Fill, Pad
using LinearAlgebra
using FindPeaks1D
using Statistics
using StatsBase
using ProgressMeter

export

    # AcousticFeatures
    Energy,
    Myriad,
    # VMyriad,
    FrequencyContours,
    SoundPressureLevel,
    ImpulseStats,
    SymmetricAlphaStableStats,
    Entropy,
    ZeroCrossingRate,
    SpectralCentroid,
    SpectralFlatness,
    PermutationEntropy,
    PSD,
    AcousticComplexityIndex,
    Score,

    # subsequences
    Subsequence,

    # utils
    spectrumflatten,
    myriadconstant,
    vmyriadconstant,
    pressure,
    envelope

include("subsequences.jl")
include("utils.jl")

abstract type AbstractAcousticFeature end

################################################################################
#
#   AcousticFeature types
#
################################################################################
struct Energy <: AbstractAcousticFeature end
name(::Energy) = ["Energy"]

struct Myriad{T<:Union{Nothing,Real}} <: AbstractAcousticFeature
    sqKscale::T
end
Myriad() = Myriad{Nothing}(nothing)
name(::Myriad) = ["Myriad"]

# struct VMyriad{T<:Union{Nothing,Real},M<:Union{Nothing,AbstractMatrix}} <: AbstractAcousticFeature
#     K²::T
#     Rₘ::M
# end
# VMyriad() = VMyriad(nothing, nothing)

struct FrequencyContours{FT<:Real,T<:Real} <: AbstractAcousticFeature
    fs::FT
    n::Int
    nv::Int # overlap
    tnorm::Union{Nothing,T} #time constant for normalization (sec)
    fd::T #frequency difference from one step to the next (Hz)
    minhprc::T
    minfdist::T
    mintlen::T
end
name(::FrequencyContours) = ["Frequency Contours"]

struct SoundPressureLevel{T<:Real} <: AbstractAcousticFeature
    ref::T
end
SoundPressureLevel() = SoundPressureLevel(1.0)
name(::SoundPressureLevel) = ["SPL"]

struct ImpulseStats{FT<:Real,T<:Real} <: AbstractAcousticFeature
    fs::FT
    k::Int
    tdist::T
    computeenvelope::Bool
end
ImpulseStats(fs) = ImpulseStats(fs, 10, 1e-3, true)
ImpulseStats(fs, k, tdist) = ImpulseStats(fs, k, tdist, true)
name(::ImpulseStats) = ["Nᵢ", "μᵢᵢ", "varᵢᵢ"]

struct SymmetricAlphaStableStats <: AbstractAcousticFeature end
name(::SymmetricAlphaStableStats) = ["α", "scale"]

struct Entropy{FT<:Real} <: AbstractAcousticFeature
    fs::FT
    n::Int
    noverlap::Int
end
name(::Entropy) = ["Temporal Entropy","Spectral Entropy","Entropy Index"]

struct ZeroCrossingRate <: AbstractAcousticFeature end
name(::ZeroCrossingRate) = ["ZCR"]

struct SpectralCentroid{FT<:Real} <: AbstractAcousticFeature
    fs::FT
end
name(::SpectralCentroid) = ["Spectral Centroid"]

struct SpectralFlatness <: AbstractAcousticFeature end
name(::SpectralFlatness) = ["Spectral Flatness"]

struct PermutationEntropy <: AbstractAcousticFeature
    m::Int
    τ::Int
    normalization::Bool
end
PermutationEntropy(m) = PermutationEntropy(m, 1, true)
name(::PermutationEntropy) = ["Permutation Entropy"]

struct PSD{FT<:Real} <: AbstractAcousticFeature
    fs::FT
    n::Int
    noverlap::Int
end
function name(f::PSD)
    "PSD-" .* string.(round.(FFTW.rfftfreq(f.n, f.fs); digits=1)) .* "Hz"
end

struct AcousticComplexityIndex{FT<:Real} <: AbstractAcousticFeature
    fs::FT
    n::Int
    noverlap::Int
    jbin::Int
end
name(::AcousticComplexityIndex) = ["Acoustic Complexity Index"]

################################################################################
#
#   Implementations
#
################################################################################
"""
    score(::Energy, x::AbstractVector{T}) where {T<:Real}

Score of `x` based on mean energy.

# Examples:
```julia-repl
julia> x = Score(Energy(), randn(9600))
2-dimensional AxisArray{Float64,2,...} with axes:
    :row, 1:1
    :col, ["Energy"]
And data, a 1×1 Array{Float64,2}:
 0.9960607967861373

julia> x = Score(Energy(), randn(9600); winlen=960, noverlap=480)
2-dimensional AxisArray{Float64,2,...} with axes:
    :row, 1:480:9121
    :col, ["Energy"]
And data, a 20×1 Array{Float64,2}:
 0.5280804987356663
 ⋮
 0.9988797206321275
```
"""
score(::Energy, x::AbstractVector{T}) where T<:Real = [mean(abs2, x)]

"""
    Score(f::Myriad{S}, x::AbstractVector{T})

Score of `x` based on myriad.

# Reference:
Mahmood et. al., "Optimal and Near-Optimal Detection in Bursty Impulsive Noise,"
IEEE Journal of Oceanic Engineering, vol. 42, no. 3, pp. 639--653, 2016.

# Examples:
```julia-repl
julia> x = Score(Myriad(), randn(9600))
2-dimensional AxisArray{Float64,2,...} with axes:
    :row, 1:1
    :col, ["Myriad"]
And data, a 1×1 Array{Float64,2}:
 27691.956992339285

julia> x = Score(Myriad(), randn(9600); winlen=960, noverlap=480)
2-dimensional AxisArray{Float64,2,...} with axes:
    :row, 1:480:9121
    :col, ["Myriad"]
And data, a 20×1 Array{Float64,2}:
 -5487.396124602646
  1977.7969182956683
  3216.396756712948
     ⋮
  2651.158251224668
  3246.7097026864853
```
"""
function score(f::Myriad{S}, x::AbstractVector{T}) where {T<:Real,S<:Real}
    sqKscale = f.sqKscale
    [sum(x -> log(sqKscale + abs2(x)), x)]
end
score(::Myriad{Nothing}, x) = score(Myriad(myriadconstant(x)), x)

# """
#     Score of `x` based on vector myriad.

#     Reference:
#     Mahmood et. al., "Optimal and Near-Optimal Detection in Bursty Impulsive Noise,"
#     IEEE Journal of Oceanic Engineering, vol. 42, no. 3, pp. 639--653, 2016.
# """
# function score(f::VMyriad{S,M}, x::AbstractVector{T}) where {T<:Real,S<:Real,M<:AbstractMatrix}
#     K², Rₘ = f.K², f.Rₘ
#     m = size(Rₘ, 1)-1
#     N = length(x)
#     N < m && throw(ArgumentError("`m` has to be larger than length of `x`"))
#     mplusonedivtwo = (m+1)/2
#     mplustwodivtwo = (m+2)/2
#     onetom = 1:m
#     invRₘ = inv(cholesky(Rₘ).L)
#     invRₘonetom = inv(cholesky(Rₘ[onetom, onetom]).L)
#     s = mplusonedivtwo*log(K²+norm²(invRₘonetom*x[onetom]))
#     for n in m+1:N
#         s += @views mplustwodivtwo*log(K²+norm²(invRₘ*x[n-m:n])) - mplusonedivtwo*log(K²+norm²(invRₘonetom*x[n-m:n-1]))
#     end
#     s
# end
# score(::VMyriad{Nothing,Nothing}, x) = score(VMyriad(vmyriadconstant(x)...), x)

"""
    score(f::FrequencyContours, x::AbstractVector{T})

Score of `x` based on frequency contours count.

# Reference:
D. Mellinger, R. Morrissey, L. Thomas, J. Yosco, "A method for detecting whistles, moans, and other frequency
contour sounds", 2011 J. Acoust. Soc. Am. 129 4055

# Examples:
```julia-repl
julia> x = Score(FrequencyContours(9600, 512, 256, 1.0, 1000.0, 99.0, 1000.0, 0.05), randn(9600))
2-dimensional AxisArray{Float64,2,...} with axes:
    :row, 1:1
    :col, ["Frequency Contours"]
And data, a 1×1 Array{Float64,2}:
 0.0038910505836575876

julia> x = Score(FrequencyContours(9600, 512, 256, 1.0, 1000.0, 99.0, 1000.0, 0.05), randn(9600); winlen=960, noverlap=480)
2-dimensional AxisArray{Float64,2,...} with axes:
    :row, 1:480:9121
    :col, ["Frequency Contours"]
And data, a 20×1 Array{Float64,2}:
 0.0
 0.0
 0.0
 0.0
 ⋮
 0.0
 0.0
 0.0
```
"""
function score(f::FrequencyContours, x::AbstractVector{T}) where T<:Real
    spec = spectrogram(x, f.n, f.nv; fs=f.fs, window=DSP.hamming)
    p  = spec.power; frequency=spec.freq; t=spec.time
    δt = t[2]-t[1]
    δf = frequency[2]-frequency[1]
    f.tnorm === nothing ? Nnorm = size(p, 2) : Nnorm = f.tnorm÷(δt) |> Int
    p    = spectrumflatten(p, Nnorm) #noise-flattened spectrogram
    # crds, _ = peakprom(Maxima(), p[:, 1], trunc(Int, f.minfdist÷δf); minprom=eps(T)+percentile(p[:, 1], f.minhprc))
    crds, _ = findpeaks1d(p[:, 1]; height=eps(T)+percentile(p[:, 1], f.minhprc), distance=trunc(Int, f.minfdist/δf))
    ctrs = [[(crd, 1)] for crd in crds]
    for (i, col) in enumerate(eachcol(p[:, 2:end]))
        col = collect(col)
        # crds,_ = peakprom(Maxima(), col, trunc(Int, f.minfdist/δf); minprom=eps(T)+percentile(col, f.minhprc))
        crds, _ = findpeaks1d(col; height=eps(T)+percentile(col, f.minhprc), distance=trunc(Int, f.minfdist/δf))
        for crd in crds
            if length(ctrs) == 0
                ctrs = [[(crd, 1)] for crd in crds]
            else
                idxselect = Int64[]
                costselect = Float64[]
                for (j, ctr) in enumerate(ctrs)
                    if (ctr[end][2] == i-1) && abs(frequency[ctr[end][1]]-frequency[crd]) <= f.fd
                        push!(idxselect, j)
                        push!(costselect, abs(frequency[ctr[end][1]]-frequency[crd]))
                    end
                end
                if isempty(idxselect)
                    push!(ctrs, [(crd, i)])
                else
                    idxopt = idxselect[argmin(costselect)]
                    push!(ctrs[idxopt], (crd, i))
                end
            end
        end
    end
    idxdelete = Int64[]
    for (i, ctr) in enumerate(ctrs)
        (length(ctr)-1)*(δt) < f.mintlen && push!(idxdelete, i)
    end
    deleteat!(ctrs, idxdelete)
    count = isempty(ctrs) ? 0 : sum(length, ctrs)
    [count/length(p)]
end

"""
    score(f::SoundPressureLevel, x::AbstractVector{T})

Score of `x` based on Sound Pressure Level (SPL). `x` is in micropascal.
In water, the common reference `ref` is 1 micropascal. In air, the
common reference `ref` is 20 micropascal.

# Examples:
```julia-repl
julia> x = Score(SoundPressureLevel(), randn(9600))
2-dimensional AxisArray{Float64,2,...} with axes:
    :row, 1:1
    :col, ["SPL"]
And data, a 1×1 Array{Float64,2}:
 -0.08307636105819256

julia> x = Score(SoundPressureLevel(), randn(9600); winlen=960, noverlap=480)
2-dimensional AxisArray{Float64,2,...} with axes:
    :row, 1:480:9121
    :col, ["SPL"]
And data, a 20×1 Array{Float64,2}:
 -2.369874304880999
  0.2795371978218069
  ⋮
 -0.04985476533237352
  0.11412307503113574
```
"""
function score(f::SoundPressureLevel, x::AbstractVector{T}) where T<:Real
    rmsx = sqrt(mean(abs2, x))
    [20*log10(rmsx/f.ref)]
end

"""
    score(f::ImpulseStats, x::AbstractVector{T})

Score of `x` based on number of impulses, mean and variance of inter-impulse intervals.
The minimum height of impulses is defined by `a+k*b` where `a` is median of the envelope
of `x` and `b` is median absolute deviation (MAD) of the envelope of `x`.

# Reference:
Matthew W Legg et al., "Analysis of impulsive biological noise due to snapping shrimp as a
point process in time", 2007.

# Examples:
```julia-repl
julia> x = Score(ImpulseStats(9600, 10, 0.01), randn(9600))
2-dimensional AxisArray{Float64,2,...} with axes:
    :row, 1:1
    :col, ["Nᵢ", "μᵢᵢ", "varᵢᵢ"]
And data, a 1×3 Array{Float64,2}:
 0.0  0.0  0.0

julia> x = Score(ImpulseStats(9600, 10, 0.01), randn(9600); winlen=960, noverlap=480)
2-dimensional AxisArray{Float64,2,...} with axes:
    :row, 1:480:9121
    :col, ["Nᵢ", "μᵢᵢ", "varᵢᵢ"]
And data, a 20×3 Array{Float64,2}:
 4.0  0.0140972  0.0621181
 0.0  0.0        0.0
 0.0  0.0        0.0
 0.0  0.0        0.0
 ⋮
 0.0  0.0        0.0
 0.0  0.0        0.0
 0.0  0.0        0.0
```
"""
function score(f::ImpulseStats, x::AbstractVector{T}) where T<:Real
    if f.computeenvelope
        x = envelope(x)
    end
    center = Statistics.median(x)
    height = center+f.k*mad(x, center=center, normalize=true)
    distance = trunc(Int, f.tdist*f.fs)
    # crds, _ = peakprom(Maxima(), x, distance; minprom=height)
    crds,_ = findpeaks1d(x; height=height, distance=distance)
    timeintervals = diff(crds)
    [convert(Float64, length(crds)), mean(timeintervals)/f.fs, var(timeintervals)/f.fs]
end

"""
    score(::SymmetricAlphaStableStats, x::AbstractVector{T})

Score of `x` based on the parameters of Symmetric Alpha Stable Distributions.
The parameter α measures the impulsiveness while the parameter scale measures
the width of the distributions.

# Reference:
https://github.com/org-arl/AlphaStableDistributions.jl

# Examples:
```julia-repl
julia> x = Score(SymmetricAlphaStableStats(), randn(9600))
2-dimensional AxisArray{Float64,2,...} with axes:
    :row, 1:1
    :col, ["α", "scale"]
And data, a 1×2 Array{Float64,2}:
 2.0  0.714388

julia> x = Score(SymmetricAlphaStableStats(), randn(9600); winlen=960, noverlap=480)
2-dimensional AxisArray{Float64,2,...} with axes:
    :row, 1:480:9121
    :col, ["α", "scale"]
And data, a 20×2 Array{Float64,2}:
 0.5      0.0
 1.90067  0.663918
 1.83559  0.614218
 ⋮
 1.80072  0.676852
 1.94506  0.677581
```
"""
function score(::SymmetricAlphaStableStats, x::AbstractVector{T}) where T<:Real
    d = fit(SymmetricAlphaStable, x)
    [d.α, d.scale]
end

"""
    score(f::Entropy, x::AbstractVector{T})

Score of `x` based on temporal entropy, spectral entropy and entropy index.

# Reference:
J. Sueur, A. Farina, A. Gasc, N. Pieretti, S. Pavoine, Acoustic Indices for Biodiversity
Assessment and Landscape Investigation, 2014.

# Examples:
```julia-repl
julia> x = Score(Entropy(9600, 96, 48), randn(9600))
2-dimensional AxisArray{Float64,2,...} with axes:
    :row, 1:1
    :col, ["Temporal Entropy", "Spectral Entropy", "Entropy Index"]
And data, a 1×3 Array{Float64,2}:
 0.984457  0.997608  0.982103

julia> x = Score(Entropy(9600, 96, 48), randn(9600); winlen=960, noverlap=480)
2-dimensional AxisArray{Float64,2,...} with axes:
    :row, 1:480:9121
    :col, ["Temporal Entropy", "Spectral Entropy", "Entropy Index"]
And data, a 20×3 Array{Float64,2}:
 0.903053  0.982299  0.887068
 0.980151  0.986018  0.966446
 0.981492  0.984845  0.966618
 0.980283  0.986635  0.967182
 0.978714  0.987383  0.966366
 ⋮
 0.97895   0.986322  0.96556
 0.980274  0.983338  0.963941
 0.980822  0.99296   0.973917
 0.979092  0.989817  0.969122
```
"""
function score(f::Entropy, x::AbstractVector{T}) where T<:Real
    sp = spectrogram(x, f.n, f.noverlap; fs=f.fs).power
    ne = normalize_envelope(x)
    n = length(ne)
    Ht = -sum(ne .* log2.(ne)) ./ log2(n)
    ns = normalize_spectrum(sp)
    N = length(ns)
    Hf = -sum(ns .* log2.(ns)) ./ log2.(N)
    H = Ht*Hf
    [Ht, Hf, H]
end

"""
    score(::ZeroCrossingRate, x::AbstractVector{T})

Score of `x` based on zero crossing rate.

# Refernce:
https://en.wikipedia.org/wiki/Zero-crossing_rate

# Examples:
```julia-repl
julia> x = Score(ZeroCrossingRate(), randn(9600))
2-dimensional AxisArray{Float64,2,...} with axes:
    :row, 1:1
    :col, ["ZCR"]
And data, a 1×1 Array{Float64,2}:
 0.5027083333333333

julia> x = Score(ZeroCrossingRate(), randn(9600); winlen=960, noverlap=480)
2-dimensional AxisArray{Float64,2,...} with axes:
    :row, 1:480:9121
    :col, ["ZCR"]
And data, a 20×1 Array{Float64,2}:
 0.2375
 0.46979166666666666
 0.4708333333333333
 0.49583333333333335
 0.5072916666666667
 ⋮
 0.5145833333333333
 0.49375
 0.5052083333333334
 0.5125
 0.4947916666666667
```
"""
function score(::ZeroCrossingRate, x::AbstractVector{T}) where T<:Real
    [count(!iszero, diff(x .> 0))/length(x)]
end

"""
    score(f::SpectralCentroid, x::AbstractVector{T})

Score of `x` based on spectral centroid.

# Reference:
https://en.wikipedia.org/wiki/Spectral_centroid

# Examples:
```julia-repl
julia> x = Score(SpectralCentroid(9600), randn(9600))
2-dimensional AxisArray{Float64,2,...} with axes:
    :row, 1:1
    :col, ["Spectral Centroid"]
And data, a 1×1 Array{Float64,2}:
 2387.4592177121676

julia> x = Score(SpectralCentroid(9600), randn(9600); winlen=960, noverlap=480)
2-dimensional AxisArray{Float64,2,...} with axes:
    :row, 1:480:9121
    :col, ["Spectral Centroid"]
And data, a 20×1 Array{Float64,2}:
 2398.6889311658415
 2362.570125358973
 2342.6919660952803
 2375.5903954079977
    ⋮
 2415.431353476017
 2453.7105902333437
 2449.222535628719
 2380.319011105224
```
"""
function score(f::SpectralCentroid, x::AbstractVector{T}) where T<:Real
    magnitudes = abs.(rfft(x))
    freqs = FFTW.rfftfreq(length(x), f.fs)
    [sum(magnitudes .* freqs) / sum(magnitudes)]
end

"""
    score(::SpectralFlatness, x::AbstractVector{T})

Score of `x` based on spectral flatness.

# Reference:
https://en.wikipedia.org/wiki/Spectral_flatness

# Examples:
```julia-repl
julia> x = Score(SpectralFlatness(), randn(9600))
2-dimensional AxisArray{Float64,2,...} with axes:
    :row, 1:1
    :col, ["Spectral Flatness"]
And data, a 1×1 Array{Float64,2}:
 0.5598932661540399

julia> x = Score(SpectralFlatness(), randn(9600); winlen=960, noverlap=480)
2-dimensional AxisArray{Float64,2,...} with axes:
    :row, 1:480:9121
    :col, ["Spectral Flatness"]
And data, a 20×1 Array{Float64,2}:
 0.5661636483227057
 0.543740942647357
 0.5854629162797854
 0.532148471407988
 ⋮
 0.5451002001200566
 0.5496170417608265
 0.537335859768483
 0.535056705285523
```
"""
function score(::SpectralFlatness, x::AbstractVector{T}) where T<:Real
    magnitudes² = (abs.(rfft(x))).^2
    [geomean(magnitudes²) / mean(magnitudes²)]
end

"""
    score(f::PermutationEntropy, x::AbstractVector{T})

Score of `x` based on permutation entropy.

# Reference:
C. Bandt, B. Pompe, "Permutation entropy: a natural complexity measure for time series",
Phys. Rev. Lett., 88 (17), 2002

# Examples:
```julia-repl
julia> x = Score(PermutationEntropy(7), randn(9600))
2-dimensional AxisArray{Float64,2,...} with axes:
    :row, 1:1
    :col, ["Permutation Entropy"]
And data, a 1×1 Array{Float64,2}:
 0.9637270776723836

julia> x = Score(PermutationEntropy(7), randn(9600), winlen=960, noverlap=480)
2-dimensional AxisArray{Float64,2,...} with axes:
    :row, 1:480:9121
    :col, ["Permutation Entropy"]
And data, a 20×1 Array{Float64,2}:
 0.4432867336969194
 0.7896679491573086
 0.7914368148634888
 0.790455877419772
 ⋮
 0.7894106035823039
 0.7886452315698513
 0.7897322855510599
 0.7884747786386084
```
"""
function score(f::PermutationEntropy, x::AbstractVector{T}) where T<:Real
    p = ordinalpatterns(x, f.m, f.τ)
    pe = -sum(p .* log2.(p))
    if f.normalization
        [pe / convert(eltype(pe), log2(factorial(big(f.m))))]
    else
        [pe]
    end
end

"""
    score(f::PSD, x::AbstractVector{T})

Score of `x` based on power spectral density in dB scale.

# Examples:
```julia-repl
julia> x = Score(PSD(96000, 64, 32), randn(9600))
2-dimensional AxisArray{Float64,2,...} with axes:
    :row, 1:1
    :col, ["PSD-0Hz", "PSD-1500Hz", "PSD-3000Hz", "PSD-4500Hz", "PSD-6000Hz", "PSD-7500Hz", "PSD-9000Hz", "PSD-10500Hz", "PSD-12000Hz", "PSD-13500Hz"  …  "PSD-34500Hz", "PSD-36000Hz", "PSD-37500Hz", "PSD-39000Hz", "PSD-40500Hz", "PSD-42000Hz", "PSD-43500Hz", "PSD-45000Hz", "PSD-46500Hz", "PSD-48000Hz"]
And data, a 1×33 Array{Float64,2}:
 -49.611  -47.1275  -46.7286  -46.4742  -46.6452  …  -47.0801  -47.2065  -46.5577  -46.4154  -50.4786

julia> x = Score(PSD(96000, 64, 32), randn(9600), winlen=960, noverlap=480)
2-dimensional AxisArray{Float64,2,...} with axes:
    :row, 1:480:9121
    :col, ["PSD-0Hz", "PSD-1500Hz", "PSD-3000Hz", "PSD-4500Hz", "PSD-6000Hz", "PSD-7500Hz", "PSD-9000Hz", "PSD-10500Hz", "PSD-12000Hz", "PSD-13500Hz"  …  "PSD-34500Hz", "PSD-36000Hz", "PSD-37500Hz", "PSD-39000Hz", "PSD-40500Hz", "PSD-42000Hz", "PSD-43500Hz", "PSD-45000Hz", "PSD-46500Hz", "PSD-48000Hz"]
And data, a 20×33 Array{Float64,2}:
 -54.3507  -49.2203  -47.9715  -50.4869  -51.4884  …  -48.8581  -49.8427  -50.804   -47.5865  -51.502
 -48.4897  -45.3653  -46.2605  -47.4234  -47.4109     -46.3487  -46.3496  -48.243   -45.3489  -49.5331
 -48.0083  -45.1101  -45.8474  -47.284   -45.2843     -46.2456  -46.5622  -47.2382  -46.6219  -47.9124
 -49.023   -45.4548  -44.5266  -45.9787  -44.7397     -47.6285  -48.4443  -47.2613  -47.7996  -48.1423
   ⋮                                               ⋱                        ⋮
 -49.9071  -46.6817  -47.1582  -45.9655  -48.3396     -46.986   -46.8983  -45.6008  -47.0211  -48.4817
 -49.0467  -47.1668  -46.9087  -47.0215  -47.8279     -46.8043  -47.2044  -45.6053  -47.0023  -48.222
 -49.7118  -47.3381  -47.219   -45.3647  -45.6587     -47.3541  -47.4126  -46.1465  -46.491   -48.1833
```
"""
function score(f::PSD, x::AbstractVector{T}) where T<:Real
    p = welch_pgram(x, f.n, f.noverlap; fs=f.fs)
    pow2db.(power(p))
end

"""
    score(f::AcousticComplexityIndex, x::AbstractVector{T})

Score of `x` based on Acoustic Complexity Index. `jbin` is the temporal size of a sub-spectrogram.

# Reference:
N. Pieretti, A. Farina, D. Morri, "A new methodology to infer the singing activity of an avian community: The Acoustic Complexity Index (ACI)", 2011.

# Examples:
```julia-repl
julia> x = Score(AcousticComplexityIndex(96000, 1024, 0, 30), randn(960000))
2-dimensional AxisArray{Float64,2,...} with axes:
    :row, 1:1
    :col, ["Acoustic Complexity Index"]
And data, a 1×1 Matrix{Float64}:
 15394.052148047322

julia> x = Score(AcousticComplexityIndex(96000, 1024, 0, 30), randn(960000), winlen=96000, noverlap=48000)
2-dimensional AxisArray{Float64,2,...} with axes:
    :row, 1:48000:912001
    :col, ["Acoustic Complexity Index"]
And data, a 20×1 Matrix{Float64}:
 998.0694761443425
 1493.1632775077805
    ⋮
 1496.4075980628913
 1486.5057244355821
```
"""
function score(f::AcousticComplexityIndex, x::AbstractVector{T}) where T<:Real
    sp = spectrogram(x, f.n, f.noverlap; fs=f.fs).power
    m, n = size(sp)
    starts = range(1, n-f.jbin+1; step=f.jbin)
    nstarts = length(starts)
    acitmp = Matrix{T}(undef, m, nstarts)
    for i ∈ 1:nstarts
        subsp = sp[:,starts[i]:starts[i]+f.jbin-1]
        D = sum(abs.(diff(subsp; dims=2)); dims=2)
        acitmp[:,i] = D ./ (sum(subsp; dims=2) .+ eps(T))
    end
    [sum(acitmp)]
end

"""
    Score(f,
          x;
          winlen=length(x),
          noverlap=0,
          padtype=:fillzeros,
          subseqtype=Float64,
          preprocess=identity,
          map=map,
          showprogress=true)

Compute acoustic feature `f` scores of a time series signal `x` using sliding windows.

By default, window length `winlen` is the length of `x`, i.e., the whole signal is used to compute
a score, and overlapping samples `noverlap` is 0. The `padtype` specifies the form of padding, and
for more information, refer to `ImageFiltering.jl`. The signal is subject to preliminary processing
`preprocess`. Acoustic feature scores of subseqences can be computed through mapping
`map`. `showprogress` is used to monitor the computations.
"""
function Score(f::AbstractAcousticFeature,
               x::AbstractVector{T};
               winlen::Int = length(x),
               noverlap::Int = 0,
               padtype::Symbol = :fillzeros,
               subseqtype::DataType = Float64,
               preprocess::Function = identity,
               map::Function = map,
               showprogress::Bool = false) where {T<:Real}
    (noverlap < 0) && (return throw(ArgumentError("`noverlap` must be larger or equal to zero.")))
    xlen = length(x)
    if winlen < xlen
        subseq = Subsequence(x, winlen, noverlap; padtype=padtype)
        if showprogress
            s = (@showprogress map(x -> score(f, preprocess(convert.(subseqtype, x))), subseq))
        else
            s = map(x -> score(f, preprocess(convert.(subseqtype, x))), subseq)
        end
        return AxisArray(vcat(reshape.(s, 1, :)...)::Matrix{subseqtype};
                         row=collect(1:step(subseq):xlen),
                         col=name(f))
    elseif winlen == xlen
        stmp = score(f, preprocess(convert.(subseqtype, x)))
        return AxisArray(reshape(stmp, (1, length(stmp))); row=ones(Int,1), col=name(f))
    else
        return throw(ArgumentError("`winlen` must be smaller or equal to the length of `x`."))
    end
end

(f::AbstractAcousticFeature)(x; kwargs...) = Score(f, x; kwargs...)

end


text:
#!/usr/bin/env julia

using Luxor

using Test

using Random
Random.seed!(42)

function getpath_1(fname)
    Drawing(400, 250, fname)

    background("white")
    origin()
    setline(0.75)
    sethue("black")

    fontsize(220)
    t = "N" # no curves on this one
    translate(-textextents(t)[3]/2, textextents(t)[4]/2)
    textpath(t)
    pathdata = getpath()
    outline = Point[]
    for i in pathdata[1:end-1]
        if length(i.points) == 2
            x = i.points[1]
            y = i.points[2]
            push!(outline, Point(x, y))
        end
    end
    poly(outline, :stroke, close=true)
    for i in 5:5:35
        poly(offsetpoly(outline, i), :stroke, close=true)
    end
    @test finish() == true
end

fname= "figures-get-path.pdf"
getpath_1(fname)
println("...finished test: output in $(fname)")

ext: jl


text:
peaks(x::AbstractArray,y::AbstractArray) = 3*(1-x).^2.*exp(-(x.^2) - (y+1).^2) - 10*(x/5 - x.^3 - y.^5).*exp(-x.^2-y.^2) - 1/3*exp(-(x+1).^2 - y.^2)
peaks(x::Number,y::Number) = 3*(1-x)^2*exp(-(x^2) - (y+1)^2) - 10*(x/5 - x^3 - y^5)*exp(-x^2-y^2) - 1/3*exp(-(x+1)^2 - y^2)
peaks(v::AbstractVector) = peaks(v[1],v[2])

Power(a,b) = a^b
E = e
grad_peaks(x::Number, y::Number) = SVector((-2*Power(E,-2*x - Power(x,2) - Power(1 + y,2))* (-(Power(E,2*y)*(1 + x)) + 9*Power(E,2*x)*(1 - 2*Power(x,2) + Power(x,3)) + Power(E,1 + 2*x + 2*y)*(3 - 51*Power(x,2) + 30*Power(x,4) + 30*x*Power(y,5))))/3., (-2*Power(E,-2*x - Power(x,2) - Power(1 + y,2))* (-(Power(E,2*y)*y) + 9*Power(E,2*x)*Power(-1 + x,2)*(1 + y) + 3*Power(E,1 + 2*x + 2*y)*y*(-2*x + 10*Power(x,3) + 5*Power(y,3)*(-5 + 2*Power(y,2)))))/3.)
grad_peaks(v::AbstractVector) = grad_peaks(v[1], v[2])

grad_scaled_peaks(x, y) = 0.1*grad_peaks(x/1000, y/1000)
grad_scaled_peaks(v::AbstractVector) = 0.1*grad_peaks(v/1000)
scaled_peaks(x,y) = 100*peaks(x/1000, y/1000)
scaled_peaks(v::AbstractVector) = 100*peaks(v/1000)

function plot_peaks(n=100)
    x = linspace(-4, 4, n)
    y = linspace(-4, 4, n)
    surface(x,y,peaks)
end

function plot_scaled_peaks(n=100)
    x = linspace(-4000, 4000, n)
    y = linspace(-4000, 4000, n)
    surface(x,y,scaled_peaks, ratio=:equal)
end

peaks2(x,y) = 3*(1-x)^2*exp(-(x^2) - (y+1)^2) - 7*(x/5 - x^2 - y^4)*exp(-x^2/2-(y)^2) - 1/3*exp(-(x+1)^2 - y^2)
peaks2(v) = peaks2(v[1], v[2])
scaled_peaks2(x,y) = 70*peaks2(x/1000, y/1000)
scaled_peaks2(v) = scaled_peaks2(v[1], v[2])

grad_peaks2(x,y) = SVector((exp(-2*x - Power(x,2) - Power(1 + y,2))*
      (10*exp(2*y)*(1 + x) - 90*exp(2*x)*(1 - 2*Power(x,2) + Power(x,3)) -
        21*exp(1 + 2*x + Power(x,2)/2. + 2*y)*
         (1 - Power(x,2) + 5*Power(x,3) + 5*x*(-2 + Power(y,4)))))/15.,
   (2*exp(-Power(1 + x,2) - Power(y,2))*y)/3. +
    28*exp(-Power(x,2)/2. - Power(y,2))*Power(y,3) -
    6*exp(-Power(x,2) - Power(1 + y,2))*Power(-1 + x,2)*(1 + y) -
    (14*exp(-Power(x,2)/2. - Power(y,2))*y*(-x + 5*Power(x,2) + 5*Power(y,4)))/5.)
grad_peaks2(v) = grad_peaks2(v[1], v[2])

grad_scaled_peaks2(v) = 0.07*grad_peaks2(v[1]/1000, v[2]/1000)

#=
grad_peaks2(x,y) = SVector((Power(E,-2*x - Power(x,2) - Power(1 + y,2))*
      (10*Power(E,2*y)*(1 + x) - 90*Power(E,2*x)*(1 - 2*Power(x,2) + Power(x,3)) -
        21*Power(E,1 + 2*x + Power(x,2)/2. + 2*y)*
         (1 - Power(x,2) + 5*Power(x,3) + 5*x*(-2 + Power(y,4)))))/15.,
   (2*Power(E,-Power(1 + x,2) - Power(y,2))*y)/3. +
    28*Power(E,-Power(x,2)/2. - Power(y,2))*Power(y,3) -
    6*Power(E,-Power(x,2) - Power(1 + y,2))*Power(-1 + x,2)*(1 + y) -
    (14*Power(E,-Power(x,2)/2. - Power(y,2))*y*(-x + 5*Power(x,2) + 5*Power(y,4)))/5.)
=#


text:
"""
   set_position!(legend, s)

Internal function to set the position of the legend using a describer like "tl" (top-left).
"""
function set_position!(o::Union{Legend,Box2D}, v::String)
   @assert get_backend() == GLE "position/only GLE backend supported"
   o.position = get(GLE_POSITION, v) do
      throw(OptionValueError("position", v))
   end
   return nothing
end

"""
   set_nobox!(o, b)

Internal function to toggle nobox on or off for an object.
"""
set_nobox!(o::Union{Legend,Box2D,Colorbar}, b::Bool) = (o.nobox = b)

"""
   set_margins!(legend, m)

Internal function to set the margins of the legend (internal distance from legend box to elements).
"""
set_margins!(o::Legend, m::T2F) = (o.margins = m)

"""
   set_offset!(obj, m)

Internal function to set the offset of the object (external distance from object to axis).
"""
set_offset!(o::Legend, v::T2F) = (o.offset = v)


text:
@doc raw"""
LevStehnoexp{T} <: Kernel
The kernel for of the levensteindistance
`` k(x, x^\prime) = \exp(-levensthein(x_i,x_j))``

"""

struct LevStehnoexp{T} <: Kernel
    s::T
end

"""
levensthein(xᵢ,xⱼ)
dynamic calculations of the model the Levenstein distance between two vectors of modules

"""
function levenshtein(xᵢ, xⱼ)
    n, m = length(xᵢ), length(xⱼ)
    lev = zeros(Int, n+1, m+1)
    lev[:,1] = 0:n
    lev[1,:] = 0:m
    for j in 1:m
        for i in 1:n
            lev[i+1, j+1] = min(
                lev[i, j+1] + 1,
                lev[i+1, j] + 1,
                xᵢ[i]==xⱼ[j] ? lev[i, j] : lev[i, j] + 1
            )
        end
    end
    return last(lev)
end
ew(k::LevStehnoexp, x::AbstractVector{N} where N) = exp.(-k.s .*zeros(length(x)))
ew(k::LevStehnoexp, x::AbstractVector{N} where N, x′::AbstractVector{N} where N) = [exp(-k.s*levenshtein(xᵢ,xⱼ)) for (xᵢ,xⱼ) in zip(x, x′)]

pw(k::LevStehnoexp, x::AbstractVector{N} where N) = reshape([exp(-k.s*levenshtein(xᵢ,xⱼ)) for xᵢ in x for xⱼ in x] ,(length(x),length(x)))
pw(k::LevStehnoexp, x::AbstractVector{N} where N, x′::AbstractVector{N} where N) = reshape([exp(-k.s*levenshtein(xᵢ,xⱼ)) for xⱼ in x′ for xᵢ in x] ,(length(x),length(x′)))

####
#Cosine kernel
####

####
#make vectors first
####

function myword2vec(word,letter2index)
    word2vec = letter2index  |> t -> length(t) |> l-> zeros(l)
    for letter in word
        word2vec[letter2index[letter]] += 1
    end
    return word2vec
end

function cossim(xi,xj)
    letter2index = Dict( letter => index for (index,letter) in enumerate(Set([xi...,xj...])))
    xi_v = myword2vec(xi,letter2index)
    xj_v = myword2vec(xj,letter2index)
    return dot(xi_v,xj_v)/(norm(xi_v) * norm(xj_v))
end

struct CosStehno <: Kernel end

ew(::CosStehno, x::AbstractVector{N} where N) = ones(length(x))
ew(::CosStehno, x::AbstractVector{N} where N, x′::AbstractVector{N} where N) = [cossim(xᵢ,xⱼ) for (xᵢ,xⱼ) in zip(x, x′)]

pw(k:: CosStehno, x::AbstractVector{N} where N) = reshape([cossim(xᵢ,xⱼ) for xᵢ in x for xⱼ in x] ,(length(x),length(x)))
pw(k:: CosStehno, x::AbstractVector{N} where N, x′::AbstractVector{N} where N) = reshape([cossim(xᵢ,xⱼ) for xⱼ in x′ for xᵢ in x] ,(length(x),length(x′)))


text:
struct LJSAFTParam <: EoSParam
    segment::SingleParam{Float64}
    b::PairParam{Float64}
    T_tilde::PairParam{Float64}
    epsilon_assoc::AssocParam{Float64}
    bondvol::AssocParam{Float64}
end

abstract type LJSAFTModel <: SAFTModel end
@newmodel LJSAFT LJSAFTModel LJSAFTParam

export LJSAFT
function LJSAFT(components; idealmodel=BasicIdeal, userlocations=String[], ideal_userlocations=String[], verbose=false)
    params,sites = getparams(components, ["SAFT/LJSAFT"]; userlocations=userlocations, verbose=verbose)
    segment = params["m"]

    k = params["k"]
    zeta = params["zeta"]

    T_tilde = epsilon_LorentzBerthelot(params["T_tilde"], k)
    params["b"].values .*= 1E-3
    params["b"].values .^= 1/3
    b = sigma_LorentzBerthelot_mod(params["b"],zeta)
    b.values .^= 3
    epsilon_assoc = params["epsilon_assoc"]
    bondvol = params["bondvol"]

    packagedparams = LJSAFTParam(segment, b, T_tilde, epsilon_assoc, bondvol)
    references = ["10.1021/ie9602320"]

    model = LJSAFT(packagedparams, sites, idealmodel; ideal_userlocations=ideal_userlocations, references=references, verbose=verbose)
    return model
end

function a_res(model::LJSAFTModel, V, T, z)
    return @f(a_seg) + @f(a_chain) + @f(a_assoc)
end

function a_seg(model::LJSAFTModel, V, T, z)
    D = LJSAFTconsts.D
    C = LJSAFTconsts.C
    C0 = LJSAFTconsts.C0
    C1 = LJSAFTconsts.C1
    C2 = LJSAFTconsts.C2
    C4 = LJSAFTconsts.C4
    x = z/∑(z)
    m = model.params.segment.values

    T̃ = @f(Tm)
    b̄ = @f(bm)

    Tst = T/T̃
    m̄ = ∑(m .* x)
    ρ = ∑(z)/V
    ρst = m̄*b̄*ρ
    η = ρst*π/6*(∑(D[i+3]*Tst^(i/2) for i ∈ -2:1)+D[end]*log(Tst))^3

    A_HS = Tst*(5/3*log(1-η)+(η*(34-33η+4η^2))/(6*(1-η)^2))
    ΔB2 = ∑(C[j+8]*Tst^(j/2) for j ∈ -7:0)
    A0 = ∑(C0[j-1]*ρst^j for j ∈ 2:5)
    A1 = ∑(C1[j-1]*Tst^(-1/2)*ρst^j for j ∈ 2:6)
    A2 = ∑(C2[j-1]*Tst^(-1)*ρst^j for j ∈ 2:6)
    A4 = ∑(C4[j-1]*Tst^(-2)*ρst^j for j ∈ 2:6)

    γ = 1.92907278
    return m̄*(A_HS+exp(-γ*ρst^2)*ρst*Tst*ΔB2+A0+A1+A2+A4)/Tst
end

function Tm(model::LJSAFTModel, V, T, z)
    x = z/∑(z)
    T̃ = model.params.T_tilde.values
    b = model.params.b.values
    m = model.params.segment.values
    comps = @comps
    return ∑(m[i]*m[j]*x[i]*x[j]*b[i,j]*T̃[i,j] for i ∈ 1:length(z) for j ∈ comps)/∑(m[i]*m[j]*x[i]*x[j]*b[i,j] for i ∈ comps for j ∈ comps)
end

function bm(model::LJSAFTModel, V, T, z)
    x = z/∑(z)
    comps = @comps
    b = model.params.b.values
    m = model.params.segment.values
    return ∑(m[i]*m[j]*x[i]*x[j]*b[i,j] for i ∈ comps for j ∈ comps)/∑(m[i]*m[j]*x[i]*x[j] for i ∈ comps for j ∈ comps)
end

function a_chain(model::LJSAFTModel, V, T, z)
    x = z/∑(z)
    m = model.params.segment.values
    return -∑(x[i]*(m[i]-1)*log(@f(g_LJ,i)) for i ∈ @comps)
end

function g_LJ(model::LJSAFTModel, V, T, z, i)
    ∑z = ∑(z)
    x = z/∑(z)
    m = model.params.segment.values
    b = model.params.b.diagvalues
    T̃ = model.params.T_tilde.diagvalues

    Tst = T/T̃[i]
    ρ = ∑z/V
    ρ̄ = b[i]*m[i]*x[i]*ρ
    a = LJSAFTconsts.a

    return (1+∑(a[i,j]*ρ̄^i*Tst^(1-j) for i ∈ 1:5 for j ∈ 1:5))
end

function a_assoc(model::LJSAFTModel, V, T, z)
    x = z/∑(z)
    n = model.sites.n_sites
    X_ = @f(X)
    return ∑(x[i]*∑(n[i][a]*(log(X_[i][a])+(1-X_[i][a])/2) for a ∈ @sites(i)) for i ∈ @comps)
end

function X(model::LJSAFTModel, V, T, z)
    _1 = one(V+T+first(z))
    ∑z = ∑(z)
    x = z/∑z
    ρ = ∑z/V
    n = model.sites.n_sites
    itermax = 500
    dampingfactor = 0.5
    error = 1.
    tol = model.absolutetolerance
    iter = 1
    X_ = [[_1 for a ∈ @sites(i)] for i ∈ @comps]
    X_old = deepcopy(X_)
    while error > tol
        iter > itermax && error("X has failed to converge after $itermax iterations")
        for i ∈ @comps, a ∈ @sites(i)
            rhs = 1/(1+∑(ρ*x[j]*∑(n[j][b]*X_old[j][b]*@f(Δ,i,j,a,b) for b ∈ @sites(j)) for j ∈ @comps))
            X_[i][a] = (1-dampingfactor)*X_old[i][a] + dampingfactor*rhs
        end
        error = sqrt(∑(∑((X_[i][a] - X_old[i][a])^2 for a ∈ @sites(i)) for i ∈ @comps))
        for i = 1:length(X_)
            X_old[i] .= X_[i]
        end
        iter += 1
    end
    return X_
end

function Δ(model::LJSAFTModel, V, T, z, i, j, a, b)
    ∑z = ∑(z)
    x = z/∑z
    m = model.params.segment.values
    _b = model.params.b.values
    T̃ = model.params.T_tilde.values

    Tst = T/T̃[i,j]
    ρ = ∑z/V
    ρ̄ = x[i]*_b[i,j]*m[i]*ρ
    ϵ_assoc = model.params.epsilon_assoc.values
    κ = model.params.bondvol.values
    b_ = LJSAFTconsts.b
    I = ∑(b_[i+1,j+1]*ρ̄^i*Tst^j for i ∈ 0:4 for j ∈ 0:4)/3.84/1e4
    return 4π*(exp(ϵ_assoc[i,j][a,b]/T)-1)*κ[i,j][a,b]*I*_b[i,j]
end

const LJSAFTconsts = (
    a = [0.49304346593882 2.1528349894745 -15.955682329017 24.035999666294 -8.6437958513990;
           -0.47031983115362 1.1471647487376  37.889828024211 -84.667121491179 39.643914108411;
            5.0325486243620 -25.915399226419 -18.862251310090 107.63707381726 -66.602649735720;
           -7.3633150434385  51.553565337453 -40.519369256098 -38.796692647218 44.605139198378;
            2.9043607296043 -24.478812869291  31.500186765040 -5.3368920371407 -9.5183440180133],
    b = [-0.03915181 0.08450471 0.06889053 -0.01034279  0.5728662e-3;
             -0.5915018  0.9838141 -0.4862279   0.1029708 -0.6919154e-2;
               1.908368  -3.415721   2.124052  -0.4298159  0.02798384;
             -0.7957312  0.7187330 -0.9678804   0.2431675 -0.01644710;
             -0.9399577   2.314054 -0.4877045  0.03932058 -0.1600850e-2],
    D = [0.011117524,-0.076383859,1.080142248, 0.000693129,-0.063920968],
    C = [-0.58544978,0.43102052,0.87361369,-4.13749995,2.90616279,-7.02181962,0.,0.0245987],
    C0 = [2.01546797,-28.17881636,28.28313847,-10.42402873],
    C1 = [-19.58371655,75.62340289,-120.70586598,93.92740328,-27.37737354],
    C2 = [29.34470520,-112.35356937,170.64908980,-123.06669187,34.42288969],
    C4 = [-13.37031968,65.38059570,-115.09233113,88.91973082,-25.62099890]

   )


text:
################################################# FILE DESCRIPTION #########################################################

# This file contains the implementation of the Vertex DataFrame.

################################################# IMPORT/EXPORT ############################################################
export listvprops, hasvprop, getvprop, setvprop!

################################################# MUTATION #################################################################

""" Add a row to the vertex dataframe """
function addvertex!(x::AbstractDataFrame)
   push!(x, @data fill(NA, ncol(x)))
   return
end

""" Remove a row from the vertex dataframe """
function rmvertex!(x::AbstractDataFrame, vs)
   if !isempty(x)
      deleterows!(x, vs)
   end
end

################################################# PROPERTY ACCESSORS #######################################################

""" List the column headers in the vertex dataframe """
listvprops(g::Graph) = names(vdata(g))

""" Check if the vertex dataframe has the input column """
hasvprop(g::Graph, prop::Symbol) = haskey(vdata(g), prop)

################################################# GETVPROP #################################################################

"""
Retrieve vertex properties.

getvprop(g::Graph, v::VertexID, vprop::Symbol) -> Fetch the value of a property for vertex v
"""
getvprop(g::Graph, v::VertexID, vprop::Symbol) = vdata(g)[vprop][v]

""" getvprop(g::Graph, vs::VertexList, vprop::Symbol) -> Fetch the value of a property for v in vs """
getvprop(g::Graph, vs::VertexList, vprop::Symbol) = vdata(g)[vprop][vs]

""" getvprop(g::Graph, ::Colon, vprop::Symbol) -> Fetch the value of a property for all verices """
getvprop(g::Graph, ::Colon, vprop::Symbol) = copy(vdata(g)[vprop])

################################################# SETVPROP #################################################################

"""
Set vertex properties.

setvprop!(g::Graph, v::VertexID, val(s), vprop::Symbol) -> Set a property for v
"""
function setvprop!(g::Graph, v::VertexID, val, vprop::Symbol)
   if hasvprop(g, vprop)
      vdata(g)[vprop][v] = val
   else
      error("Property $vprop doesn't exist. Please create it on all vertices first.")
   end
end

""" setvprop!(g::Graph, vs::VertexList, val(s), vprop::Symbol) -> Set a property for v in vs """
function setvprop!(g::Graph, vs::VertexList, val, vprop::Symbol)
   if hasvprop(g, vprop)
      vdata(g)[vprop][vs] = val
   else
      error("Property $vprop doesn't exist. Please create it on all vertices first.")
   end
end

""" setvprop!(g::Graph, ::Colon, val(s), vprop::Symbol) -> Set a property for all vertices in g """
function setvprop!(g::Graph, ::Colon, vals::AbstractVector, vprop::Symbol)
   if length(vals) == nv(g)
      vdata(g)[vprop] = vals
   else
      error("Trying to set $(length(vals)) values to $(nv(g)) properties")
   end
end

setvprop!(g::Graph, ::Colon, val, vprop::Symbol) = setvprop!(g, :, fill(val, nv(g)), vprop)

################################################# SUBGRAPHING ##############################################################

""" Return a sub-dataframe for a list of vertices """
subgraph(x::AbstractDataFrame, vs::VertexList) = x[vs,:] # TODO: Do this with 1D indexing

""" Return a sub-dataframe for a list of vertices, containing a subset of properties """
subgraph(x::AbstractDataFrame, vs::VertexList, vprops::Vector{Symbol}) = x[vs,vprops] # TODO: Do this with 1D indexing


text:
using LocalBinaryPatterns: RotationEncodingTable, rotation_encoding_table

@testset "bitrotate_encoding" begin
    @testset "RotationEncodingTable" begin
        X_UInt8 = RotationEncodingTable{UInt8}(0:255)
        X_UInt16 = RotationEncodingTable{UInt16}(0:255)

        @test length(unique(X_UInt8)) == 36
        @test length(unique(X_UInt16)) == 129
        @test X_UInt8[0b1010_1001] == 0b0011_0101
        @test X_UInt8[0b1111_1101] != X_UInt16[0b0000_0000_1111_1101]

        X_UInt8_half = RotationEncodingTable{UInt8}(128:255)
        @test all(128 .<= X_UInt8_half .<= 255)
        @test length(unique(X_UInt8_half)) == 35
        @test X_UInt8_half[0b1010_1001] == 0b1001_1010

        @test_throws ArgumentError RotationEncodingTable{UInt8}(0:65535)

        X = RotationEncodingTable{UInt32}(0:2^24-1) # 64MB storage for 24 bits
        # ensure the getindex calculation finishes in finite time
        t = @elapsed sum(X)
        @test t < 10 # is expected to be smaller than 1

        # everytime it allocates new memory to store the result
        @test RotationEncodingTable{UInt8}(0:255) !== RotationEncodingTable{UInt8}(0:255)
    end

    # repeat calls are cached in global dictionary
    lookup = rotation_encoding_table(UInt8, 0:255)
    new_lookup = rotation_encoding_table(UInt8, 0:255)
    @test lookup === new_lookup
    @test lookup === rotation_encoding_table(UInt8, 8)
end


text:
level(v::BilevelVariableRef) = v.level
mylevel(v::BilevelVariableRef) = v.level
in_level(v::BilevelVariableRef, level::Level) = (
    v.level === BOTH ||
    v.level === level ||
    (v.level === DUAL_OF_LOWER && level === UPPER_ONLY))

in_upper(v::BilevelVariableRef) = in_upper(mylevel(v))
in_lower(v::BilevelVariableRef) = in_lower(mylevel(v))
upper_ref(v::BilevelVariableRef) = v.model.var_upper[v.idx]
lower_ref(v::BilevelVariableRef) = v.model.var_lower[v.idx]

function bound_ref(v::BilevelVariableRef)
    if mylevel(v) == LOWER_ONLY
        return lower_ref(v)
    elseif mylevel(v) == UPPER_ONLY || mylevel(v) == DUAL_OF_LOWER
        return upper_ref(v)
    elseif mylevel(v) == BOTH
        m = v.model
        i = v.idx
        upper_ref = m.var_upper[i]
        # this checks actual owner:
        if haskey(m.upper_to_lower_link, upper_ref)
            return upper_ref
        else
            return m.var_lower[i]
        end
    else
        error("Unknown level")
    end
end

function solver_ref(v::BilevelVariableRef)
    m = v.model
    if mylevel(v) == LOWER_ONLY
        return m.sblm_to_solver[
            m.lower_to_sblm[JuMP.index(lower_ref(v))]]
    else
        return m.sblm_to_solver[
            m.upper_to_sblm[JuMP.index(upper_ref(v))]]
    end
end

Base.broadcastable(v::BilevelVariableRef) = Ref(v)
Base.copy(v::BilevelVariableRef) = v
Base.:(==)(v::BilevelVariableRef, w::BilevelVariableRef) =
    v.model === w.model && v.idx == w.idx && v.level == w.level
JuMP.owner_model(v::BilevelVariableRef) = v.model
JuMP.isequal_canonical(v::BilevelVariableRef, w::BilevelVariableRef) = v == w
JuMP.variable_type(::AbstractBilevelModel) = BilevelVariableRef
# add in BOTH levels
function JuMP.add_variable(inner::InnerBilevelModel, v::JuMP.AbstractVariable, name::String="")
    m = bilevel_model(inner)
    m.nextvaridx += 1
    vref = BilevelVariableRef(m, m.nextvaridx, BOTH)
    # break info so that bounds go to correct level
    var_upper, var_lower = split_variable(inner, v)
    v_upper = JuMP.add_variable(m.upper, var_upper, name)
    m.var_upper[vref.idx] = v_upper
    v_lower = JuMP.add_variable(m.lower, var_lower, name)
    m.var_lower[vref.idx] = v_lower
    m.var_level[vref.idx] = BOTH
    set_link!(inner, v_upper, v_lower)
    m.variables[vref.idx] = v # save complete data
    JuMP.set_name(vref, name)
    m.var_info[vref.idx] = empty_info(v)
    vref
end
function JuMP.add_variable(single::SingleBilevelModel, v::JuMP.AbstractVariable, name::String="")
    m = bilevel_model(single)
    m.nextvaridx += 1
    vref = BilevelVariableRef(m, m.nextvaridx, level(single))
    v_level = JuMP.add_variable(mylevel_model(single), v, name)
    mylevel_var_list(single)[vref.idx] = v_level
    m.var_level[vref.idx] = level(single)
    m.variables[vref.idx] = v
    JuMP.set_name(vref, name)
    m.var_info[vref.idx] = empty_info(v)
    vref
end
function JuMP.delete(m::AbstractBilevelModel, vref::BilevelVariableRef)
    error("No deletion on bilevel models")
    # delete!(m.variables, vref.idx)
    # delete!(m.varnames, vref.idx)
end
JuMP.is_valid(m::BilevelModel, vref::BilevelVariableRef) = vref.idx in keys(m.variables)
JuMP.is_valid(m::InnerBilevelModel, vref::BilevelVariableRef) =
    JuMP.is_valid(bilevel_model(m), vref) && in_level(vref, level(m))
JuMP.num_variables(m::BilevelModel) = length(m.variables)
JuMP.num_variables(m::UpperModel) = length(m.m.var_upper)
JuMP.num_variables(m::LowerModel) = length(m.m.var_lower)
function empty_info(::JuMP.AbstractVariable)
    return VariableInfo()
end

"""
Split variable because actual owner of the variable should be the one holding the bounds.
"""
function split_variable(::UpperModel, v::JuMP.AbstractVariable)
    var_upper = v
    var_lower = JuMP.ScalarVariable(JuMP.VariableInfo(
        false, NaN,
        false, NaN,
        false, NaN,
        v.info.has_start, v.info.start,
        false, false))
    return var_upper, var_lower
end
function split_variable(::LowerModel, v::JuMP.AbstractVariable)
    var_lower = v
    var_upper = JuMP.ScalarVariable(JuMP.VariableInfo(
        false, NaN,
        false, NaN,
        false, NaN,
        v.info.has_start, v.info.start,
        false, false))
    return var_upper, var_lower
end

# Internal function
variable_info(vref::BilevelVariableRef) = vref.model.variables[vref.idx].info
function update_variable_info(vref::BilevelVariableRef, info::JuMP.VariableInfo)
    vref.model.variables[vref.idx] = JuMP.ScalarVariable(info)
end

function JuMP.has_lower_bound(vref::BilevelVariableRef)
    if mylevel(vref) == DUAL_OF_LOWER
        return !isnan(get_dual_lower_bound(get_constrain_ref(vref)))
    end
    return variable_info(vref).has_lb
end
function JuMP.lower_bound(vref::BilevelVariableRef)
    if mylevel(vref) == DUAL_OF_LOWER
        return get_dual_lower_bound(get_constrain_ref(vref))
    end
    @assert !JuMP.is_fixed(vref)
    variable_info(vref).lower_bound
end
function JuMP.set_lower_bound(vref::BilevelVariableRef, lower)
    if mylevel(vref) == DUAL_OF_LOWER
        set_dual_lower_bound(get_constrain_ref(vref), lower)
        return vref.model.variables[vref.idx]
    end
    info = variable_info(vref)
    JuMP.set_lower_bound(bound_ref(vref), lower)
    update_variable_info(vref,
                         JuMP.VariableInfo(true, lower,
                                           info.has_ub, info.upper_bound,
                                           info.has_fix, info.fixed_value,
                                           info.has_start, info.start,
                                           info.binary, info.integer))
end
function JuMP.delete_lower_bound(vref::BilevelVariableRef)
    if mylevel(vref) == DUAL_OF_LOWER
        set_dual_lower_bound(get_constrain_ref(vref), NaN)
        return vref.model.variables[vref.idx]
    end
    info = variable_info(vref)
    JuMP.delete_lower_bound(bound_ref(vref))
    update_variable_info(vref,
                         JuMP.VariableInfo(false, info.lower_bound,
                                           info.has_ub, info.upper_bound,
                                           info.has_fix, info.fixed_value,
                                           info.has_start, info.start,
                                           info.binary, info.integer))
end
function JuMP.has_upper_bound(vref::BilevelVariableRef)
    if mylevel(vref) == DUAL_OF_LOWER
        return !isnan(get_dual_upper_bound(get_constrain_ref(vref)))
    end
    return variable_info(vref).has_ub
end
function JuMP.upper_bound(vref::BilevelVariableRef)
    if mylevel(vref) == DUAL_OF_LOWER
        return get_dual_upper_bound(get_constrain_ref(vref))
    end
    @assert !JuMP.is_fixed(vref)
    variable_info(vref).upper_bound
end
function JuMP.set_upper_bound(vref::BilevelVariableRef, upper)
    if mylevel(vref) == DUAL_OF_LOWER
        set_dual_upper_bound(get_constrain_ref(vref), upper)
        return vref.model.variables[vref.idx]
    end
    info = variable_info(vref)
    JuMP.set_upper_bound(bound_ref(vref), upper)
    update_variable_info(vref,
                         JuMP.VariableInfo(info.has_lb, info.lower_bound,
                                           true, upper,
                                           info.has_fix, info.fixed_value,
                                           info.has_start, info.start,
                                           info.binary, info.integer))
end
function JuMP.delete_upper_bound(vref::BilevelVariableRef)
    if mylevel(vref) == DUAL_OF_LOWER
        set_dual_upper_bound(get_constrain_ref(vref), NaN)
        return vref.model.variables[vref.idx]
    end
    info = variable_info(vref)
    JuMP.delete_upper_bound(bound_ref(vref))
    update_variable_info(vref,
                         JuMP.VariableInfo(info.has_lb, info.lower_bound,
                                           false, info.upper_bound,
                                           info.has_fix, info.fixed_value,
                                           info.has_start, info.start,
                                           info.binary, info.integer))
end

JuMP.is_fixed(vref::BilevelVariableRef) = variable_info(vref).has_fix
JuMP.fix_value(vref::BilevelVariableRef) = variable_info(vref).fixed_value
function JuMP.fix(vref::BilevelVariableRef, value; force::Bool=false)
    if mylevel(vref) == DUAL_OF_LOWER
        error("Dual variable cannot be fixed.")
    end
    info = variable_info(vref)
    JuMP.fix(bound_ref(vref), value; force=force)
    update_variable_info(vref,
                         JuMP.VariableInfo(info.has_lb, info.lower_bound,
                                           info.has_ub, info.upper_bound,
                                           true, value,
                                           info.has_start, info.start,
                                           info.binary, info.integer))
end
function JuMP.unfix(vref::BilevelVariableRef)
    if mylevel(vref) == DUAL_OF_LOWER
        error("Dual variable cannot be fixed.")
    end
    info = variable_info(vref)
    JuMP.unfix(bound_ref(vref))
    update_variable_info(vref,
                         JuMP.VariableInfo(info.has_lb, info.lower_bound,
                                           info.has_ub, info.upper_bound,
                                           false, info.fixed_value,
                                           info.has_start, info.start,
                                           info.binary, info.integer))
end

function JuMP.start_value(vref::BilevelVariableRef)
    if mylevel(vref) == DUAL_OF_LOWER
        return get_dual_start(get_constrain_ref(vref))
    end
    variable_info(vref).start
end
function JuMP.set_start_value(vref::BilevelVariableRef, start)
    if mylevel(vref) == DUAL_OF_LOWER
        set_dual_start(get_constrain_ref(vref), start)
        return vref.model.variables[vref.idx]
    end
    info = variable_info(vref)
    in_upper(vref) && JuMP.set_start_value(upper_ref(vref), start)
    in_lower(vref) && JuMP.set_start_value(lower_ref(vref), start)
    update_variable_info(vref,
                         JuMP.VariableInfo(info.has_lb, info.lower_bound,
                                           info.has_ub, info.upper_bound,
                                           info.has_fix, info.fixed_value,
                                           true, start,
                                           info.binary, info.integer))
end

JuMP.is_binary(vref::BilevelVariableRef) = variable_info(vref).binary
function JuMP.set_binary(vref::BilevelVariableRef)
    if mylevel(vref) == DUAL_OF_LOWER
        error("Dual variable cannot be binary.")
    end
    @assert !JuMP.is_integer(vref)
    info = variable_info(vref)
    JuMP.set_binary(bound_ref(vref))
    update_variable_info(vref,
                         JuMP.VariableInfo(info.has_lb, info.lower_bound,
                                           info.has_ub, info.upper_bound,
                                           info.has_fix, info.fixed_value,
                                           info.has_start, info.start,
                                           true, info.integer))
end
function JuMP.unset_binary(vref::BilevelVariableRef)
    if mylevel(vref) == DUAL_OF_LOWER
        error("Dual variable cannot be binary.")
    end
    info = variable_info(vref)
    JuMP.unset_binary(bound_ref(vref))
    update_variable_info(vref,
                         JuMP.VariableInfo(info.has_lb, info.lower_bound,
                                           info.has_ub, info.upper_bound,
                                           info.has_fix, info.fixed_value,
                                           info.has_start, info.start,
                                           false, info.integer))
end
JuMP.is_integer(vref::BilevelVariableRef) = variable_info(vref).integer
function JuMP.set_integer(vref::BilevelVariableRef)
    if mylevel(vref) == DUAL_OF_LOWER
        error("Dual variable cannot be integer.")
    end
    @assert !JuMP.is_binary(vref)
    info = variable_info(vref)
    JuMP.set_integer(bound_ref(vref))
    update_variable_info(vref,
                         JuMP.VariableInfo(info.has_lb, info.lower_bound,
                                           info.has_ub, info.upper_bound,
                                           info.has_fix, info.fixed_value,
                                           info.has_start, info.start,
                                           info.binary, true))
end
function JuMP.unset_integer(vref::BilevelVariableRef)
    if mylevel(vref) == DUAL_OF_LOWER
        error("Dual variable cannot be integer.")
    end
    info = variable_info(vref)
    JuMP.unset_integer(bound_ref(vref))
    update_variable_info(vref,
                         JuMP.VariableInfo(info.has_lb, info.lower_bound,
                                           info.has_ub, info.upper_bound,
                                           info.has_fix, info.fixed_value,
                                           info.has_start, info.start,
                                           info.binary, false))
end

function JuMP.value(v::BilevelVariableRef; result::Int = 1)::Float64
    m = owner_model(v)
    solver = m.solver
    ref = solver_ref(v)
    return MOI.get(solver, MOI.VariablePrimal(result), ref)
end


text:
using FileIO
using Parameters

include("random_matrices.jl")
include("network_sims.jl")
include("LS_calculations.jl")

#------------------------------------------------------------------------------
# This is the part where we find which eta and network realization to

J_seeds = collect(5000:5019)
etas = collect(-1.0:0.1:1.0)

ind = parse(Int64, ARGS[1])
ps = [(eta, J_seed) for eta in etas, J_seed in J_seeds]
eta, J_seed = ps[ind]
#------------------------------------------------------------------------------

# Where the results will be saved
dir = "results1"

N = 1000
g = 10.0
s = 1.0

h0_seed = 6000
dh_seed = 7000
Q0_seed = 8000

J = generate_J(N, eta, seed=J_seed)
h = generate_nv(N, Float64, seed=h0_seed)
dh = generate_nv(N, Float64, seed=dh_seed)
Q = generate_ONB(N, Float64, seed=Q0_seed)

dt = 0.01
t_bound = 2000.0
t_ONS = 1.0
wait = 200.0

@time Lspec, max_LE = sim_lyapunov(J, g, h, dh, Q; s=s, dt=dt, t_bound=t_bound,
		t_ONS=t_ONS, wait=wait, return_max_LE=true, printstuff=false,
		p_factor=1e-3, scale_limit=100)

dimensionality = KY_attractor_dim(Lspec)
entropy_rate = KS_entropy_rate(Lspec)

fn = "Lspecs_E$(ind)_N=$N,g=$g,eta=$eta,J_seed=$J_seed.jld2"
out = Dict{String,Any}()
@pack! out = N, g, eta, s, J_seed, h0_seed, dh_seed, Q0_seed, dt, t_bound,
		t_ONS, wait, Lspec, max_LE, dimensionality, entropy_rate

save(joinpath(dir, fn), out)


text:
# * Scene

mutable struct Scene{O<:Intersectable{<:AbstractFloat}}
    objs::Vector{O}
end
# Scene(objs::Vector{O}) where {T,O<:Intersectable{T}} = Scene{T,O}(objs)
Scene{O}() where {O} = Scene(Vector{O}())

Base.eltype(::Scene{O}) where O = O

Base.show(io::IO, scene::S) where {O,S<:Scene{O}} =
    write(io, "$S with $(length(scene.objs)) $(O)s")

function Scene(::Type{O}, filename::String) where {T,O<:Intersectable{T}}
    scene = Scene{O}()
    append!(scene, FileIO.load(filename))
    scene
end

function Base.append!(scene::Scene{Triangle{T}}, mesh::HomogenousMesh) where T
    objs = map(mesh.faces) do face
        Triangle(SVector.(mesh.vertices[face]), SVector.(mesh.normals[face]))
    end
    append!(scene.objs, objs)
    scene
end

function Base.append!(scene::Scene{O}, objs::Vector{O}) where {T,O<:Intersectable{T}}
    append!(scene.objs, objs)
    scene
end

function extents(scene::Scene)
    e = extents.(scene.objs)
    join_extents(a,b) = [(min(a[i][1],b[i][1]),max(a[i][2],b[i][2])) for i=1:3]
    reduce(join_extents, e, init=e[1])
end

Images.center(scene::Scene) = mean.(extents(scene))

# ** Random scenes
rand_point(::Type{T}) where T = SVector{3}(2rand(T)-1,2rand(T)-1,2rand(T)-1)

Random.rand(::Type{Triangle{T}}) where T =
    Triangle((rand_point(T), rand_point(T), rand_point(T)))

Random.rand(::Type{Sphere{T}}) where T =
    Sphere(rand_point(T), T(0.2)*rand(T))

Random.rand(::Type{Intersectable{T}}) where T =
    rand(rand([Triangle{T}, Sphere{T}]))

Random.rand(::Type{O}, n::Integer) where {T,O<:Intersectable{T}} =
    [rand(O) for i=1:n]

Random.rand(::Type{S}, n::Integer) where {T,O<:Intersectable{T},S<:Scene{O}} =
    Scene(rand(O, n))

# ** Scene file format

Base.write(io::IOStream, v::SVector) = write(io, "(", join(string.(v), ","), ")")

function Base.write(io::IOStream, s::Sphere)
    write(io, "Sphere ", string(eltype(s)), " ")
    write(io, s.pos, " ")
    write(io, string(s.radius), "\n")
end

function Base.write(io::IOStream, t::Triangle)
    write(io, "Triangle ", string(eltype(t)), " ")
    for v in [t.o, t.e₁, t.e₂, t.na, t.nb, t.nc]
        write(io, v, " ")
    end
    write(io, "\n")
end

function FileIO.save(f::File{format"RSC"},
                     objs::AbstractVector{O}) where {T,O<:Intersectable{T}}
    open(f, "w") do file
        write(file, magic(format"RSC"))
        write(file, '\n')
        for o in objs
            write(file, o)
        end
    end
end

FileIO.save(f::File{format"RSC"}, scene::Scene) = save(f, scene.objs)

function Base.parse(V::Type{SVector{N,T}}, s::AbstractString) where {N,T<:Number}
    s[1] == '(' && s[end] == ')' || error("Malformed $(V), $(s)")
    s = s[2:end-1]
    V([parse(T,c) for c in split(s, ",")]...)
end

read_primitive(P::Type{Sphere{T}}, pos, radius) where {T<:Number} =
    P(parse(SVector{3,T}, pos), parse(T, radius))

read_primitive(P::Type{Triangle{T}}, o, e₁, e₂, na, nb, nc) where {T<:Number} =
    P(parse(SVector{3,T}, o),
      parse(SVector{3,T}, e₁),
      parse(SVector{3,T}, e₂),
      parse(SVector{3,T}, na),
      parse(SVector{3,T}, nb),
      parse(SVector{3,T}, nc))

function FileIO.load(f::File{format"RSC"})
    open(f) do file
        skipmagic(file)
        readline(stream(file))
        FileIO.load(file)
    end
end

function FileIO.load(s::Stream{format"RSC"})
    objs = Vector{Intersectable}()
    for line in eachline(stream(s))
        d = split(strip(line))
        P = eval(Symbol(d[1]))
        T = eval(Symbol(d[2]))
        push!(objs, read_primitive(P{T}, d[3:end]...))
    end
    O = promote_type(unique(typeof.(objs))...)
    Scene{O}(objs)
end

export Scene, extents, save, load


text:
const warned_once = Ref(false)
if VERSION >= v"1.5"
    function forced_depwarn(msg, sym)
        opts = Base.JLOptions()
        if !warned_once[] && !(opts.depwarn == 1)
            @warn msg
            @info """It is recommended that you fix this now to avoid breakage when a new version is released and this warning is removed.
                     Tip: to see all deprecation warnings together with code locations, launch Julia with `--depwarn=yes` and rerun your code."""
            warned_once[] = true
        else
            Base.depwarn(msg, sym)
        end
        return nothing
    end
else
    forced_depwarn(msg, sym) = Base.depwarn(msg, sym)
end

# a perhaps "permanent" deprecation
Base.@deprecate_binding permuteddimsview PermutedDimsArray


text:
using CopositiveAnalyticCenter, Test, DelimitedFiles

@testset "6 × 6 random matrices" begin
    for i in 1:10
        A = readdlm("randmat_6x6_v$i.txt")
        @test !is_completely_positive(A)
    end
end


text:
using GeoFormatTypes, Test
using GeoFormatTypes: Geom, CRS, Extended, Unknown

@testset "Test construcors" begin
    @test_throws ArgumentError ProjString("+lat_ts=56.5 +ellps=GRS80")
    @test_throws ArgumentError ProjJSON(Dict("fype" => 1))
    @test_throws ArgumentError ProjJSON("fype")
    @test_throws ArgumentError EPSG("ERROR:4326")
    @test EPSG("EPSG:4326") == EPSG(4326)
end

@testset "Test constructors" begin
    @test ProjString("+proj=test") isa ProjString
    @test ProjJSON(Dict("type" => "GeographicCRS")) isa ProjJSON
    @test ProjJSON("type: GeographicCRS") isa ProjJSON
    @test EPSG(4326) isa EPSG
    @test WellKnownText("test") isa WellKnownText{Unknown}
    @test WellKnownBinary([1, 2, 3, 4]) isa WellKnownBinary{Unknown}
    @test WellKnownText2("test") isa WellKnownText2{Unknown}
    @test ESRIWellKnownText("test") isa ESRIWellKnownText{Unknown}
    @test WellKnownText(Extended(), "test") isa WellKnownText{Extended}
    @test WellKnownBinary(Extended(), [1, 2, 3, 4]) isa WellKnownBinary{Extended}
    @test WellKnownText2(CRS(), "test") isa WellKnownText2{CRS}
    @test ESRIWellKnownText(Geom(), "test") isa ESRIWellKnownText{Geom}
    @test GML("test") isa GML{Unknown}
    @test GML(Geom(), "test") isa GML{Geom}
    @test GML(CRS(), "test") isa GML{CRS} # Probably doesn't actually exist
    @test KML("test") isa KML
    @test GeoJSON("test") isa GeoJSON
end

@testset "Test conversion to string or int" begin
    @test convert(String, ProjString("+proj=test")) == "+proj=test"
    @test convert(String, EPSG(4326)) == "EPSG:4326"
    @test convert(Int, EPSG(4326)) == 4326
    @test convert(String, WellKnownText("test")) == "test"
    @test convert(String, WellKnownText2("test")) == "test"
    @test convert(String, ESRIWellKnownText("test")) == "test"
    @test convert(String, GML("test")) == "test"
    @test convert(String, KML("test")) == "test"
    @test convert(String, GeoJSON("test")) == "test"
end

# `convert` placeholder methods
Base.convert(target::Type{<:GeoFormat}, mode::Union{Geom,Type{Geom}}, source::GeoFormat; kwargs...) =
    (:geom, kwargs...)
Base.convert(target::Type{<:GeoFormat}, mode::Union{CRS,Type{CRS}}, source::GeoFormat; kwargs...) =
    (:crs, kwargs...)

@testset "Test convert mode allocation" begin
    @testset "Test identical type is passed through unchanged" begin
        @test convert(WellKnownText, WellKnownText(Extended(), "test")) == WellKnownText(Extended(), "test")
        @test convert(ProjString, ProjString("+proj=test")) == ProjString("+proj=test")
    end
    @testset "Test conversions are assigned to crs or geom correctly" begin
        @test convert(WellKnownText, WellKnownText2(CRS(), "test")) == (:crs,)
        @test convert(WellKnownText2, WellKnownText(CRS(), "test")) == (:crs,)
        @test convert(WellKnownBinary, WellKnownText(CRS(), "test")) == (:crs,)
        @test convert(ProjString, WellKnownText(CRS(), "test")) == (:crs,)
        @test convert(EPSG, ProjString("+proj=test")) == (:crs,)
        @test convert(CoordSys, ProjString("+proj=test")) == (:crs,)

        @test convert(GeoJSON, WellKnownText(Geom(), "test")) == (:geom,)
        @test convert(KML, WellKnownText(Geom(), "test")) == (:geom,)
        @test convert(GML, WellKnownText(Geom(), "test")) == (:geom,)
        @test convert(ESRIWellKnownText, WellKnownText(Geom(), "test")) == (:geom,)
        @test convert(WellKnownBinary, WellKnownText(Geom(), "test")) == (:geom,)
        @test convert(WellKnownText2, WellKnownText(Geom(), "test")) == (:geom,)
        @test convert(WellKnownText2, WellKnownText(Geom(), "test")) == (:geom,)
        @test convert(WellKnownText, WellKnownText2(Geom(), "test")) == (:geom,)

        @test convert(GeoJSON, WellKnownText(Extended(), "test")) == (:geom,)
        @test convert(KML, WellKnownText(Extended(), "test")) == (:geom,)
        @test convert(GML, WellKnownText(Extended(), "test")) == (:geom,)
        @test convert(ESRIWellKnownText, WellKnownText(Extended(), "test")) == (:geom,)
        @test convert(WellKnownBinary, WellKnownText(Extended(), "test")) == (:geom,)
        @test convert(WellKnownText2, WellKnownText(Extended(), "test")) == (:geom,)
        @test convert(WellKnownText2, WellKnownText(Extended(), "test")) == (:geom,)
        @test convert(WellKnownText, WellKnownText2(Extended(), "test")) == (:geom,)

        @test convert(GeoJSON, WellKnownText(Unknown(), "test")) == (:geom,)
        @test convert(KML, WellKnownText(Unknown(), "test")) == (:geom,)
        @test convert(GML, WellKnownText(Unknown(), "test")) == (:geom,)
        @test convert(ESRIWellKnownText, WellKnownText(Unknown(), "test")) == (:geom,)
        @test convert(WellKnownBinary, WellKnownText(Unknown(), "test")) == (:geom,)
        @test convert(WellKnownText2, WellKnownText(Unknown(), "test")) == (:geom,)
        @test convert(WellKnownText2, WellKnownText(Unknown(), "test")) == (:geom,)
        @test convert(WellKnownText, WellKnownText2(Unknown(), "test")) == (:geom,)
    end
    @testset "Test kargs pass through convert" begin
        @test convert(WellKnownText, WellKnownText2(CRS(), "test"); order=:trad) == (:crs, :order => :trad,)
        @test convert(GML, WellKnownText(Extended(), "test"); order=:custom) == (:geom, :order => :custom)
    end
    @testset "Test conversions that are not possible throw an error" begin
        @test_throws ArgumentError convert(KML, ProjString("+proj=test"))
        @test_throws ArgumentError convert(GeoJSON, ProjString("+proj=test"))
        @test_throws ArgumentError convert(ProjString, WellKnownText(Geom(), "test"))
        @test_throws ArgumentError convert(CoordSys, WellKnownText(Geom(), "test"))
        @test_throws ArgumentError convert(EPSG, WellKnownText(Geom(), "test"))
    end
end


text:
using ITensors
using ITensorTDVP
using Random
using Test

@testset "DMRG-X" begin
  function heisenberg(n; h=zeros(n))
    os = OpSum()
    for j in 1:(n - 1)
      os += 0.5, "S+", j, "S-", j + 1
      os += 0.5, "S-", j, "S+", j + 1
      os += "Sz", j, "Sz", j + 1
    end
    for j in 1:n
      if h[j] ≠ 0
        os -= h[j], "Sz", j
      end
    end
    return os
  end

  n = 10
  s = siteinds("S=1/2", n)

  Random.seed!(12)

  W = 12
  # Random fields h ∈ [-W, W]
  h = W * (2 * rand(n) .- 1)
  H = MPO(heisenberg(n; h), s)

  initstate = rand(["↑", "↓"], n)
  ψ = MPS(s, initstate)

  dmrg_x_kwargs = (
    nsweeps=20, reverse_step=false, normalize=true, maxdim=20, cutoff=1e-10, outputlevel=0
  )

  ϕ = dmrg_x(ProjMPO(H), ψ; dmrg_x_kwargs...)

  @test inner(ψ', H, ψ) / inner(ψ, ψ) ≈ inner(ϕ', H, ϕ) / inner(ϕ, ϕ) rtol = 1e-1
  @test inner(H, ψ, H, ψ) ≉ inner(ψ', H, ψ)^2 atol = 1e-1
  @test inner(H, ϕ, H, ϕ) ≈ inner(ϕ', H, ϕ)^2 atol = 1e-7
end

nothing


text:
module AD

using ..NiLang
using NiLangCore
using MatchCore, TupleTools

import ..NiLang: ROT, IROT, SWAP,
    chfield, value, NoGrad, INC, DEC, HADAMARD,
    AddConst, SubConst, NEG, INV
using NiLangCore: default_constructor

export GVar, grad, Loss, NoGrad, @nograd

include("vars.jl")
include("stack.jl")
include("gradfunc.jl")
include("checks.jl")

include("instructs.jl")
include("ulog.jl")
include("jacobian.jl")
include("hessian_backback.jl")
include("complex.jl")

if Base.VERSION >= v"1.4.2"
    include("precompile.jl")
    _precompile_()
end

end


text:
type SimpleOG
    B :: Int64
    M :: Int64
    alpha :: Float64
    beta :: Float64
    R :: Array{Float64}
    Q :: Array{Float64}
end

function SimpleOG(;B=10, M=5, alpha=0.5, beta=0.9)

    u(c) = c^alpha
    n = B + M + 1
    m = M + 1

    R = Array(Float64,n,m)
    Q = zeros(Float64,n,m,n)

    for a in 0:M
        Q[:, a + 1, (a:(a + B)) + 1] = 1 / (B + 1)
        for s in 0:(B + M)
            R[s + 1, a + 1] = a<=s ? u(s - a) : -Inf
        end
    end

    return SimpleOG(B, M, alpha, beta, R, Q)
end

max_stars_repo_name: mwhchen/quantecon

max_issues_repo_name: mwhchen/quantecon

max_forks_repo_name: mwhchen/quantecon


text:
import Base: IteratorSize, HasLength, length

"""
    InstantTrigger(dt::DateTime)

A trigger which should trigger job at a given instant (a given `DateTime` for example)
"""
struct InstantTrigger <: AbstractTrigger
    dt_fire_at::DateTime
end

IteratorSize(trigger::InstantTrigger) = HasLength()
length(trigger::InstantTrigger) = 1

"""
    Trigger(dt::DateTime)

Return an `InstantTrigger` which should trigger job at a given `DateTime` `dt`
"""
Trigger(dt::DateTime) = InstantTrigger(dt)

"""
    Trigger(d::Date)

Return an `InstantTrigger` which should trigger job at a given `Date` `d`(at midnight)
"""
Trigger(d::Date) = InstantTrigger(DateTime(d))

"""
    get_next_dt_fire(trigger, dt_previous_fire, dt_now)

Return instant at which a job should be run (given current instant `dt_now` and instant when job was previously run (`dt_previous_fire`)
"""
function get_next_dt_fire(trigger::InstantTrigger, dt_previous_fire, dt_now)
    if dt_previous_fire == DateTime(0)
        trigger.dt_fire_at
    else
        DateTime(0)
    end
end


text:
abstract type AbstractImplicitSurface end

mutable struct ImplicitResult{T<:CSGBuilding.AbstractImplicitSurface, F<:Real}
    surface::T
    value::F
    signint::Int
end

value(result::ImplicitResult) = result.value

function normal(result::ImplicitResult, coords)
    return result.signint * normal(result.surface, coords)
end

Base.show(io::IO, result::ImplicitResult) =
    print(io, "ImplicitResult: $(result.surface), $(result.value)")

Base.min(x::ImplicitResult, y::ImplicitResult) = ifelse(isless(x.value, y.value), x, y)

Base.max(x::ImplicitResult, y::ImplicitResult) = ifelse(isless(y.value, x.value), x, y)

## sphere

struct ImplicitSphere{T<:Real} <: AbstractImplicitSurface
    center::SVector{3,T}
    radius::T
end

function ImplicitSphere(center::AbstractArray, radius)
    T = promote_type(eltype(center), typeof(radius))
    nc = convert(SVector{3, T}, center)
    return ImplicitSphere(nc, convert(T, radius))
end

Base.show(io::IO, surface::ImplicitSphere) = print(io, "ImplSphere: R$(surface.radius)")
_name(surface::ImplicitSphere) = "Sphere"

function evaluate(surface::ImplicitSphere, coords)
    val = norm(coords-surface.center) - surface.radius
    return ImplicitResult(surface, val, 1)
end

function normal(surface::ImplicitSphere, coords)
    dist = coords-surface.center
    if ! isapprox(norm(dist), 0)
        return normalize(dist)
    else
        return convert(typeof(coords), [0,0,0])
    end
end

## plane

struct ImplicitPlane{T<:Real} <: AbstractImplicitSurface
    point::SVector{3,T}
    normal::SVector{3,T}
end

function ImplicitPlane(point::AbstractArray, normal::AbstractArray)
    T = promote_type(eltype(point), eltype(normal))
    np = convert(SVector{3, T}, point)
    nn = convert(SVector{3, T}, normal)
    return ImplicitPlane(np, nn)
end

Base.show(io::IO, surface::ImplicitPlane) = print(io, "ImplPlane: n$(surface.normal)")
_name(surface::ImplicitPlane) = "Plane"

function evaluate(surface::ImplicitPlane, coords)
    val = dot(coords-surface.point, surface.normal)
    return ImplicitResult(surface, val, 1)
end

normal(surface::ImplicitPlane, coords) = surface.normal

"""
    issame(x::AbstractImplicitSurface, y::AbstractImplicitSurface, params)

For any non-`ImplicitPlane` return false.
"""
issame(x::AbstractImplicitSurface, y::AbstractImplicitSurface, params) = false

"""
    issame(pl1::ImplicitPlane, pl2::ImplicitPlane, params)

Are the two planes the same? Parameters `α_conv, ϵ_conv`.
"""
function issame(pl1::ImplicitPlane, pl2::ImplicitPlane, params)
    @unpack α_conv, ϵ_conv = params
    abs(dot(pl1.normal, pl2.normal)) < α_conv && return false
    d = pl1.point-pl2.point
    abs(dot(d, pl1.normal)) > ϵ_conv && return false
    return true
end

## cylinder

struct ImplicitCylinder{T<:Real} <: AbstractImplicitSurface
    axis::SVector{3,T}
    center::SVector{3,T}
    radius::T
end

function ImplicitCylinder(axis::AbstractArray, center::AbstractArray, radius)
    T = promote_type(eltype(axis), eltype(center), typeof(radius))
    na = convert(SVector{3, T}, axis)
    nc = convert(SVector{3, T}, center)
    return ImplicitCylinder(na, nc, convert(T, radius))
end

"""
    vectorfromline(p, d, q)

Compute the vector directing from a line to `q` point.
The line is defined by one point `p` and it's direction vector `d`.
"""
vectorfromline(p, d, q) = q-p-d*dot(q-p, d)

function evaluate(surface::ImplicitCylinder, coords)
    c = surface.center
    a = surface.axis
    val = norm(vectorfromline(c, a, coords)) - surface.radius
    return ImplicitResult(surface, val, 1)
end

function normal(surface::ImplicitCylinder, coords)
    c = surface.center
    a = surface.axis
    tnorm = vectorfromline(c, a, coords)
    if ! isapprox(norm(tnorm), 0)
        return normalize(tnorm)
    else
        return convert(typeof(coords), [0,0,0])
    end
end

function Base.show(io::IO, surface::ImplicitCylinder)
    return print(io, "ImplCylinder: R$(surface.radius)")
end
_name(surface::ImplicitCylinder) = "Cylinder"

## cone

struct ImplicitCone{T<:Real} <: AbstractImplicitSurface
    apex::SVector{3,T}
    axis::SVector{3,T}
    opang::T
end

function ImplicitCone(apex::AbstractArray, axis::AbstractArray, opang)
    T = promote_type(eltype(apex), eltype(axis), typeof(opang))
    nap = convert(SVector{3, T}, apex)
    nax = convert(SVector{3, T}, axis)
    return ImplicitCone(nap, nax, convert(T, radius))
end

function evaluate(surface::ImplicitCone, coords)
    d, _ = project2cone(surface, coords)
    return ImplicitResult(surface, d, 1)
end

function normal(surface::ImplicitCone, coords)
    _, n = project2cone(surface, coords)
    return n
end

function Base.show(io::IO, surface::ImplicitCone)
    return print(io, "ImplCone: op:$(rad2deg(surface.opang)°)")
end
_name(surface::ImplicitCone) = "Cone"

## translational

struct ImplicitTranslational <: AbstractImplicitSurface
    coordframe
    contour
    # center of gravity
    center
    # normal of the contour is parallel to the direction
    # towards the center of the contour?
    # == should flip the computed normals to direct outwards?
    # this is used in e.g. CSGBuilding
    # true means, that the computed normals must be turned to direct outside
    outwards::Int
    # should the computed normal be flipped to match the measured points
    # this is used in this package to ensure that in/outwards is correct
    # true means that computed normals must be turned to match the measured points
    flipnormal::Int
end

function evaluate(surface::ImplicitTranslational, coords)
    proj = project2sketchplane([coords], surface.coordframe)
    d, _, _ = dn2shape_outw(proj[1], surface)
    return ImplicitResult(surface, d, 1)
end

function normal(surface::ImplicitTranslational, coords)
    proj = project2sketchplane([coords], surface.coordframe)
    _, n, _ = dn2shape_outw(proj[1], surface)
    return n
end

function Base.show(io::IO, surface::ImplicitTranslational)
    return print(io, "ImplTransl")
end
_name(surface::ImplicitTranslational) = "Transl"


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule CUTENSOR_CUDA112_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("CUTENSOR_CUDA112")
JLLWrappers.@generate_main_file("CUTENSOR_CUDA112", UUID("973b338e-f8c7-52f3-bf7b-688ce4a92402"))
end  # module CUTENSOR_CUDA112_jll


text:
using Tikz

colors = Dict( "redred" => [150.,50.,50.] )

style = tikzStyle( "style1", lw = 2.0 )

obj = tikzObj( "example2.tex", standalone = true, colors = colors, styles = [style] )

x = linspace( 0., 3., 50).^1.
y = sin.( 2*π * x) .* cos.( 1./6 * π * x )

axis1 = axisObj( obj )

plot( obj, axis1, x, y )#, smooth=false, linestyle="--", color="blueblue", style=style )

closeObj( obj, axis1 )

closeObj( obj )

max_stars_repo_name: mopg/julia-tikz

max_issues_repo_name: mopg/julia-tikz

max_forks_repo_name: mopg/julia-tikz


text:
puzzleinput = map(collect, readlines("inputs/day11.txt"))
puzzleinput = hcat(puzzleinput...)
puzzleinput = map(c -> c == 'L' ? 1 : 0, puzzleinput)
puzzleinput = permutedims(puzzleinput)

# First Half

shiftl(a) = hcat(a[:, 2:end], fill(0, (size(a)[1], 1)))
shiftr(a) = hcat(fill(0, (size(a)[1], 1)), a[:, 1:end - 1])
shiftu(a) = vcat(a[2:end, :], fill(0, (1, size(a)[2])))
shiftd(a) = vcat(fill(0, (1, size(a)[2])), a[1:end - 1, :])
# Shift left, right, up and down an array a filling the new cells with 0

function newoccupied1(o, l)
    # Starting from occupied o, returns new occupied after one round of the rules
    # The seat layout is l

    ao =
        shiftl(o) +
        shiftr(o) +
        shiftu(o) +
        shiftd(o) +
        shiftl(shiftd(o)) +
        shiftr(shiftd(o)) +
        shiftr(shiftu(o)) +
        shiftl(shiftu(o))
    # Calculates how many adjacent are occupied for all the places

    function rule(i)
        if o[i] == 0 && l[i] == 1 && ao[i] == 0
            # If a seat is empty (0) and there are no occupied seats adjacent to it,
            # the seat becomes occupied.
            1
        elseif o[i] == 1 && ao[i] > 3
            # If a seat is occupied (1) and four or more seats adjacent to it are
            # also occupied, the seat becomes empty.
            0
        else
            # Otherwise, the seat's state does not change.
            o[i]
        end
    end

    map(rule, CartesianIndices(l))
    # Apply the rule to all the places
end

occupied = fill(0, size(puzzleinput))
while true
    no = newoccupied1(occupied, puzzleinput)
    if no == occupied
        # Exit when further applications of the rules cause no seats to change state
        break
    end
    global occupied = no
end

sum(occupied)

# Second Half

function newoccupied2(o, l)
    # Starting from occupied o, returns new occupied after one round of the rules
    # The seat layout is l
    function sees(i)
        # Calculates how many first seats in each direction from place i
        # are occupied
        ij = [0; 0]
        ij[1] = i[1]
        ij[2] = i[2]
        m = [size(o)[1]; size(o)[2]]
        occ = 0
        for di = -1:1
            for dj = -1:1
                if di == 0 && dj == 0
                    continue
                end
                d = [di; dj]
                dij = ij + [di; dj]
                while dij[1] > 0 && dij[2] > 0 && dij[1] <= m[1] && dij[2] <= m[2]
                    if l[dij[1], dij[2]] == 1 # There is a seat in the layout
                        if o[dij[1], dij[2]] == 1 # And is occupied
                            occ = occ + 1
                        end
                        break
                    end
                    dij = dij + d
                end
            end
        end
        occ
    end

    ao = map(sees, CartesianIndices(l))
    # Calculates how many first seat in each direction are occupied
    # for all the places

    function rule(i)
        if o[i] == 0 && l[i] == 1 && ao[i] == 0
            # If a seat is empty (0) and there are no occupied seats adjacent to it,
            # the seat becomes occupied.
            1
        elseif o[i] == 1 && ao[i] > 4
            # If a seat is occupied (1) and five or more visible occupied seats are
            # also occupied, the seat becomes empty.
            0
        else
            # Otherwise, the seat's state does not change.
            o[i]
        end
    end

    map(rule, CartesianIndices(l))
    # Apply the rule to all the places
end

occupied = fill(0, size(puzzleinput))
while true
    no = newoccupied2(occupied, puzzleinput)
    if no == occupied
        # Exit when further applications of the rules cause no seats to change state
        break
    end
    global occupied = no
end

sum(occupied)


text:
module GadgetronExamples

include("recon_acquisitions.jl")
end


text:
using Test

include("word-count.jl")

@testset "no words" begin
    @test wordcount(" .\n,\t!^&*()~@#\$%{}[]:;'/<>") == Dict()
end

@testset "count one word" begin
    @test wordcount("word") == Dict("word" => 1)
end

@testset "count one of each word" begin
    @test wordcount("one of each") == Dict("one" => 1, "of" => 1, "each" => 1)
end

@testset "multiple occurrences of a word" begin
    @test wordcount("one fish two fish red fish blue fish") == Dict("one" => 1, "fish" => 4, "two" => 1, "red" => 1, "blue" => 1)
end

@testset "handles cramped lists" begin
    @test wordcount("one,two,three") == Dict("one" => 1, "two" => 1, "three" => 1)
end

@testset "handles expanded lists" begin
    @test wordcount("one,\ntwo,\nthree") == Dict("one" => 1, "two" => 1, "three" => 1)
end

@testset "ignore punctuation" begin
    @test wordcount("car: carpet as java: javascript!!&@\$%^&") == Dict("car" => 1, "carpet" => 1, "as" => 1, "java" => 1, "javascript" => 1)
end

@testset "include numbers" begin
    @test wordcount("testing, 1, 2 testing") == Dict("testing" => 2, "1" => 1, "2" => 1)
end

@testset "normalize case" begin
    @test wordcount("go Go GO Stop stop") == Dict("go" => 3, "stop" => 2)
end

@testset "with apostrophes" begin
    @test wordcount("First: don't laugh. Then: don't cry.") == Dict("first" => 1, "don't" => 2, "laugh" => 1, "then" => 1, "cry" => 1)
end

@testset "with quotations" begin
    @test wordcount("Joe can't tell between 'large' and large.") == Dict("joe" => 1, "can't" => 1, "tell" => 1, "between" => 1, "large" => 2, "and" => 1)
end

@testset "substrings from the beginning" begin
    @test wordcount("Joe can't tell between app, apple and a.") == Dict("joe" => 1, "can't" => 1, "tell" => 1, "between" => 1, "app" => 1, "apple" => 1, "and" => 1, "a" => 1)
end

@testset "multiple spaces not detected as a word" begin
    @test wordcount(" multiple   whitespaces") == Dict("multiple" => 1, "whitespaces" => 1)
end

@testset "alternating word separators not detected as a word" begin
    @test wordcount(",\n,one,\n ,two \n 'three'") == Dict("one" => 1, "two" => 1, "three" => 1)
end


text:
"""
    LD(A::Diagonal)
    LD(A::HBlikDiag)
    LD(A::DenseMatrix)

Return `log(det(tril(A)))` evaluated in place.
"""
LD(d::Diagonal{T}) where {T<:Number} = sum(log, d.diag)

function LD(d::UniformBlockDiagonal{T}) where {T}
    dat = d.data
    sum(log, dat[j, j, k] for j in axes(dat, 2), k in axes(dat, 3))
end

LD(d::DenseMatrix{T}) where {T} = sum(log, d[k] for k in diagind(d))

"""
    logdet(m::LinearMixedModel)

Return the value of `log(det(Λ'Z'ZΛ + I)) + m.optsum.REML * log(det(LX*LX'))`
evaluated in place.

Here LX is the diagonal term corresponding to the fixed-effects in the blocked
lower Cholesky factor.
"""
function LinearAlgebra.logdet(m::LinearMixedModel{T}) where {T}
    L = m.L
    nre = m.dims.nretrms
    s = LD(getblock(L, 1,1))
    for j in 2:nre
        s += LD(getblock(L, j, j))
    end
    if m.optsum.REML
        feindex = nre + 1
        s += LD(getblock(L, feindex, feindex))
    end
    s + s
end