text:
"""
    AbstractDictionary{I, T}

Abstract type for a dictionary between unique indices of type `I` to elements of type `T`.

At minimum, an `AbstractDictionary` must implement:

 * `getindex(::AbstractDictionary{I, T}, ::I) --> T`
 * `isassigned(::AbstractDictionary{I}, ::I) --> Bool`
 * `keys(::AbstractDictionary{I, T}) --> AbstractIndices{I}`

If values can be set/mutated, then an `AbstractDictionary` should implement:

 * `issettable(::AbstractDictionary)` (returning `true`)
 * `setindex!(dict::AbstractDictionary{I, T}, ::T, ::I}` (returning `dict`)

If arbitrary indices can be added to or removed from the dictionary, implement:

 * `isinsertable(::AbstractDictionary)` (returning `true`)
 * `insert!(dict::AbstractDictionary{I, T}, ::I, ::T}` (returning `dict`)
 * `delete!(dict::AbstractDictionary{I, T}, ::I}` (returning `dict`)
"""
abstract type AbstractDictionary{I, T}; end
abstract type AbstractIndices{I} <: AbstractDictionary{I, I}; end

Base.eltype(d::AbstractDictionary) = eltype(typeof(d))
Base.eltype(::Type{<:AbstractDictionary{I, T}}) where {I, T} = T
Base.keytype(d::AbstractDictionary) = keytype(typeof(d))
Base.keytype(::Type{<:AbstractDictionary{I, T}}) where {I, T} = I

function Base.keys(dict::AbstractDictionary)
    error("Every AbstractDictionary type must define a method for `keys`: $(typeof(dict))")
end

@propagate_inbounds function Base.getindex(dict::AbstractDictionary{I}, i::I) where {I}
    if !(istokenizable(dict))
        error("Every AbstractDictionary type must define a method for `getindex`: $(typeof(dict))")
    end

    (hastoken, token) = gettoken(dict, i)
    @boundscheck if !hastoken
        throw(IndexError("Dictionary does not contain index: $i"))
    end
    return @inbounds gettokenvalue(dict, token)
end

@inline function Base.iterate(dict::AbstractDictionary{I, T}, s...) where {I, T}
    if istokenizable(dict)
        tmp = iteratetoken(keys(dict), s...)
        tmp === nothing && return nothing
        (t, s2) = tmp
        return (gettokenvalue(dict, t), s2)
    else
        tmp = iterate(keys(dict), s...)
        tmp === nothing && return nothing
        (i, s2) = tmp
        return (@inbounds dict[i], s2)
    end
end

Base.length(d::AbstractDictionary) = length(keys(d))
Base.IteratorSize(d::AbstractDictionary) = Base.IteratorSize(keys(d))

Base.checkindex(d::AbstractDictionary{I}, i) where {I} = checkindex(d, convert(I, i))
Base.checkindex(d::AbstractDictionary{I}, i::I) where {I} = checkindex(keys(d), i)

## Comparisons

# dictionaries are isequal if they iterate in the same order
function Base.isequal(d1::AbstractDictionary, d2::AbstractDictionary)
    if d1 === d2
        return true
    end

    if sharetokens(d1, d2)
        @inbounds for t in tokens(d1)
            if !isequal(gettokenvalue(d1, t), gettokenvalue(d2, t))
                return false
            end
        end
        return true
    end

    if length(d1) != length(d2)
        return false
    end

    for (p1, p2) in zip(pairs(d1), pairs(d2))
        if !isequal(p1, p2)
            return false
        end
    end

    return true
end

# The indices must be isequal and the values ==, same ordering
function Base.:(==)(d1::AbstractDictionary, d2::AbstractDictionary)
    out = true

    if sharetokens(d1, d2)
        @inbounds for t in tokens(d1)
            out &= gettokenvalue(d1, t) == gettokenvalue(d2, t)
            if out === false
                return false
            end
        end
        return out
    end

    if length(d1) != length(d2)
        return false
    end

    for ((i1,x1), (i2,x2)) in zip(pairs(d1), pairs(d2))
        if !isequal(i1, i2)
            return false
        end
        out &= x1 == x2 # make sure it works for `missing`
        if out === false
            return false
        end
    end

    return out
end

"""
    isdictequal(d1, d2)

Determine if two dictionaries are equivalent. Dictionaries `d1` and `d2` are equivalent if
`issetequal(keys(d1), keys(d2))` and for each key `i`, `d1[i] == d2[i]`.

Example

```julia
julia> isdictequal(Dictionary(['a','b'],[1,2]), Dictionary(['b','a'],[2,1]))
true

julia> isdictequal(Dictionary(['a','b'],[1,2]), Dictionary(['b','a'],[2,3]))
false

julia> isdictequal(Dictionary(['a','b'],[1,2]), Dictionary(['b','a'],[2,missing]))
missing
```
"""
function isdictequal(d1::AbstractDictionary, d2::AbstractDictionary)
    out = true

    if sharetokens(d1, d2)
        @inbounds for t in tokens(d1)
            out &= gettokenvalue(d1, t) == gettokenvalue(d2, t)
            if out === false
                return false
            end
        end
        return out
    end

    if length(d1) != length(d2)
        return false
    end

    for (i,x1) in pairs(d1)
        (hastoken, t) = gettoken(d2, i)
        if !hastoken
            return false
        end
        out &= x1 == gettokenvalue(d2, t) # make sure it works for `missing`
        if out === false
            return false
        end
    end

    return out
end

# Lexical ordering based on iteration (of pairs - lesser key takes priority over lesser value, as implmeneted in `cmp(::Pair)`)
function Base.isless(dict1::AbstractDictionary, dict2::AbstractDictionary)
    if sharetokens(dict1, dict2)
        @inbounds for t in tokens(dict1)
            c = cmp(gettokenvalue(dict1, t), gettokenvalue(dict2, t))
            if c == -1
                return true
            elseif c == 1
                return false
            end
        end
        return false # they are isequal
    end

    # We want to iterate... until one is longer than the other (no `zip`)
    pairdict1 = pairs(dict1)
    pairdict2 = pairs(dict2)
    tmp1 = iterate(pairdict1)
    tmp2 = iterate(pairdict2)
    while tmp1 !== nothing
        if tmp2 === nothing
            return false # shorter collections are isless in lexical ordering
        end
        (p1, s1) = tmp1
        (p2, s2) = tmp2
        c = cmp(p1, p2)

        if c == -1
            return true
        elseif c == 1
            return false
        end

        tmp1 = iterate(pairdict1, s1)
        tmp2 = iterate(pairdict2, s2)
    end
    return tmp2 !== nothing # shorter collections are isless in lexical ordering
end

function Base.cmp(dict1::AbstractDictionary, dict2::AbstractDictionary)
    if sharetokens(dict1, dict2)
        @inbounds for t in tokens(dict1)
            c = cmp(gettokenvalue(dict1, t), gettokenvalue(dict2, t))
            if c == -1
                return -1
            elseif c == 1
                return 1
            end
        end
        return 0
    end

    # We want to iterate... until one is longer than the other (no `zip`)
    pairdict1 = pairs(dict1)
    pairdict2 = pairs(dict2)
    tmp1 = iterate(pairdict1)
    tmp2 = iterate(pairdict2)
    while tmp1 !== nothing
        if tmp2 === nothing
            return 1 # shorter collections are isless in lexical ordering
        end
        (p1, s1) = tmp1
        (p2, s2) = tmp2
        c = cmp(p1, p2)

        if c == -1
            return -1
        elseif c == 1
            return 1
        end

        tmp1 = iterate(pairdict1, s1)
        tmp2 = iterate(pairdict2, s2)
    end
    if tmp2 === nothing # shorter collections are isless in lexical ordering
        return 0
    else
        return -1
    end
end

## Hashing

function Base.hash(dict::AbstractDictionary, h::UInt)
    h1 = h
    h2 = h
    for (i, v) in pairs(dict)
        h1 = hash(i, h1)
        h2 = hash(v, h2)
    end

    return hash(hash(UInt === UInt64 ? 0x8955a87bc313a509 : 0xa9cff5d1, h1), h2)
end

## unique returns an Indices

function Base.unique(d::AbstractDictionary)
    out = empty(d, eltype(d)) # an AbstractIndices
    for x in d
        set!(out, x)
    end
    return out
end

function _distinct(f, ::Type{T}, itr) where T
    out = T()
    for x in itr
        i = f(x)
        (hadtoken, token) = gettoken!(out, x)
        if !hadtoken
            @inbounds settokenvalue!(out, token, i)
        end
    end
    return out
end

### Settable interface

"""
    issettable(dict::AbstractDictionary)

Return `true` if the dictionary `dict` obeys the settable interface, or `false` otherwise.

A mutable dictionary is one where the *values* can be modified (but not necessarily the
indices). The mutable interface requires the dictionary to implement:

 * `setindex!(dict::AbstractDictionary{I, T}, value::I, index::T)`

New settable dictionaries are primarily created through the `similar` function (for
unitialized values), as well as `fill`, `zeros`, `ones`, `trues` and `falses` (for
initialized values).

See also `isinsertable`.
"""
issettable(::AbstractDictionary) = false

function Base.isassigned(dict::AbstractDictionary{I}, i::I) where {I}
    if !(istokenizable(dict))
        error("Every settable AbstractDictionary type must define a method for `isassigned`: $(typeof(dict))")
    end

    (hasindex, token) = gettoken(dict, i)
    if hasindex
        return istokenassigned(dict, token)
    else
        return false
    end
end

@propagate_inbounds function Base.setindex!(dict::AbstractDictionary{I, T}, value::T, i::I) where {I, T}
    if !(istokenizable(dict))
        error("Every settable AbstractDictionary type must define a method for `setindex!`: $(typeof(dict))")
    end

    (hastoken, token) = gettoken(dict, i)
    @boundscheck if !hastoken
        throw(IndexError("Dictionary does not contain index: $i"))
    end
    settokenvalue!(dict, token, value)
    return dict
end

# similar is the primary way to construct an `issettable` dictionary, and allows us to write
# generic "factories" like `fill`, `copy`, `zeros`...

"""
    similar(d::AbstractDictionary, [T=eltype(d)])

Construct a new `issettable` dictionary with identical `keys` as `d` and an element type of
`T`. The initial values are3unitialized/undefined.
"""
Base.similar(d::AbstractDictionary) = similar(keys(d), eltype(d))
Base.similar(d::AbstractDictionary, ::Type{T}) where {T} = similar(keys(d), T)

function Base.similar(indices::AbstractIndices{I}, ::Type{T}) where {I, T}
    return similar(convert(Indices{I}, indices), T)
end

function Base.merge(d1::AbstractDictionary, d2::AbstractDictionary)
    # Note: need to copy the keys
    out = similar(copy(keys(d1)), eltype(d1))
    copyto!(out, d1)
    merge!(out, d2)
    return out
end

if isdefined(Base, :mergewith) # Julia 1.5+
    function Base.mergewith(combner, d1::AbstractDictionary, d2::AbstractDictionary)
        # Note: need to copy the keys
        out = similar(copy(keys(d1)), eltype(d1))
        copyto!(out, d1)
        mergewith!(combner, out, d2)
        return out
    end
end

# fill! and fill

function Base.fill!(d::AbstractDictionary, value)
    for t in tokens(d)
        settokenvalue!(d, t, value)
    end
    return d
end

"""
    fill(value, d::AbstractDictionary, [T = typeof(value)])

Construct a new `issettable` dictionary with identical `keys` as `d` and all elements
initialized to `value`.

An element type can optionally be provided, which can be useful for constructing containers
that accept different types of values, for example `fill(d, missing, Union{Missing, Bool})`.
"""
Base.fill(value, d::AbstractDictionary) = fill(value, d, typeof(value))
function Base.fill(value, indices::AbstractDictionary, ::Type{T}) where {T}
    out = similar(keys(indices), T)
    fill!(out, value)
    return out
end

# zeros, ones, falses, trues
Base.zeros(d::AbstractDictionary) = zeros(Float64, d)
Base.zeros(::Type{T}, d::AbstractDictionary) where {T} = fill(zero(T), d, T)

Base.ones(d::AbstractDictionary) = ones(Float64, d)
Base.ones(::Type{T}, d::AbstractDictionary) where {T} = fill(one(T), d, T)

Base.falses(d::AbstractDictionary) = falses(Bool, d)
Base.falses(::Type{T}, d::AbstractDictionary) where {T} = fill(false, d, T) # T could be `Union{Missing, Bool}`, for example

Base.trues(d::AbstractDictionary) = trues(Bool, d)
Base.trues(::Type{T}, d::AbstractDictionary) where {T} = fill(true, d, T) # T could be `Union{Missing, Bool}`, for example

# rand

# Drawing FROM a dictionary requires is to create a `Sampler(rng, dict, Val(1/Inf))` and overload `rand(rng, S)`
function Random.rand(rng::AbstractRNG, dict::AbstractDictionary)
    Random.rand(rng::AbstractRNG, Random.Sampler(rng, dict, Val(Inf)))
end

function Random.rand(rng::AbstractRNG, s::Random.SamplerTrivial{<:AbstractDictionary})
    dict = s[]
    inds = keys(dict)
    t = randtoken(rng, inds)
    return @inbounds gettokenvalue(dict, t)
end

# Creating a random dictionary

function Random.rand!(dict::AbstractDictionary; S = Random.Sampler(Random.GLOBAL_RNG, eltype(dict), Val(1)))
    map!(() -> rand(Random.GLOBAL_RNG, S), dict)
    return dict
end

function Random.rand!(rng::AbstractRNG, dict::AbstractDictionary; S = Random.Sampler(rng, eltype(dict), Val(1)))
    map!(() -> rand(rng, S), dict)
    return dict
end

sampletype(::Random.Sampler{T}) where {T} = T

Random.rand(S, dict::AbstractDictionary) = rand(Random.GLOBAL_RNG, S, dict)

function Random.rand(rng::AbstractRNG, S, dict::AbstractDictionary)
    sampler = Random.Sampler(rng, S, Val(Inf))
    out = similar(dict, sampletype(sampler))
    rand!(rng, out, S=sampler)
    return out
end

# randn
Random.randn!(dict::AbstractDictionary) = randn!(Random.GLOBAL_RNG, dict)

function Random.randn!(rng::AbstractRNG, dict::AbstractDictionary)
    map!(() -> randn(rng, eltype(dict)), dict)
    return dict
end

Random.randn(dict::AbstractDictionary) = randn(Random.GLOBAL_RNG, Float64, dict)
Random.randn(::Type{T}, dict::AbstractDictionary) where {T} = randn(Random.GLOBAL_RNG, T, dict)
Random.randn(rng::AbstractRNG, dict::AbstractDictionary) = randn(rng, Float64, dict)

function Random.randn(rng::AbstractRNG, ::Type{T}, dict::AbstractDictionary) where {T}
    out = similar(dict, T)
    randn!(rng, out)
    return out
end

# Copying - note that this doesn't necessarily copy the indices! (`copy(keys(dict))` can do that)
"""
    copy(dict::AbstractDictionary)
    copy(dict::AbstractDictionary, ::Type{T})

Create a shallow copy of the values of `dict`. Note that `keys(dict)` is not copied, and
therefore care must be taken that inserting/deleting elements. A new element type `T` can
optionally be specified.
"""
Base.copy(dict::AbstractDictionary) = copy(dict, eltype(dict))
function Base.copy(d::AbstractDictionary, ::Type{T}) where {T}
    out = similar(d, T)
    copyto!(out, d)
    return out
end

function Base.copyto!(out::AbstractDictionary, d::AbstractDictionary)
    map!(identity, out, d)
end

Base.last(d::AbstractDictionary) = first(Iterators.reverse(d))
Base.firstindex(d::AbstractDictionary) = first(keys(d))
Base.lastindex(d::AbstractDictionary) = first(Iterators.reverse(keys(d)))


text:
@testset "Data" begin

@test artificialIn_SAR_image(2) == 2*π*ones(2,2)

@test artificial_S1_slope_signal(20,0.) == repeat([-π/2],20)

@test ismissing( artificial_S1_signal(-1.) )
@test ismissing( artificial_S1_signal(2.) )
@test artificial_S1_signal(2) == [-3*π/4,-3*π/4]

# for the remainder check data types only
@test length( artificial_S1_signal(20) ) == 20

@test size( artificial_S2_whirl_image(64) ) == (64,64)
@test length(artificial_S2_whirl_image(64)[1,1]) == 3

@test size( artificial_S2_rotation_image(64) ) == (64,64)
@test length(artificial_S2_rotation_image(64)[1,1]) == 3

@test size( artificial_S2_whirl_patch(8) ) == (8,8)
@test length(artificial_S2_whirl_patch(8)[1,1]) == 3

@test size( artificial_SPD_image(8) ) == (8,8)
@test size(artificial_SPD_image(8)[1,1]) == (3,3)

@test size( artificial_SPD_image2(8) ) == (8,8)
@test size(artificial_SPD_image2(8)[1,1]) == (3,3)
@test eltype( artificial_SPD_image2(8))  == Array{Float64,2}

@test length( artificial_S2_lemniscate([0.,0.,1.], 20))  == 20
@test length( artificial_S2_lemniscate([0.,0.,1.], 20)[1]) == 3

end

ext: jl


text:
module Y
    f(x) = x + 1
    g(x) = x + 2

    export f
end

module Z
    f(x) = x + 10
    g(x) = x + 20
    h(x) = x + 30
    k(x) = x + 40
    export f, h, k
end

module X

using Reexport

@reexport using Y
# all of Y's exported symbols available here
@reexport using Z: f, g
# Z's x and y symbols available here

f(x) = x + 100;
g(x) = x + 200;
h(x) = x + 300;

export f, g

end

using X


text:
#
# Manopt.jl – The manifold of symmetric positive definite matrices
# with affine metric.
#
#
# ---
# Manopt.jl - Ronny Bergmann - 2017-07-06

import Base: exp, log, show
import LinearAlgebra: norm, dot

using Base: eps
using LinearAlgebra: svd, Diagonal, diag, diagm, eigen, eigvals, eigvecs, tr, triu, qr, cholesky, Hermitian

export SymmetricPositiveDefinite, SPDPoint, SPDTVector, show
export distance, dot, exp, norm, manifoldDimension, parallelTransport
export validateMPoint, validateTVector
export randomMPoint, randomTVector
export typeofMPoint, typeofTVector
export validateMPoint, validateTVector
export zeroTVector
# Types
# ---
@doc doc"""
    SymmetricPositiveDefinite <: Manifold

The manifold $\mathcal M = \mathcal P(n)$ of $n\times n$ symmetric positive
definite matrices.

# Fields
* `name` – representative String representing the current manifold
* `n` – size of the matrices of this manifold, i.e. $n\times n$ matrices
* `abbreviation` – short descriptor for the manifold.

# Constructor
    SymmetricPositiveDefinite(n)

construct the manifold of `n`-by-`n` matrices with affine metric.
"""
struct SymmetricPositiveDefinite <: Manifold
  name::String
  n::Int
  abbreviation::String
  SymmetricPositiveDefinite(n::Int) = new(
      "$n-by-$n symmetric positive definite matrices",n,"SPD($n) affine")
end
@doc doc"""
    SPDPoint <: MPoint

A point $x$ on the [`SymmetricPositiveDefinite`](@ref) manifold
$\mathcal M = \mathcal P(n)$ of $n\times n$,
represented in the redundant way of a symmetric positive definite matrix.
"""
struct SPDPoint{ T <: AbstractFloat} <: MPoint
	value::Matrix{T}
	SPDPoint{T}(v::Matrix{T}) where T<:AbstractFloat = new(v);
end
SPDPoint(v::Matrix{T}) where {T <: AbstractFloat} = SPDPoint{T}(v)
getValue(x::SPDPoint) = x.value
@doc doc"""
    SPDTVector <: TVector

A tangent vector

$\xi \in T_x\mathcal M = \{ x^{\frac{1}{2}}\nu x^{\frac{1}{2}}
\big| \nu\in\mathbb R^{n,n}\text{ with }\nu=\nu^{\mathrm{T}}\}$

to the [`SymmetricPositiveDefinite`](@ref) manifold
$\mathcal M = \mathcal P(n)$ at the [`SPDPoint`](@ref) `x` represented in the
redundant way of a skew symmetric matrix $\nu$, i.e. in the Lie algebra
$T_I\mathcal P(n)$, where $I\in\mathbb R^{n\times n}$ denotes the identity
matrix.
"""
struct SPDTVector{T <: AbstractFloat} <: TVector
	value::Matrix{T}
  	SPDTVector{T}(value::Matrix{T}) where {T <: AbstractFloat} = new(value);
end
SPDTVector(value::Matrix{T}) where {T<:AbstractFloat} = SPDTVector{T}(value)
getValue(ξ::SPDTVector) = ξ.value
# Traits
# ---
# (a) P(n) is a matrix manidolf
@traitimpl IsMatrixM{SymmetricPositiveDefinite}
@traitimpl IsMatrixP{SPDPoint}
@traitimpl IsMatrixTV{SPDTVector}
# Functions
# ---
@doc doc"""
    distance(M,x,y)

compute the Riemannian distance on the [`SymmetricPositiveDefinite`](@ref)
manifold $\mathcal M=\mathcal P(n)$, given by

```math
d_{\mathcal P(n)}(x,y) = \lVert \operatorname{Log}(x^{-\frac{1}{2}}yx^{-\frac{1}{2}})\rVert,
```

where $\operatorname{Log}$ denotes the matrix logarithm and the Norm is the Frobenius norm
in matrices
"""
function distance(M::SymmetricPositiveDefinite,x::SPDPoint,y::SPDPoint)
    S = log(getValue(x)\getValue(y))
    return real(sqrt(tr(transpose(S)*S)))
end
@doc doc"""
    dot(M,x,ξ,ν)

compute the innter product of the two [`SPDTVector`](@ref)s `ξ,ν` from the tangent
space $T_x\mathcal M$ of the [`SPDPoint`](@ref) `x` on the
[`SymmetricPositiveDefinite`](@ref) manifold `M` given by the formula

$ \langle \xi, \nu \rangle_x = \operatorname{tr}(x^{-1}\xi x^{-1}\nu ),$

where $\operatorname{tr}(y)$ denotes the trace of the matrix $y$.
"""
dot(M::SymmetricPositiveDefinite, x::SPDPoint, ξ::SPDTVector, ν::SPDTVector) = tr( (getValue(x)\getValue(ξ)) * (getValue(x)\getValue(ν)) )

@doc doc"""
    exp(M,x,ξ,[t=1.0])

compute the exponential map on the [`SymmetricPositiveDefinite`](@ref) manifold `M`
 `M`$=\mathcal P(n)$ with respect to the [`SPDPoint`](@ref) `x` and the
[`SPDTVector`](@ref)` ξ`, which can be shortened with `t` to `tξ`.
The formula reads

$\exp_x\xi = x^{\frac{1}{2}}\operatorname{Exp}(x^{-\frac{1}{2}}\xi x^{-\frac{1}{2}})x^{\frac{1}{2}},$

where $\operatorname{Exp}$ denotes the matrix exponential
"""
function exp(M::SymmetricPositiveDefinite, x::SPDPoint, ξ::SPDTVector, t::Float64=1.0)
    S = getValue(x)*real(exp(getValue(x)\(t*getValue(ξ))))
    return SPDPoint( 0.5*(transpose(S)+S) )
end
@doc doc"""
    log(M,x,y)

compute the logarithmic map on the [`SymmetricPositiveDefinite`](@ref) manifold
$\mathcal M=\mathcal P(n)$, i.e. the [`SPDTVector`](@ref) whose corresponding
[`geodesic`](@ref) starting from [`SPDPoint`](@ref) `x` reaches the
[`SPDPoint`](@ref)`y` after time 1.
The formula reads

$\log_x y = x^{\frac{1}{2}}\operatorname{Log}(x^{-\frac{1}{2}} y x^{-\frac{1}{2}})x^{\frac{1}{2}},$

where $\operatorname{Log}$ denotes the matrix logarithm.
"""#
function log(M::SymmetricPositiveDefinite,x::SPDPoint,y::SPDPoint)
    S = getValue(x)*real(log(getValue(x)\getValue(y)))
    return SPDTVector( 0.5*(transpose(S) + S) )
end
@doc doc"""
    manifoldDimension(M)

return the manifold dimension of the [`SymmetricPositiveDefinite`](@ref)
[`Manifold`](@ref) `M`, i.e. for $n\times n$ matrices the dimension
is $\frac{n(n+1)}{2}$.
"""
manifoldDimension(M::SymmetricPositiveDefinite) = (M.n+1)*M.n/2
@doc doc"""
    manifoldDimension(x)

returns the manifold dimension of the [`SymmetricPositiveDefinite`](@ref) manifold `M`
manifold the [`SPDPoint`](@ref) `x` belongs to,
i.e. for $n\times n$ matrices the dimension is
$\frac{n(n+1)}{2}$.
"""
manifoldDimension(x::SPDPoint) = Int( size( getValue(x), 1)*(size( getValue(x), 1)+1)/2 )
@doc doc"""
    norm(M,x,ξ)

compute the norm of the [`SPDTVector`](@ref)` ξ` from the tangent space $T_x\mathcal M$
at the [`SPDPoint`](@ref) `x` on the [`SymmetricPositiveDefinite`](@ref) manifold `M`
induced by the inner product [`dot`](@ref) as $\lVert\xi\rVert_x = \sqrt{\langle\xi,\xi\rangle_x}$.
"""
function norm(M::SymmetricPositiveDefinite,x::SPDPoint,ξ::SPDTVector)
    S = getValue(x)\getValue(ξ)
    return real(sqrt(tr( transpose(S)*S )))
end
@doc doc"""
    parallelTransport(M,x,y,ξ)

compute the paralllel transport of the [`SPDTVector`](@ref)` ξ` from
the tangent space $T_x\mathcal M$ at [`SPDPoint`](@ref) `x` to
$T_y\mathcal M$ at [`SPDPoint`](@ref)`y` on the [`SymmetricPositiveDefinite`](@ref) manifold `M`
along the [`geodesic`](@ref) $g(\cdot;x,y)$.
The formula reads

$P_{x\to y}(\xi) = x^{\frac{1}{2}}
\operatorname{Exp}\bigl(
\frac{1}{2}x^{-\frac{1}{2}}\log_x(y)x^{-\frac{1}{2}}
\bigr)
x^{-\frac{1}{2}}\xi x^{-\frac{1}{2}}
\operatorname{Exp}\bigl(
\frac{1}{2}x^{-\frac{1}{2}}\log_x(y)x^{-\frac{1}{2}}
\bigr)
x^{\frac{1}{2}},$

where $\operatorname{Exp}$ denotes the matrix exponential
and `log` the logarithmic map.
"""
function parallelTransport(M::SymmetricPositiveDefinite,x::SPDPoint,y::SPDPoint,ξ::SPDTVector)
  svd1 = svd( getValue(x) )
  if norm(getValue(x)-getValue(y))<1e-13
    return ξ
  end
  U = svd1.U
  S = svd1.S
  Ssqrt = sqrt.(S)
  SsqrtInv = Matrix(  Diagonal( 1 ./ Ssqrt )  )
  Ssqrt = Matrix(  Diagonal( Ssqrt )  )
  xSqrt = U*Ssqrt*transpose(U)
  xSqrtInv = U*SsqrtInv*transpose(U)
  tξ = xSqrtInv * getValue(ξ) * xSqrtInv
  tY = xSqrtInv * getValue(y) * xSqrtInv
  svd2 = svd(tY)
  Se = Matrix(  Diagonal( log.(svd2.S) )  )
  Ue = svd2.U
  tY2 = Ue*Se*transpose(Ue)
  eig1 = eigen(0.5*tY2)
  Sf = Matrix(  Diagonal( exp.(eig1.values) )  )
  Uf = eig1.vectors
  ν = xSqrt*Uf*Sf*transpose(Uf)*(0.5*(tξ+transpose(tξ)))*Uf*Sf*transpose(Uf)*xSqrt
  return SPDTVector(ν)
end

@doc doc"""
    randomMPoint(M,:Gaussian[, σ=1.0])

gerenate a random symmetric positive definite matrix on the
[`SymmetricPositiveDefinite`](@ref) manifold `M`.
"""
function randomMPoint(M::SymmetricPositiveDefinite,::Val{:Gaussian},σ::Float64=1.0)
    D = Diagonal( 1 .+ randn(M.n) ) # random diagonal matrix
    s = qr(σ * randn(M.n,M.n)) # random q
    return SPDPoint(s.Q*D* transpose(s.Q))
end
@doc doc"""
    randomTVector(M,x [,:Gaussian,σ = 1.0])

generate a random tangent vector in the tangent space of the [`SPDPoint`](@ref)
`x` on the [`SymmetricPositiveDefinite`](@ref) manifold `M` by using a Gaussian distribution
with standard deviation `σ` on an ONB of the tangent space.
"""
function randomTVector(M::SymmetricPositiveDefinite, x::SPDPoint, ::Val{:Gaussian}, σ::Float64 = 0.01)
    # generate ONB in TxM
    I = SPDPoint(one(zeros(M.n,M.n)))
    Ξ,κ = tangentONB(M,I,SPDTVector(one(zeros(M.n,M.n))))
    Ξx = parallelTransport.(Ref(M), Ref(I), Ref(x), Ξ)
    return sum( randn(length(Ξx)) .* Ξx )
end
@doc doc"""
    randomTVector(M,x [,:Gaussian,σ = 1.0])

generate a random tangent vector in the tangent space of the [`SPDPoint`](@ref)
`x` on the [`SymmetricPositiveDefinite`](@ref) manifold `M` by using a Rician distribution
with standard deviation `σ`.
"""
function randomTVector(M::SymmetricPositiveDefinite, x::SPDPoint, ::Val{:Rician}, σ::Real = 0.01)
    # Rician
    C = cholesky( Hermitian(getValue(x)) )
    R = sqrt(σ) * triu( randn(M.n,M.n),0)
    T = C.L * transpose(R)*R*C.U
    return log(M,x, SPDPoint(T))
end

@doc doc"""
    Ξ,κ = tangentONB(M,x,y)

compute a ONB in the tangent space of the [`SPDPoint`](@ref) `x` on the
[`SymmetricPositiveDefinite`](@ref) manifold `M` where the first vector is given by the
normed `log(M,x,y)`, i.e. the direction to the [`SPDPoint`](@ref) `y`.
"""
tangentONB(M::SymmetricPositiveDefinite, x::SPDPoint, y::SPDPoint) = tangentONB(M,x,log(M,x,y))
@doc doc"""
    Ξ,κ = tangentONB(M,x,y)

compute a ONB in the tangent space of the [`SPDPoint`](@ref) `x` on the
[`SymmetricPositiveDefinite`](@ref) manifold `M` where the first
vector is the normed tangent vector of the [`SPDTVector`](@ref) `ξ`.

The basis is computed using the eigenvectors $v_i$, $i=1,\ldots,n$, of `ξ` and
define

```math
\xi_{ij} := \begin{cases}
    \frac{1}{2} (v_i v_j^{\mathrm{T}} + v_j v_i^{\mathrm{T}})
    & \mathrm{ if } i=j	,\\
    \frac{1}{\sqrt{2}} (v_i v_j^{\mathrm{T}} + v_j v_i^{\mathrm{T}}),
    & \mathrm{ if } i \neq j,
\end{cases} \qquad\text{where } i=1,\ldots,n, k=i,\ldots,n,
```

and the correspoinding eigenvalues of the curvature tensor are given using the
eigenvalues $\lambda_i$, $i=1,\ldots,n$ of `ξ` as

```math
\kappa_{i,j} = -\frac{1}{4}(\lambda_i-\lambda_j)^2,\qquad
\text{where } i=1,\ldots,n, k=i,\ldots,n,
```
"""
function tangentONB(M::SymmetricPositiveDefinite,x::SPDPoint,ξ::SPDTVector)
    xSqrt = sqrt(getValue(x))
    v = eigvecs(getValue(ξ))
    Ξ = [
        SPDTVector( (i==j ? 1/2 : 1/sqrt(2))*
            ( v[:,i] * transpose(v[:,j])  +  v[:,j] * transpose(v[:,i]) )
        )
        for i=1:M.n for j= i:M.n
    ]
    λ = eigvals(getValue(ξ))
    κ = [ -1/4 * (λ[i]-λ[j])^2 for i=1:M.n for j= i:M.n ]
  return Ξ,κ
end

typeofTVector(::Type{SPDPoint{T}}) where T = SPDTVector{T}
typeofMPoint(::Type{SPDTVector{T}}) where T = SPDPoint{T}

@doc doc"""
    typicalDistance(M)

returns the typical distance on the [`SymmetricPositiveDefinite`](@ref)
manifold `M` $\sqrt{\frac{n(n+1)}{2}}$.
"""
typicalDistance(M::SymmetricPositiveDefinite) = sqrt(manifoldDimension(M));

@doc doc"""
    validateMPoint(M,x)

validate that the [`SPDPoint`](@ref) `x` is a valid point on the manifold
[`SymmetricPositiveDefinite`](@ref) manifold `M`, i.e. the matrix is symmetric and
positive definite.
"""
function validateMPoint(M::SymmetricPositiveDefinite, x::SPDPoint)
      if manifoldDimension(M) ≠ manifoldDimension(x)
        throw(DomainError(
            "The point $x does not lie on $M,, since the manifold dimension of $M ($(manifoldDimension(M)))does not fit the manifold dimension of $x ($(manifoldDimension(x)))."
        ))
    end
    if norm(getValue(x) - transpose(getValue(x))) > 10^(-14)
        throw(DomainError(
            "The point $x does not lie on $M, since the matrix of $x is not symmetric."
        ))
    end
    e = eigvals(getValue(x))
    if  !( (all(isreal.(e))) && (all( e .> 0)) ) # not spd
        throw(DomainError(
            "The point $x does not lie on $M, since the matrix of $x is not symmetric positive definite."
        ))
    end
    return true
end

@doc doc"""
    validateTVector(M,x,ξ)

validate, that the [`SPDTVector`](@ref)` ξ` is a tangent vector at the
[`SPDPoint`](@ref) `x` on the [`SymmetricPositiveDefinite`](@ref) `M`,
i.e. all dimensions are corrrect and the matrix is skew symmetric since
we only store the corresponding value in the Lie algebra.
"""
function validateTVector(M::SymmetricPositiveDefinite,x::SPDPoint,ξ::SPDTVector)
    if !validateMPoint(M,x)
        return false
    end
    if size(getValue(x)) != size(getValue(ξ))
        throw(DomainError(
            "The tangent vector $ξ can not be a tangent vector to $x (on $M), since the size of the matrix of ξ ($(size(getValue(ξ))) does not match the size of its base point matrix ($(size(getValue(x))))."
        ))
    end
    if norm(getValue(ξ) - transpose(getValue(ξ))) > 10^(-14)
        throw(DomainError(
            "The tangent vector $ξ does not represent a tangent vector (to $x) on $M, since the matrix of $ξ is not symmetric."
        ))
    end
    return true
end

@doc doc"""
    ξ = zeroTVector(M,x)

returns a zero vector in the tangent space $T_x\mathcal M$ of the
[`SPDPoint`](@ref) $x\in\mathcal P(n)$ on the
[`SymmetricPositiveDefinite`](@ref) manifold `M`.
"""
zeroTVector(M::SymmetricPositiveDefinite, x::SPDPoint) = SPDTVector(  zero( getValue(x) )  );
# Display
# ---
show(io::IO, M::SymmetricPositiveDefinite) = print(io, "The Manifold of $(M.name)")
show(io::IO, x::SPDPoint) = print(io, "SPD($(x.value))")
show(io::IO, ξ::SPDTVector) = print(io, "SPDT($(ξ.value))")


text:
@inline function ExaTronKernels.dtrpcg(n::Int, A::oneDeviceArray{Float64,2},
                        g::oneDeviceArray{Float64,1}, delta::Float64,
                        L::oneDeviceArray{Float64,2},
                        tol::Float64, stol::Float64, itermax::Int,
                        w::oneDeviceArray{Float64,1},
                        p::oneDeviceArray{Float64,1},
                        q::oneDeviceArray{Float64,1},
                        r::oneDeviceArray{Float64,1},
                        t::oneDeviceArray{Float64,1},
                        z::oneDeviceArray{Float64,1})
  tx = get_local_id()
  zero = 0.0
  one = 1.0

  # Initialize the iterate w and the residual r.
  fill!(w, 0.0)

  # Initialize the residual t of grad q to -g.
  # Initialize the residual r of grad Q by solving L*r = -g.
  # Note that t = L*r.
  dcopy(n,g,1,t,1)
  dscal(n,-one,t,1)
  dcopy(n,t,1,r,1)
  dnsol(n, L, r)

  # Initialize the direction p.
  dcopy(n,r,1,p,1)

  # Initialize rho and the norms of r and t.
  rho = ddot(n,r,1,r,1)
  rnorm0 = sqrt(rho)

  # Exit if g = 0.
  iters = 0
  if rnorm0 == zero
    iters = 0
    info = 1
    return info, iters
  end

  for iters=1:itermax

    # Note:
    # Q(w) = 0.5*w'Bw + h'w, where B=L^{-1}AL^{-T}, h=L^{-1}g.
    # Then p'Bp = p'L^{-1}AL^{-T}p = p'L^{-1}Az = p'q.
    # alpha = r'r / p'Bp.

    dcopy(n,p,1,z,1)
    dtsol(n, L, z)

    # Compute q by solving L*q = A*z and save L*q for
    # use in updating the residual t.
    dssyax(n, A, z, q)
    dcopy(n,q,1,z,1)
    dnsol(n, L, q)

    # Compute alpha and determine sigma such that the trust region
    # constraint || w + sigma*p || = delta is satisfied.
    ptq = ddot(n,p,1,q,1)
    if ptq > zero
      alpha = rho/ptq
    else
      alpha = zero
    end
    sigma = dtrqsol(n,w,p,delta)

    # Exit if there is negative curvature or if the
    # iterates exit the trust region.

    if (ptq <= zero) || (alpha >= sigma)
      daxpy(n,sigma,p,1,w,1)
      if ptq <= zero
        info = 3
      else
        info = 4
      end

      return info, iters
    end

    # Update w and the residuals r and t.
    # Note that t = L*r.

    daxpy(n,alpha,p,1,w,1)
    daxpy(n,-alpha,q,1,r,1)
    daxpy(n,-alpha,z,1,t,1)

    # Exit if the residual convergence test is satisfied.

    rtr = ddot(n,r,1,r,1)
    rnorm = sqrt(rtr)
    tnorm = sqrt(ddot(n,t,1,t,1))

    if tnorm <= tol
      info = 1
      return info, iters
    end

    if rnorm <= stol
      info = 2
      return info, iters
    end

    # Compute p = r + beta*p and update rho.
    beta = rtr/rho
    dscal(n,beta,p,1)
    daxpy(n,one,r,1,p,1)
    rho = rtr
  end

  iters = itermax
  info = 5
  return info, iters
end


text:
# ---
# title: 1414. Find the Minimum Number of Fibonacci Numbers Whose Sum Is K
# id: problem1414
# author: Indigo
# date: 2022-2-18
# difficulty: Medium
# categories: Array, Greedy
# link: <https://leetcode.com/problems/find-the-minimum-number-of-fibonacci-numbers-whose-sum-is-k/description/>
# hidden: true
# ---
#
# Given an integer `k`, _return the minimum number of Fibonacci numbers whose
# sum is equal to_`k`. The same Fibonacci number can be used multiple times.
#
# The Fibonacci numbers are defined as:
#
#   * `F1 = 1`
#   * `F2 = 1`
#   * `Fn = Fn-1 + Fn-2` for `n > 2.`
#
# It is guaranteed that for the given constraints we can always find such
# Fibonacci numbers that sum up to `k`.
#
#
#
# **Example 1:**
#
#
#
#     Input: k = 7
#     Output: 2
#     Explanation: The Fibonacci numbers are: 1, 1, 2, 3, 5, 8, 13, ...
#     For k = 7 we can use 2 + 5 = 7.
#
# **Example 2:**
#
#
#
#     Input: k = 10
#     Output: 2
#     Explanation: For k = 10 we can use 2 + 8 = 10.
#
#
# **Example 3:**
#
#
#
#     Input: k = 19
#     Output: 3
#     Explanation: For k = 19 we can use 1 + 5 + 13 = 19.
#
#
#
#
# **Constraints:**
#
#   * `1 <= k <= 10^9`
#
#
## @lc code=start
using LeetCode

function find_min_fibonacci_numbers(k::Int)
    fibs = [0, 1, 1]
    while fibs[end] <= k
        push!(fibs, fibs[end] + fibs[end - 1])
    end
    """
    If we need 2 or more fib_i's,
    we can always use fib_{i+1} and fib_{i-2} to replace them:
    2fib_i = fib_i + fib_{i-1} + fib_{i-2} = fib_{i+1} + fib_{i-2}.

    So greedy algorithm can be performed.
    """
    pop!(fibs)
    res = 0
    idx = length(fibs) + 1
    while k != 0
        idx -= 1
        k >= fibs[idx] || continue
        k -= fibs[idx]
        res += 1
    end
    return res
end
## @lc code=end


text:
# Preprocessing and postprocessing helpers

export normalize!, one_hot, inv_one_hot

function normalize!(x::AbstractArray{Float64})
    (d, n) = size(x)
    @assert n > 1
    μ = sum(x, dims=1) ./ n
    x .-= μ
    σ = @. sqrt($sum(x ^ 2, dims=1) / (n - 1))
    x ./= σ
end

function one_hot(y::AbstractVector{Int}, n_classes::Int)::Matrix{Float64}
    n = length(y)
    one_hot_y = zeros(n_classes, n)
    # Creates one-hot vectors. Better to make it an `AbstractVector`.
    one_hot_y[CartesianIndex.(y, 1:n)] .= 1.0
    one_hot_y
end

function inv_one_hot(y::AbstractMatrix{Float64})::Vector{Int}
    vec(mapslices(argmax, y, dims=1))
end

max_stars_repo_name: gyk/TrivialSolutions

max_issues_repo_name: gyk/TrivialSolutions

max_forks_repo_name: gyk/TrivialSolutions


text:
@testset "Scores; Guess is correct; User is given a point" begin
    # Arrange
    gameservice = TestingGameService()
    dictionary = SwedishDictionary(["DATORSPEL"])
    game = NiancatGame(dictionary, gameservice)

    team = Team(1, "defaultteam", "")
    user = User(1, "name", "displayname", team)
    gamecommand(game, user, SetPuzzle("DATORSPLE"))

    # Act
    response = gamecommand(game, user, Guess("DATORSPEL"))

    # Assert
    @test length(gameservice.scores) == 1
    (actualuser, score) = gameservice.scores[1]
    @test user == actualuser
    @test score.value == 1
end

@testset "Scores; Guess is wrong; User is not given a point" begin
    # Arrange
    gameservice = TestingGameService()
    dictionary = SwedishDictionary(["DATORSPEL"])
    game = NiancatGame(dictionary, gameservice)

    team = Team(1, "defaultteam", "")
    user = User(1, "name", "displayname", team)
    gamecommand(game, user, SetPuzzle("DATORSPLE"))

    # Act
    response = gamecommand(game, user, Guess("DATORLESP"))

    # Assert
    @test length(gameservice.scores) == 0
end

@testset "Scores; Guess same thing twice; Both scores events have the same key" begin
    # Arrange
    gameservice = TestingGameService()
    dictionary = SwedishDictionary(["DATORSPEL"])
    game = NiancatGame(dictionary, gameservice)

    team = Team(1, "defaultteam", "")
    user = User(1, "name", "displayname", team)
    gamecommand(game, user, SetPuzzle("DATORSPLE"))

    # Act
    gamecommand(game, user, Guess("DATORSPEL"))
    gamecommand(game, user, Guess("DATORSPEL"))

    # Assert
    @test length(gameservice.scores) == 2
    (_user1, score1) = gameservice.scores[1]
    (_user2, score2) = gameservice.scores[2]

    @test score1.key == score2.key
end

@testset "Scores; Guess two different solutions; The score events have different keys" begin
    # Arrange
    gameservice = TestingGameService()
    dictionary = SwedishDictionary(["DATORSPEL", "LEDARPOST"])
    game = NiancatGame(dictionary, gameservice)

    team = Team(1, "defaultteam", "")
    user = User(1, "name", "displayname", team)
    gamecommand(game, user, SetPuzzle("DATORSPLE"))

    # Act
    gamecommand(game, user, Guess("DATORSPEL"))
    gamecommand(game, user, Guess("LEDARPOST"))

    # Assert
    @test length(gameservice.scores) == 2
    (_user1, score1) = gameservice.scores[1]
    (_user2, score2) = gameservice.scores[2]

    @test score1.key != score2.key
end

@testset "Scores; Guess same thing twice with non-alphabetic differences; Both scores events have the same key" begin
    # Arrange
    gameservice = TestingGameService()
    dictionary = SwedishDictionary(["DATORSPEL"])
    game = NiancatGame(dictionary, gameservice)

    team = Team(1, "defaultteam", "")
    user = User(1, "name", "displayname", team)
    gamecommand(game, user, SetPuzzle("DATORSPLE"))

    # Act
    gamecommand(game, user, Guess("DATOR-SPEL"))
    gamecommand(game, user, Guess("DATORSPEL"))

    # Assert
    @test length(gameservice.scores) == 2
    (_user1, score1) = gameservice.scores[1]
    (_user2, score2) = gameservice.scores[2]

    @test score1.key == score2.key
end

@testset "Scores; Guess two different solutions; The score events have the same round" begin
    # Arrange
    gameservice = TestingGameService()
    dictionary = SwedishDictionary(["DATORSPEL", "LEDARPOST"])
    game = NiancatGame(dictionary, gameservice)

    team = Team(1, "defaultteam", "")
    user = User(1, "name", "displayname", team)
    gamecommand(game, user, SetPuzzle("DATORSPLE"))

    # Act
    gamecommand(game, user, Guess("DATORSPEL"))
    gamecommand(game, user, Guess("LEDARPOST"))

    # Assert
    @test length(gameservice.scores) == 2
    (_user1, score1) = gameservice.scores[1]
    (_user2, score2) = gameservice.scores[2]

    @test score1.round == score2.round
end

@testset "Scores; Score points in two different rounds; The score events have different rounds" begin
    # Arrange
    gameservice = TestingGameService()
    dictionary = SwedishDictionary(["DATORSPEL", "LEDARPOST", "ORDPUSSEL"])
    game = NiancatGame(dictionary, gameservice)

    team = Team(1, "defaultteam", "")
    user = User(1, "name", "displayname", team)

    # Act
    gamecommand(game, user, SetPuzzle("DATORSPLE"))
    gamecommand(game, user, Guess("DATORSPEL"))

    gamecommand(game, user, SetPuzzle("ORDPUSSLE"))
    gamecommand(game, user, Guess("ORDPUSSEL"))

    # Assert
    @test length(gameservice.scores) == 2
    (_user1, score1) = gameservice.scores[1]
    (_user2, score2) = gameservice.scores[2]

    @test score1.round != score2.round
end


text:
using Clang.cindex
using Clang.wrap_c
using Compat

XCB_INCLUDE = "/usr/include/xcb"

xcb_header = [joinpath(XCB_INCLUDE, "xcb.h")]

# Set up include paths
clang_includes = ASCIIString[]
push!(clang_includes, XCB_INCLUDE)

# Clang arguments
clang_extraargs = ["-v"]
# clang_extraargs = ["-D", "__STDC_LIMIT_MACROS", "-D", "__STDC_CONSTANT_MACROS"]

# Callback to test if a header should actually be wrapped (for exclusion)
function wrap_header(top_hdr::ASCIIString, cursor_header::ASCIIString)
    return startswith(dirname(cursor_header), XCB_INCLUDE)
end

lib_file(hdr::ASCIIString) = "xcb"
output_file(hdr::ASCIIString) = Pkg.dir("XCB", "src", "c_interface.jl")

function wrap_cursor(name::ASCIIString, cursor)
    return true
end

const wc = wrap_c.init(;
                        headers = xcb_header,
                        output_file = Pkg.dir("XCB", "src", "c_interface.jl"),
                        common_file = Pkg.dir("XCB", "src", "c_types.jl"),
                        clang_includes      = clang_includes,
                        clang_args          = clang_extraargs,
                        header_wrapped      = wrap_header,
                        header_library      = lib_file,
                        header_outputfile   = output_file,
                        cursor_wrapped      = wrap_cursor,
                        clang_diagnostics = true)

run(wc)


text:
__precompile__(true)

module Arrow

using CategoricalArrays, Compat, Compat.Dates

const BITMASK = UInt8[1, 2, 4, 8, 16, 32, 64, 128]
const ALIGNMENT = 8

import Base: getindex, setindex!
import Base: convert, show, unsafe_string, checkbounds, write, values, copy
import Base: length, endof, size, eltype, start, next, done, getindex, isassigned, view
import Base: IndexStyle
import Base: >, ≥, <, ≤, ==
import CategoricalArrays: levels

if VERSION ≥ v"0.7.0-"
else
    using Missings
    import Base: isnull
end

abstract type ArrowVector{T} <: AbstractVector{T} end
export ArrowVector

include("utils.jl")
include("primitives.jl")
include("lists.jl")
include("arrowvectors.jl")
include("datetime.jl")
include("dictencoding.jl")
include("bitprimitives.jl")

#=~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
TODO
    -make primitives with args re-interpreted views (to what extent possible)
    -support Bools!!!
    -better support for converting arrow formats
    -get views working properly (only available in 0.7!)
    -clean up bounds checking so it's not redundant! right now it's all fucked up
    -support nested lists (think should work but haven't tested)
    -allow users to specify ordering of sub-buffers

    -are there any ways to get views of strings???
    -ensure that convert(Array{Union{T,Missing}}, A) is efficient in 0.7
    -investigate breaking due to new reinterpret in 0.7
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~=#

end  # module Arrow


text:
#We used this function to always keep the saved pseudos updated
function with_pseudos(f::Function)
    pp = config_path("pseudos.jld2")
    pseudos = ispath(pp) ? JLD2.load(pp)["pseudos"] :
              Dict{String,Dict{Symbol,Vector{String}}}()
    t = f(pseudos)
    JLD2.jldsave(config_path("pseudos.jld2"); pseudos=pseudos)
    return t
end

function pseudos(set::String, fuzzy::String)
    with_pseudos() do all_sets
        out = Dict{Symbol,String}()
        pseudos = get(all_sets, set, nothing)
        pseudos === nothing && return nothing
        @info length(pseudos)
        for (k, ps) in pseudos
            p = getfirst(x -> occursin(fuzzy, x), ps)
            if p !== nothing
                out[k] = p
            end
        end
        return out
    end
end

function rm_pseudos!(set::String)
    with_pseudos() do pseudos
        pop!(pseudos, set, nothing)
        return nothing
    end
end

"""
    configure_pseudoset(set_name::String, dir::String)

Reads the specified `dir` and sets up the pseudos for `set`.
"""
function configure_pseudoset(set_name::String, dir::String)
    with_pseudos() do pseudos
        files = readdir(dir)
        pseudos[set_name] = Dict{Symbol,Vector{String}}([el.symbol => String[]
                                                         for el in Structures.ELEMENTS])
        for pseudo_string in files
            element = Symbol(titlecase(String(split(split(pseudo_string, ".")[1], "_")[1])))
            if haskey(pseudos[set_name], element)
                push!(pseudos[set_name][element], joinpath(dir, pseudo_string))
            end
        end
        return length(files)
    end
end

"""
    pseudo_sets()

Lists all pseudo_sets.
"""
pseudo_sets() = with_pseudos(x -> keys(x))


text:
export PSDMatrix, jacobian

struct PSDMatrix{N,T,F} <: Function where {N<:Integer,T<:DataType}
    data::Vector{T}
    initial_params::F
end

function PSDMatrix(T::DataType, N::Integer, init::Function)
    PSDMatrix{N,T,typeof(init)}( T.(init()), init )
end
PSDMatrix(N::Integer, init::Function) = PSDMatrix(Float32,N,()->vec(init()))
PSDMatrix(N::Integer) = PSDMatrix(N,()->vec(Flux.glorot_uniform(N,N)))

Flux.trainable(P::PSDMatrix) = (P.data,)
function DiffEqFlux.initial_params(P::PSDMatrix{N,T}) where {T,N}
    T.(P.initial_params())
end
DiffEqFlux.paramlength(::PSDMatrix{N}) where {N} = N*N

(m::PSDMatrix)(::AbstractVecOrMat) = posdef(m.data)
(m::PSDMatrix)(::AbstractVecOrMat, θ) = posdef(θ)
function jacobian(m::PSDMatrix{N,T}, ::AbstractVecOrMat, ::Any=nothing) where {N,T}
    collect(zeros(T,N,N) for _=1:N)
end


text:
function _create_new_shunt(sub_data::Dict, pm_data::Dict)
    for (idx, value) in pm_data["shunt"]
        if value["shunt_bus"] == sub_data["shunt_bus"] && value["source_id"][1] == sub_data["source_id"][1]
            return false, idx
        end
    end
    return true
end

function _handle_bs_bounds(shunt::Dict{String, Any})
    bsmin, bsmax = 0.0, 0.0
    for (i,el) in shunt["REACTANCE GROUPS"]
        if el["REACTANCE"] > 0.0
            bsmax += el["UNITIES"]*el["REACTANCE"]
        elseif el["REACTANCE"] < 0.0
            bsmin += el["UNITIES"]*el["REACTANCE"]
        end
    end
    return bsmin, bsmax
end

function _handle_bs(shunt::Dict{String, Any}; type = "bus")
    bs = 0
    if type == "bus" && shunt["CONTROL MODE"] == 'C'
        bs += shunt["INITIAL REACTIVE INJECTION"]
    else
        for (i,el) in shunt["REACTANCE GROUPS"]
            if el["STATUS"] == 'L'
                bs += el["OPERATING UNITIES"]*el["REACTANCE"]
            end
        end
    end
    return bs
end

function _pwf2pm_DBSH_shunt!(pm_data::Dict, pwf_data::Dict, shunt::Dict; add_control_data::Bool=false)

    shunt_bus = shunt["TO BUS"] === nothing ? shunt["FROM BUS"] : shunt["EXTREMITY"]
    n = count(x -> x["shunt_bus"] == shunt_bus, values(pm_data["shunt"]))

    sub_data = Dict{String,Any}()

    sub_data["shunt_bus"] = shunt_bus

    sub_data["gs"] = 0.0
    sub_data["bs"] = _handle_bs(shunt)

    status = pwf_data["DBAR"]["$(sub_data["shunt_bus"])"]["STATUS"]
    if status == 'L'
        sub_data["status"] = 1
    elseif status == 'D'
        sub_data["status"] = 0
    end

    if add_control_data
        sub_data["control_data"] = Dict{String,Any}()

        sub_data["control_data"]["section"] = "DBSH"

        sub_data["control_data"]["shunt_type"] = shunt["CONTROL MODE"] == 'F' ? 1 : 2
        sub_data["control_data"]["shunt_control_type"] = shunt["CONTROL MODE"] == 'F' ? 1 : shunt["CONTROL MODE"] == 'D' ? 2 : 3

        sub_data["control_data"]["vmmin"] = shunt["MINIMUM VOLTAGE"]
        sub_data["control_data"]["vmmax"] = shunt["MAXIMUM VOLTAGE"]

        sub_data["control_data"]["controlled_bus"] = shunt["CONTROLLED BUS"]
        bs_bounds = _handle_bs_bounds(shunt)
        sub_data["control_data"]["bsmin"] = bs_bounds[1]
        sub_data["control_data"]["bsmax"] = bs_bounds[2]
        @assert sub_data["control_data"]["bsmin"] <= sub_data["control_data"]["bsmax"]
        sub_data["control_data"]["inclination"] = nothing
    end

    sub_data["source_id"] = ["switched shunt", sub_data["shunt_bus"], "0$(n+1)"]
    sub_data["index"] = length(pm_data["shunt"]) + 1

    if _create_new_shunt(sub_data, pm_data)[1]
        idx = string(sub_data["index"])
        pm_data["shunt"][idx] = sub_data
    else
        idx = _create_new_shunt(sub_data, pm_data)[2]
        pm_data["shunt"][idx]["gs"] += sub_data["gs"]
        pm_data["shunt"][idx]["bs"] += sub_data["bs"]
        if add_control_data
            pm_data["shunt"][idx]["control_data"]["bsmin"] += sub_data["control_data"]["bsmin"]
            pm_data["shunt"][idx]["control_data"]["bsmax"] += sub_data["control_data"]["bsmax"]
        end
    end
end

function _pwf2pm_DBSH_shunt!(pm_data::Dict, pwf_data::Dict; add_control_data::Bool=false)
    if haskey(pwf_data, "DBSH")

        for (i,shunt) in pwf_data["DBSH"]
            _pwf2pm_DBSH_shunt!(pm_data, pwf_data, shunt, add_control_data = add_control_data)
        end
    end
end

function _pwf2pm_DCER_shunt!(pm_data::Dict, pwf_data::Dict, shunt::Dict; add_control_data::Bool=false)
    n = count(x -> x["shunt_bus"] == shunt["BUS"], values(pm_data["shunt"]))

    sub_data = Dict{String,Any}()
    sub_data["shunt_bus"] = shunt["BUS"]
    sub_data["gs"]        = 0.0
    sub_data["bs"]        = shunt["REACTIVE GENERATION"]

    status   = pwf_data["DBAR"]["$(sub_data["shunt_bus"])"]["STATUS"]

    if status == 'L'
        sub_data["status"] = 1
    elseif status == 'D'
        sub_data["status"] = 0
    end

    if add_control_data
        sub_data["control_data"] = Dict{String,Any}()

        sub_data["control_data"]["section"] = "DCER"

        sub_data["control_data"]["shunt_type"]         = 2
        sub_data["control_data"]["shunt_control_type"] = 3

        sub_data["control_data"]["bsmin"] = shunt["MINIMUM REACTIVE GENERATION"]
        sub_data["control_data"]["bsmax"] = shunt["MAXIMUM REACTIVE GENERATION"]

        @assert sub_data["control_data"]["bsmin"] <= sub_data["control_data"]["bsmax"]

        ctrl_bus = pm_data["bus"]["$(shunt["CONTROLLED BUS"])"]
        sub_data["control_data"]["vmmin"] = ctrl_bus["vm"]
        sub_data["control_data"]["vmmax"] = ctrl_bus["vm"]
        sub_data["control_data"]["controlled_bus"] = shunt["CONTROLLED BUS"]
        sub_data["control_data"]["inclination"] = shunt["INCLINATION"]
    end

    sub_data["source_id"] = ["switched shunt", sub_data["shunt_bus"], "0$(n+1)"]
    sub_data["index"]     = length(pm_data["shunt"]) + 1

    if _create_new_shunt(sub_data, pm_data)[1]
        idx = string(sub_data["index"])
        pm_data["shunt"][idx] = sub_data
    else
        idx = _create_new_shunt(sub_data, pm_data)[2]
        pm_data["shunt"][idx]["gs"] += sub_data["gs"]
        pm_data["shunt"][idx]["bs"] += sub_data["bs"]
        if add_control_data
            pm_data["shunt"][idx]["control_data"]["bsmin"] += sub_data["control_data"]["bsmin"]
            pm_data["shunt"][idx]["control_data"]["bsmax"] += sub_data["control_data"]["bsmax"]
        end
    end
end

function _pwf2pm_DCER_shunt!(pm_data::Dict, pwf_data::Dict; add_control_data::Bool=false)
    # Assumption - the reactive generation is already considering the number of unities
    if haskey(pwf_data, "DCER")

        for (i,shunt) in pwf_data["DCER"]
            _pwf2pm_DCER_shunt!(pm_data, pwf_data, shunt, add_control_data = add_control_data)
        end
    end
end

function _pwf2pm_DBAR_shunt!(pm_data::Dict, pwf_data::Dict, bus::Dict; add_control_data::Bool=false)
    sub_data = Dict{String,Any}()

    sub_data["shunt_bus"] = bus["NUMBER"]
    sub_data["gs"] = 0.0
    sub_data["bs"] =  bus["TOTAL REACTIVE POWER"]

    if bus["STATUS"] == 'L'
        sub_data["status"] = 1
    elseif bus["STATUS"] == 'D'
        sub_data["status"] = 0
    end

    if add_control_data
        sub_data["control_data"] = Dict{String,Any}()

        sub_data["control_data"]["section"] = "DBAR"

        sub_data["control_data"]["shunt_type"] = 1
        sub_data["control_data"]["shunt_control_type"] = 1

        sub_data["control_data"]["bsmin"] = sub_data["bs"]
        sub_data["control_data"]["bsmax"] = sub_data["bs"]
        @assert sub_data["control_data"]["bsmin"] <= sub_data["control_data"]["bsmax"]

        sub_data["control_data"]["vmmin"] = bus["VOLTAGE"]
        sub_data["control_data"]["vmmax"] = bus["VOLTAGE"]
        sub_data["control_data"]["controlled_bus"] = bus["CONTROLLED BUS"]
        sub_data["control_data"]["inclination"] = nothing
    end

    n = count(x -> x["shunt_bus"] == sub_data["shunt_bus"], values(pm_data["shunt"]))
    sub_data["source_id"] = ["fixed shunt", sub_data["shunt_bus"], "0$(n+1)"]
    sub_data["index"] = length(pm_data["shunt"]) + 1

    idx = string(sub_data["index"])
    pm_data["shunt"][idx] = sub_data
end

function _pwf2pm_DBAR_shunt!(pm_data::Dict, pwf_data::Dict; add_control_data::Bool=false)
     # fixed shunts specified at DBAR
     fixed_shunt_bus = findall(x -> x["TOTAL REACTIVE POWER"] != 0.0, pwf_data["DBAR"])

     for bus in fixed_shunt_bus
        _pwf2pm_DBAR_shunt!(pm_data, pwf_data, pwf_data["DBAR"][bus], add_control_data = add_control_data)
     end
end

# Assumption - if there are more than one shunt for the same bus we sum their values into one shunt (source: Organon)
# CAUTION: this might be an Organon error
function _pwf2pm_shunt!(pm_data::Dict, pwf_data::Dict; add_control_data::Bool=false)
    pm_data["shunt"] = Dict{String, Any}()

    _pwf2pm_DBAR_shunt!(pm_data, pwf_data, add_control_data = add_control_data)

    _pwf2pm_DCER_shunt!(pm_data, pwf_data, add_control_data = add_control_data)

    _pwf2pm_DBSH_shunt!(pm_data, pwf_data, add_control_data = add_control_data)
end


text:
# --------------------------------------------------------------------------
# ACE.jl and SHIPs.jl: Julia implementation of the Atomic Cluster Expansion
# Copyright (c) 2019 Christoph Ortner <christophortner0@gmail.com>
# Licensed under ASL - see ASL.md for terms and conditions.
# --------------------------------------------------------------------------

using StaticArrays, ACE, BenchmarkTools, LinearAlgebra
using ACE.SphericalHarmonics
SH = ACE.SphericalHarmonics
using ACE: alloc_temp, alloc_temp_d, alloc_B, alloc_dB, evaluate!, evaluate_d!

@info("Spherical Harmonics Evaluation")

##

suite = BenchmarkGroup()

# (@benchmark runmany!($Y, $tmp, $basis, 100)) |> display
# @btime runmany!($Y, $tmp, $basis, 100)
#
# R = ACE.rand_sphere()
# @btime evaluate!($Y, $tmp, $basis, $R)

R = JVecF(0.4, 0.7, -0.9)

for L in [5, 10, 15]
   basis = SH.SHBasis(L)
   tmp = alloc_temp(basis)
   Y = alloc_B(basis)
   tmp_d = alloc_temp_d(basis)
   dY = alloc_dB(basis)
   suite["Complex SH($L) - evaluate!"] =
         (@benchmarkable evaluate!($Y, $tmp, $basis, $R))
   suite["Complex SH($L) - evaluate_d!"] =
         (@benchmarkable evaluate_d!($Y, $dY, $tmp_d, $basis, $R))

   rbasis = SH.RSHBasis(L)
   rtmp = alloc_temp(rbasis)
   rY = alloc_B(rbasis)
   rtmp_d = alloc_temp_d(rbasis)
   rdY = alloc_dB(rbasis)
   suite["Real SH($L) - evaluate!"] =
         (@benchmarkable evaluate!($rY, $rtmp, $rbasis, $R))
   suite["Real SH($L) - evaluate_d!"] =
         (@benchmarkable evaluate_d!($rY, $rdY, $rtmp_d, $rbasis, $R))
end

if !isdefined(Main, :globalsuite)
   tune!(suite)
   results = run(suite; verbose=true)
   display(results)
else
   globalsuite["ylm"] = suite
end


text:
module Weights

import JuLIP, IPFitting
using IPFitting: LsqDB, Dat
using JuLIP: forces, Atoms, positions

function hess_weights!(cn::String, db::LsqDB; h = :auto, wE = 0.0, rscal = 3, verbose=false)

   @error("`hess_weights!` needs to be rewritten!")

   # temporarily assume h is a floating point number
   # => this needs to be automated
   @assert (h isa AbstractFloat)

   # loop through the data_groups
   for (key, dg) in db.data_groups
      # if the data_group has the correct configname
      if cn == configname(key)
         # then loop through the dats
         for d in dg
            if forces(d) == nothing
               verbose && @warn("""hess_weights!: a training configuration does
                                  not contain energy and forces => ignore""")
               d.info["W"] = Dict( "F" => zeros(3*length(d)),
                                   "E" => 0.0,
                                   "V" => zeros(6) )
               continue
            end
            # get configuration info, R vectors
            at = Atoms(d)
            X = positions(at)
            R = [ JuLIP.project_min(at, x - X[1])  for x in X ]

            # if h == :auto => COMPUTE h <<<<<<<<

            # now fix the scaling for the force weights
            r = norm.(R)
            r3 = (ones(3) * r')[:]
            wF = r3.^rscal / h
            wF[1:3] = 0.0   # this sets the weight for atom 1 to zero
                            # it contains no information...
            # assign this weight back to `d::Dat`
            d.info["W"] = Dict( "E" => wE,
                                "F" => wF,
                                "V" => zeros(6) )
         end
      end
   end

   return nothing
end

end # module

# function hess_weights_hook!(w, d::Dat)
#    at = Atoms(d)
#    if energy(d) == nothing || forces(d) == nothing
#       warn("""hess_weights_hook!: a training configuration does not contain
#               energy and forces => ignore""")
#      return w
#   end
#    # compute R vectors
#    X = positions(at)
#    h = 0.01 # norm(X[1])
#    # if h < 1e-5 || h > 0.02
#    #    warn("unexpected location of X[1] in hess_weights_hook => ignore")
#    #    @show X[1], h
#    #    return w
#    # end
#
#    # give this energy a lot of weight to make sure we match the
#    # ground state (which we assume this is)
#    w[1] = 0.0
#
#    # now fix the scaling for the force weights
#    # X[1] *= 0
#    R = [ JuLIP.project_min(at, x - X[1])  for x in X ]
#    r = norm.(R)
#    r3 = (ones(3) * r')[:]
#    If = 2:(1+3*length(R))
#    w[If] .= w[If] .* (r3.^7) / h
#    w[2] = 0.0
#    return w
# end


text:
using PyPlot
using PyCall
using Latexify
using JLD2
using Base
export plot_1D_solutions, plot_2D_phase_diagram, transform_solutions
export _set_plotting_settings, _prepare_colorbar, _prettify_label

"Set global plotting settings"
function _set_plotting_settings()
    plt.style.use("default") #reset settings
    rcParams = PyPlot.PyDict(PyPlot.matplotlib."rcParams")
    rcParams["text.usetex"]       = true
    rcParams["font.family"]       = "sans-serif"
    rcParams["font.sans-serif"]   = ["Helvetica"]
    rcParams["legend.frameon"]    = false
    rcParams["legend.fontsize"]   = "x-large"
    rcParams["patch.linewidth"]   =  0.3
    rcParams["lines.linewidth"]   = 2
    rcParams["lines.markersize"]  = 7
    rcParams["grid.linestyle"]    =  "-"
    rcParams["grid.linewidth"]    =  0.75
    rcParams["axes.axisbelow"]    =  true
    rcParams["axes.linewidth"]    =  1.25
    rcParams["axes.labelsize"]    = "large"
    rcParams["xtick.labelsize"]   = "xx-large"
    rcParams["ytick.labelsize"]   = "xx-large"
    rcParams["xtick.direction"]   = "in"
    rcParams["ytick.direction"]   = "in"
    rcParams["xtick.major.size"]  = 6
    rcParams["ytick.major.size"]  = 6
    rcParams["xtick.minor.size"]  = 3
    rcParams["ytick.minor.size"]  = 3
    rcParams["xtick.major.width"] = 1
    rcParams["ytick.major.width"] = 1
    rcParams["xtick.minor.width"] = .5
    rcParams["ytick.minor.width"] = .5
    rcParams["xtick.major.pad"]   = 7
    rcParams["ytick.major.pad"]   = 7
    rcParams["figure.autolayout"] = true
    rcParams["figure.dpi"] = 220
end

# the total number of solutions
_total_elements(array) = prod(size(array))

_add_dim!(x::Array) = reshape(x, (1,size(x)...)) #quickfix for the case where axx is not a matrix. It would require an extra singleton dimension

_squeeze!(arr) = dropdims(arr, dims = tuple(findall(size(arr) .== 1)...)) #remove singleton dimensions of an array

"""
$(TYPEDSIGNATURES)

Goes over a solution and an equally-sized array (a "mask") of booleans.
true  -> solution unchanged
false -> changed to NaN (omitted from plotting)
"""
function filter_solutions(solution::Vector,  booleans)
    _total_elements(solution) == _total_elements(booleans) || error("attempt to filter a solution using a wrongly-sized boolean array")
    rules = Dict(1 => 1., 0 => NaN)
    factors = [rules[pt] for pt in booleans]
    return solution .* factors
end

"""
$(TYPEDSIGNATURES)

Takes a `Result` object and a string `f` representing a Symbolics.jl expression.
Returns an array with the values of `f` evaluated for the respective solutions.
Additional substitution rules can be specified in `rules` in the format `("a" => val)` or `(a => val)`
"""
function transform_solutions(res::Result, f::String; rules=Dict())
    # a string is used as input - a macro would not "see" the user's namespace while the user's namespace does not "see" the variables
    transformed = [Vector{ComplexF64}(undef, length(res.solutions[1])) for k in res.solutions] # preallocate

    # define variables in rules in this namespace
    new_keys = declare_variable.(string.(keys(Dict(rules))))
    expr = f isa String ? Num(eval(Meta.parse(f))) : f

    fixed_subs = merge(res.fixed_parameters, Dict(zip(new_keys, values(Dict(rules)))))
    expr = substitute_all(expr, Dict(fixed_subs))

    vars = res.problem.variables
    all_symbols = cat(vars, collect(keys(res.swept_parameters)), dims=1)
    comp_func = build_function(expr, all_symbols)
    f = eval(comp_func)

    # preallocate an array for the numerical values, rewrite parts of it
    # when looping through the solutions
    vals = Vector{ComplexF64}(undef, length(all_symbols))
    n_vars = length(vars)
    n_pars = length(all_symbols) - n_vars

    for idx in CartesianIndices(res.solutions)
        params_values = res.swept_parameters[Tuple(idx)...]
        vals[end-n_pars+1:end] .= params_values # param values are common to all branches
        for (branch,soln) in enumerate(res.solutions[idx])
            vals[1:n_vars] .= soln
            transformed[idx][branch] = Base.invokelatest(f, vals)
        end
    end
    return transformed
end

"""
$(TYPEDSIGNATURES)

For a transformed, physical solution, evaluates multi-solution condition (e.g. Base.maximum or Base.argmax).
If real_function=true (default) the function is assume to be real and small imaginary parts are discarded.
"""
function map_multi_solutions(solution,mapping=nothing; real_function=true)
    length(vcat(solution...)[1])>1  && error("input solutions have wrong dimensions. Please transform first.")#sanity check for solutions still containing the solution tuple u1,v1,... (not transformed)
    isempty(methods(mapping)) && error("mapping is not a Julia function") #sanity check for mapping being a function
    if real_function
        @warn "imaginary parts discarded in multi-solution map"
        return mapping.(real.(filter.(!isnan,solution))) #input is a filtered solution. Remove NaN values for mapping_evaluation
    else
        return mapping.(filter.(!isnan,solution))
    end
end

"Construct a matrix of subplots from a linear subplot array, even if input length and # of rows or columns is not commensurable"
function resize_axes!(f,axs,nrows,ncols)
    gspec = pyimport("matplotlib.gridspec")
    gs = gspec.GridSpec(nrows,ncols,hspace=0.6,wspace=0.6)

    if nrows >= ncols
        for i in 1:nrows
            for j in 0:ncols
                k = i + j*ncols
                if k <= length(axs)
                    axs[k].set_position(gs[k].get_position(f))
                end
            end
        end
    else
        for i in 1:nrows
            for j in 0:ncols
                k = i + j*nrows
                if k <= length(axs)
                    axs[k].set_position(gs[k].get_position(f))
                end
            end
        end
    end
    return axs
end

"inserts underscores in variable names for prettier transformed laTeX strings"
function _prettify_label(res::Result,label::String)
    replace_rules = [string("u",k)=>string("u_",k) for k in 1:length(res.problem.variables)÷2]
    append!(replace_rules,[string("v",k)=>string("v_",k) for k in 1:length(res.problem.variables)÷2])
    return reduce(replace, replace_rules, init=label)
end

"""
    plot_1D_solutions(res::Result;
                        x::String, y::String,
                        xscale="linear",yscale="linear"
                        ,plot_only=["physical"],
                        marker_classification="stable",filename=nothing, kwargs...)

Make a 1D plot of a `Result` object.

Keyword arguments
- `x`, `y`: Expressions to plot on as independent/dependent variables (parsed into Symbolics.jl).
- `xscale`, `yscale`: scale for x/y dimensions (e.g. "linear" or "log")
- `plot_only`: a list of strings corresponding to the solution classes of `Result`. Only solutions which belong to the listed classes are plotted.
- `marker_classification`: A class of the solutions (created by `classify_solutions!`) which is distinguished with different markers. Entering an inequality creates a new class "custom_class".
- `filename`: if different from `nothing`, plotted data and parameter values are exported to `./filename.jld2`. Otherwise, data is returned as a dictionary.
- `kwargs`: any additional keywords arguments for the matplotlib plotting

The strings in `marker_classification` allows the user to stablish custom criteria for binary classification of solutions. For instance, if `marker_classification = "ω^15* sqrt(u1^2 + v1^2) < 0.1"`,
for a system with harmonic variables u1,v1, then solutions are classified as `true` according to that criterion and `false` according to its complement.
"""
function plot_1D_solutions(res::Result; x::String, y::String, xscale="linear",yscale="linear",plot_only=["physical"],marker_classification="stable",filename=nothing,kwargs...)
    _set_plotting_settings()
    length(size(res.solutions)) != 1 && error("1D plots of not-1D datasets are usually a bad idea.")

    X = transform_solutions(res, x)
    Y = transform_solutions(res, y) # first transform, then filter

    #retrieve data from specific branch
    relevant_indices = findall(x->x==true, [!all(isnan.(getindex.(Y, branch))) for branch in 1:length(Y[1])])

    f,ax = subplots(1,1,figsize=(7,4))

    #filtering of the solutions according to keyword arguments
    "stable" in plot_only && "physical" ∉ plot_only && error("Stability is not defined for unphysical solutions!")
    ~(sum([class in keys(res.classes) for class in plot_only])==length(plot_only)) && error("Class not found. Solutions are classified according to ", keys(res.classes))
    ("binary_labels" in plot_only) && print("Binary label is not Boolean and thus was ignored in plotting")

    for class in plot_only
        if  class!="binary_labels"
            Y = filter_solutions.(Y, res.classes[class])
        end
    end

    sol_type, not_sol_type = _classify_plot_data(res, marker_classification)

    Ys, Yu = filter_solutions.(Y, res.classes[sol_type]), filter_solutions.(Y, [.!el for el in res.classes[sol_type]])
    lines = ax.plot(X, Ys, "o"; kwargs...) #Nan solutions are ignored
    ax.set_prop_cycle(nothing) #reset color cycler state (NaN solutions aren't shown but the color cycler runs)
    append!(lines,ax.plot(X, Yu, "x"; kwargs...))

    xdata,ydata = [line.get_xdata() for line in lines], [line.get_ydata() for line in lines]
    markers = [line.get_marker() for line in lines]
    marker_dict = Dict("o"=>sol_type,"x"=>not_sol_type)
    save_dict= Dict(string(x) => xdata,string(y)=>ydata,"marker_dict"=>marker_dict,"markers"=>markers)

    ax.set_xlabel(latexify(x),fontsize=24)
    ax.set_ylabel(latexify(_prettify_label(res,y)),fontsize=24)

    #legend preparation
    ignored_idx = [all(isnan.(line.get_ydata())) for line in lines] # make up a legend with only non ignored entries in the plotter
    Nb   = sum(.~ignored_idx) #number of branches

    palette = plt.rcParams["axes.prop_cycle"].by_key()["color"] # the currently used palette (to match legend and plot colours)
    palette = repeat(plt.rcParams["axes.prop_cycle"].by_key()["color"],Nb ÷ length(palette) + 1)[1:Nb]

    leg_classes = [plt.Line2D([0], [0]; marker="o", color="w", label=sol_type, markerfacecolor="k", kwargs...),
            plt.Line2D([0], [0]; marker="X", color="w", label=not_sol_type,markerfacecolor="k", kwargs...)]

    leg_branches = [plt.Line2D([0], [0]; marker="o", color="w", label=k, markerfacecolor=palette[k], kwargs...) for k in 1:Nb]

    leg = cat(leg_classes, leg_branches..., dims=1)

    ax.legend(handles=leg, bbox_to_anchor=(-0.25, 0.95))
    ax.set_xscale(xscale)
    ax.set_yscale(yscale)
    f.tight_layout()

    !isnothing(filename) ? JLD2.save(_jld2_name(filename), save_dict) : nothing
    return save_dict

end

"""
plot_1D_solutions_spaghetti(res::Result;
                        z::String,
                        filename=nothing)

Produces a "spaghetti plot" of 1D `Result`object`, with harmonic variables in a given plane and the parameter `z` in the perpendicular dimension

Keyword arguments
- `z`: Parameter expression to plot on as dependent variables z (parsed into Symbolics.jl).
- `zscale`: scale for z dimension (e.g. "linear" or "log")
- `z_aspect`: aspect ratio for the parameter dimension
- `kwargs`: any additional keywords arguments for the matplotlib plotting
"""

function plot_1D_solutions_spaghetti(res,z::String,z_scale="linear",zaspect=2,kwargs...)
    var_names = [string(v) for v in res.problem.variables]
    n_dof     = length(var_names)÷2 #number of harmonic variable pairs
    nsolsmax  = sum(any.(classify_branch(res, "physical"))) #maximum number of physical solutions

    fig = plt.figure(figsize=(5*(2*n_dof),5),tight_layout=true)

    legend_elements = [plt.Line2D([1], [1], linestyle="-", color="k", label="stable", markerfacecolor="k", markersize=5),
                           plt.Line2D([1], [1], linestyle="--", color="k", label="unstable",markerfacecolor="k", markersize=5)]

    for i = 0:n_dof-1
        us = transform_solutions(res, var_names[2*i+1])
        vs = transform_solutions(res, var_names[2*i+2])
        zs  = reduce(hcat,transform_solutions(res, z))

        #stable solutions
        uS = filter_solutions.(filter_solutions.(us, res.classes["physical"]), res.classes["stable"])
        vS = filter_solutions.(filter_solutions.(vs, res.classes["physical"]), res.classes["stable"])
        #unstable solutions
        uU = filter_solutions.(filter_solutions.(us, res.classes["physical"]), [.!el for el in res.classes["stable"]])
        vU = filter_solutions.(filter_solutions.(vs, res.classes["physical"]), [.!el for el in res.classes["stable"]])

        #transform into array for easier handling
        uS = reduce(hcat,uS)
        vS = reduce(hcat,vS);

        uU = reduce(hcat,uU)
        vU = reduce(hcat,vU);

        ax = fig[:add_subplot](1,n_dof,i+1, projection="3d")
        for i in 1:nsolsmax
            ax.plot(real.(uS[i,:]),real.(vS[i,:]),zs[i,:]; kwargs...)
            ax.plot(real.(uU[i,:]),real.(vU[i,:]),zs[i,:],ls="--";kwargs...)
        end
        ax.set_xlabel( HarmonicBalance.latexify(_prettify_label(res,var_names[2*i+1])),fontsize=24,labelpad=15)
        ax.set_ylabel( HarmonicBalance.latexify(_prettify_label(res,var_names[2*i+2])),fontsize=24,labelpad=15)
        ax.set_zlabel( HarmonicBalance.latexify(z),fontsize=24,labelpad=15)

        ax.set_zscale(z_scale)
        ax.w_xaxis.set_pane_color((1.0, 1.0, 1.0, 1.0))
        ax.w_yaxis.set_pane_color((1.0, 1.0, 1.0, 1.0))
        ax.w_zaxis.set_pane_color((1.0, 1.0, 1.0, 1.0))

        ax.legend(handles=legend_elements,loc="best")
        ax.set_box_aspect((1, 1, zaspect))
    end
    fig.tight_layout()
end

""" For a plotted `marker_classification`, classify the data in `res`.
If `marker_classification` is an existing class, choose that class. If an inequality, parse it to create a new class."""
function _classify_plot_data(res, marker_classification)
    if marker_classification ∈ keys(res.classes)
        sol_type, not_sol_type = marker_classification, "not " * marker_classification
    elseif ~isnothing(marker_classification) # if an inequality is supplied, create a new class
        try
            classify_solutions!(res, marker_classification, "custom class")
        catch
            error(marker_classification * " is not a valid classification or an inequality!")
        end
        sol_type, not_sol_type = "custom class", "not custom class"
    end
    sol_type, not_sol_type
end

"""
    plot_1D_jacobian_eigenvalues(res::Result; x::String, physical=true, stable=false,marker_re="o",marker_im="X", filename=nothing,kwargs...)

Make a 1D plot of the Jacobian eigenvalues for each of the solutions in a `Result` object.

Keyword arguments

- `x`: The function on the x axis (a string parsed into Symbolics.jl).
- `physical`, `stable`: Booleans specifying whether unphysical and/or unstable solutions are shown.
- `marker_re`, `marker_im`: The markers to use for the Re and Im parts of the eigenvalues.
- `ax`: axis object from `PyCall.PyObject` setting the coordinate system where data will be plotted. If not given, it is created automatically.
- `filename`: if different from `nothing`, plotted data and parameter values are exported to `./filename.jld2`. Otherwise, data is returned as a dictionary.

"""
function plot_1D_jacobian_eigenvalues(res::Result; x::String, physical=true, stable=false,marker_re="o",marker_im="X",ax=nothing, filename=nothing, kwargs...)
    _set_plotting_settings()

    xplot = transform_solutions(res, x) #independent variable to plot

    #filtering of the solutions according to keyword arguments
    !physical && stable && error("Stability is not defined for unphysical solutions!")

    numbers = [Base.not_int.(isnan.(prod.(s))) for s in res.solutions] # true for solution which are not NaN
    to_evaluate = physical ? (stable ? res.classes["stable"] : res.classes["physical"]) : numbers
    to_evaluate = [to_evaluate[i] .* numbers[i] for (i,_) in enumerate(numbers)]

    nsolsmax  = sum(any.(classify_branch(res, "physical"))) #maximum number of physical solutions
    input_ax = ax
    if isnothing(input_ax) #create figure if axes are not provided, otherwise accept input
        f,ax = subplots(1,nsolsmax,figsize=(4*nsolsmax,4))
    end

    if nsolsmax==1
        ax = _add_dim!([ax])
    end

    lines_re,lines_im =[],[]
    for branch in 1:nsolsmax
        ax[branch].axhline(0,ls="dashdot",lw=2,alpha=0.8) #reference value for unstability test
        λs = [NaN*ones(ComplexF64, length(res.problem.variables)) for i in 1:length(res.solutions)]
        for index in 1:length(res.solutions)
            # calculate if in to_evaluate, otherwise leave the NaN
            if to_evaluate[index][branch]
                J = res.jacobian(get_single_solution(res, branch=branch, index=index))
                λs[index] = eigvals(J)
            end
        end

        X_plot = real.([xplot[idx][branch] for idx in 1:length(res.solutions)])
        append!(lines_re,ax[branch].plot(X_plot, real.(λs), string(marker_re,"r"); kwargs...))
        append!(lines_im,ax[branch].plot(X_plot, imag.(λs), string(marker_im,"g"); kwargs...))  #
        ax[branch].set_title(string("solution ", branch),fontsize=12,loc="left");
        ax[branch].set_xlabel(latexify(string(x)),fontsize=24)

    end

    save_dict =  Dict(zip(["re","im"],[Dict("data"=>Dict(),"marker"=>[]) for d in ["re","im"]]))
    for (axis,lines,marker) in zip(["re","im"],[lines_re,lines_im],[marker_re,marker_im])
        xdata = [line.get_xdata() for line in lines]
        ydata = [line.get_ydata() for line in lines]
        save_dict[axis]["data"] = Dict(string(x) => xdata,string(axis," part(eig)")=>ydata)
        save_dict[axis]["marker"] = marker
    end

    if !isnothing(filename)
        JLD2.save(_jld2_name(filename), save_dict)
    else
        return save_dict
    end

    legend_elements = [plt.Line2D([0], [0], marker=marker_re, color="w", label=L"\Re({\mathrm{eig}(J)})",
                        markerfacecolor="r", markersize=7),
                        plt.Line2D([0], [0], marker=marker_im, color="w", label=L"\Im({\mathrm{eig}(J)})",
                        markerfacecolor="g", markersize=7)]
   ax[1].legend(handles=legend_elements;)
end

"""Take a set of conditions and multi-solution maps on solutions transformed by `z`
    Returns `physical` solutions, filtered by conditions and apply multi-solution maps"""
function _preprocess_2D_solutions(res::Result,z=nothing,plot_only=nothing)
    nvar  = length(res.solutions[1,1][1]) #number of variables
    nsols = length(res.solutions[1,1]) #maximum number of solutions
    x,y = collect(values(res.swept_parameters))

    if isnothing(z)  Z = res.solutions   else  Z =  transform_solutions(res, z)  end # first transform, then filter
    Z = real.(filter_solutions.(Z,res.classes["physical"]))

    local map_s #isolate mapping from plot_only array
    if !isnothing(plot_only)
        for class in plot_only
            if  isempty(methods(class)) #class is not a multi-solution mapping
                Z = filter_solutions.(Z, res.classes[class])
            else !isempty(methods(class))
                map_s = class #unpack multisolution mapping to be applied
            end
        end
    end

    if isnothing(z) #quantity to plot depending on the input. Transform the solution structs into tensors that are easier to handle by imshow
        Z  = reshape(reduce(hcat,[reduce(hcat,sol) for sol in Z]),nvar,nsols,length(x),length(y))
        map_s = nothing
    else #transformation of solution is requested
        if Base.@isdefined map_s#class is a multi-solution mapping. Apply only after filtering
            Z = map_multi_solutions(Z,map_s; real_function=true)
        else
            Z = reshape(reduce(hcat,[reduce(hcat,sol) for sol in Z]),nsols,length(x),length(y))
            map_s = nothing
        end
    end
    Z,map_s
end

"""
    plot_2D_solutions(res::Result,ax=nothing; filename=nothing)

Make a 2D plot of each of solutions vs swept parameters for a `Result` object, obtained from a `get_steady_states` applied to a 2D parameter grid.`.

Keyword arguments

- `ax`: axis object from `PyCall.PyObject` setting the coordinate system where data will be plotted. If not given, it is created automatically.
- `filename`: if different from `nothing`, plotted data and parameter values are exported to `./filename.jld2`.  Otherwise, data is returned as a dictionary.
- `z`: The function on the z axis (a string parsed into Symbolics.jl). If `z=nothing`, raw solutions are displayed
- `plot_only`: Array of labels to filter physical solutions (e.g. "stable") and multi-solution methods if `z!=nothing` (e.g. maximum)
"""
function plot_2D_solutions(res::Result; ax=nothing, filename=nothing, z=nothing,plot_only=nothing)
    _set_plotting_settings()
    nvar  = length(res.solutions[1,1][1]) #number of variables
    nsols = length(res.solutions[1,1]) #maximum number of solutions
    x,y = collect(values(res.swept_parameters))
    extent=[x[1],x[end],y[1],y[end]]
    var_names = _get_var_name_labels(res) #variable names for plot labels

    Z,map_s = _preprocess_2D_solutions(res,z,plot_only)

    input_ax = ax #set up axes depending on input

    if isnothing(z)
        nrow,ncol = nsols,nvar
    else
        nrow =1
        if !isnothing(map_s)
            ncol = 1
        else
            ncol = nsols
        end
    end

    if isnothing(input_ax)
        f,ax = subplots(nrow,ncol,figsize=(4*ncol,4*nrow)) #create figure if axes are not provided, otherwise accept input
    else
        length(vec(input_ax)) <= nrow*ncol || error("insufficient # of panels requested, please increase #rows or #cols") #sanity check before any plot is made
    end

    if ncol==1
        ax = _add_dim!([ax])
    end

    px,py = string.(collect(keys(res.swept_parameters)))  #swept parameter strings
    if isnothing(z) #plot solution variables
        save_dict = Dict([string("panel (",m,",",l,")")=> Dict() for m in 1:nvar for l in 1:nsols])
        for m in 1:ncol
            for l in 1:nrow
                a = ax[l,m].imshow(Z[m,l,:,end:-1:1]',extent=extent,aspect="auto")
                colorbar(a,ax=ax[l,m])

                save_dict[string("panel (",m,",",l,")")]= Dict("variable"=>var_names[m],"solution #"=>l,"data"=>a.get_array(),
                            string("(",px,"_min ",px,"_max ",py,"_min ",py,"_max)")=>extent)

            end
            ax[1,m].set_title(latexify(_prettify_label(res,var_names[m])),fontsize=20)
        end
    else #plot transformed solutions or functions thereof
        save_dict = Dict([string("panel (",l,")")=> Dict() for l in 1:ncol])
        for l in 1:ncol
            if !isnothing(map_s)
                a = ax[l].imshow(Z[:,end:-1:1]',extent=extent,aspect="auto")
                ax[1].text(0.05, 0.9, string(Symbol(map_s)),c="w", transform=ax[1].transAxes,fontsize=14)
            else
                a = ax[l].imshow(Z[l,:,end:-1:1]',extent=extent,aspect="auto")
                ax[l].set_title(string("solution ",l),fontsize=18)
            end
            colorbar(a,ax=ax[l])

            save_dict[string("panel (",l,")")]= Dict("variable"=>z,"solution #"=>l,"data"=>a.get_array(),
                        string("(",px,"_min ",px,"_max ",py,"_min ",py,"_max)")=>extent)

        end
        f.suptitle(latexify(_prettify_label(res,z)),fontsize=18,y=1.01)

    end
    for ax_ in ax
        ax_.set_xlabel(latexify(px),fontsize=24);
        ax_.set_ylabel(latexify(py),fontsize=24);
    end

    f.tight_layout()

    !isnothing(filename) ? JLD2.save(_jld2_name(filename), save_dict) : nothing
    return save_dict
end

"""Discrete colorbar preparation for phase diagram"""
function _prepare_colorbar(f,ax,im,Nmax;Nmax_discrete=9)  #creates and inserts in a phase diagram figure f a discrete or continuous colorbar
    if Nmax <= Nmax_discrete
        f.colorbar(im, ax=ax,ticks=collect(1:Nmax), boundaries=collect(1:Nmax+1).-0.5)
    else
        f.colorbar(im, ax=ax)
    end
end

"""
    plot_2D_phase_diagram(res::Result; stable=false,observable="nsols",ax=nothing, filename=nothing)

Make a 2D phase diagram plot of each of solutions vs swept parameters for a `Result` object, obtained from a `get_steady_states` applied to a 2D parameter grid.

Keyword arguments

- `stable`: whether only stable solutions are depicted
- `observable`: reference observable to represent dynamical phases in the problem. If `observable="nsols"`, number of solutions for each point is shown.
   If instead `observable="binary"`, the result of classification of bistrings `[is_stable(solution_1),is_stable(solution_2),...]` is presented (see `classify_binaries!(Result)` function).
- `filename`: if different from `nothing`, plotted data and parameter values are exported to `./filename.jld2`.Otherwise, data is returned as a dictionary.

"""
function plot_2D_phase_diagram(res::Result; stable=false,observable="nsols",ax=nothing, filename=nothing)
    _set_plotting_settings()
    observable_choice = ["nsols", "binary"]
    observable ∈ observable_choice  ||error("Only the following 2D observables are allowed:  ", observable_choice)

    X,Y = collect(values(res.swept_parameters))
    x,y = collect(keys(res.swept_parameters))
    extent = [X[1], X[end], Y[1], Y[end]]

    if observable=="nsols" #number of solutions plot
        if stable
            obs_2D =  sum.(res.classes["stable"])
            pd_title = "number of stable solutions"
        else
            obs_2D =  sum.(res.classes["physical"])
            pd_title = "total number of solutions"
        end
    elseif observable=="binary"
        obs_2D = res.classes["binary_labels"]
        pd_title = "binary classified solutions"
    end

    input_ax = ax
    if isnothing(input_ax) #create figure if axes are not provided, otherwise accept input
        f,ax = subplots(1,1,figsize=(6,5))
    end
    # rearrange num_of_sols to match imshow output
    Nmax = maximum(obs_2D)
    im = ax.imshow(obs_2D[:,end:-1:1]', extent=extent, aspect="auto",vmin=minimum(obs_2D),vmax=Nmax)
    if isnothing(input_ax)
        _prepare_colorbar(f,ax,im,Nmax)
    end

    px,py =string.([x,y])#swept parameter strings

    ax.set_xlabel(latexify(px),fontsize=24)
    ax.set_ylabel(latexify(py),fontsize=24)
    ax.set_title(pd_title)

    save_dict =  Dict("observable"=>observable,"data"=>im.get_array(),
                      string("(",px,"_min ",px,"_max ",py,"_min ",py,"_max)")=>extent)

    !isnothing(filename) ? JLD2.save(_jld2_name(filename),save_dict) : nothing
    return save_dict,im,Nmax

end


text:
using IntcodeVM;
using DelimitedFiles;
using Combinatorics;
using Sockets;

const code = vec(readdlm(length(ARGS) >= 1 ? ARGS[1] : "input.txt", ',', Int, '\n'))
const port = length(ARGS) >= 2 ? parse(Int,  ARGS[2]) : 60000

# Open a server socket and launch the Intcode VM as a thread.
@async begin
    vm_listener = listen(port)
    vm_socket = accept(vm_listener)
    println("Connected!")
    IntcodeVM.run(code, in=vm_socket, out=vm_socket)
    close(vm_socket)
    close(vm_listener)
end

println("Hull painting robot is listening for commands on port $(port).")

position = 0
heading = im
const white = Set()
const socket = connect(port)
const colors = Dict(0 => "black", 1 => "white", "0" => "BLACK", "1" => "WHITE")
while true
    println("0: Black, 1: White")
    print("HULLPNTR.EXE> ")
    command = readline()
    if command in ["quit", "q", "exit", "bye"]
        break
    elseif command == ""
        continue
    end
    println(socket, command)
    (color, direction) = (readline(socket), readline(socket))
    println("This panel is $(colors[color]).")
    println("The robot should turn $(direction == "0" ? "left" : "right") 90 degrees.")
    println()
end

close(socket)


text:
"""
    QuadratureVI(;ϵ::T=1e-5, nGaussHermite::Integer=20, clipping=Inf, natural::Bool=true, optimiser=Momentum(0.0001))

Variational Inference solver by approximating gradients via numerical integration via Quadrature

## Keyword arguments
- `ϵ::T` : convergence criteria
- `nGaussHermite::Int` : Number of points for the integral estimation
- `clipping::Real` : Limit the gradients values to avoid overshooting
- `natural::Bool` : Use natural gradients
- `optimiser` : Optimiser used for the variational updates. Should be an Optimiser object from the [Flux.jl](https://github.com/FluxML/Flux.jl) library, see list here [Optimisers](https://fluxml.ai/Flux.jl/stable/training/optimisers/) and on [this list](https://github.com/theogf/AugmentedGaussianProcesses.jl/tree/master/src/inference/optimisers.jl). Default is `Momentum(0.0001)`
"""
mutable struct QuadratureVI{T} <: NumericalVI{T}
    n_points::Int # Number of points for the quadrature
    nodes::Vector{T} # Nodes locations
    weights::Vector{T} # Weights for each node
    clipping::T # Clipping value of the gradient
    ϵ::T #Convergence criteria
    n_iter::Integer #Number of steps performed
    stoch::Bool #Use of mini-batches
    batchsize::Int #Size of mini-batches
    ρ::T #Stochastic Coefficient
    NaturalGradient::Bool
    HyperParametersUpdated::Bool # Flag for updating kernel matrices
    vi_opt::NVIOptimizer

    function QuadratureVI{T}(
        ϵ::T,
        n_points::Integer,
        optimiser,
        Stochastic::Bool,
        clipping::Real,
        batchsize::Int,
        natural::Bool,
    ) where {T}
        gh = gausshermite(n_points)
        return new{T}(
            n_points,
            gh[1] .* sqrt2,
            gh[2] ./ sqrtπ,
            clipping,
            ϵ,
            0,
            Stochastic,
            batchsize,
            one(T),
            natural,
            true,
            NVIOptimizer(optimiser),
        )
    end
end

function QuadratureVI(;
    ϵ::T=1e-5,
    nGaussHermite::Integer=100,
    optimiser=Momentum(1e-5),
    clipping::Real=0.0,
    natural::Bool=true,
) where {T<:Real}
    return QuadratureVI{T}(ϵ, nGaussHermite, optimiser, false, clipping, 1, natural)
end

"""
    QuadratureSVI(nMinibatch::Int; ϵ::T=1e-5, nGaussHermite::Int=20, clipping=Inf, natural=true, optimiser=Momentum(0.0001))

Stochastic Variational Inference solver by approximating gradients via numerical integration via Gaussian Quadrature.
See [`QuadratureVI`](@ref) for a more detailed reference.

## Arguments
-`nMinibatch::Integer` : Number of samples per mini-batches

## Keyword arguments
- `ϵ::T` : convergence criteria, which can be user defined
- `nGaussHermite::Int` : Number of points for the integral estimation (for the QuadratureVI)
- `natural::Bool` : Use natural gradients
- `optimiser` : Optimiser used for the variational updates. Should be an Optimiser object from the [Flux.jl](https://github.com/FluxML/Flux.jl) library, see list here [Optimisers](https://fluxml.ai/Flux.jl/stable/training/optimisers/) and on [this list](https://github.com/theogf/AugmentedGaussianProcesses.jl/tree/master/src/inference/optimisers.jl). Default is `Momentum(0.0001)`
"""
function QuadratureSVI(
    nMinibatch::Integer;
    ϵ::T=1e-5,
    nGaussHermite::Integer=100,
    optimiser=Momentum(1e-5),
    clipping::Real=0.0,
    natural=true,
) where {T<:Real}
    return QuadratureVI{T}(ϵ, nGaussHermite, optimiser, true, clipping, nMinibatch, natural)
end

function expec_loglikelihood(
    l::AbstractLikelihood, i::QuadratureVI, y, μ::AbstractVector, diagΣ::AbstractVector
)
    return mapreduce(+, y, μ, diagΣ) do y, μ, σ²
        apply_quad(i, y, μ, σ², l)
    end
end

function apply_quad(i::QuadratureVI, y::Real, μ::Real, σ²::Real, l::AbstractLikelihood)
    xs = i.nodes * sqrt(σ²) .+ μ
    return dot(i.weights, loglikelihood.(Ref(l), y, xs))
    # return mapreduce((w, x) -> w * Distributions.loglikelihood(l, y, x), +, i.weights, xs)# loglikelihood.(l, y, x))
end

function grad_expectations!(m::AbstractGPModel{T,L,<:QuadratureVI}, state, y) where {T,L}
    for (gp, opt_state, kernel_matrices) in zip(m.f, state.opt_state, state.kernel_matrices)
        μ = mean_f(gp, kernel_matrices)
        Σ = var_f(gp, kernel_matrices)
        for i in 1:batchsize(inference(m))
            opt_state.ν[i], opt_state.λ[i] = grad_quad(
                likelihood(m), y[i], μ[i], Σ[i], inference(m)
            )
        end
    end
end

# Compute the first and second derivative of the log-likelihood using the quadrature nodes
function grad_quad(l::AbstractLikelihood, y::Real, μ::Real, σ²::Real, i::AbstractInference)
    x = i.nodes * sqrt(max(σ², zero(σ²))) .+ μ
    Edloglike = dot(i.weights, ∇loglikehood.(l, y, x))
    Ed²loglike = dot(i.weights, hessloglikehood.(l, y, x))
    if i.clipping != 0
        return (
            abs(Edloglike) > i.clipping ? sign(Edloglike) * i.clipping : -Edloglike,
            abs(Ed²loglike) > i.clipping ? sign(Ed²loglike) * i.clipping : -Ed²loglike,
        )
    else
        return -Edloglike, Ed²loglike
    end
end


text:
# BONDGRAPH
struct BondGraph <: g.AbstractGraph{Int64}
    type::Symbol
    name::Symbol
    nodes::Vector{T} where {T<:AbstractNode}
    bonds::Vector{Bond}
end
BondGraph(name = :BG) = BondGraph(:BG, Symbol(name), AbstractNode[], Bond[])

# PROPERTIES
type(bg::BondGraph) = bg.type

name(bg::BondGraph) = bg.name

nodes(bg::BondGraph) = bg.nodes

bonds(bg::BondGraph) = bg.bonds

# AbstractNode properties
parameters(bg::BondGraph) = [p for c in components(bg) for p in parameters(c)]

states(bg::BondGraph) = [x for c in components(bg) for x in states(c)]

function equations(bg::BondGraph; simplify_eqs = true)
    isempty(bg.nodes) && return Equation[]
    sys = ODESystem(bg; simplify_eqs = simplify_eqs)
    return equations(sys)
end

# Filtering
components(bg::BondGraph) = filter(x -> x isa Component, bg.nodes)
junctions(bg::BondGraph) = filter(x -> x isa Junction, bg.nodes)

getnodes(bg::BondGraph, T::DataType) = filter(x -> x isa T, bg.nodes)
getnodes(bg::BondGraph, t::Symbol) = filter(x -> type(x) == t, bg.nodes)
getnodes(bg::BondGraph, n) = filter(x -> name(x) == Symbol(n), bg.nodes)

getbonds(bg::BondGraph, t::Tuple) = getbonds(bg, t[1], t[2])
getbonds(bg::BondGraph, n1::AbstractNode, n2::AbstractNode) = filter(b -> n1 in b && n2 in b, bg.bonds)

# Base functions
show(io::IO, bg::BondGraph) = print(io, "BondGraph $(bg.type):$(bg.name) ($(g.nv(bg)) Nodes, $(g.ne(bg)) Bonds)")

# Easier referencing systems using a.b notation
function getproperty(bg::BondGraph, sym::Symbol)
    # Calling getfield explicitly avoids using "a.b" and causing a StackOverflowError
    allnodes = getfield(bg, :nodes)
    names = [getfield(n, :name) for n in allnodes]
    symnodes = allnodes[names.==sym]
    if isempty(symnodes)
        return getfield(bg, sym)
    elseif length(symnodes) == 1
        return symnodes[1]
    else
        return symnodes
    end
end

# BONDGRAPH NODE
struct BondGraphNode <: AbstractNode
    bondgraph::BondGraph
    type::Symbol
    name::Symbol
    exposed::Vector{SourceSensor}
    freeports::Vector{Bool}
    vertex::RefValue{Int}
end
function BondGraphNode(bg::BondGraph, name = name(bg); vertex::Int = 0, deepcopy = false)
    _bg = deepcopy ? deepcopy(bg) : bg

    exposed_ports = getnodes(_bg, SourceSensor)
    freeports = fill(true, length(exposed_ports))

    BondGraphNode(_bg, :BG, Symbol(name), exposed_ports, freeports, Ref(vertex))
end

# Easier referencing systems using a.b notation
function getproperty(bgn::BondGraphNode, sym::Symbol)
    bg = getfield(bgn, :bondgraph)
    try
        return getproperty(bg, sym)
    catch
        return getfield(bgn, sym)
    end
end

exposed(bgn::BondGraphNode) = bgn.exposed

max_stars_repo_name: mic-pan/BondGraphs

max_issues_repo_name: mic-pan/BondGraphs

max_forks_repo_name: mic-pan/BondGraphs


text:
# Properties of dry air
const Planet = CLIMAParameters.Planet
Planet.molmass_dryair(ps::AbstractEarthParameterSet) = 28.97e-3
Planet.R_d(ps::AbstractEarthParameterSet)            = gas_constant() / Planet.molmass_dryair(ps)
Planet.kappa_d(ps::AbstractEarthParameterSet)        = 2 / 7
Planet.cp_d(ps::AbstractEarthParameterSet)           = Planet.R_d(ps) / Planet.kappa_d(ps)
Planet.cv_d(ps::AbstractEarthParameterSet)           = Planet.cp_d(ps) - Planet.R_d(ps)

# Properties of water
Planet.ρ_cloud_liq(ps::AbstractEarthParameterSet)           = 1e3
Planet.ρ_cloud_ice(ps::AbstractEarthParameterSet)           = 916.7
Planet.molmass_water(ps::AbstractEarthParameterSet)         = 18.01528e-3
Planet.molmass_ratio(ps::AbstractEarthParameterSet)         = Planet.molmass_dryair(ps) / Planet.molmass_water(ps)
Planet.R_v(ps::AbstractEarthParameterSet)                   = gas_constant() / Planet.molmass_water(ps)
Planet.cp_v(ps::AbstractEarthParameterSet)                  = 1859
Planet.cp_l(ps::AbstractEarthParameterSet)                  = 4181
Planet.cp_i(ps::AbstractEarthParameterSet)                  = 2100
Planet.cv_v(ps::AbstractEarthParameterSet)                  = Planet.cp_v(ps) - Planet.R_v(ps)
Planet.cv_l(ps::AbstractEarthParameterSet)                  = Planet.cp_l(ps)
Planet.cv_i(ps::AbstractEarthParameterSet)                  = Planet.cp_i(ps)
Planet.T_freeze(ps::AbstractEarthParameterSet)              = 273.15
Planet.T_min(ps::AbstractEarthParameterSet)                 = 150.0
Planet.T_max(ps::AbstractEarthParameterSet)                 = 1000.0
Planet.T_icenuc(ps::AbstractEarthParameterSet)              = 233.00
Planet.T_triple(ps::AbstractEarthParameterSet)              = 273.16
Planet.T_0(ps::AbstractEarthParameterSet)                   = Planet.T_triple(ps)
Planet.LH_v0(ps::AbstractEarthParameterSet)                 = 2.5008e6
Planet.LH_s0(ps::AbstractEarthParameterSet)                 = 2.8344e6
Planet.LH_f0(ps::AbstractEarthParameterSet)                 = Planet.LH_s0(ps) - Planet.LH_v0(ps)
Planet.e_int_v0(ps::AbstractEarthParameterSet)              = Planet.LH_v0(ps) - Planet.R_v(ps) * Planet.T_0(ps)
Planet.e_int_i0(ps::AbstractEarthParameterSet)              = Planet.LH_f0(ps)
Planet.press_triple(ps::AbstractEarthParameterSet)          = 611.657
Planet.surface_tension_coeff(ps::AbstractEarthParameterSet) = 0.072

Planet.entropy_dry_air(ps::AbstractEarthParameterSet)       = 6864.8
Planet.entropy_water_vapor(ps::AbstractEarthParameterSet)   = 10513.6
Planet.entropy_reference_temperature(ps::AbstractEarthParameterSet)  = 298.15

# Properties of sea water
Planet.ρ_ocean(ps::AbstractEarthParameterSet)        = 1.035e3
Planet.cp_ocean(ps::AbstractEarthParameterSet)       = 3989.25

# Planetary parameters
Planet.planet_radius(ps::AbstractEarthParameterSet)         = 6.371e6
Planet.day(ps::AbstractEarthParameterSet)                   = 86400
Planet.Omega(ps::AbstractEarthParameterSet)                 = 7.2921159e-5
Planet.grav(ps::AbstractEarthParameterSet)                  = 9.81
Planet.year_anom(ps::AbstractEarthParameterSet)             = 365.26 * Planet.day(ps)
Planet.orbit_semimaj(ps::AbstractEarthParameterSet)         = 1 * astro_unit()
Planet.tot_solar_irrad(ps::AbstractEarthParameterSet)       = 1362.0
Planet.epoch(ps::AbstractEarthParameterSet)                 = 2451545.0 * Planet.day(ps)
Planet.mean_anom_epoch(ps::AbstractEarthParameterSet)       = deg2rad(357.52911)
Planet.obliq_epoch(ps::AbstractEarthParameterSet)           = deg2rad(23.432777778)
Planet.lon_perihelion_epoch(ps::AbstractEarthParameterSet)  = deg2rad(282.937348)
Planet.eccentricity_epoch(ps::AbstractEarthParameterSet)    = 0.016708634
Planet.lon_perihelion(ps::AbstractEarthParameterSet)        = deg2rad(282.937348)
Planet.MSLP(ps::AbstractEarthParameterSet)                  = 1.01325e5
Planet.T_surf_ref(ps::AbstractEarthParameterSet)            = 290.0
Planet.T_min_ref(ps::AbstractEarthParameterSet)             = 220.0