text:
# refinements can store parent informations
abstract type CellParents <: AbstractGridIntegerArray1D end

function ExtendableGrids.instantiate(xgrid::ExtendableGrid, ::Type{CellParents})
    ncells = num_sources(xgrid[CellNodes])
    VectorOfConstants{K}(0,ncells)
end

# functions that tell how to split one ElementGeometry into another
const _split_rule_Triangle2D = reshape([1,2,3],3,1)
const _split_rule_Tetrahedron3D = reshape([1,2,3,4],4,1)
const _split_rule_Quadrilateral2D_Triangle2D = [1 2 3;1 3 4]'
const _split_rule_Hexahedron3D_Tetrahedron3D = [1 2 4 8; 2 3 4 8; 2 7 3 8; 6 8 7 2; 6 5 8 2; 8 5 1 2]'

split_refine_rule(::Type{Triangle2D}, ::Type{Triangle2D}) = _split_rule_Triangle2D
split_refine_rule(::Type{Tetrahedron3D}, ::Type{Tetrahedron3D}) = _split_rule_Tetrahedron3D
split_refine_rule(::Type{<:Quadrilateral2D}, ::Type{Triangle2D}) = _split_rule_Quadrilateral2D_Triangle2D
split_refine_rule(::Type{<:Hexahedron3D}, ::Type{Tetrahedron3D}) = _split_rule_Hexahedron3D_Tetrahedron3D
#[1 2 3 7; 1 3 4 7; 1 5 6 7; 1 8 5 7; 1 6 2 7; 1 4 8 7]

"""
$(TYPEDSIGNATURES)

generates a new ExtendableGrid by splitting each cell into subcells of the specified targetgeometry

split rules exist for
- Quadrilateral2D into Triangle2D
- Hexahedron3D into Tetrahedron3D
"""
function split_grid_into(source_grid::ExtendableGrid{T,K}, targetgeometry::Type{<:AbstractElementGeometry}) where {T,K}
    xgrid=ExtendableGrid{T,K}()
    xgrid[Coordinates]=source_grid[Coordinates]
    oldCellGeometries = source_grid[CellGeometries]
    oldCellRegions = source_grid[CellRegions]
    EG = Base.unique(oldCellGeometries)

    split_rules = Array{Array{Int,2},1}(undef,length(EG))
    for j = 1 : length(EG)
        split_rules[j] = split_refine_rule(EG[j],targetgeometry)
    end
    singleEG::Bool = false
    if length(EG) == 1
        singleEG = true
    end
    xCellNodes = zeros(K,0)
    xCellRegions = zeros(K,0)
    oldCellNodes=source_grid[CellNodes]
    nnodes4item::Int = num_targets(oldCellNodes,1)
    ncells::Int = 0
    itemEG = EG[1]
    iEG::Int = 1
    split_rule::Array{Int,2} = split_rules[iEG]
    for cell = 1 : num_sources(oldCellNodes)
        if !singleEG
            nnodes4item = num_targets(oldCellNodes,cell)
            itemEG = oldCellGeometries[cell]
            iEG = findfirst(isequal(itemEG), EG)
            split_rule = split_rules[iEG]
        end
        for j = 1 : size(split_rule,2), k = 1 : size(split_rule,1)
            append!(xCellNodes,oldCellNodes[split_rule[k,j],cell])
        end
        for j = 1 : size(split_rule,2)
            push!(xCellRegions,oldCellRegions[cell])
        end
        ncells += size(split_rule,2)
    end
    xCellNodes = reshape(xCellNodes,num_nodes(targetgeometry),ncells)
    xgrid[CellNodes] = xCellNodes
    xgrid[CellGeometries] = VectorOfConstants(targetgeometry,ncells)
    xgrid[CoordinateSystem]=source_grid[CoordinateSystem]
    if typeof(oldCellRegions) <: VectorOfConstants
        xgrid[CellRegions] = VectorOfConstants{K}(1,ncells)
    else
        xgrid[CellRegions] = xCellRegions
    end

    # find new boundary faces (easy in 2D, not so easy in 3D)
    if dim_element(targetgeometry) == 2 # BFaces are Edge1D wich stay the same
        xgrid[BFaceNodes]=source_grid[BFaceNodes]
        xgrid[BFaceRegions]=source_grid[BFaceRegions]
        xgrid[BFaceGeometries]=VectorOfConstants(facetype_of_cellface(targetgeometry,1),num_sources(xgrid[BFaceNodes]))
    elseif dim_element(targetgeometry) == 3
        # BFaces may be split into different shapes, e.g. from Quadrilateral2D to two Triangle2D
        # and it is hard to predict how they are splitted
        # so we do something lazy here and search for new faces that lie in old bfaces
        oldBFaceNodes = source_grid[BFaceNodes]
        oldBFaceRegions = source_grid[BFaceRegions]
        newFaceNodes = xgrid[FaceNodes]
        nfaces = num_sources(newFaceNodes)
        nbfaces = num_sources(oldBFaceNodes)

        newBFaceNodes = zeros(K,0)
        newBFaceRegions = zeros(K,0)
        newnbfaces = 0
        nnodes = size(xgrid[Coordinates],2)
        flag4item = zeros(Bool,nnodes)
        nodes_per_bface::Int = 0
        nodes_per_face::Int = 0
        common_nodes::Int = 0
        for bface = 1 : nbfaces
            # flag nodes of old bface
            nodes_per_bface = num_targets(oldBFaceNodes,bface)
            for j = 1 : nodes_per_bface
                flag4item[oldBFaceNodes[j,bface]] = true
            end

            # find matching faces
            for face = 1 : nfaces
                nodes_per_face = num_targets(newFaceNodes,face)
                common_nodes = 0
                for k = 1 : nodes_per_face
                    if flag4item[newFaceNodes[k,face]] == true
                        common_nodes += 1
                    else
                        break
                    end
                end
                if common_nodes == nodes_per_face
                    append!(newBFaceNodes,newFaceNodes[:,face])
                    push!(newBFaceRegions,oldBFaceRegions[bface])
                    newnbfaces += 1
                end
            end

            # reset flags
            for j = 1 : nodes_per_bface
                flag4item[oldBFaceNodes[j,bface]] = false
            end
        end

        newBFaceNodes = reshape(newBFaceNodes,num_nodes(facetype_of_cellface(targetgeometry,1)),newnbfaces)
        xgrid[BFaceNodes]=newBFaceNodes
        xgrid[BFaceRegions]=newBFaceRegions
        xgrid[BFaceGeometries]=VectorOfConstants(facetype_of_cellface(targetgeometry,1),newnbfaces)
    end

    return xgrid
end

uniform_refine_needfacemidpoints(::Type{<:AbstractElementGeometry}) = true
uniform_refine_needcellmidpoints(::Type{<:AbstractElementGeometry}) = false

# uniform refinement rules in 1D
# first k nodes are the CellNodes
# next node is the CellMidpoint
const _uniform_rule_Edge1D = [1 3; 3 2]'
const _uniform_rule_Triangle2D = [1 4 6; 4 2 5; 6 5 3; 5 6 4]'
const _uniform_rule_Quadrilateral2D = [1 5 9 8; 2 6 9 5; 3 7 9 6; 4 8 9 7]'
const _uniform_rule_Tetrahedron3D = [   1 5 6 7;
                                        2 5 9 8;
                                        3 10 6 8;
                                        4 10 9 7;
                                        10 5 8 9;
                                        5 10 7 9;
                                        5 10 8 6;
                                        10 5 7 6]'
const _uniform_rule_Hexahedron3D = [     1   9  21  12  13  22  27  25;
                                         9   2  10  21  22  14  23  27;
                                        12  21  11   4  25  27  24  16;
                                        21  10   3  11  27  23  15  24;
                                        13  22  27  25   5  17  26  20;
                                        22  14  23  27  17   6  18  26;
                                        25  27  24  16  20  26  19   8;
                                        27  23  15  24  26  18   7  19]'

uniform_refine_rule(::Type{<:Edge1D}) = _uniform_rule_Edge1D
uniform_refine_needcellmidpoints(::Type{<:Edge1D}) = true

# uniform refinement rules in 2D
# first k nodes are the CellNodes
# next m nodes are the CellFaces midpoints
# next node is the CellMidpoint (if needed)
uniform_refine_rule(::Type{<:Triangle2D}) = _uniform_rule_Triangle2D
uniform_refine_rule(::Type{<:Quadrilateral2D}) = _uniform_rule_Quadrilateral2D
uniform_refine_needcellmidpoints(::Type{<:Quadrilateral2D}) = true

# uniform refinement rules in 3D
# first k nodes are the CellNodes
# next m nodes are the CellEdges midpoints
# next n nodes are the CellFaces midpoints
# next node is the CellMidpoint (if needed)
uniform_refine_rule(::Type{<:Tetrahedron3D}) = _uniform_rule_Tetrahedron3D
uniform_refine_needfacemidpoints(::Type{<:Tetrahedron3D}) = false
uniform_refine_rule(::Type{<:Hexahedron3D}) = _uniform_rule_Hexahedron3D
uniform_refine_needcellmidpoints(::Type{<:Hexahedron3D}) = true

"""
$(TYPEDSIGNATURES)

generates a new ExtendableGrid by uniform refinement of each cell in the given grid

uniform refinement rules are available for these AbstractElementGeometries:
- Line1D (bisection into two subsegments)
- Triangle2D (red refinement into four subtriangles)
- Quadrilateral2D (into four subquadrilaterals)
- Tetrahedron (into eight subtetrahedrons)
- Hexahedron (into eight subhexahedrons)

if multiple geometries are in the mesh uniform refinement will only work
if all refinement rules refine faces and edges (in 3D) equally
(so no hanging nodes are created)
"""
function uniform_refine(source_grid::ExtendableGrid{T,K}; store_parents = false) where {T,K}
    # @logmsg MoreInfo "Uniform refinement of $(num_sources(source_grid[CellNodes])) cells"

    xgrid = ExtendableGrid{T,K}()
    xgrid[CoordinateSystem]=source_grid[CoordinateSystem]

    # unpack stuff from source grid
    oldCoordinates = source_grid[Coordinates]
    oldCellGeometries = source_grid[CellGeometries]
    oldCellRegions = source_grid[CellRegions]
    EG = Base.unique(oldCellGeometries)

    singleEG = false
    if length(EG) == 1
        singleEG = true
    end

    # get dimension of CellGeometries
    # currently it is assumed to be the same for all cells
    dim::Int = dim_element(EG[1])

    refine_rules::Array{Array{Int,2},1} = Array{Array{Int,2},1}(undef,length(EG))
    need_facemidpoints = uniform_refine_needfacemidpoints(EG[1])
    for j = 1 : length(EG)
        refine_rules[j] = uniform_refine_rule(EG[j])
        @assert uniform_refine_needfacemidpoints(EG[j]) == need_facemidpoints
    end

    xCellNodes = VariableTargetAdjacency(K)
    xCellGeometries = Array{DataType,1}(undef,0)
    xCellRegions = zeros(K,0)
    oldCellNodes::Adjacency{K} = source_grid[CellNodes]
    oldCellFaces::Adjacency{K} = source_grid[CellFaces]
    nfaces::Int = 0
    if dim > 1
        oldFaceNodes = source_grid[FaceNodes]
        nfaces = num_sources(oldFaceNodes)
    end
    nedges::Int = 0
    if dim > 2
        oldEdgeNodes::Adjacency{K} = source_grid[EdgeNodes]
        nedges = num_sources(oldEdgeNodes)
        oldCellEdges::Adjacency{K} = source_grid[CellEdges]
    end

    # determine number of new vertices
    itemEG = Triangle2D
    newvertices::Int = 0 # in 1D no additional vertices on the faces are needed
    if dim == 2 # in 2D each face is halved
        newvertices = nfaces
    elseif dim == 3 # in 2D each face and edge is halved
        newvertices = nedges
        if need_facemidpoints
            newvertices += nfaces
        end
    end
    oldvertices::Int = size(oldCoordinates,2)
    newnode::Int = oldvertices + newvertices
    # additionally cell midpoints are needed for some refinements
    for cell = 1 : num_sources(oldCellNodes)
        itemEG = oldCellGeometries[cell]
        if uniform_refine_needcellmidpoints(itemEG) == true
            newvertices += 1
        end
    end
    xCoordinates::Array{T,2} = zeros(T,size(oldCoordinates,1),oldvertices+newvertices)
    @views xCoordinates[:,1:oldvertices] = oldCoordinates

    newvertex::Array{T,1} = zeros(T,size(xCoordinates,1))
    nnodes4item::Int = 0
    if dim > 2 # add edge midpoints to Coordinates
        for edge = 1 : nedges
            nnodes4item = num_targets(oldEdgeNodes,edge)
            fill!(newvertex,0.0)
            for k = 1 : nnodes4item, d = 1 : size(xCoordinates,1)
                newvertex[d] += xCoordinates[d,oldEdgeNodes[k,edge]]
            end
            newvertex ./= nnodes4item
            for d = 1 : size(xCoordinates,1)
                xCoordinates[d,oldvertices+edge] = newvertex[d]
            end
        end
    end
    if dim > 1 && need_facemidpoints
        # add face midpoints to Coordinates
        for face = 1 : nfaces
            nnodes4item = num_targets(oldFaceNodes,face)
            fill!(newvertex,0.0)
            for k = 1 : nnodes4item, d = 1 : size(xCoordinates,1)
                newvertex[d] += xCoordinates[d,oldFaceNodes[k,face]]
            end
            newvertex ./= nnodes4item
            for d = 1 : size(xCoordinates,1)
                xCoordinates[d,oldvertices+nedges+face] = newvertex[d]
            end
        end
    end

    # determine new cells
    nnodes4item = num_nodes(itemEG)
    nfaces4item::Int = num_faces(itemEG)
    nedges4item::Int = num_edges(itemEG)
    ncells::Int = 0
    iEG::Int = 1
    subitemnodes::Array{K,1} = zeros(K,max_num_targets_per_source(oldCellNodes)+max_num_targets_per_source(oldCellFaces)+1)
    m::Int = 0
    xCellParents::Array{K,1} = zeros(K,0)
    refine_rule::Array{Int,2} = refine_rules[iEG]
    for cell = 1 : num_sources(oldCellNodes)
        if !singleEG
            itemEG = oldCellGeometries[cell]
            nnodes4item = num_nodes(itemEG)
            nfaces4item = num_faces(itemEG)
            nedges4item = num_edges(itemEG)
            iEG = findfirst(isequal(itemEG), EG)
            refine_rule = refine_rules[iEG]
        end
        if uniform_refine_needcellmidpoints(itemEG) == true
            # add cell midpoint to Coordinates
            newnode += 1
            fill!(newvertex,0.0)
            for k = 1 : nnodes4item, d = 1 : size(xCoordinates,1)
                newvertex[d] += xCoordinates[d,oldCellNodes[k,cell]]
            end
            newvertex ./= nnodes4item
            for d = 1 : size(xCoordinates,1)
                xCoordinates[d,newnode] = newvertex[d]
            end
        end
        for j = 1 : size(refine_rule,2)
            for k = 1 : size(refine_rule,1)
                m = refine_rule[k,j]
                if dim == 1
                    if m <= nnodes4item
                        subitemnodes[k] = oldCellNodes[m,cell]
                    else
                        subitemnodes[k] = newnode
                    end
                elseif dim == 2
                    if m <= nnodes4item
                        subitemnodes[k] = oldCellNodes[m,cell]
                    elseif m <= nnodes4item + nfaces4item
                        subitemnodes[k] = oldvertices + oldCellFaces[m-nnodes4item,cell]
                    else
                        subitemnodes[k] = newnode
                    end
                elseif dim == 3
                    if m <= nnodes4item
                        subitemnodes[k] = oldCellNodes[m,cell]
                    elseif m <= nnodes4item + nedges4item
                        subitemnodes[k] = oldvertices + oldCellEdges[m-nnodes4item,cell]
                    elseif m <= nnodes4item + nedges4item + nfaces4item
                        subitemnodes[k] = oldvertices + nedges + oldCellFaces[m-nnodes4item-nedges4item,cell]
                    else
                        subitemnodes[k] = newnode
                    end
                end
            end
            append!(xCellNodes,view(subitemnodes,1:size(refine_rule,1)))
            push!(xCellRegions, oldCellRegions[cell])
            if !singleEG
                push!(xCellGeometries,itemEG)
            end
            if store_parents
                push!(xCellParents,cell)
            end
        end
        ncells += size(refine_rule,2)
    end

    # assign new cells to grid
    xgrid[Coordinates] = xCoordinates

    if singleEG
        nnodes4item = size(oldCellNodes,1)
        xgrid[CellNodes] = reshape(xCellNodes.colentries,nnodes4item,num_sources(xCellNodes))
        xgrid[CellGeometries] = VectorOfConstants(EG[1],ncells)
    else
        xgrid[CellNodes] = xCellNodes
        xgrid[CellGeometries] = xCellGeometries
    end
    if typeof(oldCellRegions) <: VectorOfConstants
        xgrid[CellRegions] = VectorOfConstants{K}(1,ncells)
    else
        xgrid[CellRegions] = xCellRegions
    end

    # determine new boundary faces
    oldBFaceNodes::Adjacency{K} = source_grid[BFaceNodes]
    oldBFaceFaces::Array{K,1} = source_grid[BFaceFaces]
    oldBFaceRegions = source_grid[BFaceRegions]
    oldBFaceGeometries = source_grid[BFaceGeometries]

    if dim == 1
        xgrid[BFaceNodes] = oldBFaceNodes
        xgrid[BFaceRegions] = oldBFaceRegions
        xgrid[BFaceGeometries] = oldBFaceGeometries
    else
        xBFaceRegions = zeros(K,0)
        xBFaceGeometries = Array{DataType,1}(undef,0)
        nbfaces = num_sources(oldBFaceNodes)
        if dim == 2 || typeof(oldBFaceNodes) == Array{K,2}
            xBFaceNodes = zeros(K,0)
        else
            xBFaceNodes = VariableTargetAdjacency(K)
        end
        if dim == 3
            xCellEdges = source_grid[CellEdges]
            xNodeEdges::Adjacency{K} = atranspose(oldEdgeNodes)
        else
            xCellEdges = []
        end

        EG = Base.unique(oldBFaceGeometries)

        refine_rules = Array{Array{Int,2},1}(undef,length(EG))
        for j = 1 : length(EG)
            refine_rules[j] = uniform_refine_rule(EG[j])
        end

        face::K = 0
        edge::Int = 0
        n1::Int = 0
        ne4n1::Int = 0
        n2::Int = 0
        ne4n2::Int = 0
        newnbfaces::Int = 0
        bface_enum_rule::Array{Int,2} = local_cellfacenodes(EG[1])
        for bface = 1 : nbfaces
            face = oldBFaceFaces[bface]
            itemEG = oldBFaceGeometries[bface]
            nnodes4item = num_nodes(itemEG)
            nfaces4item = num_faces(itemEG)
            nedges4item = num_edges(itemEG)
            iEG = findfirst(isequal(itemEG), EG)
            bface_enum_rule = local_cellfacenodes(itemEG)

            for j = 1 : size(refine_rules[iEG],2)
                for k = 1 : size(refine_rules[iEG],1)
                    m = refine_rules[iEG][k,j]
                    if dim == 2
                        if m <= nnodes4item
                            subitemnodes[k] = oldBFaceNodes[m,bface]
                        else
                            subitemnodes[k] = oldvertices + face
                        end
                    elseif dim == 3
                        if m <= nnodes4item
                            subitemnodes[k] = oldBFaceNodes[m,bface]
                        elseif m <= nnodes4item + nfaces4item
                            edge = m-nnodes4item # local number
                            # find global edge number
                            n1 = oldBFaceNodes[bface_enum_rule[1,edge],bface]
                            n2 = oldBFaceNodes[bface_enum_rule[2,edge],bface]
                            ne4n1 = num_targets(xNodeEdges,n1)
                            ne4n2 = num_targets(xNodeEdges,n2)
                            for ge = 1 : ne4n1, gf = 1 : ne4n2
                                if xNodeEdges[ge,n1] == xNodeEdges[gf,n2]
                                    edge = xNodeEdges[gf,n2]
                                    break
                                end
                            end
                         #   edge = intersect(xNodeEdges[:,oldBFaceNodes[bface_enum_rule[1,edge],bface]],xNodeEdges[:,oldBFaceNodes[bface_enum_rule[2,edge],bface]])[1]
                            subitemnodes[k] = oldvertices + edge
                        else
                            subitemnodes[k] = oldvertices + nedges + face
                        end
                    end
                end
                append!(xBFaceNodes,view(subitemnodes,1:size(refine_rules[iEG],1)))
                push!(xBFaceGeometries,itemEG)
                push!(xBFaceRegions,oldBFaceRegions[bface])
                newnbfaces += 1
            end
        end
        if dim == 2 || typeof(oldBFaceNodes) == Array{K,2}
            xgrid[BFaceNodes] = reshape(xBFaceNodes,(size(oldBFaceNodes,1),newnbfaces))
        else
            xgrid[BFaceNodes] = xBFaceNodes
        end
        xgrid[BFaceRegions] = xBFaceRegions
        xgrid[BFaceGeometries] = xBFaceGeometries
    end

    if store_parents
        xgrid[CellParents] = xCellParents
    end

    return xgrid
end

function uniform_refine(source_grid::ExtendableGrid{T,K}, nrefinements::Int; store_parents = false) where {T,K}
    xgrid = source_grid
    parents = 1:num_sources(xgrid[CellNodes])
    for j=1:nrefinements
        xgrid = uniform_refine(xgrid; store_parents = store_parents)
        if store_parents
            parents = parents[xgrid[CellParents]]
        end
    end
    if store_parents
        xgrid[CellParents] = parents
    end
    return xgrid
end

# barycentric refinement rules
# first k nodes are the CellNodes, k+1-th  node is cell midpoint
const _barycentric_rule_Triangle2D = [1 2 4; 2 3 4; 3 1 4]'
barycentric_refine_rule(::Type{<:Triangle2D}) = _barycentric_rule_Triangle2D

"""
$(TYPEDSIGNATURES)

generates a new ExtendableGrid by barycentric refinement of each cell in the source grid

barycentric refinement is available for these ElementGeometries
- Quadrilateral2D (first split into Triangle2D)
- Triangle2D
"""
function barycentric_refine(source_grid::ExtendableGrid{T,K}) where {T,K}
    # @logmsg MoreInfo "Barycentric refinement of $(num_sources(source_grid[CellNodes])) cells"

    # split first into triangles
    source_grid = split_grid_into(source_grid,Triangle2D)

    xgrid = ExtendableGrid{T,K}()
    oldCoordinates = source_grid[Coordinates]
    oldCellGeometries = source_grid[CellGeometries]
    oldCellRegions = source_grid[CellRegions]
    EG = Base.unique(oldCellGeometries)

    refine_rules = Array{Array{Int,2},1}(undef,length(EG))
    for j = 1 : length(EG)
        refine_rules[j] = barycentric_refine_rule(EG[j])
    end

    singleEG = false
    if length(EG) == 1
        singleEG = true
    end
    xCellNodes::Adjacency = VariableTargetAdjacency(K)
    xCellGeometries = Array{DataType,1}(undef,0)
    xCellRegions = zeros(K,0)

    oldCellNodes = source_grid[CellNodes]
    oldCellFaces = source_grid[CellFaces]

    # determine number of new vertices
    itemEG = Triangle2D
    newvertices = num_sources(oldCellNodes)
    oldvertices = size(oldCoordinates,2)
    xCoordinates = zeros(T,size(oldCoordinates,1),oldvertices+newvertices)
    @views xCoordinates[:,1:oldvertices] = oldCoordinates

    # determine new cells
    nnodes4item::Int = num_targets(oldCellNodes,1)
    ncells::Int = 0
    iEG::Int = 1
    subitemnodes::Array{K,1} = zeros(K,max_num_targets_per_source(oldCellNodes)+max_num_targets_per_source(oldCellFaces)+1)
    newnode::Int = oldvertices
    m::Int = 0
    newvertex::Array{T,1} = zeros(T,size(xCoordinates,1))
    refine_rule::Array{Int,2} = refine_rules[iEG]
    for cell = 1 : num_sources(oldCellNodes)
        if !singleEG
            nnodes4item = num_targets(oldCellNodes,cell)
            itemEG = oldCellGeometries[cell]
            iEG = findfirst(isequal(itemEG), EG)
            refine_rule = refine_rules[iEG]
        end

        # add cell midpoint to Coordinates
        newnode += 1
        fill!(newvertex,0.0)
        for k = 1 : nnodes4item, d = 1 : size(xCoordinates,1)
            newvertex[d] += xCoordinates[d,oldCellNodes[k,cell]]
        end
        newvertex ./= nnodes4item
        for d = 1 : size(xCoordinates,1)
            xCoordinates[d,newnode] = newvertex[d]
        end

        for j = 1 : size(refine_rule,2)
            for k = 1 : size(refine_rule,1)
                m = refine_rules[iEG][k,j]
                if m <= nnodes4item
                    subitemnodes[k] = oldCellNodes[m,cell]
                else
                    subitemnodes[k] = newnode
                end
            end
            append!(xCellNodes,view(subitemnodes,1:size(refine_rule,1)))
            if !singleEG
                push!(xCellGeometries,itemEG)
            end
            push!(xCellRegions, oldCellRegions[cell])
        end
        ncells += size(refine_rule,2)
    end

    xgrid[Coordinates] = xCoordinates
    if singleEG
        nnodes4item = size(oldCellNodes,1)
        xgrid[CellNodes] = reshape(xCellNodes.colentries,nnodes4item,num_sources(xCellNodes))
        xgrid[CellGeometries] = VectorOfConstants(EG[1],num_sources(xCellNodes))
    else
        xgrid[CellNodes] = xCellNodes
        xgrid[CellGeometries] = xCellGeometries
    end
    if typeof(oldCellRegions) <: VectorOfConstants
        xgrid[CellRegions] = VectorOfConstants{K}(1,num_sources(xCellNodes))
    else
        xgrid[CellRegions] = xCellRegions
    end
    xgrid[BFaceNodes]=source_grid[BFaceNodes]
    xgrid[BFaceRegions]=source_grid[BFaceRegions]
    xgrid[BFaceGeometries]=source_grid[BFaceGeometries]
    xgrid[CoordinateSystem]=source_grid[CoordinateSystem]
    return xgrid
end


text:
module CalliopeJuMP

export build_julia_model, model_dict

# Import external modules
using NCDatasets; using JuMP; using AxisArrays; using YAML;

include("Util.jl")
include("Variables.jl")
include("Objective.jl")

module load_constraints

    export
        load_capacity_constraints, load_costs_constraints,
        load_dispatch_constraints, load_energy_balance_constraints,
        load_milp_constraints, load_network_constraints,
        load_conversion_constraints, load_conversion_plus_constraints,
        load_export_constraints

    include(joinpath("constraints", "capacity.jl"))
    include(joinpath("constraints", "costs.jl"))
    include(joinpath("constraints", "dispatch.jl"))
    include(joinpath("constraints", "energy_balance.jl"))
    include(joinpath("constraints", "milp.jl"))
    include(joinpath("constraints", "network.jl"))
    include(joinpath("constraints", "conversion.jl"))
    include(joinpath("constraints", "conversion_plus.jl"))
    include(joinpath("constraints", "export.jl"))

end

function build_julia_model(path_to_dataset)

    # Bring in Dataset
    dataset = NCDatasets.Dataset(path_to_dataset);

    #
    ## Sets
    #
    sets = Dict()
    parameters = Dict()
    dimensions = [NCDatasets.nc_inq_dimname(dataset.ncid, i)
                  for i in NCDatasets.nc_inq_dimids(dataset.ncid, false)]

    for var in keys(dataset)
        show(var)
        if var in dimensions
            if var == "timesteps"
                show(dataset[var][1])
            end
            sets[var] = dataset[var][:]
        else
            parameters[var] = get_variable(dataset, var)
        end
    end

    # Create JuMP model, including assigning a solver

    run_solver = dataset.attrib["run.solver"]
    if run_solver == "cplex"
        using CPLEX
        solver = CplexSolver
    elseif run_solver == "gurobi"
        using Gurobi
        solver = GurobiSolver
    elseif run_solver == "glpk"
        using GLPKMathProgInterface
        if haskey(sets, "loc_techs_milp") || haskey(sets, "loc_techs_purchase")
            solver = GLPKSolverMIP
        else
            solver = GLPKSolverLP
        end
    else
        throw("unknown solver", run_solver)
    end
    backend_model = JuMP.Model(solver = solver());

    model_dict = Dict("backend_model"=>backend_model,
                      "parameters"=>parameters,
                      "defaults"=>YAML.load(dataset.attrib["defaults"]),
                      "sets"=>sets,
                      "dataset"=>dataset)

    #
    ## Variables
    #
    model_dict["variables"] = initialize_decision_variables(model_dict);

    #
    ## Expressions
    #
    model_dict["expressions"] = Dict()

    #
    ## Constraints
    #

    model_dict["constraints"] = merge(
        load_constraints.load_capacity_constraints(model_dict),
        load_constraints.load_cost_constraints(model_dict),
        load_constraints.load_dispatch_constraints(model_dict),
        load_constraints.load_energy_balance_constraints(model_dict),
        load_constraints.load_network_constraints(model_dict)
    )

    if haskey(sets, "loc_techs_conversion")
        model_dict["constraints"] = merge(
            model_dict["constraints"],
            load_constraints.load_conversion_constraints(model_dict)
        )
    end

    if haskey(sets, "loc_techs_conversion_plus")
        model_dict["constraints"] = merge(
            model_dict["constraints"],
            load_constraints.load_conversion_plus_constraints(model_dict)
        )
    end

    if haskey(sets, "loc_techs_milp") || haskey(sets, "loc_techs_purchase")
        model_dict["constraints"] = merge(
            model_dict["constraints"],
            load_constraints.load_milp_constraints(model_dict)
        )
    end

    if haskey(sets, "loc_techs_export")
        model_dict["constraints"] = merge(
            model_dict["constraints"],
            load_constraints.load_export_constraints(model_dict)
        )
    end

    #
    ## Objective
    #
    cost_minimization(model_dict)

    solve(model_dict["backend_model"])
    return model_dict

end

end

ext: jl

max_stars_repo_name: calliope-project/CalliopeJuMP

max_issues_repo_name: calliope-project/CalliopeJuMP

max_forks_repo_name: calliope-project/CalliopeJuMP


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

using Markdown
using InteractiveUtils

# ╔═╡ 804078c8-3aa1-11eb-2f4d-69380de7ab92
begin
	function ingredients(path::String)
	# this is from the Julia source code (evalfile in base/loading.jl)
	# but with the modification that it returns the module instead of the last object
	name = Symbol(basename(path))
	m = Module(name)
	Core.eval(m,
        Expr(:toplevel,
             :(eval(x) = $(Expr(:core, :eval))($name, x)),
             :(include(x) = $(Expr(:top, :include))($name, x)),
             :(include(mapexpr::Function, x) = $(Expr(:top, :include))(mapexpr, $name, x)),
             :(include($path))))
	m
	end
end;

# ╔═╡ d6e610fe-3aa2-11eb-1397-31d1f09ef073
begin
	using Serialization

	md"""
	## GAME TERMINAL

	Send commands to the game engine via a serialized string. The serialized string file will be deleted after attempting to @eval it. Press the "`" (backquote) key in-game to execute an awaiting serialized string.
	"""
end

# ╔═╡ 7650e106-3aca-11eb-11aa-67e7cd77b705
begin
	terminal = """length(gs[:group][:selected]); println("this is cool!")"""
	serialize(projectdir() * "/terminal.jls", terminal)
	"Command saved to disk @ $(projectdir() * "/terminal.jls")"
end

# ╔═╡ b2696b50-3b5a-11eb-0cd2-ed9365906359
md"""
#### Common commands
"""

# ╔═╡ c2d1f23c-3b5a-11eb-3a24-c3cd27d005d0

# ╔═╡ Cell order:
# ╟─804078c8-3aa1-11eb-2f4d-69380de7ab92
# ╟─d6e610fe-3aa2-11eb-1397-31d1f09ef073
# ╠═7650e106-3aca-11eb-11aa-67e7cd77b705
# ╟─b2696b50-3b5a-11eb-0cd2-ed9365906359
# ╠═c2d1f23c-3b5a-11eb-3a24-c3cd27d005d0


text:
"""

```
merge(r::Vector{ComplexMixutures.Result})
```

This function merges the results of MDDF calculations obtained by running the same
analysis on multiple trajectories, or multiple parts of the same trajectory. It returns
a Result structure of the same type, with all the functions and counters representing averages
of the set provided weighted by the number of frames read in each Result set.

"""
function Base.merge(r::Vector{Result})

  nr = length(r)
  nframes_read = r[1].nframes_read
  error = false
  for ir in 2:nr
    nframes_read += r[ir].nframes_read
    if r[ir].nbins != r[1].nbins
      println("ERROR: To merge Results, the number of bins of the histograms of both sets must be the same.")
    end
    if (r[ir].cutoff - r[1].cutoff) > 1.e-8
      println("ERROR: To merge Results, cutoff distance of the of the histograms of both sets must be the same.")
    end
  end
  if error
    error(" Incompatible set of results to merge. ")
  end

  # List of files and weights
  nfiles = 0
  for ir in 1:nr
    nfiles += length(r[ir].files)
  end
  files = Vector{String}(undef,nfiles)
  weights = Vector{Float64}(undef,nfiles)

  # Final resuls
  R = Result(options=r[1].options,
             nbins=r[1].nbins,
             dbulk=r[1].dbulk,
             cutoff=r[1].cutoff,
             irefatom=r[1].irefatom,
             lastframe_read=r[nr].lastframe_read,
             nframes_read=nframes_read,
             autocorrelation=r[1].autocorrelation,
             solute = r[1].solute,
             solvent = r[1].solvent,
             files=files,
             weights=weights)

  # Average results weighting the data considering the number of frames of each data set
  @. R.d = r[1].d
  ifile = 0
  for ir in 1:nr

    w = r[ir].nframes_read / nframes_read

    @. R.mddf += w*r[ir].mddf
    @. R.kb += w*r[ir].kb

    @. R.rdf += w*r[ir].rdf
    @. R.kb_rdf += w*r[ir].kb_rdf

    @. R.md_count += w*r[ir].md_count
    @. R.md_count_random += w*r[ir].md_count_random

    @. R.sum_md_count += w*r[ir].sum_md_count
    @. R.sum_md_count_random += w*r[ir].sum_md_count_random

    @. R.solute_atom += w*r[ir].solute_atom
    @. R.solvent_atom += w*r[ir].solvent_atom

    @. R.rdf_count += w*r[ir].rdf_count
    @. R.rdf_count_random += w*r[ir].rdf_count_random

    @. R.sum_rdf_count += w*r[ir].sum_rdf_count
    @. R.sum_rdf_count_random += w*r[ir].sum_rdf_count_random

    R.density.solute += w*r[ir].density.solute
    R.density.solvent += w*r[ir].density.solvent
    R.density.solvent_bulk += w*r[ir].density.solvent_bulk

    R.volume.total += w*r[ir].volume.total
    R.volume.bulk += w*r[ir].volume.bulk
    R.volume.domain += w*r[ir].volume.domain
    R.volume.shell += w*r[ir].volume.shell

    for j in 1:length(r[ir].files)
      ifile += 1
      R.files[ifile] = r[ir].files[j]
      R.weights[ifile] = w*r[ir].weights[j]
    end

  end

  return R

end


text:
__precompile__(false)
module OpenNI2

import Base: open, close, start, read, read!
using BinDeps

# Load required dependency
deps = joinpath(Pkg.dir("OpenNI2"), "deps", "deps.jl")
if isfile(deps)
    include(deps)
else
    error("OpenNI2 not properly installed. Please run Pkg.build(\"OpenNI2\")")
end

using Cxx

Libdl.dlopen(libOpenNI2, Libdl.RTLD_GLOBAL)

ni2_header_path = replace(dirname(libOpenNI2), "/lib\/ni2", "/include\/ni2")
addHeaderDir(ni2_header_path, kind=C_System)

cxx"""
#include <memory>
#include <OpenNI.h>
"""

typealias StdSharedPtr{T} cxxt"std::shared_ptr<$T>"

typealias DepthPixel cxxt"openni::DepthPixel" # uint16_t
typealias Grayscale16Pixel cxxt"openni::Grayscale16Pixel" # uint16_t
typealias RGB888Pixel cxxt"openni::RGB888Pixel" # struct

### Enums ###

for name in [
    :STATUS_OK,
    :STATUS_ERROR,
    :STATUS_NOT_IMPLEMENTED,
    :STATUS_NOT_SUPPORTED,
    :STATUS_BAD_PARAMETER,
    :STATUS_OUT_OF_FLOW,
    :STATUS_NO_DEVICE,
    :STATUS_TIME_OUT,
    :SENSOR_IR,
    :SENSOR_COLOR,
    :SENSOR_DEPTH,
    :PIXEL_FORMAT_DEPTH_1_MM,
    :PIXEL_FORMAT_DEPTH_100_UM,
    :PIXEL_FORMAT_SHIFT_9_2,
    :PIXEL_FORMAT_SHIFT_9_3,
    :PIXEL_FORMAT_RGB888,
    :PIXEL_FORMAT_YUV422,
    :PIXEL_FORMAT_GRAY8,
    :PIXEL_FORMAT_GRAY16,
    :PIXEL_FORMAT_JPEG,
    :PIXEL_FORMAT_YUYV,
    :DEVICE_STATE_OK,
    :DEVICE_STATE_ERROR,
    :DEVICE_STATE_NOT_READY,
    :DEVICE_STATE_EOF,
    :IMAGE_REGISTRATION_OFF,
    :IMAGE_REGISTRATION_DEPTH_TO_COLOR,
    ]
    cppname = string("openni::", name)
    ex = Expr(:macrocall, Symbol("@icxx_str"), string(cppname, ";"))
    @eval begin
        global const $name = $ex
        @assert isa($name, Cxx.CppEnum)
    end
end

const ANY_DEVICE = icxx"openni::ANY_DEVICE;"
const TIMEOUT_FOREVER = icxx"openni::TIMEOUT_FOREVER;"
const SensorType = Cxx.CppEnum{Symbol("openni::SensorType"),UInt32}
const ImageRegistrationMode =
    Cxx.CppEnum{Symbol("openni::ImageRegistrationMode"),UInt32}

@inline function checkStatus(rc)
    if rc != STATUS_OK
        error("[OpenNI error (code: $(rc.val))]:\n$(getExtendedError())")
    end
    rc
end

### OpenNI ###

function initialize()
    rc = icxx"openni::OpenNI::initialize();"
    checkStatus(rc)
end
shutdown() = icxx"openni::OpenNI::shutdown();"
getVersion() = icxx"openni::OpenNI::getVersion();"
function getExtendedError()
    icxx"openni::OpenNI::getExtendedError();" |> bytestring
end

# Julia-friendly getVersion
function getVersionNumber()
    version = getVersion()
    major = icxx"$version.major;"
    minor = icxx"$version.minor;"
    maintenance = icxx"$version.maintenance;"
    build = icxx"$version.build;"
    convert(VersionNumber, "$major.$minor.$maintenance-$build")
end

### DeviceInfo ###

type DeviceInfo{T}
    handle::T
end

for f in [
    :getUri,
    :getVendor,
    :getName,
    ]
    ex = Expr(:macrocall, Symbol("@icxx_str"), "\$(di.handle).$f();")
    @eval $f(di::DeviceInfo) = $ex |> bytestring
end

for f in [
    :getUsbVendorId,
    :getUsbProductId,
    ]
    ex = Expr(:macrocall, Symbol("@icxx_str"), "\$(di.handle).$f();")

    @eval $f(di::DeviceInfo) = $ex
end

### SensorInfo ###

# stream->getSensorInfo() return a const reference
type SensorInfo{T}
    handle::T
end

# but device->getSensorInfo(type) return a pointer
type SensorInfoPtr{T}
    handle::T
end

getSensorType(si::SensorInfo) = icxx"$(si.handle).getSensorType();"
getSensorType(si::SensorInfoPtr) = icxx"$(si.handle)->getSensorType();"

### Device ###

type DevicePtr{T<:StdSharedPtr}
    handle::T
end

function (::Type{DevicePtr})()
    handle =  icxx"std::shared_ptr<openni::Device>(new openni::Device);"
    DevicePtr(handle)
end

function open(device::DevicePtr, deviceURI=ANY_DEVICE)
    rc = icxx"$(device.handle)->open($deviceURI);"
    checkStatus(rc)
end

for f in [
    :setDepthColorSyncEnabled,
    ]
    ex = Expr(:macrocall, Symbol("@icxx_str"), "\$(device.handle)->$f(\$v);")
    @eval begin
        function $f(device::DevicePtr, v)
            rc = $ex
            checkStatus(rc)
        end
    end
end

for f in [
    :close,
    :getImageRegistrationMode,
    :isValid,
    :isFile,
    :getDepthColorSyncEnabled,
    ]
    ex = Expr(:macrocall, Symbol("@icxx_str"), "\$(device.handle)->$f();")
    @eval $f(device::DevicePtr) = $ex
end

function hasSensor(device::DevicePtr, sensorType::SensorType)
    icxx"$(device.handle)->hasSensor($sensorType);"
end

function getDeviceInfo(device::DevicePtr)
    handle = icxx"$(device.handle)->getDeviceInfo();"
    @assert isa(handle, Cxx.CppRef)
    DeviceInfo(handle)
end

function getSensorInfo(device::DevicePtr, sensorType::SensorType)
    handle = icxx"$(device.handle)->getSensorInfo($sensorType);"
    @assert isa(handle, Cxx.CppPtr)
    SensorInfoPtr(handle)
end

function isImageRegistrationModeSupported(device::DevicePtr,
        mode::ImageRegistrationMode)
    icxx"$(device.handle)->isImageRegistrationModeSupported($mode);"
end

function setImageRegistrationMode(device::DevicePtr,
        mode::ImageRegistrationMode)
    rc = icxx"$(device.handle)->setImageRegistrationMode($mode);"
    checkStatus(rc)
end

# workaround to dispatch on openni::Array<T>
type SupportedVideoModes{T}
    handle::T
end

function getSupportedVideoModes(si::SensorInfo)
    handle = icxx"$(si.handle).getSupportedVideoModes();"
    SupportedVideoModes(handle)
end

function getSupportedVideoModes(si::SensorInfoPtr)
    handle = icxx"$(si.handle)->getSupportedVideoModes();"
    SupportedVideoModes(handle)
end

# TODO: should dispatch on openni::Array<T>
Base.start(ar::SupportedVideoModes) = 0
Base.next(ar::SupportedVideoModes,i) = (ar[i], i+1)
Base.done(ar::SupportedVideoModes,i) = i >= length(ar)
Base.getindex(ar::SupportedVideoModes,i) = icxx"($(ar.handle))[$i];"
Base.length(ar::SupportedVideoModes) = icxx"$(ar.handle).getSize();"

### VideoStream ###

type VideoStreamPtr{T<:StdSharedPtr}
    handle::T
end

function (::Type{VideoStreamPtr})()
    handle = icxx"std::shared_ptr<openni::VideoStream>(new openni::VideoStream);"
    VideoStreamPtr(handle)
end

function create(stream::VideoStreamPtr, device::DevicePtr,
        typ::SensorType=SENSOR_DEPTH)
    rc = icxx"$(stream.handle)->create(*$(device.handle), $typ);"
    checkStatus(rc)
end

for f in [
    :setVideoMode,
    :setMirroringEnabled,
    ]
    ex = Expr(:macrocall, Symbol("@icxx_str"), "\$(stream.handle)->$f(\$v);")
    @eval begin
        function $f(stream::VideoStreamPtr, v)
            rc = $ex
            checkStatus(rc)
        end
    end
end

for f in [:start, :resetCropping]
    ex = Expr(:macrocall, Symbol("@icxx_str"), "\$(stream.handle)->$f();")
    @eval begin
        function $f(stream::VideoStreamPtr)
            rc = $ex
            checkStatus(rc)
        end
    end
end

for f in [
    :stop,
    :destroy,
    :isValid,
    :getVideoMode,
    :getMaxPixelValue,
    :getMinPixelValue,
    :isCroppingSupported,
    :getMirroringEnabled,
    :getHorizontalFieldOfView,
    :getVerticalFieldOfView,
    ]
    ex = Expr(:macrocall, Symbol("@icxx_str"), "\$(stream.handle)->$f();")
    @eval $f(stream::VideoStreamPtr) = $ex
end

function getSensorInfo(stream::VideoStreamPtr)
    handle = icxx"$(stream.handle)->getSensorInfo();"
    SensorInfo(handle)
end

typealias VideoMode cxxt"openni::VideoMode"
typealias VideoModeRef rcpp"openni::VideoMode"
const VideoModeValOrRef = Union{VideoMode, VideoModeRef}

for f in [
    :getPixelFormat,
    :getResolutionX,
    :getResolutionY,
    :getFps,
    ]
    ex = Expr(:macrocall, Symbol("@icxx_str"), "\$(videoMode).$f();")
    @eval $f(videoMode::VideoModeValOrRef) = $ex
end

function waitForAnyStream(streams; timeout::Cint=TIMEOUT_FOREVER)
    raw_streams = pcpp"openni::VideoStream"[]
    for s in streams
        push!(raw_streams, icxx"$(s.handle).get();")
    end

    readyStreamIndex = Cint[1]
    rc = icxx"""
    openni::OpenNI::waitForAnyStream($(pointer(raw_streams)),
        $(length(raw_streams)),
        $(pointer(readyStreamIndex)),
        $timeout);
    """
    checkStatus(rc)
    readyStreamIndex[1]
end

### VideoFrameRef ###

type VideoFrameRef{T}
    handle::T
end

function (::Type{VideoFrameRef})()
    handle = icxx"openni::VideoFrameRef();"
    VideoFrameRef(handle)
end

function readFrame(stream::VideoStreamPtr, frameRef::VideoFrameRef)
    rc = icxx"$(stream.handle)->readFrame(&$(frameRef.handle));"
    checkStatus(rc)
end
function read!(stream::VideoStreamPtr, frameRef::VideoFrameRef)
    readFrame(stream, frameRef)
end
function read(stream::VideoStreamPtr)
    frame = VideoFrameRef()
    read!(stream, frame)
end

for f in [
    :getDataSize,
    :getData,
    :getSensorType,
    :getVideoMode,
    :getTimestamp,
    :getFrameIndex,
    :getWidth,
    :getHeight,
    :getCroppingEnabled,
    :getCropOriginX,
    :getCropOriginY,
    :getStrideInBytes,
    :isValid,
    :release,
    ]
    ex = Expr(:macrocall, Symbol("@icxx_str"), "\$(frame.handle).$f();")
    @eval $f(frame::VideoFrameRef) = $ex
end

function Base.convert{T,N}(::Type{Array{T,N}}, frame::VideoFrameRef)
    mode = getVideoMode(frame)
    w = Int(getResolutionX(mode))
    h = Int(getResolutionY(mode))
    ar = pointer_to_array(convert(Ptr{T}, getData(frame)), w * h)
    if N == 1
        return ar
    elseif N == 2
        return reshape(ar, w, h)
    else
        error("N <= 2 is supported")
    end
end

### Nite ###

# TODO: to be configurable
const libNiTE2 = Pkg.dir("OpenNI2", "deps", "NiTE", "Redist", "libNiTE2.dylib")
if isfile(libNiTE2)
    export NiTE
    info("Enable NiTE support")
    include("NiTE.jl")
end

end # module


text:
#   This file is part of UnitSystems.jl. It is licensed under the MIT license
#   UnitSystems Copyright (C) 2021 Michael Reed

export Universe, coupling, finestructure, electronunit, protonunit, protonelectron
#const Mu,Ru,SB,hh,cc,m0,e0,ke,me,mp,mu,ee,FF,Z0,G0,Eh,a0,re,g0,lP,ϵ₀,mB = Mᵤ,Rᵤ,σ,𝘩,𝘤,μ₀,ε₀,kₑ,mₑ,mₚ,mᵤ,𝘦,𝔉,Z₀,G₀,Eₕ,a₀,rₑ,g₀,ℓP,ε₀,μB
export slug, ft, KJ1990, KJ2014, RK1990, RK2014, mₑ1990, mₑ2014, temp, units
export slugs, kilograms, lbm, meters, feet, rankine, kelvin, moles, molecules
export UnitSystem, US, SI, MKS, CGS, CGS2019, CGSm, CGSe, HLU, FFF

# == Metric is different
const κ = einstein(SI2019)
const σ = stefan(SI2019) #
const μB = magneton(SI2019) #
const ε₀ = permittivity(SI2019) #
const kₑ = coulomb(SI2019) #
const mₚ = protonmass(SI2019)
const mᵤ = atomicmass(SI2019)
const Mᵤ = molarmass(SI2019)
const 𝔉 = faraday(SI2019) #
const Φ₀ = magneticflux(SI2019) #
const Z₀ = impedance(SI2019) #
const G₀ = conductance(SI2019) #
const Eₕ = hartree(SI2019)
const a₀ = bohr(SI2019)
const rₑ = electronradius(SI2019)
const RK = klitzing(SI2019) #
const KJ = josephson(SI2019) #
const RH,Ry = R∞*mₚ/(mₑ+mₚ),𝘩*𝘤*R∞

const ℓP = length(PlanckGauss,SI2019)
const tP = time(PlanckGauss,SI2019)
const TP = temperature(PlanckGauss,SI2019)

const lS = length(Stoney,SI2019)
const tS = time(Stoney,SI2019)
const mS = mass(Stoney,SI2019)
const qS = charge(Stoney,SI2019)

const lA = length(Hartree,SI2019)
const tA = time(Hartree,SI2019)
const mA = mass(Hartree,SI2019)
const qA = charge(Hartree,SI2019)

const lQCD = length(QCD,SI2019)
const tQCD = time(QCD,SI2019)
const mQCD = mass(QCD,SI2019)

# non standard units

const BTUftlb = 3600/0.5778thermalconductivity(English) # BTU⋅ft⁻¹⋅lb⁻¹
const BTUJ = energy(English)*BTUftlb # BTU⋅J⁻¹

# constant aliases

const mpe, meu, mpu, ainv, aG = μₚₑ, μₑᵤ, μₚᵤ, αinv, αG
const Mu,Ru,SB,hh,cc,m0,e0,ke,me,mp,mu,ee,FF,Z0,G0,Eh,a0,re,g0,lP,aL,ϵ₀ = Mᵤ,Rᵤ,σ,𝘩,𝘤,μ₀,ε₀,kₑ,mₑ,mₚ,mᵤ,𝘦,𝔉,Z₀,G₀,Eₕ,a₀,rₑ,g₀,ℓP,αL,ε₀
export κ, GG, NA, kB, Rᵤ, σ, 𝘩, ħ, 𝘤, μ₀, ε₀, kₑ, mₑ, mₚ, mᵤ, 𝘦, 𝔉, Φ₀, Z₀, G₀, Eₕ, R∞, a₀, rₑ, KJ, RK, Ru, SB, hh, cc, m0, e0, ke, me, mp, mu, ee, FF, Z0, G0, Eh, a0, re, μB
export αG, αinv, μₚₑ, μₑᵤ, μₚᵤ, mpe, meu, mpu, mP, δμ₀, Mᵤ, Mu, RH, Ry, ΔνCs, Kcd, ainv
export cal, kcal, calₜₕ, kcalₜₕ, calᵢₜ, kcalᵢₜ, ℓP, g₀, g0, atm, lbm, BTUJ, BTUftlb, aG
export lP, tP, TP, lS, tS, mS, qS, lA, tA, mA, qA, lQCD, tQCD, mQCD, ϵ₀, αL, aL

# engineering unit systems docs

cgstext(US,AMP,cgs=eval(US)) = """
```Julia
julia> boltzmann($US) # erg⋅K⁻¹
$(boltzmann(cgs))

julia> planckreduced($US) # erg⋅s⋅rad⁻¹
$(planckreduced(cgs))

julia> lightspeed($US) # cm⋅s⁻¹
$(lightspeed(cgs))

julia> permeability($US) # statH⋅cm⁻¹
$(permeability(cgs))

julia> electronmass($US) # g
$(electronmass(cgs))

julia> rationalization($US)
$(rationalization(cgs))
```
"""

for U ∈ (:CGSm,:CGSe,:EMU,:ESU)
    (EU,AMP) = QuoteNode.(U ∉ (:CGSe,:ESU) ? (:EMU,:Bi) : (:ESU,:statA))
@eval @doc """
    $($(QuoteNode(U)))::UnitSystem{1e10*Rᵤ*mₑ/μₑᵤ,1e7*ħ,100𝘤,$($EU≠:EMU ? "(100𝘤)^-2" : 1),1000mₑ,4π}

Centimetre-gram-second `UnitSystem` variant based on `$($EU)` (non-rationalized).

$(cgstext($(QuoteNode(U)),$AMP))
""" $U

U ∉ (:CGSm,:CGSe) && @eval @doc """
    $(Symbol($(QuoteNode(U)),:2019))::UnitSystem{1e7*kB,1e7*ħ,100𝘤,$($EU≠:EMU ? "1e3*μ₀/𝘤^2" : "1e7*μ₀"),1000mₑ}

Centimetre-gram-second `UnitSystem` variant of tuned `SI2019` based on `$($EU)` (rationalized).

$(cgstext(Symbol($(QuoteNode(U)),:2019),$AMP))
""" $(Symbol(U,:2019))
end

@doc """
    Thomson::UnitSystem{1e10*Rᵤ*mₑ/μₑᵤ,1e7*ħ,100𝘤,1,1000mₑ,4π,1/2}

Centimetre-gram-second `UnitSystem` variant `Thomson` (EMU-Lorentz, non-rationalized).

```Julia
julia> boltzmann(Thomson) # erg⋅K⁻¹
$(boltzmann(Thomson))

julia> planckreduced(Thomson) # erg⋅s⋅rad⁻¹
$(planckreduced(Thomson))

julia> lightspeed(Thomson) # cm⋅s⁻¹
$(lightspeed(Thomson))

julia> permeability(Thomson) # abH⋅cm⁻¹
$(permeability(Thomson))

julia> electronmass(Thomson) # g
$(electronmass(Thomson))

julia> rationalization(Thomson)
$(rationalization(Thomson))

julia> lorentz(Thomson)
$(lorentz(Thomson))
```
""" Thomson

@doc """
    Gauss::UnitSystem{1e10*Rᵤ*mₑ/μₑᵤ,1e7*ħ,100𝘤,1,1000mₑ,4π,0.01/𝘤}

Centimetre-gram-second `UnitSystem` variant `CGS` (Gauss-Lorentz, non-rationalized).

```Julia
julia> boltzmann(Gauss) # erg⋅K⁻¹
$(boltzmann(Gauss))

julia> planckreduced(Gauss) # erg⋅s⋅rad⁻¹
$(planckreduced(Gauss))

julia> lightspeed(Gauss) # cm⋅s⁻¹
$(lightspeed(Gauss))

julia> permeability(Gauss) # statH⋅cm⁻¹
$(permeability(Gauss))

julia> electronmass(Gauss) # g
$(electronmass(Gauss))

julia> rationalization(Gauss)
$(rationalization(Gauss))

julia> lorentz(Gauss)
$(lorentz(Gauss))
```
""" Gauss, CGS

@doc """
    LorentzHeaviside::UnitSystem{1e10*Rᵤ*mₑ/μₑᵤ,1e7*ħ,100𝘤,1,1000mₑ,1,0.01/𝘤}

Centimetre-gram-second `UnitSystem` variant `HLU` (Heaviside-Lorentz, rationalized).

```Julia
julia> boltzmann(LorentzHeaviside) # erg⋅K⁻¹
$(boltzmann(LorentzHeaviside))

julia> planckreduced(LorentzHeaviside) # erg⋅s⋅rad⁻¹
$(planckreduced(LorentzHeaviside))

julia> lightspeed(LorentzHeaviside) # cm⋅s⁻¹
$(lightspeed(LorentzHeaviside))

julia> permeability(HLU) # hlH⋅cm⁻¹
$(permeability(LorentzHeaviside))

julia> electronmass(LorentzHeaviside) # g
$(electronmass(LorentzHeaviside))

julia> rationalization(LorentzHeaviside)
$(rationalization(LorentzHeaviside))

julia> lorentz(LorentzHeaviside)
$(lorentz(LorentzHeaviside))
```
""" LorentzHeaviside, HLU

@doc """
    MTS::UnitSystem{1e6*Rᵤ*mₑ/μₑᵤ,1000ħ,𝘤,4π/1e4,mₑ/1000}

Metre-tonne-second `UnitSystem` variant of `Metric` system.

```Julia
julia> boltzmann(MTS) # kJ⋅K⁻¹
$(boltzmann(MTS))

julia> planckreduced(MTS) # kJ⋅s⋅rad⁻¹
$(planckreduced(MTS))

julia> lightspeed(MTS) # m⋅s⁻¹
$(lightspeed(MTS))

julia> permeability(MTS) # kH⋅m⁻¹
$(permeability(MTS))

julia> electronmass(MTS) # t
$(electronmass(MTS))
```
""" MTS

@doc """
    Metric::UnitSystem{Rᵤ*mₑ/μₑᵤ/0.001,ħ,𝘤,4π*1e-7,mₑ}

Systeme International d'Unites (the SI units) adopted as the preferred `UnitSystem`.

```Julia
julia> boltzmann(Metric) # J⋅K⁻¹
$(boltzmann(Metric))

julia> planckreduced(Metric) # J⋅s⋅rad⁻¹
$(planckreduced(Metric))

julia> lightspeed(Metric) # m⋅s⁻¹
$(lightspeed(Metric))

julia> permeability(Metric) # H⋅m⁻¹
$(permeability(Metric))

julia> electronmass(Metric) # kg
$(electronmass(Metric))
```
""" Metric

@doc """
    SI2019::UnitSystem{kB,ħ,𝘤,μ₀,mₑ}

Systeme International d'Unites (the SI units) with `μ₀` for a tuned `charge` exactly.

```Julia
julia> boltzmann(SI2019) # J⋅K⁻¹
$(boltzmann(SI2019))

julia> planckreduced(SI2019) # J⋅s⋅rad⁻¹
$(planckreduced(SI2019))

julia> lightspeed(SI2019) # m⋅s⁻¹
$(lightspeed(SI2019))

julia> permeability(SI2019) # H⋅m⁻¹
$(permeability(SI2019))

julia> electronmass(SI2019) # kg
$(electronmass(SI2019))
```
""" SI2019, SI

@doc """
    CODATA::UnitSystem{Rᵤ*mₑ2014/μₑᵤ/0.001,2/RK2014/KJ2014^2/π,𝘤,2RK2014/𝘤/αinv,mₑ2014}

Metric `UnitSystem` based on Committee on Data of the International Science Council.

```Julia
julia> boltzmann(CODATA) # J⋅K⁻¹
$(boltzmann(CODATA))

julia> planckreduced(CODATA) # J⋅s⋅rad⁻¹
$(planckreduced(CODATA))

julia> lightspeed(CODATA) # m⋅s⁻¹
$(lightspeed(CODATA))

julia> permeability(CODATA) # H⋅m⁻¹
$(permeability(CODATA))

julia> electronmass(CODATA) # kg
$(electronmass(CODATA))
```
""" CODATA

@doc """
    Conventional::UnitSystem{Rᵤ*mₑ1990/μₑᵤ/0.001,2/RK1990/KJ1990^2/π,𝘤,2RK1990/𝘤/αinv,mₑ1990}

Conventional electronic `UnitSystem` with 1990 tuned `josephson` and `klitzing` constants.

```Julia
julia> boltzmann(Conventional) # J⋅K⁻¹
$(boltzmann(Conventional))

julia> planckreduced(Conventional) # J⋅s⋅rad⁻¹
$(planckreduced(Conventional))

julia> lightspeed(Conventional) # m⋅s⁻¹
$(lightspeed(Conventional))

julia> permeability(Conventional) # H⋅m⁻¹
$(permeability(Conventional))

julia> electronmass(Conventional) # kg
$(electronmass(Conventional))
```
""" Conventional

@doc """
    IAU::UnitSystem{Rᵤ*mₑ/μₑᵤ/0.001/Jₛ,ħ*day/Jₛ,day*𝘤/au,4π*1e-7*day^2/Jₛ,mₑ/mₛ}

Astronomical (solar) `UnitSystem` defined by International Astronomical Union.

```Julia
julia> boltzmann(IAU) # M⊙⋅au²⋅D⁻²⋅K⁻¹
$(boltzmann(IAU))

julia> planckreduced(IAU) # M⊙⋅au²⋅D⁻¹⋅rad⁻¹
$(planckreduced(IAU))

julia> lightspeed(IAU) # au⋅D⁻¹
$(lightspeed(IAU))

julia> permeability(IAU) # M⊙⋅au²⋅C⁻²
$(permeability(IAU))

julia> electronmass(IAU) # M⊙
$(electronmass(IAU))
```
""" IAU

@doc """
    English::UnitSystem{kB*rankine/slug/ft^2,ħ/slug/ft^2,𝘤/ft,4π,mₑ/slug}

Engineering `UnitSystem` historically used by Britain and United States.

```Julia
julia> boltzmann(English) # ft⋅lb⋅°R⁻¹
$(boltzmann(English))

julia> planckreduced(English) # ft⋅lb⋅s⋅rad⁻¹
$(planckreduced(English))

julia> lightspeed(English) # ft⋅s⁻¹
$(lightspeed(English))

julia> permeability(English) # slug⋅ft²⋅?⁻²
$(permeability(English))

julia> electronmass(English) # slugs
$(electronmass(English))
```
""" English

@doc """
    EnglishUS::UnitSystem{1000Rᵤ*mₑ/μₑᵤ*rankine/slug/ftUS^2,ħ/slug/ftUS^2,𝘤/ftUS,4π,mₑ/slug}

Engineering `UnitSystem` based on the geophysical US survey foot (1200/3937).

```Julia
julia> boltzmann(EnglishUS) # ftUS⋅lb⋅°R⁻¹
$(boltzmann(EnglishUS))

julia> planckreduced(EnglishUS) # ftUS⋅lb⋅s⋅rad⁻¹
$(planckreduced(EnglishUS))

julia> lightspeed(EnglishUS) # ftUS⋅s⁻¹
$(lightspeed(EnglishUS))

julia> permeability(EnglishUS) # slug⋅ftUS²⋅?⁻²
$(permeability(EnglishUS))

julia> electronmass(EnglishUS) # slugs
$(electronmass(EnglishUS))
```
""" EnglishUS

@doc """
    FFF::UnitSystem{Rᵤ*mₑ/μₑᵤ/0.001*rankine/Jf,ħ/14day/Jf,14day*𝘤/201.168,0,mₑ/mf}

Furlong–firkin–fortnight `FFF` is a humorous `UnitSystem` based on unusal impractical units.

```Julia
julia> boltzmann(FFF) # fir⋅fur²⋅ftn⁻²⋅F⁻¹
$(boltzmann(FFF))

julia> planckreduced(FFF) # fir⋅fur²⋅ftn⁻¹⋅rad⁻¹
$(planckreduced(FFF))

julia> lightspeed(FFF) # fur⋅ftn⁻¹
$(lightspeed(FFF))

julia> permeability(FFF) # fir⋅fur²⋅Inf⁻²
$(permeability(FFF))

julia> electronmass(FFF) # fir
$(electronmass(FFF))
```
""" FFF

@doc """
    Kennelly::UnitSystem{Rᵤ*mₑ/μₑᵤ/0.001,ħ,𝘤,1e-7,mₑ,4π}

Kennelly ? variant `UnitSystem` of the standard `Metric` units ???

```Julia
julia> boltzmann(Kennelly) # J⋅K⁻¹
$(boltzmann(Kennelly))

julia> planckreduced(Kennelly) # J⋅s⋅rad⁻¹
$(planckreduced(Kennelly))

julia> lightspeed(Kennelly) # m⋅s⁻¹
$(lightspeed(Kennelly))

julia> permeability(Kennelly) # H⋅m⁻¹
$(permeability(Kennelly))

julia> electronmass(Kennelly) # kg
$(electronmass(Kennelly))

julia> rationalization(Kennelly)
$(rationalization(Kennelly))
```
""" Kennelly

# natural unit system docs

textunits(U,S) = """
```Julia
julia> boltzmann($S)
$(boltzmann(U))

julia> planckreduced($S)
$(planckreduced(U))

julia> lightspeed($S)
$(lightspeed(U))

julia> permeability($S)
$(permeability(U))

julia> electronmass($S)
$(electronmass(U))
```
"""

@doc """
    Planck::UnitSystem{1,1,1,1,√(4π*αG)}

Planck `UnitSystem` with the `electronmass` value `√(4π*αG)` using gravitational coupling.

$(textunits(Planck,:Planck))
""" Planck

@doc """
    PlanckGauss::UnitSystem{1,1,1,4π,√αG}

Planck (Gauss) `UnitSystem` with `permeability` of `4π` and `electronmass` coupling `√αG`.

$(textunits(PlanckGauss,:PlanckGauss))

The well known `PlanckGauss` values for `length`, `time`, `mass`, and `temperature` are:
```Julia
julia> length(PlanckGauss,Metric) # ℓP
$(length(PlanckGauss,Metric))

julia> time(PlanckGauss,Metric) # tP
$(time(PlanckGauss,Metric))

julia> mass(PlanckGauss,Metric) # mP
$(mass(PlanckGauss,Metric))

julia> temperature(PlanckGauss,Metric) # TP
$(temperature(PlanckGauss,Metric))
```
""" PlanckGauss

@doc """
    Stoney::UnitSystem{1,αinv,1,4π,√(αG*αinv)}

Stoney `UnitSystem` with `permeability` of `4π` and `electronmass` coupling `√(αG*αinv)`.

$(textunits(Stoney,:Stoney))

The well known `Stoney` values for `length`, `time`, `mass`, and `charge` are:
```Julia
julia> length(Stoney,Metric) # lS
$(length(Stoney,Metric))

julia> time(Stoney,Metric) # tS
$(time(Stoney,Metric))

julia> mass(Stoney,Metric) # mS
$(mass(Stoney,Metric))

julia> charge(Stoney,Metric) # qS
$(charge(Stoney,Metric))
```
""" Stoney

@doc """
    Hartree::UnitSystem{1,1,αinv,4π/αinv^2,1}

Hartree atomic `UnitSystem` with `lightspeed` of `αinv` and `permeability` of `4π/αinv^2`.

$(textunits(Hartree,:Hartree))

The well known `Hartree` atomic unit values for `length`, `time`, `mass`, and `charge` are:
```Julia
julia> length(Hartree,Metric) # lA
$(length(Hartree,Metric))

julia> time(Hartree,Metric) # tA
$(time(Hartree,Metric))

julia> mass(Hartree,Metric) # mA
$(mass(Hartree,Metric))

julia> charge(Hartree,Metric) # qA
$(charge(Hartree,Metric))
```
""" Hartree

@doc """
    Rydberg::UnitSystem{1,1,2αinv,π/αinv^2,1/2}

Rydberg `UnitSystem` with `lightspeed` of `2αinv` and `permeability` of `π/αinv^2`.

$(textunits(Rydberg,:Rydberg))
""" Rydberg

@doc """
    Schrodinger::UnitSystem{1,1,αinv,4π/αinv^2,√(αG*αinv)}

Schrodinger `UnitSystem` with `permeability` of `4π/αinv^2` and `electronmass` of `√(αG*αinv)`.

$(textunits(Schrodinger,:Schrodinger))
""" Schrodinger

@doc """
    Electronic::UnitSystem{1,αinv,1,4π,1}

Electronic `UnitSystem` with `planckreduced` of `αinv` and `permeability` of `4π`.

$(textunits(Electronic,:Electronic))
""" Electronic

@doc """
    Natural::UnitSystem{1,1,1,1,1}

Natural `UnitSystem` with all primary constants having unit value.

$(textunits(Natural,:Natural))

The well known `Natural` values for `length`, `time`, `mass`, and `charge` are:
```Julia
julia> length(Natural,Metric)
$(length(Natural,Metric))

julia> time(Natural,Metric)
$(time(Natural,Metric))

julia> mass(Natural,Metric)
$(mass(Natural,Metric))

julia> charge(Natural,Metric)
$(charge(Natural,Metric))
```
""" Natural

@doc """
    NaturalGauss::UnitSystem{1,1,1,4π,1}

Natural (Gauss) `UnitSystem` with the Gaussian `permeability` value of `4π`.

$(textunits(NaturalGauss,:NaturalGauss))
""" NaturalGauss

@doc """
    QCD::UnitSystem{1,1,1,1,1/μₚₑ}

Qunatum chromodynamics `UnitSystem` with `electronmass` of `1/μₚₑ` or `1/$μₚₑ`.

$(textunits(QCD,:QCD))

The well known `QCD` values for `length`, `time`, `mass`, and `charge` are:
```Julia
julia> length(QCD,Metric) # lQCD
$(length(QCD,Metric))

julia> time(QCD,Metric) # tQCD
$(time(QCD,Metric))

julia> mass(QCD,Metric) # mQCD
$(mass(QCD,Metric))

julia> charge(QCD,Metric)
$(charge(QCD,Metric))
```
""" QCD

@doc """
    QCDGauss::UnitSystem{1,1,1,4π,1/μₚₑ}

Qunatum chromodynamics (Gauss) `UnitSystem` with `electronmass` of `1/μₚₑ`.

$(textunits(QCDGauss,:QCDGauss))

The well known `QCDGauss` values for `length`, `time`, `mass`, and `charge` are:
```Julia
julia> length(QCDGauss,Metric) # lQCD
$(length(QCDGauss,Metric))

julia> time(QCDGauss,Metric) # tQCD
$(time(QCDGauss,Metric))

julia> mass(QCDGauss,Metric) # mQCD
$(mass(QCDGauss,Metric))

julia> charge(QCDGauss,Metric)
$(charge(QCDGauss,Metric))
```
""" QCDGauss

@doc """
    QCDoriginal::UnitSystem{1,1,1,4π/αinv,1/μₚₑ}

Qunatum chromodynamics (original) `UnitSystem` with `permeability` of `4π/αinv`.

$(textunits(QCDoriginal,:QCDoriginal))

The well known `QCDoriginal` values for `length`, `time`, `mass`, and `charge` are:
```Julia
julia> length(QCDoriginal,Metric) # lQCD
$(length(QCDoriginal,Metric))

julia> time(QCDoriginal,Metric) # tQCD
$(time(QCDoriginal,Metric))

julia> mass(QCDoriginal,Metric) # mQCD
$(mass(QCDoriginal,Metric))

julia> charge(QCDoriginal,Metric)
$(charge(QCDoriginal,Metric))
```
""" QCDoriginal


text:
module BrownianModels

using Phylodendron2

include("simulate.jl")
export
	simulate_bm

include("rel.jl")
export
	v2brlength!,
	calc_llh!,
	init_model!,
	optimv!,
	phylip_llh

end


text:
@testset "Analysis tools" begin
    # minimalistic test
    @testset "Relative Energy regain" begin
        x = [1f0 -2; 4 5; 7 8 ]
        @test DeconvOptim.relative_energy_regain(x, x .* 1) == (Float32[0.0, 0.333, 0.5, 0.601], Float32[1.0, 1.0, 1.0, 1.0])
        @test DeconvOptim.relative_energy_regain(x, x .* 0.5) == (Float32[0.0, 0.333, 0.5, 0.601], Float32[0.75, 0.75, 0.75, 0.75])
    end

    # minimalistic test
    @testset "Normalized Cross Correlation" begin
        x = [1f0 -2; 4 5; 7 8 ]

        @test DeconvOptim.normalized_cross_correlation(x, x) == 1.0f0
    end

    @testset "Trace Deconvolution" begin
        sz = (10,10)
        x = rand(sz...)
        psf = rand(10,10)
        psf /= sum(psf)
        y = DeconvOptim.conv(x,psf)
        # test whether starting with the ground truth really yield the perfect reconstruction after 0 iterations
        opt_options, summary = DeconvOptim.options_trace_deconv(x, 0, nothing)
        res = deconvolution(y,psf; initial=x, mapping=nothing, padding=0.0, opt_options=opt_options)

        @test summary["nccs"][1] > 0.999
        @test summary["nvars"][1] < 0.001
        @test (summary["best_nvar_img"] ≈ x)
    end
end


text:
using ESGIWeatherLogistics2021
using Test

@testset "ESGIWeatherLogistics2021.jl" begin
    # Write your tests here.
end


text:
# This file was generated, do not modify it. # hide
ŷ = MLJ.predict(mtm, rows=test)
round(rms(ŷ, y[test])^2, sigdigits=4)

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials


text:
@testset "inference.jl" begin
    # Only test threading if 1.3+.
    if VERSION > v"1.2"
        @testset "threaded sampling" begin
            # Test that chains with the same seed will sample identically.
            @testset "rng" begin
                model = gdemo_default

                # multithreaded sampling with PG causes segfaults on Julia 1.5.4
                # https://github.com/TuringLang/Turing.jl/issues/1571
                samplers = @static if VERSION <= v"1.5.3" || VERSION >= v"1.6.0"
                    (
                        HMC(0.1, 7),
                        PG(10),
                        IS(),
                        MH(),
                        Gibbs(PG(3, :s), HMC(0.4, 8, :m)),
                        Gibbs(HMC(0.1, 5, :s), ESS(:m)),
                    )
                else
                    (
                        HMC(0.1, 7),
                        IS(),
                        MH(),
                        Gibbs(HMC(0.1, 5, :s), ESS(:m)),
                    )
                end
                for sampler in samplers
                    Random.seed!(5)
                    chain1 = sample(model, sampler, MCMCThreads(), 1000, 4)

                    Random.seed!(5)
                    chain2 = sample(model, sampler, MCMCThreads(), 1000, 4)

                    @test chain1.value == chain2.value
                end
            end

            # Smoke test for default sample call.
            Random.seed!(100)
            chain = sample(gdemo_default, HMC(0.1, 7), MCMCThreads(), 1000, 4)
            check_gdemo(chain)

            # run sampler: progress logging should be disabled and
            # it should return a Chains object
            sampler = Sampler(HMC(0.1, 7), gdemo_default)
            chains = sample(gdemo_default, sampler, MCMCThreads(), 1000, 4)
            @test chains isa MCMCChains.Chains
        end
    end
    @testset "chain save/resume" begin
        Random.seed!(1234)

        alg1 = HMCDA(1000, 0.65, 0.15)
        alg2 = PG(20)
        alg3 = Gibbs(PG(30, :s), HMC(0.2, 4, :m))

        chn1 = sample(gdemo_default, alg1, 5000; save_state=true)
        check_gdemo(chn1)

        chn1_resumed = Turing.Inference.resume(chn1, 1000)
        check_gdemo(chn1_resumed)

        chn1_contd = sample(gdemo_default, alg1, 5000; resume_from=chn1)
        check_gdemo(chn1_contd)

        chn1_contd2 = sample(gdemo_default, alg1, 5000; resume_from=chn1, reuse_spl_n=1000)
        check_gdemo(chn1_contd2)

        chn2 = sample(gdemo_default, alg2, 1000; save_state=true)
        check_gdemo(chn2)

        chn2_contd = sample(gdemo_default, alg2, 1000; resume_from=chn2)
        check_gdemo(chn2_contd)

        chn3 = sample(gdemo_default, alg3, 5000; save_state=true)
        check_gdemo(chn3)

        chn3_contd = sample(gdemo_default, alg3, 1000; resume_from=chn3)
        check_gdemo(chn3_contd)
    end
    @testset "Contexts" begin
        # Test LikelihoodContext
        @model function testmodel1(x)
            a ~ Beta()
            lp1 = getlogp(__varinfo__)
            x[1] ~ Bernoulli(a)
            global loglike = getlogp(__varinfo__) - lp1
        end
        model = testmodel1([1.0])
        varinfo = Turing.VarInfo(model)
        model(varinfo, Turing.SampleFromPrior(), Turing.LikelihoodContext())
        @test getlogp(varinfo) == loglike

        # Test MiniBatchContext
        @model function testmodel2(x)
            a ~ Beta()
            x[1] ~ Bernoulli(a)
        end
        model = testmodel2([1.0])
        varinfo1 = Turing.VarInfo(model)
        varinfo2 = deepcopy(varinfo1)
        model(varinfo1, Turing.SampleFromPrior(), Turing.LikelihoodContext())
        model(varinfo2, Turing.SampleFromPrior(), Turing.MiniBatchContext(Turing.LikelihoodContext(), 10))
        @test isapprox(getlogp(varinfo2) / getlogp(varinfo1), 10)
    end
    @testset "Prior" begin
        N = 5000

        # Note that all chains contain 3 values per sample: 2 variables + log probability
        Random.seed!(100)
        chains = sample(gdemo_d(), Prior(), N)
        @test chains isa MCMCChains.Chains
        @test size(chains) == (N, 3, 1)
        @test mean(chains, :s) ≈ 3 atol=0.1
        @test mean(chains, :m) ≈ 0 atol=0.1

        Random.seed!(100)
        chains = sample(gdemo_d(), Prior(), MCMCThreads(), N, 4)
        @test chains isa MCMCChains.Chains
        @test size(chains) == (N, 3, 4)
        @test mean(chains, :s) ≈ 3 atol=0.1
        @test mean(chains, :m) ≈ 0 atol=0.1

        Random.seed!(100)
        chains = sample(gdemo_d(), Prior(), N; chain_type = Vector{NamedTuple})
        @test chains isa Vector{<:NamedTuple}
        @test length(chains) == N
        @test all(length(x) == 3 for x in chains)
        @test all(haskey(x, :lp) for x in chains)
        @test mean(x[:s][1] for x in chains) ≈ 3 atol=0.1
        @test mean(x[:m][1] for x in chains) ≈ 0 atol=0.1
    end

    @testset "chain ordering" begin
        for alg in (Prior(), Emcee(10, 2.0))
            chain_sorted = sample(gdemo_default, alg, 1, sort_chain=true)
            @test names(MCMCChains.get_sections(chain_sorted, :parameters)) == [:m, :s]

            chain_unsorted = sample(gdemo_default, alg, 1, sort_chain=false)
            @test names(MCMCChains.get_sections(chain_unsorted, :parameters)) == [:s, :m]
        end
    end

    @testset "chain iteration numbers" begin
        for alg in (Prior(), Emcee(10, 2.0))
            chain = sample(gdemo_default, alg, 10)
            @test range(chain) == 1:10

            chain = sample(gdemo_default, alg, 10; discard_initial=5, thinning=2)
            @test range(chain) == range(6; step=2, length=10)
        end
    end

    # Copy-paste from integration tests in DynamicPPL.
    @testset "assume" begin
        @model function test_assume()
            x ~ Bernoulli(1)
            y ~ Bernoulli(x / 2)
            return x, y
        end

        smc = SMC()
        pg = PG(10)

        res1 = sample(test_assume(), smc, 1000)
        res2 = sample(test_assume(), pg, 1000)

        check_numerical(res1, [:y], [0.5]; atol=0.1)
        check_numerical(res2, [:y], [0.5]; atol=0.1)

        # Check that all xs are 1.
        @test all(isone, res1[:x])
        @test all(isone, res2[:x])
    end
    @testset "beta binomial" begin
        prior = Beta(2, 2)
        obs = [0, 1, 0, 1, 1, 1, 1, 1, 1, 1]
        exact = Beta(prior.α + sum(obs), prior.β + length(obs) - sum(obs))
        meanp = exact.α / (exact.α + exact.β)

        @model function testbb(obs)
            p ~ Beta(2, 2)
            x ~ Bernoulli(p)
            for i in 1:length(obs)
                obs[i] ~ Bernoulli(p)
            end
            return p, x
        end

        smc = SMC()
        pg = PG(10)
        gibbs = Gibbs(HMC(0.2, 3, :p), PG(10, :x))

        chn_s = sample(testbb(obs), smc, 1000)
        chn_p = sample(testbb(obs), pg, 2000)
        chn_g = sample(testbb(obs), gibbs, 1500)

        check_numerical(chn_s, [:p], [meanp]; atol=0.05)
        check_numerical(chn_p, [:x], [meanp]; atol=0.1)
        check_numerical(chn_g, [:x], [meanp]; atol=0.1)
    end
    @testset "forbid global" begin
        xs = [1.5 2.0]
        # xx = 1

        @model function fggibbstest(xs)
            s ~ InverseGamma(2, 3)
            m ~ Normal(0, sqrt(s))
            # xx ~ Normal(m, sqrt(s)) # this is illegal

            for i in 1:length(xs)
                xs[i] ~ Normal(m, sqrt(s))
                # for xx in xs
                # xx ~ Normal(m, sqrt(s))
            end
            return s, m
        end

        gibbs = Gibbs(PG(10, :s), HMC(0.4, 8, :m))
        chain = sample(fggibbstest(xs), gibbs, 2)
    end
    @testset "new grammar" begin
        x = Float64[1 2]

        @model function gauss(x)
            priors = TArray{Float64}(2)
            priors[1] ~ InverseGamma(2, 3)         # s
            priors[2] ~ Normal(0, sqrt(priors[1])) # m
            for i in 1:length(x)
                x[i] ~ Normal(priors[2], sqrt(priors[1]))
            end
            return priors
        end

        chain = sample(gauss(x), PG(10), 10)
        chain = sample(gauss(x), SMC(), 10)

        @model function gauss2(::Type{TV}=Vector{Float64}; x) where {TV}
            priors = TV(undef, 2)
            priors[1] ~ InverseGamma(2, 3)         # s
            priors[2] ~ Normal(0, sqrt(priors[1])) # m
            for i in 1:length(x)
                x[i] ~ Normal(priors[2], sqrt(priors[1]))
            end
            return priors
        end

        chain = sample(gauss2(; x=x), PG(10), 10)
        chain = sample(gauss2(; x=x), SMC(), 10)

        chain = sample(gauss2(Vector{Float64}; x=x), PG(10), 10)
        chain = sample(gauss2(Vector{Float64}; x=x), SMC(), 10)
    end
    @testset "new interface" begin
        obs = [0, 1, 0, 1, 1, 1, 1, 1, 1, 1]

        @model function newinterface(obs)
            p ~ Beta(2, 2)
            for i in 1:length(obs)
                obs[i] ~ Bernoulli(p)
            end
            return p
        end

        chain = sample(
            newinterface(obs), HMC{Turing.ForwardDiffAD{2}}(0.75, 3, :p, :x), 100
        )
    end
    @testset "no return" begin
        Random.seed!(5)
        @model function noreturn(x)
            s ~ InverseGamma(2, 3)
            m ~ Normal(0, sqrt(s))
            for i in 1:length(x)
                x[i] ~ Normal(m, sqrt(s))
            end
        end

        chain = sample(noreturn([1.5 2.0]), HMC(0.1, 10), 4000)
        check_numerical(chain, [:s, :m], [49 / 24, 7 / 6])
    end
    @testset "observe" begin
        Random.seed!(5)
        @model function test()
            z ~ Normal(0, 1)
            x ~ Bernoulli(1)
            1 ~ Bernoulli(x / 2)
            0 ~ Bernoulli(x / 2)
            return x
        end

        is = IS()
        smc = SMC()
        pg = PG(10)

        res_is = sample(test(), is, 10000)
        res_smc = sample(test(), smc, 1000)
        res_pg = sample(test(), pg, 100)

        @test all(isone, res_is[:x])
        @test res_is.logevidence ≈ 2 * log(0.5)

        @test all(isone, res_smc[:x])
        @test res_smc.logevidence ≈ 2 * log(0.5)

        @test all(isone, res_pg[:x])
    end
    @testset "sample" begin
        alg = Gibbs(HMC(0.2, 3, :m), PG(10, :s))
        chn = sample(gdemo_default, alg, 1000)
    end
    @testset "vectorization @." begin
        # https://github.com/FluxML/Tracker.jl/issues/119
        if Turing.Core.ADBackend() !== Turing.Core.TrackerAD
            @model function vdemo1(x)
                s ~ InverseGamma(2, 3)
                m ~ Normal(0, sqrt(s))
                @. x ~ Normal(m, sqrt(s))
                return s, m
            end

            alg = HMC(0.01, 5)
            x = randn(100)
            res = sample(vdemo1(x), alg, 250)

            @model function vdemo1b(x)
                s ~ InverseGamma(2, 3)
                m ~ Normal(0, sqrt(s))
                @. x ~ Normal(m, $(sqrt(s)))
                return s, m
            end

            res = sample(vdemo1b(x), alg, 250)

            @model function vdemo2(x)
                μ ~ MvNormal(zeros(size(x, 1)), I)
                @. x ~ $(MvNormal(μ, I))
            end

            D = 2
            alg = HMC(0.01, 5)
            res = sample(vdemo2(randn(D, 100)), alg, 250)

            # Vector assumptions
            N = 10
            setchunksize(N)
            alg = HMC(0.2, 4)

            @model function vdemo3()
                x = Vector{Real}(undef, N)
                for i in 1:N
                    x[i] ~ Normal(0, sqrt(4))
                end
            end

            t_loop = @elapsed res = sample(vdemo3(), alg, 1000)

            # Test for vectorize UnivariateDistribution
            @model function vdemo4()
                x = Vector{Real}(undef, N)
                @. x ~ Normal(0, 2)
            end

            t_vec = @elapsed res = sample(vdemo4(), alg, 1000)

            @model vdemo5() = x ~ MvNormal(zeros(N), 4 * I)

            t_mv = @elapsed res = sample(vdemo5(), alg, 1000)

            println("Time for")
            println("  Loop : ", t_loop)
            println("  Vec  : ", t_vec)
            println("  Mv   : ", t_mv)

            # Transformed test
            @model function vdemo6()
                x = Vector{Real}(undef, N)
                @. x ~ InverseGamma(2, 3)
            end

            sample(vdemo6(), alg, 1000)

            N = 3
            @model function vdemo7()
                x = Array{Real}(undef, N, N)
                @. x ~ [InverseGamma(2, 3) for i in 1:N]
            end

            sample(vdemo7(), alg, 1000)
        end
    end
    @testset "vectorization .~" begin
        @model function vdemo1(x)
            s ~ InverseGamma(2, 3)
            m ~ Normal(0, sqrt(s))
            x .~ Normal(m, sqrt(s))
            return s, m
        end

        alg = HMC(0.01, 5)
        x = randn(100)
        res = sample(vdemo1(x), alg, 250)

        @model function vdemo2(x)
            μ ~ MvNormal(zeros(size(x, 1)), I)
            return x .~ MvNormal(μ, I)
        end

        D = 2
        alg = HMC(0.01, 5)
        res = sample(vdemo2(randn(D, 100)), alg, 250)

        # Vector assumptions
        N = 10
        setchunksize(N)
        alg = HMC(0.2, 4)

        @model function vdemo3()
            x = Vector{Real}(undef, N)
            for i in 1:N
                x[i] ~ Normal(0, sqrt(4))
            end
        end

        t_loop = @elapsed res = sample(vdemo3(), alg, 1000)

        # Test for vectorize UnivariateDistribution
        @model function vdemo4()
            x = Vector{Real}(undef, N)
            return x .~ Normal(0, 2)
        end

        t_vec = @elapsed res = sample(vdemo4(), alg, 1000)

        @model vdemo5() = x ~ MvNormal(zeros(N), 4 * I)

        t_mv = @elapsed res = sample(vdemo5(), alg, 1000)

        println("Time for")
        println("  Loop : ", t_loop)
        println("  Vec  : ", t_vec)
        println("  Mv   : ", t_mv)

        # Transformed test
        @model function vdemo6()
            x = Vector{Real}(undef, N)
            return x .~ InverseGamma(2, 3)
        end

        sample(vdemo6(), alg, 1000)

        @model function vdemo7()
            x = Array{Real}(undef, N, N)
            return x .~ [InverseGamma(2, 3) for i in 1:N]
        end

        sample(vdemo7(), alg, 1000)
    end
    @testset "Type parameters" begin
        N = 10
        setchunksize(N)
        alg = HMC(0.01, 5)
        x = randn(1000)
        @model function vdemo1(::Type{T}=Float64) where {T}
            x = Vector{T}(undef, N)
            for i in 1:N
                x[i] ~ Normal(0, sqrt(4))
            end
        end

        t_loop = @elapsed res = sample(vdemo1(), alg, 250)
        t_loop = @elapsed res = sample(vdemo1(Float64), alg, 250)

        vdemo1kw(; T) = vdemo1(T)
        t_loop = @elapsed res = sample(vdemo1kw(; T=Float64), alg, 250)

        @model function vdemo2(::Type{T}=Float64) where {T<:Real}
            x = Vector{T}(undef, N)
            @. x ~ Normal(0, 2)
        end

        t_vec = @elapsed res = sample(vdemo2(), alg, 250)
        t_vec = @elapsed res = sample(vdemo2(Float64), alg, 250)

        vdemo2kw(; T) = vdemo2(T)
        t_vec = @elapsed res = sample(vdemo2kw(; T=Float64), alg, 250)

        @model function vdemo3(::Type{TV}=Vector{Float64}) where {TV<:AbstractVector}
            x = TV(undef, N)
            @. x ~ InverseGamma(2, 3)
        end

        sample(vdemo3(), alg, 250)
        sample(vdemo3(Vector{Float64}), alg, 250)

        vdemo3kw(; T) = vdemo3(T)
        sample(vdemo3kw(; T=Vector{Float64}), alg, 250)
    end
end


text:
# Utilities for demo.jl
module DemoTools

using LinearAlgebra
using TensorKit
using JLD2
using MERAKit

"""
    block_op(op, num_sites)

Take an operator `op` that defines the local term of a global operator, and block sites
together to form a new operator for which each site corresponds to `num_sites` sites of the
original operator, and which sums up to the same global operator. `num_sites` should be a
power of 2.
"""
function block_op(op::MERAKit.SquareTensorMap{2}, num_sites)
    while num_sites > 1
        VL = space(op, 1)
        VR = space(op, 2)
        eyeL = id(VL)
        eyeR = id(VR)
        opcross = eyeR ⊗ op ⊗ eyeL
        opleft = op ⊗ eyeL ⊗ eyeR
        opright = eyeL ⊗ eyeR ⊗ op
        op_new_unfused = opcross + 0.5*(opleft + opright)
        fusionspace = domain(op)
        fusetop = isomorphism(fuse(fusionspace), fusionspace)
        fusebottom = isomorphism(fusionspace, fuse(fusionspace))
        op = (fusetop ⊗ fusetop) * op_new_unfused * (fusebottom ⊗ fusebottom)
        num_sites /= 2
    end
    if num_sites != 1
        msg = "`num_sites` needs to be a power of 2"
        throw(ArgumentError(msg))
    end
    return op
end

function block_op(op::MERAKit.SquareTensorMap{1}, num_sites)
    while num_sites > 1
        V = space(op, 1)
        eye = id(V)
        opleft = op ⊗ eye
        opright = eye ⊗ op
        op_new_unfused = opleft + opright
        fusionspace = domain(op_new_unfused)
        fusetop = isomorphism(fuse(fusionspace), fusionspace)
        fusebottom = isomorphism(fusionspace, fuse(fusionspace))
        op = fusetop * op_new_unfused * fusebottom
        num_sites /= 2
    end
    if num_sites != 1
        msg = "`num_sites` needs to be a power of 2"
        throw(ArgumentError(msg))
    end
    return op
end

function goldenchain_hamiltonian(JI, Jτ; block_size = 1)
    V = GradedSpace[FibonacciAnyon](:I => 0, :τ => 1)
    H = TensorMap(zeros, Float64, V*V ← V*V)
    H.data[FibonacciAnyon(:I)] .= JI
    H.data[FibonacciAnyon(:τ)] .= Jτ
    H = block_op(H, block_size)
    H = H / block_size
    return H
end

"""
    XXZ_hamiltonian(J_xy, J_z; symmetry=:none, block_size=1)

Return the local Hamiltonian term for the XXZ model: J_xy*XX + J_xy*YY + J_z*ZZ.

`symmetry` should be `:none` or `:group`, and determmines whether the Hamiltonian should be
an explicitly U(1) symmetric `TensorMap` or a dense one. `block_size` determines how many
sites to block together, and should be a power of 2. If sites are blocked, then the
Hamiltonian is divided by a constant to make the energy per site be the same as for the
original.
"""
function XXZ_hamiltonian(J_xy, J_z; symmetry=:none, block_size=1)
    if symmetry == :U1 || symmetry == :group
        V = ℂ[U₁](-1=>1, 1=>1)
        Z = TensorMap(zeros, Float64, V ← V)
        Z.data[U₁(1)] .= 1.0
        Z.data[U₁(-1)] .= -1.0
        ZZ = Z ⊗ Z
        XXplusYY = TensorMap(zeros, Float64, V ⊗ V ← V ⊗ V)
        XXplusYY.data[U₁(0)] .= [0.0 2.0; 2.0 0.0]
    elseif symmetry == :none
        V = ℂ^2
        Z = TensorMap(zeros, Float64, V ← V)
        Z.data .= [1.0 0.0; 0.0 -1.0]
        ZZ = Z ⊗ Z
        XXplusYY = TensorMap(zeros, Float64, V ⊗ V ← V ⊗ V)
        XXplusYY.data[2,3] = 2.0
        XXplusYY.data[3,2] = 2.0
    else
        error("Unknown symmetry: $symmetry")
    end
    H = J_xy*XXplusYY + J_z*ZZ
    H = block_op(H, block_size)
    H = H / block_size
    return H
end

"""
    ising_hamiltonian(J=-1.0, h=1.0; symmetry=:none, block_size=1)

Return the local Hamiltonian term for the Ising model: J*XX + h*Z

`symmetry` should be `:none`, `:group`, or `:anyons` and determmines whether the Hamiltonian
should be an explicitly Z2 symmetric or anyonic `TensorMap`, or a dense one. `block_size`
determines how many sites to block together, and should be a power of 2. If sites are
blocked, then the Hamiltonian is divided by a constant to make the energy per site be the
same as for the original.
"""
function ising_hamiltonian(J=-1.0, h=1.0; symmetry=:none, block_size=1)
    if symmetry == :Z2
        V = ℂ[ℤ₂](0=>1, 1=>1)
        # Pauli Z
        Z = TensorMap(zeros, Float64, V ← V)
        Z.data[ℤ₂(0)] .= 1.0
        Z.data[ℤ₂(1)] .= -1.0
        eye = id(V)
        ZI = Z ⊗ eye
        IZ = eye ⊗ Z
        # Pauli XX
        XX = TensorMap(zeros, Float64, V ⊗ V ← V ⊗ V)
        XX.data[ℤ₂(0)] .= [0.0 1.0; 1.0 0.0]
        XX.data[ℤ₂(1)] .= [0.0 1.0; 1.0 0.0]
        H = J*XX + h/2 * (ZI+IZ)
    elseif symmetry == :anyons
        V = RepresentationSpace{IsingAnyon}(:I => 0, :ψ => 0, :σ => 1)
        H = TensorMap(zeros, Float64, V ⊗ V ← V ⊗ V)
        if J == 1.0 && h == 1.0
            H.data[IsingAnyon(:I)] .= 2.0
            H.data[IsingAnyon(:ψ)] .= -2.0
        elseif J == -1.0 && h == 1.0
            H.data[IsingAnyon(:I)] .= -2.0
            H.data[IsingAnyon(:ψ)] .= 2.0
        else
            msg = "Can't create an anyon symmetric Ising Hamiltonian that isn't critical."
            throw(ArgumentError(msg))
        end
    elseif symmetry == :none
        V = ℂ^2
        # Pauli matrices
        X = TensorMap(zeros, Float64, V ← V)
        Z = TensorMap(zeros, Float64, V ← V)
        eye = id(V)
        X.data .= [0.0 1.0; 1.0 0.0]
        Z.data .= [1.0 0.0; 0.0 -1.0]
        XX = X ⊗ X
        ZI = Z ⊗ eye
        IZ = eye ⊗ Z
        H = J*XX + h/2 * (ZI+IZ)
    else
        error("Unknown symmetry: $symmetry")
    end
    H = block_op(H, block_size)
    H = H / block_size
    return H
end

"""
    magop(; symmetry = :none, block_size = 1)

Return the magnetization operator for the Ising model, blocked over `block_size` sites.
"""
function magop(; symmetry = :none, block_size = 1)
    if symmetry == :none
        V = ℂ^2
    elseif symmetry == :group || symmetry == :Z2
        V = Z2Space(0 => 1, 1 => 1)
    else
        msg = "Unknown symmetry: $symmetry"
        raise(ArgumentError(msg))
    end
    X = TensorMap(zeros, Float64, V ← V)
    if symmetry == :none
        X.data .= [0.0 1.0; 1.0 0.0]
    end
    magop = block_op(X, block_size)
    return magop
end

"""
    store_mera(path, m)

Store the MERA `m` on disk at `path`.

The storage format is a JLD2 file of `pseudoserialize(m)`.  We use pseudoserialization to
both work around limitations of JLD(2) and improve compatibility between different versions
of MERAKit.jl

See also: [`load_mera`](@ref)
"""
function store_mera(path, m)
    deser = pseudoserialize(m)
    @save path deser
end

"""
    load_mera(path)

Load the MERA at stored at `path`.

The storage format is a JLD2 file of `pseudoserialize(m)`.  We use pseudoserialization to
both work around limitations of JLD(2) and improve compatibility between different versions
of MERAKit.jl

See also: [`store_mera`](@ref)
"""
function load_mera(path)
    @load path deser
    m = depseudoserialize(deser...)
    return m
end

end  # module


text:
__precompile__()

"""

# GeoEfficiency Package
introduce a fast and flexible tool to calculate in batch or individually the `geometrical efficiency`
for a set of common radiation detectors shapes (cylindrical,Bore-hole, Well-type) as seen form
a source. The source can be a point, a disc or even a cylinder.

# Quick Usage
*  geoEff()	: Calculate the geometrical efficiency for one geometrical setup return only the value of the geometrical efficiency.

*  calc() 	: Calculate the geometrical efficiency for one geometrical setup and display full information on the console.

*  calcN()	: Calculate the geometrical efficiency for geometrical setup(s) and display full information on the console until the user quit.

*  batch()	: Calculate the geometrical efficiency using data in the **`$(join(split(dataDir,"/travis")))`** folder in batch mode.

!!! note
    for more information and updates refer to the repository at [`GitHub.com`](https://github.com/DrKrar/GeoEfficiency.jl/)

"""
module GeoEfficiency

export
	about,

 # Config

 # Input_Console

 # Physics_Model
	Point,
	source,
	Detector,
	CylDetector,
	BoreDetector,
	WellDetector,

# Input_Batch
	getDetectors,
	setSrcToPoint,
	typeofSrc,

 # Calculations
	geoEff,

 # Output_Interface
 	max_batch,
 	calc,
	calcN,
	batch,
	batchInfo

include("Config.jl") # to overwrite defaults edit parameters; restore by comment out this line.
include("Error.jl")	# define error system for the package.
include("Input_Console.jl")
include("Physics_Model.jl")
include("Input_Batch.jl")
include("Calculations.jl")
include("Output_Interface.jl")

#------------------------ about ---------------------------

using Dates

const abt = """
\n
\t **************************************************
\t **            -=) GeoEfficiency (=-             **
\t **  Accurate Geometrical Efficiency Calculator  **
\t **   First Created on Fri Aug 14 20:12:01 2015  **
\t **************************************************

\t Author:        Mohamed E. Krar,  @e-mail: DrKrar@gmail.com
\t Auth_Profile:  https://www.researchgate.net/profile/Mohamed_Krar3
\t Repository:    https://github.com/DrKrar/GeoEfficiency.jl/
\t Version:       v"0.9.4-dev" - ($(Date(now()) - Date("2019-05-21")) old master)
\t Documentation: https://GeoEfficiency.GitHub.io/dev/index.html
\t PDF_Manual:    https://GeoEfficiency.GitHub.io/dev/GeoEfficiency.jl.pdf
\n
\n\tBatch Mode Calculations
\t-  read files by defaul from directory `$(join(split(dataDir, "/travis")))`
\t-  save results by default to directory `$(join(split(resultdir, "/travis")))`
\n\tfor more information see `batch`, `batchInfo`.
\n
"""

"""
$abt
"""
about() = printstyled(abt, color = :green, bold = true)
about()

end #module


text:
#---- Set up
println("\n\n\n Analysis begins:")
cd("/Users/yanivabir/Files/ShohamyLab/DecisionToExplore/exp4_replication/decExp4")
# Criteria
deck_ratio_deviation = 0.3
reward = 0.25

# Load libraries
using CSV, DataFrames, Statistics, Distributions

#---- Get data
date = r"2021-11-04|2021-11-05"
sess = "4"
subdir = "now"

# List relevant files
files = filter(x -> occursin(date, x) &&
    occursin("sess" * sess, x) &&
    !occursin("plan", x) &&
    !occursin("int", x) &&
    occursin("csv", x), readdir("../Data/"* subdir* "/"))

# Open each file, and then vcat them together
dat = ((x -> DataFrame!(CSV.File("../Data/" * subdir * "/" * x))).(files)) |> d -> vcat(d..., cols = :union)

#---- Exclude Ps who didn't finish

quality = combine(groupby(dat, :PID), :n_warnings => maximum => :n_w) # Summarize warnings each participant recieved

kick = combine(groupby(dat, :PID), :category => (x -> "kick-out" in x) => :kick)
filter!(x -> !ismissing(x.kick) && x.kick == true, kick)

filter!(r -> !(r.PID in kick.PID), dat) # Filter data according to quality table

println("\nExcluded during experiment:")
println(kick.PID)
#---- Prepare data
# Remove instructions and training
a = filter(r -> r.trial_type == "structure-quiz", dat) |>  # Find structure-quiz as beginning of interesting data
    x -> combine(groupby(x, [:PID, :sess]), :trial_index => minimum => :start_trial) # Find fist appearance of
dat = innerjoin(dat, a, on = [:PID, :sess]) |> # Join with data
    x -> filter!(r -> r.trial_index >= r.start_trial, x) # filter data

# Deal with missing block number data
for ii in 1:nrow(dat)
    if ismissing(dat[ii, :block])
        dat[ii, :block] = dat[ii - 1, :block]
    end
end

# Make deck names into numbers
dat.LCard = passmissing(x -> parse(Int64, x[8:end])).(dat.LCard)
dat.RCard = passmissing(x -> parse(Int64, x[8:end])).(dat.RCard)
dat.chosen_deck = passmissing(x -> parse(Int64, x[8:end])).(dat.chosen_deck)

# Make table names into numbers
dat.TL = passmissing(x -> parse(Int64, x[end])).(dat.TL)
dat.TR = passmissing(x -> parse(Int64, x[end])).(dat.TR)
dat.chosen_object = passmissing(x -> parse(Int64, x[end])).(dat.chosen_object)
dat.table = passmissing(x -> parse(Int64, x[end])).(dat.table)

# Which deck was chosen?
dat.chosen_deck =
    (x -> !ismissing(x.chosen_deck) ? x.chosen_deck : (ismissing(x.chosen_deck_side) ? missing :
        (x.chosen_deck_side == "L" ? x.LCard : x.RCard))).(eachrow(dat))

#---- Compute bonus payments
bonus = filter(x -> x.trial_type == "test-feedback", dat) |>
    x -> combine(groupby(x, :PID), :correct => (x -> sum((x .== "true") .| (x .== true))) => :correct)

bonus.reward = (x -> round(x, digits = 2)).(bonus.correct .* reward)
if sess == "4"
    bonus.reward = bonus.reward .+ 2
end

println("\nBonuses:")
println(sort(bonus, [:reward]))
if subdir != "raw"
    CSV.write("../Data/" * subdir * "/bonus.csv", sort(bonus, [:reward])[!, [:PID, :reward]], header = false)
end
#---- Exclude according to memory quiz
memory = filter(x -> x.trial_type == "memory-quiz", dat) # Select trials
memory = memory[:, filter(x -> !isequal(unique(memory[!, x]),
    [missing]), names(memory))] # Remove empty variables
memory.correct = memory.chosen_table .== memory.answer # Compute correct

memory = combine(groupby(memory, [:PID]), :correct => sum => :correct, :correct => length => :n)
memory.p = 1 .- cdf.(Binomial.(memory.n, 0.25), memory.correct)

pre = nrow(memory)
mem_exc = filter(x -> x.p > 0.05, memory)
filter!(x -> x.p < 0.05, memory)

println("\nRemoved ", pre - nrow(memory), " participants for low memory score")
println(mem_exc.PID)
println("\nAverage memory score remaining: ", mean(memory.correct ./ memory.n))

dat = filter(x -> x.PID in memory.PID, dat) # Exclude from data

#---- Exclude according to deck choices
# Separate learning phase
learn = filter(x -> !ismissing(x.category) &&
    x.category in ["object_choice", "deck_choice", "too-quick", "focus-blur"], dat)

# Remove aborted trials
filter!(x -> !(x.category == "object_choice" && x.rt == "null"), learn) # Remove no response object choice trials

bad_indx = []
for ii in nrow(learn):-1:1
    if learn[ii, :category] == "deck_choice" && learn[ii, :rt] == "null"
        push!(bad_indx, ii, ii - 1)
    elseif learn[ii, :category] == "too-quick" || learn[ii, :category] == "focus-blur"
        push!(bad_indx, ii, ii - 1)
        if learn[ii-1, :category] == "deck_choice"
            push!(bad_indx, ii - 2)
        end
    end
end

learn = learn[filter(x -> !(x in bad_indx), 1:nrow(learn)), :]

# Remove empty columns
learn = learn[:, filter(x -> !isequal(unique(learn[!, x]), [missing]), names(learn))]

# Make RT numeric, on seconds scale
learn.rt = (x -> parse(Float64, x)).(learn.rt) / 1000

# Transform dataframe to have one line per trial: merge table and deck choice lines
wl = filter(x -> x.category == "object_choice", learn)

rename!(wl, :rt => :table_rt)

wl[!, :rt] .= 0.

for ii in 1:nrow(wl)
    id = wl[ii, :PID]
    sess = wl[ii, :sess]
    bl = wl[ii, :block]
    ind = wl[ii, :trial_index]

    # Candidate row to pull data from - has same Participand ID, same block, matching trial index
    cand = filter(x -> x.PID == id && x.sess == sess && x.block == bl && x.trial_index == ind + 3, learn)
    if nrow(cand) > 0
        for jj in [:LCard, :RCard, :rt, :chosen_deck_side, :chosen_deck, :reward]
            wl[ii, jj] = cand[!, jj][1]
        end
    end
end

rename!(wl, :rt => :deck_rt);

wl.trial = reduce(vcat, [1:nrow(x) for x in groupby(wl, [:PID, :sess, :block])])

# Count which deck was chosen
decks = combine(groupby(wl, [:PID, :sess, :block, :chosen_object]),
    :chosen_deck => (x -> minimum(skipmissing(x))) => :min_deck)
wl = innerjoin(wl, decks, on = [:PID, :sess, :block, :chosen_object])
wl.chosen_deck_N = [x.chosen_deck == x.min_deck ? 1 : 2 for x in eachrow(wl)]

wl.running_table_n = [x.chosen_object + 4 * (x.block) for x in eachrow(wl)]

# For each subject and game, compute the absolute deviation from choosing both decks 50:50
deck_s = combine(groupby(wl, [:PID, :running_table_n]),
    :chosen_deck_N => (x -> abs(0.5 - mean(x .== 1))) => :ratio,
    :chosen_deck_N => length => :n_trials)

# Multiply each table's ratio by the number of trials for that table
deck_s.wrat = deck_s.ratio .* deck_s.n_trials

# Compute weighted average
deck_s = combine(groupby(deck_s, :PID), [:wrat, :n_trials] => ((a,b) -> sum(a) / sum(b)) => :wrat)

# Remove subjects by exclusion criterion
pre = nrow(deck_s)
deck_exc = filter(x -> x.wrat > deck_ratio_deviation, deck_s)
filter!(x -> x.wrat <= deck_ratio_deviation, deck_s)
println("\nRemoved ", pre - nrow(deck_s), " participants for low entropy deck choices")
println(deck_exc.PID)

# Remove subjects from data
filter!(x -> x.PID in deck_s.PID, dat)

# Print list of PIDs for next session
function completed_next(x)
    next_sess = filter(y -> occursin("sess" * string(parse(Int, sess) + 1), y) &&
        !occursin("plan", y) &&
        !occursin("int", y) &&
        occursin("csv", y), readdir("../Data/"* subdir* "/"))

    sum((y -> occursin(x, y)).(next_sess)) > 0
end

future = filter(x -> completed_next(x), unique(dat.PID))

println("\nValid Ps this session: n=", length(unique(dat.PID)))

println("\nPIDs for next session:")
next = filter(x -> !(x in future),unique(dat.PID))
println("\nn=", length(next))
(x -> println(x)).(next)


text:
function Kolmogorov_smirnov_test(dat::AbstractArray,d::ContinuousUnivariateDistribution,xmin::Number,xmax::Int64 = round(Int,1e5))
  data = sort(dat)
  n = float(length(data))
  max_indx = findlast(x -> x <= xmax, data)
  min_indx = findfirst(x -> x >= xmin, data)
  data = data[min_indx:max_indx]

  act_cdf = collect(0:length(data)-1) / n
  thr_cdf = cdf(d,float(data))
  D = maximum(abs(act_cdf-thr_cdf))
  return D
end

function Kolmogorov_smirnov_test(dat::AbstractArray,d::DiscreteUnivariateDistribution,xmin::Int64,xmax::Int64 = round(Int,1e5))
  alpha, xmin = params(d)
  data = round(Int,sort(dat))
  n = float(length(data))
  max_indx = findlast(x -> x <= xmax, data)
  min_indx = findfirst(x -> x >= xmin, data)
  data = data[min_indx:max_indx]

  thr_cdf = cdf(d,collect(xmin:data[end])+1)
  occurence = counts(data)[Int(xmin) - (data[1] - 1):end]
  act_cdf = occurence/sum(occurence)
  act_cdf = cumsum(act_cdf)
  D = maximum(abs(act_cdf-thr_cdf))
  return D
end

#helper function
function create_histogram(x::AbstractArray)
  h = zeros(typeof(x[1]),0)
  max = 0
  d = 0
  for i = 1:length(x)
    if max < x[i]
      for j=1:(x[i]-max)
        push!(h,0)
      end
      h[x[i]] =1
      max = x[i]
    else
      h[x[i]] +=1
    end
  end
  return h
end