text:
# The following methods are uses in VoronoiFVM.
# Do we need a more systematic approach here ?
"""
$(SIGNATURES)

Number of edges in grid.
"""
num_edges(grid::ExtendableGrid)=haskey(grid,EdgeNodes) ?  num_sources(grid[EdgeNodes]) : 0

#########################################
# VERTEX0D <: AbstractElementGeometry0D #       [1]
#########################################

const _refcoords_for_geometry_Vertex0D = reshape([0],1,1)
const _local_celledgenodes_Vertex0D = reshape([1],1,1)

"""
$(SIGNATURES)

Coordinates of reference geometry of 0D vertex
"""
refcoords_for_geometry(::Type{<:AbstractElementGeometry0D}) = _refcoords_for_geometry_Vertex0D

"""
$(SIGNATURES)

Number of nodes of 0D vertex
"""
num_nodes(::Type{<:AbstractElementGeometry0D}) = 1

"""
$(SIGNATURES)

Number of faces of 0D vertex
"""
num_faces(::Type{<:AbstractElementGeometry0D}) = 0

"""
$(SIGNATURES)

Number of edges of 0D vertex
"""
num_edges(::Type{<:AbstractElementGeometry0D}) = 0

"""
$(SIGNATURES)

Cell-edge node numbering for 1D edge
"""
local_celledgenodes(::Type{Vertex0D}) = _local_celledgenodes_Vertex0D

#######################################
# EDGE1D <: AbstractElementGeometry1D #     [1]-----[2]        [1] = 0, [2] = 1
#######################################

const _refcoords_for_geometry_Edge1D = reshape([0; 1]',1,2)
const _local_cellfacenodes_Edge1D = reshape([1; 2],1,2)
const _local_celledgenodes_Edge1D = reshape([1; 2],2,1)

"""
$(SIGNATURES)

Coordinates of reference geometry of 1D edge
"""
refcoords_for_geometry(::Type{<:AbstractElementGeometry1D}) = _refcoords_for_geometry_Edge1D

"""
$(SIGNATURES)

Number of nodes for 1D edge
"""
num_nodes(::Type{<:AbstractElementGeometry1D}) = 2

"""
$(SIGNATURES)

Number of faces for 1D edge
"""
num_faces(::Type{<:AbstractElementGeometry1D}) = 2

"""
$(SIGNATURES)

Number of edges for 1D edge
"""
num_edges(::Type{<:AbstractElementGeometry1D}) = 0

"""
$(SIGNATURES)

Cell-face node numbering for 1D edge
"""
local_cellfacenodes(::Type{<:AbstractElementGeometry1D}) = _local_cellfacenodes_Edge1D

"""
$(SIGNATURES)

Cell-edge node numbering for 1D edge
"""
local_celledgenodes(::Type{<:AbstractElementGeometry1D}) = _local_celledgenodes_Edge1D

"""
$(SIGNATURES)

Number of edges of 1D edge
"""
num_edges(::Type{Edge1D})=1

"""
$(SIGNATURES)

Geometries of faces of 1D edge
"""
facetype_of_cellface(::Type{<:AbstractElementGeometry1D}, k) = Vertex0D

xrefFACE2xrefCELL(::Type{<:AbstractElementGeometry1D}) = [ [(xref4FACE) -> [1]],
                                                           [(xref4FACE) -> [1]] ]

xrefFACE2xrefOFACE(::Type{<:AbstractElementGeometry1D}) = [(xref4FACE) -> xref4FACE, (xref4FACE) -> 1 .- xref4FACE]

#                   [3]
#                    | \
##############       |   \                    [1] = (0,0)
# Triangle2D #       |     \                  [2] = (1,0)
##############       |       \                [3] = (0,1)
#                    |         \
#                   [1]--------[2]

const _refcoords_for_geometry_Triangle2D = [0 0; 1 0; 0 1]'
const _local_cellfacenodes_Triangle2D = [1 2; 2 3; 3 1]' # local edgenodes are the same

"""
$(SIGNATURES)

Coordinates of reference geometry of 2D triangle
"""
refcoords_for_geometry(::Type{<:Triangle2D}) = _refcoords_for_geometry_Triangle2D

"""
$(SIGNATURES)

Number of nodes in 2D triangle
"""
num_nodes(::Type{<:Triangle2D}) = 3

"""
$(SIGNATURES)

Number of faces in 2D triangle
"""
num_faces(::Type{<:Triangle2D}) = 3

"""
$(SIGNATURES)

Number of edges in 2D triangle
"""
num_edges(::Type{<:Triangle2D}) = 3

"""
$(SIGNATURES)

Cell-face node numbering for 2D triangle
"""
local_cellfacenodes(::Type{<:Triangle2D}) = _local_cellfacenodes_Triangle2D

"""
$(SIGNATURES)

Cell-edge node numbering for 2D triangle
"""
local_celledgenodes(::Type{<:Triangle2D}) = _local_cellfacenodes_Triangle2D

"""
$(SIGNATURES)

Geometries of faces of 2D triangle
"""
facetype_of_cellface(::Type{<:Triangle2D}, k) = Edge1D

# maps of reference coords on cell face to reference coords in cell
xrefFACE2xrefCELL(::Type{<:Triangle2D}) = [ (xref4FACE) -> [xref4FACE[1],0],
                                            (xref4FACE) -> [1-xref4FACE[1],xref4FACE[1]],
                                            (xref4FACE) -> [0,1-xref4FACE[1]]
                                            ]

# maps of reference coords on face to reference coords in face with other orientation
xrefFACE2xrefOFACE(::Type{<:Triangle2D}) = [(xref4FACE) -> xref4FACE,                                    # orientation 1 = [1,2,3]
                                            (xref4FACE) -> [xref4FACE[1],1-xref4FACE[1]-xref4FACE[2]],   # orientation 2 = [1,3,2]
                                            (xref4FACE) -> [1-xref4FACE[1]-xref4FACE[2],xref4FACE[2]],   # orientation 3 = [3,2,1]
                                            (xref4FACE) -> [xref4FACE[2],xref4FACE[1]]                   # orientation 4 = [2,1,3]
                                            ]

#                        [4]--------[3]
#                         |          |             [1] = (0,0)
###################       |          |             [2] = (1,0)
# Quadrilateral2D #       |          |             [3] = (1,1)
###################       |          |             [4] = (0,1)
#                        [1]--------[2]

const _refcoords_for_geometry_Quadrilateral2D = [0 0; 1 0; 1 1; 0 1]'
const _local_cellfacenodes_Quadrilateral2D = [1 2; 2 3; 3 4; 4 1]' # local edgenodes are the same

"""
$(SIGNATURES)

Coordinates of reference geometry of 2D quadrilateral
"""
refcoords_for_geometry(::Type{<:Quadrilateral2D}) = _refcoords_for_geometry_Quadrilateral2D

"""
$(SIGNATURES)

Number of nodes in 2D quadrilateral
"""
num_nodes(::Type{<:Quadrilateral2D}) = 4

"""
$(SIGNATURES)

Number of faces in 2D quadrilateral
"""
num_faces(::Type{<:Quadrilateral2D}) = 4

"""
$(SIGNATURES)

Number of edges in 2D quadrilateral
"""
num_edges(::Type{<:Quadrilateral2D}) = 4

"""
$(SIGNATURES)

Cell-face node numbering for 2D quadrilateral
"""
local_cellfacenodes(::Type{<:Quadrilateral2D}) = _local_cellfacenodes_Quadrilateral2D

"""
$(SIGNATURES)

Cell-edge node numbering for 2D quadrilateral
"""
local_celledgenodes(::Type{<:Quadrilateral2D}) = _local_cellfacenodes_Quadrilateral2D

"""
$(SIGNATURES)

Geometries of faces of 2D quadrilateral
"""
facetype_of_cellface(::Type{<:Quadrilateral2D}, k) = Edge1D

# maps of reference coords on cell face to reference coords in cell
xrefFACE2xrefCELL(::Type{<:Quadrilateral2D}) = [ (xref4FACE) -> [xref4FACE[1],0],
                                                 (xref4FACE) -> [1,xref4FACE[1]],
                                                 (xref4FACE) -> [1-xref4FACE[1],1],
                                                 (xref4FACE) -> [0,1-xref4FACE[1]]
                                                 ]

#                      [4]
#                       |\\
#################       | \ \                    [1] = (0,0,0)
# Tetrahedron3D #       |  \  \                  [2] = (1,0,0)
#################       |   \   \                [3] = (0,1,0)
#                       | _-[3]-_ \              [4] = (0,0,1)
#                      [1]--------[2]

const _refcoords_for_geometry_Tetrahedron3D = [0 0 0; 1 0 0; 0 1 0; 0 0 1]'
const _local_cellfacenodes_Tetrahedron3D = [1 3 2; 1 2 4; 2 3 4; 1 4 3]'
const _local_celledgenodes_Tetrahedron3D = [1 2; 1 3; 1 4; 2 3; 2 4; 3 4]'

"""
$(SIGNATURES)

Coordinates of reference geometry of 3D tetrahedron
"""
refcoords_for_geometry(::Type{<:Tetrahedron3D}) = _refcoords_for_geometry_Tetrahedron3D

"""
$(SIGNATURES)

Number of nodes in 3D tetrahedron
"""
num_nodes(::Type{<:Tetrahedron3D}) = 4

"""
$(SIGNATURES)

Number of faces in 3D tetrahedron
"""
num_faces(::Type{<:Tetrahedron3D}) = 4

"""
$(SIGNATURES)

Number of edges in 3D tetrahedron
"""
num_edges(::Type{<:Tetrahedron3D}) = 6

"""
$(SIGNATURES)

Cell-face node numbering for 3D tetrahedron
"""
local_cellfacenodes(::Type{<:Tetrahedron3D}) = _local_cellfacenodes_Tetrahedron3D

"""
$(SIGNATURES)

Geometries of faces of 3D tetrahedron
"""
facetype_of_cellface(::Type{<:Tetrahedron3D}, k) = Triangle2D

"""
$(SIGNATURES)

Cell-edge node numbering for 3D tetrahedron
"""
local_celledgenodes(::Type{<:Tetrahedron3D}) = _local_celledgenodes_Tetrahedron3D

# maps of reference coords on cell face to reference coords in cell
xrefFACE2xrefCELL(::Type{<:Tetrahedron3D}) = [ (xref4FACE) -> [xref4FACE[2],xref4FACE[1],0],
                                               (xref4FACE) -> [xref4FACE[1],0,xref4FACE[2]],
                                               (xref4FACE) -> [1-xref4FACE[1]-xref4FACE[2],xref4FACE[1],xref4FACE[2]],
                                               (xref4FACE) -> [0,xref4FACE[2],xref4FACE[1]]
                                                ]

#
#                         [8]--------[7]
#                        / |        / |          [1] = (0,0,0)
#                     [5]--------[6]  |          [2] = (1,0,0)
################       |   |      |   |          [3] = (1,1,0)
# Hexahedron3D #       |   |      |   |          [4] = (0,1,0)
################       |   |      |   |          [5] = (0,0,1)
#                      |  [4]-----|--[3]         [6] = (1,0,1)
#                      | /        | /            [7] = (1,1,1)
#                     [1]--------[2]             [8] = (0,1,1)

const _refcoords_for_geometry_Hexahedron3D = [0 0 0; 1 0 0; 1 1 0; 0 1 0; 0 0 1; 1 0 1; 1 1 1; 0 1 1]'
const _local_cellfacenodes_Hexahedron3D = [4 3 2 1; 1 2 6 5; 2 3 7 6; 3 4 8 7; 4 1 5 8; 5 6 7 8]'
const _local_celledgenodes_Hexahedron3D = [1 2; 2 3; 3 4; 4 1; 1 5; 2 6; 3 7; 4 8; 5 6; 6 7; 7 8; 8 5]'

"""
$(SIGNATURES)

Coordinates of reference geometry of 3D hexahedron
"""
refcoords_for_geometry(::Type{<:Hexahedron3D}) = _refcoords_for_geometry_Hexahedron3D

"""
$(SIGNATURES)

Number of nodes in 3D hexahedron
"""
num_nodes(::Type{<:Hexahedron3D}) = 8
num_faces(::Type{<:Hexahedron3D}) = 6
num_edges(::Type{<:Hexahedron3D}) = 12

"""
$(SIGNATURES)

Cell-face node numbering for 3D hexahedron
"""
local_cellfacenodes(::Type{<:Hexahedron3D}) = _local_cellfacenodes_Hexahedron3D

"""
$(SIGNATURES)

Geometries of faces of 3D hexahedron
"""
facetype_of_cellface(::Type{<:Hexahedron3D}, k) = Quadrilateral2D

"""
$(SIGNATURES)

Geometries of faces of 3D parallelepiped
"""
facetype_of_cellface(::Type{<:Parallelepiped3D}, k) = Parallelogram2D

"""
$(SIGNATURES)

Cell-edge node numbering for 3D hexahedron
"""
local_celledgenodes(::Type{<:Hexahedron3D}) = _local_celledgenodes_Hexahedron3D

###############
### VOLUMES ###
###############

function Volume4ElemType(Coords, Nodes, ::Type{<:Vertex0D}, ::Type{<:ExtendableGrids.AbstractCoordinateSystem})
    function closure(item)
        return 0
    end
end

function Volume4ElemType(Coords, Nodes, ::Type{<:Edge1D}, ::Type{Cartesian1D})
    function closure(item)
        return abs(Coords[1, Nodes[2,item]] - Coords[1, Nodes[1,item]])
    end
end

function Volume4ElemType(Coords, Nodes, ::Type{<:Edge1D}, ::Type{Cartesian2D})
    function closure(item)
        return sqrt((Coords[1, Nodes[2,item]] - Coords[1, Nodes[1,item]]).^2 + (Coords[2, Nodes[2,item]] - Coords[2, Nodes[1,item]]).^2)
    end
end

function Volume4ElemType(Coords, Nodes, ::Type{<:Edge1D}, ::Type{Cartesian3D})
    function closure(item)
        return sqrt((Coords[1, Nodes[2,item]] - Coords[1, Nodes[1,item]]).^2 + (Coords[2, Nodes[2,item]] - Coords[2, Nodes[1,item]]).^2  + (Coords[3, Nodes[2,item]] - Coords[3, Nodes[1,item]]).^2)
    end
end

function Volume4ElemType(Coords, Nodes, ::Type{<:Triangle2D}, ::Type{Cartesian2D})
    function closure(item)
        return 1 // 2 * ( Coords[1, Nodes[1, item]] * (Coords[2, Nodes[2,item]] -  Coords[2, Nodes[3, item]])
                      +   Coords[1, Nodes[2, item]] * (Coords[2, Nodes[3,item]] -  Coords[2, Nodes[1, item]])
                      +   Coords[1, Nodes[3, item]] * (Coords[2, Nodes[1,item]] -  Coords[2, Nodes[2, item]]) )
    end
end

function Volume4ElemType(Coords, Nodes, ::Type{<:Quadrilateral2D}, ::Type{Cartesian2D})
    function closure(item)
        return 1//2 * (   (Coords[1, Nodes[1, item]] - Coords[1, Nodes[3, item]]) * (Coords[2, Nodes[2, item]] - Coords[2, Nodes[4, item]])
                    + (Coords[1, Nodes[4, item]] - Coords[1, Nodes[2, item]]) * (Coords[2, Nodes[1, item]] - Coords[2, Nodes[3, item]]) );

    end
end

function Volume4ElemType(Coords::Array{Tc,2}, Nodes, ::Type{<:Triangle2D}, ::Type{Cartesian3D}) where {Tc}
    d12 = zeros(Tc,3)
    d14 = zeros(Tc,3)
    function closure(item)
        # norm(cross(p(1)-p(2), p(1)-p(3)), 2)
        for k = 1 : 3
            d12[k] = Coords[k, Nodes[1, item]] - Coords[k, Nodes[2, item]]
            d14[k] = Coords[k, Nodes[1, item]] - Coords[k, Nodes[3, item]]
        end
        return sqrt((d12[2]*d14[3]-d12[3]*d14[2])^2 + (d12[3]*d14[1]-d12[1]*d14[3])^2 + (d12[1]*d14[2]-d12[2]*d14[1])^2) / 2;
    end
end

function Volume4ElemType(Coords::Array{Tc,2}, Nodes, ::Type{<:Parallelogram2D}, ::Type{Cartesian3D}) where {Tc}
    d12 = zeros(Tc,3)
    d14 = zeros(Tc,3)
    function closure(item)
        # norm(cross(p(1)-p(2), p(1)-p(3)), 2)
        for k = 1 : 3
            d12[k] = Coords[k, Nodes[1, item]] - Coords[k, Nodes[2, item]]
            d14[k] = Coords[k, Nodes[1, item]] - Coords[k, Nodes[4, item]]
        end
        return sqrt((d12[2]*d14[3]-d12[3]*d14[2])^2 + (d12[3]*d14[1]-d12[1]*d14[3])^2 + (d12[1]*d14[2]-d12[2]*d14[1])^2);
    end
end

function Volume4ElemType(Coords, Nodes, ::Type{<:Parallelepiped3D}, ::Type{Cartesian3D})
    function closure(item)
        return    ((Coords[1, Nodes[5, item]] - Coords[1, Nodes[1, item]]) * ( (Coords[2, Nodes[2, item]] - Coords[2, Nodes[1, item]]) * (Coords[3, Nodes[4, item]] - Coords[3, Nodes[1, item]]) - (Coords[2, Nodes[4, item]] - Coords[2, Nodes[1, item]]) * (Coords[3, Nodes[2, item]] - Coords[3, Nodes[1, item]]))
        + (Coords[2, Nodes[5, item]] - Coords[2, Nodes[1, item]]) * ( (Coords[3, Nodes[2, item]] - Coords[3, Nodes[1, item]]) * (Coords[1, Nodes[4, item]] - Coords[1, Nodes[1, item]]) - (Coords[1, Nodes[2, item]] - Coords[1, Nodes[1, item]]) * (Coords[3, Nodes[4, item]] - Coords[3, Nodes[1, item]]))
        + (Coords[3, Nodes[5, item]] - Coords[3, Nodes[1, item]]) * ( (Coords[1, Nodes[2, item]] - Coords[1, Nodes[1, item]]) * (Coords[2, Nodes[4, item]] - Coords[2, Nodes[1, item]]) - (Coords[2, Nodes[2, item]] - Coords[2, Nodes[1, item]]) * (Coords[1, Nodes[4, item]] - Coords[1, Nodes[1, item]])));

    end
end

function Volume4ElemType(Coords, Nodes, ::Type{<:Tetrahedron3D}, ::Type{Cartesian3D})
    function closure(item)
        return    1 // 6 * ((Coords[1, Nodes[4, item]] - Coords[1, Nodes[1, item]]) * ( (Coords[2, Nodes[2, item]] - Coords[2, Nodes[1, item]]) * (Coords[3, Nodes[3, item]] - Coords[3, Nodes[1, item]]) - (Coords[2, Nodes[3, item]] - Coords[2, Nodes[1, item]]) * (Coords[3, Nodes[2, item]] - Coords[3, Nodes[1, item]]))
        + (Coords[2, Nodes[4, item]] - Coords[2, Nodes[1, item]]) * ( (Coords[3, Nodes[2, item]] - Coords[3, Nodes[1, item]]) * (Coords[1, Nodes[3, item]] - Coords[1, Nodes[1, item]]) - (Coords[1, Nodes[2, item]] - Coords[1, Nodes[1, item]]) * (Coords[3, Nodes[3, item]] - Coords[3, Nodes[1, item]]))
        + (Coords[3, Nodes[4, item]] - Coords[3, Nodes[1, item]]) * ( (Coords[1, Nodes[2, item]] - Coords[1, Nodes[1, item]]) * (Coords[2, Nodes[3, item]] - Coords[2, Nodes[1, item]]) - (Coords[2, Nodes[2, item]] - Coords[2, Nodes[1, item]]) * (Coords[1, Nodes[3, item]] - Coords[1, Nodes[1, item]])));
    end
end

###############
### NORMALS ###
###############

function Normal4ElemType!(normal, Coords, Nodes, item, ::Type{<:Vertex0D}, ::Type{Cartesian2D})
    normal[1] = 0.0
    normal[2] = 0.0
end

function Normal4ElemType!(normal, Coords, Nodes, item, ::Type{<:Vertex0D}, ::Type{Cartesian1D})
    normal[1] = 1.0
end

function Normal4ElemType!(normal, Coords, Nodes, item, ::Type{<:Edge1D}, ::Type{Cartesian2D})
    # rotate tangent
    normal[1] = Coords[2, Nodes[2,item]] - Coords[2, Nodes[1,item]]
    normal[2] = Coords[1, Nodes[1,item]] - Coords[1, Nodes[2,item]]
    # divide by length
    normal ./= sqrt(normal[1]^2+normal[2]^2)
end

function Normal4ElemType!(normal, Coords, Nodes, item, ::Type{<:Quadrilateral2D}, ::Type{Cartesian3D})
    # cross(p(1)-p(2), p(1)-p(4)) / length
    normal[1]  = (Coords[2, Nodes[1, item]] - Coords[2, Nodes[2, item]]) * (Coords[3, Nodes[1, item]] - Coords[3, Nodes[4, item]])
    normal[1] -= (Coords[3, Nodes[1, item]] - Coords[3, Nodes[2, item]]) * (Coords[2, Nodes[1, item]] - Coords[2, Nodes[4, item]])
    normal[2]  = (Coords[3, Nodes[1, item]] - Coords[3, Nodes[2, item]]) * (Coords[1, Nodes[1, item]] - Coords[1, Nodes[4, item]])
    normal[2] -= (Coords[1, Nodes[1, item]] - Coords[1, Nodes[2, item]]) * (Coords[3, Nodes[1, item]] - Coords[3, Nodes[4, item]])
    normal[3]  = (Coords[1, Nodes[1, item]] - Coords[1, Nodes[2, item]]) * (Coords[2, Nodes[1, item]] - Coords[2, Nodes[4, item]])
    normal[3] -= (Coords[2, Nodes[1, item]] - Coords[2, Nodes[2, item]]) * (Coords[1, Nodes[1, item]] - Coords[1, Nodes[4, item]])
    # divide by length
    normal ./= sqrt(normal[1]^2+normal[2]^2+normal[3]^2)

    ## old version
    #d12 = @views Coords[:, Nodes[1, item]] - Coords[:, Nodes[2, item]]
    #d14 = @views Coords[:, Nodes[1, item]] - Coords[:, Nodes[4, item]]
    #normal[1] = d12[2]*d14[3]-d12[3]*d14[2]
    #normal[2] = d12[3]*d14[1]-d12[1]*d14[3]
    #normal[3] = d12[1]*d14[2]-d12[2]*d14[1]
end

function Normal4ElemType!(normal, Coords, Nodes, item, ::Type{<:Triangle2D}, ::Type{Cartesian3D})
    # cross(p(1)-p(2), p(1)-p(3)) / length
    normal[1]  = (Coords[2, Nodes[1, item]] - Coords[2, Nodes[2, item]]) * (Coords[3, Nodes[1, item]] - Coords[3, Nodes[3, item]])
    normal[1] -= (Coords[3, Nodes[1, item]] - Coords[3, Nodes[2, item]]) * (Coords[2, Nodes[1, item]] - Coords[2, Nodes[3, item]])
    normal[2]  = (Coords[3, Nodes[1, item]] - Coords[3, Nodes[2, item]]) * (Coords[1, Nodes[1, item]] - Coords[1, Nodes[3, item]])
    normal[2] -= (Coords[1, Nodes[1, item]] - Coords[1, Nodes[2, item]]) * (Coords[3, Nodes[1, item]] - Coords[3, Nodes[3, item]])
    normal[3]  = (Coords[1, Nodes[1, item]] - Coords[1, Nodes[2, item]]) * (Coords[2, Nodes[1, item]] - Coords[2, Nodes[3, item]])
    normal[3] -= (Coords[2, Nodes[1, item]] - Coords[2, Nodes[2, item]]) * (Coords[1, Nodes[1, item]] - Coords[1, Nodes[3, item]])
    # divide by length
    normal ./= sqrt(normal[1]^2+normal[2]^2+normal[3]^2)
end

################
### TANGENTS ###
################

function Tangent4ElemType!(tangent, Coords, Nodes, item, ::Type{<:Edge1D}, ::Type{Cartesian2D})
    tangent[1] = Coords[1,Nodes[2,item]] - Coords[1, Nodes[1,item]]
    tangent[2] = Coords[2,Nodes[2,item]] - Coords[2, Nodes[1,item]]
    # divide by length
    tangent ./= sqrt(tangent[1]^2+tangent[2]^2)
end

function Tangent4ElemType!(tangent, Coords, Nodes, item, ::Type{<:Edge1D}, ::Type{Cartesian3D})
    tangent[1] = Coords[1,Nodes[2,item]] - Coords[1, Nodes[1,item]]
    tangent[2] = Coords[2,Nodes[2,item]] - Coords[2, Nodes[1,item]]
    tangent[3] = Coords[3,Nodes[2,item]] - Coords[3, Nodes[1,item]]
    # divide by length
    tangent ./= sqrt(tangent[1]^2+tangent[2]^2+tangent[3]^2)
end


text:
"""
A `PointCloud` is a container for points sharing a common set of per-point
attributes.  Points within the cloud can be accessed by index, and vectors of
attributes may be accessed by name.  There is one required attribute
`position`, which is used to build a datastructure for fast spatial lookup of
k-nearest neighbours, or points within a fixed radius.

### Constructors

```julia
PointCloud(positions::Vector{FixedSizeArrays.Vec{3, AbstractFloat}})
PointCloud(positions::Matrix{AbstractFloat})
```

### Example

```julia
using FixedSizeArrays
using PointClouds

# Create positions of 3D Vec coordinates
positions = [Vec(i,-i,1.0) for i = 1:10]

# Create a PointCloud from positions
cloud = PointCloud(positions)

# Create new intensity attributes
cloud[:intensity] = collect(1:10)

# Create capture time for each point sample
cloud[:time] = [time() for i = 1:10]

# Find low intensity points
low_intensity_cloud = cloud[cloud[:intensity] .< 5]

# Find vector of time stamps per point
gps_time = low_intensity_cloud[:time]

# Find points within 5 units of [1,1,1]
nearby = cloud[inrange(cloud, [1,1,1], 5.0)]
```
"""

type PointCloud{Dim,T,SIndex}
    positions::Vector{SVector{Dim,T}}
    spatial_index::SIndex
    attributes::Dict{Symbol,Vector}
end

# Create a `PointCloud` for FixedSizeArrays Vec using a KDTree for spatial
# indexing based on `positions`.
# TODO should be able to change the spatial index
function PointCloud{Dim, T <: AbstractFloat}(positions::Vector{SVector{Dim, T}})
    PointCloud(positions, KDTree(positions), Dict{Symbol,Vector}())
end

# Create a `PointCloud` from an 3xN array of points, using a KDTree for spatial
# indexing based on `positions`.
PointCloud{T <: AbstractFloat}(points::Matrix{T}) = PointCloud(convert_positions(points))

"""
    convert_positions{T <: Real}(points::Matrix{T}) -> positions::Vector{FixedSizeArrays.Vec{3, Float64}}

Convert a standard julia array of points into a FixedSizeArray
"""
function convert_positions{T}(points::Matrix{T})
    ndims = size(points, 1)
    npoints = size(points, 2)
    if isbits(T)
        data = reinterpret(SVector{ndims, T}, points, (npoints, ))
    else
        data = SVector{ndims, T}[SVector{ndims, T}(points[:, i]) for i in 1:npoints]
    end
    return data
end

function show(io::IO, cloud::PointCloud)
    print(io, "PointCloud(N=$(length(cloud)), attributes=$(keys(cloud.attributes)))")
end

#------------------------
# Container functionality acting on columns
keys(cloud::PointCloud) = keys(cloud.attributes)
haskey(cloud::PointCloud, attrname::Symbol) = haskey(cloud.attributes, attrname)

function getindex(cloud::PointCloud, attrname::Symbol)
    haskey(cloud, attrname) || error("Point cloud $cloud has no attribute $attrname")
    cloud.attributes[attrname]
end

function setindex!(cloud::PointCloud, value::Vector, attrname::Symbol)
    if length(value) != length(cloud)
        error("length($attrname) = $(length(value)) not equal to number of points = $(length(cloud))")
    end
    cloud.attributes[attrname] = value
end

# Access to columns with "blessed" type for efficiency.
# TODO: Figure out how we can remove these, so that eg, cloud[:normal] is the
# definitive way to access normals.
# Note that the `positions()` is used rather than `position()`, to avoid a
# clash with `Base.position()` which means IO stream position... ugh.

"Return vector of point cloud positions"
positions(cloud::PointCloud)   = cloud.positions

"Return array of point cloud normals as fixed size vectors"
function normals{Dim,T,SIndex}(cloud::PointCloud{Dim,T,SIndex})
    cloud.attributes[:normal]::Vector{Vec{Dim,Float32}}
end

#--------------------------
# Container functionality acting on rows

length(cloud::PointCloud) = length(cloud.positions)
endof(cloud::PointCloud) = length(cloud.positions)

# Subset of a point cloud
function getindex{Dim,T,SIndex}(cloud::PointCloud{Dim,T,SIndex}, row_inds::AbstractVector)
    pos = positions(cloud)[row_inds]
    tree = KDTree(pos)  # TODO NEED TO GET SPATIAL INDEX FROM TYPE
    attrs = Dict{Symbol,Vector}()
    for (k,v) in cloud.attributes
        attrs[k] = v[row_inds]
    end
    PointCloud{Dim,T,SIndex}(pos, tree, attrs)
end

# Concatenate point clouds
function vcat{Dim,T,SIndex}(cloud1::PointCloud{Dim,T,SIndex}, clouds::PointCloud{Dim,T,SIndex}...)
    # Workaround for a weird compiler bug, which causes the deepcopy type assert to fail.
    # See https://github.com/JuliaLang/julia/issues/19041
    #attrs = deepcopy(cloud1.attributes)
    attrs = Dict{Symbol,Vector}()
    for (k,v) in cloud1.attributes
        attrs[k] = deepcopy(v)
    end
    pos = positions(cloud1)
    ks = Set(keys(attrs))
    for cloud in clouds
        pos = [pos; positions(cloud)]
        ki = Set(keys(cloud))
        ki == ks || error("Cannot join clouds with inconsistent attributes")
        for (k,v) in attrs
            append!(v, cloud.attributes[k])
        end
    end
    spatial_index = KDTree(pos)  # TODO NEED TO GET SPATIAL INDEX FROM TYPE
    PointCloud{Dim,T,SIndex}(pos, spatial_index, attrs)
end

"""
Group points into `M = unique(scanid)` point clouds, each containing points
with the same scan id.  Returns a dictionary mapping the unique scanid to each
cloud.
"""
function split_cloud(allpoints, scanid; min_points::Int=10000)
    unique_scans = unique(scanid)
    scans = Dict{eltype(unique_scans), typeof(allpoints)}()
    for id in unique_scans
        scans[id] = allpoints[id .== scanid]
    end
    scans
end

#--------------------------
# Spatial index lookup (knn and ball queries)
knn(cloud::PointCloud, points, k) = knn(cloud.spatial_index, points, k)
inrange(cloud::PointCloud, points, radius) = inrange(cloud.spatial_index, points, radius)

#------------------------
# Utilities for adding columns

"""
    pca_normals!(normals, query_points, cloud; neighbours::Int=10)

Compute normals of a point cloud via PCA

For each point in the 3xN array `query_points`, compute an approximate surface
normal using PCA on the local neighbourhood, using the closest `neighbours`
points.  `cloud` should support spatial lookup via the `knn()` function.
Normals are returned in the 3xN array `normals`.

TODO: Need a generalized version of this for the PCA covariance computation.
"""
function pca_normals!{V<:AbstractVector}(normals::Vector{V}, query_points, cloud; neighbours::Int=10)
    position = positions(cloud)
    for i = 1:length(query_points)
        @inbounds (inds, _) = knn(cloud, query_points[i], neighbours)
        isempty(inds) && continue
        cov_tmp = zeros(SMatrix{3,3,eltype(V)})
        for j = 1:length(inds)
            @inbounds p = position[inds[j]] - query_points[i]
            cov_tmp += p * p'
        end
        # FIXME: Avoid converting to Matrix when eig(SMatrix) is solid
        # See https://github.com/JuliaArrays/StaticArrays.jl/issues/80
        evals, evecs = eig(Symmetric(Matrix(cov_tmp)))
        # Select evec of smallest eval
        normals[i] = evecs[:,1]
    end
    return normals
end

"""
Add normals to point cloud in Float32 precision, compute using PCA.

See pca_normals!() for details.
"""
function add_normals!{Dim}(cloud::PointCloud{Dim}; kwargs...)
    # Use Float32 for estimated normals - heaps of precision at half the memory
    query_points = positions(cloud)
    normals = similar(query_points, SVector{Dim,Float32})
    pca_normals!(normals, query_points, cloud; kwargs...)
    cloud[:normals] = normals
    return cloud
end

ext: jl


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct IoK8sApiExtensionsV1beta1ReplicaSet <: SwaggerModel
    apiVersion::Any # spec type: Union{ Nothing, String } # spec name: apiVersion
    kind::Any # spec type: Union{ Nothing, String } # spec name: kind
    metadata::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApisMetaV1ObjectMeta } # spec name: metadata
    spec::Any # spec type: Union{ Nothing, IoK8sApiExtensionsV1beta1ReplicaSetSpec } # spec name: spec
    status::Any # spec type: Union{ Nothing, IoK8sApiExtensionsV1beta1ReplicaSetStatus } # spec name: status

    function IoK8sApiExtensionsV1beta1ReplicaSet(;apiVersion=nothing, kind=nothing, metadata=nothing, spec=nothing, status=nothing)
        o = new()
        validate_property(IoK8sApiExtensionsV1beta1ReplicaSet, Symbol("apiVersion"), apiVersion)
        setfield!(o, Symbol("apiVersion"), apiVersion)
        validate_property(IoK8sApiExtensionsV1beta1ReplicaSet, Symbol("kind"), kind)
        setfield!(o, Symbol("kind"), kind)
        validate_property(IoK8sApiExtensionsV1beta1ReplicaSet, Symbol("metadata"), metadata)
        setfield!(o, Symbol("metadata"), metadata)
        validate_property(IoK8sApiExtensionsV1beta1ReplicaSet, Symbol("spec"), spec)
        setfield!(o, Symbol("spec"), spec)
        validate_property(IoK8sApiExtensionsV1beta1ReplicaSet, Symbol("status"), status)
        setfield!(o, Symbol("status"), status)
        o
    end
end # type IoK8sApiExtensionsV1beta1ReplicaSet

const _property_map_IoK8sApiExtensionsV1beta1ReplicaSet = Dict{Symbol,Symbol}(Symbol("apiVersion")=>Symbol("apiVersion"), Symbol("kind")=>Symbol("kind"), Symbol("metadata")=>Symbol("metadata"), Symbol("spec")=>Symbol("spec"), Symbol("status")=>Symbol("status"))
const _property_types_IoK8sApiExtensionsV1beta1ReplicaSet = Dict{Symbol,String}(Symbol("apiVersion")=>"String", Symbol("kind")=>"String", Symbol("metadata")=>"IoK8sApimachineryPkgApisMetaV1ObjectMeta", Symbol("spec")=>"IoK8sApiExtensionsV1beta1ReplicaSetSpec", Symbol("status")=>"IoK8sApiExtensionsV1beta1ReplicaSetStatus")
Base.propertynames(::Type{ IoK8sApiExtensionsV1beta1ReplicaSet }) = collect(keys(_property_map_IoK8sApiExtensionsV1beta1ReplicaSet))
Swagger.property_type(::Type{ IoK8sApiExtensionsV1beta1ReplicaSet }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiExtensionsV1beta1ReplicaSet[name]))}
Swagger.field_name(::Type{ IoK8sApiExtensionsV1beta1ReplicaSet }, property_name::Symbol) =  _property_map_IoK8sApiExtensionsV1beta1ReplicaSet[property_name]

function check_required(o::IoK8sApiExtensionsV1beta1ReplicaSet)
    true
end

function validate_property(::Type{ IoK8sApiExtensionsV1beta1ReplicaSet }, name::Symbol, val)
end


text:
# ENV["JULIA_DEBUG"] = Metatheory
using Metatheory

fibo = @theory x y n begin
  x::Int + y::Int => x + y
  fib(n::Int) => (n < 2 ? n : :(fib($(n - 1)) + fib($(n - 2))))
end

params = SaturationParams(timeout = 60)
g = EGraph(:(fib(10)))
@time saturate!(g, fibo, params)

z = EGraph(:(fib(10)))
@time saturate!(z, fibo, params)

@testset "Fibonacci" begin
  @test 55 == extract!(g, astsize)
end


text:
using AbstractPlotting
using Test

# does this machine have a OPENGL?
const OPENGL = haskey(ENV, "OPENGL") || haskey(ENV, "GITLAB_CI") || try success(pipeline(`glxinfo`, `grep version`)) catch; false end # if it's Gitlab, it must be JuliaGPU

OPENGL && begin @info "OpenGL detected"; using GLMakie end
OPENGL || @warn "No OpenGL detected!  Software tests only."

# does this machine have FFMPEG?  We'll take it on faith if you tell us...
const fmp = haskey(ENV, "FFMPEG") || try; AbstractPlotting.@ffmpeg_env success(`$ffmpeg -version`); catch; false end;

const _MINIMAL = get(ENV, "ABSTRACTPLOTTING_MINIMAL", "false")

scene = scatter(rand(4))

include("conversions.jl")
include("quaternions.jl")
include("projection_math.jl")

@testset "basic functionality" begin
    scene = scatter(rand(4))
    @test scene[Axis].ticks.title_gap[] == 3
    scene[Axis].ticks.title_gap = 4
    @test scene[Axis].ticks.title_gap[] == 4
end

# if get(ENV, "IS_TRAVIS_CI", "false") == "false"
#   exit(0);
# end

using MakieGallery

if _MINIMAL == "false"

    # Load all entries in the database

    database = MakieGallery.load_database()

   ## Exceptions are made on basis of title and not index,
    # because index may change as MakieGallery changes.

   # All these require FFMpeg and need to save,
   # and are therefore ignored
   ffmpeg_exs  = [
               "Animation", "Lots of Heatmaps","Animated Scatter",
               "Chess Game", "Record Video",  "Animated surface and wireframe",
               "Moire", "pong", "pulsing marker", "Travelling wave",
               "Type recipe for molecule simulation"
   ]

   # All these require GLMakie and so are ignored
   glmakie_exs = [
               "Textured Mesh", "Load Mesh", "Wireframe of a Mesh",
               "FEM mesh 3D", "Normals of a Cat", "Line GIF"
   ]

   # Requires GDAL (a GL package) so ignored
   gdal_exs = [
           "WorldClim visualization"
   ]

    # Requires GLMakie and ModernGL, so ignored
   moderngl_exs = [
               "Explicit frame rendering"
   ]

   # use Stepper plots, which save, so ignored
   save_exs = [
             "Axis theming", "Labels", "Color Legend",
             "Stepper demo"
   ]

   # hopefullly fixed by next tag of MakieGallery (already in AbstractPlotting
   color_exs = ["colormaps"]

   # curl fails with this for some reason, so it has been ignored.
   curl_exs = ["Earth & Ships"]

   if OPENGL
       if fmp
           # all test infrastructure in place - no exceptions!
           exc_str = Set()
       else
           exc_str = Set(ffmpeg_exs)
       end
   else
       # combine all exceptions into a single Set
       exc_str = exc_str = union(ffmpeg_exs, glmakie_exs, gdal_exs, moderngl_exs, save_exs, color_exs, curl_exs)
   end

else

    # Load only short tests

    database = MakieGallery.load_tests()

    # This one is broken, all the others are fine.  I don't know why this fails.  It doesn't fail on my machine :P
    exc_str = Set(["Comparing contours, image, surfaces and heatmaps"])

end

if !OPENGL # run software only tests...

    @testset "Gallery short tests" begin

        # iterate over database
        @testset "$(database[i].title) (#$i)" for i in 1:length(database)

               # skip if the title is in the list of exceptions
              if database[i].title ∈ exc_str

                 print(
                 "Skipping "
                 * database[i].title *
                 "\n(removed from tests explicitly)\n"
                 )

                 continue

               end

               @debug("Running " * database[i].title * "\n(index $i)\n")
               # evaluate the entry
               @test_nowarn MakieGallery.eval_example(database[i]);

       end

    end

else # full MakieGallery comparisons here

    using GLMakie

    @info("Running full tests - artifacts will be stored!")

    for exc in exc_str

        printstyled("Excluded ", color = :yellow, bold = true)
        println(exc)

    end

    filter!(entry -> !(entry.title in exc_str), database)

    tested_diff_path = joinpath(@__DIR__, "tested_different")
    test_record_path = joinpath(@__DIR__, "test_recordings")
    rm(tested_diff_path, force = true, recursive = true)
    mkpath(tested_diff_path)
    rm(test_record_path, force = true, recursive = true)
    mkpath(test_record_path)

    examples = MakieGallery.record_examples(test_record_path)

    MakieGallery.run_comparison(test_record_path, tested_diff_path)

end


text:
# Clayton, N. S. and Dickinson, A. (1999), Motivational control of caching
# behaviour in the scrub jay (Aphelocoma coerulescens), Animal Behaviour,
# 2:435-444

function _summarize(::Experiment{:Clayton99C_exp1}, results)
    combine(groupby(results, [:group, :foodtype, :action, :stage]),
       df -> DataFrame(μ = mean(df.counts), sem = sem(df.counts),
                       n = length(df.counts)))
end
function statistical_tests(exp::Experiment{:Clayton99C_exp1}, data)
    tests = Test[]
    push!(tests, Test("prefeding effect on stones cache",
                      ttest(@where(data, :action .== "cache",
                                         :foodtype .== "stone",
                                         :group .== "S",
                                         :stage .== 1).counts,
                            @where(data, :action .== "cache",
                                         :foodtype .== "stone",
                                         (:group .== "pP") .| (:group .== "pP+S")).counts,
                            paired = false), # why is the degree 1,20 in the paper?
                      exp.tests))
    push!(tests, Test("prefeding effect on peanuts cache",
                      ttest(@where(data, :action .== "cache",
                                         :foodtype .== "peanut",
                                         :group .== "P",
                                         :stage .== 1).counts,
                            @where(data, :action .== "cache",
                                         :foodtype .== "peanut",
                                         (:group .== "pP") .| (:group .== "pP+S")).counts,
                            paired = false), # why is the degree 1,20 in the paper?
                      exp.tests))
    push!(tests, Test("precaching effect on stones cache",
                      ttest(@where(data, :action .== "cache",
                                         :foodtype .== "stone",
                                         :group .== "S", :trial .== "S",
                                         :stage .== 1).counts,
                            @where(data, :action .== "cache",
                                         :foodtype .== "stone",
                                         :group .== "S",
                                         :stage .== 2).counts,
                            paired = true),
                      exp.tests))
    push!(tests, Test("precaching effect on peanuts cache",
                      ttest(@where(data, :action .== "cache",
                                         :foodtype .== "peanut",
                                         (:group .== "S") .| (:group .== "pP+S"),
                                         :stage .== 2).counts,
                            @where(data, :action .== "cache",
                                         :foodtype .== "peanut",
                                         :group .== "P",
                                         :stage .== 1).counts,
                            paired = false),
                      exp.tests))
    id = "interactions"
    push!(tests, Test(id,
                      splitplotanova(@where(data, :stage .== 2, :action .== "cache"),
                                     exp.tests[id].locals,
                                     betweenfactors = [:group, :order],
                                     withinfactors = :foodtype,
                                     id = :id, data = :counts),
                      exp.tests))
    id = "stones cache in stage 2"
    push!(tests, Test(id,
                      anova(@where(data, :stage .== 2, :action .== "cache",
                                         :foodtype .== "stone"),
                            @formula(counts ~ group), exp.tests[id].locals),
                      exp.tests))
    id = "peanuts cache in stage 2"
    push!(tests, Test(id,
                      anova(@where(data, :stage .== 2, :action .== "cache",
                                         :foodtype .== "peanut"),
                            @formula(counts ~ group), exp.tests[id].locals),
                      exp.tests))
    TestCollection(tests)
end
function run!(::Experiment{:Clayton99C_exp1}, models)
    results = DataFrame(group = String[], trial = String[], order = String[],
                        foodtype = String[], action = String[],
                        stage = Int[], id = Int[], counts = Int[])
    id = 0
    for group in ["pP", "S", "P", "pP+S"]
        for i in 1:6
            id += 1
            m = models[id]
            trials = ["S", "P"]
            if i > 3
                reverse!(trials)
                order = "PS"
            else
                order = "SP"
            end
            for trial in trials
                remove!(m, Any)
                wait!(m, 4u"hr")
                tray = Tray()
                add!(m, tray) # stage 1
                if group == "pP"
                    add!(m, PowderedPeanut, 50)
                elseif group == "S"
                    add!(m, Stone, 50)
                elseif group == "P"
                    add!(m, Peanut, 50)
                else
                    add!(m, PowderedPeanut, 50)
                    add!(m, Stone, 50)
                end
                wait!(m, 1u"hr")
                if group == "P"
                    cache = countcache(tray, Peanut)
                    push!(results, [group, trial, order, "peanut", "cache", 1, id, cache])
                    push!(results, [group, trial, order, "peanut", "eat", 1, id,
                                    50 - countfooditems(m, Peanut)])
                elseif group == "S"
                    push!(results, [group, trial, order, "stone", "cache", 1, id,
                                    countcache(tray, Stone)])
                end
                remove!(m, Any)
                pilfer!(tray)
                add!(m, tray) # stage 2
                foodtype = trial == "S" ? Stone : Peanut
                add!(m, foodtype, 50)
                wait!(m, 1u"hr")
                cache = countcache(tray, foodtype)
                push!(results, [group, trial, order, lowercase(string(foodtype)), "cache", 2,
                                id, cache])
                if foodtype == Peanut
                    push!(results, [group, trial, order, "peanut", "eat", 2, id,
                                    50 - countfooditems(m, Peanut)])
                end
                remove!(m, Any)
                add!(m, MaintenanceDiet)
                wait!(m, 18u"hr" + 1u"d")
            end
        end
    end
    results
end

function _summarize(::Experiment{:Clayton99C_exp2}, results)
    combine(groupby(results, [:group, :action, :stage]),
       df -> DataFrame(μ = mean(df.counts), sem = sem(df.counts),
                       n = length(df.counts)))
end
function statistical_tests(exp::Experiment{:Clayton99C_exp2}, data)
    tests = Test[]
    for action in ("cache", "eat")
        push!(tests, Test(action,
                          splitplotanova(@where(data, :action .== action),
                                         exp.tests[action].locals,
                                         betweenfactors = [:group, :foodtype2, :traychanged],
                                         withinfactors = :stage,
                                         id = :id, data = :counts),
                          exp.tests))
        for stage in 1:2
            id = "$action stage $stage"
            push!(tests, Test(id,
                              anova(@where(data, :action .== action,
                                           :stage .== stage),
                                    @formula(counts ~ group * foodtype2 * traychanged), exp.tests[id].locals),
                              exp.tests))
        end
    end
    for group in ("same", "different")
        for action in ("cache", "eat")
#   with stage as within factor I don't get the same Df as in the paper
            id = "$group group stage 1 vs stage 2 $action"
            push!(tests, Test(id,
                              splitplotanova(@where(data, :action .== action,
                                             :group .== group),
                                             exp.tests[id].locals,
                                             betweenfactors = [:foodtype2, :traychanged],
                                             withinfactors = :stage,
                                             id = :id, data = :counts),
#                               anova(@where(data, :action .== action,
#                                            :group .== group),
#                               @formula(counts ~ stage*foodtype2*traychanged)),
                              exp.tests))
        end
    end
    TestCollection(tests)
end
function run!(::Experiment{:Clayton99C_exp2}, models)
    results = DataFrame(group = String[], action = String[], stage = Int[],
                        foodtype2 = String[], traychanged = Bool[],
                        id = Int[], counts = Int[])
    id = 0
    for group in ["same", "different"]
        for i in 1:12
            id += 1
            m = models[id]
            foodtypes = [Peanut, Kibble]
            i > 6 && reverse!(foodtypes)
            remove!(m, Any)
            wait!(m, 4u"hr")
            tray = Tray(1)
            foodtype2 = "same" == group ? string(foodtypes[1]) : string(foodtypes[2])
            traychanged = i % 2 == 0
            for stage in 1:2
                foodtypeindex = 1
                if stage == 2
                    traychanged && (tray = Tray(2)) # Did they put it at the same place as the first tray or not?
                    foodtypeindex = "same" == group ? 1 : 2
                end
                add!(m, tray)
                add!(m, foodtypes[foodtypeindex], 50)
                wait!(m, 1u"hr")
                cache = countcache(tray)
                push!(results, [group, "cache", stage, foodtype2, traychanged, id, cache])
                push!(results, [group, "eat", stage, foodtype2, traychanged, id,
                                50 - countfooditems(m, foodtypes[foodtypeindex])])
                remove!(m, Any)
                pilfer!(tray)
            end
        end
    end
    results
end

function _summarize(::Experiment{:Clayton99C_exp3}, results)
    combine(groupby(results, [:group, :action]),
       df -> DataFrame(μ = mean(df.counts), sem = sem(df.counts),
                       n = length(df.counts)))
end
function statistical_tests(exp::Experiment{:Clayton99C_exp3}, data)
    tests = Test[]
    for action in ("cache", "eat")
        push!(tests, Test(action,
                          anova(@where(data, :exp .== "a", :action .== action),
                                @formula(counts ~ group * foodtype),
                                exp.tests[action].locals),
                          exp.tests))
    end
    push!(tests, Test("within-subject prefeeding",
                      is_significantly_bigger(@where(data, :exp .== "b"),
                                              conditionkey = :group,
                                              condbigger = "nothing",
                                              condsmaller = "prefed",
                                              valuekey = :counts,
                                              paired = true),
                      exp.tests))
    TestCollection(tests)
end
function run!(::Experiment{:Clayton99C_exp3}, models; samebirds = true)
    results = DataFrame(exp = String[],
                        group = String[], action = String[], foodtype = String[],
                        id = Int[], counts = Int[])
    id = 0
    for group in ["same", "different", "none"]
        for i in 1:8
            id += 1
            m = models[id]
            foodtypes = [Peanut, Kibble]
            i > 6 && reverse!(foodtypes)
            remove!(m, Any)
            wait!(m, 4u"hr")
            tray = Tray(1)
            add!(m, tray)
            group != "none" && add!(m, foodtypes[1] == Peanut ? PowderedPeanut :
                                    PowderedKibble, 50)
            wait!(m, 1u"hr")
            remove!(m, Any)
            add!(m, tray)
            foodtypeindex = group == "same" ? 1 : 2
            add!(m, foodtypes[foodtypeindex], 50)
            wait!(m, 1u"hr")
            cache = countcache(tray)
            push!(results, ["a", group, "cache",
                            string(foodtypes[foodtypeindex]), id, cache])
            push!(results, ["a", group, "eat", string(foodtypes[foodtypeindex]), id,
                            50 - countfooditems(m, foodtypes[foodtypeindex])])
            remove!(m, Any)
        end
    end
    samebirds && (id = 0)
    for i in 1:23
        id += 1
        m = models[id]
        add!(m, MaintenanceDiet)
        wait!(m, 100u"d")
        for day in 1:2
            remove!(m, Any)
            wait!(m, 4u"hr")
            prefed = (day == 1 && i <= 11) || (day == 2 && i > 11)
            prefed && add!(m, PowderedPeanut, 50)
            wait!(m, 1u"hr") # not mentioned in paper
            remove!(m, Any)
            tray = Tray()
            add!(m, tray)
            add!(m, Peanut, 50) # amount not mentioned in paper
            wait!(m, 1u"hr") # not mentioned in paper
            cache = countcache(tray)
            push!(results, ["b", prefed ? "prefed" : "nothing",
                            "cache", i <= 11 ? "first" : "second", id, cache])
            remove!(m, Any)
            add!(m, MaintenanceDiet)
            wait!(m, 18u"hr")
        end
    end
    results
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

@doc raw"""NodeConfigStatus describes the status of the config assigned by Node.Spec.ConfigSource.

    IoK8sApiCoreV1NodeConfigStatus(;
        active=nothing,
        assigned=nothing,
        error=nothing,
        lastKnownGood=nothing,
    )

    - active::IoK8sApiCoreV1NodeConfigSource : Active reports the checkpointed config the node is actively using. Active will represent either the current version of the Assigned config, or the current LastKnownGood config, depending on whether attempting to use the Assigned config results in an error.
    - assigned::IoK8sApiCoreV1NodeConfigSource : Assigned reports the checkpointed config the node will try to use. When Node.Spec.ConfigSource is updated, the node checkpoints the associated config payload to local disk, along with a record indicating intended config. The node refers to this record to choose its config checkpoint, and reports this record in Assigned. Assigned only updates in the status after the record has been checkpointed to disk. When the Kubelet is restarted, it tries to make the Assigned config the Active config by loading and validating the checkpointed payload identified by Assigned.
    - error::String : Error describes any problems reconciling the Spec.ConfigSource to the Active config. Errors may occur, for example, attempting to checkpoint Spec.ConfigSource to the local Assigned record, attempting to checkpoint the payload associated with Spec.ConfigSource, attempting to load or validate the Assigned config, etc. Errors may occur at different points while syncing config. Earlier errors (e.g. download or checkpointing errors) will not result in a rollback to LastKnownGood, and may resolve across Kubelet retries. Later errors (e.g. loading or validating a checkpointed config) will result in a rollback to LastKnownGood. In the latter case, it is usually possible to resolve the error by fixing the config assigned in Spec.ConfigSource. You can find additional information for debugging by searching the error message in the Kubelet log. Error is a human-readable description of the error state; machines can check whether or not Error is empty, but should not rely on the stability of the Error text across Kubelet versions.
    - lastKnownGood::IoK8sApiCoreV1NodeConfigSource : LastKnownGood reports the checkpointed config the node will fall back to when it encounters an error attempting to use the Assigned config. The Assigned config becomes the LastKnownGood config when the node determines that the Assigned config is stable and correct. This is currently implemented as a 10-minute soak period starting when the local record of Assigned config is updated. If the Assigned config is Active at the end of this period, it becomes the LastKnownGood. Note that if Spec.ConfigSource is reset to nil (use local defaults), the LastKnownGood is also immediately reset to nil, because the local default config is always assumed good. You should not make assumptions about the node's method of determining config stability and correctness, as this may change or become configurable in the future.
"""
mutable struct IoK8sApiCoreV1NodeConfigStatus <: SwaggerModel
    active::Any # spec type: Union{ Nothing, IoK8sApiCoreV1NodeConfigSource } # spec name: active
    assigned::Any # spec type: Union{ Nothing, IoK8sApiCoreV1NodeConfigSource } # spec name: assigned
    error::Any # spec type: Union{ Nothing, String } # spec name: error
    lastKnownGood::Any # spec type: Union{ Nothing, IoK8sApiCoreV1NodeConfigSource } # spec name: lastKnownGood

    function IoK8sApiCoreV1NodeConfigStatus(;active=nothing, assigned=nothing, error=nothing, lastKnownGood=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1NodeConfigStatus, Symbol("active"), active)
        setfield!(o, Symbol("active"), active)
        validate_property(IoK8sApiCoreV1NodeConfigStatus, Symbol("assigned"), assigned)
        setfield!(o, Symbol("assigned"), assigned)
        validate_property(IoK8sApiCoreV1NodeConfigStatus, Symbol("error"), error)
        setfield!(o, Symbol("error"), error)
        validate_property(IoK8sApiCoreV1NodeConfigStatus, Symbol("lastKnownGood"), lastKnownGood)
        setfield!(o, Symbol("lastKnownGood"), lastKnownGood)
        o
    end
end # type IoK8sApiCoreV1NodeConfigStatus

const _property_map_IoK8sApiCoreV1NodeConfigStatus = Dict{Symbol,Symbol}(Symbol("active")=>Symbol("active"), Symbol("assigned")=>Symbol("assigned"), Symbol("error")=>Symbol("error"), Symbol("lastKnownGood")=>Symbol("lastKnownGood"))
const _property_types_IoK8sApiCoreV1NodeConfigStatus = Dict{Symbol,String}(Symbol("active")=>"IoK8sApiCoreV1NodeConfigSource", Symbol("assigned")=>"IoK8sApiCoreV1NodeConfigSource", Symbol("error")=>"String", Symbol("lastKnownGood")=>"IoK8sApiCoreV1NodeConfigSource")
Base.propertynames(::Type{ IoK8sApiCoreV1NodeConfigStatus }) = collect(keys(_property_map_IoK8sApiCoreV1NodeConfigStatus))
Swagger.property_type(::Type{ IoK8sApiCoreV1NodeConfigStatus }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiCoreV1NodeConfigStatus[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1NodeConfigStatus }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1NodeConfigStatus[property_name]

function check_required(o::IoK8sApiCoreV1NodeConfigStatus)
    true
end

function validate_property(::Type{ IoK8sApiCoreV1NodeConfigStatus }, name::Symbol, val)
end


text:
using Statistics
using StatsBase:sample
using EvoTrees
using BenchmarkTools
# using BSON: @save, @load
using Plots
using NetworkLayout

# prepare a dataset
features = rand(Int(1.25e5), 100)
# features = rand(100, 10)
X = features
Y = rand(size(X, 1))
𝑖 = collect(1:size(X, 1))

# train-eval split
𝑖_sample = sample(𝑖, size(𝑖, 1), replace=false)
train_size = 0.8
𝑖_train = 𝑖_sample[1:floor(Int, train_size * size(𝑖, 1))]
𝑖_eval = 𝑖_sample[floor(Int, train_size * size(𝑖, 1)) + 1:end]

X_train, X_eval = X[𝑖_train, :], X[𝑖_eval, :]
Y_train, Y_eval = Y[𝑖_train], Y[𝑖_eval]

#############################
# CPU - linear
#############################
params1 = EvoTreeRegressor(T=Float32,
    loss=:linear, metric=:none,
    nrounds=10,
    λ=1.0, γ=0.5, η=0.2,
    max_depth=5, min_weight=1.0,
    rowsample=0.5, colsample=0.5, nbins=32,
    device="cpu")

# for 100k 10 rounds: 410.477 ms (44032 allocations: 182.68 MiB)
# for 100k 100 rounds: 2.177 s (404031 allocations: 626.45 MiB)
# for 1.25e6 no eval: 6.244 s (73955 allocations: 2.18 GiB)
# for 1.25e6 mse with eval data:  6.345 s (74009 allocations: 2.18 GiB)
@time model = fit_evotree(params1, X_train, Y_train);
EvoTrees.save(model, "data/model_linear_cpu.bson")
plot(model, 2)

params1 = EvoTreeRegressor(T=Float32,
    loss=:linear, metric=:none,
    nrounds=10,
    λ=1.0, γ=0.5, η=0.2,
    max_depth=5, min_weight=1.0,
    rowsample=0.5, colsample=0.5, nbins=32,
    device="gpu")

# for 100k 10 rounds: 410.477 ms (44032 allocations: 182.68 MiB)
# for 100k 100 rounds: 2.177 s (404031 allocations: 626.45 MiB)
# for 1.25e6 no eval: 6.244 s (73955 allocations: 2.18 GiB)
# for 1.25e6 mse with eval data:  6.345 s (74009 allocations: 2.18 GiB)
@time model = fit_evotree(params1, X_train, Y_train);
EvoTrees.save(model, "data/model_linear_gpu.bson")

#############################
# CPU - Gaussian
#############################
params1 = EvoTreeGaussian(T=Float64,
    loss=:gaussian, metric=:none,
    nrounds=10,
    λ=1.0, γ=0.1, η=0.1,
    max_depth=5, min_weight=1.0,
    rowsample=0.5, colsample=0.5, nbins=64,
    device="gpu")

@time model = fit_evotree(params1, X_train, Y_train);
EvoTrees.save(model, "data/model_gaussian_gpu.bson")


text:
const AdaBoostRegressor_ = sken(:AdaBoostRegressor)
@sk_reg mutable struct AdaBoostRegressor <: MMI.Deterministic
    base_estimator::Any    = nothing
    n_estimators::Int      = 50::(_ > 0)
    learning_rate::Float64 = 1.0::(_ > 0)
    loss::String           = "linear"::(_ in ("linear","square","exponential"))
    random_state::Any      = nothing
end
MMI.fitted_params(model::AdaBoostRegressor, (f, _, _)) = (
    estimators           = f.estimators_,
    estimator_weights    = f.estimator_weights_,
    estimator_errors     = f.estimator_errors_,
    feature_importances  = f.feature_importances_
    )
add_human_name_trait(AdaBoostRegressor, "AdaBoost ensemble regression")

# ----------------------------------------------------------------------------
const AdaBoostClassifier_ = sken(:AdaBoostClassifier)
@sk_clf mutable struct AdaBoostClassifier <: MMI.Probabilistic
    base_estimator::Any    = nothing
    n_estimators::Int      = 50::(_ > 0)
    learning_rate::Float64 = 1.0::(_ > 0)
    algorithm::String      = "SAMME.R"::(_ in ("SAMME", "SAMME.R"))
    random_state::Any      = nothing
end
MMI.fitted_params(m::AdaBoostClassifier, (f, _, _)) = (
    estimators        = f.estimators_,
    estimator_weights = f.estimator_weights_,
    estimator_errors  = f.estimator_errors_,
    classes           = f.classes_,
    n_classes         = f.n_classes_
    )
meta(AdaBoostClassifier,
    input   = Table(Continuous),
    target  = AbstractVector{<:Finite},
    weights = false,
    )

# ============================================================================
const BaggingRegressor_ = sken(:BaggingRegressor)
@sk_reg mutable struct BaggingRegressor <: MMI.Deterministic
    base_estimator::Any      = nothing
    n_estimators::Int        = 10::(_>0)
    max_samples::Union{Int,Float64}  = 1.0::(_>0)
    max_features::Union{Int,Float64} = 1.0::(_>0)
    bootstrap::Bool          = true
    bootstrap_features::Bool = false
    oob_score::Bool          = false
    warm_start::Bool         = false
    n_jobs::Option{Int}      = nothing
    random_state::Any        = nothing
    verbose::Int             = 0
end
MMI.fitted_params(model::BaggingRegressor, (f, _, _)) = (
    estimators          = f.estimators_,
    estimators_samples  = f.estimators_samples_,
    estimators_features = f.estimators_features_,
    oob_score           = model.oob_score ? f.oob_score_ : nothing,
    oob_prediction      = model.oob_score ? f.oob_prediction_ : nothing
    )
add_human_name_trait(BaggingRegressor, "bagging ensemble regressor")

# ----------------------------------------------------------------------------
const BaggingClassifier_ = sken(:BaggingClassifier)
@sk_clf mutable struct BaggingClassifier <: MMI.Probabilistic
    base_estimator::Any      = nothing
    n_estimators::Int        = 10::(_>0)
    max_samples::Union{Int,Float64}  = 1.0::(_>0)
    max_features::Union{Int,Float64} = 1.0::(_>0)
    bootstrap::Bool          = true
    bootstrap_features::Bool = false
    oob_score::Bool          = false
    warm_start::Bool         = false
    n_jobs::Option{Int}      = nothing
    random_state::Any        = nothing
    verbose::Int             = 0
end
MMI.fitted_params(m::BaggingClassifier, (f, _, _)) = (
    base_estimator        = f.base_estimator_,
    estimators            = f.estimators_,
    estimators_samples    = f.estimators_samples_,
    estimators_features   = f.estimators_features_,
    classes               = f.classes_,
    n_classes             = f.n_classes_,
    oob_score             = m.oob_score ? f.oob_score_ : nothing,
    oob_decision_function = m.oob_score ? f.oob_decision_function_ : nothing
    )
meta(BaggingClassifier,
    input   = Table(Continuous),
    target  = AbstractVector{<:Finite},
    weights = false,
    human_name   = "bagging ensemble classifier"
    )

# ============================================================================
const GradientBoostingRegressor_ = sken(:GradientBoostingRegressor)
@sk_reg mutable struct GradientBoostingRegressor <: MMI.Deterministic
    loss::String                    = "ls"::(_ in ("ls","lad","huber","quantile"))
    learning_rate::Float64          = 0.1::(_>0)
    n_estimators::Int               = 100::(_>0)
    subsample::Float64              = 1.0::(_>0)
    criterion::String               = "friedman_mse"::(_ in ("mse","mae","friedman_mse"))
    min_samples_split::Union{Int,Float64} = 2::(_>0)
    min_samples_leaf::Union{Int,Float64}  = 1::(_>0)
    min_weight_fraction_leaf::Float64     = 0.0::(_≥0)
    max_depth::Int                 = 3::(_>0)
    min_impurity_decrease::Float64 = 0.0::(_≥0)
    init::Any                      = nothing
    random_state::Any              = nothing
    max_features::Union{Int,Float64,String,Nothing} = nothing::(_===nothing || (isa(_,String) && (_ in ("auto","sqrt","log2"))) || _ > 0)
    alpha::Float64                 = 0.9::(_>0)
    verbose::Int                   = 0
    max_leaf_nodes::Option{Int}    = nothing::(_===nothing || _>0)
    warm_start::Bool               = false
#    presort::Union{Bool,String}    = "auto"::(isa(_, Bool) || _ == "auto")
    validation_fraction::Float64   = 0.1::(_>0)
    n_iter_no_change::Option{Int}  = nothing
    tol::Float64                   = 1e-4::(_>0)
end
MMI.fitted_params(m::GradientBoostingRegressor, (f, _, _)) = (
    feature_importances = f.feature_importances_,
    train_score         = f.train_score_,
    loss                = f.loss_,
    init                = f.init_,
    estimators          = f.estimators_,
    oob_improvement     = m.subsample < 1 ? f.oob_improvement_ : nothing
    )
add_human_name_trait(GradientBoostingRegressor, "gradient boosting ensemble regression")

# ----------------------------------------------------------------------------
const GradientBoostingClassifier_ = sken(:GradientBoostingClassifier)
@sk_clf mutable struct GradientBoostingClassifier <: MMI.Probabilistic
    loss::String                    = "deviance"::(_ in ("deviance","exponential"))
    learning_rate::Float64          = 0.1::(_>0)
    n_estimators::Int               = 100::(_>0)
    subsample::Float64              = 1.0::(_>0)
    criterion::String               = "friedman_mse"::(_ in ("mse","mae","friedman_mse"))
    min_samples_split::Union{Int,Float64} = 2::(_>0)
    min_samples_leaf::Union{Int,Float64}  = 1::(_>0)
    min_weight_fraction_leaf::Float64     = 0.0::(_≥0)
    max_depth::Int                 = 3::(_>0)
    min_impurity_decrease::Float64 = 0.0::(_≥0)
    init::Any                      = nothing
    random_state::Any              = nothing
    max_features::Union{Int,Float64,String,Nothing} = nothing::(_===nothing || (isa(_,String) && (_ in ("auto","sqrt","log2"))) || _ > 0)
    verbose::Int                   = 0
    max_leaf_nodes::Option{Int}    = nothing::(_===nothing || _>0)
    warm_start::Bool               = false
#    presort::Union{Bool,String}    = "auto"::(isa(_, Bool) || _ == "auto")
    validation_fraction::Float64   = 0.1::(_>0)
    n_iter_no_change::Option{Int}  = nothing
    tol::Float64                   = 1e-4::(_>0)
end
MMI.fitted_params(m::GradientBoostingClassifier, (f, _, _)) = (
    n_estimators        = f.n_estimators_,
    feature_importances = f.feature_importances_,
    train_score         = f.train_score_,
    loss                = f.loss_,
    init                = f.init_,
    estimators          = f.estimators_,
    oob_improvement     = m.subsample < 1 ? f.oob_improvement_ : nothing
    )
meta(GradientBoostingClassifier,
    input   = Table(Continuous),
    target  = AbstractVector{<:Finite},
    weights = false
    )

# ============================================================================
const RandomForestRegressor_ = sken(:RandomForestRegressor)
@sk_reg mutable struct RandomForestRegressor <: MMI.Deterministic
    n_estimators::Int              = 100::(_ > 0)
    criterion::String              = "mse"::(_ in ("mae", "mse"))
    max_depth::Option{Int}         = nothing::(_ === nothing || _ > 0)
    min_samples_split::Union{Int,Float64} = 2::(_ > 0)
    min_samples_leaf::Union{Int,Float64}  = 1::(_ > 0)
    min_weight_fraction_leaf::Float64     = 0.0::(_ ≥ 0)
    max_features::Union{Int,Float64,String,Nothing} = "auto"::(_ === nothing || (isa(_, String) && (_ in ("auto","sqrt","log2"))) || _ > 0)
    max_leaf_nodes::Option{Int}    = nothing::(_ === nothing || _ > 0)
    min_impurity_decrease::Float64 = 0.0::(_ ≥ 0)
    bootstrap::Bool                = true
    oob_score::Bool                = false
    n_jobs::Option{Int}            = nothing
    random_state::Any              = nothing
    verbose::Int                   = 0
    warm_start::Bool               = false
    ccp_alpha::Float64             =0.0::(_ ≥ 0)
    max_samples::Union{Nothing,Float64,Int} =
        nothing::(_ === nothing || (_ ≥ 0 && (_ isa Integer || _ ≤ 1)))
end
MMI.fitted_params(model::RandomForestRegressor, (f, _, _)) = (
    estimators          = f.estimators_,
    feature_importances = f.feature_importances_,
    n_features          = f.n_features_,
    n_outputs           = f.n_outputs_,
    oob_score           = model.oob_score ? f.oob_score_ : nothing,
    oob_prediction      = model.oob_score ? f.oob_prediction_ : nothing
    )
meta(RandomForestRegressor,
    input   = Table(Count,Continuous),
    target  = AbstractVector{Continuous},
    weights = false
    )

# ----------------------------------------------------------------------------
const RandomForestClassifier_ = sken(:RandomForestClassifier)
@sk_clf mutable struct RandomForestClassifier <: MMI.Probabilistic
    n_estimators::Int              = 100::(_ > 0)
    criterion::String              = "gini"::(_ in ("gini","entropy"))
    max_depth::Option{Int}         = nothing::(_ === nothing || _ > 0)
    min_samples_split::Union{Int,Float64} = 2::(_ > 0)
    min_samples_leaf::Union{Int,Float64}  = 1::(_ > 0)
    min_weight_fraction_leaf::Float64     = 0.0::(_ ≥ 0)
    max_features::Union{Int,Float64,String,Nothing} = "auto"::(_ === nothing || (isa(_, String) && (_ in ("auto","sqrt","log2"))) || _ > 0)
    max_leaf_nodes::Option{Int}    = nothing::(_ === nothing || _ > 0)
    min_impurity_decrease::Float64 = 0.0::(_ ≥ 0)
    bootstrap::Bool                = true
    oob_score::Bool                = false
    n_jobs::Option{Int}            = nothing
    random_state::Any              = nothing
    verbose::Int                   = 0
    warm_start::Bool               = false
    class_weight::Any              = nothing
    ccp_alpha::Float64             =0.0::(_ ≥ 0)
    max_samples::Union{Nothing,Float64,Int} =
        nothing::(_ === nothing || (_ ≥ 0 && (_ isa Integer || _ ≤ 1)))
end
MMI.fitted_params(m::RandomForestClassifier, (f, _, _)) = (
    estimators            = f.estimators_,
    classes               = f.classes_,
    n_classes             = f.n_classes_,
    n_features            = f.n_features_,
    n_outputs             = f.n_outputs_,
    feature_importances   = f.feature_importances_,
    oob_score             = m.oob_score ? f.oob_score_ : nothing,
    oob_decision_function = m.oob_score ? f.oob_decision_function_ : nothing
    )
meta(RandomForestClassifier,
    input   = Table(Count,Continuous),
    target  = AbstractVector{<:Finite},
    weights = false
    )

const ENSEMBLE_REG = Union{Type{<:AdaBoostRegressor}, Type{<:BaggingRegressor}, Type{<:GradientBoostingRegressor}}

MMI.input_scitype(::ENSEMBLE_REG)  = Table(Continuous)
MMI.target_scitype(::ENSEMBLE_REG) = AbstractVector{Continuous}

# ============================================================================
const ExtraTreesRegressor_ = sken(:ExtraTreesRegressor)
@sk_reg mutable struct ExtraTreesRegressor <: MMI.Deterministic
    n_estimators::Int              = 100::(_>0)
    criterion::String              = "mse"::(_ in ("mae", "mse"))
    max_depth::Option{Int}         = nothing::(_ === nothing || _ > 0)
    min_samples_split::Union{Int,Float64}  = 2::(_ > 0)
    min_samples_leaf::Union{Int,Float64}   = 1::(_ > 0)
    min_weight_fraction_leaf::Float64      = 0.0::(_ ≥ 0)
    max_features::Union{Int,Float64,String,Nothing} = "auto"::(_ === nothing || (isa(_, String) && (_ in ("auto","sqrt","log2"))) || _ > 0)
    max_leaf_nodes::Option{Int}    = nothing::(_ === nothing || _ > 0)
    min_impurity_decrease::Float64 = 0.0::(_ ≥ 0)
    bootstrap::Bool                = true
    oob_score::Bool                = false
    n_jobs::Option{Int}            = nothing
    random_state::Any              = nothing
    verbose::Int                   = 0
    warm_start::Bool               = false
end
MMI.fitted_params(m::ExtraTreesRegressor, (f, _, _)) = (
    estimators          = f.estimators_,
    feature_importances = f.feature_importances_,
    n_features          = f.n_features_,
    n_outputs           = f.n_outputs_,
    oob_score           = m.oob_score ? f.oob_score_ : nothing,
    oob_prediction      = m.oob_score ? f.oob_prediction_ : nothing,
    )
meta(ExtraTreesRegressor,
    input   = Table(Continuous),
    target  = AbstractVector{Continuous},
    weights = false
    )

"""
$(MMI.doc_header(ExtraTreesRegressor))

Extra trees regressor, fits a number of randomized decision trees on
various sub-samples of the dataset and uses averaging to improve the
predictive accuracy and control over-fitting.

"""
ExtraTreesRegressor

# ----------------------------------------------------------------------------
const ExtraTreesClassifier_ = sken(:ExtraTreesClassifier)
@sk_clf mutable struct ExtraTreesClassifier <: MMI.Probabilistic
    n_estimators::Int              = 100::(_>0)
    criterion::String              = "gini"::(_ in ("gini", "entropy"))
    max_depth::Option{Int}         = nothing::(_ === nothing || _ > 0)
    min_samples_split::Union{Int,Float64}  = 2::(_ > 0)
    min_samples_leaf::Union{Int,Float64}   = 1::(_ > 0)
    min_weight_fraction_leaf::Float64      = 0.0::(_ ≥ 0)
    max_features::Union{Int,Float64,String,Nothing} = "auto"::(_ === nothing || (isa(_, String) && (_ in ("auto","sqrt","log2"))) || _ > 0)
    max_leaf_nodes::Option{Int}    = nothing::(_ === nothing || _ > 0)
    min_impurity_decrease::Float64 = 0.0::(_ ≥ 0)
    bootstrap::Bool                = true
    oob_score::Bool                = false
    n_jobs::Option{Int}            = nothing
    random_state::Any              = nothing
    verbose::Int                   = 0
    warm_start::Bool               = false
    class_weight::Any              = nothing
end
MMI.fitted_params(m::ExtraTreesClassifier, (f, _, _)) = (
    estimators            = f.estimators_,
    classes               = f.classes_,
    n_classes             = f.n_classes_,
    feature_importances   = f.feature_importances_,
    n_features            = f.n_features_,
    n_outputs             = f.n_outputs_,
    oob_score             = m.oob_score ? f.oob_score_ : nothing,
    oob_decision_function = m.oob_score ? f.oob_decision_function_ : nothing,
    )
meta(ExtraTreesClassifier,
    input   = Table(Continuous),
    target  = AbstractVector{<:Finite},
    weights = false
    )

"""
$(MMI.doc_header(ExtraTreesClassifier))

Extra trees classifier, fits a number of randomized decision trees on
various sub-samples of the dataset and uses averaging to improve the
predictive accuracy and control over-fitting.

"""
ExtraTreesClassifier


text:
function compress(file)
        to = dirname(file)
        dest = joinpath(to, string(basename(file), ".7z"))
        try
                cmd = "7z a -y -t7z -- $(dest) $(file)"
                run(pipeline(ifelse(Sys.iswindows(), `cmd /c $cmd`, `sh -c $cmd`); stdout=joinpath(to, "7zip.log")))
        catch e
                rm(dest; force=true)
                throw(e)
        end
        return dest
end

function decompress(file)
        to = dirname(file)
        source = joinpath(dirname(file), basename(file)[1:end-3])
        try
                cmd = "7z e -y -o$(to) $(file)"
                run(pipeline(ifelse(Sys.iswindows(), `cmd /c $cmd`, `sh -c $cmd`); stdout=joinpath(to, "7zip.log")))
        catch e
                rm(source; force=true)
                throw(e)
        end
        return source
end


text:
## ------------------------------------------------------------------
# const _META_SECTION_KEY = "Meta"
# function new_document!(dockey::String; kwargs...)

#     # book
#     book = currbook()
#     isnothing(book) && error("No Book selected. See `openbook`.")

#     # Add doc
#     dockey = format_idkey(dockey)
#     doc = RBDoc(book, dockey)
#     book[dockey] = doc
#     currdoc!(doc)

#     # Add Meta section
#     sec = add_section!(doc, _META_SECTION_KEY)
#     for (key, value) in kwargs
#         add_pair!(sec, string(key), value)
#     end

#     return doc
# end

# # Use a kwargs as input for the dockey
# function new_document!(dockey::Symbol; kwargs...)
#     !haskey(kwargs, dockey) && error("refkey ':$(dockey)' missing!")
#     dockey = string(kwargs[dockey])
#     new_document!(dockey::String; kwargs...)
# end


text:
module RelSaturate
using Hecke

Hecke.add_verbose_scope(:Saturate)

export saturate!

function dlog(dl::Dict, x, p::Int)
  if iszero(x)
    throw(Hecke.BadPrime(1))
  end
  if haskey(dl, x)
    return dl[x]
  end
#  println("difficult for ", parent(x))
  i = 2
  y = x*x
  while !haskey(dl, y)
    y *= x
    i += 1
    @assert i <= p
  end
  #OK: we know x^i = g^dl[y] (we don't know g)
  v = dl[y]
  g = gcd(p, i)
  r = div(p, g)
  @assert v % g == 0
  e = invmod(div(i, g), r)*div(v, g) % r
  if e == 0
    e = r
  end
  dl[x] = e
  y = x*x
  f = (e*2) % p
  while !isone(y)
    if haskey(dl, y)
      @assert dl[y] == f
    end
    dl[y] = f
    y *= x
    f = (f+e) % p
  end
  g = [ a for (a,b) = dl if b == 1]
  @assert length(g) == 1
  @assert g[1]^dl[x] == x
  return dl[x]
end

function Hecke.matrix(R::Hecke.Ring, M::MatElem)
  return matrix(R, nrows(M), ncols(M), elem_type(R)[R(M[i,j]) for i=1:nrows(M) for j=1:ncols(M)])
end

function mod_p(R, Q::NfOrdIdl, p::Int, T)
  Zk = order(Q)
  F, mF = Hecke.ResidueFieldSmall(Zk, Q)
  mF1 = Hecke.NfToFqNmodMor_easy(mF, number_field(Zk))
  @assert size(F) % p == 1
  pp, e = Hecke.ppio(Int(size(F)-1), p)
  dl = Dict{elem_type(F), Int}()
  dl[F(1)] = 0
  lp = factor(p)
  while true
    x = rand(F)
    if iszero(x)
      continue
    end
    x = x^e
    if any(i-> x^div(pp, Int(i)) == 1, keys(lp.fac))
      continue
    else
      dlog(dl, x, pp)
      @assert length(dl) == pp
      break
    end
  end
  #TODO: in the image of mF1, if the input is a FacElem, the exponents should be reduced by pp.
  #This avoids some inverses.
  return matrix(T, 1, length(R), Int[dlog(dl, image(mF1, x, pp)^e, pp) % p for x = R])
end

Hecke.lift(A::fmpz_mat) = A

#= idea
  let G = <a_1, ... a_n> and p a prime
    for prime ideals Q s.th. N(Q) = 1 (p) we do
        log_Q(a_i)
        nullspace() to get comb. that looks locally like an p-th power
        change a_i
=#

function lift_nonsymmetric(a::nmod_mat)
  z = fmpz_mat(nrows(a), ncols(a))
  z.base_ring = FlintZZ
  ccall((:fmpz_mat_set_nmod_mat_unsigned, :libflint), Nothing,
          (Ref{fmpz_mat}, Ref{nmod_mat}), z, a)
  return z
end

function lift_nonsymmetric(a::gfp_mat)
  z = fmpz_mat(nrows(a), ncols(a))
  z.base_ring = FlintZZ
  ccall((:fmpz_mat_set_nmod_mat_unsigned, :libflint), Nothing,
          (Ref{fmpz_mat}, Ref{gfp_mat}), z, a)
  return z
end

function saturate_exp(c::Hecke.ClassGrpCtx, p::Int, stable = 1.5)
  ZK = order(c.FB.ideals[1])
  T, mT = torsion_unit_group(ZK)
  sT = Int(order(T))

  R = vcat(c.R_gen, c.R_rel)
  K = nf(ZK)
  zeta = mT(T[1])
  if gcd(sT, p) != 1 && !(hash(zeta) in c.RS) # && order is promising...
    push!(R, K(zeta))
#  else
#    println("NOT doint zeta")
  end
  T = GF(p, cached = false)
  A = identity_matrix(T, length(R))
  cA = ncols(A)
  i = 1

  S = Hecke.PrimesSet(Hecke.p_start, -1, Int(p), 1)
  for q in S
    @vprint :Saturate 1 "Finding primes for saturation: $i/$(stable*ncols(A))\n"
    if isdefining_polynomial_nice(K) && isindex_divisor(ZK, q)
      continue
    end
    if discriminant(ZK) % q == 0
      continue
    end
    #if gcd(div(q-1, Int(pp)), pp) > 1 # not possible if cond(k) is involved
    #  continue
    #end
    @vtime :Saturate 3 lq = prime_decomposition(ZK, q, 1)
    if length(lq) == 0
      continue
    end
    for Q in lq
      try
        @vtime :Saturate 1 z = mod_p(R, Q[1], Int(p), T)
        z = z*A
        rrz, z = nullspace(z)
        if iszero(rrz)
          return zero_matrix(FlintZZ, 0, length(R))
        end
        A = A*sub(z, 1:nrows(z), 1:rrz)
        if cA == ncols(A)
          break #the other ideals are going to give the same info
                #for multi-quad as the field is normal
        end
      catch e
        if !isa(e, Hecke.BadPrime)
          rethrow(e)
        end
      end
    end
    if cA == ncols(A)
      i += 1
    else
      i = 0
      cA = ncols(A)
    end
    if i> stable*ncols(A)
      break
    end
  end
  return lift_nonsymmetric(A)
end

fe(a::FacElem) = a
fe(a::nf_elem) = FacElem(a)

function elems_from_sat(c::Hecke.ClassGrpCtx, z)
  res = []#Tuple{FacElem{nf_elem, AnticNumberField}, }[]
  for i=1:ncols(z)
    a = fe(c.R_gen[1])^FlintZZ(z[1, i])
    b = FlintZZ(z[1, i]) * c.M.bas_gens[1]
    for j=2:length(c.R_gen)
      a *= fe(c.R_gen[j])^FlintZZ(z[j, i])
      b += FlintZZ(z[j, i]) * c.M.bas_gens[j]
    end
    for j=1:length(c.R_rel)
      a *= fe(c.R_rel[j])^FlintZZ(z[j + length(c.R_gen), i])
      b += FlintZZ(z[j + length(c.R_gen), i]) * c.M.rel_gens[j]
    end
    push!(res, (a, b))
  end
  return res
end

function saturate!(d::Hecke.ClassGrpCtx, U::Hecke.UnitGrpCtx, n::Int, stable = 3.5)
  @assert isprime(n)
  @vtime :Saturate 1 c = simplify(d, U)
  success = false
  while true
    @vprint :Saturate 1 "Computing candidates for the saturation ...\n"
    @vtime :Saturate 1 e = saturate_exp(c, n, stable)
    if nrows(e) == 0
      @vprint :Saturate 1  "sat yielded nothing new at ", stable, success, "\n"
      return success
    end
    se = sparse_matrix(e)'

    A = sparse_matrix(FlintZZ)
    K = nf(c)
    t, mt = torsion_unit_group(maximal_order(K))
    zeta = K(mt(t[1]))

    @vprint :Saturate 1 "sat: (Hopefully) enlarging by $(ncols(e)) elements\n"

    wasted = false
    for i=1:ncols(e)
      r = e[:, i]
      @assert content(r) == 1
      a = FacElem(K(1))
      fac_a = SRow(FlintZZ)
      for j = 1:length(c.R_gen)
        if !iszero(r[j, 1])
          mul!(a, a, fe(c.R_gen[j])^r[j, 1])
          fac_a += r[j, 1] * c.M.bas_gens[j]
        end
      end
      for j=1:length(c.R_rel)
        if !iszero(r[j + length(c.R_gen), 1])
          mul!(a, a, fe(c.R_rel[j])^r[j + length(c.R_gen), 1])
          fac_a += r[j + length(c.R_gen), 1] * c.M.rel_gens[j]
        end
      end
      if nrows(e) > length(c.R_gen) + length(c.R_rel) && !iszero(r[nrows(e), 1])
        @assert length(c.R_gen) + length(c.R_rel) + 1 == nrows(e)
        a *= fe(zeta)^r[nrows(e), 1]
      end

      decom = Dict((c.FB.ideals[k], v) for (k,v) = fac_a)

      #for (k, v) in a.fac
      #  if iszero(v)
      #    delete!(a.fac, k)
      #  end
      #end

      @vtime :Saturate 1 fl, x = ispower(a, n, decom = decom)

      if fl
        @assert isa(x, FacElem)
        success = true
        fac_a = divexact(fac_a, n)
        Hecke.class_group_add_relation(d, x, fac_a)
        @vprint :Saturate 1  "sat added new relation\n"
        if iszero(fac_a) #to make sure the new unit is used!
          #find units can be randomised...
          #maybe that should also be addressed elsewhere
          @vprint :Saturate 2  "sat added new unit\n"
          Hecke._add_dependent_unit(U, x)
        end
      else
        @vprint :Saturate 2  "sat wasted time, local power wasn't global\n"
        wasted = true
      end
    end
    if wasted
      stable *= 2
    else
      @vprint :Saturate  1 "sat success at ", stable, "\n"
      return success
    end
  end
end

function simplify(c::Hecke.ClassGrpCtx, U::Hecke.UnitGrpCtx)
  d = Hecke.class_group_init(c.FB, SMat{fmpz}, add_rels = false)

  Hecke.module_trafo_assure(c.M)
  trafos = c.M.trafo

  for i=1:length(c.FB.ideals)
    c.M.basis.rows[i].values[1] == 1 && continue
    @assert all(x -> x > 0, c.M.basis.rows[i].values)
    x = zeros(fmpz, length(c.R_gen) + length(c.R_rel))
    x[i] = 1
    for j in length(trafos):-1:1
      Hecke.apply_right!(x, trafos[j])
    end
    y = FacElem(vcat(c.R_gen, c.R_rel), x)
    fl = Hecke.class_group_add_relation(d, y, deepcopy(c.M.basis.rows[i]))
    @assert fl
  end
  for i=1:length(U.units)
    Hecke.class_group_add_relation(d, U.units[i], SRow(FlintZZ))
  end
  return d
end

end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct RoutePropertiesFormat2 <: SwaggerModel
    addressPrefix::Any # spec type: Union{ Nothing, String } # spec name: addressPrefix
    nextHopType::Any # spec type: Union{ Nothing, RouteNextHopType2 } # spec name: nextHopType
    nextHopIpAddress::Any # spec type: Union{ Nothing, String } # spec name: nextHopIpAddress
    provisioningState::Any # spec type: Union{ Nothing, ProvisioningState } # spec name: provisioningState

    function RoutePropertiesFormat2(;addressPrefix=nothing, nextHopType=nothing, nextHopIpAddress=nothing, provisioningState=nothing)
        o = new()
        validate_property(RoutePropertiesFormat2, Symbol("addressPrefix"), addressPrefix)
        setfield!(o, Symbol("addressPrefix"), addressPrefix)
        validate_property(RoutePropertiesFormat2, Symbol("nextHopType"), nextHopType)
        setfield!(o, Symbol("nextHopType"), nextHopType)
        validate_property(RoutePropertiesFormat2, Symbol("nextHopIpAddress"), nextHopIpAddress)
        setfield!(o, Symbol("nextHopIpAddress"), nextHopIpAddress)
        validate_property(RoutePropertiesFormat2, Symbol("provisioningState"), provisioningState)
        setfield!(o, Symbol("provisioningState"), provisioningState)
        o
    end
end # type RoutePropertiesFormat2

const _property_map_RoutePropertiesFormat2 = Dict{Symbol,Symbol}(Symbol("addressPrefix")=>Symbol("addressPrefix"), Symbol("nextHopType")=>Symbol("nextHopType"), Symbol("nextHopIpAddress")=>Symbol("nextHopIpAddress"), Symbol("provisioningState")=>Symbol("provisioningState"))
const _property_types_RoutePropertiesFormat2 = Dict{Symbol,String}(Symbol("addressPrefix")=>"String", Symbol("nextHopType")=>"RouteNextHopType2", Symbol("nextHopIpAddress")=>"String", Symbol("provisioningState")=>"ProvisioningState")
Base.propertynames(::Type{ RoutePropertiesFormat2 }) = collect(keys(_property_map_RoutePropertiesFormat2))
Swagger.property_type(::Type{ RoutePropertiesFormat2 }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_RoutePropertiesFormat2[name]))}
Swagger.field_name(::Type{ RoutePropertiesFormat2 }, property_name::Symbol) =  _property_map_RoutePropertiesFormat2[property_name]

function check_required(o::RoutePropertiesFormat2)
    (getproperty(o, Symbol("nextHopType")) === nothing) && (return false)
    true
end

function validate_property(::Type{ RoutePropertiesFormat2 }, name::Symbol, val)
end


text:
# # Clocked thermostat

# ## Model

module clocked_thermostat #jl

using ReachabilityAnalysis, ModelingToolkit

const A1 = 40.0
const A2 = 30.0
const B = 0.5
const t1 = 75
const t2 = 65

const var′ = @variables t t′

function thermostat_on′()
    invariant = HalfSpace(t <= t1, var′)
    @system(x' = [-B 0; 0 0]*x + [A1, 1], x ∈ invariant)
end

function thermostat_off′()
    invariant = HalfSpace(t >= t2, var′)
    @system(x' = [-B 0; 0 0]*x + [A2, 1], x ∈ invariant)
end

function thermostat_hybrid′()
    automaton = LightAutomaton(2)
    add_transition!(automaton, 1, 2, 1)
    add_transition!(automaton, 2, 1, 2)

    mode1 = thermostat_on′()
    mode2 = thermostat_off′()
    modes = [mode1, mode2]

    ## transition on -> off
    guard = HalfSpace(t >= t1, var′)
    trans1 = ConstrainedIdentityMap(1, guard)
    ## transition off -> on
    guard = HalfSpace(t <= t2, var′)
    trans2 = ConstrainedIdentityMap(1, guard)
    resetmaps = [trans1, trans2]

    return HybridSystem(automaton, modes, resetmaps, [AutonomousSwitching()])
end

function model(X0)
    H = thermostat_hybrid′()
return IVP(H, X0)
end

end  #jl

# ## References

max_stars_repo_name: mforets/ReachabilityLibrary

max_issues_repo_name: mforets/ReachabilityLibrary

max_forks_repo_name: mforets/ReachabilityLibrary