text:
###############################################################################
#
# Leaf optical parameters
#
###############################################################################
"""
    mutable struct LeafOpticals{FT}

Struct for leaf optical properties

# Fields
$(TYPEDFIELDS)
"""
Base.@kwdef mutable struct LeafOpticals{FT}
    # TODO Add explanations to each field
    nr    ::Array{FT,1} = zeros(FT, 2)
    Km    ::Array{FT,1} = zeros(FT, 2)
    Kab   ::Array{FT,1} = zeros(FT, 2)
    Kant  ::Array{FT,1} = zeros(FT, 2)
    Kcar  ::Array{FT,1} = zeros(FT, 2)
    Kw    ::Array{FT,1} = zeros(FT, 2)
    KBrown::Array{FT,1} = zeros(FT, 2)
    phi   ::Array{FT,1} = zeros(FT, 2)
    KcaV  ::Array{FT,1} = zeros(FT, 2)
    KcaZ  ::Array{FT,1} = zeros(FT, 2)
    "Wave length `[nm]`, same as `WL` in [`WaveLengths`](@ref)`"
    lambda::Array{FT,1} = zeros(FT, 2)
end

max_stars_repo_name: climate-machine/Land


text:
function Base.show(io::IO, ::MIME"text/plain", m::Metabolite)
    for fname in fieldnames(Metabolite)
        if fname == :charge
            c = isnothing(getfield(m, fname)) ? nothing : string(getfield(m, fname))
            _pretty_print_keyvals(io, "Metabolite.$(string(fname)): ", c)
        else
            _pretty_print_keyvals(io, "Metabolite.$(string(fname)): ", getfield(m, fname))
        end
    end
end

ext: jl


text:
using Test
using DelimitedFiles, LinearAlgebra
import MLJBase
const Mlj = MLJBase
using StableRNGs
rng = StableRNG(123)
using BetaML.Clustering

println("*** Additional testing for the Clustering algorithms...")

#println("Testing MLJ interface for Clustering models....")
# evaluate seem not supported for unsupervised models


text:
@testset "func: vec()" begin
    A = TensorMap(randn, ComplexF64, ℂ^2 * ℂ^2 * ℂ^2, ℂ^2 * ℂ^2 * ℂ^2)
    Aarr = convert(Array, A)
    vec_Aarr = permutedims(
        reshape(permutedims(Aarr, (1, 2, 4, 5, 3, 6)), (2, 2, 2, 2, 4)), (1, 2, 5, 3, 4)
    )
    @test convert(Array, PEPOKit._vec(A)[1]) ≈ vec_Aarr
end

max_stars_repo_name: jack-dunham/PEPOKit

max_issues_repo_name: jack-dunham/PEPOKit

max_forks_repo_name: jack-dunham/PEPOKit


text:
#  ======================================================
#  Constructors
#  ======================================================

mutable struct SimpleValueGraph{V <: Integer,
                                E_VAL,
                                E_VAL_C <: Union{Nothing, EdgeValContainer{E_VAL}}
                               } <: AbstractSimpleValueGraph{V, E_VAL}
    ne::Int
    fadjlist::Adjlist{V}
    edgevals::E_VAL_C
end

const ConstValueGraph{V, E_VAL} = SimpleValueGraph{V, E_VAL, Nothing}
const InOutValueGraph{V, E_VAL} = SimpleValueGraph{V, E_VAL, <: EdgeValContainer{E_VAL}}

function SimpleValueGraph(nv::V,
                          E_VAL::Type=default_edgeval_type;
                          const_edgevals::Bool = false
                         ) where {V<:Integer}
    fadjlist = Adjlist{V}(nv)
    if const_edgevals
        return ConstValueGraph{V, E_VAL, Nothing}(nv, fadjlist, nothing)
    end
    edgevals = create_edgeval_list(nv, E_VAL)
    return InOutValueGraph{V, E_VAL, typeof(edgevals)}(V(0), fadjlist, edgevals)
end

SimpleValueGraph{V}(nv::Integer, E_VAL::Type=default_edgeval_type; args...) where {V} = SimpleValueGraph(V(nv), E_VAL, args...)
SimpleValueGraph{V, E_VAL}(nv::Integer, args...) where {V, E_VAL} = SimpleValueGraph(V(nv), E_VAL, args...)

# TODO weights are not symmetric
function SimpleValueGraph(g::SimpleGraph{V},
                          E_VAL::Type=default_edgeval_type;
                          const_edgevals::Bool=false,
                          edgeval_initializer = (s, d) -> default_edgeval(E_VAL)) where {V}
    n = nv(g)
    fadjlist = deepcopy_adjlist(g.fadjlist)
    if const_edgevals
        return ConstValueGraph{V, E_VAL}(nv, fadjlist, nothing)
    end
    edgevals = Vector{Vector{E_VAL}}(undef, n)
    for s in Base.OneTo(n)
        edgevals[s] = [edgeval_initializer(V(s), d) for d in fadjlist[s]]
    end
    InOutValueGraph{V, E_VAL, typeof(edgevals)}(ne(g), fadjlist, edgevals)
end

# TODO this function has some issues with typesafety
# TODO weights are not symmetric
# TODO const edgevals
function SimpleValueGraph(g::SimpleGraph{V},
                          E_VAL::Type{<:TupleOrNamedTuple};
                          const_edgevals::Bool=false,
                          edgeval_initializer = (s, d) -> default_edgeval(E_VAL)) where {V}
    n = nv(g)
    fadjlist = deepcopy_adjlist(g.fadjlist)
    if const_edgevals
        return ConstValueGraph{V, E_VAL}(nv, fadjlist, nothing)
    end

    E_VAL_C = edgevals_container_type(Val(E_VAL))
    edgevals = E_VAL_C( T(undef, n) for T in E_VAL_C.types )
    for s in Base.OneTo(n)
        for (i, T) in enumerate(E_VAL.types)
            edgevals[i][s] = Vector{T}(undef, length(fadjlist[s]))
        end
    end
    for s in Base.OneTo(n)
        for (j, d) in enumerate(fadjlist[s])
            w = edgeval_initializer(V(s), d)
            for (i, T) in enumerate(E_VAL.types)
                edgevals[i][s][j] = w[i]
            end
        end
    end
    InOutValueGraph{V, E_VAL, E_VAL_C}(ne(g), fadjlist, edgevals)
end

# =========================================================
# Interface
# =========================================================

function add_edge!(g::InOutValueGraph{V, E_VAL},
                   s::Integer,
                   d::Integer,
                   value=default_edgeval(E_VAL)) where {V, E_VAL}
    verts = vertices(g)
    (s in verts && d in verts) || return false # edge out of bounds
    edgevals = g.edgevals
    @inbounds list = g.fadjlist[s]
    index = searchsortedfirst(list, d)
    @inbounds if index <= length(list) && list[index] == d
        # edge already there, replace value, but return false
        set_value_for_index!(edgevals, s, index, value)
        # TODO maybe shortcircuit if it is a self-loop
        index = searchsortedfirst(g.fadjlist[d], s)
        set_value_for_index!(edgevals, d, index, value)
        return false
    end

    insert!(list, index, d)
    insert_value_for_index!(edgevals, s, index, value)
    g.ne += 1

    s == d && return true # selfloop

    @inbounds list = g.fadjlist[d]
    index = searchsortedfirst(list, s)
    insert!(list, index, s)
    insert_value_for_index!(edgevals, d, index, value)
    return true # edge successfully added
end

function add_edge!(g::ConstValueGraph,
                   s::Integer,
                   d::Integer)
    verts = vertices(g)
    (s in verts && d in verts) || return false # edge out of bounds
    @inbounds list = g.fadjlist[s]
    index = searchsortedfirst(list, d)
    @inbounds if index <= length(list) && list[index] == d
        # edge already there, return false
        return false
    end

    insert!(list, index, d)
    g.ne += 1

    s == d && return true # selfloop

    @inbounds list = g.fadjlist[d]
    index = searchsortedfirst(list, s)
    insert!(list, index, s)
    return true # edge successfully added
end

add_edge!(g::SimpleValueGraph, e::SimpleEdge)      = add_edge!(g, src(e), dst(e))
add_edge!(g::InOutValueGraph, e::SimpleEdge, u)   = add_edge!(g, src(e), dst(e), u)
add_edge!(g::InOutValueGraph, e::SimpleValueEdge) = add_edge!(g, src(e), dst(e), val(e))

function rem_edge!(g::InOutValueGraph, s::Integer, d::Integer)
    verts = vertices(g)
    (s in verts && d in verts) || return false # edge out of bounds
    edgevals = g.edgevals
    @inbounds list = g.fadjlist[s]
    index = searchsortedfirst(list, d)
    @inbounds (index <= length(list) && list[index] == d) || return false
    deleteat!(list, index)
    delete_value_for_index!(edgevals, s, index)

    g.ne -= 1
    s == d && return true # self-loop

    @inbounds list = g.fadjlist[d]
    index = searchsortedfirst(list, s)
    deleteat!(list, index)
    delete_value_for_index!(edgevals, d, index)

    return true
end

function rem_edge!(g::ConstValueGraph, s::Integer, d::Integer)
    verts = vertices(g)
    (s in verts && d in verts) || return false # edge out of bounds
    edgevals = g.edgevals
    @inbounds list = g.fadjlist[s]
    index = searchsortedfirst(list, d)
    @inbounds (index <= length(list) && list[index] == d) || return false
    deleteat!(list, index)

    g.ne -= 1
    s == d && return true # self-loop

    @inbounds list = g.fadjlist[d]
    index = searchsortedfirst(list, s)
    deleteat!(list, index)

    return true
end

rem_edge!(g::SimpleValueGraph, e::SimpleEdge) = rem_edge!(g, src(e), dst(e))
rem_edge!(g::InOutValueGraph, e::SimpleValueEdge) = rem_edge!(g, src(e), dst(e))

# TODO rem_vertex!, rem_vertices!

function has_edge(g::SimpleValueGraph, s::Integer, d::Integer)
    verts = vertices(g)
    (s in verts && d in verts) || return false # edge out of bounds
    @inbounds list_s = g.fadjlist[s]
    @inbounds list_d = g.fadjlist[d]
    if length(list_s) > length(list_d)
        d = s
        list_s = list_d
    end
    return LightGraphs.insorted(d, list_s)
end

function has_edge(g::SimpleValueGraph{V, E_VAL}, s::Integer, d::Integer, value) where {V, E_VAL}
    verts = vertices(g)
    (s in verts && d in verts) || return false # edge out of bounds
    @inbounds list_s = g.fadjlist[s]
    @inbounds list_d = g.fadjlist[d]
    if length(list_s) > length(list_d)
        s, d = d, s
        list_s = list_d
    end
    index = searchsortedfirst(list_s, d)
    @inbounds return (index <= length(list_s) && list_s[index] == d && value_for_index(g.edgevals, E_VAL, s, index) == value)
end

has_edge(g::SimpleValueGraph, e::SimpleEdge)      = has_edge(g, src(e), dst(e))
has_edge(g::SimpleValueGraph, e::SimpleEdge, u)   = has_edge(g, src(e), dst(e), u)
has_edge(g::SimpleValueGraph, e::SimpleValueEdge) = has_edge(g, src(e), dst(e), val(e))

# TODO rest methods for get_value
# TODO lots of duplicated code
function get_edgeval(g::InOutValueGraph{V, E_VAL}, s::Integer, d::Integer) where {V, E_VAL}
    verts = vertices(g)
    (s in verts && d in verts) || return nothing # TODO may raise bounds error?
    @inbounds list_s = g.fadjlist[s]
    @inbounds list_d = g.fadjlist[d]
    if length(list_s) > length(list_d)
        s, d = d, s
        list_s = list_d
    end
    index = searchsortedfirst(list_s, d)
    @inbounds if index <= length(list_s) && list_s[index] == d
        return value_for_index(g.edgevals, E_VAL, s, index)
    end
    return nothing
end

function get_edgeval(g::InOutValueGraph{V, E_VAL}, s::Integer, d::Integer, key) where {V <: Integer, E_VAL <: TupleOrNamedTuple}
    verts = vertices(g)
    (s in verts && d in verts) || return nothing # TODO may raise bounds error?
    @inbounds list_s = g.fadjlist[s]
    @inbounds list_d = g.fadjlist[d]
    if length(list_s) > length(list_d)
        s, d = d, s
        list_s = list_d
    end
    index = searchsortedfirst(list_s, d)
    @inbounds if index <= length(list_s) && list_s[index] == d
        return value_for_index(g.edgevals, E_VAL, s, index, key)
    end
    return nothing
end

function get_edgeval(g::ConstValueGraph{V, E_VAL}, s::Integer, d::Integer) where {V, E_VAL}
    if has_edge(g, s, d)
        return default_edgeval(E_VAL)
    end
    return nothing
end

function get_edgeval(g::ConstValueGraph{V, E_VAL}, s::Integer, d::Integer, key) where {V, E_VAL}
    if has_edge(g, s, d)
        return default_edgeval(E_VAL)[key]
    end
    return nothing
end

get_edgeval!(g::SimpleValueGraph, e::SimpleEdge)   = get_edgeval(g, src(e), dst(e))
get_edgeval!(g::SimpleValueGraph, e::SimpleEdge, key)   = get_edgeval(g, src(e), dst(e), key)

function set_edgeval!(g::InOutValueGraph, s::Integer, d::Integer, value)
    verts = vertices(g)
    (s in verts && d in verts) || return false
    @inbounds list = g.fadjlist[s]
    index = searchsortedfirst(list, d)
    @inbounds index <= length(list) && list[index] == d || return false
    edgevals = g.edgevals
    set_value_for_index!(edgevals, s, index, value)

    @inbounds list = g.fadjlist[d]
    index = searchsortedfirst(list, s)
    set_value_for_index!(edgevals, d, index, value)

    return true
end

# TODO create generated methods for NamedTuple (and also Tuple?)
function set_edgeval!(g::InOutValueGraph{V, <: TupleOrNamedTuple}, s::Integer, d::Integer, key, value) where {V}
    verts = vertices(g)
    (s in verts && d in verts) || return false
    @inbounds list = g.fadjlist[s]
    index = searchsortedfirst(list, d)
    @inbounds index <= length(list) && list[index] == d || return false
    edgevals = g.edgevals
    set_value_for_index!(edgevals, s, index, key, value)

    @inbounds list = g.fadjlist[d]
    index = searchsortedfirst(list, s)
    set_value_for_index!(edgevals, d, index, key, value)

    return true
end

set_edgeval!(g::InOutValueGraph, e::SimpleEdge, u)   = set_value!(g, src(e), dst(e), u)
set_edgeval!(g::InOutValueGraph, e::SimpleValueEdge) = set_value!(g, src(e), dst(e), edgeval(e))

is_directed(::Type{<:SimpleValueGraph}) = false
is_directed(g::SimpleValueGraph)        = false

outneighbors(g::SimpleValueGraph, v::Integer) = g.fadjlist[v]
inneighbors(g::SimpleValueGraph,  v::Integer) = outneighbors(g, v)

outedgevals(g::InOutValueGraph{V, E_VAL, <: Adjlist}, v::Integer) where {V, E_VAL} =
    g.edgevals[v]
outedegevals(g::ConstValueGraph{V, E_VAL}, v::Integer) where {V, E_VAL} = Base.Iterators.repeated(default_edgeval(E_VAL), length(g.fadjlist[v]))

inedgevals(g::SimpleValueGraph, v::Integer) = outedgevals(g, v)
all_edgevals(g::SimpleValueGraph, v::Integer) = outedgevals(g, v)

# TODO implement these with iterators instead of generators
outedgevals(g::InOutValueGraph{V, E_VAL, <: Tuple}, v::Integer) where {V, E_VAL} =
    ( Tuple( adjlist[v][i] for adjlist in g.edgevals ) for i in OneTo(length(g.fadjlist[v])) )
outedgevals(g::InOutValueGraph{V, E_VAL, T}, v::Integer) where {V, E_VAL, T <: NamedTuple} =
    ( NamedTuple{Tuple(T.names)}( adjlist[v][i] for adjlist in g.edgevals ) for i in OneTo(length(g.fadjlist[v])) )

outedgevals(g::InOutValueGraph{V, E_VAL, <: TupleOrNamedTuple}, v::Integer, key) where {V, E_VAL} =
    g.edgevals[key][v]
outedegevals(g::ConstValueGraph{V, E_VAL}, v::Integer, key) where {V, E_VAL} = Base.Iterators.repeated(default_edgeval(E_VAL)[key], length(g.fadjlist[v]))

inedgevals(g::SimpleValueGraph{V, E_VAL}, v::Integer, key) where {V, E_VAL} =
    outedgevals(g, v, key)
all_edgevals(g::SimpleValueGraph{V, E_VAL}, v::Integer, key) where {V, E_VAL} =
    outedgevals(g, v, key)

# TODO maybe add a sizehint kwarg
function add_vertex!(g::InOutValueGraph{V, E_VAL, E_VAL_C}) where {V, E_VAL, E_VAL_C}
    # TODO There are overflow checks in Julia Base, use these
    (nv(g) + one(V) <= nv(g)) && return false # overflow
    push!(g.fadjlist, V[])
    if E_VAL_C <: TupleOrNamedTuple
        for (i, T) in enumerate(E_VAL.types)
            push!(g.edgevals[i], T[])
        end
    else
        push!(g.edgevals, E_VAL[])
    end
    return true
end

function add_vertex!(g::ConstValueGraph{V}) where {V}
    # TODO There are overflow checks in Julia Base, use these
    (nv(g) + one(V) <= nv(g)) && return false # overflow
    push!(g.fadjlist, V[])
    return true
end

# ====================================================================
# Iterators
# ====================================================================

function iterate(iter::SimpleValueEdgeIter{<:InOutValueGraph}, state=(one(eltype(iter.g)), 1) )
    g = iter.g
    fadjlist = g.fadjlist
    V = eltype(g)
    E_VAL = edgeval_type(g)
    n::V = nv(g)
    u::V, i = state

    @inbounds while u < n
        if i > length(fadjlist[u])
            u += V(1)
            i = searchsortedfirst(fadjlist[u], u)
            continue
        end
        e = SimpleValueEdge(u, fadjlist[u][i], value_for_index(g.edgevals, E_VAL, u, i))
        return e, (u, i + 1)
    end

    (n == 0 || i > length(fadjlist[n])) && return nothing

    e = SimpleValueEdge(n, n, value_for_index(g.edgevals, E_VAL, n, 1))
    return e, (u, i + 1)
end

function iterate(iter::SimpleValueEdgeIter{<:ConstValueGraph}, state=(one(eltype(iter.g)), 1) )
    g = iter.g
    fadjlist = g.fadjlist
    V = eltype(g)
    n::V = nv(g)
    u::V, i = state

    @inbounds while u < n
        if i > length(fadjlist[u])
            u += V(1)
            i = searchsortedfirst(fadjlist[u], u)
            continue
        end
        e = SimpleEdge(u, fadjlist[u][i])
        return e, (u, i + 1)
    end

    (n == 0 || i > length(fadjlist[n])) && return nothing

    e = SimpleEdge(n, n, default_edgeval(E_VAL))
    return e, (u, i + 1)
end


text:
export AbstractSolution

"""
    AbstractSolution

Abstract type representing the solution of a rechability problem.
"""
abstract type AbstractSolution end


text:
using Plots, SNN

Ne = 800;      Ni = 200
E = SNN.IZ(;N = Ne, param = SNN.IZParameter(;a = 0.02, b = 0.2, c = -65, d = 8))
I = SNN.IZ(;N = Ni, param = SNN.IZParameter(;a = 0.1, b = 0.2, c = -65, d = 2))

EE = SNN.SpikingSynapse(E, E, :v; σ = 0.5,  p = 0.8)
EI = SNN.SpikingSynapse(E, I, :v; σ = 0.5,  p = 0.8)
IE = SNN.SpikingSynapse(I, E, :v; σ = -1.0, p = 0.8)
II = SNN.SpikingSynapse(I, I, :v; σ = -1.0, p = 0.8)
P = [E, I]
C = [EE, EI, IE, II]

SNN.monitor([E, I], [:fire])
for t = 1:1000
    E.I .= 5randn(Ne)
    I.I .= 2randn(Ni)
    SNN.sim!(P, C, 1ms)
end
SNN.raster(P) |> display


text:
using DeformableRegistration
using Images
using Logging
import DeformableRegistration.Distance.ngfDistance

# Normalized gradient field distance
function ngfDistance(referenceImage::regImage,templateImage::regImage,
                     transformedGrid::Array{Float64,1};
                     doDerivative::Bool=false,doHessian::Bool=false,options::regOptions=regOptions(), centeredGrid::Array{Float64,1}=zeros(1)
                     )

    opt = options

    edgeParameterR=options.edgeParameterR
    edgeParameterT=options.edgeParameterT
    debug = false
    useEdgeParameterInNumerator=options.useEdgeParameterInNumerator
    parametricOnly=options.parametricOnly

    if(size(centeredGrid)[1]==1)
      centeredGrid = getCellCenteredGrid(referenceImage).data
    end

    if(checkForOddNumberOfGridPoints(referenceImage.data,transformedGrid))
      error("StaggeredGrids are not supported. Please transform to cell centered first.")
    end

    #c,dTtuple = linearInter2D(templateImage.data, ΩT, mT, deformedGrid,doDerivative=doDerivative)
    if doDerivative
      transformedImage, dX_transformedImage, dY_transformedImage =
          interpolateImage(templateImage,transformedGrid,doDerivative=true)
      dT = spdiagm((dX_transformedImage, dY_transformedImage),[0,prod(size(transformedImage))])

    else
      transformedImage =
          interpolateImage(templateImage,transformedGrid,doDerivative=false)
    end

    Logging.debug("\n\n edgeParameterR: %5e    edgeParameterT %5e \n", edgeParameterR, edgeParameterT)

    m = size(referenceImage.data)
    h = referenceImage.voxelsize
    Rc = Array(referenceImage.data)
    shift = referenceImage.shift

    shortFiniteDiffX = spdiagm((-ones(m[1]-1,1),ones(m[1]-2,1)),[0, 1],m[1]-1,m[1])/(2*h[1])
    #shortFiniteDiffX[1] = shortFiniteDiffX[2]
    #shortFiniteDiffX[end] = shortFiniteDiffX[end-1]

    shortFiniteDiffY = spdiagm((-ones(m[2]-1,1),ones(m[2]-2,1)),[0, 1],m[2]-1,m[2])/(2*h[2])
    #shortFiniteDiffY[1] = shortFiniteDiffY[2]
    #shortFiniteDiffY[end] = shortFiniteDiffY[end-1]

    averageX = spdiagm((0.5*ones(m[1]-2,1),0.5*ones(m[1]-1,1)),[-1, 0],m[1],m[1]-1)
    averageY = spdiagm((ones(m[2]-2,1)*0.5,ones(m[2]-1,1)*0.5),[-1, 0],m[2],m[2]-1)

    G1 = sparse(kron(speye(m[2]),shortFiniteDiffX));
    G2 = sparse(kron(shortFiniteDiffY,speye(m[1])));

    AvgX = sparse(kron(speye(m[2]),averageX));
    AvgY = sparse(kron(averageY,speye(m[1])));

    gradTx  = G1*transformedImage[:]
    gradTy  = G2*transformedImage[:]

    gradRx  = G1*Rc[:]
    gradRy =  G2*Rc[:]

    lengthGT = sqrt.(sum(AvgX * (gradTx.^2) + AvgY * (gradTy.^2), 2) + edgeParameterT^2) #epsilon-norm
    lengthGR = sqrt.(sum(AvgX * (gradRx.^2) + AvgY * (gradRy.^2), 2) + edgeParameterR^2)

    if useEdgeParameterInNumerator
      r1  = sum(AvgX * (gradRx.*gradTx) +
		AvgY * (gradRy.*gradTy) ,2) .+ edgeParameterR*edgeParameterT    # numerator
    else
      r1  = sum(AvgX * (gradRx.*gradTx) +
		AvgY * (gradRy.*gradTy) ,2)    # numerator
    end
    r2  = 1./(lengthGT.*lengthGR)             #  ... and denominator
    rc  = r1 .* r2                            # combine things and finalize

    dFunctionValue = 0
    drc = 0

    if doDerivative
      dr1 = (AvgX * spdiagm(gradRx[:]) * G1) + (AvgY * spdiagm(gradRy[:]) * G2) #(AvgX * gradRx) * AvgX * G1 + (AvgY * gradRy) * AvgX * G2
      dr2 = ( spdiagm(-r2.^2)
              * spdiagm(lengthGR)
              * (spdiagm(1./(2 .* lengthGT)))
              * ((2 * AvgX * spdiagm(gradTx) * G1) + (2 * AvgY * spdiagm(gradTy) * G2))
            )

      if(parametricOnly)
          N = prod(m)
          Q = sparse([centeredGrid[1:N] centeredGrid[N+1:end] ones(N)])
          Q = [Q spzeros(size(Q)[1],size(Q)[2]); spzeros(size(Q)[1],size(Q)[2]) Q]
          dT = dT * Q
      end
      #dr2Partial = reshape(  r1 .* -1 ./ (lengthGR .* lengthGT.^3)  ,(size(r1)))

      #drc = (    spdiag((r2 .* (AvgX * gradRx) + dr2Partial .* (AvgX * gradTx) )[:]) * AvgX * G1
      #         + spdiag((r2 .* (AvgY * gradRy) + dr2Partial .* (AvgY * gradTy) )[:]) * AvgY * G2 )
      drc = spdiagm(r1)*dr2 + spdiagm(r2) * dr1
      drc_times_dT = drc*dT;
      dFunctionValue  = Array(-2.0*prod(h)*rc'*drc_times_dT)

      if doHessian
        d2FunctionValue(x) =  (2*prod(h) .* drc_times_dT' * drc_times_dT) * x   # note the missing minus sign is not a bug!
			else
			  d2FunctionValue = 0
		  end
	else
		d2FunctionValue = 0
	end

	functionValue  = (length(rc)*prod(h) - prod(h) * (rc'*rc))[1]

	if(ndims(dFunctionValue)>1)
		dFunctionValue = vec(dFunctionValue)
	end
	#return functionValue[1], dFunctionValue', d2FunctionValue, r1, 1.-(rc.^2)
	return [functionValue,dFunctionValue,d2FunctionValue] #, drc

end


text:
# import TrajectoryOptimization
include("N_plots.jl")

# Read in the system
model,obj = Dynamics.pendulum!
obj_c = Dynamics.pendulum_constrained![2]
obj_c.u_min[1] = -2
obj_c.u_max[1] = 2

# Params
N = 51
method = :hermite_simpson

# Initial Trajectory
U0 = ones(1,N)
X0 = line_trajectory(obj.x0,obj.xf,N)

X0_rollout = copy(X0)
solver = Solver(model,obj_c,N=N)
rollout!(X0_rollout,U0,solver)

# Knot point test params
Ns = [6,11,15,21,41,51,81,101,201,401,501,801,1001]
dt_truth = 1e-3

obj.tf ./ (Ns.-1)
N_truth = Int(obj.tf / dt_truth + 1)

#####################################
#          UNCONSTRAINED            #
#####################################
# Solver Options
opts = SolverOptions()
opts.use_static = false
opts.cost_tolerance = 1e-4
opts.outer_loop_update_type = :default
opts.iterations = 500

group  = "pendulum/unconstrained"
solver = Solver(model, obj, N=N)
res, stats = solve(solver, U0)

solver.state.constrained
using Juno, Profile
Profile.init(n = 10^7, delay = 0.01)
@profile solve_dircol(solver,X0_rollout,U0)
Juno.profiler()
res_d, stats_d  = solve_dircol(solver,X0,U0)

solver_truth, res_truth,  = run_dircol_truth(model, obj, dt_truth, X0_rollout, U0, group)
time_truth = get_time(solver_truth)
plot(time_truth, res_truth.X')

err_mid, eterm_mid, stats_mid = run_Ns(model, obj, Ns, :midpoint)
err_rk3, eterm_rk3, stats_rk3 = run_Ns(model, obj, Ns, :rk3)
err_foh, eterm_foh, stats_foh = run_Ns(model, obj, Ns, :rk3_foh)
err_rk4, eterm_rk4, stats_rk4 = run_Ns(model, obj, Ns, :rk4)

save_data(group)

plot_stat("error",group)
plot_stat("error_final",group,yscale=:log10)
plot_stat("runtime",group)
plot_stat("iterations",group)

#####################################
#            CONSTRAINED            #
#####################################
opts = SolverOptions()
opts.cost_tolerance = 1e-4
opts.cost_tolerance_intermediate = 1e-4
opts.constraint_tolerance = 1e-3
opts.outer_loop_update_type = :individual
opts.use_static = false
opts.constraint_decrease_ratio = .25

group = "pendulum/constrained"
solver_truth, res_truth,  = run_dircol_truth(model, obj_c, dt_truth, X0_rollout, U0, group)
time_truth = get_time(solver_truth)
plot(res_truth.X')

Ns = [21,41,51,81,101,201,401,501,801,1001]
err_mid, eterm_mid, stats_mid = run_Ns(model, obj_c, Ns, :midpoint)
err_rk3, eterm_rk3, stats_rk3 = run_Ns(model, obj_c, Ns, :rk3)
err_foh, eterm_foh, stats_foh = run_Ns(model, obj_c, Ns, :rk3_foh)
err_rk4, eterm_rk4, stats_rk4 = run_Ns(model, obj_c, Ns, :rk4)

save_data(group)

plot_stat("error",group)
plot_stat("error_final",group,yscale=:log10)
plot_stat("runtime",group,legend=:bottomright)
plot_stat("iterations",group)
plot_stat("c_max",group,yscale=:log10)

opts = SolverOptions()
opts.cost_tolerance = 1e-16
opts.cost_tolerance_intermediate = 1e-16
opts.constraint_tolerance = 1e-16
opts.use_static = false
opts.iterations = 250
opts.iterations_outerloop = 10

solver = solver(model, obj_c, opts=opts, N=51, integration=:rk3)

#####################################
#            INFEASIBLE             #
#####################################
# Solver Options
opts = SolverOptions()
opts.cost_tolerance = 1e-4
opts.cost_tolerance_intermediate = 1e-1
opts.constraint_tolerance = 1e-3
opts.outer_loop_update_type = :individual
opts.constraint_decrease_ratio = .25
opts.use_static = false
opts.resolve_feasible = false

solver_truth, res_truth,  = run_dircol_truth(model, obj_c, dt_truth, X0, U0, "pendulum/infeasible")
time_truth = get_time(solver_truth)
plot(res_truth.X')

obj_c

Ns = [21,51,81,101,201,401,501,801,1001]
err_mid, eterm_mid, stats_mid = run_Ns(model, obj_c, Ns, :midpoint, infeasible=true)
err_rk3, eterm_rk3, stats_rk3 = run_Ns(model, obj_c, Ns, :rk3, infeasible=true)
err_foh, eterm_foh, stats_foh = run_Ns(model, obj_c, Ns, :rk3_foh, infeasible=true)
err_rk4, eterm_rk4, stats_rk4 = run_Ns(model, obj_c, Ns, :rk4, infeasible=true)

save_data("pendulum/infeasible")

plot_stat("error",group)
plot_stat("error_final",group,yscale=:log10)
plot_stat("runtime",group,legend=:bottomright)
plot_stat("iterations",group)
plot_stat("c_max",group,yscale=:log10)

solver = Solver(model, obj_c, opts=opts, N=51, integration=:rk3_foh)
disable_logging(Logging.Debug)
solver.opts.verbose = true
solver.opts.live_plotting = true
solver.opts.cost_tolerance_intermediate = 1e-6
solve(solver,X0,U0)


text:
using Zygote, LinearAlgebra
using GenericLinearAlgebra: _eigvals!

function test_eigmin(A::AbstractMatrix)
    vals = _eigvals!(copy(A))
    ev = minimum(real.(vals))
    return ev
end

A = rand(4, 4); A = A + A';

@show eigmin(A)
@show test_eigmin(A)
@show Zygote.gradient(test_eigmin, A)

max_stars_repo_name: ericphanson/Attempt-diff-eigmin

max_issues_repo_name: ericphanson/Attempt-diff-eigmin

max_forks_repo_name: ericphanson/Attempt-diff-eigmin


text:
#####
##### Gradient AD implementation using Flux
#####

import .Flux

struct FluxGradientLogDensity{L} <: ADGradientWrapper
    ℓ::L
end

"""
$(SIGNATURES)

Gradient using algorithmic/automatic differentiation via Flux.
"""
ADgradient(::Val{:Flux}, ℓ) = FluxGradientLogDensity(ℓ)

Base.show(io::IO, ∇ℓ::FluxGradientLogDensity) = print(io, "Flux AD wrapper for ", ∇ℓ.ℓ)

function logdensity_and_gradient(∇ℓ::FluxGradientLogDensity, x::AbstractVector{T}) where {T}
    @unpack ℓ = ∇ℓ
    y, back = Flux.Tracker.forward(x -> logdensity(ℓ, x), x)
    yval = Flux.Tracker.data(y)
    # work around https://github.com/FluxML/Flux.jl/issues/497
    z = T <: Real ? zero(T) : 0.0
    S = typeof(z + 0.0)
    S(yval)::S, (S.(first(Flux.Tracker.data.(back(1)))))::Vector{S}
end


text:
isdefined(Base, :__precompile__) && __precompile__(false)

module StateSpaceRoutines

    using Distributions: Distribution, MvNormal, pdf, Weights, sample, logpdf
    using Tracker
    using LinearAlgebra, Statistics
    using Roots: fzero
    using HDF5, JLD2
    using Distributed
    using SharedArrays

    export

        # filters/kalman_filter.jl
        init_stationary_states, kalman_filter, chand_recursion, kalman_likelihood,
        kalman_filtered_states,

        # filters/tempered_particle_filter
        tempered_particle_filter, initialize_state_draws,
        solve_inefficiency, mutation, weight_kernel!, next_φ,correction!,
        selection!, mutation!, update_c,

        # smoothers/
        hamilton_smoother, koopman_smoother, koopman_disturbance_smoother,
        carter_kohn_smoother, durbin_koopman_smoother

    const VERBOSITY = Dict(:none => 0, :low => 1, :high => 2)

    include("filters/util.jl")
    include("filters/kalman_filter.jl")
    include("filters/chand_recursion.jl")
    include("filters/tempered_particle_filter/util.jl")
    include("filters/tempered_particle_filter/initialization.jl")
    include("filters/tempered_particle_filter/correction.jl")
    include("filters/tempered_particle_filter/selection.jl")
    include("filters/tempered_particle_filter/mutation.jl")
    include("filters/tempered_particle_filter/tempered_particle_filter.jl")

    include("smoothers/util.jl")
    include("smoothers/hamilton_smoother.jl")
    include("smoothers/koopman_smoother.jl")
    include("smoothers/carter_kohn_smoother.jl")
    include("smoothers/durbin_koopman_smoother.jl")

end


text:
module TracePrecompiles

function run_julia(cmd)
    current_proj = unsafe_string(Base.JLOptions().project)
    run(`$(Base.julia_cmd()) --project=$(current_proj) --startup-file=no $(cmd)`)
end

function trace_compiles(package, trace_file, outputfile)
    tdir = mktempdir(; cleanup=true)
    trace_out = joinpath(tdir, "precompiles.jl")
    run_julia(`--trace-compile=$(trace_out) $trace_file`)
    prs = joinpath(@__DIR__, "process-prs.jl")
    run_julia(` $prs $(package) $(trace_out) $(outputfile)`)
end

end


text:
export BrukerDatasetStore, BrukerStore, findBrukerFiles, getNewCalibNum

struct BrukerDatasetStore <: DatasetStore
  path::String
end

if ispath("/opt/mpidata")
  const BrukerStore = BrukerDatasetStore("/opt/mpidata")
end

path(e::Experiment{BrukerDatasetStore}) = joinpath( path(e.study), string(e.num) )
path(s::Study{BrukerDatasetStore}, numExp::Integer) = joinpath(path(s),string(numExp))
iscalib(e::Experiment{BrukerDatasetStore}) = _iscalib(path(e))

###  Implementations of abstract interfaces ###

readonly(::BrukerDatasetStore) = true

studydir(d::BrukerDatasetStore) = d.path
calibdir(d::BrukerDatasetStore) = error("BrukerDatasetStore has no calibdir")

function getStudy(d::BrukerDatasetStore, studyfolder::String)
  study = nothing
  studypath = joinpath(d.path,studyfolder)
  if !ishidden(studypath) && isdir(studypath)
    w = split(studyfolder,'_')
    if length(w) >= 5 && length(w[1])==8 # only these can be study folders
      # w_ = w[1:end-2]
      # date = w[1]
      # date = string(date[1:4],"/",date[5:6],"/",date[7:8])

      w = split(studyfolder,'_')
      dateStr = w[1]
      timeStr = w[2]
      date = DateTime(string(dateStr[1:4],"-",dateStr[5:6],"-",dateStr[7:8],"T",
			   timeStr[1:2],":",timeStr[3:4],":",timeStr[5:6]))

      j = JcampdxFile()
      subjfile = string(studypath,"/subject")
      if isfile(subjfile)
        read(j,string(studypath,"/subject"),maxEntries=14) #magic number...
        name = latin1toutf8(j["SUBJECT_study_name"])
        # name = string(latin1toutf8(j["SUBJECT_name_string"]),
        #              "_",latin1toutf8(j["SUBJECT_study_name"]),
        #              "_",latin1toutf8(j["SUBJECT_study_nr"]))
        s1 = latin1toutf8(j["SUBJECT_id"])
	      s2 = latin1toutf8(j["SUBJECT_name_string"])
        subject = (s1 == s2) ? s1 : s1*s2
      else
        # Workaround if no subject file is present => use first dataset
        # and derive the study from the Brukerfile
        r = readdir(studypath)

        found = false
        for file in r

          if tryparse(Int64,file) != nothing
            b = BrukerFileFast(joinpath(studypath, file ))
            name = studyName(b)
            subject = experimentSubject(b)
            found = true
            break
          end
        end
        if !found
          return nothing
        end
      end
      study = Study(d, name, studyfolder, date, subject)
    end
  end
  return study
end

function getExperiments(s::Study{BrukerDatasetStore})
  @static if Sys.isunix()
    files = findBrukerFiles(path(s),1,1)
  else
    files = findBrukerFiles(path(s))
  end

  experiments = Experiment[]

  for file in files
    prefix, ext = splitext(splitpath(file)[end])
    num = tryparse(Int64, prefix)
    if num != nothing
      #try
      exp = getExperiment(s,num)
      push!(experiments, exp)
      #catch e
      #  @debug "" e
      #end
    end
  end
  return experiments
end

@static if Sys.isunix()
  function findBrukerFiles(path::AbstractString, mindepth::Int=1, maxdepth::Int=2)
    candidatePaths = split(read(`find $path -maxdepth $maxdepth -mindepth $mindepth -type d`,String),"\n")[1:end-1]
    mask = zeros(Bool,length(candidatePaths))
    for (i,candidatePath) in enumerate(candidatePaths)
      if isfile(joinpath(candidatePath,"acqp")) &&
         isfile(joinpath(candidatePath,"method")) &&
         isfile(joinpath(candidatePath,"visu_pars"))
        mask[i] = true
      end
    end
    return String.(candidatePaths[mask])
  end
else
  function findBrukerFiles(path::AbstractString)
    files = readdir(path)
    bfiles = String[]
    for file in files
      if isdir(joinpath(path,file))
       try
        if isfile(joinpath(path,file,"acqp")) &&
           isfile(joinpath(path,file,"method")) &&
           isfile(joinpath(path,file,"visu_pars"))
          push!(bfiles, joinpath(path,file))
        else
          rfiles = findBrukerFiles(joinpath(path,file))
          if rfiles != nothing && length(rfiles) > 0
            push!(bfiles, rfiles...)
          end
        end
       catch
        continue
       end
      end
    end
  return bfiles
  end
end


text:
result = [
    2.32551,
    -0.155293,
    0.469214,
    -0.0870457,
    0.271364,
    -0.222398,
    1.01423,
    -0.179009,
    1.01724,
    -0.152972,
    0.216039,
    -0.251637,
    1.05034,
    -0.231289,
    0.245388,
    -0.231289,
    1.03371,
    -0.258872,
    1.03256,
    -0.262519,
    1.04748,
    -0.259143,
    1.0535,
    -0.266484,
    1.04711,
    -0.267177,
    1.02131,
    -0.280381,
]

p_gen_matpower_3bus = [20.3512373930753, 100.0, 100.0]
q_gen_matpower_3bus = [45.516916781567232, 10.453799727283879, -31.992561631394636]

pf_sys5_re = PSB.build_system(PSB.PSITestSystems, "c_sys5_re"; add_forecasts=false)
remove_component!(Line, pf_sys5_re, "1")
remove_component!(Line, pf_sys5_re, "2")
br = get_component(Line, pf_sys5_re, "6")
PSY.set_x!(br, 20.0)
PSY.set_r!(br, 2.0)

@testset "NLsolve Power Flow testing" begin
    # This is a negative test. The data passed for sys5_re is known to be infeasible.
    @test_logs(
        (:error, "The powerflow solver returned convergence = false"),
        match_mode = :any,
        @test !run_powerflow!(pf_sys5_re, finite_diff=true)
    )
    #Compare results between finite diff methods and Jacobian method
    res_finite_diff = run_powerflow(
        PSB.build_system(PSB.PSITestSystems, "c_sys14"; add_forecasts=false),
        finite_diff=true,
    )
    res_jacobian =
        run_powerflow(PSB.build_system(PSB.PSITestSystems, "c_sys14"; add_forecasts=false))
    @test LinearAlgebra.norm(
        res_finite_diff["bus_results"].Vm - res_jacobian["bus_results"].Vm,
    ) <= 1e-6
    @test run_powerflow!(
        PSB.build_system(PSB.PSITestSystems, "c_sys14"; add_forecasts=false),
        finite_diff=true,
        method=:newton,
    )

    sys = PSB.build_system(PSB.PSITestSystems, "c_sys14"; add_forecasts=false)
    branch = first(PSY.get_components(Line, sys))
    dyn_branch = DynamicBranch(branch)
    add_component!(sys, dyn_branch)
    @test dyn_pf = run_powerflow!(sys)
    dyn_pf = run_powerflow(sys)
    @test LinearAlgebra.norm(dyn_pf["bus_results"].Vm - res_jacobian["bus_results"].Vm) <=
          1e-6

    sys = PSB.build_system(PSB.PSITestSystems, "c_sys14"; add_forecasts=false)
    line = get_component(Line, sys, "Line4")
    PSY.set_available!(line, false)
    run_powerflow!(sys)
    @test PSY.get_active_power_flow(line) == 0.0
    test_bus = get_component(PSY.Bus, sys, "Bus 4")
    @test isapprox(PSY.get_magnitude(test_bus), 1.002; atol=1e-3)

    sys = PSB.build_system(PSB.PSITestSystems, "c_sys14"; add_forecasts=false)
    line = get_component(Line, sys, "Line4")
    PSY.set_available!(line, false)
    res = run_powerflow(sys)
    @test res["flow_results"].P_from_to[4] == 0.0
    @test res["flow_results"].P_to_from[4] == 0.0

    sys = PSB.build_system(PSB.PSSETestSystems, "psse_240_case_renewable_sys")
    @test run_powerflow!(sys)

    sys_3bus = PSB.build_system(PSB.PSSETestSystems, "psse_3bus_gen_cls_sys")
    bus_103 = get_component(PSY.Bus, sys_3bus, "BUS 3")
    fix_shunt = PSY.FixedAdmittance("FixAdmBus3", true, bus_103, 0.0 + 0.2im)
    add_component!(sys_3bus, fix_shunt)
    df = run_powerflow(sys_3bus)
    @test isapprox(df["bus_results"].P_gen, p_gen_matpower_3bus, atol=1e-4)
    @test isapprox(df["bus_results"].Q_gen, q_gen_matpower_3bus, atol=1e-4)
end

@testset "Test 240 Case PSS/e results" begin
    file = joinpath(
        TEST_FILES_DIR,
        "test_data",
        "WECC240_v04_DPV_RE20_v33_6302_xfmr_DPbuscode_PFadjusted_V32_noRemoteVctrl.raw",
    )
    system = System(
        file,
        bus_name_formatter=x -> strip(string(x["name"])) * "-" * string(x["index"]),
        runchecks=false,
    )

    pf_bus_result_file = joinpath(TEST_FILES_DIR, "test_data", "pf_bus_results.csv")
    pf_gen_result_file = joinpath(TEST_FILES_DIR, "test_data", "pf_gen_results.csv")

    pf = run_powerflow!(system)
    @test pf
    pf_result_df = run_powerflow(system)

    v_diff, angle_diff, number = psse_bus_results_compare(pf_bus_result_file, pf_result_df)
    p_diff, q_diff, names = psse_gen_results_compare(pf_gen_result_file, system)

    base_power = get_base_power(system)
    @test norm(v_diff, Inf) < DIFF_INF_TOLERANCE
    @test norm(v_diff, 2) / length(v_diff) < DIFF_L2_TOLERANCE
    @test norm(angle_diff, Inf) < DIFF_INF_TOLERANCE
    @test norm(angle_diff, 2) / length(angle_diff) < DIFF_L2_TOLERANCE
    @test norm(p_diff, Inf) < DIFF_INF_TOLERANCE * base_power
    @test norm(p_diff, 2) / length(p_diff) < DIFF_L2_TOLERANCE
    @test sum(q_diff) < DIFF_INF_TOLERANCE * base_power
    @test norm(q_diff, 2) / length(q_diff) < DIFF_L2_TOLERANCE
end


text:
"""
    normalizedpath(ht, vocab_hash, in_word)

Generate the path through a Huffman encoding for in_word.  Returns nodepath
(the series of nodes from the root of the Huffman tree
to the word leaf), and branchpath (the binary representation of whether to
branch left or right from the node).
"""
function normalizedpath(ht::HuffmanTree, vocab_hash::Dict, in_word::String)
    nodepath, branchpath = _rootpath(ht, vocab_hash, in_word)
    nodepath = nodepath .- length(vocab_hash)
    return nodepath, branchpath
end

# generate the path from leaf to root; because leaf nodes all have index less
# than internal nodes, the path is offset relative to 1 by the number of leaves
function _rootpath(ht::HuffmanTree, vocab_hash::Dict, in_word::String)
    nodepath = Array{Int, 1}()
    branchpath = Array{Int, 1}()
    leaf = vocab_hash[in_word]
    lparent = ht.nodeparent[leaf]
    lbin = ht.branch[leaf]
    while lparent > 0
        push!(nodepath, lparent)
        push!(branchpath, lbin)
        lbin = ht.branch[lparent]
        lparent = ht.nodeparent[lparent]
    end
    return nodepath, branchpath
end

"""
    allpaths(ht, vocab_hash, vocab)

Generate all paths through a Huffman encoding for every word in a vocabulary
hash, by using a vocabulary array.  Returns nodepaths (the series of nodes
from the root of the Huffman tree to the word leaf), and branchpaths (the binary
representation of whether to branch left or right from the node).
"""
function allpaths(ht::HuffmanTree, vocab_hash::Dict, vocab::Array)
    nodepaths = Array{Array{Int, 1}, 1}()
    branchpaths = Array{Array{Int, 1}, 1}()
    for key in 1:length(vocab_hash)
        wordnodes, wordbranches = normalizedpath(ht, vocab_hash, vocab[key])
        push!(nodepaths, wordnodes)
        push!(branchpaths, wordbranches)
    end
    return nodepaths, branchpaths
end


text:
# Container data structure (structure-of-arrays style) for DG elements
struct ElementContainer3D{NVARS, POLYDEG} <: AbstractContainer
  u::Array{Float64, 5}                   # [variables, i, j, k, elements]
  u_t::Array{Float64, 5}                 # [variables, i, j, k, elements]
  u_tmp2::Array{Float64, 5}              # [variables, i, j, k, elements]
  u_tmp3::Array{Float64, 5}              # [variables, i, j, k, elements]
  inverse_jacobian::Vector{Float64}      # [elements]
  node_coordinates::Array{Float64, 5}    # [orientation, i, j, k, elements]
  surface_ids::Matrix{Int}               # [direction, elements]
  surface_flux_values::Array{Float64, 5} # [variables, i, j, direction, elements]
  cell_ids::Vector{Int}                  # [elements]
end

function ElementContainer3D{NVARS, POLYDEG}(capacity::Integer) where {NVARS, POLYDEG}
  # Initialize fields with defaults
  n_nodes = POLYDEG + 1
  u = fill(NaN, NVARS, n_nodes, n_nodes, n_nodes, capacity)
  u_t = fill(NaN, NVARS, n_nodes, n_nodes, n_nodes, capacity)
  # u_rungakutta is initialized to non-NaN since it is used directly
  u_tmp2 = fill(0.0, NVARS, n_nodes, n_nodes, n_nodes, capacity)
  u_tmp3 = fill(0.0, NVARS, n_nodes, n_nodes, n_nodes, capacity)
  inverse_jacobian = fill(NaN, capacity)
  node_coordinates = fill(NaN, 3, n_nodes, n_nodes, n_nodes, capacity) # 3 = ndims
  surface_ids = fill(typemin(Int), 2 * 3, capacity) # 3 = ndims
  surface_flux_values = fill(NaN, NVARS, n_nodes, n_nodes, 2 * 3, capacity) # 3 = ndims
  cell_ids = fill(typemin(Int), capacity)

  elements = ElementContainer3D{NVARS, POLYDEG}(u, u_t, u_tmp2, u_tmp3, inverse_jacobian, node_coordinates,
                                                surface_ids, surface_flux_values, cell_ids)

  return elements
end

# Return number of elements
nelements(elements::ElementContainer3D) = length(elements.cell_ids)

# Container data structure (structure-of-arrays style) for DG interfaces
struct InterfaceContainer3D{NVARS, POLYDEG} <: AbstractContainer
  u::Array{Float64, 5}      # [leftright, variables, i, j, interfaces]
  neighbor_ids::Matrix{Int} # [leftright, interfaces]
  orientations::Vector{Int} # [interfaces]
end

function InterfaceContainer3D{NVARS, POLYDEG}(capacity::Integer) where {NVARS, POLYDEG}
  # Initialize fields with defaults
  n_nodes = POLYDEG + 1
  u = fill(NaN, 2, NVARS, n_nodes, n_nodes, capacity)
  neighbor_ids = fill(typemin(Int), 2, capacity)
  orientations = fill(typemin(Int), capacity)

  interfaces = InterfaceContainer3D{NVARS, POLYDEG}(u, neighbor_ids, orientations)

  return interfaces
end

# Return number of interfaces
ninterfaces(interfaces::InterfaceContainer3D) = length(interfaces.orientations)

# Container data structure (structure-of-arrays style) for DG boundaries
struct BoundaryContainer3D{NVARS, POLYDEG} <: AbstractContainer
  u::Array{Float64, 5}                # [leftright, variables, i, j, boundaries]
  neighbor_ids::Vector{Int}           # [boundaries]
  orientations::Vector{Int}           # [boundaries]
  neighbor_sides::Vector{Int}         # [boundaries]
  node_coordinates::Array{Float64, 4} # [orientation, i, j, elements]
end

function BoundaryContainer3D{NVARS, POLYDEG}(capacity::Integer) where {NVARS, POLYDEG}
  # Initialize fields with defaults
  n_nodes = POLYDEG + 1
  u = fill(NaN, 2, NVARS, n_nodes, n_nodes, capacity)
  neighbor_ids = fill(typemin(Int), capacity)
  orientations = fill(typemin(Int), capacity)
  neighbor_sides = fill(typemin(Int), capacity)
  node_coordinates = fill(NaN, 3, n_nodes, n_nodes, capacity) # 3 = ndims

  boundaries = BoundaryContainer3D{NVARS, POLYDEG}(u, neighbor_ids, orientations, neighbor_sides,
                                                   node_coordinates)

  return boundaries
end

# Return number of boundaries
nboundaries(boundaries::BoundaryContainer3D) = length(boundaries.orientations)

# Container data structure (structure-of-arrays style) for DG L2 mortars
# Positions/directions for orientations = 1, large_sides = 2:
# mortar is orthogonal to x-axis, large side is in positive coordinate direction wrt mortar
#   /----------------------------\  /----------------------------\
#   |             |              |  |                            |
#   | upper, left | upper, right |  |                            |
#   |      3      |      4       |  |                            |
#   |             |              |  |           large            |
#   |-------------|--------------|  |             5              |
# z |             |              |  |                            |
#   | lower, left | lower, right |  |                            |
# ^ |      1      |      2       |  |                            |
# | |             |              |  |                            |
# | \----------------------------/  \----------------------------/
# |
# ⋅----> y
# Left and right are always wrt to a coordinate direction:
# * left is always the negative direction
# * right is always the positive direction
#
# Left and right are used *both* for the numbering of the mortar faces *and* for the position of the
# elements with respect to the axis orthogonal to the mortar.
struct L2MortarContainer3D{NVARS, POLYDEG} <: AbstractContainer
  u_upper_left ::Array{Float64, 5} # [leftright, variables, i, j, mortars]
  u_upper_right::Array{Float64, 5} # [leftright, variables, i, j, mortars]
  u_lower_left ::Array{Float64, 5} # [leftright, variables, i, j, mortars]
  u_lower_right::Array{Float64, 5} # [leftright, variables, i, j, mortars]
  neighbor_ids::Matrix{Int}  # [position, mortars]
  # Large sides: left -> 1, right -> 2
  large_sides::Vector{Int}   # [mortars]
  orientations::Vector{Int}  # [mortars]
end

function L2MortarContainer3D{NVARS, POLYDEG}(capacity::Integer) where {NVARS, POLYDEG}
  # Initialize fields with defaults
  n_nodes = POLYDEG + 1
  u_upper_left  = fill(NaN, 2, NVARS, n_nodes, n_nodes, capacity)
  u_upper_right = fill(NaN, 2, NVARS, n_nodes, n_nodes, capacity)
  u_lower_left  = fill(NaN, 2, NVARS, n_nodes, n_nodes, capacity)
  u_lower_right = fill(NaN, 2, NVARS, n_nodes, n_nodes, capacity)
  neighbor_ids = fill(typemin(Int), 5, capacity)
  large_sides = fill(typemin(Int), capacity)
  orientations = fill(typemin(Int), capacity)

  l2mortars = L2MortarContainer3D{NVARS, POLYDEG}(u_upper_left, u_upper_right,
                                                  u_lower_left, u_lower_right,
                                                  neighbor_ids, large_sides, orientations)

  return l2mortars
end

# Return number of L2 mortars
nmortars(l2mortars::L2MortarContainer3D) = length(l2mortars.orientations)

# Allow printing container contents
function Base.show(io::IO, c::L2MortarContainer3D{NVARS, POLYDEG}) where {NVARS, POLYDEG}
  println(io, '*'^20)
  for idx in CartesianIndices(c.u_upper_left)
    println(io, "c.u_upper_left[$idx] = $(c.u_upper_left[idx])")
  end
  for idx in CartesianIndices(c.u_upper_right)
    println(io, "c.u_upper_right[$idx] = $(c.u_upper_right[idx])")
  end
  for idx in CartesianIndices(c.u_lower_left)
    println(io, "c.u_lower_left[$idx] = $(c.u_lower_left[idx])")
  end
  for idx in CartesianIndices(c.u_lower_right)
    println(io, "c.u_lower_right[$idx] = $(c.u_lower_right[idx])")
  end
  println(io, "transpose(c.neighbor_ids) = $(transpose(c.neighbor_ids))")
  println(io, "c.large_sides = $(c.large_sides)")
  println(io, "c.orientations = $(c.orientations)")
  println(io, '*'^20)
end


text:
module Polymer

using LinearAlgebra
using ArgCheck
using REPL: symbol_latex
using LaTeXStrings
using YAML

include("utils.jl")
export
    unicodesymbol2string,
    reverse_dict

include("parameters.jl")
export
    AbstractParameter,
    PolymerParameter
export
    χParam,
    NParam,
    χNParam,
    fParam,
    ϕParam,
    RgParam,
    CParam,
    bParam,
    αParam,
    τParam,
    DEFAULT_PARAMETERS
export
    description,
    value_type,
    variable_symbol,
    as_variable_name,
    as_ascii_label,
    as_plot_label

include("chiN.jl")
export
    AbstractχNMatrix,
    χNMatrix

include("types.jl")
export
    PolymerSystemType,
    NeatPolymer, PolymerBlend,
    PolymerSolution
export
    SpaceDimension,
    D1, D2, D3
export
    ConfinementType,
    BulkConfinement,
    BoxConfinement,
    SlabConfinement,
    DiskConfinement,
    SphereConfinement,
    CylinderConfinement
export
    ChargedType,
    Neutral,
    SmearedCharge,
    DiscreteCharge
export
    BlockEnd,
    FreeEnd,
    BranchPoint,
    PolymerBlock
export
    AbstractSpecie,
    KuhnSegment
export
    AbstractMolecule,
    SmallMolecule,
    AbstractPolymer,
    BlockCopolymer,
    RandomCopolymer,
    AlternatingCopolymer,
    Particle,
    GiantMolecule
export
    AbstractComponent,
    Component,
    AbstractSystem,
    PolymerSystem
export
    isconfined,
    ischarged,
    isfreeblockend,
    multicomponent,
    ncomponents,
    specie,
    species,
    nspecies,
    nblocks,
    systemtype

include("systems.jl")
export
    homopolymer_chain,
    diblock_chain,
    solvent,
    AB_system,
    AB_A_system,
    AB_S_system

include("make.jl")
export
    ObjType
export
    load_config,
    make

end # module


text:
# Macros used to express ISPC concepts in Julia code.

using Base.Meta

function ispc_esc(expr)
    # println(expr)
    return esc(expr)
    # All ISPC entry points are escaped as @fastmath.
    # This is to prevent Julia from inserting a number of runtime
    # checks into the lowered AST that are either untranslateable
    # or unnecessary because they'll be done by the ISPC compiler.
    # Note: this has NO influence on whether the ISPC code will be
    # compiled with the '--fastmath' option!
    # return Base.FastMath.make_fastmath(esc(expr))
end

"""
The ISPC `foreach` construct:

@foreach(array) do index
    ...
end

@foreach(:active) do index
    ...
end
"""
macro foreach(func_expr, array_expr...)
    signature, block = func_expr.args
    # Turn the lambda into a let block:
    foreach_id = gensym(:foreach)
    code = ispc_esc(quote
        $(Expr(:meta, :ispc, foreach_id, :foreach, array_expr))
        let $([:($a = ISPC.foreachindex($i, $s))
                    for (i, (s, a)) in enumerate(
                            zip(array_expr, signature.args))]...)
            # We use the iteration expression as an argument to
            # foreachindex() to pull the necessary variables into
            # the closure. The :meta statement does not do that.
            $block
        end
        $(Expr(:meta, :ispc, foreach_id))
        nothing # otherwise the :meta above becomes a return value
    end)
end

# That function will not be called, its sole purpose
# is to give a type to foreach indices so that type
# inference can work on them:
@noinline foreachindex(args...) = 1

"""
The ISPC `unmasked` construct:

@unmasked begin
    ...
end
"""
macro unmasked(block)
    unmasked_id = gensym(:unmasked)
    ispc_esc(quote
        $(Expr(:meta, :ispc, unmasked_id, :unmasked))
        $(block)
        $(Expr(:meta, :ispc, unmasked_id))
        nothing # otherwise the :meta above becomes a return value
    end)
end

"""
The ISPC `coherent` construct used to mark a test expression
as likely to evaluate uniformly over a whole gang:

if x > 0 @coherent
    ...
elseif x == 0 @coherent
    ...
else
    ...
end

while x >= 0 @coherent
    ...
end
"""
macro coherent()
    return Expr(:meta, :ispc, :coherent)
end

"""
Marks a function that may contain ISPC kernels:

@ispc function foo(...)
    ...
    @kernel() do
        ...
    end
    ...
end
"""
macro ispc(fdef::Expr)

    # What kind of function definition do we have?
    func = macroexpand(fdef)

    if func.head == :function
        # That's a plain function
        return esc(ISPC.make_ispc_main(expand(func)))

    elseif func.head == :stagedfunction
        # That's a @generated function.

        # This case is a bit harder to handle. We need to call the
        # original @generated function in order to know the actual
        # function body where we want to extract the kernels, and
        # to do that we need to know the argument types.

        # So we'll have to declare the original @generated function,
        # declare a new @generated function to catch argument types,
        # obtain the actual body from the original function, and run
        # the extraction routine in the new @generated function.

        # First let's escape the name of the @generated block:
        stagedname, stageddef = rename_fdef(func)

        # Now build another @generated function to catch the argument types:
        funcspec = func.args[1]
        argspec = func.args[1].args[2:end]
        argnames = map(get_argname, argspec)

        ispcbody = quote
            # println($stagedname)
            # println(methods($stagedname))
            # m = methods($stagedname).defs
            # println(m)
            # println(methods($stagedname).defs.func)
            # println(methods($stagedname).defs.func.env)
            # sleep(0.5)
            # ISPC.print_lambda(methods($stagedname).defs.func.code.ast)
            # sleep(0.5)

            flowered = code_typed($stagedname, Tuple{$(argnames...)})
            body = flowered[1]

            # # Get the @generated body for the argument types we have:
            # body = methods($stagedname).defs.func($(argnames...))
            # println(body)

            # # Create a plain function with the same signature and the
            # # body we just got, otherwise the code won't be lowered
            # # right because the arguments aren't there:
            # funcspec = $(Expr(:quote, funcspec))

            # # Run it through make_ispc_main():
            # func = Expr(:function, funcspec, body)
            # println(func)
            ispcmain = ISPC.make_ispc_main(body)
            # println(ispcmain)

            # Return the *body* of the transformed function:
            body = ispcmain.args[3]
            quote
                $(body.args...) # splatting avoids "misplaced return" error
            end
        end
        ispcfunc = Expr(:stagedfunction, funcspec, ispcbody)
        # The @ispc macro returns both the original, escaped @generated function
        # and the ispc @generated function:
        esc(quote
            import ISPC
            $stageddef
            @Base.__doc__($ispcfunc)
        end)
    end
end

"""
Defines the entry point to an ISPC kernel.
Must be used inside an @ispc function.

@kernel([compile-opts]) do
    ...
end
"""
macro kernel(args...)
    if length(args) == 2
        lambda, opts = args
    elseif length(args) == 1
        lambda, = args
        opts = :(``)
    end

    # Generate a unique identifier for this kernel.
    # Wrap it in a Val{} so that it can be used for
    # multiple dispatch, and also to prevent it from
    # being turned into a GlobalRef when the AST is
    # lowered.
    identifier = Val{gensym()}

    # Kernels don't return anything:
    tree = lambda.args[2]
    push!(tree.args, :nothing)

    # insert!(tree.args, 1, Expr(:boundscheck, false))
    # push!(tree.args, Expr(:boundscheck, :pop))

    return ispc_esc(quote
        ISPC.ispc_kernel($identifier, $lambda, $opts)
    end)
end

@noinline function ispc_kernel(args...)
    # Note: this should not be called, it's just a tag that
    # will be picked up when the kernels are extracted.
    # Provide a helpful error message:
    error("Functions containing ISPC kernels must be marked with @ispc")
end


text:
module TestStats
    using Base.Test
    using DataArrays

    autocor(DataArray([1, 2, 3, 4, 5]))

    @assert isequal(xtabs([1, 2, 2, 2, 3]), Dict([(2, 3), (3, 1), (1, 1)]))
end


text:
### A Pluto.jl notebook ###
# v0.17.1

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local iv = try Base.loaded_modules[Base.PkgId(Base.UUID("6e696c72-6542-2067-7265-42206c756150"), "AbstractPlutoDingetjes")].Bonds.initial_value catch; b -> missing; end
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : iv(el)
        el
    end
end

# ╔═╡ dc0623d0-3438-4bf3-9bd0-ea2b35589b7e
begin
    import Pkg
    # activate the shared project environment
    Pkg.activate(Base.current_project())
	#Pkg.add("Plots")
    using Omega, RCall, Plots, Distributions
end

# ╔═╡ d5bac81e-e0e1-40f0-b889-529f8ed7719c
using PlutoUI

# ╔═╡ 01271358-45ae-11ec-30f1-8521b73aa741
md"# Constrained Flows

Deep neural networks can be used to cosntruct very flexible families of probability distributions.  However, these families might be __too__ flexible in the sense that they contain distributions which are inconsistent with facts that we know.

We can use the distributional inference features of Omega to incorporate certain kinds of domain knowledge into these families.  For example, we might know that the distribution has exactly two modes."

# ╔═╡ 6353b3f7-e6e5-46e4-80e1-272d8b9286bd
md"First let's setup a Bayesian normalizing flow from.  A normalizing flow is a real values probability distribution, which among other things, allows us to easily compute the pdf of any value"

# ╔═╡ 836d9f48-5f7a-4963-955d-9f4e4a7c8833
function g(k, y)
	if k == 1
		1
	elseif k == 2
		log(y/(1-y))
	elseif iseven(k)
		(y - 0.5)
	else
		(y - 0.5) * log(y/(1-y))
	end
end

# ╔═╡ e7b83751-0720-4cda-8efa-3f14ee0c98e3
@bind x1 Slider(-10:0.1:10)

# ╔═╡ b565412f-3e41-44f3-879e-fd7aace7a1b6
@bind x2 Slider(-10:0.1:10)

# ╔═╡ 2a0e9e4f-cdf7-4ac7-9dae-9a3249c61dd7
@bind x3 Slider(-10:0.1:10)

# ╔═╡ af3ebb43-a593-4a07-8b63-869c139f6f05
@bind x4 Slider(-10:0.1:10)

# ╔═╡ f21c55ce-4539-437a-a667-2f489134a4f9
@bind x5 Slider(-10:0.1:10)

# ╔═╡ a1f683a2-f9d3-4829-abb0-6fe6b831b1e8
@bind x6 Slider(-10:0.1:10)

# ╔═╡ 97b74c99-d460-485a-a297-85735d4cb68a
@bind x7 Slider(-10:0.1:10)

# ╔═╡ 4ba5a65b-c239-4c57-80df-305148c3d832
@bind x8 Slider(-10:0.1:10)

# ╔═╡ 7d78fea3-e4d9-4948-9fa0-7e2a156bc2ec
@bind x9 Slider(-10:0.1:10)

# ╔═╡ 75bdc477-2d3b-43a2-bfb7-899b51aa6d25
@bind x10 Slider(-10:0.1:10)

# ╔═╡ 817599f0-d59f-40e4-8509-b21abcaf9d02

# ╔═╡ 93adbc21-47c7-4050-bf74-b9d782b42930
PlutoUI.RangeSlider

# ╔═╡ d2efde36-bfe2-49e6-9fde-cdd4651c1c34
#m = Metalog([x1, x2, x3, x4, x5, x6, x7, x8, x9, x10])

# ╔═╡ 3e674ded-4ec1-4019-a678-1ecafaa75330
md"Now let's use the random conditional distribution"

# ╔═╡ 1c7e60a0-eba8-40a6-a1ca-5f85025ec451
function nmodes(xs::Vector{Float64})
	@show rcopy(R"multimode::nmodes($xs, 3.0)")
end

# ╔═╡ 19729f65-9185-407e-ae33-f5678cc3157e
nmodes_dist(ml; nsamples = 10000) = nmodes(randsample(ml, nsamples))

# ╔═╡ 4869a4b0-c5a1-453d-8697-b5fb3a394a60
α = manynth(Normal(0, 3), 1:10)

# ╔═╡ ad4473fd-9ddc-49b9-bfc9-79139fe825de
function quantile(m::Metalog, y)
	ncoeff = length(m.α)

	sum(αi * g(i, y) for (i, αi) in enumerate(m.α))
end

# ╔═╡ 8ca63ca5-9bb3-4dd0-a44c-a468fd922837
m = Metalog(rand(3))

# ╔═╡ 0bb1ac69-53a2-4bc4-86a1-fbe04c61bf00
sample() = quantile(m, rand())

# ╔═╡ d2c3ed3e-c337-43fb-9645-ffd605116317
data = [sample() for i = 1:10000]

# ╔═╡ e80fd58a-e123-4303-af20-68d1582d102b
histogram(data)

# ╔═╡ 29c67690-ae3f-42d8-8e97-9ddd9de1474c
nmodes(data)

# ╔═╡ 6ffb21fa-0c48-4ade-a56f-331acca729f2
M = pw(Metalog, α)

# ╔═╡ 8e2da9d5-3e9b-4efc-a190-698235e87ef5

# ╔═╡ 9424a06b-fd04-40e2-a6a2-cbfb6f665e4c
X(ω) = quantile(M(ω), (@~ StdUniform{Float64}())(ω))

# ╔═╡ 64de684f-e38e-48e3-9eb3-88cc2f67b39e
x⃗_α = rid(X, α)

# ╔═╡ 2f9c9660-cbb8-4bf4-8cc0-08a18cfced9e
randsample(x⃗_α)

# ╔═╡ 9ebfa368-5053-49a3-a5cf-f3a3e55c96ab
x⃗_α_nmodes = pw(nmodes_dist, x⃗_α)

# ╔═╡ 583b1289-4a88-4885-b6ab-3aa325369dda
ω = defω()

# ╔═╡ a35a3b27-7203-4b2d-9076-2e777283034d
x⃗_α_nmodes(ω)

# ╔═╡ 954f6d88-c734-4038-86b5-923cddd3b798
datadata = randsample(x⃗_α(ω), 10000)

# ╔═╡ c069e0eb-e8c5-47e1-a362-d7a8764fde41
histogram(datadata)

# ╔═╡ bcc32ecd-7479-4be3-8e3c-a7a913a5cbfb
nmodes(datadata)

# ╔═╡ 3226c9ad-9959-4ede-97ed-924d506fa2e5
nmodes([X(ω) for i = 1:10000])

# ╔═╡ a711f8d2-b87b-47ca-991a-7af29980c14e
M_ = randsample(cnd(M, x⃗_α_nmodes ==ₚ 3), 1; alg = RejectionSample)

# ╔═╡ 3cf857e8-fd28-482a-a664-17a16c5fc958
sim(m::Metalog; nsamples = 10000) = [quantile(m, rand()) for i = 1:nsamples]

# ╔═╡ cd8c89b2-0ad4-474b-9fc8-9b33aae02704
dx = sim(M_[1])

# ╔═╡ 7c321d1f-61c9-4d5c-bca7-4887a96742e7
nmodes(dx)

# ╔═╡ 1c391c6e-9519-4e01-bb73-7e57e8ceb7ea

# ╔═╡ c79791f2-ce51-4308-98b6-58f3b08ae6f9
struct Metalog{A}
	α::A  # Coefficients
end

# ╔═╡ 28fe7962-d5b2-4ec9-a08e-a2d32a7d789c
(m::Metalog)(ω) = quantile(m(ω), (@~ StdUniform{Float64}())(ω))

# ╔═╡ Cell order:
# ╠═01271358-45ae-11ec-30f1-8521b73aa741
# ╠═dc0623d0-3438-4bf3-9bd0-ea2b35589b7e
# ╠═6353b3f7-e6e5-46e4-80e1-272d8b9286bd
# ╠═c79791f2-ce51-4308-98b6-58f3b08ae6f9
# ╠═836d9f48-5f7a-4963-955d-9f4e4a7c8833
# ╠═ad4473fd-9ddc-49b9-bfc9-79139fe825de
# ╠═d5bac81e-e0e1-40f0-b889-529f8ed7719c
# ╠═e7b83751-0720-4cda-8efa-3f14ee0c98e3
# ╠═b565412f-3e41-44f3-879e-fd7aace7a1b6
# ╠═2a0e9e4f-cdf7-4ac7-9dae-9a3249c61dd7
# ╠═af3ebb43-a593-4a07-8b63-869c139f6f05
# ╠═f21c55ce-4539-437a-a667-2f489134a4f9
# ╠═a1f683a2-f9d3-4829-abb0-6fe6b831b1e8
# ╠═97b74c99-d460-485a-a297-85735d4cb68a
# ╠═4ba5a65b-c239-4c57-80df-305148c3d832
# ╠═7d78fea3-e4d9-4948-9fa0-7e2a156bc2ec
# ╠═75bdc477-2d3b-43a2-bfb7-899b51aa6d25
# ╠═817599f0-d59f-40e4-8509-b21abcaf9d02
# ╠═93adbc21-47c7-4050-bf74-b9d782b42930
# ╠═d2efde36-bfe2-49e6-9fde-cdd4651c1c34
# ╠═8ca63ca5-9bb3-4dd0-a44c-a468fd922837
# ╠═0bb1ac69-53a2-4bc4-86a1-fbe04c61bf00
# ╠═d2c3ed3e-c337-43fb-9645-ffd605116317
# ╠═e80fd58a-e123-4303-af20-68d1582d102b
# ╠═29c67690-ae3f-42d8-8e97-9ddd9de1474c
# ╠═3e674ded-4ec1-4019-a678-1ecafaa75330
# ╠═1c7e60a0-eba8-40a6-a1ca-5f85025ec451
# ╠═19729f65-9185-407e-ae33-f5678cc3157e
# ╠═4869a4b0-c5a1-453d-8697-b5fb3a394a60
# ╠═6ffb21fa-0c48-4ade-a56f-331acca729f2
# ╠═28fe7962-d5b2-4ec9-a08e-a2d32a7d789c
# ╠═8e2da9d5-3e9b-4efc-a190-698235e87ef5
# ╠═9424a06b-fd04-40e2-a6a2-cbfb6f665e4c
# ╠═64de684f-e38e-48e3-9eb3-88cc2f67b39e
# ╠═2f9c9660-cbb8-4bf4-8cc0-08a18cfced9e
# ╠═9ebfa368-5053-49a3-a5cf-f3a3e55c96ab
# ╠═583b1289-4a88-4885-b6ab-3aa325369dda
# ╠═a35a3b27-7203-4b2d-9076-2e777283034d
# ╠═954f6d88-c734-4038-86b5-923cddd3b798
# ╠═c069e0eb-e8c5-47e1-a362-d7a8764fde41
# ╠═bcc32ecd-7479-4be3-8e3c-a7a913a5cbfb
# ╠═3226c9ad-9959-4ede-97ed-924d506fa2e5
# ╠═a711f8d2-b87b-47ca-991a-7af29980c14e
# ╠═3cf857e8-fd28-482a-a664-17a16c5fc958
# ╠═cd8c89b2-0ad4-474b-9fc8-9b33aae02704
# ╠═7c321d1f-61c9-4d5c-bca7-4887a96742e7
# ╠═1c391c6e-9519-4e01-bb73-7e57e8ceb7ea

max_stars_repo_name: zenna/expect

max_issues_repo_name: zenna/expect

max_forks_repo_name: zenna/expect


text:
abstract type  AbstractTaxon end

struct Taxon <: AbstractTaxon
    taxid::Int
    name::String
    rank::Symbol
    db::DB
end

Base.show(io::IO, taxon::Taxon) = print(io, "$(taxon.taxid) [$(String(taxon.rank))] $(taxon.name)")
AbstractTrees.printnode(io::IO, taxon::Taxon) = print(io, taxon)
AbstractTrees.nodetype(::Taxon) = Taxon

function Taxon(taxid::Int, db::DB)
    name = db.names[taxid]
    rank = get(db.ranks, taxid, Symbol("no rank"))
    return Taxon(taxid, name, rank, db)
end

function Taxon(name::String, db::DB)
    taxid_canditates = findall(isequal(name), db.names)
    length(taxid_canditates) == 0 && error("There is no candidates for ",name)
    length(taxid_canditates) == 1 && return Taxon(taxid_canditates[1],db)
    length(taxid_canditates) > 1 && error("There are several candidates for ",name)
end

"""
    get(db::DB, taxid::Int, default)

Return the `Taxon` object stored for the given taxid, or the given default value if no mapping for the taxid is present.
"""

function Base.get(db::DB, taxid::Int, default)
    try
        return Taxon(taxid, db)
    catch
        return default
    end
end

"""
    get(db::DB, name::String, default)

Return the `Taxon` object stored for the given name, or the given default value if no mapping for the name is present.
"""

function Base.get(db::DB, name::String, default)
    try
        return Taxon(name, db)
    catch
        return default
    end
end

"""
    taxid(taxon::Taxon)

Return the taxid of the given `Taxon` object.
"""

taxid(taxon::Taxon) = taxon.taxid

"""
    parent(taxon::Taxon)

Return the `Taxon` object that is the parent of the given `Taxon` object.
"""

function Base.parent(taxon::Taxon)
   parent_taxid = get(taxon.db.parents, taxon.taxid, nothing)
   if parent_taxid === nothing
        return nothing
   end
   parent = Taxon(parent_taxid, taxon.db)
   return parent
end

"""
    children(taxon::Taxon)

Return the vector of `Taxon` objects that are children of the given `Taxon` object.
"""

function AbstractTrees.children(taxon::Taxon)
    children_taxid = findall(isequal(taxon.taxid), taxon.db.parents)
    children_taxon = map(x -> Taxon(x, taxon.db), children_taxid)
    return children_taxon
end

struct UnclassifiedTaxon <:AbstractTaxon
    name::String
    rank::Symbol
    source::Taxon
end

function UnclassifiedTaxon(rank, source)
    name = "unclassified " * source.name * " " * String(rank)
    UnclassifiedTaxon(name, rank, source)
end

Base.show(io::IO, taxon::UnclassifiedTaxon) = print(io, "Unclassified [$(String(taxon.rank))] $(taxon.name)")

"""
    rank(taxon::AbstractTaxon)

Return the rank of the given `Taxon` object.
It also works for an `UnclassifiedTaxon` object.
"""

function rank(taxon::AbstractTaxon)
    taxon.rank
end

function name(taxon::AbstractTaxon)
    taxon.name
end


text:
function get_additional_state_dimension(obj::DSSTJacobiansMapper)
    return jcall(obj, "getAdditionalStateDimension", jint, ())
end

function get_name(obj::AbstractJacobiansMapper)
    return jcall(obj, "getName", JString, ())
end

function get_parameters(obj::AbstractJacobiansMapper)
    return jcall(obj, "getParameters", jint, ())
end

function get_parameters_jacobian(obj::DSSTJacobiansMapper, arg0::SpacecraftState, arg1::Vector{Vector{jdouble}})
    return jcall(obj, "getParametersJacobian", void, (SpacecraftState, Vector{Vector{jdouble}}), arg0, arg1)
end

function get_state_jacobian(obj::DSSTJacobiansMapper, arg0::SpacecraftState, arg1::Vector{Vector{jdouble}})
    return jcall(obj, "getStateJacobian", void, (SpacecraftState, Vector{Vector{jdouble}}), arg0, arg1)
end

function set_initial_jacobians(obj::DSSTJacobiansMapper, arg0::SpacecraftState, arg1::Vector{Vector{jdouble}}, arg2::Vector{Vector{jdouble}}, arg3::Vector{jdouble})
    return jcall(obj, "setInitialJacobians", void, (SpacecraftState, Vector{Vector{jdouble}}, Vector{Vector{jdouble}}, Vector{jdouble}), arg0, arg1, arg2, arg3)
end

function set_short_period_jacobians(obj::DSSTJacobiansMapper, arg0::SpacecraftState)
    return jcall(obj, "setShortPeriodJacobians", void, (SpacecraftState,), arg0)
end


text:
using Geodesy
using Geodesy: ED50
using Base.Test

@test_approx_eq distance(ENU(1, 1, 1), ENU(2, 2, 2)) sqrt(3)
@test_approx_eq distance(ECEF(1, 1, 1), ECEF(3, 3, 3)) sqrt(12)
let lla = LLA(1, 2, 3)
    @test_throws MethodError distance(lla, lla)
end

randLLA() = LLA{WGS84}((rand() - .5) * 360,
                       (rand() - .5) * 180,
                       (rand() - .5) * 18_000)

# Test ECEF/ENU distance equivalence

for _ = 1:1_000
    lla1 = randLLA()
    lla2 = randLLA()

    enu1a = ENU(lla1, lla1)
    enu2a = ENU(lla2, lla1)

    ecef1a = ECEF(lla1)
    ecef2a = ECEF(lla2)

    @test_approx_eq distance(enu1a, enu2a) distance(ecef1a, ecef2a)
end

# Test WGS84 data
# Source: http://geographiclib.sourceforge.net/1.1/geodesic.html

for (lat1, lon1, azi1, lat2, lon2, azi2, s12, a12, m12, area) in [
    # random
    (13.34131310843, 0, 64.893732588321, 27.167030166114807887, 81.147882318022755456, 98.180388726997955923, 8484008.9511816, 76.426044312192519403, 6185338.0042860860151, 23494817394476.983874),
    (29.35936791137, 0, 3.562776653705, 67.931453540325162437, 170.572112180718477742, 171.729063831994540061, 9188093.1894032, 82.602801498601116, 6334987.8142568269251, 119120891536718.654989),
    (17.365583427331, 0, 106.832754755529, -23.719740766843660611, 143.94225684272702158, 86.009322008401023546, 16216991.9060447, 146.134657498478139171, 3545186.5038365272107, -14695430107490.419363),
    (48.225427017676, 0, 119.239152833356, 4.900279381681411244, 49.439329609674356182, 144.230393300911335776, 6700530.0165299, 60.344013877711501923, 5530443.3028408957319, 17653442365603.374915),

    # nearly antipodal
    (38.310693597209, 0, 13.055421665051, -37.99474655739563392, 179.800407280950849166, 167.00176496813340991, 19966877.947874, 179.675966711820268126, 75499.3174931018213, 109009980753799.724835),
    (22.977144517038, 0, 122.603891541977, -23.169021570280095852, 179.206837703954580828, 57.523483269982844513, 19944201.7818916, 179.644110241838873311, 55984.8030212716797, -45965171372686.893463),
    (32.617903598529, 0, 117.001691155741, -32.783714078432253875, 179.160224074954599178, 63.207164700860993621, 19944331.9513927, 179.633995938746164263, 50407.0843857503313, -38005514593031.257346),
    (36.25014491329, 0, 104.397824516396, -36.265811672674164215, 179.452719450823176582, 75.646812747826007655, 19976387.5305775, 179.936961729692998254, 9702.9269601585606, -20311056951502.674737),

    # short distances
    (7.512857811123, 0, 19.179996654393, 7.512970768929170776, .0000393712808054, 19.180001802174479956, 13.2265639, .000119209131652008, 13.2265638999905, 3630790.473767),
    (60.577316649216, 0, 51.477685773428, 60.581278059486801371, .010115459170752622, 51.486496706339778805, 708.7627697, .006372099166879661, 708.7627682463562, 6235117278.20435),
    (25.004926378698, 0, 111.712849557686, 25.004832186058303729, .000259567207005928, 111.712959275561025562, 28.203346, .000254055203171289, 28.2033459999077, 77441164.466799),
    (84.118481951125, 0, 94.447662813765, 84.118029655588002885, .056382918243493123, 94.503748905242322208, 647.3047293, .005815043081559776, 647.3047281960948, 39730883653.342566),

    # one end near pole
    (89.99441357489, 0, 108.409374600533, 67.187701448774233223, 71.578016303409870693, 179.986321804301608569, 2546481.7111054, 22.879378547944501141, 2479770.9556252980024, 50706854674009.113617),
    (89.993683218483, 0, 137.249500453314, -8.627003845730393534, 42.751125116580170258, 179.995648792331395129, 10955444.4092702, 98.593858084909860974, 6306522.9078558324466, 30282404215629.211912),
    (89.992173317056, 0, 13.90689177244, 27.740193102384975039, 166.092112190149224805, 179.997869019139657646, 6933163.6633116, 62.346626656988626695, 5649575.4547274351368, 117662908002758.730672),
    (89.99100400556, 0, 21.35479737404, -33.934802175065847398, 158.647383210433154775, 179.996042526461038959, 13759331.1308878, 123.854141018333117523, 5296783.2872557291945, 112385332108250.260572),

    # ends near opposite poles
    (89.997595195216, 0, 179.619593949157, -89.993757086836343885, .526939485610036088, 179.853466395429310428, 20002965.5687943, 179.991323264018997889, 965.8895184347734, 165680964874.397503),
    (89.995765698917, 0, 153.12232383545, -89.996211465670864213, 176.527666913598374567, 30.349989028255825932, 20003874.7682409, 179.999490741949810376, 56.690601014786, -86974926782513.13143),
    (89.993659799817, 0, 116.603332371716, -89.994016930268323357, 134.74824745276636986, 108.648360077793119385, 20003400.6490436, 179.995231659367964821, 530.80913864348, -5635497067621.528193),
    (89.990644140796, 0, 30.46588676412, -89.989017576528539514, 175.124277494839408787, 154.40978557307445724, 20003725.8286419, 179.998152796509465907, 205.6313961813705, 87804891390950.526562),

    # nearly meridional
    (64.715000600285, 0, 179.996224985326, 8.300822939700970608, .003170819821410212, 179.99836617721385843, 6261642.7021158, 56.367229161264380247, 5301555.0052719458082, 1513999906.274607),
    (9.942453894759, 0, 179.997872956141, -83.118967952518969703, .017659335996671428, 179.982568934704861551, 10332922.954675, 93.005823560078472002, 6365254.4421798373908, -10836007339.812858),
    (35.084724228461, 0, 179.999375888592, -17.59161194552863581, .000517914072294356, 179.999463803115894174, 5829774.6537315, 52.5305200980786251, 5046870.9808586923786, 62063370.750408),
    (74.869415682946, 0, .003912050298, 50.158127543633727267, 179.994997645955753609, 179.998404330082463559, 6133133.7314448, 55.115697853366717881, 5235493.1516311444341, 127512500697476.134139),

    # nearly equatorial
    (.000781908619, 0, 89.996315754589, .00113877572073136, 5.579164251062796394, 89.996409342160800773, 621069.7246462, 5.597933071232296147, 620082.0999033472963, 66003372.470098),
    (.00359314356, 0, 89.996052792958, -.003009733334200574, 171.458850316526555958, 90.004405319663036586, 19086711.9083442, 172.035653279020831801, 880771.6586722240014, 5890685275.699263),
    (.004051952929, 0, 90.00047847819, -.000831645030925742, 94.684815255707904411, 90.00398124980729247, 10540265.4251863, 95.003343523917344558, 6332530.6112875111861, 2470357284.501839),
    (.003167017824, 0, 90.000698533704, .002970767627028171, 11.153642898368961157, 90.001297854876781738, 1241617.8463507, 11.191164740509070739, 1233738.0764828187615, 422675979.680255),

    # running between vertices
    (77.168961116231, 0, 90, -77.168961116231, 179.865654057283834783, 90.000000000000000002, 20002265.499854625753, 180, .0000000000001, .000002),
    (68.172063421589, 0, 90, -68.172063421589, 179.775116649629150722, 90.000000000000000001, 19999263.6453246101296, 180, .0000000000001, .000001),
    (72.350857979384, 0, 90, -72.350857979384, 179.816605713764011176, 90.000000000000000003, 20000827.0495409551884, 179.999999999999999999, .0000000000001, .000002),
    (47.474134062996, 0, 90, -47.474134062996, 179.591518711643492433, 90, 19988532.7238678314003, 180, 0, 0),

    # ending close to vertices
    (24.04314613745, 0, 90.045671956752, -24.043146137492923579, 179.448625072251282743, 89.954328067136969735, 19975879.31117, 179.999999946272898281, .0416508759595, -64468770242.133552),
    (43.407004212773, 0, 90.041907806109, -43.407004212805402022, 179.561037570819915051, 89.958092235644803931, 19986149.3817786, 179.999999955708889468, .0239608044566, -59236521272.663064),
    (60.893175111551, 0, 90.030887879944, -60.893175111602473374, 179.705866045884665835, 89.969112291280908219, 19995946.5267377, 179.99999990435001247, .0152818050346, -43716837097.405323),
    (57.297318482011, 0, 90.006736496503, -57.297318482011214432, 179.673356341455815136, 89.993263506332003021, 19994084.1111944, 179.999999998173649615, .000475461712, -9532089026.025471)
]
    p1 = LL{WGS84}(lon1, lat1)
    p2 = LL{WGS84}(lon2, lat2)

    try
        vi = Geodesy.vicentys_inverse(p1, p2)
        @test abs(vi[1] - s12) < 1e-4
    catch e
        @test contains(e.msg, "antipodal")
        @test 179 < abs(lon2 - lon1) < 181
    end
end

# Test data over the Hayford ellipsoid, "exact" distances this time
# Source: http://geographiclib.sourceforge.net/1.1/geodesic.html

function todec(x)
    r = match(r"^([\-0-9]+)d([0-9]+)'([0-9\.]+)$", x)
    Geodesy.dms2decimal(map(float, r.captures)...)
end
function testvals(t)
    lat1, lat2, dlon, azi1, azi2, s12 = t
    todec(lat1), todec(lat2), todec(dlon), todec(azi1), todec(azi2), float(s12)
end

for t in [
    ("37d19'54.95367" ,"26d07'42.83946"  ,"41d28'35.50729" ,"95d27'59.63088905553491" ,"118d05'58.96160858886728","4085966.7025902201825" )
    ("35d16'11.24862" ,"67d22'14.77638"  ,"137d47'28.31435","15d44'23.74849770324814" ,"144d55'39.92147266777383","8084823.8382961415712" )
    ("1d00'00"        ,"-0d59'53.83076"  ,"179d17'48.02997","88d59'59.99897053696689" ,"91d00'06.11835637659787" ,"19959999.9998034962728")
    ("1d00'00"        ,"1d01'15.18952"   ,"179d46'17.84244","4d59'59.9999565312182"   ,"174d59'59.88480004914149","19780006.5587880182731")
    ("41d41'45.88"    ,"-41d41'46.20"    ,"179d59'59.44"   ,"179d58'49.16247860972594","0d01'10.83761892414861"  ,"20004566.7228054132931")
    ("0d00'00"        ,"0d00'00"         ,"179d41'49.78063","30d00'00.00002084815889" ,"149d59'59.99997915184111","19996147.416826781925" )
    ("30d00'00"       ,"-30d00'00"       ,"179d40'00"      ,"39d24'51.80601183179472" ,"140d35'08.19398816820528","19994364.6068583984182")
    ("60d00'00"       ,"-59d59'00"       ,"179d50'00"      ,"29d11'51.07006518417487" ,"150d49'06.86792886226858","20000433.9629039632049")
    ("30d00'00"       ,"-29d50'00"       ,"179d48'00"      ,"16d02'28.33895348073478" ,"163d59'10.33689436714454","19983420.1535833515076")
    ("30d00'00"       ,"-29d55'00"       ,"179d48'00"      ,"18d38'12.55689701007199" ,"161d22'45.43724069005652","19992241.7634404403113")
    ("34d28'44.76421" ,"-34d28'44.76421" ,"179d30'00"      ,"89d59'59.78608904242135" ,"90d00'00.21391095757865" ,"19981603.2781440234735")
    ("0d00'00"        ,"0d00'00"         ,"179d23'38.18182","89d59'51.5863776476884"  ,"90d00'08.4136223523116"  ,"19970827.8695289752144")
    ("34d28'44.764213","-34d28'44.764213","179d30'00"      ,"90d00'00.00000000479721" ,"90d00'00.00000000479721" ,"19981603.2781440234735")
    ("56d41'58.297496","-56d41'58.297496","179d40'00"      ,"90d00'00.0000006754891"  ,"90d00'00.0000006754891"  ,"19994364.6068583984183")
    ("29d45'00"       ,"29d45'00.23848"  ,"0d00'00.25626"  ,"43d09'29.0634310812721"  ,"43d09'29.19059165112047" ,"10.0665488514951"      )
    ("41d41'45.88"    ,"41d41'46.2"      ,"0d00'00.56"     ,"52d40'39.39067110974434" ,"52d40'39.76317180931766" ,"16.2839750636094"      )
    ("46d00'00"       ,"46d00'01"        ,"0d00'01.816"    ,"51d41'11.54166708300919" ,"51d41'12.8479912183896"  ,"49.8037586266135"      )
    ("40d00'00"       ,"40d00'02"        ,"0d00'04.75"     ,"61d18'00.21784992705284" ,"61d18'03.27110871411274" ,"128.4581417556041"     )
    ("38d00'00"       ,"38d00'04.765"    ,"0d00'05.554"    ,"42d41'10.35602545797828" ,"42d41'13.77545984493403" ,"199.8717376568129"     )
    ("30d00'00"       ,"37d53'32.46584"  ,"116d19'16.68843","45d00'00.00000438165649" ,"129d08'12.32600911217221","10002499.999860115911" )
    ("30d19'54.95367" ,"-30d11'50.15681" ,"179d58'17.84244","2d23'52.10812966032674"  ,"177d36'19.67010864697937","19989590.5480170316779")
    ("0d39'49.12586"  ,"-0d45'14.13112"  ,"179d58'17.84244","177d39'39.01010365081902","2d20'21.15303599834387"  ,"19994529.4454322340309")
    ("0d00'54.95367"  ,"0d00'42.83946"   ,"179d28'17.84244","54d08'27.73161935117439" ,"125d51'32.27232710890937","19977290.7711390609949")
    ("40d00'00"       ,"-40d00'05.75932" ,"179d55'15.59578","170d15'10.8812277902834" ,"9d44'49.94565915545135"  ,"20003827.8511392345191")
    ("37d00'00"       ,"28d15'36.69535"  ,"2d37'39.52918"  ,"164d59'59.99997936397788","166d25'16.2593819478948" ,"1000000.0001515200393" )
    ("38d30'45"       ,"-35d25'35"       ,"179d45'00"      ,"3d22'19.5694402667547"   ,"176d45'41.43721792855263","19661438.0251956080829")
    ("60d00'00"       ,"-60d00'00"       ,"179d41'47"      ,"90d00'00.00891264563132" ,"90d00'00.00891264563132" ,"19996104.3689008382234")
]
    lat1, lat2, dlon, azi1, azi2, s12 = testvals(t)

    p1 = LL{ED50}(0, lat1)
    p2 = LL{ED50}(dlon, lat2)

    try
        vi = Geodesy.vicentys_inverse(p1, p2)
        @test abs(vi[1] - s12) < 1e-4
    catch e
        @test contains(e.msg, "antipodal")
        @test 179 < dlon < 181
    end
end


text:
psttable = readdlm(pkgdir(@__MODULE__)*"/src/transposedpst.txt", Int)

function pst(p::Piece, sq::Square, table=psttable)::Real
    sqval = (pcolor(p) == WHITE ? sq.val : ((sq.val - 1) ⊻ 7) + 1)::Int
    table[ptype(p).val, sqval]
end

pstEval(g::Union{Game,SimpleGame}, d = 1, table=psttable, c::PieceColor=WHITE)::Real =
    isdraw(g) ? 0 : pstEval(board(g), d, table, c)

function pstEval(b::Board, d::Integer = 1, table=psttable, c::PieceColor=WHITE)
    score = 0
    ps = Piece.((1:6) .+ (8*(c==BLACK)))
    for p = ps
        for s = pieces(b, p)
            score += pst(p, s, table)
        end
    end
    score
end

function pstFull(b::Board, d::Integer=1, table=psttable)
    pstEval(b, d, table, sidetomove(b)) - pstEval(b, d, table, -sidetomove(b))
end

function pstFull(g::Union{Game,SimpleGame}, d::Integer=1, table=psttable)
    c = sidetomove(board(g))
    pstEval(g, d, table, c) - pstEval(g, d, table, -c)
end

function pstAndMate(g, d::Integer=1, table=psttable; matescore=10000)
    b = typeof(g) == Board ? g : board(g)
    if ischeckmate(b)
        return -max(d, 1) * matescore
    end
    if isdraw(g) #|| isrepetitiondraw(g) # too slow
        return 0
    end
    any(m->LichessBot.move_is_mate(b,m), moves(b)) && return max(d, 1) * matescore
    return pstFull(g, d, table)
end


text:
module Day16

_inp = "59773419794631560412886746550049210714854107066028081032096591759575145680294995770741204955183395640103527371801225795364363411455113236683168088750631442993123053909358252440339859092431844641600092736006758954422097244486920945182483159023820538645717611051770509314159895220529097322723261391627686997403783043710213655074108451646685558064317469095295303320622883691266307865809481566214524686422834824930414730886697237161697731339757655485312568793531202988525963494119232351266908405705634244498096660057021101738706453735025060225814133166491989584616948876879383198021336484629381888934600383957019607807995278899293254143523702000576897358"

inp = parse.(Int, split(Day16._inp, ""))

pat(i) = IterTools.imap(x -> (x for _ in 1:i), [0, 1, 0, -1]) |> Iterators.flatten

function simulate(v, reps)
    v = v |> copy |> x -> repeat(x, reps)
    nextv = v |> similar
    for j in 1:100
        for i in 1:length(v)
            pattern = i |> pat |> Iterators.cycle |> x -> Iterators.drop(x, 1)
            nextv[i] = abs(sum(prod, zip(v, pattern)) % 10)
        end

        v = nextv
    end

    v
end

solvea(v) = simulate(v, 1)[1:8] |> join

function solveb(v)
    offset = sum(v[k] * 10^(k - 1) for k in 1:7)
    simulate(v, 10000)[offset:offset + 8] |> join
end

end

max_stars_repo_name: AnAverageHuman/competitive

max_issues_repo_name: AnAverageHuman/competitive

max_forks_repo_name: AnAverageHuman/competitive


text:
### Getting and Setting Properties
#Mutable not defined
Mutable = include("MutableTypes.jl")
include("GLib.jl")
include("gtype.jl")
include("../gen/gtk_get_set_gen.jl")
#=
requires GV -that requires-> GValue (gvalues.jl) -that requires->  gtk_get_set_gen.jl ("../gen/gtk_get_set_gen.jl")
=#
struct GValue
    g_type::GType
    field2::UInt64
    field3::UInt64
    GValue() = new(0, 0, 0)
end
const GV = Union{Mutable{GValue},Ptr{GValue}} #Mutable not #try:added MutableTypes.jl
Base.zero(::Type{GValue}) = GValue()
function gvalue(::Type{T}) where {T}
    v = mutable(GValue())
    v[] = T
    v
end
function gvalue(x)
    T = typeof(x)
    v = gvalue(T)
    v[T] = x
    v
end
function gvalues(xs...)
    v = zeros(GValue, length(xs))
    for (i, x) in enumerate(xs)
        T = typeof(x)
        gv = mutable(v, i)
        gv[] = T  # init type
        gv[T] = x # init value
    end
    finalizer((v) -> for i = 1:length(v)
            ccall((:g_value_unset, libgobject), Nothing, (Ptr{GValue},), pointer(v, i))
        end, v)
    v
end

function setindex!(dest::GV, src::GV)
    ccall((:g_value_transform, libgobject), Cint, (Ptr{GValue}, Ptr{GValue}), src, dest) != 0
    src
end # GV not defined

setindex!(::Type{Nothing}, v::GV) = v #GV not defined
setindex!(v::GLib.GV, x) = setindex!(v, x, typeof(x)) #GLib (include ) #GV not defined
setindex!(gv::GV, x, i::Int) = setindex!(mutable(gv, i), x)# GV not defined

getindex(gv::GV, i::Int, ::Type{T}) where {T} = getindex(mutable(gv, i), T)#GV not defined
getindex(gv::GV, i::Int) = getindex(mutable(gv, i)) # GV not defined
getindex(v::GV, i::Int, ::Type{Nothing}) = nothing#GV not defined

let handled = Set()
    global make_gvalue, getindex
    function make_gvalue(pass_x, as_ctype, to_gtype, with_id, cm::Module, allow_reverse::Bool = true, fundamental::Bool = false)
        with_id === :error && return
        if isa(with_id, Tuple)
            with_id = with_id::Tuple{Symbol,Any}
            with_id = :(ccall($(Expr(:tuple, Meta.quot(Symbol(string(with_id[1], "_get_type"))), with_id[2])), GType, ()))
        end
        if pass_x !== Union{} && !(pass_x in handled)
            Core.eval(cm, quote
                function Base.setindex!(v::GLib.GV, ::Type{T}) where {T<:$pass_x}
                    ccall((:g_value_init, GLib.libgobject), Nothing, (Ptr{GLib.GValue}, Csize_t), v, $with_id)
                    v
                end
                function Base.setindex!(v::GLib.GV, x, ::Type{T}) where {T<:$pass_x}
                    $(
                        if to_gtype == :string
                            :(x = GLib.bytestring(x))
                        elseif to_gtype == :pointer || to_gtype == :boxed
                            :(x = GLib.mutable(x))
                        elseif to_gtype == :gtype
                            :(x = GLib.g_type(x))
                        end
                    )
                    ccall(($(string("g_value_set_", to_gtype)), GLib.libgobject), Nothing, (Ptr{GLib.GValue}, $as_ctype), v, x)
                    if isa(v, GLib.MutableTypes.MutableX)
                        finalizer((v::GLib.MutableTypes.MutableX) -> ccall((:g_value_unset, GLib.libgobject), Nothing, (Ptr{GLib.GValue},), v), v)
                    end
                    v
                end
            end)
        end
        if to_gtype == :static_string
            to_gtype = :string
        end
        if pass_x !== Union{} && !(pass_x in handled)
            push!(handled, pass_x)
            Core.eval(cm, quote
                function Base.getindex(v::GLib.GV, ::Type{T}) where {T<:$pass_x}
                    x = ccall(($(string("g_value_get_", to_gtype)), GLib.libgobject), $as_ctype, (Ptr{GLib.GValue},), v)
                    $(
                        if to_gtype == :string
                            :(x = GLib.bytestring(x))
                        elseif pass_x == Symbol
                            :(x = Symbol(x))
                        end
                    )
                    return Base.convert(T, x)
                end
            end)
        end
        if fundamental || allow_reverse
            fn = Core.eval(cm, quote
                function (v::GLib.GV)
                    x = ccall(($(string("g_value_get_", to_gtype)), GLib.libgobject), $as_ctype, (Ptr{GLib.GValue},), v)
                    $(
                        if to_gtype == :string
                            :(x = GLib.bytestring(x))
                        end
                    )
                    $(
                        if pass_x !== Union{}
                            :(return Base.convert($pass_x, x))
                        else
                            :(return x)
                        end
                    )
                end
            end)
            allow_reverse && pushfirst!(gvalue_types, [pass_x, Core.eval(cm, :(() -> $with_id)), fn])
            return fn
        end
        return nothing
    end
end #let

macro make_gvalue(pass_x, as_ctype, to_gtype, with_id, opt...)
    esc(:(make_gvalue($pass_x, $as_ctype, $to_gtype, $with_id, $__module__, $(opt...))))
end

function make_gvalue_from_fundamental_type(i, cm)
    (name, ctype, juliatype, g_value_fn) = fundamental_types[i]
    return make_gvalue(juliatype, ctype, g_value_fn, fundamental_ids[i], cm, false, true)
end

const gvalue_types = Any[]
const fundamental_fns = tuple(Function[make_gvalue_from_fundamental_type(i, @__MODULE__) for
                                       i in 1:length(fundamental_types)]...)#ERROR:fundamental_types not defined
@make_gvalue(Symbol, Ptr{UInt8}, :static_string, :(g_type(AbstractString)), false) #GV
@make_gvalue(Type, GType, :gtype, (:g_gtype, :libgobject))#GType
@make_gvalue(Ptr{GBoxed}, Ptr{GBoxed}, :gboxed, :(g_type(GBoxed)), false)
#=WARNING: redefinition of constant gvalue_types. This may fail, cause incorrect answers, or produce other errors.
Any[]=#

function getindex(gv::GV, ::Type{Any})
    gtyp = unsafe_load(gv).g_type
    if gtyp == 0
        error("Invalid GValue type")
    end
    if gtyp == g_type(Nothing)
        return nothing
    end
    # first pass: fast loop for fundamental types
    for (i, id) in enumerate(fundamental_ids)
        if id == gtyp  # if g_type == id
            return fundamental_fns[i](gv)
        end
    end
    # second pass: user defined (sub)types
    for (typ, typefn, getfn) in gvalue_types
        if g_isa(gtyp, typefn())
            return getfn(gv)
        end
    end
    # last pass: check for derived fundamental types (which have not been overridden by the user)
    for (i, id) in enumerate(fundamental_ids)
        if g_isa(gtyp, id)
            return fundamental_fns[i](gv)
        end
    end
    typename = g_type_name(gtyp)
    error("Could not convert GValue of type $typename to Julia type")
end
#end

get_gtk_property(w::GObject, name::AbstractString, ::Type{T}) where {T} = get_gtk_property(w, String(name)::String, T)
get_gtk_property(w::GObject, name::Symbol, ::Type{T}) where {T} = get_gtk_property(w, String(name), T)
function get_gtk_property(w::GObject, name::String, ::Type{T}) where {T}
    v = gvalue(T)
    ccall((:g_object_get_property, libgobject), Nothing,
        (Ptr{GObject}, Ptr{UInt8}, Ptr{GValue}), w, name, v)
    val = v[T]
    ccall((:g_value_unset, libgobject), Nothing, (Ptr{GValue},), v)
    return val
end

set_gtk_property!(w::GObject, name, ::Type{T}, value) where {T} = set_gtk_property!(w, name, convert(T, value))
set_gtk_property!(w::GObject, name::AbstractString, value) = set_gtk_property!(w::GObject, String(name)::String, value)
set_gtk_property!(w::GObject, name::Symbol, value) = set_gtk_property!(w::GObject, String(name), value)
function set_gtk_property!(w::GObject, name::String, value)
    ccall((:g_object_set_property, libgobject), Nothing,
        (Ptr{GObject}, Ptr{UInt8}, Ptr{GValue}), w, name, gvalue(value))

    w
end

struct FieldRef{T}
    obj::T
    field::Symbol

    global function getproperty(obj::T, field::Symbol) where {T<:GObject}
        isdefined(obj, field) && return getfield(obj, field)
        new{T}(obj, field)
    end

    FieldRef(obj::T, field::Symbol) where {T} = new{T}(obj, field)
end

#= UncommentMe
getindex(f::FieldRef, ::Type{T}) where {T} = get_gtk_property(f.obj, f.field,T)

function setindex!(f::FieldRef, value::T, ::Type{T}) where {T}
    isdefined(f.obj,f.field) && return setfield!(f.obj, f.field, value)
    set_gtk_property!(f.obj, f.field, value)
    return f
end #FieldRef  undefined
=#
setindex!(f::FieldRef, value::K, ::Type{T}) where {K,T} = setindex!(f, convert(T, value), T)

function show(io::IO, w::GObject)
    READABLE = 0x00000001
    DEPRECATED = 0x80000000
    print(io, typeof(w), '(')
    if unsafe_convert(Ptr{GObject}, w) == C_NULL
        print(io, "<NULL>)")
        return
    end
    n = mutable(Cuint)
    props = ccall((:g_object_class_list_properties, libgobject), Ptr{Ptr{GParamSpec}},
        (Ptr{Nothing}, Ptr{Cuint}), G_OBJECT_GET_CLASS(w), n)
    v = gvalue(String)
    first = true
    for i = 1:unsafe_load(n)
        param = unsafe_load(unsafe_load(props, i))
        if !first
            print(io, ", ")
        else
            first = false
        end
        print(io, GLib.bytestring(param.name))
        if (param.flags & READABLE) != 0 &&
           (param.flags & DEPRECATED) == 0 &&
           (ccall((:g_value_type_transformable, libgobject), Cint,
               (Int, Int), param.value_type, g_type(AbstractString)) != 0)
            ccall((:g_object_get_property, libgobject), Nothing,
                (Ptr{GObject}, Ptr{UInt8}, Ptr{GValue}), w, param.name, v)
            str = ccall((:g_value_get_string, libgobject), Ptr{UInt8}, (Ptr{GValue},), v)
            value = (str == C_NULL ? "NULL" : GLib.bytestring(str))
            if param.value_type == g_type(AbstractString) && str != C_NULL
                print(io, "=\"", value, '"')
            else
                print(io, '=', value)
            end
        end
    end
    print(io, ')')
    ccall((:g_value_unset, libgobject), Ptr{Nothing}, (Ptr{GValue},), v)
end

#immutable GTypeQuery
#  g_type::Int
#  type_name::Ptr{UInt8}
#  class_size::Cuint
#  instance_size::Cuint
#  GTypeQuery() = new(0, 0, 0, 0)
#end
#function gsizeof(gtyp)
#    q = mutable(GTypeQuery)
#    ccall((:g_type_query, libgobject), Nothing, (Int, Ptr{GTypeQuery},), gtyp, q)
#    q[].instance_size
#end

max_stars_repo_name: adamwillisMastery/MiniGtk

max_issues_repo_name: adamwillisMastery/MiniGtk

max_forks_repo_name: adamwillisMastery/MiniGtk


text:
const ev_to_K = 11604.52500617u"K"
const e_statC = 4.80320425e-10u"cm^(3/2) * g^(1/2) * s^−1"

"""
    electron_gyrofrequency(B)

The electron gyrofrequency, or the angular frequency of the circular motion of an
electron in the plane perpendicular to the magnetic field.

# Arguments
- `B`: magnetic flux density, typically in units of tesla
"""
function electron_gyrofrequency(B)
    return (e_statC * B) / (ElectronMass * SpeedOfLightInVacuum)
end

"""
    ion_gyrofrequency(B)

The ion gyrofrequency, or the angular frequency of the circular motion of an ion in
the plane perpendicular to the magnetic field.

# Arguments
 - `B`: magnetic flux density, typically in units of tesla
 - `Z`: charge state, unitless
 - `m_i`: ion mass, typically in grams
"""
function ion_gyrofrequency(B, Z, m_i)
    return (Z * e_statC * B) / (m_i * SpeedOfLightInVacuum)
end

"""
    electron_plasma_frequency(n_e)

Electron plasma frequency, or the frequency at which the electrons in the plasma oscillate

# Arguments
 - `n_e`: electron number density, typically in units of 1/cm^3
"""
function electron_plasma_frequency(n_e)
    return sqrt((4π * n_e * e_statC^2) / (ElectronMass |> u"g"))
end

"""
    ion_plasma_frequency(n_i, m_i, Z)

Ion plasma frequency, or the frequency at which the ions in the plasma oscillate

# Arguments
 - `n_i`: ion number density, typically in units of 1/cm^3
 - `Z`: charge state, dimensionless
 - `m_i`: ion mass, e.g., grams
"""
function ion_plasma_frequency(n_i, m_i, Z)
    return sqrt((4π * n_i * Z^2 * e_statC^2) / m_i)
end

"""
# Arguments
 - ``:
"""
function electron_trapping_rate()
end

"""
# Arguments
 - ``:
"""
function ion_trapping_rate()
end

"""
# Arguments
 - ``:
"""
function electron_collision_rate()
end

"""
# Arguments
 - ``:
"""
function ion_collision_rate()
end

"""
# Arguments
 - ``:
"""
function electron_debroglie_length()
end

"""
# Arguments
 - ``:
"""
function classical_minimum_approach_dist()
end

"""
# Arguments
 - ``:
"""
function electron_gryoradius()
end

"""
# Arguments
 - ``:
"""
function ion_gyroradius()
end

"""
# Arguments
 - ``:
"""
function Debye_length()
end

"""
# Arguments
 - ``:
"""
function magnetic_debye_length()
end

"""
    electron_thermal_velocity(T_e)

# Arguments
 - `T_e`: electron temperature in eV
"""
function electron_thermal_velocity(T_e::Unitful.Energy)
    T = ustrip(u"eV", T_e) * ev_to_K
	electron_thermal_velocity(T)
end

function electron_thermal_velocity(T_e::Unitful.Temperature)
	sqrt(T_e * BoltzmannConstant / ElectronMass) |> u"cm/s"
end

"""
# Arguments
 - ``:
"""
function ion_thermal_velocity(T_i, m_i)
    return sqrt((BoltzmannConstant * T_i) / m_i)
end

"""
    ion_sound_velocity(γ, Z, T_e, m_i)

Ion sound speed, the speed of the longitudinal waves resulting from the
mass of the ions and the pressure of the electrons.

# Arguments
 - `γ`: adiabatic index
 - `Z`: charge state
 - `T_e`: electron temperature in eV
 - `m_i`: ion mass
"""
function ion_sound_velocity(γ, Z, T_e, m_i)
    return sqrt(γ * Z * k_B * T_e / m_i)
end

"""
    alfven_velocity(B, n_i, m_i)

The Alfven velocity, the speed of the waves resulting from the mass of the
ions and the restoring force of the magnetic field

# Arguments
 - `B`: magnetic flux, typically in teslas
 - `n_i`: ion number density, typically in 1/cm^3
 - `m_i`: ion mass, typically in grams
"""
function alfven_velocity(B, n_i, m_i)
    return B / sqrt(4π * n_i * m_i)
end