text:
module LDB
using Base: AbstractFloat
export
    # energy map
    EnergyMap,
    TimeFrequency,
    ProbabilityDensity,
    Signatures,
    energy_map,
    # discriminant measures
    DiscriminantMeasure,
    ProbabilityDensityDM,
    SignaturesDM,
    AsymmetricRelativeEntropy,
    SymmetricRelativeEntropy,
    HellingerDistance,
    LpDistance,
    EarthMoverDistance,
    discriminant_measure,
    # discriminant power
    DiscriminantPower,
    BasisDiscriminantMeasure,
    FishersClassSeparability,
    RobustFishersClassSeparability,
    discriminant_power,
    # local discriminant basis
    LocalDiscriminantBasis,
    fit!,
    fit_transform,
    transform,
    inverse_transform,
    change_nfeatures

using AverageShiftedHistograms,
      LinearAlgebra,
      Wavelets,
      Parameters,
      Statistics,
      StatsBase

using ..Utils,
      ..DWT,
      ..BestBasis

## ENERGY MAPS
"""
Energy map for Local Discriminant Basis. Current available types are:
- [`TimeFrequency`](@ref)
- [`ProbabilityDensity`](@ref)
"""
abstract type EnergyMap end

@doc raw"""
    TimeFrequency <: EnergyMap

An energy map based on time frequencies, a measure based on the differences of
derived quantities from projection ``Z_i``, such as mean class energies or
cumulants.

**See also:** [`EnergyMap`](@ref), [`ProbabilityDensity`](@ref),
    [`Signatures`](@ref)
"""
struct TimeFrequency <: EnergyMap end

@doc raw"""
    ProbabilityDensity <: EnergyMap

An energy map based on probability density, a measure based on the differences
among the pdfs of ``Z_i``. Since we do not know the true density functions of
the coefficients, the PDFs are estimated using the Average Shifted Histogram
(ASH).

**See also:** [`EnergyMap`](@ref), [`TimeFrequency`](@ref), [`Signatures`](@ref)
"""
struct ProbabilityDensity <: EnergyMap end

@doc raw"""
    Signatures <: EnergyMap

An energy map based on signatures, a measure that uses the Earth Mover's
Distance (EMD) to compute the discriminating  power of a coordinate. Signatures
provide us with a fully data-driven representation, which can be efficiently
used with EMD. This representation is more efficient than a histogram and is
able to represent complex data structure with fewer samples.

Here, a signature for the coefficients in the ``j``-th level, ``k``-th node,
``l``-th index of class ``c`` is defined as

``s_{j,k,l}^{(c)} = \{(\alpha_{i;j,k,l}^{(c)}, w_{i;j,k,l}^{(c)})\}_{i=1}^{N_c}``

where ``\alpha_{i;j,k,l}^{(c)}`` and ``w_{i;j,k,l}^{(c)}`` are the expansion
coefficients and weights at location ``(j,k,l)`` for signal ``i`` of class ``c``
respectively. Currently, the two valid types of weights are `:equal` and `:pdf`.

# Argumemts
- `weight::Symbol`: Type of weight to be used to compute ``w_{i;j,k,l}^{(c)}``.
    Available methods are `:equal` and `pdf`. Default is set to `:equal`.

**See also:** [`EnergyMap`](@ref), [`TimeFrequency`](@ref),
    [`ProbabilityDensity`](@ref)
"""
struct Signatures <: EnergyMap
    weight::Symbol
    Signatures(weight=:equal) = weight ∈ [:equal, :pdf] ? new(weight) :
        error("Invalid weight type. Valid weight types are :equal and :pdf.")
end

"""
    energy_map(Xw, y, method)

Returns the Time Frequency Energy map or the Probability Density Energy map
depending on the input `method` (`TimeFrequency()` or `ProbabilityDensity()`).

**See also:** [`EnergyMap`](@ref). [`TimeFrequency`](@ref),
    [`ProbabilityDensity`](@ref)
"""
function energy_map(Xw::AbstractArray{S,3}, y::AbstractVector{T},
        method::TimeFrequency) where {S<:Number, T}

    # basic summary of data
    c = unique(y)       # unique classes
    nc = length(c)      # number of classes
    Ny = length(y)
    n, L, Nx = size(Xw)

    # parameter checking
    @assert Nx == Ny
    @assert nc > 1
    @assert isdyadic(n)
    @assert 1 <= L-1 <= maxtransformlevels(n)

    # construct normalized energy map for each class
    Γ = Array{Float64, 3}(undef, (n,L,nc))
    @inbounds begin
        for (i,cᵢ) in enumerate(c)
            idx = findall(yᵢ -> yᵢ==cᵢ, y)
            x = Xw[:,1,idx]
            norm_sum = sum(mapslices(xᵢ -> norm(xᵢ,2)^2, x, dims = 1))
            en = sum(Xw[:,:,idx].^2, dims=3)
            Γ[:,:,i] = en / norm_sum
        end
    end

    return Γ
end

function energy_map(Xw::AbstractArray{S,3}, y::AbstractVector{T},
        method::ProbabilityDensity) where {S<:Number, T}

    # basic summary of data
    c = unique(y)       # unique classes
    nc = length(c)      # number of classes
    Ny = length(y)
    n, L, Nx = size(Xw)

    # parameter checking
    @assert Nx == Ny
    @assert nc > 1
    @assert isdyadic(n)
    @assert 1 <= L-1 <= maxtransformlevels(n)

    # construct empirical probability density for each coefficent of each class
    nbins = ceil(Int, (30*Nx)^(1/5)) # number of bins/histogram
    mbins = ceil(Int, 100/nbins)     # number of histograms M/nbins, M=100 is arbitrary
    Γ = Array{Float64,4}(undef, (n, L, (nbins+1)*mbins, nc))
    @inbounds begin
        for (i,cᵢ) in enumerate(c)      # iterate over each class
            idx = findall(yᵢ -> yᵢ==cᵢ, y)
            xw = Xw[:,:,idx]            # wavelet packet for class cᵢ
            for j in axes(xw,1)
                for k in axes(xw,2)
                    z = @view Xw[j,k,:] # coefficients at (j,k) for all signals
                    zᵢ = xw[j,k,:]      # coefficients at (j,k) for cᵢ signals

                    # ash parameter setup
                    σ = std(z)
                    s = 0.5
                    δ = (maximum(z)-minimum(z)+σ)/((nbins+1)*mbins-1)
                    rng = range(minimum(z)-s*σ, step=δ, length=(nbins+1)*mbins)

                    # empirical pdf
                    epdf = ash(zᵢ, rng=rng, m=mbins, kernel=Kernels.triangular)
                    _, Γ[j,k,:,i] = xy(epdf)
                end
            end
        end
    end

    return Γ
end

function energy_map(Xw::AbstractArray{S,3}, y::AbstractVector{T},
        method::Signatures) where {S<:Number, T}

    # basic summary of data
    c = unique(y)       # unique classes
    nc = length(c)      # number of classes
    Ny = length(y)
    n, L, Nx = size(Xw)

    # parameter checking
    @assert Nx == Ny
    @assert nc > 1
    @assert isdyadic(n)
    @assert 1 <= L-1 <= maxtransformlevels(n)

    # form signatures in a structure of a named tuple
    Γ = method.weight==:equal ?
        Array{NamedTuple{(:coef, :weight), Tuple{Array{S}, Float64}},1}(undef, nc) :      # equal weights
        Array{NamedTuple{(:coef, :weight), Tuple{Array{S}, Array{Float64}}},1}(undef, nc) # pdf-based weights
    for (i, cᵢ) in enumerate(c)
        idx = findall(yᵢ -> yᵢ==cᵢ, y)
        xw = Xw[:,:,idx]            # wavelet packet for class cᵢ
        if method.weight == :equal
            Nc = length(idx)
            w = 1/Nc
        else
            Nc = length(idx)
            nbins = ceil(Int, (30*Nx)^(1/5)) # number of bins/histogram
            mbins = ceil(Int, 100/nbins)     # number of histograms M/nbins, M=100 is arbitrary

            # compute weights
            w = Array{Float64,3}(undef, (n,L,Nc))
            for j in axes(xw,1)
                for k in axes(xw,2)
                    z = @view xw[j,k,:] # coefficients at (j,k) for cᵢ signals

                    # ash parameter setup
                    σ = std(z)
                    s = 0.5
                    δ = (maximum(z)-minimum(z)+σ)/((nbins+1)*mbins-1)
                    rng = range(minimum(z)-s*σ, step=δ, length=(nbins+1)*mbins)

                    # empirical pdf
                    epdf = ash(z, rng=rng, m=mbins, kernel=Kernels.triangular)
                    for l in 1:Nc
                        w[j,k,l] = AverageShiftedHistograms.pdf(epdf, z[l])
                    end
                end
            end
        end
        Γ[i] = (coef = xw, weight = w)
    end
    return Γ
end

## DISCRIMINANT MEASURES
"""
Discriminant measure for Local Discriminant Basis. Current available subtypes
are:
- [`ProbabilityDensityDM`](@ref)
- [`SignaturesDM`](@ref)
"""
abstract type DiscriminantMeasure end

"""
Discriminant measure for Probability Density and Time Frequency based energy
maps. Current available measures are:
- [`AsymmetricRelativeEntropy`](@ref)
- [`SymmetricRelativeEntropy`](@ref)
- [`LpDistance`](@ref)
- [`HellingerDistance`](@ref)
"""
abstract type ProbabilityDensityDM <: DiscriminantMeasure end

"""
Discriminant measure for Signatures based energy maps. Current available
measures are:
- [`EarthMoverDistance`](@ref)
"""
abstract type SignaturesDM <: DiscriminantMeasure end

@doc raw"""
    AsymmetricRelativeEntropy <: ProbabilityDensityDM

Asymmetric Relative Entropy discriminant measure for the Probability Density and
Time Frequency based energy maps. This measure is also known as cross entropy
and Kullback-Leibler divergence.

Equation: ``D(p,q) = \sum p(x) \log \frac{p(x)}{q(x)}``
"""
struct AsymmetricRelativeEntropy <: ProbabilityDensityDM end

@doc raw"""
    SymmetricRelativeEntropy <: ProbabilityDensityDM

Symmetric Relative Entropy discriminant measure for the Probability Density and
Time Frequency energy maps. Similar idea to the Asymmetric Relative Entropy, but
this aims to make the measure more symmetric.

Equation: Denote the Asymmetric Relative Entropy as ``D_A(p,q)``, then

``D(p,q) = D_A(p,q) + D_A(q,p) = \sum p(x) \log \frac{p(x)}{q(x)} + q(x) \log \frac{q(x)}{p(x)}``
"""
struct SymmetricRelativeEntropy <: ProbabilityDensityDM end

@doc raw"""
    LpDistance <: ProbabilityDensityDM

``\ell^p`` Distance discriminant measure for the Probability Density and Time
Frequency based energy maps. The default ``p`` value is set to 2.

Equation: ``W(q,r) = ||q-r||_p^p = \sum_{i=1}^n (q_i - r_i)^p``
"""
@with_kw struct LpDistance <: ProbabilityDensityDM
    p::Number = 2
end

@doc raw"""
    HellingerDistance <: ProbabilityDensityDM

Hellinger Distance discriminant measure for the Probability Density energy
map.

Equation: ``H(p,q) = \sum_{i=1}^n (\sqrt{p_i} - \sqrt{q_i})^2``
"""
struct HellingerDistance <: ProbabilityDensityDM end

@doc raw"""
    EarthMoverDistance <: SignaturesDM

Earth Mover Distance discriminant measure for the Signatures energy map.

Equation:
``E(P,Q) = \frac{\sum_{k=1}^{m+n+1} |\hat p_k - \hat q_k| (r_{k+1} - r_k)}{w_\Sigma}``

where ``r_1, r_2, \ldots, r_{m+n}`` is the sorted list of ``p_1, \ldots, p_m,
q_1, \ldots, q_n`` and ``\hat p_k = \sum_{p_i \leq r_k} w_{p_i}``,
``\hat q_k = \sum_{q)i \leq r_k} w_{q_i}``.
"""
struct EarthMoverDistance <: SignaturesDM end

"""
    discriminant_measure(Γ, dm)

Returns the discriminant measure of each node calculated from the energy maps.
"""
function discriminant_measure(Γ::AbstractArray{T},
                              dm::ProbabilityDensityDM) where T<:Number

    # basic summary of data
    @assert 3 <= ndims(Γ) <= 4
    n = size(Γ,1)
    L = size(Γ,2)
    C = size(Γ)[end]
    @assert C > 1       # ensure more than 1 class

    D = zeros(Float64, (n,L))
    @inbounds begin
        for i in 1:(C-1)
            for j in (i+1):C
                if ndims(Γ) == 3
                    D += discriminant_measure(Γ[:,:,i], Γ[:,:,j], dm)
                else
                    D += discriminant_measure(Γ[:,:,:,i], Γ[:,:,:,j], dm)
                end
            end
        end
    end

    return D
end

# discriminant measure for EMD
function discriminant_measure(Γ::AbstractArray{NamedTuple{(:coef, :weight), Tuple{S1,S2}},1},
                              dm::SignaturesDM) where
                             {S1<:Array{<:Number},
                              S2<:Union{AbstractFloat,Array{<:AbstractFloat}}}
    # Basic summary of data
    C = length(Γ)
    @assert C > 1
    n = size(Γ[1][:coef], 1)
    L = size(Γ[1][:coef], 2)

    D = zeros(Float64, (n,L))
    @inbounds begin
        for i in 1:(C-1)
            for j in (i+1):C
                D += discriminant_measure(Γ[i], Γ[j], dm)
            end
        end
    end

    return D
end

# discriminant measure between 2 energy maps
function discriminant_measure(Γ₁::AbstractArray{T}, Γ₂::AbstractArray{T},
        dm::ProbabilityDensityDM) where T<:Number

    # parameter checking and basic summary
    @assert 2 <= ndims(Γ₁) <= 3
    @assert size(Γ₁) == size(Γ₂)
    n = size(Γ₁,1)
    L = size(Γ₁,2)

    D = Array{T, 2}(undef, (n,L))
    @inbounds begin
        for i in axes(D,1)
            for j in axes(D,2)
                if ndims(Γ₁) == 2       # time frequency energy map case
                    D[i,j] = discriminant_measure(Γ₁[i,j], Γ₂[i,j], dm)
                else                    # probability density energy map case
                    for k in axes(Γ₁,3)
                        D[i,j] += discriminant_measure(Γ₁[i,j,k], Γ₂[i,j,k], dm)
                    end
                end
            end
        end
    end

    return D
end

# discriminant measure between 2 nergy maps for EMD
function discriminant_measure(Γ₁::NamedTuple{(:coef, :weight), Tuple{S1,S2}},
        Γ₂::NamedTuple{(:coef, :weight), Tuple{S1,S2}},
        dm::SignaturesDM) where
        {S1<:Array{T} where T<:Number,
        S2<:Union{AbstractFloat,Array{<:AbstractFloat}}}

    # parameter checking and basic summary
    n = size(Γ₁[:coef],1)
    L = size(Γ₁[:coef],2)
    @assert n == size(Γ₁[:coef],1) == size(Γ₂[:coef],1)
    @assert L == size(Γ₁[:coef],2) == size(Γ₂[:coef],2)

    D = Array{Float64,2}(undef, (n,L))
    for i in 1:n
        for j in 1:L
            # signatures
            if typeof(Γ₁[:weight]) <: AbstractFloat # equal weight
                P = (coef=Γ₁[:coef][i,j,:], weight=Γ₁[:weight])
                Q = (coef=Γ₂[:coef][i,j,:], weight=Γ₂[:weight])
            else                                    # probability density weight
                P = (coef=Γ₁[:coef][i,j,:], weight=Γ₁[:weight][i,j,:])
                Q = (coef=Γ₂[:coef][i,j,:], weight=Γ₂[:weight][i,j,:])
            end
            D[i,j] = discriminant_measure(P, Q, dm)
        end
    end
    return D
end

# Asymmetric Relative Entropy
function discriminant_measure(p::T, q::T, dm::AsymmetricRelativeEntropy) where
        T<:Number

    # parameter checking
    @assert p >= 0 && q >= 0

    if p == 0 || q == 0
        return 0
    else
        return p * log(p/q)
    end
end

# Symmetric Relative Entropy
function discriminant_measure(p::T, q::T, dm::SymmetricRelativeEntropy) where
        T<:Number

    return discriminant_measure(p, q, AsymmetricRelativeEntropy()) +
        discriminant_measure(q, p, AsymmetricRelativeEntropy())
end

# Hellinger Distance
function discriminant_measure(p::T, q::T, dm::HellingerDistance) where T<:Number
    return (sqrt(p) - sqrt(q))^2
end

# Lᵖ Distance
function discriminant_measure(p::T, q::T, dm::LpDistance) where T<:Number
    return (p - q)^dm.p
end

# Earth Mover Distance
function discriminant_measure(P::NamedTuple{(:coef, :weight), Tuple{S1,S2}},
        Q::NamedTuple{(:coef, :weight), Tuple{S1, S2}},
        dm::EarthMoverDistance) where
        {S1<:Vector{T} where T<:Number,
        S2<:Union{AbstractFloat,Vector{<:AbstractFloat}}}

    # assigning tuple signatures into coef and weight
    p, w_p = P
    q, w_q = Q

    # sort signature values
    p_order = sortperm(p)
    p = p[p_order]
    w_p = typeof(w_p)<:AbstractFloat ? repeat([w_p], length(p)) : w_p[p_order]
    q_order = sortperm(q)
    q = q[q_order]
    w_q = typeof(w_q)<:AbstractFloat ? repeat([w_q], length(q)) : w_q[q_order]

    # concatenate p and q, then sort them
    r = [p; q]
    sort!(r)

    # compute emd
    n = length(r)
    emd = 0
    for i in 1:(n-1)
        # get total weight of p and q less than or equal to r[i]
        p_less = p .<= r[i]
        ∑w_p = sum(w_p[p_less])
        q_less = q .<= r[i]
        ∑w_q = sum(w_q[q_less])
        # add into emd
        emd += abs(∑w_p - ∑w_q) * (r[i+1] - r[i])
    end
    emd /= (sum(w_p) + sum(w_q))
    return emd
end

## DISCRIMINATION POWER
"""
Discriminant Power measure for the Local Discriminant Basis. Current available
measures are
- [`BasisDiscriminantMeasure`](@ref)
- [`FishersClassSeparability`](@ref)
- [`RobustFishersClassSeparability`](@ref)
"""
abstract type DiscriminantPower end

"""
    BasisDiscriminantMeasure <: DiscriminantPower

This is the discriminant measure of a single basis function computed in a
previous step to construct the energy maps.
"""
struct BasisDiscriminantMeasure <: DiscriminantPower end

@doc raw"""
    FishersClassSeparability <: DiscriminantPower

The Fisher's class separability of the expansion coefficients in the basis
function.

Equation: ``\frac{\sum_{c=1}^C \pi_c({\rm mean}_i(\alpha_{\lambda,i}^{(c)}) - {\rm mean}_c({\rm mean}_i(\alpha_{\lambda,i}^{(c)})))^2}{\sum_{c=1}^C \pi_c {\rm var}_i(\alpha_{\lambda,i}^{(c)})}``
"""
struct FishersClassSeparability <: DiscriminantPower end

@doc raw"""
    RobustFishersClassSeparability <: DiscriminantPower

The robust version of Fisher's class separability of the expansion coefficients
in the basis function.

Equation: ``\frac{\sum_{c=1}^C \pi_c({\rm med}_i(\alpha_{\lambda,i}^{(c)}) - {\rm med}_c({\rm med}_i(\alpha_{\lambda,i}^{(c)})))^2}{\sum_{c=1}^C \pi_c {\rm mad}_i(\alpha_{\lambda,i}^{(c)})}``
"""
struct RobustFishersClassSeparability <: DiscriminantPower end

"""
    discriminant_power(D, tree, dp)

Returns the discriminant power of each leaf from the local discriminant basis
(LDB) tree.
"""
function discriminant_power(D::AbstractArray{T,2}, tree::BitVector,
        dp::BasisDiscriminantMeasure) where T<:Number

    @assert length(tree) == size(D,1) - 1

    power = getbasiscoef(D, tree)
    order = sortperm(power, rev = true)

    return (power, order)
end

function discriminant_power(coefs::AbstractArray{T,2},
                            y::AbstractVector{S},
                            dp::FishersClassSeparability) where {T<:Number, S}

    n = size(coefs,1)                             # signal length

    classes = unique(y)
    C = length(coefs)                             # number of classes

    Nᵢ = Array{T,1}(undef, C)
    Eαᵢ = Array{T,2}(undef, (n,C))                # mean of each entry
    Varαᵢ = Array{T,2}(undef, (n,C))              # variance of each entry
    @inbounds begin
        for (i, c) in enumerate(classes)
            idx = findall(yᵢ -> yᵢ == c, y)
            Nᵢ[i] = length(idx)
            Eαᵢ[:,i] = mean(coefs[:, idx], dims = 2)
            Varαᵢ[:,i] = var(coefs[:, idx], dims = 2)
        end
    end
    Eα = mean(Eαᵢ, dims = 2)                      # overall mean of each entry
    pᵢ = Nᵢ / sum(Nᵢ)                             # proportions of each class

    power = ((Eαᵢ - (Eα .* Eαᵢ)).^2 * pᵢ) ./ (Varαᵢ * pᵢ)
    order = sortperm(power, rev = true)

    return (power, order)
end

function discriminant_power(coefs::AbstractArray{T,2},
                            y::AbstractVector{S},
                            dp::RobustFishersClassSeparability) where {T<:Number,S}

    n = size(coefs,1)                            # signal length

    classes = unique(y)
    C = length(classes)

    Nᵢ = Array{T,1}(undef, C)
    Medαᵢ = Array{T,2}(undef, (n,C))             # mean of each entry
    Madαᵢ = Array{T,2}(undef, (n,C))             # variance of each entry
    @inbounds begin
        for (i, c) in enumerate(classes)
            idx = findall(yᵢ -> yᵢ == c, y)
            Nᵢ[i] = length(idx)
            Medαᵢ[:,i] = median(coefs[:, idx], dims = 2)
            Madαᵢ[:,i] = mapslices(x -> mad(x, normalize = false), coefs[:, idx],
                dims = 2)
        end
    end
    Medα = median(Medαᵢ, dims = 2)               # overall mean of each entry
    pᵢ = Nᵢ / sum(Nᵢ)                            # proportions of each class

    power = ((Medαᵢ - (Medα .* Medαᵢ)).^2 * pᵢ) ./ (Madαᵢ * pᵢ)
    order = sortperm(power, rev = true)

    return (power, order)
end

## LOCAL DISCRIMINANT BASIS
"""
    LocalDiscriminantBasis

Class type for the Local Discriminant Basis (LDB), a feature selection algorithm
developed by N. Saito and R. Coifman in "Local Discriminant Bases and Their
Applications" in the Journal of Mathematical Imaging and Vision, Vol 5, 337-358
(1995). This struct contains the following field values:

# Parameters and Attributes:
- `wt::DiscreteWavelet`: a discrete wavelet for transform purposes
- `max_dec_level::Union{Integer, Nothing}`: max level of wavelet packet
    decomposition to be computed.
- `dm::DiscriminantMeasure`: the discriminant measure for the LDB algorithm.
    Supported measures are the `AsymmetricRelativeEntropy()`, `LpDistance()`,
    `SymmetricRelativeEntropy()`, and `HellingerDistance()`
- `en::EnergyMap`: the type of energy map used. Supported maps are
    `TimeFrequency()`, `ProbabilityDensity()`, and `Signatures()`.
- `dp::DiscriminantPower()`: the measure of discriminant power among expansion
    coefficients. Supported measures are `BasisDiscriminantMeasure()`,
    `FishersClassSeparability()`, and `RobustFishersClassSeparability()`.
- `top_k::Union{Integer, Nothing}`: the top-k coefficients used in each node to
    determine the discriminant measure.
- `n_features::Union{Integer, Nothing}`: the dimension of output after
    undergoing feature selection and transformation.
- `n::Union{Integer, Nothing}`: length of signal
- `Γ::Union{AbstractArray{<:AbstractFloat},
AbstractArray{NamedTuple{(:coef, :weight), Tuple{S1, S2}}} where {S1<:Array{T}
    where T<:AbstractFloat, S2<:Union{AbstractFloat, Array{<:AbstractFloat}}},
    Nothing}`: computed energy map
- `DM::Union{AbstractArray{<:AbstractFloat}, Nothing}`: computed discriminant
    measure
- `cost::Union{AbstractVector{<:AbstractFloat}, Nothing}`: computed wavelet
    packet decomposition (WPD) tree cost based on the discriminant measure `DM`.
- `tree::Union{BitVector, Nothing}`: computed best WPD tree based on the
    discriminant measure `DM`.
- `DP::Union{AbstractVector{<:AbstractFloat}, Nothing}`: computed discriminant
    power
- `order::Union{AbstractVector{Integer}, Nothing}`: ordering of `DP` by
    descending order.
"""
mutable struct LocalDiscriminantBasis
    # to be declared by user
    wt::DiscreteWavelet
    max_dec_level::Union{Integer, Nothing}
    dm::DiscriminantMeasure
    en::EnergyMap
    dp::DiscriminantPower
    top_k::Union{Integer, Nothing}
    n_features::Union{Integer, Nothing}
    # to be computed in fit!
    n::Union{Integer, Nothing}
    Γ::Union{AbstractArray{<:AbstractFloat},
             AbstractArray{NamedTuple{(:coef, :weight), Tuple{S1, S2}}} where
                {S1<:Array{T} where T<:AbstractFloat,
                 S2<:Union{AbstractFloat, Array{<:AbstractFloat}}},
             Nothing}
    DM::Union{AbstractArray{<:AbstractFloat}, Nothing}
    cost::Union{AbstractVector{<:AbstractFloat}, Nothing}
    tree::Union{BitVector, Nothing}
    DP::Union{AbstractVector{<:AbstractFloat}, Nothing}
    order::Union{AbstractVector{Integer}, Nothing}
end

"""
    LocalDiscriminantBasis([;
        wt=wavelet(WT.haar),
        max_dec_level=nothing,
        dm=AsymmetricRelativeEntropy(), em=TimeFrequency(),
        dp=BasisDiscriminantMeasure(), top_k=nothing,
        n_features=nothing]
    )

Class constructor for `LocalDiscriminantBasis`.

# Arguments:
- `wt::DiscreteWavelet`: Wavelet used for decomposition of signals. Default is
    set to be `wavelet(WT.haar)`.
- `max_dec_level::Union{Integer, Nothing}`: max level of wavelet packet
    decomposition to be computed. When `max_dec_level=nothing`, the maximum
    transform levels will be used. Default is set to be `nothing`.
- `dm::DiscriminantMeasure`: the discriminant measure for the LDB algorithm.
    Supported measures are the `AsymmetricRelativeEntropy()`, `LpDistance()`,
    `SymmetricRelativeEntropy()`, and `HellingerDistance()`. Default is set to
    be `AsymmetricRelativeEntropy()`.
- `en::EnergyMap`: the type of energy map used. Supported maps are
    `TimeFrequency()` and `ProbabilityDensity()`. Default is set to be
    `TimeFrequency()`.
- `dp::DiscriminantPower=BasisDiscriminantMeasure()`: the measure of
    discriminant power among expansion coefficients. Supported measures are
    `BasisDiscriminantMeasure()`, `FishersClassSeparability()`, and
    `RobustFishersClassSeparability()`. Default is set to be `BasisDiscriminantMeasure()`.
- `top_k::Union{Integer, Nothing}`: the top-k coefficients used in each node to
    determine the discriminant measure. When `top_k=nothing`, all coefficients
    are used to determine the discriminant measure. Default is set to be
    `nothing`.
- `n_features::Union{Integer, Nothing}`: the dimension of output after
    undergoing feature selection and transformation. When `n_features=nothing`,
    all features will be returned as output. Default is set to be `nothing`.
"""
function LocalDiscriminantBasis(; wt::DiscreteWavelet=wavelet(WT.haar),
        max_dec_level::Union{Integer, Nothing}=nothing,
        dm::DiscriminantMeasure=AsymmetricRelativeEntropy(),
        en::EnergyMap=TimeFrequency(),
        dp::DiscriminantPower=BasisDiscriminantMeasure(),
        top_k::Union{Integer, Nothing}=nothing,
        n_features::Union{Integer, Nothing}=nothing)

    return LocalDiscriminantBasis(
        wt, max_dec_level, dm, en, dp, top_k, n_features,
        nothing, nothing, nothing, nothing, nothing, nothing, nothing
    )
end

"""
    fit!(f, X, y)

Fits the Local Discriminant Basis feature selection algorithm `f` onto the
signals `X` (or the decomposed signals `Xw`) with labels `y`.

**See also:** [`LocalDiscriminantBasis`](@ref), [`fit_transform`](@ref),
    [`transform`](@ref), [`inverse_transform`](@ref), [`change_nfeatures`](@ref)
"""
function fit!(f::LocalDiscriminantBasis, X::AbstractArray{S,2},
        y::AbstractVector{T}) where {S<:Number, T}

    # basic summary of data
    n, N = size(X)

    # change LocalDiscriminantBasis parameters if necessary
    f.max_dec_level = f.max_dec_level === nothing ?
        maxtransformlevels(n) : f.max_dec_level
    @assert 1 <= f.max_dec_level <= maxtransformlevels(n)

    # wavelet packet decomposition
    Xw = Array{S, 3}(undef, (n,f.max_dec_level+1,N))
    @inbounds begin
        for i in axes(Xw,3)
            Xw[:,:,i] = wpd(X[:,i], f.wt, f.max_dec_level)
        end
    end

    # fit local discriminant basis
    fit!(f, Xw, y)
    return nothing
end

function fit!(f::LocalDiscriminantBasis, Xw::AbstractArray{S,3},
        y::AbstractVector{T}) where {S<:Number, T}

    # basic summary of data
    c = unique(y)       # unique classes
    nc = length(c)      # number of classes
    Ny = length(y)
    f.n, L, Nx = size(Xw)

    # change LocalDiscriminantBasis parameters if necessary
    f.top_k = f.top_k === nothing ? f.n : f.top_k
    f.n_features = f.n_features === nothing ? f.n : f.n_features
    f.max_dec_level = f.max_dec_level === nothing ? L-1 : f.max_dec_level

    # parameter checking
    @assert Nx == Ny
    @assert 1 <= f.top_k <= f.n
    @assert 1 <= f.n_features <= f.n
    @assert f.max_dec_level+1 == L
    @assert 1 <= f.max_dec_level <= maxtransformlevels(f.n)
    @assert nc > 1
    @assert isdyadic(f.n)

    # construct energy map for each class
    f.Γ = energy_map(Xw, y, f.en)

    # compute discriminant measure D and obtain tree cost
    f.DM = discriminant_measure(f.Γ, f.dm)
    f.cost = Vector{Float64}(undef, 1<<L-1)
    @inbounds begin
        for i in eachindex(f.cost)
            lθ = floor(Integer, log2(i))    # node level
            θ = i - 1<<lθ                   # node
            nθ = nodelength(f.n, lθ)        # node length
            rng = (θ*nθ+1):((θ+1)*nθ)
            if f.top_k < nθ
                DMθ = f.DM[rng,lθ+1]        # discriminant measure for node θ
                sort!(DMθ, rev=true)
                f.cost[i] = sum(DMθ[1:f.top_k])
            else
                f.cost[i] = sum(f.DM[rng,lθ+1])
            end
        end
    end

    # select best tree and best set of expansion coefficients
    f.tree = bestbasis_treeselection(f.cost, f.n, :max)
    Xc = getbasiscoefall(Xw, f.tree)

    # obtain and order basis functions by power of discrimination
    (f.DP, f.order) = f.dp == BasisDiscriminantMeasure() ?
        discriminant_power(f.DM, f.tree, f.dp) :
        discriminant_power(Xc, y, f.dp)

    return nothing
end

"""
    transform(f, X)

Extract the LDB features on signals `X`.

**See also:** [`LocalDiscriminantBasis`](@ref), [`fit!`](@ref),
    [`fit_transform`](@ref), [`inverse_transform`](@ref), [`change_nfeatures`](@ref)
"""
function transform(f::LocalDiscriminantBasis, X::AbstractArray{T,2}) where T
    # check necessary measurements
    n, N = size(X)
    @assert f.max_dec_level !== nothing
    @assert f.top_k !== nothing
    @assert f.n_features !== nothing
    @assert f.n !== nothing
    @assert f.Γ !== nothing
    @assert f.DM !== nothing
    @assert f.cost !== nothing
    @assert f.tree !== nothing
    @assert f.DP !== nothing
    @assert f.order !== nothing
    @assert n == f.n
    @assert 1 <= f.max_dec_level <= maxtransformlevels(f.n)
    @assert 1 <= f.top_k <= f.n
    @assert 1 <= f.n_features <= f.n

    # wpt on X based on given f.tree
    Xc = Array{T, 2}(undef, (n,N))
    @inbounds begin
        for i in axes(Xc,2)
            Xc[:,i] = wpt(X[:,i], f.wt, f.tree)
        end
    end
    return Xc[f.order[1:f.n_features],:]
end

"""
    fit_transform(f, X, y)

Fit and transform the signals `X` with labels `y` based on the LDB class `f`.

**See also:** [`LocalDiscriminantBasis`](@ref), [`fit!`](@ref),
    [`transform`](@ref), [`inverse_transform`](@ref), [`change_nfeatures`](@ref)
"""
function fit_transform(f::LocalDiscriminantBasis, X::AbstractArray{S,2},
        y::AbstractVector{T}) where {S<:Number, T}

    # get necessary measurements
    n, N = size(X)
    f.max_dec_level = f.max_dec_level===nothing ?
        maxtransformlevels(n) : f.max_dec_level
    @assert 1 <= f.max_dec_level <= maxtransformlevels(n)

    # wpd on X
    Xw = Array{S, 3}(undef, (n,f.max_dec_level+1,N))
    @inbounds begin
        for i in axes(Xw,3)
            Xw[:,:,i] = wpd(X[:,i], f.wt, f.max_dec_level)
        end
    end

    # fit LDB and return best features
    fit!(f, Xw, y)
    Xc = getbasiscoefall(Xw, f.tree)
    return Xc[f.order[1:f.n_features],:]
end

"""
    inverse_transform(f, x)

Compute the inverse transform on the feature matrix `x` to form the original
signal based on the LDB class `f`.

**See also:** [`LocalDiscriminantBasis`](@ref), [`fit!`](@ref),
    [`fit_transform`](@ref), [`transform`](@ref), [`change_nfeatures`](@ref)
"""
function inverse_transform(f::LocalDiscriminantBasis,
        x::AbstractArray{T,2}) where T<:Number

    # get necessary measurements
    @assert size(x,1) == f.n_features
    N = size(x,2)

    # insert the features x into the full matrix X padded with 0's
    Xc = zeros(f.n, N)
    Xc[f.order[1:f.n_features],:] = x

    # iwpt on X
    X = Array{T,2}(undef, (f.n, N))
    @inbounds begin
        for i in axes(Xc, 2)
            X[:,i] = iwpt(Xc[:,i], f.wt, f.tree)
        end
    end
    return X
end

"""
    change_nfeatures(f, x, n_features)

Change the number of features from `f.n_features` to `n_features`.

Note that if the input `n_features` is larger than `f.n_features`, it results in
the regeneration of signals based on the current `f.n_features` before
reselecting the features. This will cause additional features to be less
accurate and effective.

**See also:** [`LocalDiscriminantBasis`](@ref), [`fit!`](@ref), [`fit_transform`](@ref),
    [`transform`](@ref), [`inverse_transform`](@ref)
"""
function change_nfeatures(f::LocalDiscriminantBasis, x::AbstractArray{T,2},
        n_features::Integer) where T<:Number

    # check measurements
    @assert f.n_features !== nothing
    @assert size(x,1) == f.n_features ||
        throw(ArgumentError("f.n_features and number of rows of x do not match!"))
    @assert 1 <= n_features <= f.n

    # change number of features
    if f.n_features >= n_features
        f.n_features = n_features
        y = x[1:f.n_features,:]
    else
        @warn "Proposed n_features larger than currently saved n_features. Results will be less accurate since inverse_transform and transform is involved."
        X = inverse_transform(f, x)
        f.n_features = n_features
        y = transform(f, X)
    end

    return y
end

end # end module


text:
# CUDA device pointer type

import Base: eltype, convert, unsafe_convert, cconvert, isnull

export
    DevicePtr

# This wrapper type contains a pointer, but prevents all conversions from and to
# regular pointers. This ensures we don't mix host and device pointers.

immutable DevicePtr{T}
    inner::Ptr{T}

    DevicePtr() = new(C_NULL)
    DevicePtr{T}(::Ptr{T}) = throw(InexactError())
    DevicePtr{T}(::Bool, x::Ptr{T}) = new(x)
    # NOTE: this hidden constructor is not meant to be used -- use unsafe_convert instead
end

# Simple conversions between Ptr and DevicePtr are disallowed
convert{T}(::Type{Ptr{T}}, p::DevicePtr{T}) = throw(InexactError())
convert{T}(::Type{DevicePtr{T}}, p::Ptr{T}) = throw(InexactError())

# Unsafe ones are tolerated
unsafe_convert{T}(::Type{Ptr{T}}, p::DevicePtr{T}) = p.inner
unsafe_convert{T}(::Type{DevicePtr{T}}, p::Ptr{T}) = DevicePtr{T}(true, p)

cconvert{P<:DevicePtr}(::Type{P}, x) = x # defer conversions to DevicePtr to unsafe_convert

# Some convenience methods (which are already defined for Ptr{T},
# but due to the disallowed conversions we can't use those)
isnull{T}(p::DevicePtr{T}) = (p.inner == C_NULL)
eltype{T}(::Type{DevicePtr{T}}) = T

ext: jl


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

module FMI

include("FMI1.jl")
include("FMI2_comp.jl")
include("FMI3.jl")
include("assertions.jl")

# ToDo: Submodules
include("FMI2_sim.jl")
#include("FMI_neural.jl")
include("FMI_plot.jl")

### EXPORTING LISTS START ###

# FMI.jl
export fmiLoad, fmiSimulate, fmiSimulateCS, fmiSimulateME, fmiUnload
export fmiGetNumberOfStates, fmiGetTypesPlatform, fmiGetVersion, fmiInstantiate!, fmiFreeInstance!
export fmiSetDebugLogging, fmiSetupExperiment, fmiEnterInitializationMode, fmiExitInitializationMode, fmiTerminate , fmiReset
export fmiGetReal, fmiSetReal, fmiGetInteger, fmiSetInteger, fmiGetBoolean, fmiSetBoolean, fmiGetString, fmiSetString, fmiGetReal!, fmiGetInteger!, fmiGetBoolean!, fmiGetString!
export fmiGetFMUstate, fmiSetFMUstate, fmiFreeFMUstate, fmiSerializedFMUstateSize, fmiSerializeFMUstate, fmiDeSerializeFMUstate
export fmiGetDirectionalDerivative, fmiDoStep, fmiSetTime, fmiSetContinuousStates, fmi2EnterEventMode, fmiNewDiscreteStates
export fmiEnterContinuousTimeMode, fmiCompletedIntegratorStep, fmiGetDerivatives, fmiGetEventIndicators, fmiGetContinuousStates, fmiGetNominalsOfContinuousStates

# FMI2.jl
export FMU2, fmi2True, fmi2False
export fmi2SimulationResult, fmi2SimulationResultGetValuesAtIndex, fmi2SimulationResultGetTime, fmi2SimulationResultGetValues
export fmi2String2ValueReference, fmi2ValueReference2String
export fmi2Unzip, fmi2Load, fmi2Unload
export fmi2Simulate, fmi2SimulateCS, fmi2SimulateME
export fmi2GetTypesPlatform, fmi2GetVersion
export fmi2Instantiate!, fmi2FreeInstance!, fmi2SetDebugLogging
export fmi2SetupExperiment, fmi2EnterInitializationMode, fmi2ExitInitializationMode, fmi2Terminate, fmi2Reset
export fmi2GetReal, fmi2SetReal, fmi2GetInteger, fmi2SetInteger, fmi2GetBoolean, fmi2SetBoolean, fmi2GetString, fmi2SetString
export fmi2GetFMUstate, fmi2SetFMUstate, fmi2FreeFMUstate, fmi2SerializedFMUstateSize, fmi2SerializeFMUstate, fmi2DeSerializeFMUstate
export fmi2GetDirectionalDerivative
export fmi2SetRealInputDerivatives, fmi2GetRealOutputDerivatives
export fmi2DoStep, fmi2CancelStep
export fmi2GetStatus, fmi2GetRealStatus, fmi2GetIntegerStatus, fmi2GetBooleanStatus, fmi2GetStringStatus
export fmi2SetTime, fmi2SetContinuousStates
export fmi2EnterEventMode, fmi2NewDiscreteStates, fmi2EnterContinuousTimeMode, fmi2CompletedIntegratorStep, fmi2GetDerivatives, fmi2GetEventIndicators, fmi2GetContinuousStates, fmi2GetNominalsOfContinuousStates

# FMI2_comp.jl
export fmi2SetDebugLogging, fmi2SetupExperiment
export fmi2GetReal, fmi2SetReal, fmi2GetInteger, fmi2SetInteger, fmi2GetBoolean, fmi2SetBoolean, fmi2GetString, fmi2SetString
export fmi2GetFMUstate, fmi2FreeFMUstate, fmi2SerializedFMUstateSize, fmi2SerializeFMUstate, fmi2DeSerializeFMUstate
export fmi2GetDirectionalDerivative
export fmi2DoStep
export fmi2SetTime, fmi2SetContinuousStates
export fmi2NewDiscreteStates, fmi2CompletedIntegratorStep, fmi2GetDerivatives, fmi2GetEventIndicators, fmi2GetContinuousStates, fmi2GetNominalsOfContinuousStates

# FMI2_c.jl
export fmi2Instantiate, fmi2FreeInstance, fmi2SetDebugLogging
export fmi2GetTypesPlatform, fmi2GetVersion
export fmi2SetupExperiment, fmi2EnterInitializationMode, fmi2ExitInitializationMode, fmi2Terminate, fmi2Reset
export fmi2GetReal!, fmi2SetReal, fmi2GetInteger!, fmi2SetInteger, fmi2GetBoolean!, fmi2SetBoolean, fmi2GetString!, fmi2SetString
export fmi2GetFMUstate, fmi2SetFMUstate, fmi2FreeFMUstate, fmi2SerializedFMUstateSize, fmi2SerializeFMUstate, fmi2DeSerializeFMUstate
export fmi2GetDirectionalDerivative!
export fmi2SetRealInputDerivatives, fmi2GetRealOutputDerivatives
export fmi2DoStep, fmi2CancelStep
export fmi2GetStatus, fmi2GetRealStatus, fmi2GetIntegerStatus, fmi2GetBooleanStatus, fmi2GetStringStatus
export fmi2SetTime, fmi2SetContinuousStates
export fmi2EnterEventMode, fmi2NewDiscreteStates, fmi2EnterContinuousTimeMode, fmi2CompletedIntegratorStep!, fmi2GetDerivatives, fmi2GetEventIndicators, fmi2GetContinuousStates, fmi2GetNominalsOfContinuousStates

# FMI2_sim.jl
export fmi2SimulateME

# FMI_plot.jl
export fmiPlot

# FMI2_md.jl
# nothing to export

### EXPORTING LISTS END ###

fmi2Struct = Union{FMU2, fmi2Component}
fmi2Reference = Union{String, Array{String}, fmi2ValueReference, Array{fmi2ValueReference}}

""" Multiple Dispatch variants for FMUs with version 2.0.X """
function fmiLoad(pathToFMU::String)
    fmi2Load(pathToFMU)
end

"""Simulate an fmu according to its standard from 0.0 to t_stop"""
function fmiSimulate(fmu::fmi2Struct, dt::Real, t_start::Real = 0.0, t_stop::Real = 1.0, recordValues::Union{Array{fmi2ValueReference}, Array{String}} = [])
    fmi2Simulate(fmu, dt, t_start, t_stop, recordValues)
end

"""Simulate an CoSimulation fmu according to its standard from 0.0 to t_stop"""
function fmiSimulateCS(fmu::fmi2Struct, dt::Real, t_start::Real = 0.0, t_stop::Real = 1.0, recordValues::Union{Array{fmi2ValueReference}, Array{String}} = [])
    fmi2SimulateCS(fmu, dt, t_start, t_stop, recordValues)
end

"""Simulate an ModelExchange fmu according to its standard from 0.0 to t_stop"""
function fmiSimulateME(fmu::fmi2Struct, dt::Real, t_start::Real = 0.0, t_stop::Real = 1.0)
    fmi2SimulateME(fmu, dt, t_start, t_stop)
end

"""Unloads the FMU and all its instances and frees the allocated memory"""
function fmiUnload(fmu::FMU2)
    fmi2Unload(fmu)
end

"""Returns the number of states of the FMU"""
function fmiGetNumberOfStates(fmu::FMU2)
    length(fmu.modelDescription.stateValueReferences)
end

"""Returns the header file used to compile the FMU. By default returns "default" """
function fmiGetTypesPlatform(fmu::FMU2)
    fmi2GetTypesPlatform(fmu)
end

"""Returns the version of the FMU"""
function fmiGetVersion(fmu::FMU2)
    fmi2GetVersion(fmu)
end

"""Creates a new instance of the FMU"""
function fmiInstantiate!(fmu::FMU2; visible::Bool = false, loggingOn::Bool = false)

    version = fmiGetVersion(fmu)

    if version == "2.0"
        return fmi2Instantiate!(fmu; visible, loggingOn)
    else
        @assert false ["fmiInstantiate!(...): Unknwon FMI version $version !"]
    end

    nothing
end

"""Frees the allocated memory of the last instance of the FMU"""
function fmiFreeInstance!(fmu::FMU2)
    fmi2FreeInstance!(fmu)
end

"""Control the use of the logging callback function"""
function fmiSetDebugLogging(s::fmi2Struct)
    fmi2SetDebugLogging(s)
end

"""Initialize the Simulation boundries"""
function fmiSetupExperiment(fmu::fmi2Struct,
    toleranceDefined::Bool,
                tolerance::Real,
                startTime::Real,
                stopTimeDefined::Bool,
                stopTime::Real)

    fmi2SetupExperiment(fmu, toleranceDefined, tolerance, startTime, stopTimeDefined, stopTime)
end
"""Initialize the Simulation boundries"""
function fmiSetupExperiment(fmu::fmi2Struct, startTime::Real = 0.0, stopTime::Real = startTime; tolerance::Real = 0.0)

    toleranceDefined = (tolerance > 0.0)
    stopTimeDefined = (stopTime > startTime)

    fmi2SetupExperiment(fmu, toleranceDefined, tolerance, startTime, stopTimeDefined, stopTime)
end

"""Informs the FMU to enter initializaton mode"""
function fmiEnterInitializationMode(s::fmi2Struct)
    fmi2EnterInitializationMode(s)
end

"""Informs the FMU to exit initialization mode"""
function fmiExitInitializationMode(s::fmi2Struct)
    fmi2ExitInitializationMode(s)
end

"""Informs the FMU that the simulation run is terminated"""
function fmiTerminate(s::fmi2Struct)
    fmi2Terminate(s)
end

"""Resets the FMU after a simulation run"""
function fmiReset(s::fmi2Struct)
    fmi2Reset(s)
end

"""Returns the real values of an array of variables"""
function fmiGetReal(fmu::fmi2Struct, vr::fmi2Reference)
    fmi2GetReal(fmu, vr)
end

"""Writes the real values of an array of variables in the given field"""
function fmiGetReal!(fmu::fmi2Struct, vr::Union{Array{fmi2ValueReference}, Array{String}}, values::Array{<:Real})
    fmi2GetReal!(fmu, vr, values)
end

"""Set the values of an array of real variables"""
function fmiSetReal(fmu::fmi2Struct, vr::fmi2Reference, value::Array{<:Real})
    fmi2SetReal(fmu, vr, Array{Real}(value))
end

"""Set the value of a real variable"""
function fmiSetReal(fmu::fmi2Struct, vr::fmi2Reference, value::Real)
    fmi2SetReal(fmu, vr, Real(value))
end
"""Returns the integer values of an array of variables"""
function fmiGetInteger(fmu::fmi2Struct, vr::fmi2Reference)
    fmi2GetInteger(fmu, vr)
end

"""Writes the integer values of an array of variables in the given field"""
function fmiGetInteger!(fmu::fmi2Struct, vr::Union{Array{fmi2ValueReference}, Array{String}}, values::Array{<:Integer})
    fmi2GetInteger!(fmu, vr, values)
end

"""Set the values of an array of integer variables"""
function fmiSetInteger(fmu::fmi2Struct, vr::fmi2Reference, value::Array{<:Integer})
    fmi2SetInteger(fmu, vr, Array{Integer}(value))
end

"""Set the value of a integer variable"""
function fmiSetInteger(fmu::fmi2Struct, vr::fmi2Reference, value::Integer)
    fmi2SetInteger(fmu, vr, Integer(value))
end
"""Returns the boolean values of an array of variables"""
function fmiGetBoolean(fmu::fmi2Struct, vr::fmi2Reference)
    fmi2GetBoolean(fmu, vr)
end

"""Writes the boolean values of an array of variables in the given field"""
function fmiGetBoolean!(fmu::fmi2Struct, vr::Union{Array{fmi2ValueReference}, Array{String}}, values::Array{Bool})
    fmi2GetBoolean!(fmu, vr, values)
end

"""Set the values of an array of boolean variables"""
function fmiSetBoolean(fmu::fmi2Struct, vr::fmi2Reference, value::Array{Bool})
    fmi2SetBoolean(fmu, vr, value)
end

"""Set the value of a boolean variable"""
function fmiSetBoolean(fmu::fmi2Struct, vr::fmi2Reference, value::Bool)
    fmi2SetBoolean(fmu, vr, value)
end
"""Returns the string values of an array of variables"""
function fmiGetString(fmu::fmi2Struct, vr::fmi2Reference)
    fmi2GetString(fmu, vr)
end

"""Writes the string values of an array of variables in the given field"""
function fmiGetString!(fmu::fmi2Struct, vr::Union{Array{fmi2ValueReference}, Array{String}}, values::Array{String})
    fmi2GetString!(fmu, vr, values)
end

"""Set the values of an array of string variables"""
function fmiSetString(fmu::fmi2Struct, vr::fmi2Reference, value::Array{String})
    fmi2SetString(fmu, vr, value)
end

"""Set the value of a string variable"""
function fmiSetString(fmu::fmi2Struct, vr::fmi2Reference, value::String)
    fmi2SetString(fmu, vr, value)
end
"""Returns the FMU state of the fmu"""
function fmiGetFMUstate(fmu2::fmi2Struct)
    fmi2GetFMUstate(fmu2)
end
"""Sets the FMU to the given state"""
function fmiSetFMUstate(fmu2::fmi2Struct, state::fmi2FMUstate)
    fmi2SetFMUstate(fmu2, state)
end
"""Free the memory for the allocated FMU state"""
function fmiFreeFMUstate(fmu2::fmi2Struct, state::fmi2FMUstate)
    fmi2FreeFMUstate(fmu2, state)
end
"""Returns the size of the byte vector the FMU can be stored in"""
function fmiSerializedFMUstateSize(c::fmi2Struct, state::fmi2FMUstate)
    fmi2SerializedFMUstateSize(c, state)
end
"""Serialize the data in the FMU state pointer"""
function fmiSerializeFMUstate(c::fmi2Struct, state::fmi2FMUstate)
    fmi2SerializeFMUstate(c, state)
end
"""Deserialize the data in the FMU state pointer"""
function fmiDeSerializeFMUstate(c::fmi2Struct, serializedState::Array{fmi2Byte})
    fmi2DeSerializeFMUstate(c, serializedState)
end
"""Returns the values of the directional derivatives"""
function fmiGetDirectionalDerivative(fmu::fmi2Struct,
                                     vUnknown_ref::Array{Cint},
                                     vKnown_ref::Array{Cint},
                                     dvKnown::Array{Real} = nothing,
                                     dvUnknown::Array{Real} = nothing)
    fmi2GetDirectionalDerivative(fmu,
                                 Array{fmi2ValueReference}(vUnknown_ref),
                                 Array{fmi2ValueReference}(vKnown_ref),
                                 Array{Real}(dvKnown),
                                 Array{Real}(dvUnknown))
end

""" Wrapper for single directional derivative """
function fmiGetDirectionalDerivative(fmu::fmi2Struct, vUnknown_ref::Cint, vKnown_ref::Cint, dvKnown::Real = 1.0, dvUnknown::Real = 1.0)
    fmi2GetDirectionalDerivative(fmu, vUnknown_ref, vKnown_ref, dvKnown, dvUnknown)
end

"""Does one step in the CoSimulation FMU"""
function fmiDoStep(fmu::fmi2Struct, currentCommunicationPoint::Real, communicationStepSize::Real, noSetFMUStatePriorToCurrentPoint::Bool = true)
    fmi2DoStep(fmu, currentCommunicationPoint, communicationStepSize, noSetFMUStatePriorToCurrentPoint)
end
"""Does one step in the CoSimulation FMU"""
function fmiDoStep(c::fmi2Struct, communicationStepSize::Real)
    fmi2DoStep(c, communicationStepSize)
end
"""Set a time instant"""
function fmiSetTime(fmu2::fmi2Struct, time::Real)
    fmi2SetTime(fmu, time)
end
"""Set a new (continuous) state vector"""
function fmiSetContinuousStates(c::fmi2Struct, x::Union{Array{Float32}, Array{Float64}})
    fmi2SetContinuousStates(c, x)
end
"""The model enters Event Mode"""
function fmi2EnterEventMode(fmu2::fmi2Struct)
    fmi2EnterEventMode(fmu2)
end
"""Returns the next discrete states"""
function fmiNewDiscreteStates(c::fmi2Struct)
    fmi2NewDiscreteStates(c)
end
"""The model enters Continuous-Time Mode"""
function fmiEnterContinuousTimeMode(fmu2::fmi2Struct)
    fmi2EnterContinuousTimeMode(fmu2)
end
"""This function must be called by the environment after every completed step"""
function fmiCompletedIntegratorStep(fmu2::fmi2Struct,
                                     noSetFMUStatePriorToCurrentPoint::fmi2Boolean)
    fmi2CompletedIntegratorStep(fmu2,noSetFMUStatePriorToCurrentPoint)
end
"""Compute state derivatives at the current time instant and for the current states"""
function  fmiGetDerivatives(c::fmi2Struct)
    fmi2GetDerivatives(c)
end
"""Returns the event indicators of the FMU"""
function fmiGetEventIndicators(fmu2::fmi2Struct)
    fmi2GetEventIndicators(fmu2)
end
"""Return the new (continuous) state vector x"""
function fmiGetContinuousStates(fmu2::fmi2Struct)
    fmi2GetContinuousStates(fmu2)
end
"""Return the new (continuous) state vector x"""
function fmiGetNominalsOfContinuousStates(fmu2::fmi2Struct)
    fmi2GetNominalsOfContinuousStates(fmu2)
end

""" Multiple Dispatch fallback for FMUs with unsupported versions """
unsupportedFMUs = Union{FMU1,FMU3}
function fmiDoStep(fmu::unsupportedFMUs, currentCommunicationPoint::Real, communicationStepSize::Real, noSetFMUStatePriorToCurrentPoint::Bool = true)
    error(unsupportedFMU::errorType)
end

function fmiUnload(fmu::unsupportedFMUs)
    error(unsupportedFMU::errorType)
end

function fmiGetTypesPlatform(fmu::unsupportedFMUs)
    error(unsupportedFMU::errorType)
end

function fmiGetVersion(fmu::unsupportedFMUs)
    error(unsupportedFMU::errorType)
end

function fmiInstantiate!(fmu::unsupportedFMUs; visible::Bool = false, loggingOn::Bool = false)
    error(unsupportedFMU::errorType)
end

function fmiFreeInstance!(fmu::unsupportedFMUs)
    error(unsupportedFMU::errorType)
end

function fmiSetDebugLogging(fmu::unsupportedFMUs)
    error(unsupportedFMU::errorType)
end

function fmiSetupExperiment(fmu::unsupportedFMUs,
    toleranceDefined::Bool,
                tolerance::Real,
                startTime::Real,
                stopTimeDefined::Bool,
                stopTime::Real)
    error(unsupportedFMU::errorType)
end

function fmiEnterInitializationMode(fmu::unsupportedFMUs)
    error(unsupportedFMU::errorType)
end

function fmiExitInitializationMode(fmu::unsupportedFMUs)
    error(unsupportedFMU::errorType)
end

function fmiTerminate(fmu::unsupportedFMUs)
    error(unsupportedFMU::errorType)
end

function fmiReset(fmu::unsupportedFMUs)
    error(unsupportedFMU::errorType)
end

function fmiGetReal(fmu::unsupportedFMUs, vr::Array{fmi2Reference})
    error(unsupportedFMU::errorType)
end

function fmiGetReal(fmu::unsupportedFMUs, vr::fmi2Reference)
    error(unsupportedFMU::errorType)
end

function fmiSetReal(fmu::unsupportedFMUs, vr::Array{fmi2Reference}, value::Array{fmi2Real})
    error(unsupportedFMU::errorType)
end

function fmiSetReal(fmu::unsupportedFMUs, vr::fmi2ValueReference, value::fmi2Real)
    error(unsupportedFMU::errorType)
end

function fmiSimulate(fmu::unsupportedFMUs, dt::Real, t_start::Real = 0.0, t_stop::Real = 1.0)
    error(unsupportedFMU::errorType)
end

function fmiSimulateCS(fmu::unsupportedFMUs, dt::Real, t_start::Real = 0.0, t_stop::Real = 1.0)
    error(unsupportedFMU::errorType)
end

function fmiSimulateME(fmu::unsupportedFMUs, dt::Real, t_start::Real = 0.0, t_stop::Real = 1.0)
    error(unsupportedFMU::errorType)
end

function fmiSetReal(fmu::unsupportedFMUs, vr::fmi2Reference, value::Real)
    error(unsupportedFMU::errorType)
end

function fmiGetInteger(fmu::unsupportedFMUs, vr::fmi2Reference)
    error(unsupportedFMU::errorType)
end

function fmiSetInteger(fmu::unsupportedFMUs, vr::fmi2Reference, value::Array{Int64})
    error(unsupportedFMU::errorType)
end

function fmiSetInteger(fmu::unsupportedFMUs, vr::fmi2Reference, value::Int64)
    error(unsupportedFMU::errorType)
end

function fmiGetBoolean(fmu::unsupportedFMUs, vr::fmi2Reference)
    error(unsupportedFMU::errorType)
end

function fmiSetBoolean(fmu::unsupportedFMUs, vr::fmi2Reference, value::Array{Bool})
    error(unsupportedFMU::errorType)
end

function fmiSetBoolean(fmu::unsupportedFMUs, vr::fmi2Reference, value::Bool)
    error(unsupportedFMU::errorType)
end

function fmiGetString(fmu::unsupportedFMUs, vr::fmi2Reference)
    error(unsupportedFMU::errorType)
end

function fmiSetString(fmu::unsupportedFMUs, vr::fmi2Reference, value::Array{String})
    error(unsupportedFMU::errorType)
end

function fmiSetString(fmu::unsupportedFMUs, vr::fmi2Reference, value::String)
    error(unsupportedFMU::errorType)
end

function fmiGetFMUstate(fmu::unsupportedFMUs)
    error(unsupportedFMU::errorType)
end

function fmiSetFMUstate(fmu::unsupportedFMUs, state::fmi2FMUstate)
    error(unsupportedFMU::errorType)
end

function fmiFreeFMUstate(fmu::unsupportedFMUs)
    error(unsupportedFMU::errorType)
end

function fmiSerializedFMUstateSize(c::unsupportedFMUs, size::Int64)
    error(unsupportedFMU::errorType)
end

function fmiSerializeFMUstate(c::unsupportedFMUs, serializedState::fmi2Byte, size::Int64)
    error(unsupportedFMU::errorType)
end

function fmiDeSerializeFMUstate(c::unsupportedFMUs, serializedState::fmi2Byte, size::Int64)
    error(unsupportedFMU::errorType)
end

function fmiGetDirectionalDerivative(fmu::unsupportedFMUs,
                                     vUnknown_ref::Array{Cint},
                                     vKnown_ref::Array{Cint},
                                     dvKnown::Array{Real} = nothing,
                                     dvUnknown::Array{Real} = nothing)
    error(unsupportedFMU::errorType)
end

function fmiGetDirectionalDerivative(fmu::unsupportedFMUs, vUnknown_ref::Cint, vKnown_ref::Cint, dvKnown::Real = 1.0, dvUnknown::Real = 1.0)
    error(unsupportedFMU::errorType)
end

function fmiDoStep(c::unsupportedFMUs, communicationStepSize::Real)
    error(unsupportedFMU::errorType)
end

function fmiSetTime(fmu::unsupportedFMUs, time::Real)
    error(unsupportedFMU::errorType)
end

function fmiSetContinuousStates(c::unsupportedFMUs, x::Union{Array{Float32}, Array{Float64}})
    error(unsupportedFMU::errorType)
end

function fmi2EnterEventMode(fmu::unsupportedFMUs)
    error(unsupportedFMU::errorType)
end

function fmiNewDiscreteStates(c::unsupportedFMUs)
    error(unsupportedFMU::errorType)
end

function fmiEnterContinuousTimeMode(fmu::unsupportedFMUs)
    error(unsupportedFMU::errorType)
end

function fmiCompletedIntegratorStep(fmu::unsupportedFMUs,
                                     noSetFMUStatePriorToCurrentPoint::fmi2Boolean)
    error(unsupportedFMU::errorType)
end

function  fmiGetDerivatives(c::unsupportedFMUs)
    error(unsupportedFMU::errorType)
end

function fmiGetEventIndicators(fmu::unsupportedFMUs)
    error(unsupportedFMU::errorType)
end

function fmiGetContinuousStates(fmu::unsupportedFMUs)
    error(unsupportedFMU::errorType)
end

function fmiGetNominalsOfContinuousStates(fmu::unsupportedFMUs)
    error(unsupportedFMU::errorType)
end

end # module FMI


text:
using GPUAcceleratedTracking, DrWatson, DataFrames
@quickactivate "GPUAcceleratedTracking"

raw_data_df = collect_results(datadir("benchmarks/kernel/atbenchmark"))

using Query, CairoMakie

gnss_string = Dict(
    "GPSL1" => "GPS L1 C/A",
    "GPSL5" => "GPS L5 P"
)
os = "linux"

fig = Figure(
    resolution = (1000, 1000),
    font = "Times New Roman",
    # fontsize = 10,
)

sort!(raw_data_df)

samples = unique(Vector{Int64}(raw_data_df[!, :num_samples]))
sort!(samples)
freqs_gpsl1 = samples / 0.001
freqs_gpsl5 = samples[5:end] / 0.001
algorithm_names = unique(Vector{String}(raw_data_df[!, :algorithm]))
sort!(algorithm_names)
# algorithm_names = algorithm_names[2:end] #disregard 1 alg

# PLOT GPS L1 C/A
gnss = "GPSL1"
num_ants = 1; num_correlators = 3

gpsl1_1ms_1_3_ax = Axis(fig,
    xlabel = "Sampling Frequency [Hz]",
    ylabel = "Processing Time [s]",
    xscale = log10,
    yscale = log10,
    xmajorgridvisible = true,
    xminorgridvisible = true,
    xminorticksvisible = true,
    xminorticks = IntervalsBetween(9),
    xticks = collect(10 .^ (6:9)),
    # yticks = collect(10.0 .^ (-6:-2)),
    title = "$(gnss_string[gnss]) T=1 ms, M=$num_ants, L=$num_correlators"
)
# xlims!(gpsl1_1ms_1_3_ax, [10^(6),  3 * 10^(8)])
# ylims!(gpsl1_1ms_1_3_ax, [10^(-6), 2*10^(-3)])
markers = [:rect, :utriangle, :circle, :diamond, :dtriangle]
lin = Array{Lines}(undef, length(algorithm_names));
sca = Array{Scatter}(undef, length(algorithm_names));
for (idx, algorithm) in enumerate(algorithm_names)
    times = Float64.(
        (
            raw_data_df |>
                @filter(
                    _.processor == "GPU" &&
                    _.algorithm == algorithm_names[idx] &&
                    _.num_ants == num_ants &&
                    _.num_correlators == num_correlators &&
                    _.GNSS == gnss
                ) |>
                @map(
                    {
                        _.Minimum
                    }
                ) |> DataFrame
        ).Minimum
    )
    lin[idx] = lines!(
        gpsl1_1ms_1_3_ax,
        freqs_gpsl1,
        times .* 10 ^ (-9)
    )
    sca[idx] = scatter!(
        gpsl1_1ms_1_3_ax,
        freqs_gpsl1,
        times .* 10 ^ (-9),
        marker = markers[idx],
        markersize = 15
    )
end
cputimes = Float64.(
        (
            raw_data_df |>
                @filter(
                    _.processor == "CPU" &&
                    _.num_ants == num_ants &&
                    _.num_correlators == num_correlators &&
                    _.GNSS == gnss
                ) |>
                @map(
                    {
                        _.Minimum
                    }
                ) |> DataFrame
        ).Minimum
    )
cpulin = lines!(
    gpsl1_1ms_1_3_ax,
    freqs_gpsl1,
    cputimes .* 10 ^ (-9),
)
cpusca = scatter!(
    gpsl1_1ms_1_3_ax,
    freqs_gpsl1,
    cputimes .* 10 ^ (-9),
    marker = :star5,
    markersize = 15
)
realtime = lines!(gpsl1_1ms_1_3_ax, [10^(6), 3 * 10^(8)], [10^(-3), 10^(-3)], linestyle = :dashdot, color =:dimgrey )
fig[1,1] = gpsl1_1ms_1_3_ax
fig
# PLOT GPS L1 C/A
gnss = "GPSL1"
num_ants = 4; num_correlators = 3

gpsl1_1ms_4_3_ax = Axis(fig,
    xlabel = "Sampling Frequency [Hz]",
    ylabel = "Processing Time [s]",
    xscale = log10,
    yscale = log10,
    xmajorgridvisible = true,
    xminorgridvisible = true,
    xminorticksvisible = true,
    xminorticks = IntervalsBetween(9),
    xticks = collect(10 .^ (6:9)),
    # yticks = collect(10.0 .^ (-6:-2)),
    title = "$(gnss_string[gnss]) T=1 ms, M=$num_ants, L=$num_correlators"
)
# xlims!(gpsl1_1ms_4_3_ax, [10^(6),  3 * 10^(8)])
# ylims!(gpsl1_1ms_4_3_ax, [10^(-6), 2*10^(-3)])
markers = [:rect, :utriangle, :circle, :diamond, :dtriangle]
lin = Array{Lines}(undef, length(algorithm_names));
sca = Array{Scatter}(undef, length(algorithm_names));
for (idx, algorithm) in enumerate(algorithm_names)
    times = Float64.(
        (
            raw_data_df |>
                @filter(
                    _.processor == "GPU" &&
                    _.algorithm == algorithm_names[idx] &&
                    _.num_ants == num_ants &&
                    _.num_correlators == num_correlators &&
                    _.GNSS == gnss
                ) |>
                @map(
                    {
                        _.Minimum
                    }
                ) |> DataFrame
        ).Minimum
    )
    lin[idx] = lines!(
        gpsl1_1ms_4_3_ax,
        freqs_gpsl1,
        times .* 10 ^ (-9)
    )
    sca[idx] = scatter!(
        gpsl1_1ms_4_3_ax,
        freqs_gpsl1,
        times .* 10 ^ (-9),
        marker = markers[idx],
        markersize = 15
    )
end
cputimes = Float64.(
        (
            raw_data_df |>
                @filter(
                    _.processor == "CPU" &&
                    _.num_ants == num_ants &&
                    _.num_correlators == num_correlators &&
                    _.GNSS == gnss
                ) |>
                @map(
                    {
                        _.Minimum
                    }
                ) |> DataFrame
        ).Minimum
    )
cpulin = lines!(
    gpsl1_1ms_4_3_ax,
    freqs_gpsl1,
    cputimes .* 10 ^ (-9),
)
cpusca = scatter!(
    gpsl1_1ms_4_3_ax,
    freqs_gpsl1,
    cputimes .* 10 ^ (-9),
    marker = :star5,
    markersize = 15
)
realtime = lines!(gpsl1_1ms_4_3_ax, [10^(6), 3 * 10^(8)], [10^(-3), 10^(-3)], linestyle = :dash, color = :dimgrey )
fig[1,2] = gpsl1_1ms_4_3_ax
fig
# PLOT GPS L1 C/A
gnss = "GPSL1"
num_ants = 4; num_correlators = 7

gpsl1_1ms_4_7_ax = Axis(fig,
    xlabel = "Sampling Frequency [Hz]",
    ylabel = "Processing Time [s]",
    xscale = log10,
    yscale = log10,
    xmajorgridvisible = true,
    xminorgridvisible = true,
    xminorticksvisible = true,
    xminorticks = IntervalsBetween(9),
    xticks = collect(10 .^ (6:9)),
    # yticks = collect(10.0 .^ (-6:-2)),
    title = "$(gnss_string[gnss]) T=1 ms, M=$num_ants, L=$num_correlators"
)
# xlims!(gpsl1_1ms_4_7_ax, [10^(6),  3 * 10^(8)])
# ylims!(gpsl1_1ms_4_7_ax, [10^(-6), 2*10^(-3)])
markers = [:rect, :utriangle, :circle, :diamond, :dtriangle]
lin = Array{Lines}(undef, length(algorithm_names));
sca = Array{Scatter}(undef, length(algorithm_names));
for (idx, algorithm) in enumerate(algorithm_names)
    times = Float64.(
        (
            raw_data_df |>
                @filter(
                    _.processor == "GPU" &&
                    _.algorithm == algorithm_names[idx] &&
                    _.num_ants == num_ants &&
                    _.num_correlators == num_correlators &&
                    _.GNSS == gnss
                ) |>
                @map(
                    {
                        _.Minimum
                    }
                ) |> DataFrame
        ).Minimum
    )
    lin[idx] = lines!(
        gpsl1_1ms_4_7_ax,
        freqs_gpsl1,
        times .* 10 ^ (-9)
    )
    sca[idx] = scatter!(
        gpsl1_1ms_4_7_ax,
        freqs_gpsl1,
        times .* 10 ^ (-9),
        marker = markers[idx],
        markersize = 15
    )
end
cputimes = Float64.(
        (
            raw_data_df |>
                @filter(
                    _.processor == "CPU" &&
                    _.num_ants == num_ants &&
                    _.num_correlators == num_correlators &&
                    _.GNSS == gnss
                ) |>
                @map(
                    {
                        _.Minimum
                    }
                ) |> DataFrame
        ).Minimum
    )
cpulin = lines!(
    gpsl1_1ms_4_7_ax,
    freqs_gpsl1,
    cputimes .* 10 ^ (-9),
)
cpusca = scatter!(
    gpsl1_1ms_4_7_ax,
    freqs_gpsl1,
    cputimes .* 10 ^ (-9),
    marker = :star5,
    markersize = 15
)
realtime = lines!(gpsl1_1ms_4_7_ax, [10^(6), 3 * 10^(8)], [10^(-3), 10^(-3)], linestyle = :dash, color = :dimgrey )
fig[1,3] = gpsl1_1ms_4_7_ax
fig

# PLOT GPS L5 C/A
gnss = "GPSL5"
num_ants = 1; num_correlators = 3
algorithm_names = unique(Vector{String}(raw_data_df[!, :algorithm])
)
sort!(algorithm_names)
gpsl5_1ms_1_3_ax = Axis(fig,
    xlabel = "Sampling Frequency [Hz]",
    ylabel = "Processing Time [s]",
    xscale = log10,
    yscale = log10,
    xmajorgridvisible = true,
    xminorgridvisible = true,
    xminorticksvisible = true,
    xminorticks = IntervalsBetween(9),
    # xticks = collect(10 .^ (6:9)),
    # yticks = collect(10.0 .^ (-6:-2)),
    title = "$(gnss_string[gnss]) T=1 ms, M=$num_ants, L=$num_correlators"
)
# xlims!(gpsl5_1ms_1_3_ax, [10^(7), 3 * 10^(8)])
# ylims!(gpsl5_1ms_1_3_ax,[1e-6, 5e-3])
markers = [:rect, :utriangle, :circle, :diamond, :dtriangle]
lin = Array{Lines}(undef, length(algorithm_names));
sca = Array{Scatter}(undef, length(algorithm_names));
for (idx, algorithm) in enumerate(algorithm_names)
    times = Float64.(
        (
            raw_data_df |>
                @filter(
                    _.processor == "GPU" &&
                    _.algorithm == algorithm_names[idx] &&
                    _.num_ants == num_ants &&
                    _.num_correlators == num_correlators &&
                    _.GNSS == gnss
                ) |>
                @map(
                    {
                        _.Minimum
                    }
                ) |> DataFrame
        ).Minimum
    )
    lin[idx] = lines!(
        gpsl5_1ms_1_3_ax,
        freqs_gpsl5,
        times .* 10 ^ (-9)
    )
    sca[idx] = scatter!(
        gpsl5_1ms_1_3_ax,
        freqs_gpsl5,
        times .* 10 ^ (-9),
        marker = markers[idx],
        markersize = 15
    )
end
cputimes = Float64.(
        (
            raw_data_df |>
                @filter(
                    _.processor == "CPU" &&
                    _.num_ants == num_ants &&
                    _.num_correlators == num_correlators &&
                    _.GNSS == gnss
                ) |>
                @map(
                    {
                        _.Minimum
                    }
                ) |> DataFrame
        ).Minimum
    )
cpulin = lines!(
    gpsl5_1ms_1_3_ax,
    freqs_gpsl5,
    cputimes .* 10 ^ (-9),
)
cpusca = scatter!(
    gpsl5_1ms_1_3_ax,
    freqs_gpsl5,
    cputimes .* 10 ^ (-9),
    marker = :star5,
    markersize = 15
)
realtime = lines!(gpsl5_1ms_1_3_ax, [2 * 10^(7), 3 * 10^(8)], [10^(-3), 10^(-3)], linestyle = :dash, color = :dimgrey )
fig[2,1] = gpsl5_1ms_1_3_ax
fig

# PLOT GPS L5 C/A
gnss = "GPSL5"
num_ants = 4; num_correlators = 3

gpsl5_1ms_4_3_ax = Axis(fig,
    xlabel = "Sampling Frequency [Hz]",
    ylabel = "Processing Time [s]",
    xscale = log10,
    yscale = log10,
    xmajorgridvisible = true,
    xminorgridvisible = true,
    xminorticksvisible = true,
    xminorticks = IntervalsBetween(9),
    # xticks = collect(10 .^ (6:9)),
    # yticks = collect(10.0 .^ (-6:-2)),
    title = "$(gnss_string[gnss]) T=1 ms, M=$num_ants, L=$num_correlators"
)
# xlims!(gpsl5_1ms_4_3_ax, [10^(7), 3 * 10^(8)])
# ylims!(gpsl5_1ms_4_3_ax,[1e-6, 5e-3])
markers = [:rect, :utriangle, :circle, :diamond, :dtriangle]
lin = Array{Lines}(undef, length(algorithm_names));
sca = Array{Scatter}(undef, length(algorithm_names));
for (idx, algorithm) in enumerate(algorithm_names)
    times = Float64.(
        (
            raw_data_df |>
                @filter(
                    _.processor == "GPU" &&
                    _.algorithm == algorithm_names[idx] &&
                    _.num_ants == num_ants &&
                    _.num_correlators == num_correlators &&
                    _.GNSS == gnss
                ) |>
                @map(
                    {
                        _.Minimum
                    }
                ) |> DataFrame
        ).Minimum
    )
    lin[idx] = lines!(
        gpsl5_1ms_4_3_ax,
        freqs_gpsl5,
        times .* 10 ^ (-9)
    )
    sca[idx] = scatter!(
        gpsl5_1ms_4_3_ax,
        freqs_gpsl5,
        times .* 10 ^ (-9),
        marker = markers[idx],
        markersize = 15
    )
end
cputimes = Float64.(
        (
            raw_data_df |>
                @filter(
                    _.processor == "CPU" &&
                    _.num_ants == num_ants &&
                    _.num_correlators == num_correlators &&
                    _.GNSS == gnss
                ) |>
                @map(
                    {
                        _.Minimum
                    }
                ) |> DataFrame
        ).Minimum
    )
cpulin = lines!(
    gpsl5_1ms_4_3_ax,
    freqs_gpsl5,
    cputimes .* 10 ^ (-9),
)
cpusca = scatter!(
    gpsl5_1ms_4_3_ax,
    freqs_gpsl5,
    cputimes .* 10 ^ (-9),
    marker = :star5,
    markersize = 15
)
realtime = lines!(gpsl5_1ms_4_3_ax, [2 * 10^(7), 3 * 10^(8)], [10^(-3), 10^(-3)], linestyle = :dash, color = :dimgrey )
fig[2,2] = gpsl5_1ms_4_3_ax
fig

# PLOT GPS L5 C/A
gnss = "GPSL5"
num_ants = 4; num_correlators = 7

gpsl5_1ms_4_7_ax = Axis(fig,
    xlabel = "Sampling Frequency [Hz]",
    ylabel = "Processing Time [s]",
    xscale = log10,
    yscale = log10,
    xmajorgridvisible = true,
    xminorgridvisible = true,
    xminorticksvisible = true,
    xminorticks = IntervalsBetween(9),
    # xticks = collect(10 .^ (6:9)),
    # yticks = collect(10.0 .^ (-6:-2)),
    title = "$(gnss_string[gnss]) T=1 ms, M=$num_ants, L=$num_correlators"
)
# xlims!(gpsl5_1ms_4_7_ax, [10^(7), 3 * 10^(8)])
# ylims!(gpsl5_1ms_4_7_ax,[1e-6, 5e-3])
markers = [:rect, :utriangle, :circle, :diamond, :dtriangle]
lin = Array{Lines}(undef, length(algorithm_names));
sca = Array{Scatter}(undef, length(algorithm_names));
for (idx, algorithm) in enumerate(algorithm_names)
    times = Float64.(
        (
            raw_data_df |>
                @filter(
                    _.processor == "GPU" &&
                    _.algorithm == algorithm_names[idx] &&
                    _.num_ants == num_ants &&
                    _.num_correlators == num_correlators &&
                    _.GNSS == gnss
                ) |>
                @map(
                    {
                        _.Minimum
                    }
                ) |> DataFrame
        ).Minimum
    )
    lin[idx] = lines!(
        gpsl5_1ms_4_7_ax,
        freqs_gpsl5,
        times .* 10 ^ (-9)
    )
    sca[idx] = scatter!(
        gpsl5_1ms_4_7_ax,
        freqs_gpsl5,
        times .* 10 ^ (-9),
        marker = markers[idx],
        markersize = 15
    )
end
cputimes = Float64.(
        (
            raw_data_df |>
                @filter(
                    _.processor == "CPU" &&
                    _.num_ants == num_ants &&
                    _.num_correlators == num_correlators &&
                    _.GNSS == gnss
                ) |>
                @map(
                    {
                        _.Minimum
                    }
                ) |> DataFrame
        ).Minimum
    )
cpulin = lines!(
    gpsl5_1ms_4_7_ax,
    freqs_gpsl5,
    cputimes .* 10 ^ (-9),
)
cpusca = scatter!(
    gpsl5_1ms_4_7_ax,
    freqs_gpsl5,
    cputimes .* 10 ^ (-9),
    marker = :star5,
    markersize = 15
)
realtime = lines!(gpsl5_1ms_4_7_ax, [2 * 10^(7), 3 * 10^(8)], [10^(-3), 10^(-3)], linestyle = :dash, color = :dimgrey )
fig[2,3] = gpsl5_1ms_4_7_ax
fig

labels = [algorithm_names; "SIMD CPU"]
entries = [];
for entry in 1:6
    entry == 6 ? entries = vcat(entries, [cpulin; cpusca]) : entries = vcat(entries, [lin[entry]; sca[entry]])
end

lfig =
Legend(
    fig[3,1:3],
    [
        [lin[1], sca[1]],
        [lin[2], sca[2]],
        [lin[3], sca[3]],
        [lin[4], sca[4]],
        [lin[5], sca[5]],
        [cpulin, cpusca],
    ],
    labels,
    framevisible = false
)

supertitle = Label(fig[0, :],
"Desktop PC, Intel Core i5-9600K @ 3.70 GHz, NVIDIA GeForce GTX 1050 Ti (Pascal)")
fig
save(plotsdir("desktop_allplots.pdf"), fig)

max_stars_repo_name: mfkiwl/GPUAcceleratedTracking

max_issues_repo_name: mfkiwl/GPUAcceleratedTracking

max_forks_repo_name: mfkiwl/GPUAcceleratedTracking


text:
############################################################################
# Copyright (c) Johan Mabille and Sylvain Corlay                           #
# Copyright (c) QuantStack                                                 #
#                                                                          #
# Distributed under the terms of the BSD 3-Clause License.                 #
#                                                                          #
# The full license is in the file LICENSE, distributed with this software. #
############################################################################

using Test

module TensorFunctions
using CxxWrap
using Xtensor
@wrapmodule Xtensor._l_tensors
end

# Test example functions from the TensorFunctions module

# example1
arr = [[1.0 2.0]
       [3.0 4.0]]
@test TensorFunctions.test_access(arr) == 4.0
@test TensorFunctions.example1([4., 5., 6.]) == 4.0

# example2
x = [[0. 1.]
     [2. 3.]]
res = [[2. 3.]
       [4. 5.]]
y = TensorFunctions.example2(x)
@test y ≈ res atol=1e-12

# vectorize_example1
x1 = [[0 1]
      [2 3]]
x2 = [0, 1]
res = [[0 2]
       [2 4]]
y = TensorFunctions.vectorize_example1(x1, x2)
@test y == res

# readme_example1
v = [[0.0  3.0  6.0   9.0  12.0]
     [1.0  4.0  7.0  10.0  13.0]
     [2.0  5.0  8.0  11.0  14.0]]

y = TensorFunctions.readme_example1(v)
@test y ≈ 1.2853996391883833 atol=1e-12

# readme_example2
x = [[ 0.0  1.0  2.0  3.0  4.0]
     [ 5.0  6.0  7.0  8.0  9.0]
     [10.0 11.0 12.0 13.0 14.0]]
y =  [1.0, 2.0, 3.0, 4.0, 5.0]
z = TensorFunctions.readme_example2(x, y)
expected = [[-0.540302  1.257618  1.89929   0.794764 -1.040465]
            [-1.499227  0.136731  1.646979  1.643002  0.128456]
            [-1.084323 -0.583843  0.45342   1.073811  0.706945]]
@test z ≈ expected atol=1e-5

# compare_shapes
x = ones(4, 4)
y = ones(5, 5)
z = zeros(4, 4)
@assert !TensorFunctions.compare_shapes(x, y)
@assert TensorFunctions.compare_shapes(x, z)


text:
function run_unbalanced(A::AbstractMatrix{T}) where T
    n = size(A,1)
    fu = eigen(A, scale=false)
    if T <: BlasFloat
        # use our code for balancing
        Abal,  B = balance!(copy(A))
        fb1 = eigen(Abal)
        fb = LinearAlgebra.Eigen(fb1.values, lmul!(B, fb1.vectors))
    else
        fb = eigen(A, scale=true)
    end
    resids(E) = [norm(A*E.vectors[:,j]-E.values[j]*E.vectors[:,j])/abs(E.values[j]) for j in 1:n]
    ru = maximum(resids(fu))
    rb = maximum(resids(fb))
    # println("residuals: bare $ru balanced $rb")
    ru, rb
end
@testset "balancing $T" for T in (Float64, ComplexF64,
                                  BigFloat, Complex{BigFloat})
    m = round(-log10(eps(real(T)))/3)
    A = T.([1 0 10.0^(-2m); 1 1 10.0^(-m); 10.0^(2m) 10.0^m 1])
    ru, rb = run_unbalanced(A)
    @test ru > rb
end


text:
"""
    POMDPFile

Writes a pomdpx file from a pomdp defined using the explicit interface of POMDPs.jl

# Constructors
- `POMDPFile(filename)` just stores the file name if the file exists
- `POMDPFile(pomdp::POMDP, filename="model.pomdpx"; verbose=true)` writes a .pomdpx file from the pomdp definition
"""
struct POMDPFile
    filename::AbstractString

    function POMDPFile(filename)
        @assert isfile(filename) "Pomdpx file $(filename) does not exist"
        return new(filename)
    end
end

function POMDPFile(pomdp::POMDP, filename="model.pomdpx"; verbose=true)
    pomdpx = POMDPXFile(filename)
    if verbose
        println("Generating a pomdpx file: $(filename)")
    end
    write(pomdp, pomdpx)
    return POMDPFile(filename)
end


text:
"""
    KeyAlignmentError

Is thrown when the constrained dimensions of components in a `KeyedDataset` have misaligned
key values.

# Fields
* constraint::Pattern - Constraint pattern describing all dimensions that must align
* groups - An iterator of paths and keys for each non-matching group
"""
struct KeyAlignmentError <: Exception
    constraint::Pattern
    groups
end

function Base.showerror(io::IO, err::KeyAlignmentError)
    lines = ["KeyAlignmentError: Misaligned dimension keys on constraint $(err.constraint)"]

    for (paths, keyvals) in err.groups
        push!(lines, string("  $paths ∈ ", sprint(summary, keyvals)))
    end

    println(io, join(lines, "\n"))
end

"""
    KeyedDataset

A `KeyedDataset` describes an associative collection of component `KeyedArray`s with constraints
on their shared dimensions.

# Fields
- `constraints::OrderedSet{Pattern}` - Constraint [`Pattern`](@ref)s on shared dimensions.
- `data::LittleDict{Tuple, KeyedArray}` - Flattened key paths as tuples component keyed arrays.
"""
@auto_hash_equals struct KeyedDataset
    # Our constraints are a collection of pseudo path tuples typically with 1 or
    # more `:_` wildcard components
    constraints::OrderedSet{Pattern}
    # Data lookup can be by any type, but typically it'll either be symbol or tuple.
    data::LittleDict{Tuple, KeyedArray}

    function KeyedDataset(
        constraints::OrderedSet{Pattern},
        data::LittleDict,
        check=true
    )
        ds = new(constraints, data)
        check && validate(ds)
        return ds
    end
end

function KeyedDataset(pairs::Pair...; constraints=Pattern[])
    # Convert any non-tuple keys to tuples
    tupled_pairs = map(pairs) do (k, v)
        k isa Tuple ? k => v : (k,) => v
    end

    data = LittleDict(tupled_pairs)

    # If no constraints have been specified then we default to (:__, dimname)
    constraint_set = if isempty(constraints)
        OrderedSet{Pattern}(
            Pattern(:__, d) for d in Iterators.flatten(dimnames.(values(data)))
        )
    else
        OrderedSet{Pattern}(constraints)
    end

    result = KeyedDataset(constraint_set, data)
    return result
end

# Utility kwargs and empty constructor.
function KeyedDataset(; constraints=Pattern[], kwargs...)
    if isempty(kwargs)
        return KeyedDataset(OrderedSet{Pattern}(constraints), LittleDict{Tuple, KeyedArray}())
    else
        return KeyedDataset(kwargs...; constraints=constraints)
    end
end

function Base.show(io::IO, ds::KeyedDataset)
    n = length(ds.data)
    m = length(ds.constraints)

    # Extract the constraints as a vector for indexing
    constraints = collect(ds.constraints)

    lines = String["KeyedDataset with:", "  $n components"]
    for (k, v) in ds.data
        # Identify shared dimensions where appropriate
        dimensions = map(dimnames(v)) do dimname
            cidx = findall(c -> (k..., dimname) in c, constraints)
            isempty(cidx) ? string(dimname) : string(dimname, "[", _only(cidx), "]")
        end

        s = string(
            "    $k => ",
            join(size(v), "x"),
            " $(nameof(typeof(v))){$(eltype(v))}",
            " with dimension ",
            join(dimensions, ", ")
        )

        push!(lines, s)
    end

    push!(lines, "  $m constraints")

    for (i, c) in enumerate(constraints)
        _axiskeys = axiskeys(ds, c)
        key_summary = isempty(_axiskeys) ? "NA" : sprint(summary, _only(_axiskeys))
        push!(lines, "    [$i] $(c.segments) ∈ $key_summary")
    end

    print(io, join(lines, "\n"))
end

"""
    dimpaths(ds, [pattern]) -> Vector{<:Tuple}

Return a list of all dimension paths in the [`KeyedDataset`](@ref).
Optionally, you can filter the results using a [`Pattern`](@ref).

# Example
```jldoctest
julia> using AxisKeys; using AxisSets: KeyedDataset, dimpaths;

julia> ds = KeyedDataset(
           :val1 => KeyedArray(rand(4, 3, 2); time=1:4, loc=-1:-1:-3, obj=[:a, :b]),
           :val2 => KeyedArray(rand(4, 3, 2) .+ 1.0; time=1:4, loc=-1:-1:-3, obj=[:a, :b]),
       );

julia> dimpaths(ds)
6-element Vector{Tuple{Symbol, Symbol}}:
 (:val1, :time)
 (:val1, :loc)
 (:val1, :obj)
 (:val2, :time)
 (:val2, :loc)
 (:val2, :obj)
```
"""
dimpaths(ds::KeyedDataset, pattern::Pattern) = filter(in(pattern), dimpaths(ds))
function dimpaths(ds::KeyedDataset)
    paths = Iterators.flatten(((k..., d) for d in dimnames(v)) for (k, v) in ds.data)
    return collect(paths)
end

"""
    constraintmap(ds)

Returns a mapping of constraint patterns to specific dimension paths.
The returned dictionary has keys of type [`Pattern`](@ref) and the values are sets of
`Tuple`.

# Example
```jldoctest
julia> using AxisKeys; using AxisSets: KeyedDataset, constraintmap;

julia> ds = KeyedDataset(
           :val1 => KeyedArray(rand(4, 3, 2); time=1:4, loc=-1:-1:-3, obj=[:a, :b]),
           :val2 => KeyedArray(rand(4, 3, 2) .+ 1.0; time=1:4, loc=-1:-1:-3, obj=[:a, :b]),
       );

julia> collect(constraintmap(ds))
3-element Vector{Pair{AxisSets.Pattern, Set{Tuple}}}:
 Pattern((:__, :time)) => Set([(:val2, :time), (:val1, :time)])
  Pattern((:__, :loc)) => Set([(:val1, :loc), (:val2, :loc)])
  Pattern((:__, :obj)) => Set([(:val2, :obj), (:val1, :obj)])
```
"""
function constraintmap(ds::KeyedDataset)
    items = dimpaths(ds)
    return LittleDict{Pattern, Set{Tuple}}(
        c => Set(filter(in(c), items)) for c in ds.constraints
    )
end

"""
    dimnames(ds)

Returns a list of the unique dimension names within the [`KeyedDataset`](@ref).

# Example
```jldoctest
julia> using AxisKeys; using NamedDims; using AxisSets: KeyedDataset;

julia> ds = KeyedDataset(
           :val1 => KeyedArray(rand(4, 3, 2); time=1:4, loc=-1:-1:-3, obj=[:a, :b]),
           :val2 => KeyedArray(rand(4, 3, 2) .+ 1.0; time=1:4, loc=-1:-1:-3, obj=[:a, :b]),
       );

julia> dimnames(ds)
3-element Vector{Symbol}:
 :time
 :loc
 :obj
```
"""
function NamedDims.dimnames(ds::KeyedDataset)
    return unique(Iterators.flatten(dimnames(a) for a in values(ds.data)))
end

"""
    axiskeys(ds)
    axiskeys(ds, dimname)
    axiskeys(ds, pattern)
    axiskeys(ds, dimpath)

Returns a list of unique axis keys within the [`KeyedDataset`](@ref).
A `Tuple` will always be returned unless you explicitly specify the `dimpath` you want.

# Example
```jldoctest
julia> using AxisKeys; using AxisSets: KeyedDataset;

julia> ds = KeyedDataset(
           :val1 => KeyedArray(rand(4, 3, 2); time=1:4, loc=-1:-1:-3, obj=[:a, :b]),
           :val2 => KeyedArray(rand(4, 3, 2) .+ 1.0; time=1:4, loc=-1:-1:-3, obj=[:a, :b]),
       );

julia> axiskeys(ds)
(1:4, -1:-1:-3, [:a, :b])

julia> axiskeys(ds, :time)
(1:4,)

julia> axiskeys(ds, (:val1, :time))
1:4
```
"""
function AxisKeys.axiskeys(ds::KeyedDataset)
    return Tuple(unique(Iterators.flatten(axiskeys(a) for a in values(ds.data))))
end

function AxisKeys.axiskeys(ds::KeyedDataset, dimpath::Tuple)
    key, dim = dimpath[1:end-1], dimpath[end]
    component = ds.data[key]
    return axiskeys(component, dim)
end

function AxisKeys.axiskeys(ds::KeyedDataset, pattern::Pattern)
    return Tuple(unique(axiskeys(ds, p) for p in dimpaths(ds, pattern)))
end

AxisKeys.axiskeys(ds::KeyedDataset, dim::Symbol) = axiskeys(ds, Pattern(:__, dim))

"""
    validate(ds, [constraint])

Validate that all constrained dimension paths within a [`KeyedDataset`](@ref) have matching key values.
Optionally, you can test an explicit constraint [`Pattern`](@ref).

# Returns
- `true` if an error isn't thrown

# Throws
- `ArgumentError`: If the constraints are not respected
"""
function validate(ds::KeyedDataset)
    for (k, v) in constraintmap(ds::KeyedDataset)
        validate(ds, k, v)
    end
    return true
end

function validate(ds::KeyedDataset, constraint::Pattern)
    paths = filter(in(constraint), dimpaths(ds))
    validate(ds, constraint, paths)
    return true
end

function validate(ds::KeyedDataset, constraint::Pattern, paths)
    if isempty(paths)
        @debug("No dimensions match the constraint $constraint")
    else
        keys = unique(axiskeys(ds, p) for p in paths)
        n = length(keys)
        if n > 1
            groups = Tuple(Vector{Tuple}() for i in 1:n)
            for p in paths
                for i in 1:n
                    if axiskeys(ds, p) == keys[i]
                        push!(groups[i], p)
                    end
                end
            end
            throw(KeyAlignmentError(constraint, zip(groups, keys)))
        end
    end
    return true
end


text:
#####################################################################
# Base functions.
#####################################################################
function sphere(x)
  sum(x.^2)
end

# Schwefel's ellipsoid.
function ellipsoid(x)
  res = 0
  for(i in 1:length(x))
    res += sum(x[1:i])^2
  end
  res
end

function elliptic(x)
  D = length(x)
  condition = 1e+6
  coefficients = condition .^ linspace(0, 1, D)
  sum(coefficients .* x.^2)
end

function rastrigin(x)
  D = length(x)
  10 * D + sum( x.^2 ) - 10 * sum( cos( 2 * π * x ) )
end

function ackley(x)
  D = length(x)
  try
    20 - 20.*exp(-0.2.*sqrt(sum(x.^2)/D)) - exp(sum(cos(2 * π * x))/D) + e
  catch e
    # Sometimes we have gotten a DomainError from the cos function so we protect this call
    println(e)
    println("For input x = ", x)
    # Return a very large fitness value to indicate that this is NOT the solution we want.
    # TODO: Fix better way to handle this!
    9.99e100
  end
end

function schwefel1_2(x)
  D = length(x)
  partsums = zeros(D)
  partsum = 0
  for i in 1:D
    partsum += x[i]
    partsums[i] = partsum
  end
  sum(partsums.^2)
end

function rosenbrock(x)
  n = length(x)
  return( sum( 100*( x[2:n] - x[1:(n-1)].^2 ).^2 + ( x[1:(n-1)] - 1 ).^2 ) )
end

function step(x)
  sum(ceil(x + 0.5))
end

function griewank(x)
  n = length(x)
  1 + (1/4000)*sum(x.^2) - prod(cos(x ./ sqrt(1:n)))
end

function schwefel2_22(x)
  ax = abs(x)
  sum(ax) + prod(ax)
end

function schwefel2_21(x)
  maximum(abs(x))
end

# I'm unsure about this one since it does not return the expected minima at
# [1.0, 1.0].
function schwefel2_26(x)
  D = length(x)
  418.98288727243369 * D - sum(x .* sin(sqrt(abs(x))))
end

function cigar(x)
  x[1]^2 + 1e6 * sum(x[2:end].^2)
end

function cigtab(x)
  x[1]^2 + 1e8 * x[end]^2 + 1e4 * sum(x[2:(end-1)].^2)
end

# Shekel10 is a 4D, multi-minima, non-separable test problem. Our implementation
# is based on the C code in:
#   http://www.math.ntu.edu.tw/~wwang/cola_lab/test_problems/multiple_opt/multiopt_prob/Shekel10/Shekel10.c
Shekel10_A = [4 4 4 4; 1 1 1 1; 8 8 8 8; 6 6 6 6; 3 7 3 7; 2 9 2 9; 5 5 3 3; 8 1 8 1; 6 2 6 2; 7 3.6 7 3.6]
Shekel10_C = [0.1, 0.2, 0.2, 0.4, 0.4, 0.6, 0.3, 0.7, 0.5, 0.5]
function shekel(x, a, c)
  sum = 0.0
  for i in 1:length(c)
    den = 0.0
    for j in 1:size(a, 2)
      den += (x[j] - a[i,j])^2
    end
    sum = sum - 1 / (den + c[i])
  end
  return sum
end
shekel10(x) = shekel(x, Shekel10_A, Shekel10_C)

# Shekel7 is a 4D, multi-minima, non-separable test problem. Our implementation
# is based on the C code in:
#   http://www.math.ntu.edu.tw/~wwang/cola_lab/test_problems/multiple_opt/multiopt_prob/Shekel7/Shekel7.c
Shekel7_A = [4 4 4 4; 1 1 1 1; 8 8 8 8; 6 6 6 6; 3 7 3 7; 2 9 2 9; 5 5 3 3]
Shekel7_C = [0.1, 0.2, 0.2, 0.4, 0.4, 0.6, 0.3]
shekel7(x) = shekel(x, Shekel7_A, Shekel7_C)

# Shekel5 is a 4D, multi-minima, non-separable test problem. Our implementation
# is based on the C code in:
#   http://www.math.ntu.edu.tw/~wwang/cola_lab/test_problems/multiple_opt/multiopt_prob/Shekel5/Shekel5.c
Shekel5_A = [4 4 4 4; 1 1 1 1; 8 8 8 8; 6 6 6 6; 3 7 3 7]
Shekel5_C = [0.1, 0.2, 0.2, 0.4, 0.4]
shekel5(x) = shekel(x, Shekel5_A, Shekel5_C)

# Hartman 6D is a multi-minima, non-separable test problem. Our implementation is based on:
#  http://www.sfu.ca/~ssurjano/hart6.html
Hartman6_alpha = [1.0 1.2 3.0 3.2]
Hartman6_A = [10 3 17 3.50 1.7 8; 0.05 10 17 0.1 8 14; 3 3.5 1.7 10 17 8; 17 8 0.05 10 0.1 14]
Hartman6_P = 1e-4 * [1312 1696 5569 124 8283 5886; 2329 4135 8307 3736 1004 9991; 2348 1451 3522 2883 3047 6650; 4047 8828 8732 5743 1091 381]
function hartman(x, alpha, A, P)
  sum = 0.0
  for i in 1:length(alpha)
    isum = 0.0
    for j in 1:size(A, 2)
      isum += A[i,j] * (x[j] - P[i,j])^2
    end
    sum -= alpha[i] * exp(-isum)
  end
  sum
end
hartman6(x) = hartman(x, Hartman6_alpha, Hartman6_A, Hartman6_P)

# Hartman 3D is a multi-minima, non-separable test problem. Our implementation is based on:
#  http://www.sfu.ca/~ssurjano/hart3.html
# However, we get a different global minima than the one stated on that page.
Hartman3_alpha = [1.0 1.2 3.0 3.2]
Hartman3_A = [3.0 10 30; 0.1 10 35; 3.0 10 30; 0.1 10 36]
Hartman3_P = 1e-4 * [3689 1170 2673; 4699 4387 7470; 1091 8732 5547; 381 5743 8828]
hartman3(x) = hartman(x, Hartman3_alpha, Hartman3_A, Hartman3_P)
# The minima should be -3.86278 but we get a different one so use that in the problem spec:

#####################################################################
# S2 functions in addition to the base functions above. As stated
# in table II of the JADE paper: http://150.214.190.154/EAMHCO/pdf/JADE.pdf
#####################################################################

function quartic(x)
  D = length(x)
  sum( (1:D) .* x.^4 )
end

function noisy_quartic(x)
  quartic(x) + rand()
end

function s2_step(x)
  sum( ceil(x + 0.5).^2 )
end

#####################################################################
# Misc other interesting optimization functions and families.
#####################################################################

# Schwefel 2.13 is a hard one...
#function f=schwefel_213(x)
#global initial_flag
#persistent a b A alpha
#[ps,D]=size(x);
#if initial_flag==0
#    initial_flag=1;
#    load schwefel_213_data
#    if length(alpha)>=D
#        alpha=alpha(1:D);a=a(1:D,1:D);b=b(1:D,1:D);
#    else
#        alpha=-3+6*rand(1,D);
#        a=round(-100+200.*rand(D,D));
#        b=round(-100+200.*rand(D,D));
#    end
#    alpha=repmat(alpha,D,1);
#    A=sum(a.*sin(alpha)+b.*cos(alpha),2);
#end
#for i=1:ps
#    xx=repmat(x(i,:),D,1);
#    B=sum(a.*sin(xx)+b.*cos(xx),2);
#    f(i,1)=sum((A-B).^2,1);
#end

# This is a generator for the family of deceptive functions from the
# Cuccu2011 paper on novelty-based restarts. We have vectorized it to allow
# more than 1D versions. The Cuccu2011 paper uses the following values for
# (l, w) = [(5, 0),  (15, 0),  (30, 0),
#           (5, 2),  (15, 2),  (30, 2),
#           (5, 10), (15, 10), (30, 10)]
# and notes that (15, 2) and (30, 2) are the most difficult instances.
function deceptive_cuccu2011(l, w)
  (x) -> begin
    absx = abs(x)
    sumabsx = sum(absx)
    if sumabsx <= 1
      return sum(x.^2)
    elseif sumabsx >= l+1
      return sum((absx - l).^2)
    else
      return (1 - 0.5 * sum(sin( (π * w * (absx - 1)) / l ).^2))
    end
  end
end

# Deceptive/hardest instances:
deceptive_cuccu2011_15_2 = deceptive_cuccu2011(15, 2)
deceptive_cuccu2011_30_2 = deceptive_cuccu2011(30, 2)


text:
module QuDynamics
    using QuBase
    using Compat
    using ODE
    using Expokit
    using ExpmV
    VERSION < v"0.4-" && using Docile
    include("quequations.jl")
    include("propmachinery.jl")
    include("propstepsolvers.jl")
    include("propodesolvers.jl")
    include("propexpmsolvers.jl")
    include("propmcwf.jl")
    module QuTiP
    using ..QuBase
    using ..QuDynamics
    include("qutipinterface.jl")
    end
    export propagate
end


text:
# function fit(::Type{Exponential}, m::AbstractMoments)
#     # https://en.wikipedia.org/wiki/Exponential_distribution
#     n_moments(m) >= 1 || error("Need mean to estimate exponential")
#     return Exponential(mean(m))
# end

function fit(::Type{Weibull}, lower::QuantilePoint, upper::QuantilePoint)
    # https://www.johndcook.com/quantiles_parameters.pdf
    gamma = (log(-log(1-upper.p)) - log(-log(1-lower.p)))/(log(upper.q) -log(lower.q))
    beta = lower.q / (-log(1-lower.p))^(1/gamma)
    Weibull(gamma, beta)
end

# function fit_mean_quantile(::Type{Exponential}, mean::Real, qp::QuantilePoint)
#     # only fit to mean
#     fit(Type{Exponential}, AbstractMoments(mean))
# end

# function fit_mode_quantile(::Type{Exponential}, mode::Real, qp::QuantilePoint)
#     # ignore mode (its always at 0)
#     θ = -qp.q/log(1-qp.p)
#     Exponential(θ)
# end


text:
"""
`r2_score(x, y)`

Calculate the coefficient of determination as described by Nash and Sutcliffe
(https://doi.org/10.1016/0022-1694(70)90255-6).

**Arguments**
- `x`: Array of target values (data set).
- `y`: Array of predicted (fitted, modeled) values.

**Returns**
- `R2`: The coefficient of determination.
"""
function r2_score(x, y)
    μ = mean(x)
    SSres = sum((x .- y).^2)
    SStot = sum((x .- μ).^2)
    R2 = 1 - SSres/SStot
    return R2
end

"""
`ptp(x; kwargs...)`

Return the peak-to-peak range (maximum - minimum).

**Arguments**
- `x`: Array of values.

**Keyword Arguments**
- Passed to `maximum`/`minimum` functions. See their documentation.

**Returns**
- Peak-to-peak range.
"""
function ptp(x; kwargs...)
    return maximum(x; kwargs...) - minimum(x; kwargs...)
end

max_stars_repo_name: vkumpost/stoosc

max_issues_repo_name: vkumpost/stoosc

max_forks_repo_name: vkumpost/stoosc


text:
export glscreen

let screen_list = Screen[]
    global current_screen, add_screen, get_screens, empty_screens!
    current_screen() = last(screen_list)
    add_screen(screen) = push!(screen_list, screen)
    get_screens() = copy(screen_list)
    empty_screens!() = empty!(screen_list)
end

function cleanup_old_screens()
    for screen in get_screens()
        destroy!(screen)
    end
    empty!(timer_signal_dict)
    empty_screens!()
    GLFW.Terminate()
    GLFW.Init()
    reset_texture_atlas!()
    GLAbstraction.empty_shader_cache!()
end

function glscreen(name="GLVisualize";
        resolution = GLWindow.standard_screen_resolution(),
        debugging = false,
        background = RGBA(1,1,1,1)
    )

    cleanup_old_screens()

    screen = Screen(name, resolution=resolution, debugging=debugging, color=background)
    add_screen(screen)
    GLWindow.add_complex_signals!(screen) #add the drag events and such
    preserve(map(screen.inputs[:window_open]) do open
        if !open
            empty_screens!()
        end
        nothing
    end)
    screen
end

const timer_signal_dict = Dict{Int, WeakRef}()
"""
Creates a timer signal with `updates_per_second` while `window` is open.
It's reusing timer signals with the same update rate and registering the updates
with GLFW.
"""
function get_timer_signal(updates_per_second, window=current_screen())
    signal = get!(timer_signal_dict, updates_per_second) do
        # because this is a function, it'll only get executed if needed
        WeakRef(fpswhen(window.inputs[:window_open], updates_per_second))
    end.value
    signal
end

function fold_loop(v0, _)
    val, range, index = v0
    val = range[index]
    index += 1
    index>length(range) && (index = 1)
    (val, range, index)
end

function loop(range::Range, rate=60)
    t = get_timer_signal(rate)
    map(first, foldp(fold_loop, (first(range), range, 1), t))
end

function fold_bounce(v0, _)
    val, range, index, direction = v0
    val = range[index]
    index += direction
    if index in (length(range)+1, 0)
        direction = -direction
        index += 2direction
    end
    (val, range, index, direction)
end

function bounce{T}(range::Range{T}, rate=60)
    t = get_timer_signal(rate)
    map(first, foldp(fold_bounce, (first(range), range, 1, 1), t))
end

export bounce, loop


text:
struct GaussianCopula <: Copula
    correlation::Matrix{<:Real}
end

const Φ = Normal()

function sample(c::GaussianCopula, n::Integer=1)
    L = cholesky(c.correlation).L

    Z = rand(Φ, n, size(c.correlation, 2))
    X = transpose(L * transpose(Z))

    return u = cdf.(Φ, X)
end

function dimensions(c::GaussianCopula)
    return size(c.correlation, 1)
end

function to_standard_normal_space(c::GaussianCopula, u)
    L = cholesky(c.correlation).L
    return quantile.(Φ, u) * inv(L)'
end

function to_copula_space(c::GaussianCopula, s)
    L = cholesky(c.correlation).L
    return cdf.(Φ, s * L')
end

max_stars_repo_name: FriesischScott/julia-uq

max_issues_repo_name: FriesischScott/julia-uq

max_forks_repo_name: FriesischScott/julia-uq


text:
module NNlib

using Requires

export σ, sigmoid, relu, leakyrelu, elu, swish, selu, softplus, softsign, logσ, logsigmoid,
  softmax, logsoftmax, maxpool, meanpool

include("numeric.jl")
include("activation.jl")
include("softmax.jl")
include("logsoftmax.jl")
include("linalg.jl")
include("conv.jl")
include("cubroadcast.jl")

end # module


text:
using mytest
using Test

@testset "mytest.jl" begin
    # Write your own tests here.
end


text:
include("utils_dg.jl")

K = 7  # Number of elements

n = 4
nn = 1
r = jacobiGL(0, 0, n)
V = vandermonde(r, 0, 0, n)
fr = collect(range(-1,1,length=100))
fV = vandermonde(fr, 0, 0, n)
refine = fV * inv(V)
# cell average
avg = Diagonal(zeros(n+1))
avg[1:nn] .= 1
p0 = V * avg * inv(V)
# polyone
nn = 2
tmp = zeros(n+1)
tmp[1:nn] .= 1
avg = Diagonal(tmp)
p1 = V * avg * inv(V)
# poly two
nn = 2+1
tmp = zeros(n+1)
tmp[1:nn] .= 1
avg = Diagonal(tmp)
p2 = V * avg * inv(V)

# poly 6
Ω = (; a=0, b=2π, periodic = true)
a = Ω.a
b = Ω.b
mesh = Mesh(K, n, Ω.a, Ω.b, periodic = Ω.periodic)
x = mesh.x
weights = sum(mesh.M, dims = 1)
Δx = x[2] - x[1]

# inexact integration
DM = Diagonal(sum(mesh.M, dims = 1)[:])
mesh.M .= DM
mesh.Mi .= inv(DM)
mesh.lift[:,1] .= mesh.Mi[1,:]
mesh.lift[:,end] .= mesh.Mi[end,:]

b = Ω.b
a = Ω.a

u_first  = @. sin(2π/(b-a) * x) *( 1+ 0.9 * sin(8 * 2π/(b-a) * x) )
using Random
Random.seed!(123)
uo = refine * u_first
x = refine * x

Duu  = refine * compute_volume_terms(u_first .* u_first, mesh)
uDu  = 2.0 .* refine * ( u_first .* compute_volume_terms(u_first, mesh) )
ESDu = 2/3 .* Duu + 1/3 .* uDu

## Plotting
fig = Figure(resolution = (1700+600, 1000+400))

lims = (extrema(x)..., extrema(uo)...)
axo = Axis(fig[1,1], xlabel = L"x", ylabel = L"y", ylabelsize = 22,
    xlabelsize= 22, xgridstyle=:dash, ygridstyle=:dash, xtickalign = 1,
    xticksize=10, ytickalign=1, yticksize=10,  xlabelpadding = -10)

axo.limits = (extrema(x)..., (extrema(uo) .* 1.1)...)
axo.title = "Function"
axo.titlesize = 32
for i in 1:K
    lines!(axo, x[:,i], uo[:,i], linewidth = 5)
end

ax0 = Axis(fig[1,2], xlabel = L"x", ylabel = L"y", ylabelsize = 22,
    xlabelsize= 22, xgridstyle=:dash, ygridstyle=:dash, xtickalign = 1,
    xticksize=10, ytickalign=1, yticksize=10,  xlabelpadding = -10, ylims = (-1.1, 1.1))
lims2 = (extrema(x)..., (extrema(Duu) .* 1.1)...)
ax0.limits = lims2
ax0.title = "Conservative "
ax0.titlesize = 32
for i in 1:K
    lines!(ax0, x[:,i], Duu[:,i], linewidth = 5)
end

ax1 = Axis(fig[2,1], xlabel = L"x", ylabel = L"y", ylabelsize = 22,
    xlabelsize= 22, xgridstyle=:dash, ygridstyle=:dash, xtickalign = 1,
    xticksize=10, ytickalign=1, yticksize=10,  xlabelpadding = -10, ylims = (-1.1, 1.1))
ax1.limits = lims2
ax1.title = "Advective"
ax1.titlesize = 32
for i in 1:K
    lines!(ax1, x[:,i], uDu[:,i], linewidth = 5)
end

ax3 = Axis(fig[2,2], xlabel = L"x", ylabel = L"y", ylabelsize = 22,
    xlabelsize= 22, xgridstyle=:dash, ygridstyle=:dash, xtickalign = 1,
    xticksize=10, ytickalign=1, yticksize=10,  xlabelpadding = -10, ylims = (-1.1, 1.1))
ax3.limits = lims2
ax3.title = "Entropy Stable"
ax3.titlesize = 32
for i in 1:K
    lines!(ax3, x[:,i], ESDu[:,i], linewidth = 5)
end

max_stars_repo_name: sandreza/HeldSuarezVisualizationScripts

max_issues_repo_name: sandreza/HeldSuarezVisualizationScripts

max_forks_repo_name: sandreza/HeldSuarezVisualizationScripts