text:
using AmbientForcing
using PowerDynamics, Distributions
using OrdinaryDiffEq
using Test

@testset "PD_example" begin
    # Lets create a small test power grid with three nodes
    line_list = []
    append!(line_list, [StaticLine(from=1, to=2, Y=-1im / 0.02 + 4)])
    append!(line_list, [StaticLine(from=1, to=3, Y=-1im / 0.02 + 4)])

    # The PQAlgebraic Node is our constraint. The Power Output of node 3 is fixed.
    node_list = []
    append!(node_list, [SlackAlgebraic(U=1, Y_n=0)])
    append!(node_list, [FourthOrderEq(H=3.318, P=-0.6337, D=0.1, Ω=50, E_f=0.5, T_d_dash=0.1, T_q_dash=8.690, X_d_dash=0.111, X_q_dash=0.103, X_d=0.1, X_q=0.6)])
    append!(node_list, [PQAlgebraic(P=-0.6337, Q=0.0)])

    # Using the right hand side as our ODEFunction
    pg = PowerGrid(node_list, line_list)
    rpg = rhs(pg)

    # Accessing the constraint equation g of the power grid
    g = constraint_equations(rpg)

    # The operation point is a fixed point and naturally lies on the manifold
    # Thus it can be used as the initial condition for our differential equation
    op = find_operationpoint(pg)

    # Lets check if g(op) ≈ 0 meaning that the constraint is fulfilled
    @test isapprox(sum(g(op.vec)), 0.0, atol=1e-8)

    # Lets generate a random vector from the ambient space
    # First we want to perturb all variables in the grid
    Frand = random_force(rpg, [0.0, 1], Uniform)
    afoprob = ambient_forcing_problem(rpg, op.vec, 2.0, Frand, method=:ForwardDiff)
    z_new_all = ambient_forcing(afoprob, op.vec, 2.0, Frand) # Our new valid initial condition
    @test isapprox(sum(g(z_new_all)), 0.0, atol=1e-8)

    ## Next: let's perturb just the voltage at node 3!
    # Getting the index of the real and imaginary part of the voltage

    idx = idx_exclusive(rpg, ["u_r_3", "u_i_3"])
    # generate a vector only with non-vanishing components the voltage at node 3
    Frand = random_force(rpg, [0.0, 1.0], Uniform, idx)
    z_new_node_3 = ambient_forcing(afoprob, Frand)
    @test isapprox(sum(g(z_new_node_3)), 0.0, atol=1e-8)

    # It is also possible to perturb the variable using different distributions:
    # Typically SNBS the angle θ and ω are perturbed differently
    idx = idx_exclusive(rpg, ["u_r_2", "u_i_2", "θ_2", "ω_2"])
    τ = 2.0 # the integration time

    # When you want to sample eg. the angle θ from a box of [0, 2π]
    # You have to make sure to dived the distribution argument by the integration time τ
    dist_vec = [[0, 1], [0, 2], [0, 2π] ./ τ, [-5, 5] ./ τ]

    Frand = random_force(rpg, dist_vec, Uniform, idx)
    z_new_node_2 = ambient_forcing(afoprob, op.vec, τ, Frand)
    @test isapprox(sum(g(z_new_node_2)), 0.0, atol=1e-8)
end

@testset "OrdinaryDiffEq_example" begin
    # Adopting the Robertson Example from the DifferentialEquation Docs
    # https://diffeq.sciml.ai/stable/tutorials/dae_example/
    function rober(du, u, p, t)
        y₁, y₂, y₃ = u
        du[1] = -0.04 * y₁ + 1e4 * y₂ * y₃
        du[2] = 0.04 * y₁ - 1e4 * y₂ * y₃ - 3e7 * y₂^2
        du[3] = y₁ + y₂ + y₃ - 1
        nothing
    end

    # Creating the mass matrix M
    # The last row depicts our constraint
    M = [1.0 0 0
        0 1.0 0
        0 0 0]

    # Setting up the DAE as an ODE in mass matrix form
    ode_rober = ODEFunction(rober, mass_matrix=M)

    # Choosing an inital condition which fulfills the constraint
    u0 = [1.0, 0.0, 0.0]

    g_rober = constraint_equations(ode_rober)
    @test isapprox(sum(g_rober(u0)), 0.0, atol=1e-8)

    # Randomly perturbing all variables
    Frand = random_force(ode_rober, [0.0, 1], Uniform)
    afoprob_rober = ambient_forcing_problem(ode_rober, u0, 2.0, Frand)
    z_new_all = ambient_forcing(afoprob_rober, Frand)
    @test isapprox(sum(g_rober(z_new_all)), 0.0, atol=1e-8)

    # only perturbing the second variable y₂
    h = [0, 1, 0]
    z_new_2 = ambient_forcing(afoprob_rober, h)
    @test isapprox(sum(g_rober(z_new_2)), 0.0, atol=1e-8)
end


text:
# System module.
module System

using RCall
using Conda

export LIB_SKL_AVAILABLE,
       LIB_CRT_AVAILABLE

import PyCall: pyimport_conda, pycall

function check_py_dep()
  is_available = true
  try
    pyimport_conda("sklearn", "scikit-learn")
  catch
    @info "scikitlearn not available"
    is_available = false
  end
  return is_available
end

function check_r_dep()
  is_available = true
  try
    R"library(caret)"
  catch
    @info "caret not available"
    is_available = false
  end
  return is_available
end

## Check system for python dependencies.
LIB_SKL_AVAILABLE = check_py_dep()
LIB_CRT_AVAILABLE = check_r_dep()

end # module

ext: jl

max_stars_repo_name: ppalmes/TSMLextra

max_issues_repo_name: ppalmes/TSMLextra

max_forks_repo_name: ppalmes/TSMLextra


text:
function obojti(r,y)
    tmp=0
    while isborder(r,Nord)
        move!(r,West)
        tmp+=1
    end
    if y>1
        move!(r,Nord)
        move!(r,Nord)
        y-=2
    end
    for i in 1:tmp
        move!(r,Ost)
    end
end
function todefault(r)
    x=0
    y=0
    while isborder(r,Sud)==false || isborder(r,West)==false
        if isborder(r,Sud)==false
            y+=1
            move!(r,Sud)
        end
        if isborder(r,West)==false
            x+=1
            move!(r,West)
        end
    end
    return x,y
end
function tostart(r,y)
    while y>0
        if isborder(r,Nord)==false
            while isborder(r,Nord)==false && y>0
                move!(r,Nord)
                y-=1
            end
        end
        if isborder(r,Nord)==true && y>1
            obojti(r,y)
            y-=2
        end
    end
end

function calc_toboarder(r,side)
    k=0
    while isborder(r,side)==false
        move!(r,side)
        k+=1
    end
    return k
end
function task11(r)
    x,y=todefault(r)
    for i in 1:y
        move!(r,Nord)
    end
    putmarker!(r)
    k=calc_toboarder(r,Nord)
    for i in 1:x
        move!(r,Ost)
    end
    putmarker!(r)
    l=calc_toboarder(r,Ost)
    for i in 1:k
        move!(r,Sud)
    end
    putmarker!(r)
    calc_toboarder(r,Sud)
    for i in 1:l
        move!(r,West)
    end
    putmarker!(r)
    tostart(r,y)
end


text:
# Local indices
"""
A `FaceIndex` wraps an (Int, Int) and defines a face by pointing to a (cell, face).
"""
struct FaceIndex
    cellidx::Int
    idx::Int
end

"""
A `FacetIndex` wraps an (Int, Int) and defines a facet by pointing to a (cell, facet).
"""
struct FacetIndex
    cellidx::Int
    idx::Int
end

"""
A `EdgeIndex` wraps an (Int, Int) and defines an edge by pointing to a (cell, edge).
"""
struct EdgeIndex
    cellidx::Int
    idx::Int
end

function Base.getindex(e::EdgeIndex, n::Int64)
  if n == 1
    return e.cellidx
  elseif n ==2
    return e.idx
  else
    throw(BoundsError(e,n))
  end
end

function get_Normal(mesh, edge_idx::EdgeIndex)
  coords = getverticescoords(mesh, edge_idx)
  v1 = coords[2] - coords[1]
  n1 = Tensors.Vec{2}((v1[2], -v1[1]))
  return n1/norm(n1)
end

function get_vertices_matrix(mesh)
  vertices_m = Matrix{eltype(mesh.vertices[1])}(undef,getnvertices(mesh),getdim(mesh))
  for k in 1:getnvertices(mesh)
      vertices_m[k,:] = getvertexcoords(mesh,k)
  end
  vertices_m
end

function cell_volume(mesh, cell_idx::Int)
  d = getdim(mesh)
  if  d == 2
    N = getncellvertices(mesh, cell_idx)
    verts = getverticescoords(mesh, cell_idx)
    return 0.5*abs(sum(verts[j][1]*verts[mod1(j+1,N)][2]-verts[mod1(j+1,N)][1]*verts[j][2] for j ∈ 1:N))
  else
    error("cell volume not available for cells of dim $d")
  end
end

function cell_centroid(mesh, cell_idx::Int)
  d = getdim(mesh)
  if  d == 2
    verts = getverticescoords(mesh, cell_idx)
    vertices = [StaticArrays.SVector(x[1],x[2]) for x in verts]
    chull = PlanarConvexHulls.ConvexHull{PlanarConvexHulls.CCW}(vertices)
    return Tensors.Vec{2}(PlanarConvexHulls.centroid(chull))
  else
    error("cell centroid not available for cells of dim $d")
  end
end

function cell_diameter(mesh, cell_idx::Int)
  d = getdim(mesh)
  if  d == 2
    verts = getverticescoords(mesh, cell_idx)
    vertices = [StaticArrays.SVector(x[1],x[2]) for x in verts]
    chull = PlanarConvexHulls.ConvexHull{PlanarConvexHulls.CCW}(vertices)
    h = 0.0
    for vert in chull.vertices
        h = max(h,maximum([norm(vert-vert2) for vert2 in chull.vertices]))
    end
    h
  else
    error("cell diameter not available for cells of dim $d")
  end
end

max_stars_repo_name: Paulms/jFEMToools

max_issues_repo_name: Paulms/jFEMToools

max_forks_repo_name: Paulms/jFEMToools


text:
const SimpleGraphEdge = SimpleEdge

"""
    SimpleGraph{T}

A type representing an undirected graph.
"""
mutable struct SimpleGraph{T <: Integer} <: AbstractSimpleGraph{T}
    ne::Int
    fadjlist::Vector{Vector{T}} # [src]: (dst, dst, dst)
end

eltype(x::SimpleGraph{T}) where T = T

# Graph{UInt8}(6), Graph{Int16}(7), Graph{UInt8}()
"""
    SimpleGraph{T}(n=0)

Construct a `SimpleGraph{T}` with `n` vertices and 0 edges.
If not specified, the element type `T` is the type of `n`.

## Examples
```jldoctest
julia> SimpleGraph(UInt8(10))
{10, 0} undirected simple UInt8 graph
```
"""
function SimpleGraph{T}(n::Integer=0) where T <: Integer
    fadjlist = [Vector{T}() for _ = one(T):n]
    return SimpleGraph{T}(0, fadjlist)
end

# SimpleGraph(6), SimpleGraph(0x5)
SimpleGraph(n::T) where T <: Integer = SimpleGraph{T}(n)

# SimpleGraph()
SimpleGraph() = SimpleGraph{Int}()

# SimpleGraph(UInt8)
"""
    SimpleGraph(::Type{T})

Construct an empty `SimpleGraph{T}` with 0 vertices and 0 edges.

## Examples
```jldoctest
julia> SimpleGraph(UInt8)
{0, 0} undirected simple UInt8 graph
```
"""
SimpleGraph(::Type{T}) where T <: Integer = SimpleGraph{T}(zero(T))

# SimpleGraph(adjmx)
"""
    SimpleGraph{T}(adjm::AbstractMatrix)

Construct a `SimpleGraph{T}` from the adjacency matrix `adjm`.
If `adjm[i][j] != 0`, an edge `(i, j)` is inserted. `adjm` must be a square and symmetric matrix.
The element type `T` can be omitted.

## Examples
```jldoctest
julia> A1 = [false true; true false]
julia> SimpleGraph(A1)
{2, 1} undirected simple Int64 graph

julia> A2 = [2 7; 7 0]
julia> SimpleGraph{Int16}(A2)
{2, 2} undirected simple Int16 graph
```
"""
SimpleGraph(adjmx::AbstractMatrix) = SimpleGraph{Int}(adjmx)

# Graph{UInt8}(adjmx)
function SimpleGraph{T}(adjmx::AbstractMatrix) where T <: Integer
    dima, dimb = size(adjmx)
    isequal(dima, dimb) || throw(ArgumentError("Adjacency / distance matrices must be square"))
    issymmetric(adjmx) || throw(ArgumentError("Adjacency / distance matrices must be symmetric"))

    g = SimpleGraph(T(dima))
    @inbounds for i in findall(triu(adjmx) .!= 0)
        add_edge!(g, i[1], i[2])
    end
    return g
end

# SimpleGraph of a SimpleGraph
"""
    SimpleGraph{T}(g::SimpleGraph)

Construct a copy of g.
If the element type `T` is specified, the vertices of `g` are converted to this type.
Otherwise the element type is the same as for `g`.

## Examples
```jldoctest
julia> g = CompleteGraph(5)
julia> SimpleGraph{UInt8}(g)
{5, 10} undirected simple UInt8 graph
```
"""
SimpleGraph(g::SimpleGraph) = copy(g)

# converts Graph{Int} to Graph{Int32}
function SimpleGraph{T}(g::SimpleGraph) where T <: Integer
    h_fadj = [Vector{T}(x) for x in fadj(g)]
    return SimpleGraph(ne(g), h_fadj)
end

# SimpleGraph(digraph)
"""
    SimpleGraph(g::SimpleDiGraph)

Construct an undirected `SimpleGraph` from a directed `SimpleDiGraph`.
Every directed edge in `g` is added as an undirected edge.
The element type is the same as for `g`.

## Examples
```jldoctest
julia> g = PathDiGraph(Int8(5))
julia> SimpleGraph(g)
{5, 4} undirected simple Int8 graph
```
"""
function SimpleGraph(g::SimpleDiGraph)
    gnv = nv(g)
    edgect = 0
    newfadj = deepcopy(g.fadjlist)
    @inbounds for i in vertices(g)
        for j in badj(g, i)
            index = searchsortedfirst(newfadj[i], j)
            if index <= length(newfadj[i]) && newfadj[i][index] == j
                edgect += 1     # this is an existing edge - we already have it
                if i == j
                    edgect += 1 # need to count self loops
                end
            else
                insert!(newfadj[i], index, j)
                edgect += 2      # this is a new edge only in badjlist
            end
        end
    end
    iseven(edgect) || throw(AssertionError("invalid edgect in graph creation - please file bug report"))
    return SimpleGraph(edgect ÷ 2, newfadj)
end

@inbounds function cleanupedges!(fadjlist::Vector{Vector{T}}) where T <: Integer
    neg = 0
    for v in 1:length(fadjlist)
        if !issorted(fadjlist[v])
            sort!(fadjlist[v])
        end
        unique!(fadjlist[v])
        neg += length(fadjlist[v])
        # self-loops should count as one edge
        for w in fadjlist[v]
            if w == v
                neg += 1
                break
            end
        end
    end
    return neg ÷ 2
end

"""
    SimpleGraph(edge_list::Vector)

Construct a `SimpleGraph` from a vector of edges.
The element type is taken from the edges in `edge_list`.
The number of vertices is the highest that is used in an edge in `edge_list`.

### Implementation Notes
This constructor works the fastest when `edge_list` is sorted
by the lexical ordering and does not contain any duplicates.

### See also
[`SimpleGraphFromIterator`](@ref)

## Examples
```jldoctest

julia> el = Edge.([ (1, 2), (1, 5) ])
julia> SimpleGraph(el)
{5, 2} undirected simple Int64 graph
```
"""
function SimpleGraph(edge_list::Vector{SimpleGraphEdge{T}}) where T <: Integer
    nvg = zero(T)
    @inbounds(
    for e in edge_list
        nvg = max(nvg, src(e), dst(e))
    end)

    list_sizes = ones(Int, nvg)
    degs = zeros(Int, nvg)
    @inbounds(
    for e in edge_list
        s, d = src(e), dst(e)
        (s >= 1 && d >= 1) || continue
        degs[s] += 1
        if s != d
            degs[d] += 1
        end
    end)

    fadjlist = Vector{Vector{T}}(undef, nvg)
    @inbounds(
    for v in 1:nvg
        fadjlist[v] = Vector{T}(undef, degs[v])
    end)

    @inbounds(
    for e in edge_list
        s, d = src(e), dst(e)
        (s >= 1 && d >= 1) || continue
        fadjlist[s][list_sizes[s]] = d
        list_sizes[s] += 1
        if s != d
            fadjlist[d][list_sizes[d]] = s
            list_sizes[d] += 1
        end
    end)

    neg = cleanupedges!(fadjlist)
    g = SimpleGraph{T}()
    g.fadjlist = fadjlist
    g.ne = neg

    return g
end

@inbounds function add_to_fadjlist!(fadjlist::Vector{Vector{T}}, s::T, d::T) where T <: Integer
    nvg = length(fadjlist)
    nvg_new = max(nvg, s, d)
    for v = (nvg + 1):nvg_new
        push!(fadjlist, Vector{T}())
    end

    push!(fadjlist[s], d)
    if s != d
        push!(fadjlist[d], s)
    end
end

function _SimpleGraphFromIterator(iter)::SimpleGraph
    T = Union{}
    fadjlist = Vector{Vector{T}}()
    @inbounds(
    for e in iter
        typeof(e) <: SimpleGraphEdge ||
                        throw(ArgumentError("iter must be an iterator over SimpleEdge"))
        s, d = src(e), dst(e)
        (s >= 1 && d >= 1) || continue
        if T != eltype(e)
            T = typejoin(T, eltype(e))
            fadjlist = convert(Vector{Vector{T}}, fadjlist)
        end
        add_to_fadjlist!(fadjlist, s, d)
    end)

    T == Union{} && return SimpleGraph(0)
    neg  = cleanupedges!(fadjlist)
    g = SimpleGraph{T}()
    g.fadjlist = fadjlist
    g.ne = neg

    return g
end

function _SimpleGraphFromIterator(iter, ::Type{SimpleGraphEdge{T}}) where T <: Integer
    fadjlist = Vector{Vector{T}}()
    @inbounds(
    for e in iter
        s, d = src(e), dst(e)
        (s >= 1 && d >= 1) || continue
        add_to_fadjlist!(fadjlist, s, d)
    end)

    neg  = cleanupedges!(fadjlist)
    g = SimpleGraph{T}()
    g.fadjlist = fadjlist
    g.ne = neg

    return g
end

"""
    SimpleGraphFromIterator(iter)

Create a [`SimpleGraph`](@ref) from an iterator `iter`. The elements in iter must
be of `type <: SimpleEdge`.

# Examples
```jldoctest
julia> using LightGraphs

julia> g = SimpleGraph(3);

julia> add_edge!(g, 1, 2);

julia> add_edge!(g, 2, 3);

julia> h = SimpleGraphFromIterator(edges(g));

julia> collect(edges(h))
2-element Array{LightGraphs.SimpleGraphs.SimpleEdge{Int64},1}:
 Edge 1 => 2
 Edge 2 => 3
```
"""
function SimpleGraphFromIterator(iter)::SimpleGraph
    if Base.IteratorEltype(iter) == Base.EltypeUnknown()
        return _SimpleGraphFromIterator(iter)
    end
    if eltype(iter) <: SimpleGraphEdge && isconcretetype(eltype(iter))
        return _SimpleGraphFromIterator(iter, eltype(iter))
    end
    return _SimpleGraphFromIterator(iter)
end

edgetype(::SimpleGraph{T}) where T <: Integer = SimpleGraphEdge{T}

"""
    badj(g::SimpleGraph[, v::Integer])

Return the backwards adjacency list of a graph. If `v` is specified,
return only the adjacency list for that vertex.

###Implementation Notes
Returns a reference, not a copy. Do not modify result.
"""
badj(g::SimpleGraph) = fadj(g)
badj(g::SimpleGraph, v::Integer) = fadj(g, v)

"""
    adj(g[, v])

Return the adjacency list of a graph. If `v` is specified, return only the
adjacency list for that vertex.

### Implementation Notes
Returns a reference, not a copy. Do not modify result.
"""
adj(g::SimpleGraph) = fadj(g)
adj(g::SimpleGraph, v::Integer) = fadj(g, v)

copy(g::SimpleGraph) =  SimpleGraph(g.ne, deepcopy(g.fadjlist))

==(g::SimpleGraph, h::SimpleGraph) =
vertices(g) == vertices(h) &&
ne(g) == ne(h) &&
fadj(g) == fadj(h)

"""
    is_directed(g)

Return `true` if `g` is a directed graph.
"""
is_directed(::Type{SimpleGraph}) = false
is_directed(::Type{SimpleGraph{T}}) where T = false
is_directed(g::SimpleGraph) = false

function has_edge(g::SimpleGraph{T}, s, d) where T
    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 insorted(d, list_s)
end

function has_edge(g::SimpleGraph{T}, e::SimpleGraphEdge{T}) where T
    s, d = T.(Tuple(e))
    return has_edge(g, s, d)
end

"""
    add_edge!(g, e)

Add an edge `e` to graph `g`. Return `true` if edge was added successfully,
otherwise return `false`.

# Examples
```jldoctest
julia> using LightGraphs

julia> g = SimpleGraph(2);

julia> add_edge!(g, 1, 2)
true

julia> add_edge!(g, 2, 3)
false
```
"""
function add_edge!(g::SimpleGraph{T}, e::SimpleGraphEdge{T}) where T
    s, d = T.(Tuple(e))
    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 (index <= length(list) && list[index] == d) && return false  # edge already in graph
    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

"""
    rem_edge!(g, e)

Remove an edge `e` from graph `g`. Return `true` if edge was removed successfully,
otherwise return `false`.

### Implementation Notes
If `rem_edge!` returns `false`, the graph may be in an indeterminate state, as
there are multiple points where the function can exit with `false`.

# Examples
```jldoctest
julia> using LightGraphs

julia> g = SimpleGraph(2);

julia> add_edge!(g, 1, 2);

julia> rem_edge!(g, 1, 2)
true

julia> rem_edge!(g, 1, 2)
false
```
"""
function rem_edge!(g::SimpleGraph{T}, e::SimpleGraphEdge{T}) where T
    s, d = T.(Tuple(e))
    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 (index <= length(list) && list[index] == d) || return false  # edge not in graph
    deleteat!(list, index)

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

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

"""
    add_vertex!(g)

Add a new vertex to the graph `g`. Return `true` if addition was successful.

# Examples
```jldoctest
julia> using LightGraphs

julia> g = SimpleGraph(Int8(typemax(Int8) - 1))
{126, 0} undirected simple Int8 graph

julia> add_vertex!(g)
true

julia> add_vertex!(g)
false
```
"""
function add_vertex!(g::SimpleGraph{T}) where T
    (nv(g) + one(T) <= nv(g)) && return false       # test for overflow
    push!(g.fadjlist, Vector{T}())
    return true
end

"""
    rem_vertices!(g, vs, keep_order=false) -> vmap

Remove all vertices in `vs` from `g`.
Return a vector `vmap` that maps the vertices in the modified graph to the ones in
the unmodified graph.
If `keep_order` is `true`, the vertices in the modified graph appear in the same
order as they did in the unmodified graph. This might be slower.

### Implementation Notes
This function is not part of the official LightGraphs API and is subject to change/removal between major versions.

# Examples
```jldoctest
julia> using LightGraphs

julia> g = CompleteGraph{5}
{5, 10} undirected simple Int64 graph

julia> vmap = rem_vertices!(g, [2, 4], keep_order=true);

julia> vmap
3-element Array{Int64,1}:
 1
 3
 5

julia> g
{3, 3} undirected simple Int64 graph
```
"""
function rem_vertices!(g::SimpleGraph{T},
                       vs::AbstractVector{<:Integer};
                       keep_order::Bool=false
                      ) where {T <: Integer}
    # TODO There might be some room for performance improvements.
    # At the moment, we check for all edges if they stay in the graph.
    # If some vertices keep their position, this might be unnecessary.

    n = nv(g)
    isempty(vs) && return collect(Base.OneTo(n))

    # Sort and filter the vertices that we want to remove
    remove = sort(vs)
    unique!(remove)
    (1 <= remove[1] && remove[end] <= n) ||
            throw(ArgumentError("Vertices to be removed must be in the range 1:nv(g)."))

    # Create a vmap that maps vertices to their new position
    # vertices that get removed are mapped to 0
    vmap = Vector{T}(undef, n)

    if keep_order
        # traverse the vertex list and shift if a vertex gets removed
        i = 1
        @inbounds for u in vertices(g)
            if i <= length(remove) && u == remove[i]
                vmap[u] = 0
                i += 1
            else
                vmap[u] = u - (i - 1)
            end
        end
    else
        # traverse the vertex list and replace vertices that get removed
        # with the furthest one to the back that does not get removed
        i = 1
        j = length(remove)
        v = n
        @inbounds for u in vertices(g)
            u > v && break
            if i <= length(remove) && u == remove[i]
                while v == remove[j] && v > u
                   vmap[v] = 0
                   v -= one(T)
                   j -= 1
                end
                # v > remove[j] || u == v
                vmap[v] = u
                vmap[u] = 0
                v -= one(T)
                i += 1
            else
                vmap[u] = u
            end
        end
    end

    fadjlist = g.fadjlist

    # count the number of edges that will be removed
    # for an edge that gets removed we have to ensure that
    # such an edge does not get counted twice when both endpoints
    # get removed. That's why we relay on the ordering >= on the vertices.
    num_removed_edges = 0
    @inbounds for u in remove
        for v in fadjlist[u]
            if v >= u || vmap[v] != 0
                num_removed_edges += 1
            end
        end
    end
    g.ne -= num_removed_edges

    # move the lists in the adjacency list to their new position
    # The order of traversal is very important here, as otherwise we
    # could overwrite lists, that we want to keep!
    @inbounds for u in (keep_order ? (one(T):1:n) : (n:-1:one(T)))
        if vmap[u] != 0
            fadjlist[vmap[u]] = fadjlist[u]
        end
    end
    resize!(fadjlist, n - length(remove))

    # remove vertices from the lists in fadjlist
    @inbounds for list in fadjlist
        Δ = 0
        for (i, v) in enumerate(list)
            if vmap[v] == 0
                Δ += 1
            else
                list[i - Δ] = vmap[v]
            end
        end
        resize!(list, length(list) - Δ)
        if !keep_order
            sort!(list)
        end
    end

    # we create a reverse vmap, that maps vertices in the result graph
    # to the ones in the original graph. This resembles the output of
    # induced_subgraph
    reverse_vmap = Vector{T}(undef, nv(g))
    @inbounds for (i, u) in enumerate(vmap)
        if u != 0
            reverse_vmap[u] = i
        end
    end

    return reverse_vmap
end


text:
# Generates the dataset to be used in GMM example for EM course

using Plots, Distributions, Random, DataFrames, CSV
Random.seed!(123)

# Set GMM parameters
n = 200
α = 0.4
μ₀ = [-1, 2]
μ₁ = [3.0, 7.0]
Σ₀ = [2.0 1.0; 1.0 2.0]
Σ₁ = [1.5 0; 0 0.5]
𝒩₀ = MvNormal(μ₀, Σ₀)
𝒩₁ = MvNormal(μ₁, Σ₁)
X₁ = zeros(n)
X₂ = zeros(n)

for i=1:n
    X₁[i],X₂[i] = rand() < α ? rand(𝒩₁) : rand(𝒩₀)
end

# Store in data frame
df = DataFrame(X1 = X₁, X2 = X₂)

# Save to csv
CSV.write("gmm_observations.csv", df)

max_stars_repo_name: seanrattana/courses

max_issues_repo_name: seanrattana/courses

max_forks_repo_name: seanrattana/courses


text:
function juliacon2021()
    if myid() == 1
        return println(
            "Welcome to JuliaCon 2021! Find more information on https://juliacon.org/2021/."
        )
    else
        return println("Greetings from ", rand(countries), "!")
    end
    return nothing
end


text:
#=
This file is auto-generated. Do not edit.
=#

#! format: off

"""
    mutable struct EXST1 <: AVR
        Tr::Float64
        Vi_lim::NamedTuple{(:min, :max), Tuple{Float64, Float64}}
        Tc::Float64
        Tb::Float64
        Ka::Float64
        Ta::Float64
        Vr_lim::NamedTuple{(:min, :max), Tuple{Float64, Float64}}
        Kc::Float64
        Kf::Float64
        Tf::Float64
        V_ref::Float64
        ext::Dict{String, Any}
        states::Vector{Symbol}
        n_states::Int
        internal::InfrastructureSystemsInternal
    end

IEEE Type ST1 Excitation System (PTI version)

# Arguments
- `Tr::Float64`: Voltage Measurement Time Constant in s, validation range: `(0, nothing)`, action if invalid: `warn`
- `Vi_lim::NamedTuple{(:min, :max), Tuple{Float64, Float64}}`: Voltage input limits (Vi_min, Vi_max)
- `Tc::Float64`: Numerator lead-lag (lead) time constant in s, validation range: `(0, nothing)`, action if invalid: `warn`
- `Tb::Float64`: Denominator lead-lag (lag) time constant in s, validation range: `(0, nothing)`, action if invalid: `warn`
- `Ka::Float64`: Amplifier Gain, validation range: `(0, nothing)`, action if invalid: `warn`
- `Ta::Float64`: Amplifier Time Constant in s, validation range: `(0, nothing)`, action if invalid: `warn`
- `Vr_lim::NamedTuple{(:min, :max), Tuple{Float64, Float64}}`: Voltage regulator limits (regulator output) (Vr_min, Vr_max)
- `Kc::Float64`: Current field constant limiter multiplier, validation range: `(0, nothing)`
- `Kf::Float64`: Excitation control system stabilizer gain, validation range: `(eps(), 0.3)`, action if invalid: `warn`
- `Tf::Float64`: Excitation control system stabilizer time constant, validation range: `(eps(), nothing)`, action if invalid: `error`
- `V_ref::Float64`: Reference Voltage Set-point, validation range: `(0, nothing)`
- `ext::Dict{String, Any}`
- `states::Vector{Symbol}`: The states are:
	Vm: Sensed Terminal Voltage,
	Vrll: Lead-Lag state,
	Vr: Regulator Output,
	Vfb: Feedback state
- `n_states::Int`: The EXST1 has 4 states
- `internal::InfrastructureSystemsInternal`: power system internal reference, do not modify
"""
mutable struct EXST1 <: AVR
    "Voltage Measurement Time Constant in s"
    Tr::Float64
    "Voltage input limits (Vi_min, Vi_max)"
    Vi_lim::NamedTuple{(:min, :max), Tuple{Float64, Float64}}
    "Numerator lead-lag (lead) time constant in s"
    Tc::Float64
    "Denominator lead-lag (lag) time constant in s"
    Tb::Float64
    "Amplifier Gain"
    Ka::Float64
    "Amplifier Time Constant in s"
    Ta::Float64
    "Voltage regulator limits (regulator output) (Vr_min, Vr_max)"
    Vr_lim::NamedTuple{(:min, :max), Tuple{Float64, Float64}}
    "Current field constant limiter multiplier"
    Kc::Float64
    "Excitation control system stabilizer gain"
    Kf::Float64
    "Excitation control system stabilizer time constant"
    Tf::Float64
    "Reference Voltage Set-point"
    V_ref::Float64
    ext::Dict{String, Any}
    "The states are:
	Vm: Sensed Terminal Voltage,
	Vrll: Lead-Lag state,
	Vr: Regulator Output,
	Vfb: Feedback state"
    states::Vector{Symbol}
    "The EXST1 has 4 states"
    n_states::Int
    "power system internal reference, do not modify"
    internal::InfrastructureSystemsInternal
end

function EXST1(Tr, Vi_lim, Tc, Tb, Ka, Ta, Vr_lim, Kc, Kf, Tf, V_ref=1.0, ext=Dict{String, Any}(), )
    EXST1(Tr, Vi_lim, Tc, Tb, Ka, Ta, Vr_lim, Kc, Kf, Tf, V_ref, ext, [:Vm, :Vrll, :Vr, :Vfb], 4, InfrastructureSystemsInternal(), )
end

function EXST1(; Tr, Vi_lim, Tc, Tb, Ka, Ta, Vr_lim, Kc, Kf, Tf, V_ref=1.0, ext=Dict{String, Any}(), states=[:Vm, :Vrll, :Vr, :Vfb], n_states=4, internal=InfrastructureSystemsInternal(), )
    EXST1(Tr, Vi_lim, Tc, Tb, Ka, Ta, Vr_lim, Kc, Kf, Tf, V_ref, ext, states, n_states, internal, )
end

# Constructor for demo purposes; non-functional.
function EXST1(::Nothing)
    EXST1(;
        Tr=0,
        Vi_lim=(min=0.0, max=0.0),
        Tc=0,
        Tb=0,
        Ka=0,
        Ta=0,
        Vr_lim=(min=0.0, max=0.0),
        Kc=0,
        Kf=0,
        Tf=0,
        V_ref=0,
        ext=Dict{String, Any}(),
    )
end

"""Get [`EXST1`](@ref) `Tr`."""
get_Tr(value::EXST1) = value.Tr
"""Get [`EXST1`](@ref) `Vi_lim`."""
get_Vi_lim(value::EXST1) = value.Vi_lim
"""Get [`EXST1`](@ref) `Tc`."""
get_Tc(value::EXST1) = value.Tc
"""Get [`EXST1`](@ref) `Tb`."""
get_Tb(value::EXST1) = value.Tb
"""Get [`EXST1`](@ref) `Ka`."""
get_Ka(value::EXST1) = value.Ka
"""Get [`EXST1`](@ref) `Ta`."""
get_Ta(value::EXST1) = value.Ta
"""Get [`EXST1`](@ref) `Vr_lim`."""
get_Vr_lim(value::EXST1) = value.Vr_lim
"""Get [`EXST1`](@ref) `Kc`."""
get_Kc(value::EXST1) = value.Kc
"""Get [`EXST1`](@ref) `Kf`."""
get_Kf(value::EXST1) = value.Kf
"""Get [`EXST1`](@ref) `Tf`."""
get_Tf(value::EXST1) = value.Tf
"""Get [`EXST1`](@ref) `V_ref`."""
get_V_ref(value::EXST1) = value.V_ref
"""Get [`EXST1`](@ref) `ext`."""
get_ext(value::EXST1) = value.ext
"""Get [`EXST1`](@ref) `states`."""
get_states(value::EXST1) = value.states
"""Get [`EXST1`](@ref) `n_states`."""
get_n_states(value::EXST1) = value.n_states
"""Get [`EXST1`](@ref) `internal`."""
get_internal(value::EXST1) = value.internal

"""Set [`EXST1`](@ref) `Tr`."""
set_Tr!(value::EXST1, val) = value.Tr = val
"""Set [`EXST1`](@ref) `Vi_lim`."""
set_Vi_lim!(value::EXST1, val) = value.Vi_lim = val
"""Set [`EXST1`](@ref) `Tc`."""
set_Tc!(value::EXST1, val) = value.Tc = val
"""Set [`EXST1`](@ref) `Tb`."""
set_Tb!(value::EXST1, val) = value.Tb = val
"""Set [`EXST1`](@ref) `Ka`."""
set_Ka!(value::EXST1, val) = value.Ka = val
"""Set [`EXST1`](@ref) `Ta`."""
set_Ta!(value::EXST1, val) = value.Ta = val
"""Set [`EXST1`](@ref) `Vr_lim`."""
set_Vr_lim!(value::EXST1, val) = value.Vr_lim = val
"""Set [`EXST1`](@ref) `Kc`."""
set_Kc!(value::EXST1, val) = value.Kc = val
"""Set [`EXST1`](@ref) `Kf`."""
set_Kf!(value::EXST1, val) = value.Kf = val
"""Set [`EXST1`](@ref) `Tf`."""
set_Tf!(value::EXST1, val) = value.Tf = val
"""Set [`EXST1`](@ref) `V_ref`."""
set_V_ref!(value::EXST1, val) = value.V_ref = val
"""Set [`EXST1`](@ref) `ext`."""
set_ext!(value::EXST1, val) = value.ext = val


text:
toradians(degree::Float64) = degree * π / 180.0
todegrees(radian::Float64) = radian * 180.0 / π
convertunits(units::Symbol = :kilometers) = 3960.0*(units == :miles) + 6373.0*(units == :kilometers) + 57.2957795*(units == :degrees) + (units == :radians)

function distance(point1::Vector{Float64}, point2::Vector{Float64}, units::Symbol = :kilometers)
    dLat = toradians(point2[2] - point1[2])
    dLon = toradians(point2[1] - point1[1])
    lat1 = toradians(point1[2])
    lat2 = toradians(point2[2])
    a = sin(dLat/2)^2 + sin(dLon/2)^2 * cos(lat1) * cos(lat2)
    2.0 * atan2(sqrt(a), sqrt(1-a)) * convertunits(units)
end
distance(point1::GeoJSON.Point, point2::GeoJSON.Point, units::Symbol = :kilometers) = distance(point1.coordinates, point2.coordinates, units)
distance(point1::GeoJSON.Feature, point2::GeoJSON.Feature, units::Symbol = :kilometers) = distance(point1.geometry, point2.geometry, units)

function ringarea(coords::Vector{Vector{Float64}})
    f(p1::Vector{Float64},p2::Vector{Float64}) = toradians(p2[1]-p1[1]) * (2+sin(toradians(p1[2])) + sin(toradians(p2[2])))
    area_ = 0.0
    if length(coords) > 2
        for i=1:length(coords)-1
            area_ += f(coords[1], coords[i+1])
        end
        area_ *= RADIUS^2 / 2
    end
    area_
end

function polygonarea(coords::Vector{Vector{Vector{Float64}}})
    if length(coords) > 0
        area_ = abs(ringarea(coords[0]))
        for coord in coords[2:end]
            area_ -= abs(ringarea(coords))
        end
        return area_
    else
        return 0.0
    end
end
polygonarea(obj::Void) = 0.0

area(obj::GeoJSON.Point) = 0.0
area(obj::GeoJSON.MultiPoint) = 0.0
area(obj::GeoJSON.LineString) = 0.0
area(obj::GeoJSON.MultiLineString) = 0.0
area(obj::GeoJSON.Polygon) = polygonarea(obj.coordinates)
area(obj::GeoJSON.MultiPolygon) = sum(map(polygonarea, obj.coordinates))
area(obj::GeoJSON.GeometryCollection) = sum(map(area, obj.geometries))
area(obj::GeoJSON.Feature) = area(obj.geometry)
area(obj::GeoJSON.FeatureCollection) = sum(map(area, obj.features))

function nearest(target::GeoJSON.Point, points::GeoJSON.FeatureCollection)
    if length(points.features) > 0
        nearest_feature = points.features[1]
        nearest_dist = distance(target, nearest_feature.geometry)
        for feature in points.features[2:end]
            dist = distance(target, feature.geometry)
            if dist < nearest_dist
                nearest_feature = feature
                nearest_dist = dist
            end
        end
        return deepcopy(nearest_feature)
    end
end

function nearest(target::GeoJSON.Point, points::GeoJSON.GeometryCollection)
    if length(points.geometries) > 0
        nearest = points.geometries[1]
        nearest_dist = distance(target, nearest.geometry)
        for geometry in points.geometries[2:end]
            dist = distance(target, geometry)
            if dist < nearest_dist
                nearest = geometry
                nearest_dist = dist
            end
        end
        nearest = deepcopy(nearest)
    end
end
nearest(target::GeoJSON.Feature, points::GeoJSON.FeatureCollection) = nearest(target.geometry, points)

function bboxpolygon(bbox::Vector{Float64})
    lowleft = Float64[bbox[1], bbox[2]]
    topleft = Float64[bbox[1], bbox[4]]
    topright = Float64[bbox[3], bbox[4]]
    lowright = Float64[bbox[3], bbox[2]]
    GeoJSON.Polygon(Vector{Vector{Float64}}[Vector{Float64}[lowleft, lowright, topright, topleft, lowleft]])
end

topright(obj::Vector{Vector{Float64}}) = reduce(max, obj)
topright(obj::Vector{Vector{Vector{Float64}}}) = topright(map(topright, obj)::Vector{Vector{Float64}})
topright(obj::GeoJSON.Point) = obj.coordinates
topright(obj::GeoJSON.MultiPoint) = topright(obj.coordinates)
topright(obj::GeoJSON.LineString) = topright(obj.coordinates)
topright(obj::GeoJSON.MultiLineString) = topright(obj.coordinates)
topright(obj::GeoJSON.Polygon) = topright(obj.coordinates)
topright(obj::GeoJSON.MultiPolygon) = topright(map(topright, obj.coordinates))
topright(obj::GeoJSON.GeometryCollection) = topright(map(topright, obj.geometries))
topright(obj::GeoJSON.Feature) = topright(obj.geometry)
topright(obj::GeoJSON.FeatureCollection) = topright(map(topright, obj.features))

lowleft(obj::Vector{Vector{Float64}}) = reduce(min, obj)
lowleft(obj::Vector{Vector{Vector{Float64}}}) = lowleft(map(lowleft, obj)::Vector{Vector{Float64}})
lowleft(obj::GeoJSON.Point) = obj.coordinates
lowleft(obj::GeoJSON.MultiPoint) = lowleft(obj.coordinates)
lowleft(obj::GeoJSON.LineString) = lowleft(obj.coordinates)
lowleft(obj::GeoJSON.MultiLineString) = lowleft(obj.coordinates)
lowleft(obj::GeoJSON.Polygon) = lowleft(obj.coordinates)
lowleft(obj::GeoJSON.MultiPolygon) = lowleft(map(lowleft, obj.coordinates))
lowleft(obj::GeoJSON.GeometryCollection) = lowleft(map(lowleft, obj.geometries))
lowleft(obj::GeoJSON.Feature) = lowleft(obj.geometry)
lowleft(obj::GeoJSON.FeatureCollection) = lowleft(map(lowleft, obj.features))

extent(obj::GeoJSON.Point) = obj.coordinates
extent(obj::GeoJSON.MultiPoint) = [topright(obj); lowleft(obj)]
extent(obj::GeoJSON.LineString) = [topright(obj); lowleft(obj)]
extent(obj::GeoJSON.MultiLineString) = [topright(obj); lowleft(obj)]
extent(obj::GeoJSON.Polygon) = [topright(obj); lowleft(obj)]
extent(obj::GeoJSON.MultiPolygon) = [topright(obj); lowleft(obj)]
extent(obj::GeoJSON.Feature) = [topright(obj); lowleft(obj)]
extent(obj::GeoJSON.FeatureCollection) = [topright(obj); lowleft(obj)]

envelope(obj::GeoJSON.AbstractGeoJSON) = bboxpolygon(extent(obj))

function midpoint(p1::Vector{Float64}, p2::Vector{Float64})
    x1, y1, x2, y2 = p1[1], p1[2], p2[1], p2[2]
    Float64[(x1+x2)/2, (y1+y2)/2]
end
midpoint(p1::GeoJSON.Point, p2::GeoJSON.Point) = GeoJSON.Point(midpoint(p1.coordinates, p2.coordinates))
midpoint(p1::GeoJSON.Feature, p2::GeoJSON.Feature) = GeoJSON.Feature(midpoint(p1.geometry, p2.geometry))

function square(bbox::Vector{Float64})
    lowleft = Float64[bbox[1],bbox[2]]
    topleft = Float64[bbox[1],bbox[4]]
    topright = Float64[bbox[3],bbox[4]]
    lowright = Float64[bbox[3],bbox[2]]
    hdist = distance(lowleft, lowright)
    vdist = distance(lowleft, topleft)
    if hdist >= vdist
        mid = midpoint(lowleft, topleft)
        return Float64[ bbox[1],
                        mid[2] - (bbox[3] - bbox[1])/2,
                        bbox[3],
                        mid[1] + (bbox[3] - bbox[1])/2 ]
    else
        mid = midpoint(lowleft, lowright)
        return Float64[ mid[1] - (bbox[4] - bbox[2])/2,
                        bbox[2],
                        mid[1] + (bbox[4] - bbox[2])/2,
                        bbox[4] ]
    end
end

function size(bbox::Vector{Float64}, factor::Real)
    curr_xdist = bbox[3] - bbox[1]
    curr_ydist = bbox[4] - bbox[2]
    new_xdist = factor * curr_xdist
    new_ydist = factor * curr_ydist
    dx = new_xdist - curr_xdist
    dy = new_ydist - curr_ydist

    Float64[bbox[1] - dx/2, # low_x
            bbox[2] - dy/2, # low_y
            bbox[3] + dx/2, # high_x
            bbox[4] + dy/2] # high_y
end

function center(obj::GeoJSON.AbstractGeoJSON)
    ext = extent(obj)
    GeoJSON.Point(Float64[(ext[1] + ext[3])/2, (ext[2] + ext[4])/2])
end

centroid(obj::Vector{Vector{Float64}}) = reduce(+, obj) / length(obj)
function centroid(obj::Vector{Vector{Vector{Float64}}})
    points = vcat(obj...) # not memory efficient
    reduce(+, points) / length(points)
end
centroid(obj::GeoJSON.Point) = obj
centroid(obj::GeoJSON.MultiPoint) = GeoJSON.Point(centroid(obj.coordinates))
centroid(obj::GeoJSON.LineString) = GeoJSON.Point(centroid(obj.coordinates))
centroid(obj::GeoJSON.MultiLineString) = GeoJSON.Point(centroid(obj.coordinates))
centroid(obj::GeoJSON.Polygon) = GeoJSON.Point(centroid(obj.coordinates)) # what about the holes?
centroid(obj::GeoJSON.MultiPolygon) = GeoJSON.Point(centroid(vcat(obj.coordinates...))) # what about the holes?
centroid(obj::GeoJSON.GeometryCollection) = map(centroid, obj.geometries) # different from Turf.js
centroid(obj::GeoJSON.Feature) = centroid(obj.geometry)
centroid(obj::GeoJSON.FeatureCollection) = map(centroid, obj.features) # different from Turf.js

function bearing(p1::Vector{Float64}, p2::Vector{Float64})
    lon1 = toradians(p1[1])
    lon2 = toradians(p2[1])
    lat1 = toradians(p1[2])
    lat2 = toradians(p2[2])
    a = sin(lon2 - lon1) * cos(lat2)
    b = cos(lat1) * sin(lat2) - sin(lat1) * cos(lat2) * cos(lon2 - lon1)
    todegrees(atan2(a, b))
end
bearing(p1::GeoJSON.Point, p2::GeoJSON.Point) = bearing(p1.coordinates, p2.coordinates)
bearing(p1::GeoJSON.Feature, p2::GeoJSON.Feature) = bearing(p1.geometry, p2.geometry)

function destination(point::Vector{Float64}, distance::Float64, bearing::Float64, units::Symbol = :kilometers)
    lon1 = toradians(point[1])
    lat1 = toradians(point[2])
    bearing_rad = toradians(bearing)
    dist = distance / convertunits(units)
    lat2 = asin(sin(lat1)*cos(dist) + cos(lat1)*sin(dist)*cos(bearing_rad))
    lon2 = lon1 + atan2(sin(bearing_rad)*sin(dist)*cos(lat1), cos(dist) - sin(lat1)*sin(lat2))
    GeoJSON.Point(Float64[todegrees(lon2), todegrees(lat2)])
end
destination(point::GeoJSON.Point, distance::Float64, bearing::Float64, units::Symbol = :kilometers) = destination(point.coordinates, distance, bearing, units)
destination(point::GeoJSON.Feature, distance::Float64, bearing::Float64, units::Symbol = :kilometers) = destination(point.geometry, distance, bearing, units)

function linedistance(line::Vector{Vector{Float64}}, units::Symbol = :kilometers)
    dist = 0.0
    for i=1:length(line)-1
        dist += distance(line[i], line[i+1], units)
    end
    dist
end
linedistance(line::GeoJSON.LineString, units::Symbol = :kilometers) = linedistance(line.coordinates, units)
function linedistance(multiline::GeoJSON.MultiLineString, units::Symbol = :kilometers)
    dist = 0.0
    for line in multiline.coordinates
        dist += linedistance(line, units)
    end
    dist
end
linedistance(line::GeoJSON.Feature, units::Symbol = :kilometers) = linedistance(line.geometry, units)

function along(line::Vector{Vector{Float64}}, dist::Real, units::Symbol = :kilometers)
    travelled = 0.0
    n_points = length(line)
    for i=1:n_points
        if travelled <= dist && i == n_points
            break
        elseif travelled == dist
            return GeoJSON.Point(line[i])
        elseif travelled > dist
            overshot = dist - travelled
            direction = bearing(line[i], line[i-1]) - 180.0
            return destination(line[i], overshot, direction, units)
        else # travelled == dist
            travelled += distance(line[i], line[i+1], units)
        end
    end
    GeoJSON.Point(line[end])
end
along(line::GeoJSON.LineString, dist::Real, units::Symbol = :kilometers) = along(line.coordinates, dist, units)
along(line::GeoJSON.Feature, dist::Real, units::Symbol = :kilometers) = along(line.geometry, dist, units)


text:
struct CardinalDict{K,V} <: Associative{K,V}
    valued::BitArray{1}
    values::Vector{V}

    function CardinalDict{V}(n::K) where K<:Integer where V
        valued = falses(n)
        values = Vector{V}(n)
        T = type_for_indexing(n)
        return new{T,V}(valued, values)
    end
end

function CardinalDict(values::Vector{T}) where T
    n = length(values)
    dict = CardinalDict{T}(n)
    for i in 1:n
        @inbounds dict[i] = values[i]
    end
    return dict
end

const SInt = Union{Int8, Int16, Int32, Int64, Int128}

# these two support eval(parse(string(dict::CardinalDict)))

function CardinalDict(zipped::Base.Iterators.Zip2{Array{I,1},Array{V,1}}, maxkey::SInt=length(zipped)) where {I<:SInt,V}
    thekeys = map(first, zipped)
    thevals = map(last, zipped)
    maxkey  = max(maxkey, maximum(thekeys))
    result = CardinalDict{V}(maxkey)
    for (k,v) in zipped
        result[k] = v
    end
    return result
end

function CardinalDict(pairs::Vector{Pair{I,V}}, maxkey::SInt=length(pairs)) where {I<:SInt,V}
    thekeys = map(first, pairs)
    thevals = map(last, pairs)
    maxkey  = max(maxkey, maximum(thekeys))
    result = CardinalDict{V}(maxkey)
    for (k,v) in zip(thekeys, thevals)
        result[k] = v
    end
    return result
end

# direct manipulations

@inline function Base.haskey(dict::CardinalDict{K,V}, key::K) where {K,V}
   return getindex(dict.valued, key)
end
@inline Base.haskey(dict::CardinalDict{K,V}, key::J) where {J,K,V} =
    haskey(dict, key%K)

function Base.getindex(dict::CardinalDict{K,V}, key::K) where {K,V}
    haskey(dict, key) || throw(ErrorException("Key (index) $(key) has not been given a value"))
    @inbounds return getindex(dict.values, key)
end
@inline Base.getindex(dict::CardinalDict{K,V}, key::J) where {J,K,V} =
    getindex(dict, key%K)

function Base.setindex!(dict::CardinalDict{K,V}, value::V, key::K) where {K,V}
    0 < key <= keymax(dict) || throw(ErrorException("Key (index) $(key) is outside of the domain 1:$(keymax(dict))."))
    @inbounds begin
        setindex!(dict.valued, true, key)
        setindex!(dict.values, value, key)
    end
    return dict
end
@inline Base.setindex!(dict::CardinalDict{K,V}, value::V, key::J) where J where K where V =
    setindex!(dict, value, key%K)


text:
abstract type AbstractDiffn end

mutable struct DiffnPoly2nd<: AbstractDiffn
end
function (obj::DiffnPoly2nd)(
                             xy::AbstractArray{Tmsh},
		 lambda::AbstractArray{Tsol, 3}) where {Tmsh, Tsol}
	# @assert(size(lambda, 1) == Tdim)
	# @assert(size(lambda, 2) == Tdim)
	# the 3rd dimension should be dof per node
	for dof = 1:size(lambda, 3)
		lambda[1, 1, dof] = xy[1]*xy[1] + 1.0
		lambda[1, 2, dof] = xy[1]*xy[2]
		lambda[2, 2, dof] = xy[2]*xy[2] + 1.0
		lambda[2, 1, dof] = lambda[1, 2, dof]
	end
	return nothing
end

mutable struct DiffnPoly6th<: AbstractDiffn
end
function (obj::DiffnPoly6th)(
                             xy::AbstractArray{Tmsh},
		 lambda::AbstractArray{Tsol, 3}) where {Tmsh, Tsol}
	# @assert(size(lambda, 1) == Tdim)
	# @assert(size(lambda, 2) == Tdim)
	# the 3rd dimension should be dof per node
	for dof = 1:size(lambda, 3)
		lambda[1, 1, dof] = xy[1]^6 + 1.0
		lambda[1, 2, dof] = xy[1]^3*xy[2]^3
		lambda[2, 2, dof] = xy[2]^6 + 1.0
		lambda[2, 1, dof] = lambda[1, 2, dof]
	end
	return nothing
end

mutable struct DiffnPoly0th<: AbstractDiffn
end

function (obj::DiffnPoly0th)(
                              xy::AbstractArray{Tmsh},
		  lambda::AbstractArray{Tsol, 3}) where {Tmsh, Tsol}
	# @assert(size(lambda, 1) == Tdim)
	# @assert(size(lambda, 2) == Tdim)
	# the 3rd dimension should be dof per node
	for dof = 1:size(lambda, 3)
		lambda[1, 1, dof] = 1.0e1
		lambda[2, 2, dof] = 1.0e1
		lambda[1, 2, dof] = 0.0
		lambda[2, 1, dof] = 0.0
	end
	return nothing
end

global const DiffnDict = Dict{String, AbstractDiffn}(
	"poly0th" => DiffnPoly0th(),
	"poly2nd" => DiffnPoly2nd(),
	"poly6th" => DiffnPoly6th()
)

function getDiffnFunc(mesh::AbstractMesh,
					  sbp::AbstractOperator,
					  eqn::EllipticData,
					  opts)
	eqn.diffusion_func = DiffnDict[opts["Diffusion"]]
end

function calcDiffn(mesh::AbstractMesh{Tmsh},
				   sbo::AbstractOperator,
				   eqn::EllipticData{Tsol, Tres, Tdim},
				   diffusion_func::AbstractDiffn,
				   lambda::AbstractArray{Tres, 5} ) where {Tmsh, Tsol, Tres, Tdim}

	numElems = mesh.numEl
	numNodesPerElement = mesh.numNodesPerElement

	for elem = 1:numElems
		for n = 1:numNodesPerElement
			xy = sview(mesh.coords, :, n, elem)

			diffusion_func(xy, sview(lambda, :, :, :, n, elem))
		end
	end
end


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Tests related to the string width computation.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

@testset "Printable text width" verbose = true begin
    # Some examples here were obtained from:
    #   https://www.lihaoyi.com/post/BuildyourownCommandLinewithANSIescapecodes.html

    str = "\e[38;5;231;48;5;243mTes\e[3mt Test\e[1m Test\e[0m"
    @test printable_textwidth(str) == 14

    str = "\e[38;5;231;48;5;243m😃\e[3m😄\e[1m😊\e[0m"
    @test printable_textwidth(str) == 6

    str = "\e[1m😃\e[0m\e[4m😅\e[0m\e[7m🥳\e[0m"
    @test printable_textwidth(str) == 6

    str = """
    \e[38;5;231;48;5;243mTes\e[3mt Test\e[1m Test\e[0m
    \e[38;5;231;48;5;243m😃\e[3m😄\e[1m😊\e[0m
    \e[1m😃\e[0m\e[4m😅\e[0m\e[7m🥳\e[0m
    \u001b[30;1m A \u001b[31;1m B \u001b[32;1m C \u001b[33;1m D \u001b[0m
    \u001b[44;1m A \u001b[45;1m B \u001b[46;1m C \u001b[47;1m D \u001b[0m
    \u001b[1m BOLD \u001b[0m\u001b[4m Underline \u001b[0m\u001b[7m Reversed \u001b[0m
    \u001b[1m\u001b[4m\u001b[7m BOLD Underline Reversed \u001b[0m"""

    lines_width = printable_textwidth_per_line(str)

    @test length(lines_width) == 7
    @test lines_width[1] == 14
    @test lines_width[2] == 6
    @test lines_width[3] == 6
    @test lines_width[4] == 12
    @test lines_width[5] == 12
    @test lines_width[6] == 27
    @test lines_width[7] == 25
end


text:
# TODO: add tests for PMmap

GEN_EQ_FIELDS = [
    "gen_bus",
    "mbase",
    "pmax",
    "pmin",
    "startup",
    "shutdown",
    "ncost",
    "cost",
    "gen_status",
    "model",
]#"pg", "qg", "vg", "qmax", "qmin"]
BUS_EQ_FIELDS = ["bus_i", "index", "base_kv", "vmin", "vmax"]
BRANCH_EQ_FIELDS = ["f_bus", "t_bus", "br_status", "br_r", "br_x", "b_fr", "b_to"]#"angmin", "angmax",]#"g_to", "g_fr", "shift", "transformer", "tap", "rate_a", "rate_b", "rate_c"]
COMPARISON_CATS = [
    "bus",
    "branch",
    "baseMVA",
    "per_unit",
    "storage",
    "dc_line",
    "gen",
    #"switch", #not represented in PSY
    "areas",
    "shunt",
    "load",
    "name",
    "source_type",
    "source_version",
]

function compare_devices(pm_devices, pmi_devices, fields, match_expr)
    for (ix, pm_device) in pm_devices
        pmi_device_match = [d for (k, d) in pmi_devices if match_expr(d, pm_device)]
        @assert length(pmi_device_match) == 1
        for pmi_device in pmi_device_match
            match = true
            for field in fields
                if haskey(pm_device, field)
                    #@show ix, pmi_device["index"], field, pm_device[field] == pmi_device[field]
                    match = isapprox(pm_device[field], pmi_device[field], atol = 0.01)
                    !match && @show ix, field, pm_device[field], pmi_device[field]
                    !match && break
                end
            end
            @test match
        end
    end
end

function get_pm_br_ids(pmi_devices)
    pm_br_ids = Dict()
    for (k, br) in pmi_devices
        br_type = first(br["source_id"])
        if br_type == "transformer"
            id = last(split(br["source_id"][end - 1], "_"))
        else
            id = replace(last(split(br["source_id"][end], "_")), "Branch " => "")
        end
        pm_br_ids[k] = id
    end
    return pm_br_ids
end

function compare_pm_pmi(fname::String)
    pm_data = try
        PM.parse_file(fname)
    catch
        @warn "skipping testing of $fname due to PowerModels parsing error"
        return
    end
    sys = System(fname)
    pmi_data = PMI.get_pm_data(sys)
    pm_map = PMI.get_pm_map(sys)

    # there is no guarantee that we represent all data categories that PM does, so filter by the ones we care about
    pm_cats = intersect(COMPARISON_CATS, keys(pm_data))
    endswith(fname, joinpath("matpower", "case3.m")) && setdiff!(pm_cats, ["areas"]) #dont compare areas for case3.m, see: https://github.com/lanl-ansi/PowerModels.jl/pull/788

    @test isempty(setdiff(pm_cats, keys(pmi_data)))
    for cat in pm_cats
        pm_devices = pm_data[cat]
        @testset "$cat" begin
            @test haskey(pmi_data, cat)
            pmi_devices = get(pmi_data, cat, nothing)
            if !(typeof(pm_devices) <: AbstractString)
                @test isempty(setdiff(keys(pm_devices), keys(pmi_devices)))
                if cat == "gen"
                    compare_devices(
                        pm_devices,
                        pmi_devices,
                        GEN_EQ_FIELDS,
                        (d, pmd) ->
                            join(strip.(string.(d["source_id"])), "-") ==
                            join(strip.(string.(pmd["source_id"])), "-"),
                    )
                    @test isempty(setdiff(keys(pmi_devices), keys(pm_map[cat])))
                elseif cat == "bus"
                    compare_devices(
                        pm_devices,
                        pmi_devices,
                        BUS_EQ_FIELDS,
                        (d, pmd) -> d["index"] == pmd["index"],
                    )
                    @test isempty(setdiff(keys(pmi_devices), keys(pm_map[cat])))
                elseif cat == "branch"
                    pm_br_ids = get_pm_br_ids(pmi_devices)
                    compare_devices(
                        pm_devices,
                        pmi_devices,
                        BRANCH_EQ_FIELDS,
                        (d, pmd) -> pm_br_ids[string(d["index"])] == string(pmd["index"]),
                    )
                    @test isempty(setdiff(keys(pmi_devices), keys(pm_map[cat])))
                end
            end
        end
    end
end

@testset "Test PMI and PM data equivelance" begin
    for (root, dirs, files) in walkdir(PM_DATA_DIR)
        for file in files
            if !occursin("parser_test", file)
                @testset "$file" begin
                    compare_pm_pmi(joinpath(root, file))
                end
            end
        end
    end
end

@testset "Test RTS-GMLC tabular data" begin
    sys = PSB.build_system(PSITestSystems, "RTS_GMLC_sys")
    pm_data = PMI.get_pm_data(sys)
    @test length(pm_data["bus"]) == length(PSY.get_components(PSY.Bus, sys))
    @test length(pm_data["branch"]) == length(PSY.get_components(PSY.ACBranch, sys))
    @test length(pm_data["gen"]) == length(PSY.get_components(PSY.Generator, sys))
    @test length(pm_data["storage"]) == length(PSY.get_components(PSY.Storage, sys))
    @test length(pm_data["dcline"]) == length(PSY.get_components(PSY.DCBranch, sys))
    @test length(pm_data["areas"]) == length(PSY.get_components(PSY.Area, sys))
    @test length(pm_data["load"]) == length(PSY.get_components(PSY.PowerLoad, sys))
    @test length(pm_data["shunt"]) == length(PSY.get_components(PSY.FixedAdmittance, sys))
end

max_stars_repo_name: NREL-SIIP/PowerModelsInterface

max_issues_repo_name: NREL-SIIP/PowerModelsInterface

max_forks_repo_name: NREL-SIIP/PowerModelsInterface


text:
using Printf

import Base: getindex, ndims, length, *
import LinearAlgebra: ×

abstract type AbstractDomain end

"""
    IntervalDomain
"""
struct IntervalDomain{T} <: AbstractDomain
    min::T
    max::T
    periodic::Bool
end

function IntervalDomain(; min, max, periodic = false)
    @assert min < max
    return IntervalDomain(min, max, periodic)
end

"""
    ProductDomain
"""
struct ProductDomain{T} <: AbstractDomain
    domains::T
end

function ProductDomain(; domains)
    if length(domains) < 2
        error("A product domain needs more than one subdomain!")
    else
        return reduce(×, domains)
    end
end

"""
    SphericalShell
"""
struct SphericalShell{T} <: AbstractDomain
    radius::T
    height::T
end

function SphericalShell(; radius, height)
    @assert radius > 0 && height > 0
    return SphericalShell(radius, height)
end

"""
    Nary Operations
"""
×(domain1::AbstractDomain, domain2::AbstractDomain) = ProductDomain((domain1, domain2))
×(product::ProductDomain, domain::AbstractDomain) = ProductDomain((product.domains..., domain))
×(domain::AbstractDomain, product::ProductDomain) = ProductDomain((domain, product.domains...))
×(product1::ProductDomain, product2::ProductDomain) = ProductDomain((product1.domains..., product2.domains...))

"""
    Extensions
"""
Base.ndims(domain::IntervalDomain) = 1
Base.ndims(domain::ProductDomain) = +(ndims.(domain.domains)...)
Base.ndims(domain::SphericalShell) = 3

Base.length(domain::IntervalDomain) = domain.max - domain.min
Base.length(domain::ProductDomain) = length.(domain.domains)
Base.length(domain::SphericalShell) = (domain.radius, domain.radius, domain.height)

Base.getindex(domain::ProductDomain, i::Int) = domain.domains[i]

function Base.getproperty(value::ProductDomain, name::Symbol)
    if name == :periodicity
        periodicity = ones(Bool, ndims(value))
        for i in 1:ndims(value)
            periodicity[i] = value[i].periodic
        end
        return Tuple(periodicity)
    else
        # default
        return getfield(value, name)
    end
end

function Base.propertynames(::ProductDomain)
    return (:domains, :periodicity)
end

function Base.show(io::IO, domain::IntervalDomain)
    min = domain.min
    max = domain.max
    printstyled(io, "[", color = 226)
    astring = @sprintf("%0.2f", min)
    bstring = @sprintf("%0.2f", max)
    printstyled(astring, ", ", bstring, color = 7)
    domain.periodic ? printstyled(io, ")", color = 226) : printstyled(io, "]", color = 226)
end

function Base.show(io::IO, product::ProductDomain)
    for (i, domain) in enumerate(product.domains)
        print(domain)
        if i < ndims(product)
            printstyled(io, "×", color = 118)
        end
    end
end

function info(domain::ProductDomain)
    println("This is a ", ndims(domain), "-dimensional tensor product domain.")
    print("The domain is ")
    println(domain, ".")
    for (i, domain) in enumerate(domain.domains)
        domain_string = domain.periodic ? "periodic" : "wall-bounded"
        length = @sprintf("%.2f ", domain.max - domain.min)
        println("The dimension $i domain is ", domain_string, " with length ≈ ", length)
    end
end

max_stars_repo_name: CliMA/CouplerMachine

max_issues_repo_name: CliMA/CouplerMachine


text:
# Prints out a list of outliers for best-fit model:
using JLD2
using PyPlot
using Printf

include("../src/NbodyGradient/src/ttv.jl")

#@load "output_files/T1_run_hmc_student_ecc_lgndof_V1exp2nuinv_nstep2000_eps0.1_nleap20_201.jld2"
@load "T1_timing_posterior.jld2"
@load "T1_likelihood_profile_student_planetc_3.0sig.jld2"

# Print deviations of each transit:
chi_planet = zeros(Float64,n-1)
nobs_planet = zeros(Int64,n-1)
dev = zeros(Float64,ntrans)
elements = elements_grid[1,:,:,jgrid[1]]
dq0 = ttv_elements!(n,t0,h,tmax,elements,tt1,count1,0.0,0,0,rstar)
println("| Transit number | Planet | Model time | Observed time | Obs - mod | (Obs-mod)/sig | (Obs-mod)^2/sig^2 |")
println("| ---| --- | --- | --- | ---- | --- | --- |")
dur = [36.06,42.03,48.87,55.76,62.85,68.24,76.16] ./(24*60)
for j=1:ntrans
  dev[j] = (tobs_tot[j]-tt1[iplanet[j],indx[j]])/sobs_tot[j]
  chi_planet[iplanet[j]-1] += dev[j]^2
  nobs_planet[iplanet[j]-1] += 1
  if abs(dev[j]) > 3
    overlap = false
    for k=1:ntrans
      if abs(tt1[iplanet[k],indx[k]] - tt1[iplanet[j],indx[j]]) < (dur[iplanet[j]-1]+dur[iplanet[k]-1])/2 && k != j
        overlap = true
      end
    end
    println("|",@sprintf("%4i",j),"| ",@sprintf("%2i",iplanet[j]-1),"| ",@sprintf("%4i",indx[j]),"| ",
      @sprintf("%12.5f",tt1[iplanet[j],indx[j]]),"| ",@sprintf("%12.5f",tobs_tot[j]),"| ",
      @sprintf("%8.5f",tobs_tot[j]-tt1[iplanet[j],indx[j]]),"| ",
      @sprintf("%8.5f",sobs_tot[j]),"| ",@sprintf("%5.1f",dev[j]),"| ",@sprintf("%5.1f",dev[j]^2),"|",overlap)
  end
end


text:
# Kinematic julia code for the purposes of analyzing telemtric data pertaining to speed,
# acceleration and other measurements dealing of the kinematic nature.

using Dates, KernelDensity, RecipesBase, Rmath, DataFrames, Query

struct KMH <: SpeedType end;
struct MPH <: SpeedType end;

"""
Kinematic Event Structure

A basic structure for velocity, acceleration and timestamp.

Velocity must be positive!

Constructor:

	KinematicEvent(vel::Float64, aec::Float64, t::DateTime, speedtype::SpeedType)

"""
mutable struct KinematicEvent <: TelemetricEvent
    vel::Float64 # velocity
    aec::Float64 # acceleration
    t::DateTime # Time stamp
    speedtype::SpeedType # speed type
    function KinematicEvent(vel::Float64, aec::Float64, t::DateTime, speedtype::SpeedType)
        # Check to make sure velocity is positive
        if( vel < 0)
            error("Velocity must be positive")
        else
            ke = new(vel,aec,t, speedtype);
            return(ke);
        end
    end
end

function Base.show(io::IO, ke::KinematicEvent)
    println(io, "---------------------------")
    println(io, "KinematicEvent object");
    println(io, "---------------------------")
    println(io, "Velocity: ", ke.vel);
    println(io, "Acceleration: ", ke.aec);
    println(io, "Timestamp: ", ke.t);
    println(io, "SpeedType: ", ke.speedtype);
end

"""
Convert!(from::KMH,to::MPH,ke::KinematicEvent)

example:

    a = KinematicEvent(0.0,1.0,DateTime(2019,07,4),KMH())

    Convert!(MPH(),KMH(),a) # throws error
    Convert!(KMH(),MPH(),a) # throws error

"""
# Conversion function from KMH to MPH
function Convert!(from::KMH,to::MPH,ke::KinematicEvent)
    if( ke.speedtype == MPH())
        error("Speed type is already in MPH")
    end
    ke.vel = ke.vel*0.6213711922;
    ke.aec = ke.aec*0.6213711922;
    ke.speedtype = MPH();
end

"""
Convert!(from::MPH,to::KMH,ke::KinematicEvent)

example:

    a = KinematicEvent(0.0,1.0,DateTime(2019,07,4),KMH())

    Convert!(MPH(),KMH(),a) # throws error
    Convert!(KMH(),MPH(),a) # throws error

"""
# Conversion function from MPH to KMH
function Convert!(from::MPH,to::KMH,ke::KinematicEvent)
    if( ke.speedtype == KMH())
        error("Speed type is already in KMH")
    end
    ke.vel = ke.vel*1.6093440000966945;
    ke.aec = ke.aec*1.6093440000966945;
    ke.speedtype = KMH();
end

"""
Kinematic Events Structure

A basic structure for velocity, acceleration and timestamp.

Velocity must be positive!

extension of the single Kinematic Event to incorproate multiple kinematic events

Constructor:

	KinematicEvents(vels::Array{Float64,1},aecs::Array{Float64,1} ,speedtype::SpeedType,tstamps::Array{DateTime,1})

"""
mutable struct KinematicEvents<:TelemetricEvents
    vels::Array{Float64,1}
    aecs::Array{Float64,1}
    speedtype::SpeedType
    tstamps::Array{DateTime,1}
    function KinematicEvents(vels::Array{Float64,1},aecs::Array{Float64,1} ,speedtype::SpeedType,tstamps::Array{DateTime,1})
        if(size(vels,1) != size(aecs,1))
            error("Vels and Aecs must be the same size ")
        end
        # Check velocities
        for vel in vels
            if vel < 0.0
                error("vels cannot be negative, one of the velocities is negative")
            end
        end

        # return KinematicEvents
        kes = new(vels,aecs,speedtype,tstamps);
        return (kes);
    end
end

function Base.show(io::IO, kes::KinematicEvents)
    println(io, "---------------------------")
    println(io, "KinematicEvents object");
    println(io, "---------------------------")
    println(io, "No. Events: ", (kes.vels |> size)[1] );
    println(io, "SpeedType: ", kes.speedtype);
end

@recipe function plot(ev::KinematicEvents)
    vels::Array{Float64,1} = ev.vels;
    aecs::Array{Float64,1} = ev.aecs;
    seriestype := :scatter
    vels, aecs ;
end

"""
Kinematic Map Structure

A basic structure for velocity and acceleration heat map

vel bounds , lower and upper bound for velocity

aec bounds , lower and upper bound acceleration

kernel , BivariateKDE

Constructor:

	KinematicMap(kes::KinematicEvents,vel_bounds::Array{Float64,1}, aec_bounds::Array{Float64,1})

"""
mutable struct KinematicMap <: TelemetricMap
    vel_bounds::Array{Float64,1}
    aec_bounds::Array{Float64,1}
    kernel::BivariateKDE
    events::KinematicEvents
    function KinematicMap(kes::KinematicEvents,vel_bounds::Array{Float64,1}, aec_bounds::Array{Float64,1})
        # Check errors
        if(size(vel_bounds,1) != 2)
            error("vel bounds must have exactly two bounds")
        elseif(size(aec_bounds,1) != 2)
            error("aec bounds must have exactly two bounds")
        end
        # Start of try catch statement, many possible errors
        try
            vels::Array{Float64,1} = kes.vels;
            aecs::Array{Float64,1} = kes.aecs;
            dat::DataFrame = DataFrame(vels = vels, aecs = aecs, tstamps = kes.tstamps);
            dat_filtered::DataFrame = @from i in dat begin
                            @where i.vels > vel_bounds[1] && i.vels < vel_bounds[2]
                            @where i.aecs > aec_bounds[1] && i.aecs < aec_bounds[2]
                            @select {i.vels , i.aecs, i.tstamps}
                            @collect DataFrame
                        end

            vels_filtered::Array{Float64,1} = convert(Array{Float64,1},dat_filtered.vels);
            aecs_filtered::Array{Float64,1} = convert(Array{Float64,1},dat_filtered.aecs);

            kernel::BivariateKDE = kde((vels_filtered,aecs_filtered)) # BivariateKDE

            events = KinematicEvents(vels_filtered,
                                    aecs_filtered,
                                    kes.speedtype,
                                    convert(Array{DateTime,1},dat_filtered.tstamps));

            km::KinematicMap = new(vel_bounds,aec_bounds,kernel,events);
            return (km);
        catch e
            println("ERROR: $e")
            error("KinematicMap construction failed, try using debugger and entering the function")
        end # End of try catch statement
    end # End of Constructor
end # End of mutable struct

function Base.show(io::IO, map::KinematicMap)
    println(io, "---------------------------")
    println(io, "KinematicMap object");
    println(io, "---------------------------")
    println(io, "Velocity from $(map.vel_bounds[1]) to $(map.vel_bounds[2])");
    println(io, "Acceleration from $(map.aec_bounds[1]) to $(map.aec_bounds[2])")
    println(io, "No. Events: ", (map.events.vels |> size)[1] );
    println(io, "SpeedType: ", map.events.speedtype);
end
"""
plots the KinematicMap

# example

    using Rmath

    x = rnorm(100) |> x - > abs.(x)
    y = rnorm(100)

    kes = KinematicEvents(x,y KMH())
    tele_map = KinematicMap(kes,[0.0,1.0],[-1.0,1.0])

    plot(tele_map)
"""
@recipe function plot(km::KinematicMap)
    x_h::Float64 = (km.vel_bounds[2] - km.vel_bounds[1])/100
    y_h::Float64 = (km.aec_bounds[2] - km.aec_bounds[1])/100
    xy_grid = (km.vel_bounds[1]:x_h:km.vel_bounds[2],
                km.aec_bounds[1]:y_h:km.aec_bounds[2])
    # Interpolate
    ik = InterpKDE(km.kernel);

    seriestype := :heatmap
    color --> :plasma
    fill := true

    z = pdf(ik,xy_grid[1],xy_grid[2])'

    xy_grid[1],xy_grid[2],z
end
#=
"""
Overlays plot KinematicMap, only one input so far.

# example

    using Rmath

    x = rnorm(100) |> x - > abs.(x)
    y = rnorm(100)

    kes = KinematicEvents(x,y KMH())
    tele_map = KinematicMap(kes,[0.0,1.0],[-1.0,1.0])

    plot!(tele_map)
"""
function plot!(km::KinematicMap)
    x_h::Float64 = (km.vel_bounds[2] - km.vel_bounds[1])/100
    y_h::Float64 = (km.aec_bounds[2] - km.aec_bounds[1])/100
    xy_grid = (km.vel_bounds[1]:x_h:km.vel_bounds[2],
                km.aec_bounds[1]:y_h:km.aec_bounds[2])
    # Interpolate
    ik = InterpKDE(km.kernel);
    contour!(xy_grid[1],xy_grid[2],pdf(ik,xy_grid[1],xy_grid[2]))
end
=#

"""
DeviationEvents structure

Identify the deviational telemtric events at a power level α

contains:
events, original events
map, KinematicMap
deviants, deviant KinematicEvents
p_kernel, UnivariateKDE used to calculate

"""
mutable struct DeviationEvents <: TelemetricEvents
    α::Float64 # power level alpha, default is usually 0.05
    events::KinematicEvents
    bikernel::BivariateKDE
    deviants::KinematicEvents
    p_kernel::UnivariateKDE
    p_is::Array{Float64,1}
    p_star::Float64
    function DeviationEvents(events::KinematicEvents,bikernel::BivariateKDE,α = 0.05)
        ev::KinematicEvents = events
        ik = InterpKDE(bikernel);
        p_is::Array{Float64,1} = map( (v,a) -> pdf(ik,v,a), ev.vels, ev.aecs);
        try
            p_kernel = kde(p_is,npoints = 2^13);
            # Fix for numerical issues
            p_kernel.density = p_kernel.density ./ (p_kernel.density |> sum)
            # offset range due to weird left bound
            offset_p = abs(p_kernel.x[1])::Float64
	    # check bound and push

        domain_x = convert(Array{Float64,1},p_kernel.x)

	    p_kernel.x[1] < 0 ? (domain_x = (domain_x)) : nothing ;

            # Get devient events by solving for the function
            function GetStar()
                for i in 1:size(p_kernel.density,1)
                    # Solve for this
                    if( (p_kernel.density[1:i] |> sum) > α )
                        p_zero_star = domain_x[(i-1)]
                        return(p_zero_star);
                    end
                end
                return nothing
            end
            p_zero_star = GetStar()

            # get all p_is less than the zero_star
            # offset p_is
            p_is = p_is .+ offset_p
            p_zero_star += offset_p
            p_kernel.x = p_kernel.x .+ offset_p
            checks_p = p_is .< p_zero_star;
            deviants = KinematicEvents(ev.vels[checks_p],ev.aecs[checks_p],ev.speedtype,ev.tstamps[checks_p])
            de = new(α,events,bikernel,deviants,p_kernel,p_is,p_zero_star)
            return(de)
        catch e
            print("ERROR: $e")
            error("\n deviational events constructuion failed, try debugging, most likely KDE")
        end
    end
end

function Base.show(io::IO, dev::DeviationEvents)
    println(io, "---------------------------")
    println(io, "DeviationEvents object");
    println(io, "---------------------------")
    println(io, "α-level: $(dev.α)");
    println(io, "No. Deviant Events: ", (dev.deviants.vels |> size)[1] );
    println(io, "SpeedType: ", dev.deviants.speedtype);
end

"""
Plots deviation events, only one input so far.
"""
@recipe function plot(dev::DeviationEvents)
    vels::Array{Float64,1} = dev.deviants.vels;
    aecs::Array{Float64,1} = dev.deviants.aecs;
    seriestype := :scatter
    markersize --> 2
    vels, aecs
end

#=
function plot!(dev::DeviationEvents)
    vels::Array{Float64,1} = dev.deviants.vels;
    aecs::Array{Float64,1} = dev.deviants.aecs;
    scatter!(vels,aecs,markersize = 2);
end
=#

#=

abstract type TelemetricMap end;
abstract type TelemetricEvent end;
abstract type TelemetricEvents end;
abstract type SpeedType end;

ev = ke
ik = InterpKDE(te.kernel)
p_is = map( (v,a) -> pdf(ik,v,a), ev.vels, ev.aecs);
p_kernel = kde(p_is);
deviants = ev
=#
#=
using Rmath,Plots

a = rnorm(10000) |> x -> abs.(x)
b = rnorm(10000)
ke = KinematicEvents(a,b,KMH(),repeat([DateTime(2019,1)],size(a,1)))
te = KinematicMap(ke,[0.0,1],[-2.0,2.0])

de = DeviationEvents(te.events,te.kernel,0.15)

plot(te)
scatter!(de.d)
=#

#=
# Example with real data
using CSV, Plots

te = CSV.read("/Users/pocucan/Downloads/2011-08-25.csv")
te = CSV.read("/Users/nik/Downloads/2011-08-25.csv")
vels = te.speed_mph |> x -> convert(Array{Float64,1},x)
aecs = te.accel_meters_ps |> x -> convert(Array{Float64,1},x)

ke = KinematicEvents(vels,aecs,MPH(),repeat([DateTime(2019,1)],size(vels,1)))
#mp = KinematicMap(ke,[30.0,50.0],[-5.0,5.0])
mp = KinematicMap(ke,[40.0,90.0],[-2.0,2.0])

dev = DeviationEvents(mp.events,mp.kernel,0.05)

dev_1 = DeviationEvents(mp.events,mp.kernel,0.01)

=#


text:
using MLBase
# using LIBSVM
# using Mocha

include("core.jl")

function svm_simple()
    # Xfull, yfull, mask = load_image_dataset(FaceDatasets.CKMaxDataset; da)
    Xfull = load_images(FaceDatasets.CKMaxDataset, DATA_DIR_CK)
    X = Xfull[:, yfull .!= -1]
    y = yfull[yfull .!= -1]
    # train_idxs = collect(StratifiedKfold(y, 10))[1]
    # test_idxs = setdiff(1:length(y), train_idxs)
    train_idxs = collect(1:200)
    test_idxs = collect(201:length(y))
    @time model = svmtrain(y[train_idxs], X[:, train_idxs])
    @time yhat = svmpredict(model, X[:, test_idxs])

    @printf "%.2f%%\n" mean(yhat .== y[test_idxs])*100
end

function svm_rbm()
    Xfull, yfull, mask = load_image_dataset(:ck_max)
    X = Xfull[:, yfull .!= -1]
    y = yfull[yfull .!= -1]
    # train_idxs = collect(StratifiedKfold(y, 10))[1]
    # test_idxs = setdiff(1:length(y), train_idxs)
    train_idxs = collect(1:200)
    test_idxs = collect(201:length(y))

    rbm = RBM(Degenerate, Bernoulli, 1, 1)
    h5open(joinpath(MODEL_DIR, "rbm_512.h5")) do h5
        load_params(h5, rbm, "model")
    end
    Xt = Boltzmann.transform(rbm, X)

    @time model = svmtrain(y[train_idxs], Xt[:, train_idxs])
    @time yhat = svmpredict(model, Xt[:, test_idxs])

    @printf "%.2f%%\n" mean(yhat .== y[test_idxs])*100
end

function mocha_simple()
    Xfull, yfull, mask = load_image_dataset(:ck_max)
    X = Xfull[:, yfull .!= -1]
    y = yfull[yfull .!= -1]
    # train_idxs = collect(StratifiedKfold(y, 10))[1]
    # test_idxs = setdiff(1:length(y), train_idxs)
    train_idxs = collect(1:200)
    test_idxs = collect(201:length(y))

    train_data_path = joinpath(MODEL_DIR, "train_data_simple.h5")
    test_data_path = joinpath(MODEL_DIR, "test_data_simple.h5")
    h5open(train_data_path, "w") do h5
        write(h5, "data", X[:, train_idxs])
        write(h5, "label", y[train_idxs])
    end
    h5open(test_data_path, "w") do h5
        write(h5, "data", X[:, test_idxs])
        write(h5, "label", y[test_idxs])
    end

    data  = HDF5DataLayer(name="train-data", source="train-data.txt",
                          batch_size=64)
    ## conv  = ConvolutionLayer(name="conv1",n_filter=20,kernel=(5,5),bottoms=[:data],tops=[:conv])
    ## pool  = PoolingLayer(name="pool1",kernel=(2,2),stride=(2,2),bottoms=[:conv],tops=[:pool])
    ## conv2 = ConvolutionLayer(name="conv2",n_filter=50,kernel=(5,5),bottoms=[:pool],tops=[:conv2])
    ## pool2 = PoolingLayer(name="pool2",kernel=(2,2),stride=(2,2),bottoms=[:conv2],tops=[:pool2])
    fc1   = InnerProductLayer(name="ip1",output_dim=500,neuron=Neurons.ReLU(),bottoms=[:data],
                              tops=[:ip1])
    ## fc2   = InnerProductLayer(name="ip2",output_dim=10,bottoms=[:ip1],tops=[:ip2])
    loss  = SoftmaxLossLayer(name="loss",bottoms=[:ip1,:label])

    backend = DefaultBackend()
    init(backend)

    net = Net("MNIST-train", backend, [data, fc1, loss])

    exp_dir = "snapshots"
    solver_method = SGD()
    params = make_solver_parameters(solver_method, max_iter=10000, regu_coef=0.0005,
                                    mom_policy=MomPolicy.Fixed(0.9),
                                    lr_policy=LRPolicy.Inv(0.01, 0.0001, 0.75),
                                    load_from=exp_dir)
    solver = Solver(solver_method, params)

    setup_coffee_lounge(solver, save_into="$exp_dir/statistics.jld", every_n_iter=1000)

    # report training progress every 100 iterations
    add_coffee_break(solver, TrainingSummary(), every_n_iter=100)

    # save snapshots every 5000 iterations
    add_coffee_break(solver, Snapshot(exp_dir), every_n_iter=5000)

    # show performance on test data every 1000 iterations
    data_test = HDF5DataLayer(name="test-data",source="test-data.txt",batch_size=100)
    accuracy = AccuracyLayer(name="test-accuracy",bottoms=[:ip1, :label])
    test_net = Net("MNIST-test", backend, [data_test, fc1, accuracy])
    add_coffee_break(solver, ValidationPerformance(test_net), every_n_iter=1000)

    solve(solver, net)

    destroy(net)
    destroy(test_net)

end


text:
#  extension of N {\displaystyle N} N by H {\displaystyle H} H).
#  If any of these statements holds (and hence all of them hold, by their equivalence), we say G is the semidirect product of N and H, written
#      G = N ⋊ H {\displaystyle G=N\rtimes H} {\displaystyle G=N\rtimes H} or G = H ⋉ N , {\displaystyle G=H\ltimes N,} {\displaystyle G=H\ltimes N,}
#  or that G splits over N; one also says that G is a semidirect product of H acti

export SymmorphicSymmetry
export SymmorphicIrrepComponent
export ⋊, ⋉
export symmetry_product
export group, group_order, group_multiplication_table,
       element, elements, element_name, element_names

export fractional_momentum
export generator_indices, generator_elements

"""
    SymmorphicSymmetry{S1, S2, E}

Symmorphic symmetry of symmetry `S1` ⋊ `S2`, whose elements are of type `E`.
For symmorphic space group, the normal symmetry `S1` is translation symmetry, and `S2` is
point symmetry. `S1` and `S2` can either be `AbstractSymmetry` or `SymmetryEmbedding`.

# Fields
* `normal::S1`: normal subgroup
* `rest::S2`
* `elements::Array{E}`
* `element_names::Array{String}`
"""
struct SymmorphicSymmetry{
    S1<:SymmetryOrEmbedding,
    S2<:SymmetryOrEmbedding,
    E<:AbstractSpaceSymmetryOperation
}<:AbstractSymmetry{E}

    normal::S1   # e.g.) Translation Symmetry
    rest::S2     # e.g.) Point Symmetry
    elements::Array{E}
    element_names::Array{String}

    @doc """
        SymmorphicSymmetry(normal, rest)

    Construct a symmorphic symmetry `normal ⋊ rest`.
    """
    function SymmorphicSymmetry(normal::S1, rest::S2) where {S1<:SymmetryOrEmbedding, S2<:SymmetryOrEmbedding}
        if !iscompatible(normal, rest)
            throw(ArgumentError("symmetries $normal and $rest are not compatible"))
        end
        E = promote_type(eltype(normal), eltype(rest))
        if !isconcretetype(E)
            throw(ArgumentError("element type $E is not a concrete type"))
        end
        for g in generator_elements(rest), n in generator_elements(normal)
            n2 = g * n * inv(g)
            if n2 ∉ normal
                throw(ArgumentError("first argument must be normal"))
            end
        end
        elems = E[x*y for y in elements(normal), x in elements(rest)]
        elem_names = String["{ $y | $x }"  for y in element_names(normal), x in element_names(rest)]
        @assert allunique(elems)

        return new{S1, S2, E}(normal, rest, elems, elem_names)
    end
end

"""
    ⋊(normal::AbstractSymmetry, rest::AbstractSymmetry)

Return symmorphic symmetry `normal ⋊ rest`.
"""
function ⋊(normal::AbstractSymmetry, rest::AbstractSymmetry)
    return SymmorphicSymmetry(normal, rest)
end

"""
    ⋉(rest::AbstractSymmetry, normal::AbstractSymmetry)

Return symmorphic symmetry `normal ⋊ rest`.
"""
function ⋉(rest::AbstractSymmetry, normal::AbstractSymmetry)
    return SymmorphicSymmetry(normal, rest)
end

function Base.:(==)(lhs::SymmorphicSymmetry{S1, S2, E}, rhs::SymmorphicSymmetry{S1, S2, E}) where {S1, S2, E}
    return lhs.normal == rhs.normal && lhs.rest == rhs.rest
end

"""
    symmetry_product(sym)

# Arguments
* `sym::SymmorphicSymmetry{TranslationSymmetry,PointSymmetry,S}) where {S<:SpaceOperation{<:Integer}}`
"""
function symmetry_product(
    sym::SymmorphicSymmetry{TranslationSymmetry,PointSymmetry,S},
) where {S<:SpaceOperation{<:Integer}}
    function product(lhs::SpaceOperation{<:Integer}, rhs::SpaceOperation{<:Integer})
        foo = lhs * rhs
        return S(foo.matrix, sym.normal.orthocube.wrap(foo.displacement)[2])
    end
    return product
end

"""
    group(sym::SymmorphicSymmetry)

Return `FiniteGroup` of the elements.
"""
function group(sym::SymmorphicSymmetry)
    return FiniteGroup(group_multiplication_table(vec(sym.elements), symmetry_product(sym)))
end

"""
    group_order(sym::SymmorphicSymmetry)

Group order of `sym`, which is the product of the subgroup and the normal subgroup.
"""
function group_order(sym::SymmorphicSymmetry)
    return group_order(sym.normal) * group_order(sym.rest)
end

"""
    group_multiplication_table(sym::SymmorphicSymmetry)

Return group multiplication table of `sym`.
"""
function group_multiplication_table(sym::SymmorphicSymmetry)
    return group_multiplication_table(vec(sym.elements), symmetry_product(sym))
end

Base.eltype(sym::SymmorphicSymmetry{S1, S2, E}) where {S1, S2, E} = E
Base.valtype(sym::SymmorphicSymmetry{S1, S2, E}) where {S1, S2, E} = E

"""
    element(sym::SymmorphicSymmetry, g...)

Return the `i`th element of the symmorphic symmetry.
"""
element(sym::SymmorphicSymmetry, g...) = sym.elements[g...]

"""
    elements(sym::SymmorphicSymmetry)

Return the elements of the symmorphic symmetry.
"""
elements(sym::SymmorphicSymmetry) = sym.elements

"""
    element_names(sym::SymmorphicSymmetry)

Return the element names of the symmorphic symmetry.
"""
element_names(sym::SymmorphicSymmetry) = sym.element_names

"""
    element_name(sym::SymmorphicSymmetry, g...)

Return the name of the `g`th element of the symmorphic symmetry.
"""
element_name(sym::SymmorphicSymmetry, g...) = sym.element_names[g...]

# character_table(sym::SymmorphicSymmetry) = sym.character_table
# irreps(sym::SymmorphicSymmetry) = sym.irreps
# irrep(sym::SymmorphicSymmetry, idx) = sym.irreps[idx]
# num_irreps(sym::SymmorphicSymmetry) = length(sym.irreps)
# irrep_dimension(sym::SymmorphicSymmetry, idx::Integer) = 1 # size(first(sym.irreps[idx]), 1)

"""
    fractional_momentum(sym, g)

Return the `g`th fractional momentum of the normal (translation) symmetry.

# Arguments
* `sym::SymmorphicSymmetry{<:TranslationSymmetry, S2, E}`
"""
function fractional_momentum(
    sym::SymmorphicSymmetry{<:TranslationSymmetry, S2, E},
    args...
) where {S2, E}
    return fractional_momentum(sym.normal, args...)
end

"""
    generator_indices(sym::SymmorphicSymmetry)

Return the indices of the generators, which is a union of the generators of normal and rest.
"""
function generator_indices(sym::SymmorphicSymmetry)
    gn, gr = sym.normal.generators, sym.rest.generators
    return vcat(
        [CartesianIndex(x..., 1) for x in gn],
        [CartesianIndex(1, x...) for x in gr],
    )
end

"""
    generator_elements(sym::SymmorphicSymmetry)

Return the generating elements
"""
function generator_elements(sym::SymmorphicSymmetry)
    g = generator_indices(sym)
    return element(sym, g)
end


text:
"""
DuckDB type
"""
mutable struct LogicalType
    handle::duckdb_logical_type

    function LogicalType(type::DUCKDB_TYPE)
        handle = duckdb_create_logical_type(type)
        result = new(handle)
        finalizer(_destroy_type, result)
        return result
    end
    function LogicalType(handle::duckdb_logical_type)
        result = new(handle)
        finalizer(_destroy_type, result)
        return result
    end
end

function _destroy_type(type::LogicalType)
    if type.handle != C_NULL
        duckdb_destroy_logical_type(type.handle)
    end
    type.handle = C_NULL
    return
end

CreateLogicalType(::Type{T}) where {T <: String} = DuckDB.LogicalType(DuckDB.DUCKDB_TYPE_VARCHAR)
CreateLogicalType(::Type{T}) where {T <: Bool} = DuckDB.LogicalType(DuckDB.DUCKDB_TYPE_BOOLEAN)
CreateLogicalType(::Type{T}) where {T <: Int8} = DuckDB.LogicalType(DuckDB.DUCKDB_TYPE_TINYINT)
CreateLogicalType(::Type{T}) where {T <: Int16} = DuckDB.LogicalType(DuckDB.DUCKDB_TYPE_SMALLINT)
CreateLogicalType(::Type{T}) where {T <: Int32} = DuckDB.LogicalType(DuckDB.DUCKDB_TYPE_INTEGER)
CreateLogicalType(::Type{T}) where {T <: Int64} = DuckDB.LogicalType(DuckDB.DUCKDB_TYPE_BIGINT)
CreateLogicalType(::Type{T}) where {T <: UInt8} = DuckDB.LogicalType(DuckDB.DUCKDB_TYPE_UTINYINT)
CreateLogicalType(::Type{T}) where {T <: UInt16} = DuckDB.LogicalType(DuckDB.DUCKDB_TYPE_USMALLINT)
CreateLogicalType(::Type{T}) where {T <: UInt32} = DuckDB.LogicalType(DuckDB.DUCKDB_TYPE_UINTEGER)
CreateLogicalType(::Type{T}) where {T <: UInt64} = DuckDB.LogicalType(DuckDB.DUCKDB_TYPE_UBIGINT)
CreateLogicalType(::Type{T}) where {T <: Float32} = DuckDB.LogicalType(DuckDB.DUCKDB_TYPE_FLOAT)
CreateLogicalType(::Type{T}) where {T <: Float64} = DuckDB.LogicalType(DuckDB.DUCKDB_TYPE_DOUBLE)
CreateLogicalType(::Type{T}) where {T <: Date} = DuckDB.LogicalType(DuckDB.DUCKDB_TYPE_DATE)
CreateLogicalType(::Type{T}) where {T <: Time} = DuckDB.LogicalType(DuckDB.DUCKDB_TYPE_TIME)
CreateLogicalType(::Type{T}) where {T <: DateTime} = DuckDB.LogicalType(DuckDB.DUCKDB_TYPE_TIMESTAMP)
CreateLogicalType(::Type{T}) where {T <: AbstractString} = DuckDB.LogicalType(DuckDB.DUCKDB_TYPE_VARCHAR)

function CreateLogicalType(::Type{T}) where {T}
    throw(NotImplementedException("Unsupported type for CreateLogicalType"))
end

function GetTypeId(type::LogicalType)
    return duckdb_get_type_id(type.handle)
end

function GetInternalTypeId(type::LogicalType)
    type_id = GetTypeId(type)
    if type_id == DUCKDB_TYPE_DECIMAL
        type_id = duckdb_decimal_internal_type(type.handle)
    elseif type_id == DUCKDB_TYPE_ENUM
        type_id = duckdb_enum_internal_type(type.handle)
    end
    return type_id
end

function GetDecimalScale(type::LogicalType)
    return duckdb_decimal_scale(type.handle)
end

function GetEnumDictionary(type::LogicalType)
    dict::Vector{String} = Vector{String}()
    dict_size = duckdb_enum_dictionary_size(type.handle)
    for i in 1:dict_size
        val = duckdb_enum_dictionary_value(type.handle, i)
        str_val = String(unsafe_string(val))
        push!(dict, str_val)
        duckdb_free(val)
    end
    return dict
end

function GetListChildType(type::LogicalType)
    return LogicalType(duckdb_list_type_child_type(type.handle))
end

function GetStructChildCount(type::LogicalType)
    return duckdb_struct_type_child_count(type.handle)
end

function GetStructChildName(type::LogicalType, index::UInt64)
    val = duckdb_struct_type_child_name(type.handle, index)
    result = unsafe_string(val)
    duckdb_free(val)
    return result
end

function GetStructChildType(type::LogicalType, index::UInt64)
    return LogicalType(duckdb_struct_type_child_type(type.handle, index))
end


text:
using JSON

Endpoint("/examples/requests") do request::HTTP.Request
    readstring(joinpath(dirname(@__FILE__),"requests.html"))
end

Endpoint("/examples/requests/echo") do request::HTTP.Request
    if request.method == "GET"
        params = HTTP.queryparams(HTTP.URI(request.target).query)
        println("Body: $(params)")
        response = JSON.json(params)
    elseif request.method == "POST"
        data = String(request.body)
        println("Parameters: $(data)")
        response = JSON.json(Dict(:data => data))
    end
    response
end

# Endpoint("/api/math/rand") do request::HTTP.Request
#     uri = URI(request.resource)
#     if request.method == "GET"
#         param = query_params(uri)
#         if isempty(param)
#             response = Dict(:method => "GET", :value => rand())
#         else
#             n = min(1000,parse(Int,param["n"]))
#             response = Dict(:method => "GET", :value => rand(n))
#         end
#     end
#     if request.method == "POST"
#         param = query_params(uri)
#         response = Dict(:method => "POST", :value => rand())
#     end
#     JSON.json(response)
# end


text:
# Test of companion linearization

using NonlinearEigenproblemsTest
using NonlinearEigenproblems
using Test,Random
using LinearAlgebra

@bench @testset "Companion Linearization" begin

    #####################
    # Dense matrix test #
    #####################
    Random.seed!(0)
    pep = nep_gallery("pep0");
    deg = size(get_fv(pep),1) -1;
    n = size(pep,1);
    tolerance = 1e-14;

    λa,xa = newton(pep, maxit=30, logger=0, tol = tolerance,v=ones(size(pep,1)));
    Dc,Vc = polyeig(pep, DefaultEigSolver);

    ind = 1;
    for i = 1:(deg*n)
        if(abs(λa-Dc[i])/abs(λa) < tolerance*10)
            ind = i;
            break
        end
    end
    @test (abs(λa-Dc[ind])/abs(λa) < tolerance*10)

    ######################
    # Sparse matrix test #
    ######################
    pep = nep_gallery("pep0_sparse", 200, 0.03);
    deg = size(get_fv(pep),1) -1;
    n = size(pep,1);
    tolerance = 1e-14;

    λa,xa =newton(pep, maxit=30, λ=-0.75, v=ones(n), logger=0, tol = tolerance);
    Dc,Vc = polyeig(pep,DefaultEigSolver);

    ind = 1;
    nn = min(deg*n, length(Dc))
    for i = 1:nn
        if(abs(λa-Dc[i])/abs(λa) < tolerance*50)
            ind = i;
            break
        end
    end
    @test (abs(λa-Dc[ind])/abs(λa) < tolerance*50)

    ########################
    # BigFloat matrix test #
    ########################
    A0 = (Array{BigFloat,2}([1 2; 3 4]));
    A1 = (Array{BigFloat,2}([1 44; 3 4]));
    A2 = (Array{BigFloat,2}([1 44; -3 100]));
    pep = PEP([A0,A1,A2])

    # Power method for testing (since eig does not work for BigFloats)
    E,A = companion(pep);
    z=ones(BigFloat,size(A,1));
    local evp::BigFloat
    local evp_old::BigFloat
    local d::BigFloat
    tolerance = big"1e-50"
    d=Inf; evp=Inf
    k=0;
    while abs(d) > tolerance
        k=k+1
        z=z/norm(z);
        z2=E\A*z
        evp_old=evp
        evp=dot(z,z2)
        d=evp-evp_old
        @debug "Iteration $k, λ=$(Float64(evp)), Δ = $(Float64(d))"
        z=z2
    end
    @info "Solving same problem with resinv"
    λ,v = resinv(BigFloat, pep, λ = (BigFloat(evp)+0.1), v = z[1:size(pep,1)], logger=0, tol = tolerance, errmeasure=ResidualErrmeasure(pep))

    @test (abs(λ-evp)/abs(λ) < tolerance*10)

end


text:
using DataStructures

const out_dir = "../../docs/src/diagnostics/poincare_2nd/guiding_center_4d/"

# create list of figures to put in list
figures = (
    "_poincare_2nd_q",
    "_poincare_2nd_p",
    "_poincare_2nd_l"
)

# suffix of figure files
fig_suff = ".png"

# collect simulation parameters
function get_params(Δt, ntime, problem)
    params = OrderedDict()
    params["Δt"] = Δt
    params["nt"] = ntime
    params["nloop"] = nloop
    params["nsave"] = nsave
    params["nplot"] = nplot
    params["GeometricIntegrators.jl"] = package_version("GeometricIntegrators")
    params["ChargedParticleDynamics.jl"] = package_version("ChargedParticleDynamics")
    params["MagneticEquilibria.jl"] = package_version("MagneticEquilibria")
    params
end

function run_tableau_list(tableaus, pinv, prefix, out_dir, fig_dir, nplot)
    # create output directory
    if !isdir(out_dir * "/" * fig_dir)
        mkpath(out_dir * "/" * fig_dir)
    end

    for tab in tableaus
        integrator = tab[1]
        tableau    = tab[2]
        runid      = tab[3]
        filename   = out_dir * "/" * fig_dir * "/" * prefix*runid*"_nx"*string(pinv.nx)*"_ny"*string(pinv.ny)*"_nt"*string(pinv.ntime)

        sim = Simulation(pinv.equ, integrator, tableau, pinv.Δt, runid, filename * "_solution.h5", pinv.ntime, pinv.nsave)
        sol = run!(sim)

        I, J, K, L = evaluate_poincare_invariant(pinv, sol)

        write_to_hdf5(pinv, sol, filename * "_poincare_2nd.h5")

        info("Plotting results...")

        plot_integral_error(sol.t.t, I, filename * "_poincare_2nd_q" * fig_suff;
                            plot_title=L"$\Delta \int_S \omega_{ij} (q) \, dq^i \wedge dq^j$")

        if isdefined(sol, :p)
            plot_integral_error(sol.t.t, J, filename * "_poincare_2nd_p" * fig_suff;
                                plot_title=L"$\Delta \int_S dp_i \, dq^i$")
        end

        if isdefined(sol, :λ)
            plot_integral_error(sol.t.t, L, filename * "_poincare_2nd_l" * fig_suff;
                                plot_title=L"$\Delta \int_S ( \omega_{ij} (q) \, dq^i \wedge dq^j - \Delta t^2 \, \omega_{ij} (q) \, d\lambda^i \wedge d\lambda^j$")
        end

        plot_poincare_surface(sol, nplot, filename * "_area" * fig_suff; xmin=XMIN, xmax=XMAX, ymin=YMIN, ymax=YMAX)

    end
end

run_tableau_list(tableau_list, pinv, prefix, out_dir, fig_dir, nplot)

# generate HTML file


text:
using Wflow

tomlpath = joinpath(@__DIR__, "sbm_simple.toml")
Wflow.run(tomlpath)

#=
# test whether restarted runs get the same results as continuus ones, i.e. state is captured
tomlpath = joinpath(@__DIR__, "sbm_config.toml")
config = Wflow.Config(tomlpath)

# January at once, cold start
config.starttime = DateTime("2000-01-01T00:00:00")
config.endtime   = DateTime("2000-02-01T00:00:00")
config.model.reinit = true  # cold start
# note that this needs to be relative to the tomlpath
config.state.path_output = joinpath(dirname(tomlpath), "data/state-test/outstates-moselle-january.nc")
config.output.path = joinpath(dirname(tomlpath), "data/state-test/output-moselle-january.nc")
model = Wflow.initialize_sbm_model(config)
Wflow.run(model)

# first half of January, cold start
config.starttime = DateTime("2000-01-01T00:00:00")
config.endtime   = DateTime("2000-01-15T00:00:00")
config.model.reinit = true  # cold start
config.state.path_output = joinpath(dirname(tomlpath), "data/state-test/outstates-moselle-january-1of2.nc")
config.output.path = joinpath(dirname(tomlpath), "data/state-test/output-moselle-january-1of2.nc")
model = Wflow.initialize_sbm_model(config)
Wflow.run(model)

# second half of January, warm start
config.starttime = DateTime("2000-01-15T00:00:00")
config.endtime   = DateTime("2000-02-01T00:00:00")
config.model.reinit = false  # cold start
config.state.path_input = joinpath(dirname(tomlpath), "data/state-test/outstates-moselle-january-1of2.nc")
config.state.path_output = joinpath(dirname(tomlpath), "data/state-test/outstates-moselle-january-2of2.nc")
config.output.path = joinpath(dirname(tomlpath), "data/state-test/output-moselle-january-2of2.nc")
model = Wflow.initialize_sbm_model(config)
Wflow.run(model)
=#