text:
# This file is a part of ValueShapes.jl, licensed under the MIT License (MIT).

@inline _varoffset_cumsum_impl(s, x, y, rest...) = (s, _varoffset_cumsum_impl(s+x, y, rest...)...)
@inline _varoffset_cumsum_impl(s,x) = (s,)
@inline _varoffset_cumsum_impl(s) = ()
@inline _varoffset_cumsum(x::Tuple) = _varoffset_cumsum_impl(0, x...)

"""
    NamedTupleShape{names,...} <: AbstractValueShape

Defines the shape of a `NamedTuple` (resp.  set of variables, parameters,
etc.).

Constructors:

    NamedTupleShape(name1 = shape1::AbstractValueShape, ...)
    NamedTupleShape(named_shapes::NamedTuple)

Example:

```julia
shape = NamedTupleShape(
    a = ScalarShape{Real}(),
    b = ArrayShape{Real}(2, 3),
    c = ConstValueShape(42)
)

data = VectorOfSimilarVectors{Float64}(shape)
resize!(data, 10)
rand!(flatview(data))
table = shape.(data)
fill!(table.a, 4.2)
all(x -> x == 4.2, view(flatview(data), 1, :))
```

See also the documentation of [`AbstractValueShape`](@ref).
"""
struct NamedTupleShape{names,AT<:(NTuple{N,ValueAccessor} where N)} <: AbstractValueShape
    _accessors::NamedTuple{names,AT}
    _flatdof::Int

    @inline function NamedTupleShape(shape::NamedTuple{names,<:NTuple{N,AbstractValueShape}}) where {names,N}
        labels = keys(shape)
        shapes = values(shape)
        shapelengths = map(totalndof, shapes)
        offsets = _varoffset_cumsum(shapelengths)
        accessors = map(ValueAccessor, shapes, offsets)
        # acclengths = map(x -> x.len, accessors)
        # @assert shapelengths == acclengths
        n_flattened = sum(shapelengths)
        named_accessors = NamedTuple{labels}(accessors)
        new{names,typeof(accessors)}(named_accessors, n_flattened)
    end
end

export NamedTupleShape

@inline NamedTupleShape(;named_shapes...) = NamedTupleShape(values(named_shapes))

@inline _accessors(x::NamedTupleShape) = getfield(x, :_accessors)
@inline _flatdof(x::NamedTupleShape) = getfield(x, :_flatdof)

function Base.show(io::IO, shape::NamedTupleShape)
    print(io, Base.typename(typeof(shape)).name, "(")
    show(io, (;shape...))
    print(io, ")")
end

@inline totalndof(shape::NamedTupleShape) = _flatdof(shape)

@inline Base.keys(shape::NamedTupleShape) = keys(_accessors(shape))

@inline Base.values(shape::NamedTupleShape) = values(_accessors(shape))

@inline Base.getindex(d::NamedTupleShape, k::Symbol) = _accessors(d)[k]

@inline function Base.getproperty(shape::NamedTupleShape, p::Symbol)
    # Need to include internal fields of NamedTupleShape to make Zygote happy (ToDo: still true?):
    if p == :_accessors
        getfield(shape, :_accessors)
    elseif p == :_flatdof
        getfield(shape, :_flatdof)
    else
        getproperty(_accessors(shape), p)
    end
end

@inline function Base.propertynames(shape::NamedTupleShape, private::Bool = false)
    names = propertynames(_accessors(shape))
    if private
        (names..., :_flatdof, :_accessors)
    else
        names
    end
end

@inline Base.length(shape::NamedTupleShape) = length(_accessors(shape))

@inline Base.getindex(shape::NamedTupleShape, i::Integer) = getindex(_accessors(shape), i)

@inline Base.map(f, shape::NamedTupleShape) = map(f, _accessors(shape))

function Base.merge(a::NamedTuple, shape::NamedTupleShape{names}) where {names}
    merge(a, NamedTuple{names}(map(x -> valshape(x), values(shape))))
end

Base.merge(a::NamedTupleShape) = a
Base.merge(a::NamedTupleShape, b::NamedTupleShape, cs::NamedTupleShape...) = merge(NamedTupleShape(;a..., b...), cs...)

import Base.<=

function <=(a::NamedTupleShape{names}, b::NamedTupleShape{names}) where {names}
    all(map((a, b) -> a.offset == b.offset && a.shape <= b.shape, values(a), values(b)))
end

@inline <=(a::NamedTupleShape, b::NamedTupleShape) = false

valshape(x::NamedTuple) = NamedTupleShape(map(valshape, x))

(shape::NamedTupleShape)(::UndefInitializer) = map(x -> valshape(x)(undef), shape)

Base.@propagate_inbounds (shape::NamedTupleShape)(data::AbstractVector{<:Real}) = ShapedAsNT(data, shape)

@inline _multi_promote_type() = Nothing
@inline _multi_promote_type(T::Type) = T
@inline _multi_promote_type(T::Type, U::Type, rest::Type...) = promote_type(T, _multi_promote_type(U, rest...))

@inline default_unshaped_eltype(shape::NamedTupleShape) =
    _multi_promote_type(map(default_unshaped_eltype, values(shape))...)

@inline shaped_type(shape::NamedTupleShape{names}, ::Type{T}) where {names,T<:Real} =
    NamedTuple{names,Tuple{map(acc -> shaped_type(acc.shape, T), values(_accessors(shape)))...}}

function unshaped(x::NamedTuple{names}, shape::NamedTupleShape{names}) where names
    # ToDo: Improve performance of return type inference
    T = default_unshaped_eltype(shape)
    U = default_unshaped_eltype(valshape(x))
    R = promote_type(T, U)

    x_unshaped = Vector{R}(undef, totalndof(shape)...)
    shape(x_unshaped)[] = x
    x_unshaped
end

function unshaped(x::AbstractArray{<:NamedTuple{names},0}, shape::NamedTupleShape{names}) where names
    unshaped(x[], shape)
end

replace_const_shapes(f::Function, shape::NamedTupleShape) = NamedTupleShape(map(s -> replace_const_shapes(f, s), (;shape...)))

"""
    ShapedAsNT{T<:NamedTuple,...} <: AbstractArray{T,0}

View of an `AbstractVector{<:Real}` as a zero-dimensional Array containing a
`NamedTuple`, according to a specified [`NamedTupleShape`](@ref).

Constructors:

    ShapedAsNT(data::AbstractVector{<:Real}, shape::NamedTupleShape)

    shape(data)

The resulting `ShapedAsNT` shares memory with `data`. It takes the form of a
(virtual) zero-dimensional Array to make the contents as editable as `data`
itself (compared to a standard immutable NamedTuple):

```julia
x = (a = 42, b = rand(1:9, 2, 3))
shape = valshape(x)
data = Vector{Int}(undef, shape)
y = shape(data)
@assert y isa ShapedAsNT
y[] = x
@assert y[] == x
y.a = 22
y.a[] = 33
@assert shape(data) == y
@assert unshaped(y) === data
```

Use `unshaped(x)` to access `data` directly.

See also [`ShapedAsNTArray`](@ref).
"""
struct ShapedAsNT{T<:NamedTuple,D<:AbstractVector{<:Real},S<:NamedTupleShape} <: AbstractArray{T,0}
    __internal_data::D
    __internal_valshape::S
end

export ShapedAsNT

Base.@propagate_inbounds function ShapedAsNT(data::D, shape::S) where {T<:Real,D<:AbstractVector{T},S<:NamedTupleShape}
    @boundscheck _checkcompat(shape, data)
    NT_T = shaped_type(shape, T)
    ShapedAsNT{NT_T,D,S}(data, shape)
end

@inline _data(A::ShapedAsNT) = getfield(A, :__internal_data)
@inline _valshape(A::ShapedAsNT) = getfield(A, :__internal_valshape)

@inline valshape(A::ShapedAsNT) = _valshape(A)
@inline unshaped(A::ShapedAsNT) = _data(A)

function unshaped(x::ShapedAsNT{<:NamedTuple{names}}, shape::NamedTupleShape{names}) where names
    valshape(x) <= shape || throw(ArgumentError("Shape of value not compatible with given shape"))
    unshaped(x)
end

# ToDo: Move index calculation to separate function with no-op custom pullback to increase performance?
Base.@propagate_inbounds function Base.getproperty(A::ShapedAsNT, p::Symbol)
    # Need to include internal fields of ShapedAsNT to make Zygote happy (ToDo: still true?):
    if p == :__internal_data
        getfield(A, :__internal_data)
    elseif p == :__internal_valshape
        getfield(A, :__internal_valshape)
    else
        data = _data(A)
        shape = _valshape(A)
        va = getproperty(_accessors(shape), p)
        view(data, va)
    end
end

Base.@propagate_inbounds function Base.setproperty!(A::ShapedAsNT, p::Symbol, x)
    data = _data(A)
    shape = _valshape(A)
    va = getproperty(_accessors(shape), p)
    setindex!(data, x, va)
    A
end

@inline function Base.propertynames(A::ShapedAsNT, private::Bool = false)
    names = Base.propertynames(_valshape(A))
    if private
        (names..., :__internal_data, :__internal_valshape)
    else
        names
    end
end

@inline Base.size(A::ShapedAsNT) = ()
@inline Base.IndexStyle(::Type{<:ShapedAsNT}) = IndexLinear()

Base.@propagate_inbounds function _apply_ntshape_copy(data::AbstractVector{<:Real}, shape::NamedTupleShape)
    accessors = _accessors(shape)
    map(va -> getindex(data, va), accessors)
end

Base.@propagate_inbounds Base.getindex(A::ShapedAsNT) = _apply_ntshape_copy(_data(A), _valshape(A))

Base.@propagate_inbounds function Base.getindex(A::ShapedAsNT, i::Integer)
    @boundscheck Base.checkbounds(A, i)
    getindex(A)
end

Base.@propagate_inbounds function Base.getindex(A::ShapedAsNT, i::Union{AbstractArray,Colon})
    @boundscheck Base.checkbounds(A, i)
    [getindex(A)]
end

Base.@propagate_inbounds _apply_ntshape_view(A::AbstractVector{<:Real}, shape::NamedTupleShape) =
    ShapedAsNT(A, shape)

Base.@propagate_inbounds Base.view(A::ShapedAsNT) = A

Base.@propagate_inbounds function Base.view(A::ShapedAsNT, i::Integer)
    @boundscheck Base.checkbounds(A, i)
    view(A)
end

Base.@propagate_inbounds function Base.view(A::ShapedAsNT, i::Union{AbstractArray,Colon})
    @boundscheck Base.checkbounds(A, i)
    ShapedAsNTArray(view([_data(A)], :), _valshape(A))
end

Base.@propagate_inbounds function Base.setindex!(A::ShapedAsNT{<:NamedTuple{names}}, x::NamedTuple{names}) where {names}
    if @generated
        Expr(:block, map(p -> :(A.$p = x.$p), names)...)
    else
        # Shouldn't be used, ideally
        @assert false
        map(n -> setproperty!(A, n, getproperty(x, n)), names)
    end

    A
end

Base.@propagate_inbounds Base.setindex!(A::ShapedAsNT{T}, x) where {T} = setindex!(A, convert(T, x))

Base.@propagate_inbounds function Base.setindex!(A::ShapedAsNT, x, i::Integer)
    @boundscheck Base.checkbounds(A, i)
    setindex!(A, x)
end

Base.similar(A::ShapedAsNT{T}, ::Type{T}, ::Tuple{}) where T =
    ShapedAsNT(similar(_data(A)), _valshape(A))

Base.show(io::IO, ::MIME"text/plain", A::ShapedAsNT) = show(io, A)

function Base.show(io::IO, A::ShapedAsNT)
    print(io, Base.typename(typeof(A)).name, "(")
    show(io, A[])
    print(io, ")")
end

Base.copy(A::ShapedAsNT) = ShapedAsNT(copy(_data(A)), _valshape(A))

function _shaped_nt_ΔΩ(ΔΩ::NamedTuple{names}, result::NamedTuple{names}, x::ShapedAsNT{<:NamedTuple{names}}) where names
    vs = gradient_shape(valshape(x))
    ΔΩ_fixed_values = map(values(result), values(ΔΩ), values(vs)) do result, df, acc
        isnothing(df) || valshape(acc) isa ConstValueShape ? ValueShapes.const_zero(result) : df
    end
    ΔΩ_fixed = NamedTuple{names}(ΔΩ_fixed_values)
    ΔΩ_shaped = vs(similar(unshaped(x)))
    ΔΩ_shaped[] = ΔΩ_fixed
    return ΔΩ_shaped
end

# Zygote will currently ignore this, see Zygote.jl issue #811:
function ChainRulesCore.rrule(::typeof(Base.getindex), x::ShapedAsNT)
    result = x[]
    shapedasnt_getindex_pullback(ΔΩ) = (NoTangent(), _shaped_nt_ΔΩ(unthunk(ΔΩ), result, x))
    return result, shapedasnt_getindex_pullback
end
#
# So need to use
ZygoteRules.@adjoint Base.getindex(x::ShapedAsNT) = begin
    result = x[]
    shapedasnt_getindex_pullback(ΔΩ::NamedTuple) = begin
        (_shaped_nt_ΔΩ(ΔΩ, result, x),)
    end
    (result, shapedasnt_getindex_pullback)
end

function ChainRulesCore.rrule(::Type{ShapedAsNT}, A::AbstractVector{<:Real}, vs::NamedTupleShape{names}) where names
    result = ShapedAsNT(A, vs)
    function shapedasnt_pullback_impl(ΔΩ::Union{ShapedAsNT{<:NamedTuple{names}},NamedTuple{names}})
        (NoTangent(), unshaped(ΔΩ, gradient_shape(vs)), nothing)
    end
    function shapedasnt_pullback_impl(ΔΩ_c::Tangent{Any,<:NamedTuple{names}})
        ΔΩ = NamedTuple{names}((ΔΩ_c...,))
        shapedasnt_pullback_impl(ΔΩ)
    end
    function shapedasnt_pullback_impl(ΔΩ_c::Tangent{Any,<:NamedTuple{(:__internal_data, :__internal_valshape)}})
        @assert ΔΩ_c.__internal_valshape == NoTangent() || ΔΩ_c.__internal_valshape == ZeroTangent()
        (NoTangent(), ΔΩ_c.__internal_data, nothing)
    end
    shapedasnt_pullback(ΔΩ) = shapedasnt_pullback_impl(unthunk(ΔΩ))
    return result, shapedasnt_pullback
end

function ChainRulesCore.rrule(::typeof(unshaped), A::ShapedAsNT)
    result = unshaped(A)
    vs = valshape(A)
    unshaped_nt_pullback(ΔΩ::AbstractArray{<:Real}) = (NoTangent(), gradient_shape(vs)(ΔΩ))
    return result, unshaped_nt_pullback
end

function ChainRulesCore.rrule(::typeof(unshaped), A::ShapedAsNT, vs::NamedTupleShape)
    result = unshaped(A, vs)
    unshaped_nt_pullback(ΔΩ::AbstractArray{<:Real}) = (NoTangent(), gradient_shape(vs)(ΔΩ), nothing)
    return result, unshaped_nt_pullback
end

function ChainRulesCore.rrule(::typeof(unshaped), A::NamedTuple, vs::NamedTupleShape)
    result = unshaped(A, vs)
    unshaped_nt_pullback(ΔΩ::AbstractArray{<:Real}) = (NoTangent(), gradient_shape(vs)(ΔΩ)[], nothing)
    return result, unshaped_nt_pullback
end

"""
    ShapedAsNTArray{T<:NamedTuple,...} <: AbstractArray{T,N}

View of an `AbstractArray{<:AbstractVector{<:Real},N}` as an array of
`NamedTuple`s, according to a specified [`NamedTupleShape`](@ref).

`ShapedAsNTArray` implements the `Tables` API. Semantically, it acts a
broadcasted [`ShapedAsNT`](@ref).

Constructors:

    ShapedAsNTArray(
        data::AbstractArray{<:AbstractVector{<:Real},
        shape::NamedTupleShape
    )

    shape.(data)

The resulting `ShapedAsNTArray` shares memory with `data`:

```julia
using ArraysOfArrays, Tables, TypedTables

X = [
    (a = 42, b = rand(1:9, 2, 3))
    (a = 11, b = rand(1:9, 2, 3))
]

shape = valshape(X[1])
data = nestedview(Array{Int}(undef, totalndof(shape), 2))
Y = shape.(data)
@assert Y isa ShapedAsNTArray
Y[:] = X
@assert Y[1] == X[1] == shape(data[1])[]
@assert Y.a == [42, 11]
Tables.columns(Y)
@assert unshaped.(Y) === data
@assert Table(Y) isa TypedTables.Table
```

Use `unshaped.(Y)` to access `data` directly.

`Tables.columns(Y)` will return a `NamedTuple` of columns. They will contain
a copy the data, using a memory layout as contiguous as possible for each
column.
"""
struct ShapedAsNTArray{T<:NamedTuple,N,D<:AbstractArray{<:AbstractVector{<:Real},N},S<:NamedTupleShape} <: AbstractArray{T,N}
    __internal_data::D
    __internal_elshape::S
end

export ShapedAsNTArray

function ShapedAsNTArray(data::D, shape::S) where {N,T<:Real,D<:AbstractArray{<:AbstractVector{T},N},S<:NamedTupleShape}
    NT_T = shaped_type(shape, T)
    ShapedAsNTArray{NT_T,N,D,S}(data, shape)
end

# Specialize (::NamedTupleShape).(::AbstractVector{<:AbstractVector}):
Base.copy(instance::VSBroadcasted1{1,<:NamedTupleShape,AbstractVector{<:AbstractVector{<:Real}}}) =
    ShapedAsNTArray(instance.args[1], instance.f)

@inline _data(A::ShapedAsNTArray) = getfield(A, :__internal_data)
@inline _elshape(A::ShapedAsNTArray) = getfield(A, :__internal_elshape)

@inline elshape(A::ShapedAsNTArray) = _elshape(A)

@inline _bcasted_unshaped(A::ShapedAsNTArray) = _data(A)

Base.copy(instance::VSBroadcasted1{N,typeof(unshaped),ShapedAsNTArray{T,N}}) where {T,N} =
    _bcasted_unshaped(instance.args[1])

@inline function Base.getproperty(A::ShapedAsNTArray, p::Symbol)
    # Need to include internal fields of ShapedAsNTArray to make Zygote happy (ToDo: still true?):
    if p == :__internal_data
        getfield(A, :__internal_data)
    elseif p == :__internal_elshape
        getfield(A, :__internal_elshape)
    else
        data = _data(A)
        shape = _elshape(A)
        va = getproperty(_accessors(shape), p)
        view.(data, Ref(va))
    end
end

@inline function Base.propertynames(A::ShapedAsNTArray, private::Bool = false)
    names = Base.propertynames(_elshape(A))
    if private
        (names..., :__internal_data, :__internal_elshape)
    else
        names
    end
end

@inline Base.size(A::ShapedAsNTArray) = size(_data(A))
@inline Base.axes(A::ShapedAsNTArray) = axes(_data(A))
@inline Base.IndexStyle(::Type{<:ShapedAsNTArray{T,N,D}}) where {T,N,D} = IndexStyle(D)

Base.@propagate_inbounds _apply_ntshape_copy(data::AbstractArray{<:AbstractVector{<:Real}}, shape::NamedTupleShape) =
    ShapedAsNTArray(data, shape)

Base.getindex(A::ShapedAsNTArray, idxs...) = _apply_ntshape_copy(getindex(_data(A), idxs...), _elshape(A))

Base.@propagate_inbounds _apply_ntshape_view(data::AbstractArray{<:AbstractVector{<:Real}}, shape::NamedTupleShape) =
    ShapedAsNTArray(data, shape)

Base.@propagate_inbounds _apply_ntshape_view(data::AbstractArray{<:AbstractVector{<:Real}, 0}, shape::NamedTupleShape) =
    ShapedAsNT(data[], shape)

Base.view(A::ShapedAsNTArray, idxs...) = _apply_ntshape_view(view(_data(A), idxs...), _elshape(A))

function Base.setindex!(A::ShapedAsNTArray, x, idxs::Integer...)
    A_idxs = ShapedAsNT(getindex(_data(A), idxs...), _elshape(A))
    setindex!(A_idxs, x)
end

function Base.similar(A::ShapedAsNTArray{T}, ::Type{T}, dims::Dims) where T
    data = _data(A)
    U = eltype(data)
    newdata = similar(data, U, dims)
    # In case newdata is not something like an ArrayOfSimilarVectors:
    if !isempty(newdata) && !isdefined(newdata, firstindex(newdata))
        for i in eachindex(newdata)
            newdata[i] = similar(data[firstindex(data)])
        end
    end
    ShapedAsNTArray(newdata, _elshape(A))
end

Base.empty(A::ShapedAsNTArray{T,N,D,S}) where {T,N,D,S} =
    ShapedAsNTArray{T,N,D,S}(empty(_data(A)), _elshape(A))

Base.show(io::IO, ::MIME"text/plain", A::ShapedAsNTArray) = TypedTables.showtable(io, A)

Base.copy(A::ShapedAsNTArray) = ShapedAsNTArray(copy(_data(A)), _elshape(A))

Base.pop!(A::ShapedAsNTArray) = _elshape(A)(pop!(_data(A)))

# Base.push!(A::ShapedAsNTArray, x::Any)  # ToDo

Base.popfirst!(A::ShapedAsNTArray) = _elshape(A)(popfirst!(_data(A)))

# Base.pushfirst!(A::ShapedAsNTArray, x::Any)  # ToDo

function Base.append!(A::ShapedAsNTArray, B::ShapedAsNTArray)
    _elshape(A) == _elshape(B) || throw(ArgumentError("Can't append ShapedAsNTArray instances with different element shapes"))
    append!(_data(A), _data(B))
    A
end

# Base.append!(A::ShapedAsNTArray, B::AbstractArray)  # ToDo

function Base.prepend!(A::ShapedAsNTArray, B::ShapedAsNTArray)
    _elshape(A) == _elshape(B) || throw(ArgumentError("Can't prepend ShapedAsNTArray instances with different element shapes"))
    prepend!(_data(A), _data(B))
    A
end

# Base.prepend!(A::ShapedAsNTArray, B::AbstractArray)  # ToDo

function Base.deleteat!(A::ShapedAsNTArray, i)
    deleteat!(_data(A), i)
    A
end

# Base.insert!(A::ShapedAsNTArray, i::Integer, x::Any)  # ToDo

Base.splice!(A::ShapedAsNTArray, i) = _elshape(A)(splice!(_data(A), i))

# Base.splice!(A::ShapedAsNTArray, i, replacement)  # ToDo

function Base.vcat(A::ShapedAsNTArray, B::ShapedAsNTArray)
    _elshape(A) == _elshape(B) || throw(ArgumentError("Can't vcat ShapedAsNTArray instances with different element shapes"))
    ShapedAsNTArray(vcat(_data(A), _data(B)), _elshape(A))
end

# Base.vcat(A::ShapedAsNTArray, B::AbstractArray)  # ToDo

# Base.hcat(A::ShapedAsNTArray, B) # ToDo

Base.vec(A::ShapedAsNTArray{T,1}) where T = A
Base.vec(A::ShapedAsNTArray) = ShapedAsNTArray(vec(_data(A)), _elshape(A))

Tables.istable(::Type{<:ShapedAsNTArray}) = true
Tables.rowaccess(::Type{<:ShapedAsNTArray}) = true
Tables.columnaccess(::Type{<:ShapedAsNTArray}) = true
Tables.schema(A::ShapedAsNTArray{T}) where {T} = Tables.Schema(T)

function Tables.columns(A::ShapedAsNTArray)
    data = _data(A)
    accessors = _accessors(_elshape(A))
    # Copy columns to make each column as contiguous in memory as possible:
    map(va -> getindex.(data, Ref(va)), accessors)
end

@inline Tables.rows(A::ShapedAsNTArray) = A


text:
using PyPlot
export plotModel

@doc doc"""
High-level plotting routine for model simulations

plotModel(ModelName, "simulationName")

Optional arguments:
- *t* _in stochastic simulations using Gillespie's algorithm XPP, the actual simulation time has to be calculated as an auxilliary variable_
- *pars* Boolean, _plot parameter overview?_
- *vars* Array, _which variable are plotted?_
- *xlim*
- *ylim*
- *linewidth*
- *alpha*
- *colors* Dict, _[varname => "HEXVALUE"] for custom colors_
- *fig* _default true, if true make new figure, else pass existing figure handle_
- *sp* _subplot handle, only in conjunction with figure = false_

"""->
function plotModel(M, simName; t = "t", pars = false, vars = false, xlim = false, ylim =false, colors = false, linewidth = 2, fig = true, alpha = 1, sp = false)
    if fig == true
        figure()
    end
    if pars
        ax1 = subplot2grid([4,1], [0,0], 3, 1)
        plotModel_core(ax1, M, simName; t = t, vars = vars, xlim = xlim, ylim = ylim, colors = colors, linewidth = linewidth, alpha = alpha)
        ax2 = subplot2grid([4,1], [3,0], 1, 1)
        plotPars(ax2, M, simName)
    else
        if sp == false
          ax = subplot(111)
        else
          ax = sp
        end
        plotModel_core(ax, M, simName; t = t, vars = vars, xlim = xlim, ylim = ylim, colors = colors, linewidth = linewidth, alpha = alpha)
    end
    if fig == true
        suptitle(simName)
    end
end

function plotPars(sp, M, simName)
    P = M.sims[simName].P
    pars = [P[p] for p in keys(P)]
    names = [p for p in keys(P)]
    ixs = [rand() for p in pars]
    scatter(pars, ixs, marker = "o",  color = "k", alpha = 0.5, label = names)
    sp[:set_xscale]("log")
    sp[:set_xlim]([0.0001, 10000])
    for i in [1:length(pars)]
        if pars[i] != 0
            annotate(names[i], xy = (pars[i]*1.2, ixs[i]))
        end
    end
    sp[:set_yticks]([])
    grid(b = true)
    grid(b = true, which = "minor", alpha = 0.2)
    # sp[:set_yticklabels](names)
end

function plotModel_core(subplot, M, simName; t = "t", vars = false, xlim = false, ylim =false, colors = false, linewidth = 2, alpha = 1)
    D = M.sims[simName].D
    if vars == false
        vars = M.vars[2:end]
    end
    for v in vars
        if colors == false
            subplot[:plot](D[t], D[v], linewidth = linewidth, alpha = alpha, label = v)
        else
            subplot[:plot](D[t], D[v], color = colors[v], linewidth = linewidth, alpha = alpha, label = v)
        end
    end
    subplot[:set_xlabel]("t (min)")
    subplot[:set_ylabel]("")
    if xlim != false
        subplot[:set_xlim](xlim)
    end
    if ylim != false
        subplot[:set_ylim](ylim)
    end
end

ext: jl

max_stars_repo_name: novakgroupoxford/XPPjl

max_issues_repo_name: novakgroupoxford/XPPjl

max_forks_repo_name: novakgroupoxford/XPPjl


text:
"""
    SpinDoctor

Diffusion MRI simulation toolbox.
"""
module SpinDoctor

using Arpack: eigs
using DiffEqCallbacks: PresetTimeCallback, FunctionCallingCallback
using Expokit: expmv!
using GLPK: Optimizer
using LinearAlgebra
using Makie
using MiniQhull: delaunay
using OrdinaryDiffEq: OrdinaryDiffEq
using OrdinaryDiffEq:
    DiffEqArrayOperator, DiscreteCallback, ODEProblem, ODEFunction, QNDF, MagnusGL6
using Polynomials: fit
using Printf
using QuadGK
using Roots: find_zeros
using SparseArrays
using SpecialFunctions
using StaticArrays
using Statistics: mean
using TetGen: RawTetGenIO, facetlist!, tetrahedralize
using Triangulate: TriangulateIO, triangulate
using WriteVTK

# Gradidents
export PGSE, CosOGSE, SinOGSE, DoublePGSE, GeneralGradient, ScalarGradient
export integral, int_F², intervals, isconstant, echotime

export assemble_matrices
export split_mesh
export split_field
export get_cmpt_volumes

export Model
export initial_conditions

export savefield

export compute_adc_sta
export fit_adc
export fit_tensors
export plot_mesh, plot_field

export IntervalConstantBTPDE, GeneralBTPDE, HADC, Karger
export Laplace, MatrixFormalism, AnalyticalLaplace, AnalyticalMatrixFormalism
export solve, solve_multigrad

# Matrix formalism
export length2eig, eig2length, limit_lengthscale, compute_mf_diffusion_tensor

# Utils
export unitcircle, unitsphere
export compute_signal

# Callbacks
export Printer, VTKWriter, Plotter

# Recipes
export AbstractSetup, PlateSetup, CylinderSetup, SphereSetup, NeuronSetup
export coefficients, analytical_coefficients
export create_geometry

# Utils
include("utils/fresnel.jl")
include("utils/unitcircle.jl")
include("utils/unitsphere.jl")
include("utils/compute_signal.jl")

# Magnetic field gradient sequences
include("gradients/sequences.jl")
include("gradients/gradient.jl")
include("gradients/integral.jl")
include("gradients/int_F2.jl")
include("gradients/intervals.jl")
include("gradients/echotime.jl")
include("gradients/isconstant.jl")
include("gradients/diffusion_time.jl")
include("gradients/diffusion_time_sta.jl")
include("gradients/j_integral.jl")

# Datatypes
include("datatypes/femesh.jl")
include("datatypes/model.jl")

# Problems
include("problems/problems.jl")
include("problems/output_type.jl")
include("problems/solve.jl")
include("problems/solve_multigrad.jl")

# Geometry
include("geometry/call_tetgen.jl")
include("geometry/convexhull.jl")
include("geometry/create_cells.jl")
include("geometry/create_fibonacci_sphere.jl")
include("geometry/deform_domain.jl")
include("geometry/get_volumes.jl")
include("geometry/get_mesh_surface.jl")
include("geometry/get_mesh_surfacenormals.jl")
include("geometry/gmesh2fem.jl")
include("geometry/read_cells.jl")
include("geometry/read_surfaces.jl")
include("geometry/read_tetgen.jl")
include("geometry/save_cells.jl")
include("geometry/save_surfaces.jl")
include("geometry/save_tetgen.jl")
include("geometry/split_mesh.jl")
include("geometry/split_field.jl")

# Matrix assembly
include("matrix_assembly/assemble_mass_matrix.jl")
include("matrix_assembly/assemble_stiffness_matrix.jl")
include("matrix_assembly/assemble_flux_matrices.jl")
include("matrix_assembly/assemble_flux_matrix.jl")
include("matrix_assembly/couple_flux_matrix.jl")
include("matrix_assembly/assemble_matrices.jl")

# Matrix formalism
include("matrix_formalism/eig2length.jl")
include("matrix_formalism/limit_lengthscale.jl")
include("matrix_formalism/solve_laplace.jl")
include("matrix_formalism/solve_mf.jl")
include("matrix_formalism/compute_mf_diffusion_tensor.jl")

# Analytical
include("analytical/alpha_func.jl")
include("analytical/compute_bc.jl")
include("analytical/compute_beta.jl")
include("analytical/compute_int_I.jl")
include("analytical/compute_int_J.jl")
include("analytical/compute_int_K.jl")
include("analytical/compute_JY.jl")
include("analytical/compute_nu.jl")
include("analytical/compute_v.jl")
include("analytical/find_alpha.jl")
include("analytical/solve_analytical_laplace.jl")
include("analytical/solve_analytical_mf.jl")

# BTPDE
include("btpde/initial_conditions.jl")
include("btpde/solve_btpde.jl")
include("btpde/solve_btpde_midpoint.jl")
include("btpde/solve_karger.jl")

# ADC
include("adc/compute_adc_sta.jl")
include("adc/solve_hadc.jl")

# Postprocessing
include("postprocess/fit_adc.jl")
include("postprocess/fit_tensors.jl")
include("postprocess/savefield.jl")

# Plot
include("plot/plot_mesh.jl")
include("plot/plot_field.jl")

# Callbacks
include("callbacks/callbacks.jl")
include("callbacks/initialize.jl")
include("callbacks/update.jl")
include("callbacks/finalize.jl")

# Recipes
include("recipes/setups.jl")
include("recipes/coefficients.jl")
include("recipes/create_surfaces.jl")
include("recipes/create_surfaces_plates.jl")
include("recipes/create_surfaces_cylinder.jl")
include("recipes/create_surfaces_sphere.jl")
include("recipes/create_surfaces_neuron.jl")
include("recipes/create_geometry.jl")
include("recipes/radial_dimension.jl")
include("recipes/analytical_coefficients.jl")

end


text:
adj = spzeros(6,6)
adj[1,3] = 1
adj[1,5] = 1
adj[2,3] = 1
adj[2,4] = 1
adj[2,5] = 1
adj[2,6] = 1
labels = [1, 2, 1, 2, 1, 2]
counts = RE.count_neighbor_groups(adj, labels)
@test counts[1,1] == 2
@test counts[2,1] == 0
@test counts[1,2] == 2
@test counts[2,2] == 2
@test counts[1,3] == 0
@test counts[2,3] == 0
@test counts[1,4] == 0
@test counts[2,4] == 0

pairs = [(i,j) for i in [1,2] for j in [1,2]]
ns, ms = RE.compile_count_dicts(counts, labels, pairs)
@test ms[1,1] == 2
@test ms[1,2] == 0
@test ms[2,1] == 2
@test ms[2,2] == 2
@test ns[1] == 3
@test ns[2] == 3

pr = [0.1 0.2; 0.3 0.4]
ll = calculate_ll(ns, ms, pr)
@test_approx_eq_eps ll -7.6713 10e-4

i = 1
count_deltas = counts[:,i][:]
old_label = 1
new_label = 2
new_ns, new_ms = RE.adjust_counts(ns, ms, count_deltas, old_label, new_label)
@test new_ms[1,1] == 0
@test new_ms[1,2] == 0
@test new_ms[2,1] == 4
@test new_ms[2,2] == 2
@test new_ns[1] == 2
@test new_ns[2] == 4

new_ll = RE.calculate_ll(new_ns, new_ms, pr)
@test_approx_eq_eps new_ll -7.57119 10e-4

new_pr = zeros(2,2)
RE.maximize_parameters!(ns,ms,new_pr)
@test new_pr[1,1] == 2/9
@test new_pr[1,2] == 0
@test new_pr[2,1] == 2/9
@test new_pr[2,2] == 2/9

# test bipartite
labels = [1, 2, 3, 4, 3, 4]
counts = RE.count_neighbor_groups(adj, labels)
@test counts[1,1] == 0
@test counts[2,1] == 0
@test counts[3,1] == 2
@test counts[4,1] == 0
@test counts[1,2] == 0
@test counts[2,2] == 0
@test counts[3,2] == 2
@test counts[4,2] == 2
@test counts[1,3] == 0
@test counts[2,3] == 0
@test counts[3,3] == 0
@test counts[4,3] == 0
@test counts[1,4] == 0
@test counts[2,4] == 0
@test counts[3,4] == 0
@test counts[4,4] == 0
@test counts[1,5] == 0
@test counts[2,5] == 0
@test counts[3,5] == 0
@test counts[4,5] == 0
@test counts[1,6] == 0
@test counts[2,6] == 0
@test counts[3,6] == 0
@test counts[4,6] == 0
pairs = [(i,j) for i in [1,2] for j in [3,4]]
ns, ms = RE.compile_count_dicts(counts, labels, pairs)
@test ms[1,3] == 2
@test ms[1,4] == 0
@test ms[2,3] == 2
@test ms[2,4] == 2
@test ns[1] == 1
@test ns[2] == 1
@test ns[3] == 2
@test ns[4] == 2
@test RE.swap_group(1) == 2
@test RE.swap_group(2) == 1
@test RE.swap_group(3) == 4
@test RE.swap_group(4) == 3


text:
"""
    struct CartesianSpace <: EuclideanSpace{ℝ}

A real Euclidean space `ℝ^d`, which is therefore self-dual. `CartesianSpace` has no
additonal structure and is completely characterised by its dimension `d`. This is the
vector space that is implicitly assumed in most of matrix algebra.
"""
struct CartesianSpace <: EuclideanSpace{ℝ}
    d::Int
end
CartesianSpace(d::Integer = 0; dual = false) = CartesianSpace(Int(d))
function CartesianSpace(dim::Pair; dual = false)
    if dim.first === Trivial()
        return CartesianSpace(dim.second; dual = dual)
    else
        msg = "$(dim) is not a valid dimension for CartesianSpace"
        throw(SectorMismatch(msg))
    end
end
function CartesianSpace(dims::AbstractDict; kwargs...)
    if length(dims) == 0
        return CartesianSpace(0; kwargs...)
    elseif length(dims) == 1
        return CartesianSpace(first(dims); kwargs...)
    else
        msg = "$(dims) is not a valid dimension dictionary for CartesianSpace"
        throw(SectorMismatch(msg))
    end
end

# convenience constructor
Base.getindex(::RealNumbers) = CartesianSpace
Base.:^(::RealNumbers, d::Int) = CartesianSpace(d)

# Corresponding methods:
#------------------------
dim(V::CartesianSpace) = V.d
Base.axes(V::CartesianSpace) = Base.OneTo(dim(V))

Base.oneunit(::Type{CartesianSpace}) = CartesianSpace(1)
⊕(V1::CartesianSpace, V2::CartesianSpace) = CartesianSpace(V1.d+V2.d)
fuse(V1::CartesianSpace, V2::CartesianSpace) = CartesianSpace(V1.d*V2.d)
flip(V::CartesianSpace) = V

infimum(V1::CartesianSpace, V2::CartesianSpace) = CartesianSpace(min(V1.d, V2.d))
supremum(V1::CartesianSpace, V2::CartesianSpace) = CartesianSpace(max(V1.d, V2.d))

Base.show(io::IO, V::CartesianSpace) = print(io, "ℝ^$(V.d)")


text:
"""
App level functionality -- loading and managing app-wide components like configs, models, initializers, etc.
"""
module App

using Revise
using YAML
using Genie

const ASSET_FINGERPRINT = ""

function bootstrap()
  if haskey(ENV, "GENIE_ENV") && isfile(joinpath(Genie.ENV_PATH, ENV["GENIE_ENV"] * ".jl"))
    isfile(joinpath(Genie.CONFIG_PATH, "global.jl")) && include(joinpath(Genie.CONFIG_PATH, "global.jl"))
    include(joinpath(Genie.ENV_PATH, ENV["GENIE_ENV"] * ".jl"))
  else
    ENV["GENIE_ENV"] = Configuration.DEV
    eval(@__MODULE__, Meta.parse("config = Configuration.Settings(app_env = Configuration.DEV)"))
  end
end

end

### Main

using Genie.Loggers, Genie.Configuration

"""
    load_libs(dir = Genie.LIB_PATH) :: Nothing

Recursively adds subfolders of lib to LOAD_PATH.
"""
function load_libs(root_dir = Genie.LIB_PATH) :: Nothing
  push!(LOAD_PATH, root_dir)
  for (root, dirs, files) in walkdir(root_dir)
    for dir in dirs
      p = joinpath(root, dir)
      in(p, LOAD_PATH) || push!(LOAD_PATH, p)
    end
  end

  nothing
end

"""
    load_resources(dir = Genie.RESOURCES_PATH) :: Nothing

Recursively adds subfolders of resources to LOAD_PATH.
"""
function load_resources(root_dir = Genie.RESOURCES_PATH) :: Nothing
  push!(LOAD_PATH, root_dir)

  for (root, dirs, files) in walkdir(root_dir)
    for dir in dirs
      p = joinpath(root, dir)
      in(p, LOAD_PATH) || push!(LOAD_PATH, joinpath(root, dir))
    end
  end

  nothing
end

function load_helpers(root_dir = Genie.HELPERS_PATH) :: Nothing
  push!(LOAD_PATH, root_dir)

  for (root, dirs, files) in walkdir(root_dir)
    for dir in dirs
      p = joinpath(root, dir)
      in(p, LOAD_PATH) || push!(LOAD_PATH, joinpath(root, dir))
    end
  end

  nothing
end

"""
    load_acl(dir::String) :: Dict{Any,Any}

Loads the ACL file associated with the invoked `controller` and returns the rules.
"""
function load_acl(dir::String) :: Dict{Any,Any}
  file_path = joinpath(dir, Genie.GENIE_AUTHORIZATOR_FILE_NAME)
  isfile(file_path) ? YAML.load(open(file_path)) : Dict{Any,Any}()
end

"""
    load_configurations() :: Nothing

Loads (includes) the framework's configuration files.
"""
function load_configurations() :: Nothing
  loggers_path = abspath("$(Genie.CONFIG_PATH)/loggers.jl")
  if isfile(loggers_path)
    include(loggers_path)
    is_dev() && Revise.track(loggers_path)
  end

  secrets_path = abspath("$(Genie.CONFIG_PATH)/secrets.jl")
  if isfile(secrets_path)
    include(secrets_path)
    is_dev() && Revise.track(secrets_path)
  end

  nothing
end

"""
    load_initializers() :: Nothing

Loads (includes) the framework's initializers.
"""
function load_initializers() :: Nothing
  dir = abspath(joinpath(Genie.CONFIG_PATH, "initializers"))

  if isdir(dir)
    f = readdir(dir)
    for i in f
      fi = joinpath(dir, i)
      if endswith(fi, ".jl")
        include(fi)
        is_dev() && Revise.track(fi)
      end
    end
  end

  nothing
end

"""
    load_routes_definitions() :: Nothing

Loads the routes file.
"""
function load_routes_definitions(fail_on_error = is_dev()) :: Nothing
  try
    if isfile(Genie.ROUTES_FILE_NAME)
      include(Genie.ROUTES_FILE_NAME)
      is_dev() && Revise.track(Genie.ROUTES_FILE_NAME)
    end
  catch ex
    log(ex, :warn)

    fail_on_error && rethrow(ex)
  end

  nothing
end

"""
    load_channels_definitions() :: Nothing

Loads the channels file.
"""
function load_channels_definitions(fail_on_error = is_dev()) :: Nothing
  try
    if isfile(Genie.CHANNELS_FILE_NAME)
      include(Genie.CHANNELS_FILE_NAME)
      is_dev() && Revise.track(Genie.CHANNELS_FILE_NAME)
    end
  catch ex
    log(ex, :err)

    fail_on_error && rethrow(ex)
  end

  nothing
end

"""
    secret_token() :: String

Wrapper around /config/secrets.jl SECRET_TOKEN `const`.
"""
function secret_token() :: String
  if @isdefined SECRET_TOKEN
    SECRET_TOKEN
  else
    log("SECRET_TOKEN not configured - please make sure that you have a valid secrets.jl file.
          You can generate a new secrets.jl file with a random SECRET_TOKEN using Genie.REPL.write_secrets_file()
          or use the included /app/config/secrets.jl.example file as a model.", :warn)
    st = Genie.REPL.secret_token()
    Core.eval(Genie, Meta.parse("""const SECRET_TOKEN = "$st" """))

    st
  end
end

function load() :: Nothing
  App.bootstrap()

  load_configurations()

  Genie.Loggers.log_path!()
  Genie.Loggers.empty_log_queue()

  load_initializers()
  load_helpers()

  load_libs()
  load_resources()

  load_routes_definitions()
  load_channels_definitions()

  nothing
end


text:
function molecule(::Molecule"methylene")
  return """
    C            0.000000000000     0.000000000000    -0.062823577797
    H            0.000000000000    -0.976631669800     0.374014789003
    H           -0.000000000000     0.976631669800     0.374014789003
  """
end


text:
###############################################################################
#
# Soil optical parameters
#
###############################################################################
"""
    mutable struct SoilOpticals{FT}

A struct of soil optical parameters

# Fields
$(TYPEDFIELDS)
"""
mutable struct SoilOpticals{FT}
    "Soil surface temperature"
    T::FT
    "Soil color class"
    color::Int

    # broad band showrtwave albedo
    "Shortwave albedo for NIR"
    ρ_NIR::FT
    "Shortwave albedo for PAR"
    ρ_PAR::FT

    # hyperspectral shortwave albedo
    "Shortwave albedo that matches `WL` from [`Wavelengths`](@ref)"
    ρ_SW::Vector{FT}
    "Shortwave albedo that matches `WLF` from [`Wavelengths`](@ref)"
    ρ_SW_SIF::Vector{FT}
    "Shortwave absorption that equals `1 - ρ_SW`"
    ε_SW::Vector{FT}

    # hyperspectral soil albedo
    "Shortwave albedo matrix from 4 bands, wavelengths are 400:10:2500 nm"
    SW_mat_raw_4::Matrix{FT}
    "Shortwave albedo matrix from 2 bands, wavelengths are 400:10:2500 nm"
    SW_mat_raw_2::Matrix{FT}
    "Shortwave albedo matrix from 4 bands with `WL` from [`Wavelengths`](@ref)"
    SW_mat_4::Matrix{FT}
    "Shortwave albedo matrix from 2 bands with `WL` from [`Wavelengths`](@ref)"
    SW_mat_2::Matrix{FT}
    "Shortwave albedo weight from 4 bands"
    SW_vec_4::Vector{FT}
    "Shortwave albedo weight from 2 bands"
    SW_vec_2::Vector{FT}
    "Shortwave albedo, wavelengths are 400:10:2500 nm"
    ρ_SW_raw::Vector{FT}

    # hyperspectral longwave albedo
    "Longtwave albedo"
    ρ_LW::Vector{FT}

    # cache that stores mean band values (used to speed up calculations)
    "Mean value for day band 1 in NIR region"
    dry_NIR::FT
    "Mean value for day band 1 in PAR region"
    dry_PAR::FT
    "Mean value for day band 1 in NIR region"
    wet_NIR::FT
    "Mean value for day band 1 in PAR region"
    wet_PAR::FT

    # constructor
    function SoilOpticals{FT}(wls::WaveLengths{FT}) where {FT<:AbstractFloat}
        @unpack nWL, nWLF, sWL = wls;

        # read data that has a 10 nm stepping
        _dat = read_csv(SOIL_GSV);
        _raw_4 = Matrix{FT}(_dat[:,2:end-1]);

        # extend the data to 1 nm stepping by interpolating the data
        _wlr = _dat.WL;
        _wle = collect(400:2500.1);
        _ext_4 = Matrix{FT}(undef, length(_wle), 4);
        for _ie in 1:size(_ext_4,1)-1
            _wl= _wle[_ie];
            _ir = Int(fld(_wl - 400, 10)) + 1;
            _a = (_wl - _wlr[_ir]) * 0.1;
            _ext_4[_ie,:] .= (1-_a) .* _raw_4[_ir,:] .+ _a .* _raw_4[_ir+1,:];
        end
        _ext_4[end,:] = _raw_4[end,:];

        # rescale the data to match the steppings of wavelength set
        _res_4 = Matrix{FT}(undef, nWL, 4);

        # fill in the arrays
        for _i_res in 1:nWL
            _wo = findall( (_wle .>= sWL[_i_res]) .& (_wle .< sWL[_i_res+1]) )
            if length(_wo) == 0
                @warn "Some wavelengths out of bounds $(string(sWL[_i_res]))";
            end;
            _res_4[_i_res,1] = mean( _ext_4[_wo,1] );
            _res_4[_i_res,2] = mean( _ext_4[_wo,2] );
            _res_4[_i_res,3] = mean( _ext_4[_wo,3] );
            _res_4[_i_res,4] = mean( _ext_4[_wo,4] );
        end

        # 2 band values
        _raw_2 = [_raw_4[:,1] _raw_4[:,end]];
        _res_2 = [_res_4[:,1] _res_4[:,end]];

        # mean values for 2 bands
        _dry_nir = mean( _raw_2[31:end,1] );
        _dry_par = mean( _raw_2[1:31,1]   );
        _wet_nir = mean( _raw_2[31:end,2] );
        _wet_par = mean( _raw_2[1:31,2]   );

        return new(T_25(FT), 1, FT(0.2), FT(0.2), ones(FT,nWL)*FT(0.2),
                   ones(FT,nWLF)*FT(0.2), ones(FT,nWL)*FT(0.8), _raw_4, _raw_2,
                   _res_4, _res_2, ones(FT,4), ones(FT,2), ones(FT,211),
                   FT[0.1], _dry_nir, _dry_par, _wet_nir, _wet_par)
    end
end


text:
function raw_to_df(raw_data)

    # Argument ":auto" is required to generate column names in latest version of DataFrames.jl (v1.2.1)
    df = DataFrame(raw_data[1], :auto)
    df_names = Symbol.(vcat(raw_data[2]...))
    df = DataFrames.rename(df, df_names)

    timestamps = df[!,:timestamp]

    select!(df, Not([:timestamp]))

    for col in eachcol(df)
        col = Float64.(col)
    end

    df[!,:Date] = Date.(timestamps)
    return df
end

function average_price_df(currency::String, df_in::DataFrame)

    df_out_price, df_out_candle  = [DataFrame() for i = 1:2]

    df_out_price[!,:Date] = df_in[!,:Date]

    df_out_price[!,Symbol("$currency")] = (df_in[!,Symbol("open (EUR)")] + df_in[!,Symbol("high (EUR)")] + df_in[!,Symbol("low (EUR)")] + df_in[!,Symbol("close (EUR)")])/4

    candle_col = Any[]
    for i in 1:size(df_in)[1]
        push!(candle_col, (df_in[!,Symbol("open (EUR)")][i], df_in[!,Symbol("high (EUR)")][i], df_in[!,Symbol("low (EUR)")][i], df_in[!,Symbol("close (EUR)")][i]))
    end

    df_out_candle[!,:Date] = df_in[!,:Date]

    df_out_candle[!,Symbol("$currency")] = candle_col

    return df_out_price, df_out_candle
end

function vol_df(currency::String, df_in::DataFrame)

    df_out_vol = DataFrame()

    df_out_vol[!,:Date] = df_in[!,:Date]
    df_out_vol[!,Symbol("$currency")] = df_in[!,:volume]

    return df_out_vol
end

function moving_averages(Price_df::DataFrame, duration::Int64, window::Int64)

    # Price_df should have date order - oldest to latest
    Price_col = Price_df[end-duration+1-window+1:end,2]
    rows1 = length(Price_col)
    Price_SMA, Price_WMA, Price_EMA = [Float64[] for i = 1:3]

    weights = collect(1:window)/(window*(window+1)/2)
    k = 2/(window+1)

    # Calculate different Moving Averages
    for i = 1:rows1-(window-1)
        # Simple Moving Average (SMA)
        push!(Price_SMA, mean(Price_col[i:i+(window-1)]))

        # Weighted Moving Average
        push!(Price_WMA, sum(Price_col[i:i+(window-1)].*weights))

        # Exponential Moving Average
        SMA = mean(Price_col[i:i+(window-1)])
        EMA = (Price_col[i+(window-1)]*k)+(SMA*(1-k))
        push!(Price_EMA, EMA)
    end

    return Price_SMA, Price_WMA, Price_EMA
end

function calculate_ema(Price_col::Vector{Float64}, window::Int64)

    k(window) = 2/(window+1)
	Price_EMA = Float64[]
	rows = length(Price_col)

	for i = 1:rows-(window-1)

		SMA = mean(Price_col[i:i+(window-1)])
        EMA = (Price_col[i+(window-1)]*k(window)) +
		      (SMA*(1-k(window)))
        push!(Price_EMA, EMA)
	end

	return Price_EMA
end

function calculate_macd(Price_df::DataFrame, window_long::Int64 = 26,
                        window_short::Int64 = 12, window_signal::Int64 = 9)

	# Price_df should have date order - oldest to latest
	Price_col = Price_df[!, 2]

	# Calculate 26 (long) and 12 (short) period EMA
	Price_EMA_short = calculate_ema(Price_col, window_short)
	Price_EMA_long = calculate_ema(Price_col, window_long)

	# Make EMA_long and EMA_short equal
	EMA_short_col = Price_EMA_short[window_long - window_short + 1:end]

	# Calculate MACD = EMA_short - EMA_long
	MACD_col = EMA_short_col - Price_EMA_long

	# Calculate signal line (9 period EMA of MACD)
	Signal_col = calculate_ema(MACD_col, window_signal)

	df_ema = DataFrame(date = Price_df[window_long+window_signal-1:end, :Date],
                      Raw = Price_df[window_long+window_signal-1:end, 2],
                      EMA_long = Price_EMA_long[window_signal:end],
                      EMA_short = EMA_short_col[window_signal:end],
                      MACD = MACD_col[window_signal:end], Signal = Signal_col)

	return df_ema
end


text:
abstract type AbstractVisualizer end

# Interface methods
function setobject! end
function settransform! end
function delete! end
function setprop! end
function setcontrol! end

# Convenient shortcuts
setobject!(vis::AbstractVisualizer, geom::GeometryLike) = setobject!(vis, Object(geom))
setobject!(vis::AbstractVisualizer, geom::GeometryLike, material::AbstractMaterial) = setobject!(vis, Object(geom, material))
setvisible!(vis::AbstractVisualizer, visible::Bool) = setprop!(vis, "visible", visible)


text:
@testset "SIP Res" begin
    # Define semi-infinite program
    f(x) = (1/3)*x[1]^2 + x[2]^2 + x[1]/2
    gSIP(x,p) = (1.0 - (x[1]^2)*(p[1]^2))^2 - x[1]*p[1]^2 - x[2]^2 + x[2]

    x_l = Float64[-1000.0, -1000.0]
    x_u = Float64[1000.0, 1000.0]
    p_l = Float64[0.0]
    p_u = Float64[1.0]

    sip_result = sip_solve(SIPRes(), x_l, x_u, p_l, p_u, f, Any[gSIP], abs_tolerance = 1E-3)

    @test isapprox(sip_result.lower_bound, 0.19446619886176916, atol = 1E-3)
    @test isapprox(sip_result.upper_bound, 0.19500611503848964, atol = 1E-3)
    @test isapprox(sip_result.xsol[1], -0.7500000115038946, atol = 1E-2)
    @test isapprox(sip_result.xsol[2], -0.6184706298867955, atol = 1E-2)
end
#=
@testset "SIP ResRev" begin
    # Define semi-infinite program
    f(x) = (1/3)*x[1]^2 + x[2]^2 + x[1]/2
    gSIP(x,p) = (1.0 - (x[1]^2)*(p[1]^2))^2 - x[1]*p[1]^2 - x[2]^2 + x[2]

    x_l = Float64[-1000.0, -1000.0]
    x_u = Float64[1000.0, 1000.0]
    p_l = Float64[0.0]
    p_u = Float64[1.0]

    sip_result = sip_solve(SIPResRev(), x_l, x_u, p_l, p_u, f, Any[gSIP], abs_tolerance = 1E-3)

    @test isapprox(sip_result.lower_bound, 0.19446619886176916, atol = 1E-3)
    @test isapprox(sip_result.upper_bound, 0.19500611503848964, atol = 1E-3)
    @test isapprox(sip_result.xsol[1], -0.7500000115038946, atol = 1E-2)
    @test isapprox(sip_result.xsol[2], -0.6184706298867955, atol = 1E-2)
end
=#
@testset "SIP Hybrid" begin
    # Define semi-infinite program
    f(x) = (1/3)*x[1]^2 + x[2]^2 + x[1]/2
    gSIP(x,p) = (1.0 - (x[1]^2)*(p[1]^2))^2 - x[1]*p[1]^2 - x[2]^2 + x[2]

    x_l = Float64[-1000.0, -1000.0]
    x_u = Float64[1000.0, 1000.0]
    p_l = Float64[0.0]
    p_u = Float64[1.0]

    sip_result = sip_solve(SIPHybrid(), x_l, x_u, p_l, p_u, f, Any[gSIP], abs_tolerance = 1E-3)

    @test isapprox(sip_result.lower_bound, 0.19446619886176916, atol = 1E-3)
    @test isapprox(sip_result.upper_bound, 0.19500611503848964, atol = 1E-3)
    @test isapprox(sip_result.xsol[1], -0.7500000115038946, atol = 1E-2)
    @test isapprox(sip_result.xsol[2], -0.6184706298867955, atol = 1E-2)
end

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
module Day24

export get_inputs, get_solution1, get_solution2

using StaticArrays, OffsetArrays

## Input getting
function get_inputs()
    test_input1 = test_input2 = readlines(joinpath(@__DIR__, "test_input1.txt"))
    test_output1 = 10
    test_output2 = 2208
    data = readlines(joinpath(@__DIR__, "input.txt"))
    return (; test_input1, test_input2, test_output1, test_output2, data)
end

const DIRECTION_MAP = Dict(
    "ne" => ( 0,  1),
    "nw" => (-1,  1),
    "sw" => ( 0, -1),
    "se" => ( 1, -1),
    "w"  => (-1,  0),
    "e"  => ( 1,  0),
)

function parse_directions(directions)
    directions = getindex.(directions, findall(r"(ne|nw|sw|se|w|e)", directions))
    return [DIRECTION_MAP[direction] for direction in directions]
end

## Solution functions
const BLACK = true
const WHITE = false

function initialize_tiles(directions, n)
    flip_tiles = [(sum(first.(direction)), sum(last.(direction)))
        for direction in parse_directions.(directions)
    ]
    # flip_tiles = sum.(parse_directions.(directions))
    black_tiles = Set(flip_tiles)
    positions = [tile for tile in black_tiles if isodd(count(==(tile), flip_tiles))]
    bounding_box = ((:)(extrema(first(positions))...), (:)(extrema(first(positions))...))
    bounding_box = grow(bounding_box, n+1)
    tiles = zeros(Bool, bounding_box...)
    tiles[CartesianIndex.(positions)] .= true
    return (; tiles, bounding_box)
end

grow(bounding_box, n=1) = map(range -> first(range)-n:last(range)+n, bounding_box)

# Part 1
get_solution1(data) = initialize_tiles(data, 10).tiles |> count

# Part 2
function get_solution2(data)
    tile_colors, bounding_box = initialize_tiles(data, 100)
    temp = copy(tile_colors)
    adj_inds = values(DIRECTION_MAP)

    for turn in 1:100
        temp .= tile_colors

        for pos in CartesianIndices(tile_colors)
            count_black = 0
            for adj_pos in (Tuple(pos) .+ i for i in adj_inds)
                !checkbounds(Bool, tile_colors, adj_pos...) && break
                if tile_colors[adj_pos...] == BLACK
                    count_black += 1
                end
            end

            color = tile_colors[pos]
            if color==BLACK && (count_black==0 || count_black>2)
                temp[pos] = WHITE
            elseif color==WHITE && count_black==2
                temp[pos] = BLACK
            end
        end

        tile_colors .= temp
    end

    return count(values(tile_colors))
end

end


text:
using Documenter
using PooksoftAlphaVantageDataStore

makedocs(sitename="Pooksoft",
    pages = [
        "index.md",

        "Utility" =>[
            "utility.md"
        ],

        "Stock Time Series (STS)" => [
            "execute_sts_daily_api_call.md",
            "execute_sts_weekly_api_call.md",
            "execute_sts_adjusted_weekly_api_call.md",
            "execute_sts_monthly_api_call.md",
            "execute_sts_adjusted_monthly_api_call.md",
            "execute_quote_api_call.md",
            "execute_search_api_call.md"
        ],

        "Technical Indicators (TI)" => [
            "execute_ti_simple_moving_average_api_call.md",
            "execute_ti_exponential_moving_average_api_call.md",
            "execute_ti_rsi_api_call.md"
        ]
    ]
)

# Deploy -
# Documenter can also automatically deploy documentation to gh-pages.
# See "Hosting Documentation" and deploydocs() in the Documenter manual
# for more information.
deploydocs(
    repo = "github.com/Pooksoft/PooksoftAlphaVantageDataStore.jl.git",
    devbranch = "master",
    devurl = "dev",
)


text:
# Quick Structs for checking calculations
import Impero:compute
struct Wrapper{T, S} <: AbstractExpression
    data::T
    meta_data::S
end
# Struct for MetaData
struct WrapperMetaData{T}
    io_name::T
end

function Base.show(io::IO, field::Wrapper{T, S}) where {T <: Char, S}
    color = 230
    printstyled(io, field.data, color = color)
end
function Base.show(io::IO, field::Wrapper{T, S}) where {T, S <: WrapperMetaData}
    color = 230
    printstyled(io, field.meta_data.io_name, color = color)
end

compute(a::Wrapper) = a.data

macro wrapper(expr)
    rewritten_expr = _wrapper(expr)
    return rewritten_expr
end

function _wrapper(expr::Expr)
    symb = expr.args[1]
    val  = expr.args[2]
    string_symb = String(symb)
    new_expr = :($(esc(symb)) =  Wrapper($val, WrapperMetaData($string_symb)))
    return new_expr
end

macro wrapper(exprs...)
    rewritten_exprs = [_wrapper(expr) for expr in exprs]
    return Expr(:block, rewritten_exprs...)
end

## Add directional derivatives
struct DirectionalDerivative{𝒮} <: AbstractExpression
    direction::𝒮
end
struct GradientMetaData{𝒮}
    direction::𝒮
end
function (p::DirectionalDerivative)(expr::AbstractExpression)
    return Gradient(expr, GradientMetaData(p.direction))
end
function Base.show(io::IO, p::DirectionalDerivative{S}) where S <: String
    print(io, Char(0x02202) * p.direction)
end
function Base.show(io::IO, p::Gradient{S, T}) where {S, T <: GradientMetaData{String}}
    printstyled(io, Char(0x02202) * p.metadata.direction, "(", color = 165)
    print(io, p.operand)
    printstyled(io, ")", color = 165)
end

∂x = DirectionalDerivative("x")
∂y = DirectionalDerivative("y")
∂z = DirectionalDerivative("z")
∂t = DirectionalDerivative("t")


text:
function set_bond_values_reactive(;
    session::ServerSession, notebook::Notebook,
    bound_sym_names::AbstractVector{Symbol},
    is_first_values::AbstractVector{Bool}=[true for x in bound_sym_names],
    initiator=nothing,
    kwargs...
)::Union{Task,TopologicalOrder}
    # filter out the bonds that don't need to be set
    to_set = filter(zip(bound_sym_names, is_first_values) |> collect) do (bound_sym, is_first_value)
        new_value = notebook.bonds[bound_sym].value

        variable_exists = is_assigned_anywhere(notebook, notebook.topology, bound_sym)
        if !variable_exists
            # a bond was set while the cell is in limbo state
            # we don't need to do anything
            return false
        end

        # TODO: Not checking for any dependents now
        # any_dependents = is_referenced_anywhere(notebook, notebook.topology, bound_sym)

        # fix for https://github.com/fonsp/Pluto.jl/issues/275
        # if `Base.get` was defined to give an initial value (read more about this in the Interactivity sample notebook), then we want to skip the first value sent back from the bond. (if `Base.get` was not defined, then the variable has value `missing`)
        # Check if the variable does not already have that value.
        # because if the initial value is already set, then we don't want to run dependent cells again.
        eq_tester = :(try !ismissing($bound_sym) && ($bound_sym == Main.PlutoRunner.transform_bond_value($(QuoteNode(bound_sym)), $(new_value))) === true catch; false end) # not just a === comparison because JS might send back the same value but with a different type (Float64 becomes Int64 in JS when it's an integer. The `=== true` check handles cases like `[missing] == [123]`, which returns `missing`, not `true` or `false`.)
        if is_first_value && will_run_code(notebook) && WorkspaceManager.eval_fetch_in_workspace((session, notebook), eq_tester)
            return false
        end
        return true
    end .|> first

    if isempty(to_set) || !will_run_code(notebook)
        send_notebook_changes!(ClientRequest(; session, notebook, initiator))
        return TopologicalOrder(notebook.topology, Cell[], Dict{Cell, ReactivityError}())
    end

    new_values = [notebook.bonds[bound_sym].value for bound_sym in to_set]

    function custom_deletion_hook((session, notebook)::Tuple{ServerSession,Notebook}, old_workspace_name, new_workspace_name, to_delete_vars::Set{Symbol}, methods_to_delete::Set{Tuple{UUID,FunctionName}}, to_reimport::Set{Expr}; to_run::AbstractVector{Cell})
        to_delete_vars = Set([to_delete_vars..., to_set...]) # also delete the bound symbols
        WorkspaceManager.move_vars((session, notebook), old_workspace_name, new_workspace_name, to_delete_vars, methods_to_delete, to_reimport)
        for (bound_sym, new_value) in zip(to_set, new_values)
            WorkspaceManager.eval_in_workspace((session, notebook), :($(bound_sym) = Main.PlutoRunner.transform_bond_value($(QuoteNode(bound_sym)), $(new_value))))
        end
    end
    to_reeval = where_referenced(notebook, notebook.topology, Set{Symbol}(to_set))

    run_reactive_async!(session, notebook, to_reeval; deletion_hook=custom_deletion_hook, user_requested_run=false, run_async=false, kwargs...)
end

"""
Returns the set of all possible values for the binded variable `n` as returned by the widget implementation using `AbstractPlutoDingetjes.possible_bond_values(element)`. This API is meant to be used by PlutoSliderServer.
"""
possible_bond_values(session::ServerSession, notebook::Notebook, name::Symbol) = WorkspaceManager.possible_bond_values((session,notebook), name)

"""
Optimized version of `length ∘ possible_bond_values`.
"""
possible_bond_values_length(session::ServerSession, notebook::Notebook, name::Symbol) = WorkspaceManager.possible_bond_values((session,notebook), name; get_length=true)


text:
using Test, YAAD, YAAD.TestUtils

include("math.jl")
include("linalg.jl")
include("cat.jl")
include("array.jl")


text:
module SimplexQuad

using LinearAlgebra

export simplexquad

"""
    simplexquad([::Type{T},] N::Int, n::Int)
    simplexquad(N::Int, vert::AbstractMatrix)

Construct Gauss points and weights for an `n`-dimensional simplex
domain with vertices specified by the `n*(n-1)` matrix `vert`, where
each row contains the coordinates `(x1,...,xn)` for a vertex. The
order of the quadrature scheme in each dimension must be the same in
this implementation.

# Arguments
- `N`: Number of quadrature points
- `T`: Return type (defaults to Float64)
- `n`: Number of dimensions
- `vert`: matrix of simplex vertex coordinates

# Result
- `X`: Matrix of quadrature points, where each row contains the
  coordinates for a point
- `W`: Array of quadrature weights

# Sample Usage:

    X, W = simplexquad(n, vert)     # Specify the vertices
    X, W = simplexquad(T, n, dim)   # Specify the dimension and use unit simplex

`X` will be a matrix for which the `j`th column are the grid points in
each coordinate `xj`.

To integrate a function `f`, use e.g.

    sum(W[i] * f(X[i,:]) for i in 1:length(W))

Note: The standard `n`-dimensional simplex has vertices specified
`vert = diagm(n+1, n, ones(n))
"""
simplexquad

function simplexquad(::Type{T}, N::Int, n::Int) where {T}
    m = n + 1
    vert = diagm(m, n, ones(T, n))
    return simplexquad(N, vert)
end
simplexquad(N::Int, n::Int) = simplexquad(float(Int), N, n)

function simplexquad(N::Int, vert::AbstractMatrix{T}) where {T}
    N <= 0 && @error "Number of Gauss points must be a natural number"

    m, n = size(vert)
    m != n + 1 &&
        @error "The matrix of vertices must have n+1 rows and n columns"

    Nn = N^n
    if n == 1
        # The 1-D simplex is only an interval
        q, w = rquad(T, N, 0)
        len = diff(vert; dims = 1)
        X = vert[1] .+ len[1] * q
        W = abs.(len[1]) * w
        X = reshape(X, :, 1)
    else
        # Find quadrature rules for higher dimensional domains
        q = Array{Vector}(undef, n)
        w = Array{Vector}(undef, n)
        for k in 1:n
            q[k], w[k] = rquad(T, N, n - k)
        end
        Q = ndgrid(q...)
        q = reshape(cat(Q...; dims = n), Nn, n)
        W = ndgrid(w...)
        w = reshape(cat(W...; dims = n), Nn, n)
        map = zeros(T, m, m) + I
        map[2:m, 1] .= -1
        c = map * vert
        W = abs(det(c[2:m, :])) * prod(w; dims = 2)
        qp = cumprod(q; dims = 2)
        e = ones(T, Nn, 1)
        X = [e [(1 .- q[:, 1:(n - 1)]) e] .* [e qp[:, 1:(n - 2)] qp[:, n]]] * c
        @assert size(W, 2) == 1
        W = reshape(W, :)
    end

    X::Matrix{T}
    @assert size(X, 2) == n
    W::Vector{T}
    @assert length(W) == size(X, 1)
    return X, W
end

function rquad(::Type{T}, N::Int, k::Int) where {T}
    k1 = k + 1
    k2 = k + 2

    n = collect(1:N)'
    nnk = 2 * n .+ k
    A = [T(k) / k2 repeat([T(k^2)], 1, N) ./ (nnk .* (nnk .+ 2))]

    n = collect(2:N)'
    nnk = nnk[n]
    B1 = T(4) * k1 / (k2 * k2 * (k + 3))

    nk = n .+ k
    nnk2 = nnk .* nnk
    B = T(4) * (n .* nk) .^ 2 ./ (nnk2 .* nnk2 - nnk2)

    ab = [A' [T(2^k1) / k1; B1; B']]
    s = sqrt.(ab[2:N, 2])
    X, V = eigen(SymTridiagonal(ab[1:N, 1], s))
    I = sortperm(X)
    X = X[I]
    x = (X .+ 1) / 2
    w = (T(1) / 2)^k1 * ab[1, 2] * V[1, I] .^ 2

    x::Vector{T}
    @assert length(x) == N
    w::Vector{T}
    @assert length(w) == N
    return x, w
end

# These functions were formerly a part of Julia. License is MIT:
# https://julialang.org/license

ndgrid(v::AbstractVector) = copy(v)

function ndgrid(v1::AbstractVector{T}, v2::AbstractVector{T}) where {T}
    m, n = length(v1), length(v2)
    v1 = reshape(v1, m, 1)
    v2 = reshape(v2, 1, n)
    return (repeat(v1, 1, n), repeat(v2, m, 1))
end

function ndgrid_fill(a, v, s, snext)
    return for j in 1:length(a)
        a[j] = v[div(rem(j - 1, snext), s) + 1]
    end
end

function ndgrid(vs::AbstractVector{T}...) where {T}
    n = length(vs)
    sz = map(length, vs)
    out = ntuple(i -> Array{T}(undef, sz), n)
    s = 1
    for i in 1:n
        a = out[i]::Array
        v = vs[i]
        snext = s * size(a, i)
        ndgrid_fill(a, v, s, snext)
        s = snext
    end
    return out
end

integrate(f, X, W) = sum(W[i] * f(X[i, :]) for i in 1:length(W))

end


text:
# Copyright (c) 2013: Clp.jl contributors
#
# Use of this source code is governed by an MIT-style license that can be found
# in the LICENSE.md file or at https://opensource.org/licenses/MIT.

function warmup()
    bridge = MOI.instantiate(Optimizer; with_bridge_type=Float64)
    no_bridge = MOI.instantiate(Optimizer)
    cache = MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}())
    MOI.optimize!(bridge, cache)
    MOI.optimize!(no_bridge, cache)
    return
end

function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    precompile(warmup, ())
    return
end

_precompile_()


text:
using HDF5
using Test

@testset "mmap" begin

# Create a new file
fn = tempname()
f = h5open(fn, "w")
@test isopen(f)

# Write  HDF5 file
hdf5_A = d_create(f,"A",datatype(Int64),dataspace(3,3));
A = rand(Int64,3,3)
hdf5_A[:,:] = A
flush(f)
close(f)
# Read HDF5 file & MMAP
f = h5open(fn,"r")
A_mmaped = readmmap(f["A"])
@test all(A .== A_mmaped)
# Check that it is read only
@test_throws ReadOnlyMemoryError A_mmaped[1,1] = 33
close(f)
# Now check if we can write
f = h5open(fn,"r+")
A_mmaped = readmmap(f["A"])
A_mmaped[1,1] = 33
close(f)

end


text:
using Distributions, SpecialFunctions
import Distributions: mean

""" Exact Bayes (Adams and MacKay, 2007) for the Gaussian task (section 2.3.1)
"""
struct GBayesFilter{D, A} <: BayesFilter
    posteriors::Array{D, 1}
    prior::D
    probabilities::Array{Float64, 1}
    sigma::Float64
    changeprobability::Float64
    cutoff::Float64
    approximation::A
    mu_n::Array{Float64, 1}
end
function GBayesFilter(;sigma = .1, mu_0 = 0., sigma_0 = 1.,
                    changeprobability = 0.01, cutoff = 1e-16)

    prior = GaussKnownVariance(mu_0, sigma_0, sigma)
    approximation = CutoffApprox(cutoff)
    GBayesFilter([deepcopy(prior)], prior, [1.], sigma, changeprobability,
                cutoff, approximation, [0.])
end
meanparameter(f::Union{GBayesFilter, GMPN, GSOR}) = sum(postmean.(f.posteriors) .* f.probabilities)
function update!(f::Union{GBayesFilter, GMPN, GSOR}, x)
    update_posterior!(f, x)
    approximate!(f.approximation, f)
    parameter_estimate = meanparameter(f)       # E[θ | x_{1:t}]
    update_prior!(f)
    f.mu_n[1] = copy(parameter_estimate)
end

mutable struct GaussKnownVariance
    mu_0::Float64
    sigma_0::Float64
    sigma::Float64
    ν::Float64
    χ::Float64
end
function GaussKnownVariance(mu_0, sigma_0, sigma)
    ν_0 = calcν(sigma, sigma_0)
    χ_0 = calcχ(mu_0, ν_0, sigma)

    GaussKnownVariance(mu_0, sigma_0, sigma, ν_0, χ_0)
end
dataprobability(d::GaussKnownVariance, x) = gauss_χ_ν(d, x)
function gauss_χ_ν(d::GaussKnownVariance, x)
    expterm = - x^2/(2. * d.sigma^2)
    expterm -= d.χ^2 / (2. * d.ν)
    expterm += (d.χ + x/d.sigma)^2 / (2. * (d.ν + 1.))
    if isnan(expterm)
        error("expterm is NaN")
    end
    sqrtterm = sqrt(d.ν/(d.ν + 1.)) * (1. / (sqrt(2π) * d.sigma) )
    dataprobability = exp(expterm) * sqrtterm
    if isinf(dataprobability)
         error("dataprobability is Inf")
    end
    dataprobability
end

function update!(d::Vector{<:GaussKnownVariance}, x)
    @inbounds for i in 1:length(d)
        d[i].ν += 1.
        d[i].χ += x/d[i].sigma
    end
end
vargauss(d::GaussKnownVariance) = 1/(1/(d.sigma_0^2) + d.ν/(d.sigma^2))
mean(d::GaussKnownVariance) = vargauss(d) * (d.mu_0/(d.sigma_0^2) + d.χ/(d.sigma^2))

function postmean(d::GaussKnownVariance)
    mu = d.sigma * d.χ / d.ν
end


text:
@testset "How to stop NLP" begin
###############################################################################
#
# The Stopping structure eases the implementation of algorithms and the
# stopping criterion.
#
# We illustrate here the basic features of NLPStopping, which is a
# specialized version of the Stopping to the case where:
# pb is an AbstractNLPModel
# state shares the structure of the NLPAtX.
#
# NLPModels is a package to handle non-linear (constrained) optimization.
# NLPStopping is following this approach.
#
###############################################################################
#using Test, NLPModels, Stopping

#We first create a toy problem
f(x) = sum(x.^2)
x0 = zeros(5)
nlp = ADNLPModel(f, x0)
#Warning: see https://github.com/JuliaSmoothOptimizers/NLPModels.jl/blob/master/src/autodiff_model.jl
#for the proper way of defining an ADNLPModel
meta = NLPModelMeta(5, x0=x0, lvar=zeros(5), uvar = Inf * ones(5))
nlp2 = ADNLPModel(meta, Counters(), ADNLPModels.ForwardDiffAD(5), f,  x->[])
nlp_at_x = NLPAtX(x0)
x1 = ones(5)

###############################################################################
#1) Initialize the NLPStopping.
#The specificity here is that the NLPStopping requires another mandatory input:
#optimality_check which is the function later used to compute the score.
#Recall that the score is then tested at 0, to declare optimality.
#Stopping provides a default KKT function.
stop_nlp = NLPStopping(nlp, nlp_at_x, optimality_check = (x,y) -> KKT(x,y))

#Another approach is to use the lazy way:
stop_nlp_lazy = NLPStopping(nlp2) #use nlp.meta.x0 as initial point
@test stop_nlp_lazy.current_state.x == nlp2.meta.x0

###############################################################################
#2) Fill in
#Before calling start! and stop! one should fill in current information in the
#State. -> the optimality_check then exploits this knowledge.
#As seen before, we by hand use update_and_start and update_and_stop.
#Another way is to call the fill_in! function:
fill_in!(stop_nlp, x1, matrix_info = false)
@test stop_nlp.current_state.x  == x1
@test stop_nlp.current_state.fx == 5.
@test stop_nlp.current_state.Hx == zeros(0,0)
#Note that since there are no constraints, c(x) and J(x) are not called:
@test stop_nlp.current_state.Jx == zeros(0,0)
@test stop_nlp.current_state.cx == zeros(0)
#Since there are no bounds on x, the Lagrange multiplier is not updated:
@test stop_nlp.current_state.mu == zeros(0)

#would give Hx if matrix_info = true
fill_in!(stop_nlp_lazy, x1)
@test stop_nlp_lazy.current_state.Hx != zeros(0,0)
#stop_nlp_lazy.pb has bounds, so mu is a vector of size x
@test size(x0) == size(stop_nlp_lazy.current_state.mu)

###############################################################################
#3) Evaluations
#Another particularity is that the NLPModels has a counter keeping track of
#the evaluations of each function.
#Similarly the NLPStopping has a dictionary keeping all the maximum number of
#evaluations:
@test typeof(stop_nlp.meta.max_cntrs) <: Dict
#For instance the limit in evaluations of objective and gradient:
@test stop_nlp.meta.max_cntrs[:neval_obj] == typemax(Int)
@test stop_nlp.meta.max_cntrs[:neval_grad] == typemax(Int)

#Limit can be set using init_max_counters function:
stop_nlp.meta.max_cntrs = init_max_counters(obj = 3, grad = 0, hess = 0)
@test stop_nlp.meta.max_cntrs[:neval_obj] == 3
@test stop_nlp.meta.max_cntrs[:neval_grad] == 0

OK = update_and_stop!(stop_nlp, evals = stop_nlp.pb.counters)
@test OK == true
@test stop_nlp.meta.resources == true
@test status(stop_nlp) == :EvaluationLimit

###############################################################################
#4) Unbounded problem
#An additional feature of the NLPStopping is to provide an _unbounded_problem_check
#whenever \|c(x)\| or -f(x) become too large.
stop_nlp.meta.unbounded_threshold = -6.0 #by default 1.0e50
stop!(stop_nlp)
@test stop_nlp.meta.unbounded_pb == true
@test stop_nlp.current_state.fx > stop_nlp.meta.unbounded_threshold
@test stop_nlp.meta.resources == true #still true as the state has not changed

###############################################################################
#An advanced feature is the possibility to send keywords to optimality_check:
optimality_fct_test = (x,y;a = 1.0) -> a

#In this case, the optimality_check function used to compute the score may
#depend on a parameter (algorithm-dependent for instance)
stop_nlp_2 = NLPStopping(nlp, nlp_at_x, optimality_check = optimality_fct_test)
fill_in!(stop_nlp_2, x0)

OK = stop!(stop_nlp_2, a = 0.0)
@test OK == true
@test stop_nlp_2.meta.optimal == true

#However, note that the same cannot be achieved with update_and_stop!:
reinit!(stop_nlp_2)
OK = update_and_stop!(stop_nlp_2, a = 0.0)
@test OK == false

end


text:
function get_ae(obj::SphericalHarmonicsProvider)
    return jcall(obj, "getAe", jdouble, ())
end

function get_max_degree(obj::SphericalHarmonicsProvider)
    return jcall(obj, "getMaxDegree", jint, ())
end

function get_max_order(obj::SphericalHarmonicsProvider)
    return jcall(obj, "getMaxOrder", jint, ())
end

function get_mu(obj::SphericalHarmonicsProvider)
    return jcall(obj, "getMu", jdouble, ())
end

function get_offset(obj::SphericalHarmonicsProvider, arg0::AbsoluteDate)
    return jcall(obj, "getOffset", jdouble, (AbsoluteDate,), arg0)
end

function get_reference_date(obj::SphericalHarmonicsProvider)
    return jcall(obj, "getReferenceDate", AbsoluteDate, ())
end

function get_tide_system(obj::TideSystemProvider)
    return jcall(obj, "getTideSystem", TideSystem, ())
end

function on_date(obj::NormalizedSphericalHarmonicsProvider, arg0::AbsoluteDate)
    return jcall(obj, "onDate", NormalizedSphericalHarmonicsProvider_NormalizedSphericalHarmonics, (AbsoluteDate,), arg0)
end


text:
"""
    SteepestDescent

Descent method which minimizes the objective function in the direction
of antigradient at each step.
"""
mutable struct SteepestDescent{T<:AbstractFloat,V<:AbstractVector{T}} <: OptBuffer
    x::V
    xpre::V
    g::V
    dir::V
    y::T
    ypre::T
    α::T
end

function SteepestDescent(x::AbstractVector, α::Real)
    F = float(eltype(x))
    return SteepestDescent(
        similar(x, F),
        similar(x, F),
        similar(x, F),
        similar(x, F),
        F(NaN),
        F(NaN),
        convert(F, α),
    )
end

SteepestDescent(x::AbstractVector) = SteepestDescent(x, 1)

"""
`optfn!` must be the 3-arg closure that computes fdf(x + α*d) and overwrites `M`'s gradient
"""
function init!(optfn!, ::SteepestDescent{T}, x0; kw...) where {T}
    optfn!(x0, zero(T), x0)
    return
end

@inline reset!(::SteepestDescent) = nothing

function reset!(M::SteepestDescent, x0, α = M.α)
    if length(M.x) != length(x0)
        foreach((M.x, M.xpre, M.g, M.dir)) do v
            resize!(v, length(x0))
        end
    end
    M.α = α
    return
end

function callfn!(fdf, M::SteepestDescent, x, α, d)
    __update_arg!(M, x, α, d)
    y, g = fdf(M.x, M.g)
    __update_grad!(M, g)
    M.y = y
    return y, g
end

@inline function __descent_dir!(M::SteepestDescent)
    M.dir .= .-M.g
    return M.dir
end

function step!(optfn!, M::SteepestDescent; constrain_step = infstep)
    M.x, M.xpre = M.xpre, M.x
    M.ypre = M.y
    d = __descent_dir!(M)
    xpre = M.xpre
    αmax = constrain_step(xpre, d)
    α = strong_backtracking!(optfn!, xpre, d, M.ypre, M.g, α=M.α, αmax=αmax, β=1e-4, σ=0.1)
    return α
end

@inline function __update_arg!(M::SteepestDescent, x, α, d)
    map!(M.x, d, x) do a, b
        muladd(α, a, b)
    end
    return
end

@inline function __update_arg!(M::SteepestDescent, x)
    if x !== M.x
        copy!(M.x, x)
    end
    return
end

@inline function __update_xpre!(M::SteepestDescent, x)
    if x !== M.xpre
        copy!(M.xpre, x)
    end
    return
end

@inline function __update_grad!(M::SteepestDescent, g)
    if M.g !== g
        copy!(M.g, g)
    end
    return
end


text:
#=
jim.jl
jiffy image display
2019-02-23 Jeff Fessler, University of Michigan
=#

export jim

using Plots: heatmap, ColorGradient
using LaTeXStrings
using MosaicViews: mosaicview
using FFTViews: FFTView
using Colors: HSV
using Test: @test

# global default key/values
jim_def = Dict([
	:aspect_ratio => :equal,
	:clim => nothing,
	:color => :grays,
	:line3plot => true, # lines around sub image for 3d mosaic?
	:line3type => (:yellow),
	:ncol => 0,
	:padval => nothing,
	:mosaic_npad => 1,
	:title => "",
	:xlabel => "",
	:ylabel => "",
	:fft0 => false,
	:yflip => nothing, # defer to minimum value of y - see default below
	:abswarn => true, # warn when taking abs of complex images?
	])

"""
`nothing_else(x, y)`
return `y` if `x` is nothing, else return `x`
"""
function nothing_else(x, y)
	return x == nothing ? y : x
end

"""
`jim(z, ...)`

jiffy image display of `x` using `heatmap`

in
- `z` image, can be 2D or higher, if higher then it uses `mosaicviews`

option
- `aspect_ratio` for heatmap; default `:equal`
- `clim` for heatmap; default `(minimum(z),maximum(z))`
- `color` colormap; default `:grays`
- `ncol` for mosaicview for 3D and higher arrays; default `0` does auto select
- `padval` padding value for mosaic view; default `minimum(z)`
- `line3plot` lines around sub image for 3d mosaic; default `false`
- `line3type` line type around sub image for 3d mosaic; default `(:yellow)`
- `mosaic_npad` # of pixel padding for mosaic view; default `1`
- `fft0` if true use FFTView to display; default `false`
- `title` for heatmap; default `""`
- `xlabel` for heatmap; default `""`
- `ylabel` for heatmap; default `""`
- `yflip` for heatmap; default `true` if `minimum(y) >= 0`
- `x` for x axis; default `1:size(z,1)`
- `y` for y axis; default `1:size(z,2)`
- `xtick` for heatmap; default `[minimum(x),maximum(x)]`
- `ytick` for heatmap; default `[minimum(y),maximum(y)]`

out
- returns plot handle

2019-02-23 Jeff Fessler, University of Michigan
"""
function jim(z::AbstractArray{<:Real} ;
		aspect_ratio = jim_def[:aspect_ratio],
		clim = nothing_else(jim_def[:clim], (minimum(z), maximum(z))),
		color = jim_def[:color],
		line3plot = jim_def[:line3plot],
		line3type = jim_def[:line3type],
		ncol::Int = jim_def[:ncol],
		padval = nothing_else(jim_def[:padval], minimum(z)),
		mosaic_npad::Int = jim_def[:mosaic_npad],
		title::String = jim_def[:title],
		xlabel::String = jim_def[:xlabel],
		ylabel::String = jim_def[:ylabel],
		fft0::Bool = jim_def[:fft0],
		x = fft0 ? Int.(-size(z,1)/2:size(z,1)/2-1) : (1:size(z,1)),
		y = fft0 ? Int.(-size(z,2)/2:size(z,2)/2-1) : (1:size(z,2)),
		xtick = (minimum(x) < 0 && maximum(x) > 0) ?
			[minimum(x),0,maximum(x)] : [minimum(x),maximum(x)],
		ytick = (minimum(y) < 0 && maximum(y) > 0) ?
			[minimum(y),0,maximum(y)] : [minimum(y),maximum(y)],
		yflip::Bool = nothing_else(jim_def[:yflip], minimum(y) >= 0),
		abswarn::Bool = jim_def[:abswarn], # ignored here
	)

	n1,n2,n3 = size(z,1), size(z,2), size(z,3)
	xy = (x, y)
	if ndims(z) > 2
		n1 += mosaic_npad
		n2 += mosaic_npad
		n3 = size(z,3)
		if ncol == 0
			ncol = Int(floor(sqrt(prod(size(z)[3:end]))))
		end
		z = mosaicview(z, padval, ncol=ncol, npad=mosaic_npad)
		xy = () # no x,y for mosaic
	elseif fft0
		z = FFTView(z)[x,y]
	end

	# attempt an HSV colormap for phase images
	if color === :hsv
		color = ColorGradient([HSV(h,1,1) for h=LinRange(0,350,351)])
	end

	heatmap(xy..., z', transpose=false,
		aspect_ratio=aspect_ratio,
		clim=clim,
		color=color,
		title=title,
		yflip=yflip,
		xlabel=xlabel,
		ylabel=ylabel,
		xtick=xtick,
		ytick=ytick)

	if n3 > 1 && line3plot # lines around each subimage
		m1 = (1+size(z,1)) / n1 # add one because of mosaicview non-edge
		m2 = (1+size(z,2)) / n2
		plot_box! = (ox, oy) -> plot!(
			ox .+ [0,1,1,0,0]*n1 .+ 0.0,
			oy .+ [0,0,1,1,0]*n2 .+ 0.0,
			line=jim_def[:line3type], label="")
		for ii=0:n3-1
			i1 = mod(ii, m1)
			i2 = floor(ii / m1)
			plot_box!(i1*n1, i2*n2)
		end
	end
	plot!()

end # jim

# handle case of complex image
function jim(z::AbstractArray{<:Number} ;
			abswarn::Bool = jim_def[:abswarn],
			kwargs...
		)

	if !(eltype(z) <: Real)
		abswarn && (@warn "magnitude at $(caller_name())")
		z = abs.(z)
	end
	jim(z ; kwargs...)
end

"""
`jim(z, title ; kwargs...)`
"""
function jim(z::AbstractArray{<:Number}, title::Union{String,LaTeXString}
		; kwargs...)
	return jim(z ; title=title, kwargs...)
end

"""
`jim(x, y, z ; kwargs...)`
"""
function jim(x, y, z ; kwargs...)
	return jim(z ; x=x, y=y, kwargs...)
end

"""
`jim(x, y, z, title ; kwargs...)`
"""
function jim(x, y, z, title::Union{String,LaTeXString} ; kwargs...)
	return jim(z ; x=x, y=y, title=title, kwargs...)
end

"""
`jim()`
return docstring if user calls `jim()` with no arguments
"""
function jim()
	@doc jim
end

#"""
#`jim(abswarn=false)`
#"""
#function jim( ; abswarn::Bool=jim_state_abswarn)
#	global jim_state_abswarn = abswarn
#end

"""
`jim(key, value)` set default value for one of the keys
"""
function jim(key::Symbol, value)
	global jim_def
	!haskey(jim_def, key) && throw(ArgumentError("no key $key"))
	jim_def[key] = value
end

"""
`jim(:test)`

`jim(:keys)` return default keys

`jim(:defs)` return `Dict` of default keys / vals

`jim(:key)` return `Dict[key]` if possible
"""
function jim(test::Symbol)
	global jim_def
	if test === :keys
		return keys(jim_def)
	end
	if test === :defs
		return jim_def
	end
	if haskey(jim_def, test)
		return jim_def[test]
	end
	test != :test && throw("symbol $test")
	jim()
	jim(:keys)
	jim(:clim)
	@test typeof(jim(:defs)) <: Dict

	jim(ones(4,3), title="test2")
	jim(ones(4,3,5), title="test3")
	jim(1:4, 5:9, zeros(4,5), title="test3")
	jim(zeros(4,5), x=1:4, y=5:9, title="test3")
	jim(zeros(4,6), fft0=true)
	jim(x=1:4, y=5:9, rand(4,5), title="test4")
	jim(rand(4,5), color=:hsv)
	jim(:abswarn, false)
	jim(complex(rand(4,3)))
	jim(complex(rand(4,3)), "complex")
	jim(:abswarn, true)
	true
end


text:
module Transformer

# 1. first try EncoderLayer on its own:
include("layers.jl")
d = 4 # model size
t = 5 # timesteps
layer = EncoderLayer(d, 2d)
layer = Flux.mapleaves(Flux.Tracker.data, layer) # disable Tracker
x = rand(d, t)
layer(x)

# 2. then try Encoder
include("posenc.jl")
v = 10 # vocab size
n = 3 # number of layers
encoder = Encoder(d, 2d, n, v)
encoder = Flux.mapleaves(Flux.Tracker.data, encoder)
x = rand(1:v, t)
encoder(x)

# 3. try with real data and Minibatch.jl
include("data.jl")
encoder = Encoder(d, 2d, n, length(vocab_en))
encoder = Flux.mapleaves(Flux.Tracker.data, encoder)
x = en[1]
encoder(x)

end # module


text:
using DiffEqPhysics

function generate_bodies_in_cell_nodes(n::Int, m::Real, v_dev::Real, L::Real)

    rng = MersenneTwister(n);
    velocities = v_dev * randn(rng, Float64, (3, n))
    bodies = MassBody[]

    count = 1
    dL = L / (ceil(n^(1 / 3)))
    for x = dL / 2:dL:L, y = dL / 2:dL:L, z = dL / 2:dL:L
        if count > n
            break
        end
        r = SVector(x, y, z)
        v = SVector{3}(velocities[:,count])
        body = MassBody(r, v, m)
        push!(bodies, body)
        count += 1
    end
    return bodies
end

const T = 120.0 # °K
const kb = 1.38e-23 # J/K
const ϵ = T * kb
const σ = 3.4e-10 # m
const ρ = 1374 # kg/m^3
const m = 39.95 * 1.6747 * 1e-27 # kg
const N = 216#floor(Int, ρ * L^3 / m)
const L = (m*N/ρ)^(1/3)#10.229σ
const R = 2.25σ
const v_dev = sqrt(kb * T / m)
const τ = 0.5e-15 # σ/v, fs
const t1 = 0.0
const t2 = 3000τ

const _L = L / σ
const _σ = 1.0
const _ϵ = 1.0
const _m = 1.0
const _v = v_dev / sqrt(ϵ / m)
const _R = R / σ
const _τ = τ * sqrt(ϵ / m) / σ
const _t1 = t1 * sqrt(ϵ / m) / σ
const _t2 = t2 * sqrt(ϵ / m) / σ
bodies = generate_bodies_in_cell_nodes(N, _m, _v, _L)
parameters = LennardJonesParameters(_ϵ, _σ, _R)
lj_system = PotentialNBodySystem(bodies, Dict(:lennard_jones => parameters));
simulation = NBodySimulation(lj_system, (_t1, _t2), CubicPeriodicBoundaryConditions(_L), _ϵ/T);
#result = run_simulation(simulation, Tsit5())
result = @time run_simulation(simulation, VelocityVerlet(), dt=_τ)

#(rs, grf) = rdf(result)
#(ts, dr2) = msd(result)

#using Plots
#plot(rs, grf, xlim=[0, 0.4999_L], label=["Radial distribution function"],ylabel="g(r)", xlabel="r")
#plot(rs/_σ, grf, xlim=[0, 0.4999_L/_σ], label=["Radial distribution function"],ylabel="g(r)", xlabel="r/sigma")

time_now = Dates.format(now(), "yyyy_mm_dd_HH_MM_SS")
Nactual = length(bodies)
timesteps = round(length(result.solution.t))

@time save_to_pdb(result, "D:/liquid argon simulation $Nactual molecules and $timesteps steps $time_now.pdb" )


text:
function addSavingRAM_!(p::Program_,platform::String)

    #List of nonempty arrays
    l = []

    if length(p.agent.declaredSymbols["Local"])>0
        push!(l,:(Core.Array(localV)[1:N,:]))
    else
        push!(l,:(Core.Array{Float64,2}(undef,0,2)))
    end

    if length(p.agent.declaredSymbols["Identity"])>0
        push!(l,:(Core.Array(identityV)[1:N,:]))
    else
        push!(l,:(Core.Array{Int,2}(undef,0,2)))
    end

    if length(p.agent.declaredSymbols["Global"])>0
        push!(l,:(Core.Array(globalV)))
    else
        push!(l,:(Core.Array{Float64,1}()))
    end

    if length(p.agent.declaredSymbols["GlobalArray"]) > 0
        list = string("[copy(",p.agent.declaredSymbols["GlobalArray"][1],"")
        list = string(list,")]")
        push!(l,Meta.parse(string("Core.Array(",list,")")))
    else
        push!(l,:(Core.Array{Core.Array{Float64},1}()))
    end

    if length(p.agent.declaredSymbols["Medium"]) > 0
        if p.agent.dims == 1
            push!(l,:(Core.Array(mediumV)[2:end-1,:]))
        elseif p.agent.dims == 2
            push!(l,:(Core.Array(mediumV)[2:end-1,2:end-1,:]))
        elseif p.agent.dims == 3
            push!(l,:(Core.Array(mediumV)[2:end-1,2:end-1,2:end-1,:]))
        end
    else
        push!(l,:(Core.Array{Float64,1}()))
    end

    push!(p.declareVar.args,:(commRAM_ = CommunityInTime()))

    push!(p.execInit.args,
        :(begin
            ob = Community($(p.agent.dims),t,N,com.mediumN,com.simulationBox,com.radiusInteraction,com.declaredSymbols_,$(l...))
            push!(commRAM_,ob)
        end)
    )
    push!(p.execInloop.args,
        :(begin
            if t >= tSave
                tSave += dtSave
                ob = Community($(p.agent.dims),t+dt,N,com.mediumN,com.simulationBox,com.radiusInteraction,com.declaredSymbols_,$(l...))
                push!(commRAM_,ob)
            end
        end)
    )

    # push!(p.execAfter.args,
    # :(begin
    #     ob = Community(t+dt,N,com.declaredSymbols_,$(l...))
    #     push!(commRAM_,ob)
    #     end)
    # )

    push!(p.returning.args,:(commRAM_))

    return
end


text:
# Julia wrapper for header: /usr/include/libavutil/audio_fifo.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

export
    av_audio_fifo_free,
    av_audio_fifo_alloc,
    av_audio_fifo_realloc,
    av_audio_fifo_write,
    av_audio_fifo_peek,
    av_audio_fifo_peek_at,
    av_audio_fifo_read,
    av_audio_fifo_drain,
    av_audio_fifo_reset,
    av_audio_fifo_size,
    av_audio_fifo_space

function av_audio_fifo_free(af)
    ccall((:av_audio_fifo_free, libavutil), Cvoid, (Ptr{AVAudioFifo},), af)
end

function av_audio_fifo_alloc(sample_fmt::AVSampleFormat, channels::Integer, nb_samples::Integer)
    ccall((:av_audio_fifo_alloc, libavutil), Ptr{AVAudioFifo}, (AVSampleFormat, Cint, Cint), sample_fmt, channels, nb_samples)
end

function av_audio_fifo_realloc(af, nb_samples::Integer)
    ccall((:av_audio_fifo_realloc, libavutil), Cint, (Ptr{AVAudioFifo}, Cint), af, nb_samples)
end

function av_audio_fifo_write(af, data, nb_samples::Integer)
    ccall((:av_audio_fifo_write, libavutil), Cint, (Ptr{AVAudioFifo}, Ptr{Ptr{Cvoid}}, Cint), af, data, nb_samples)
end

function av_audio_fifo_peek(af, data, nb_samples::Integer)
    ccall((:av_audio_fifo_peek, libavutil), Cint, (Ptr{AVAudioFifo}, Ptr{Ptr{Cvoid}}, Cint), af, data, nb_samples)
end

function av_audio_fifo_peek_at(af, data, nb_samples::Integer, offset::Integer)
    ccall((:av_audio_fifo_peek_at, libavutil), Cint, (Ptr{AVAudioFifo}, Ptr{Ptr{Cvoid}}, Cint, Cint), af, data, nb_samples, offset)
end

function av_audio_fifo_read(af, data, nb_samples::Integer)
    ccall((:av_audio_fifo_read, libavutil), Cint, (Ptr{AVAudioFifo}, Ptr{Ptr{Cvoid}}, Cint), af, data, nb_samples)
end

function av_audio_fifo_drain(af, nb_samples::Integer)
    ccall((:av_audio_fifo_drain, libavutil), Cint, (Ptr{AVAudioFifo}, Cint), af, nb_samples)
end

function av_audio_fifo_reset(af)
    ccall((:av_audio_fifo_reset, libavutil), Cvoid, (Ptr{AVAudioFifo},), af)
end

function av_audio_fifo_size(af)
    ccall((:av_audio_fifo_size, libavutil), Cint, (Ptr{AVAudioFifo},), af)
end

function av_audio_fifo_space(af)
    ccall((:av_audio_fifo_space, libavutil), Cint, (Ptr{AVAudioFifo},), af)
end


text:
sig = [2 3; 3 2; 1 1]
S = random_S0Graph(sig)

w = random_S0_density(S.sig)
w /= tr(w)

#@show S
#display(w)

λ, x1, Z = dsw_schur2(S)
problem = maximize(trace_logm(x1, w), [ λ <= 1 ])
@time solve!(problem, () -> SCS.Optimizer(verbose=0, eps=solver_eps))
h1=problem.optval
x1=Hermitian(evaluate(x1))

λ, x2, Z = dsw_schur2(complement(S))
problem = maximize(trace_logm(x2, w), [ λ <= 1 ])
@time solve!(problem, () -> SCS.Optimizer(verbose=0, eps=solver_eps))
h2=problem.optval
x2=Hermitian(evaluate(x2))

@test x1*x2 ≈ x2*x1  atol=tol
@test x1*x2 ≈ Ψ(S, w)  atol=tol
# follows from the above
@test h1+h2 ≈ tr(w*log(Ψ(S, w)))  atol=tol

max_stars_repo_name: dstahlke/NoncommutativeGraphs

max_issues_repo_name: dstahlke/NoncommutativeGraphs

max_forks_repo_name: dstahlke/NoncommutativeGraphs


text:
# DX = [tx,ty]
function odoAdd(X::Array{Float64,1}, DX::Array{Float64,1})
    len = length(X)
    A = eye(3)
    #@show A[2,3] = X[2]
    A[1:len,3] = X

    B = eye(3)
    B[1:len,3] = DX
    #B[2,3] = DX[2]
    retval = zeros(len,1)
    AB = A*B
    retval[1:len,1] = AB[1:len,3]
    #retval[2,1] = AB[2,3]
    return retval
end

function evalPotential(odom::Odo, Xi::Array{Graphs.ExVertex,1}, Xid::Int64)
    rz,cz = size(odom.Zij)
    # implicit equation portion -- bi-directional pairwise function
    if Xid == Xi[1].index #odom.
        #Z = (odom.Zij\eye(rz)) # this will be used for group operations
        Z = - odom.Zij
        Xval = getVal(Xi[2])
    elseif Xid == Xi[2].index
        Z = odom.Zij
        Xval = getVal(Xi[1])
    else
        error("Bad evalPairwise Odo")
    end

    r,c = size(Xval)
    RES = zeros(r,c*cz)
    # increases the number of particles based on the number of modes in the measurement Z
    for i in 1:(c*cz)
        ent = diag(odom.Cov).*randn(size(vec(Z[:,floor(Int,i/(c+1)+1)])))
        RES[:,i] = odoAdd(Xval[:,i], ent+Z[:,floor(Int,i/(c+1)+1)])
    end
    return RES
end

function evalPotential(odom::OdoMM, Xi::Array{Graphs.ExVertex,1}, Xid::Int64)
    rz,cz = size(odom.Zij)
    # implicit equation portion -- bi-directional pairwise function
    if Xid == Xi[1].index #odom.
        #Z = (odom.Zij\eye(rz)) # this will be used for group operations
        Z = -odom.Zij

        XvalMM = Array{Array{Float64,2},1}(2)
        for i in 2:3
            XvalMM[i-1] = getVal(Xi[i])
        end

        len1 = size(XvalMM[1],2)
        len2 = size(XvalMM[2],2)
        if len1 != len2
            error("odoMM -- must be the same size $(len1), $(len2)")
        #    tmp = kde!(XvalMM[2][:,:], "lcv")
        #    XvalMM[2], = sample(p,len1)
        #elseif len2 > len1
        #    tmp = kde!(XvalMM[1][:,:], "lcv")
        #    XvalMM[1], = sample(p,len2)
        end

        len = len1
        s = rand(len).>0.5 # assume fixed weight for now
        p = (1:len)[s]
        np = (1:len)[!s]
        Xval = zeros(size(XvalMM[1]))
        Xval[:,p] = XvalMM[1][:,p]
        Xval[:,np] = XvalMM[2][:,np]
    elseif Xid == Xi[2].index || Xid == Xi[3].index
        Z = odom.Zij
        Xval = getVal(Xi[1])
    else
        error("Bad evalPairwise OdoMM")
    end

    r,c = size(Xval)
    RES = zeros(r,c*cz)
    # increases the number of particles based on the number of modes in the measurement Z
    for i in 1:(c*cz)
        ent = diag(odom.Cov).*randn(size(vec(Z[:,floor(Int,i/(c+1)+1)])))
        RES[:,i] = odoAdd(Xval[:,i], ent+Z[:,floor(Int,i/(c+1)+1)])
    end
    return RES
end

function rangeAdd(X::Array{Float64,1}, DX::Array{Float64,1})
    A = eye(2)
    A[1,2] = X[1]
    B = eye(2)
    B[1,2] = DX[1]
    return [(A*B)[1,2]]
end

function evalPotential(rang::Ranged, Xi::Array{Graphs.ExVertex,1}, Xid::Int64)
    if Xid == Xi[1].index #rang.
        Z = -rang.Zij
        Xval = getVal(Xi[2])
    elseif Xid == Xi[2].index
        Z = rang.Zij
        Xval = getVal(Xi[1])
    else
        error("Bad evalPairwise Ranged")
    end

    r,c = size(Xval)
    cz = size(rang.Zij,1)
    RES = zeros(r,c*cz)

    for i in 1:(c*cz) # for each mode in the measurement
        ent = rang.Cov[1]*randn(size(vec(Z[:,floor(Int,i/(c+1)+1)])))
        RES[:,i] = rangeAdd(Xval[:,i], ent+Z[floor(Int,i/(c+1)+1)])
    end
    return RES
end

function evalPotential(obs::Obsv2, Xi::Array{Graphs.ExVertex,1}; N::Int64=300)#, from::Int64)
    # @show obs.bws, typeof(obs.bws)
    pd = kde!(obs.pts, obs.bws[:,1])
    return KernelDensityEstimate.sample(pd,N)[1]
    # return obs.Cov[1]*randn()+obs.Zi
end

function evalPotentialSpecific(Xi::Array{Graphs.ExVertex,1}, typ::Singleton, solvefor::Int64; N::Int64=100)
  outpts = evalPotential(typ, Xi, N=N) # , solvefor
  return outpts
end

function evalPotentialSpecific(Xi::Array{Graphs.ExVertex,1}, typ::Pairwise, solvefor::Int64; N::Int64=100)
  return evalPotential(typ, Xi, solvefor)
end

# Multiple dispatch occurs internally, resulting in factor graph potential evaluations
function evalFactor2(fgl::FactorGraph, fct::Graphs.ExVertex, solvefor::Int64; N::Int64=100)
  # return evalPotential(fct.attributes["data"].fnc, solvefor) #evalPotential(fct.attributes["fnc"], solvefor)
  Xi = Graphs.ExVertex[]
  for id in fct.attributes["data"].fncargvID
    push!(Xi, dlapi.getvertex(fgl,id)) # TODO -- should use local mem only for this part, update after ## fgl.v[id]
  end
  return evalPotentialSpecific(Xi, fct.attributes["data"].fnc, solvefor, N=N) #evalPotential(fct.attributes["fnc"], solvefor)
end

function findRelatedFromPotential(fg::FactorGraph, idfct::Graphs.ExVertex, vertid::Int64, N::Int64) # vert
    # if vert.index == vertid
        ptsbw = evalFactor2(fg, idfct, vertid, N=N); # idfct[2] # assuming it is properly initialized TODO
        sum(abs(ptsbw)) < 1e-14 ? error("findRelatedFromPotential -- an input is zero") : nothing

        Ndim = size(ptsbw,1)
        Npoints = size(ptsbw,2)
        # Assume we only have large particle population sizes, thanks to addNode!
        p = kde!(ptsbw, "lcv")
        if Npoints != N # this is where we control the overall particle set size
            p = resample(p,N)
        end
        return p
    # end
    # return Union{}
end