text:
### copied from https://github.com/sswatson/Graphics2D.jl/blob/master/src/Graphics2D.jl
### at 08.12.2017 (commit: e6e2e1e7646d458b1db859721eb2c1253ced322b)

import Base.show,
       Base.*,
       Base.-,
       Cairo

COLORS = Dict(
"red" =>    [1.0, 0.0, 0.0, 1.0],
"green" =>  [0.0, 1.0, 0.0, 1.0],
"blue" =>   [0.0, 0.0, 1.0, 1.0],
"yellow" => [1.0, 1.0, 0.0, 1.0],
"purple" => [0.5, 0.0, 0.5, 1.0],
"orange" => [1.0, 0.4, 0.0, 1.0],
"white" =>  [1.0, 1.0, 1.0, 1.0],
"black" =>  [0.0, 0.0, 0.0, 1.0],
"gray" =>   [0.5, 0.5, 0.5, 1.0]
)

abstract type GraphicElement end

struct GraphicsBoundingBox
    xmin::Real
    xmax::Real
    ymin::Real
    ymax::Real
end

struct GraphicPoint <: GraphicElement
    x::Real
    y::Real
    pointsize::Real
    color::Array{Float64,1}
end

struct GraphicLine{T<:Real} <: GraphicElement
    coords::Array{T,2}
    linewidth::Real
    color::Array{Float64,1}
    fill::Bool
    fillcolor::Array{Float64,1}
end

struct GraphicArrow{T<:Real} <: GraphicElement
    coords::Array{T,2}
    linewidth::Real
    color::Array{Float64,1}
    arrowsize::Real
    arrowloc::Real
end

struct GraphicCircle <: GraphicElement
    center::Array{Float64,1}
    radius::Float64
    linewidth::Real
    color::Array{Float64,1}
    fill::Bool
    fillcolor::Array{Float64,1}
end

struct GraphicArc <: GraphicElement
    center::Array{Float64,1}
    radius::Real
    theta1::Real
    theta2::Real
    linewidth::Real
    color::Array{Float64,1}
    fill::Bool
    fillcolor::Array{Float64,1}
end

struct GraphicText <: GraphicElement
    location::Array{Float64,1}
    text::AbstractString
    textsize::Real
    color::Array{Float64,1}
end

*(x::Real,s::String) = x * COLORS[s]

tocolor(string_or_rgb) = isa(string_or_rgb,AbstractString) ?
                            COLORS[string_or_rgb] : string_or_rgb
function GraphicPoint(x::Real,
               y::Real;
               pointsize::Real=0.005,
               color=[0.0,0.0,0.0,1.0])
    return GraphicPoint(x,y,pointsize,tocolor(color))
end

function GraphicPoint(z::Complex{T}; args...) where T <: Real
    return GraphicPoint(real(z),imag(z); args...)
end

function GraphicPoint(p::Tuple{T,T}; args...) where T <: Real
    return GraphicPoint(p[1],p[2]; args...)
end

function GraphicPoint(a::Array{T,1}; args...) where T <: Real
    return GraphicPoint(a[1],a[2]; args...)
end

"""
    GraphicPoint(args...)

Draws a point in the gui.

# Optinal Arguments

- `pointsize::Real=0.005`
- `color=[0.0,0.0,0.0,1.0]`
"""
function GraphicPoint(; args...)
    return x -> GraphicPoint(x; args...)
end

function GraphicLine(A::Array{T,2};
                       linewidth::Real=1.0,
                       color=[0.0,0.0,0.0,1.0],
                       fill::Bool=false,
                       fillcolor=[1.0,1.0,1.0,1.0]) where T<:Real
    return GraphicLine(A,linewidth,tocolor(color),fill,tocolor(fillcolor))
end

function GraphicLine(a::Array{T,2};args...) where T
    return GraphicLine(map(float,a);args...)
end

function GraphicLine(a::Array{Complex{T},1};args...) where T
    return GraphicLine(hcat(real(a),imag(a));args...)
end

function GraphicLine(a::Array{Array{T,1},1};args...) where T
    return GraphicLine(hcat(T[a[k][1] for k=1:length(a)],T[a[k][2] for k=1:length(a)]);args...)
end

function GraphicLine(a::Array{Tuple{T,T},1};args...) where T
    return GraphicLine(hcat(T[a[k][1] for k=1:length(a)],T[a[k][2] for k=1:length(a)]);args...)
end

function GraphicLine(a::GraphicArrow)
    return GraphicLine(a.coords;linewidth=a.linewidth,color=a.color)
end

"""
    GraphicLine(args...)

Draws a line in the gui. Needs a two element vector with the start and the endpoint.

# Optinal Arguments

- `linewidth::Real=1.0`
- `color=[0.0,0.0,0.0,1.0]`
- `fill::Bool=false`
- `fillcolor=[1.0,1.0,1.0,1.0])`
"""
function GraphicLine(; args...)
    return x -> GraphicLine(x; args...)
end

function GraphicArrow(A::Array{T,2};
                        linewidth::Real=1.0,
                        color=[0.0,0.0,0.0,1.0],
                        arrowsize::Real=0.05,
                        arrowloc::Real=0.5) where T<:Real
    return GraphicArrow(A,linewidth,tocolor(color),arrowsize,arrowloc)
end

function GraphicArrow(a::Array{T,2};args...) where T
    return GraphicArrow(map(float,a);args...)
end

function GraphicArrow(a::Array{Complex{T},1};args...) where T
    return GraphicArrow(hcat(real(a),imag(a));args...)
end

function GraphicArrow(a::Array{Array{T,1},1};args...) where T
    return GraphicArrow(hcat(T[a[k][1] for k=1:length(a)],T[a[k][2] for k=1:length(a)]);args...)
end

function GraphicArrow(a::Array{Tuple{T,T},1};args...) where T
    return GraphicArrow(hcat(T[a[k][1] for k=1:length(a)],T[a[k][2] for k=1:length(a)]);args...)
end

"""
    GraphicArrow(args...)

Draws an arrow in the gui. Needs a two element vector with the start and the endpoint. The arrow head lies in the endpoint.

# Optinal Arguments

- `linewidth::Real=1.0`
- `color=[0.0,0.0,0.0,1.0]`
- `arrowsize::Real=0.05`
- `arrowloc::Real=0.5`
"""
function GraphicArrow(; args...)
    return x -> GraphicArrow(x; args...)
end

function -(a::GraphicArrow)
    return GraphicArrow(a.coords[end:-1:1,1:end],
                 a.linewidth,
		 a.color,
		 a.arrowsize,
		 a.arrowloc)
end

function *(i::Integer,a::GraphicArrow)
    return signbit(i) ? -a : a
end

function GraphicCircle(center::Array{T,1},
                     radius::U;
                     linewidth::Real=1.0,
                     color=[0.0,0.0,0.0, 1.0],
                     fill::Bool=false,
                     fillcolor=[1.0,1.0,1.0, 1.0]
                     ) where {T,U}
    return GraphicCircle(map(float, center), float(radius), linewidth,
        tocolor(color), fill, tocolor(fillcolor))
end

function GraphicCircle(x::T,
                     y::T,
                     radius::U;
                     args...
                     )  where {T,U}
    return GraphicCircle(T[x,y],radius;args...)
end

function GraphicCircle(p::Tuple{T,T},
                     radius::U;
                     args...
                     ) where {T,U}
    return GraphicCircle(p[1],p[2],radius;args...)
end

"""
    GraphicCircle(args...)

Draws a circle in the gui. Needs a tuple with the center and the radius.

# Optinal Arguments

- `linewidth::Real=1.0`
- `color=[0.0,0.0,0.0,1.0]`
- `fill::Bool=false`
- `fillcolor=[1.0,1.0,1.0,1.0])`
"""
function GraphicCircle(; args...)
    return (x, r) -> GraphicCircle(x, r; args...)
end

function GraphicArc(center::Array{T,1},
                radius,
                theta1,
                theta2;
                linewidth::Real=1.0,
                color=[0.0,0.0,0.0, 1.0],
                fill::Bool=false,
                fillcolor=[1.0,1.0,1.0, 1.0]
                ) where T
    return GraphicArc(map(float,center),
               float(radius),
               float(theta1),
               float(theta2),
               linewidth,
               tocolor(color),
               fill,
               tocolor(fillcolor))
end

function GraphicArc(; args...)
    return (center, radius, theta1, theta2) -> GraphicArc(center, radius, theta1, theta2; args...)
end

function GraphicText(location::Array{T,1},
                text::AbstractString;
                textsize::Float64=0.2,
                color=[1.0,1.0,1.0, 1.0]
                ) where T
    return GraphicText(location,text,textsize,tocolor(color))
end

"""
    GraphicText(args...)

Draws text in the gui. Needs a tuple with the startposition and the text.

# Optinal Arguments

- `textsize::Real=0.2`
- `color=[1.0,1.0,1.0, 1.0]`
"""
function GraphicText(; args...)
    return (x, text) -> GraphicText(x, text; args...)
end

function arrowhead(a::GraphicArrow)
    A = a.coords[end-1,:]
    B = a.coords[end,:]
    z = (B-A)[1] + im*((B-A)[2])
    θ = angle(z)
    ψ = 20π/180
    l = a.arrowsize
    tip = A + a.arrowloc*(B-A)
    C = tip-abs(z)*l*[cos(θ+ψ),sin(θ+ψ)]
    D = tip-abs(z)*l*[cos(θ-ψ),sin(θ-ψ)]

    return GraphicLine([tip,C,D,tip],fill=true,color=a.color,fillcolor=a.color)
end

function show(io::IO,line::GraphicLine)
    print(io,"GraphicLine(coordinates =\n")
    show(io,line.coords)
    print(io,",\nlinewidth = ",line.linewidth)
    if line.fill
        print(io,",\nfillcolor = ")
        print(io,"[",line.fillcolor[1],",",
              line.fillcolor[2],",",line.fillcolor[3],",",line.fillcolor[4],"]")
    else
        print(io,",\nfill = ",line.fill)
    end
    print(io,")")
end

function show(io::IO,arrow::GraphicArrow)
    print(io,"GraphicArrow(coordinates =\n")
    show(io,arrow.coords)
    print(io,",\nlinewidth = ",arrow.linewidth)
    print(io,",\narrowsize = ",arrow.arrowsize)
    print(io,",\narrowloc = ",arrow.arrowloc)
    print(io,")")
end

function show(io::IO,point::GraphicPoint)
    print(io,"GraphicPoint(")
    print(io,point.x)
    print(io,",")
    print(io,point.y)
    print(io,")")
end

function boundingbox(l::GraphicLine)
    # Coordinates are returned in the order xmin, xmax, ymax, ymin
    return GraphicsBoundingBox(minimum(l.coords[:,1]), maximum(l.coords[:,1]),
        minimum(l.coords[:,2]),maximum(l.coords[:,2]))
end

function boundingbox(a::GraphicArrow)
    # Coordinates are returned in the order xmin, xmax, ymax, ymin
    return boundingbox(map(boundingbox,[GraphicLine(a),arrowhead(a)]))
end

function boundingbox(c::GraphicCircle)
    return GraphicsBoundingBox(c.center[1] - c.radius, c.center[1] + c.radius,
        c.center[2] - c.radius, c.center[2] + c.radius)
end

function boundingbox(p::GraphicPoint)
    return GraphicsBoundingBox(p.x,p.x,p.y,p.y)
end

function boundingbox(a::GraphicArc)
    return GraphicsBoundingBox(a.center[1] - a.radius,
        a.center[1] + a.radius, a.center[2] - a.radius, a.center[2] + a.radius)
end

function boundingbox(t::GraphicText)
    return GraphicsBoundingBox(t.location[1],t.location[1],t.location[2],
        t.location[2])
end

function boundingbox(A::Array)
    a = minimum([bb.xmin for bb in A])
    b = maximum([bb.xmax for bb in A])
    c = maximum([bb.ymax for bb in A])
    d = minimum([bb.ymin for bb in A])
    if a == b || c == d
        error(logger, "Bounding box has empty interior")
    end
    return GraphicsBoundingBox(a - 0.05*(b-a), b + 0.05*(b-a), d - 0.05*(c-d),
        c + 0.05*(c-d))
end

function sketch(cr::Cairo.CairoContext,
                l::GraphicLine,
                bb::GraphicsBoundingBox=boundingbox(l))

    Cairo.set_line_width(cr,l.linewidth)

    Cairo.move_to(cr,l.coords[1,1],l.coords[1,2]);
    for i=1:size(l.coords)[1]-1
        Cairo.line_to(cr,l.coords[i+1,1],l.coords[i+1,2]);
    end

    if l.coords[end,:] == l.coords[1,:]
        Cairo.close_path(cr)
    end

    if l.coords[end,:] == l.coords[1,:] && l.fill
        Cairo.set_source_rgba(cr,l.fillcolor...)
        Cairo.fill_preserve(cr)
        Cairo.set_source_rgba(cr,l.color...)
        Cairo.stroke(cr)
    else
        Cairo.set_source_rgba(cr,l.color...)
        Cairo.stroke(cr)
    end
end

function sketch(cr::Cairo.CairoContext,
                a::GraphicArrow,
                bb::GraphicsBoundingBox=boundingbox(a))

    Cairo.set_line_width(cr,a.linewidth)

    Cairo.move_to(cr,a.coords[1,1],a.coords[1,2]);
    for i=1:size(a.coords)[1]-1
        Cairo.line_to(cr,a.coords[i+1,1],a.coords[i+1,2]);
    end

    Cairo.set_source_rgba(cr,a.color...)
    Cairo.stroke(cr)

    sketch(cr,arrowhead(a),bb)

end

function sketch(cr::Cairo.CairoContext,
                c::GraphicCircle,
                bb::GraphicsBoundingBox=boundingbox(c))

    Cairo.set_line_width(cr,c.linewidth)
    Cairo.arc(cr, c.center[1], c.center[2], c.radius, 0, 2*pi)

    if c.fill
        Cairo.set_source_rgba(cr,c.fillcolor...)
        Cairo.fill_preserve(cr)
    end

    Cairo.set_source_rgba(cr,c.color...)
    Cairo.stroke(cr)
end

function sketch(cr::Cairo.CairoContext,
                a::GraphicArc,
                bb::GraphicsBoundingBox=boundingbox(a))
    Cairo.arc(cr, a.center[1], a.center[2], a.radius, a.theta1, a.theta2);
    Cairo.set_source_rgba(cr,a.color...)
    Cairo.stroke(cr)
end

function sketch(cr::Cairo.CairoContext,
                p::GraphicPoint,
                bb::GraphicsBoundingBox=boundingbox(p))

    Cairo.arc(cr, p.x, p.y,
              p.pointsize*max(bb.xmax-bb.xmin,bb.ymax-bb.ymin), 0, 2*pi);
    Cairo.set_source_rgba(cr,p.color...)
    Cairo.fill(cr)
end

function sketch(cr::Cairo.CairoContext,
                t::GraphicText,
                bb::GraphicsBoundingBox=boundingbox(t))

    Cairo.save(cr)

    Cairo.set_font_size(cr,t.textsize)
    Cairo.set_source_rgba(cr,t.color...)

    x_bearing, y_bearing, width, height = Cairo.text_extents(cr,t.text)

    Cairo.move_to(cr,t.location[1] - width/2 - x_bearing, t.location[2] + (3/2)*height + y_bearing)

    # Cairo.scale(cr,1,-1)
    # Cairo.translate(cr, 0.0, -height)
    Cairo.show_text(cr,t.text)
    Cairo.restore(cr)
end


text:
"""
    License(; name="MIT", path=nothing, destination="LICENSE")

Creates a license file.

## Keyword Arguments
- `name::AbstractString`: Name of a license supported by PkgTemplates.
  Available licenses can be seen
  [here](https://github.com/invenia/PkgTemplates.jl/tree/master/templates/licenses).
- `path::Union{AbstractString, Nothing}`: Path to a custom license file.
  This keyword takes priority over `name`.
- `destination::AbstractString`: File destination, relative to the repository root.
  For example, `"LICENSE.md"` might be desired.
"""
struct License <: BasicPlugin
    path::String
    destination::String
end

function License(;
    name::AbstractString="MIT",
    path::Union{AbstractString, Nothing}=nothing,
    destination::AbstractString="LICENSE",
)
    if path === nothing
        path = default_file("licenses", name)
        isfile(path) || throw(ArgumentError("License '$(basename(path))' is not available"))
    end
    return License(path, destination)
end

source(p::License) = p.path
destination(p::License) = p.destination
view(::License, t::Template, ::AbstractString) = Dict(
    "AUTHORS" => join(t.authors, ", "),
    "YEAR" => year(today()),
)

ext: jl


text:
######################### Gradient Computation #######################
#TODO: LSQ Gradient
# Non-functional
function leastSqGrad(mesh::Mesh, matrix::AbstractArray{Float64, 2}, stencil=zeros(2,2))
    # Stencil should be a list of lists, with each sublist containing the cells contained in the stencil of the main cell

    #cells, cVols, cCenters, faces, fAVecs, fCenters, boundaryFaces = mesh
    #bdryFaceIndices = Array(nFaces-nBdryFaces:nFaces)

    nCells, nFaces, nBoundaries, nBdryFaces = unstructuredMeshInfo(mesh)
    nVars = size(matrix, 2) # Returns the length of the second dimension of "matrix"

    grad = zeros(nCells, nVars, 3)

    L11 = zeros(nCells,nVars)
    L12 = zeros(nCells,nVars)
    L13 = zeros(nCells,nVars)
    L22 = zeros(nCells,nVars)
    L23 = zeros(nCells,nVars)
    L33 = zeros(nCells,nVars)
    L1f = zeros(nCells,nVars)
    L2f = zeros(nCells,nVars)
    L3f = zeros(nCells,nVars)

    @fastmath for f in 1:nFaces-nBdryFaces

        c1 = mesh.faces[f][1]
        c2 = mesh.faces[f][2]

        dx = mesh.cCenter[c2][1] - mesh.cCenter[c1][1]
        dy = mesh.cCenter[c2][2] - mesh.cCenter[c1][2]
        dz = mesh.cCenter[c2][3] - mesh.cCenter[c1][3]

        weight = 1 / sqrt(dx*dx + dy*dy + dz*dz)

        wdx = weight * dx
        wdy = weight * dy
        wdz = weight * dz

        for v in 1:nVars
            dv = matrix[c2,v] - matrix[c1,v]
            wdv = weight * dv

            L11[c1,v] += wdx^2
            L12[c1,v] += wdx * wdy
            L13[c1,v] += wdx * wdz
            L22[c1,v] += wdy^2
            L23[c1,v] += wdy * wdz
            L33[c1,v] += wdz^2

            L1f[c1,v] += wdx * wdv
            L2f[c1,v] += wdy * wdv
            L3f[c1,v] += wdz * wdv

            L11[c2,v] += wdx^2
            L12[c2,v] += wdx * wdy
            L13[c2,v] += wdx * wdz
            L22[c2,v] += wdy^2
            L23[c2,v] += wdy * wdz
            L33[c2,v] += wdz^2

            L1f[c2,v] += wdx * wdv
            L2f[c2,v] += wdy * wdv
            L3f[c2,v] += wdz * wdv

        end

    end

    # Deal with boundary faces, and add them to the matrix vectors

    return grad
end

#=
    Takes the gradient of (scalar) data provided in matrix form (passed into arugment 'matrix'):
    Cell      x1      x2      x3
    Cell 1    x1_1    x2_1    x3_1
    Cell 2    x1_2    x2_2    x3_2
    ...
    (where x1, x2, and x3 are arbitrary scalars)
    and output a THREE-DIMENSIONAL gradient arrays of the following form
    Cell      x1          x2          x3
    Cell 1    grad(x1)_1  grad(x2)_1  grad(x3)_1
    Cell 2    grad(x1)_2  grad(x2)_2  grad(x3)_2
    ...
    Where each grad(xa)_b is made up of THREE elements for the (x,y,z) directions

    Ex. Gradient @ cell 1 of P would be: greenGaussGrad(mesh, P)[1, 1, :]
=#
function greenGaussGrad(mesh::Mesh, matrix::AbstractArray{Float64, 2}, valuesAtFaces=false)
    nCells, nFaces, nBoundaries, nBdryFaces = unstructuredMeshInfo(mesh)
    nVars = size(matrix, 2)

    #Interpolate values to faces if necessary
    if valuesAtFaces
        faceVals = matrix
    else
        faceVals = linInterp_3D(mesh, matrix)
    end

    # Create matrix to hold gradients
    grad = zeros(nCells, nVars, 3)

    # Integrate fluxes from each face
    @fastmath for f in eachindex(mesh.faces)
        ownerCell = mesh.faces[f][1]
        neighbourCell = mesh.faces[f][2]

        @inbounds for v in 1:nVars
            for d in 1:3
                # Every face has an owner
                grad[ownerCell, v, d] += mesh.fAVecs[f][d] * faceVals[f, v]

                # Boundary faces don't - could split into two loops
                if neighbourCell > -1
                    grad[neighbourCell, v, d] -= mesh.fAVecs[f][d] * faceVals[f, v]
                end
            end
        end
    end

    # Divide integral by cell volume to obtain gradients
    @simd for c in 1:nCells
        for v in 1:nVars
            for d in 1:3
                grad[c,v,d] /= mesh.cVols[c]
            end
        end
    end

    return grad
end

####################### Face value interpolation ######################
#=
    Interpolates to all INTERIOR faces
    Arbitrary value matrix interpolation
    Example Input Matrix:
    Cell      x1       x2       x3
    Cell 1    x1_c1    x2_c1    x3_c1
    Cell 2    x1_c2    x2_c2    x3_c2
    ...
    Outputs a matrix of the following form
    Cell      x1       x2       x3
    Face 1    x1_f1    x2_f1    x3_f1
    Face 2    x1_f2    x2_f2    x3_f2
    ...
=#
function linInterp_3D(mesh::Mesh, matrix::AbstractArray{Float64, 2}, faceVals=nothing)
    nCells, nFaces, nBoundaries, nBdryFaces = unstructuredMeshInfo(mesh)
    nVars = size(matrix, 2)

    # Make a new matrix if one is not passed in
    if faceVals == nothing
        faceVals = zeros(nFaces, nVars)
    end

    # Boundary face fluxes must be set separately
    for f in 1:nFaces-nBdryFaces
        # Find value at face using linear interpolation
        c1 = mesh.faces[f][1]
        c2 = mesh.faces[f][2]

        #TODO: Precompute these distances?
        c1Dist = 0
        c2Dist = 0
        for i in 1:3
            c1Dist += (mesh.cCenters[c1][i] - mesh.fCenters[f][i])^2
            c2Dist += (mesh.cCenters[c2][i] - mesh.fCenters[f][i])^2
        end
        totalDist = c1Dist + c2Dist

        for v in 1:nVars
            faceVals[f, v] = matrix[c1, v]*(c2Dist/totalDist) + matrix[c2, v]*(c1Dist/totalDist)
        end
    end

    return faceVals
end

# Similar to linInterp_3D. Instead of linearly interpolating, selects the maximum value of the two adjacent cell centers as the face value
function maxInterp(mesh::Mesh, vars...)
    nCells, nFaces, nBoundaries, nBdryFaces = unstructuredMeshInfo(mesh)

    nVars = size(vars, 1)
    faceVals = zeros(nFaces, nVars)

    # Boundary face fluxes must be set separately
    @inbounds @fastmath for f in 1:nFaces-nBdryFaces
        c1 = mesh.faces[f][1]
        c2 = mesh.faces[f][2]

        for v in 1:nVars
            faceVals[f, v] = max(vars[v][c1], vars[v][c2])
        end
    end

    return faceVals
end

# Similar to linInterp_3D. Instead of linearly interpolating, calculates the change (delta) of each variable across each face (required for JST method)
function faceDeltas(mesh::Mesh, sln::SolutionState)
    nCells, nFaces, nBoundaries, nBdryFaces = unstructuredMeshInfo(mesh)
    nVars = size(sln.cellState, 2)

    faceDeltas = zeros(nFaces, nVars)

    # Boundary face fluxes must be set separately (faceDelta is zero at all possible boundary conditions right now)
    @inbounds @fastmath for f in 1:nFaces-nBdryFaces
        ownerCell = mesh.faces[f][1]
        neighbourCell = mesh.faces[f][2]

        for v in 1:nVars
            faceDeltas[f, v] = sln.cellState[neighbourCell, v] - sln.cellState[ownerCell, v]
        end
    end

    return faceDeltas
end

######################### TimeStepping #################################
# Calculate sln.cellPrimitives and sln.cellFluxes from sln.cellState
function decodeSolution_3D(sln::SolutionState, fluid::Fluid)
    nCells = size(sln.cellState, 1)
    for c in 1:nCells
        # Updates cell primitives
        @views decodePrimitives3D!(sln.cellPrimitives[c,:], sln.cellState[c,:], fluid)
        # Updates mass, xMom, eV2 x,y,z-direction fluxes
        @views calculateFluxes3D!(sln.cellFluxes[c, :], sln.cellPrimitives[c,:], sln.cellState[c,:])
    end
end

#=
    Calculates sln.fluxResiduals from sln.faceFluxes
=#
function integrateFluxes_unstructured3D(mesh::Mesh, sln::SolutionState)
    nCells, nFaces, nBoundaries, nBdryFaces = unstructuredMeshInfo(mesh)

    # Recomputing flux balances, so wipe existing values
    fill!(sln.fluxResiduals, 0)
    nVars = size(sln.fluxResiduals, 2)

    #### Flux Integration ####
    @inbounds @fastmath for f in eachindex(mesh.faces)
        ownerCell = mesh.faces[f][1]
        neighbourCell = mesh.faces[f][2]

        for v in 1:nVars
            i1 = (v-1)*3 + 1
            i2 = i1+2
            @views flow = dot(sln.faceFluxes[f, i1:i2], mesh.fAVecs[f])

            # Subtract from owner cell
            sln.fluxResiduals[ownerCell, v] -= flow
            # Add to neighbour cell
            if neighbourCell > -1
                sln.fluxResiduals[neighbourCell, v] += flow
            end
        end
    end

    # Divide by cell volume
    for c in 1:nCells
        for v in 1:nVars
            sln.fluxResiduals[c,v] /= mesh.cVols[c]
        end
    end

    return sln.fluxResiduals
end

max_stars_repo_name: henrystoldt/fvCFD


text:
module PkgTest

# include("gaussian_vc_benchmark.jl")
# include("dyestuff.jl")
# include("dyestuff2.jl")
# include("sleepstudy.jl")
include("randtest.jl")

end


text:
function read_kpts( filname )
    file = open(filname)
    str = readline(file)
    Nkpts = parse( Int, str )
    kpts = zeros( Float64, 3,Nkpts )
    for ik = 1:Nkpts
        str = split(readline(file))
        kpts[1,ik] = parse( Float64, str[1] )
        kpts[2,ik] = parse( Float64, str[2] )
        kpts[3,ik] = parse( Float64, str[3] )
    end
    close(file)
    return kpts
end


text:
using ShortStrings
using Test, Random

# write your own tests here
r = randstring.(1:15)
@test all(ShortString15.(r) .== r)
a = ShortString15.(r)
@test fsort(a) |> issorted

r = randstring.(1:7)
@test all(ShortString7.(r) .== r)
a = ShortString7.(r)
@test fsort(a) |> issorted

r = randstring.(1:3)
@test all(ShortString{UInt32}.(r) .== r)
a = ShortString{UInt32}.(r)
@test fsort(a) |> issorted

@test all(ShortString3.(r) .== r)
a = ShortString3.(r)
@test fsort(a) |> issorted

@test collect(ShortString3("abc")) == ['a', 'b', 'c']


text:
module Interfaces

VERSION < v"0.4-" && using Docile

#MUTABLE_UNION_LIB_PATH = "deps/usr/lib/lib_mutable_union.so"
FILE_PATH = @__FILE__
PKG_PATH = dirname(FILE_PATH)
MUTABLE_UNION_LIB = Libdl.dlopen("$PKG_PATH/../deps/usr/lib/lib_mutable_union.so")

export @interface, @implements, implements, methods_exist, MUTABLE_UNION_LIB

include("core.jl")

end


text:
using Documenter, CoupledFields, Cairo

makedocs(
    modules = [CoupledFields],
    clean = false,
    sitename = "CoupledFields.jl",
    pages = [
        "Home" => "index.md",
        "Guide" => ["man/Example1.md", "man/Example2.md"],
        "Library" => "lib/library.md"
    ],
)

deploydocs(
    repo = "github.com/Mattriks/CoupledFields.jl.git",
)

max_stars_repo_name: Mattriks/CoupledFields

max_forks_repo_name: Mattriks/CoupledFields


text:
######################################################################
# counting various spell types
######################################################################
using AMDB
# you may need to Pkg.add these...
using DataStructures
using GZip
using UnicodePlots

# you only need to read the data once per session
records = GZip.open(deserialize, "/Data/AMDB/subsample/AMDB_subsample.jls.gz", "r")

"""
Demo for counting the distribution of the number of person spells.
"""
function count_person_lengths(records)
    c = counter(Int)
    for (id,record) in records
        push!(c, length(record.person_spells))
    end
    sort(collect(c), by=first)
end

count_person_lengths(records)   # looks sensible

# a plot for the number of AMP spells
UnicodePlots.histogram(filter(l->l ≤ 100,
                              [length(r.AMP_spells) for (_,r) in records]))

# which individuals have above 1000 spells?
filter((_,r) -> length(r.AMP_spells) > 1000, records)


text:
using LoopVectorization
"""
    Fermi.ConfigurationInteraction.DetOperations.Determinant(α::String, β::String)

Constructor function for Determinant object using strings.
"""
function Determinant(α::String, β::String)

    αint = parse(Int, reverse(α); base=2)
    βint = parse(Int, reverse(β); base=2)

    Determinant(αint, βint)
end

"""
    Fermi.ConfigurationInteraction.DetOperations.αlist(D::Determinant)

Return the alpha string of a Determinant as a list
"""
function αlist(D::Determinant)

    return [parse(Int, ss) for ss in reverse(bitstring(D.α))]
end

"""
    Fermi.ConfigurationInteraction.DetOperations.βlist(D::Determinant)

Return the beta string of a Determinant as a list
"""
function βlist(D::Determinant)

    return [parse(Int, ss) for ss in reverse(bitstring(D.β))]
end

"""
    Fermi.ConfigurationInteraction.DetOperations.αindexi!(D::Determinant, Out::Array{Int64,1})

Write the indexes of the occupid alpha electrons of the Determinant to a given list
"""
function αindex!(D::Determinant, Out::Array{Int64,1})

    i = 1
    e = 1

    # Loop until 'e' electrons are found. Be careful! If 'e' is greater than
    # the number of electrons you will get stuck!
    while e ≤ length(Out)
        if 1<<(i-1) & D.α ≠ 0
            Out[e] = i
            e += 1
        end
        i += 1
    end
end

"""
    Fermi.ConfigurationInteraction.DetOperations.βindex!(D::Determinant, Out::Array{Int64.1})

Write the indexes of the occupid beta electrons of the Determinant to a given list
"""
function βindex!(D::Determinant, Out::Array{Int64,1})

    i = 1
    e = 1

    # Loop until 'e' electrons are found. Be careful! If 'e' is greater than
    # the number of electrons you will get stuck!
    while e ≤ length(Out)
        if 1<<(i-1) & D.β ≠ 0
            Out[e] = i
            e += 1
        end
        i += 1
    end
end

"""
    Fermi.ConfigurationInteraction.DetOperations.αexcitation_level(D1::Determinant, D2::Determinant)

Compare two determinants to return the excitation level of the alpha electrons
"""
@inline function αexcitation_level(D1::Determinant, D2::Determinant;bail=false)
    count_ones(D1.α ⊻ D2.α) >> 1
end

"""
    Fermi.ConfigurationInteraction.DetOperations.βexcitation_level(D1::Determinant, D2::Determinant)

Compare two determinants to return the excitation level of the beta electrons
"""
@inline function βexcitation_level(D1::Determinant, D2::Determinant;bail=false)
    count_ones(D1.β ⊻ D2.β) >> 1
end

"""
    Fermi.ConfigurationInteraction.DetOperations.excitation_level(D1::Determinant, D2::Determinant)

Compare two determinants to return the excitation level between them
"""
@inline function excitation_level(D1::Determinant, D2::Determinant)
    return αexcitation_level(D1, D2) + βexcitation_level(D1,D2)
end

"""
    Fermi.ConfigurationInteraction.DetOperations.αexclusive(D1::Determinant, D2::Determinant)

Compare two determinants to return the indexes of the alpha electrons present in the first, but
not in the second
"""
function αexclusive(D1::Determinant, D2::Determinant)

    αexcl = D1.α ⊻ D2.α & D1.α

    out = []
    i = 1
    # Save alphas exclusives, in crescent order
    while 1<<((i-1)&63) ≤ αexcl
        if 1<<((i-1)&63) & αexcl ≠ 0
            push!(out, i)
        end
        i += 1
    end

    return out
end

"""
    Fermi.ConfigurationInteraction.DetOperations.first_αexclusive(D1::Determinant, D2::Determinant)

Compare two determinants to return the index of the first alpha electron present in the first, but
not in the second
"""
function first_αexclusive(D1::Determinant, D2::Determinant)

    x = D1.α ⊻ D2.α & D1.α
    return 64 - leading_zeros(x&(~(x-1)))
end

"""
    Fermi.ConfigurationInteraction.DetOperations.second_αexclusive(D1::Determinant, D2::Determinant)

Compare two determinants to return the index of the second alpha electron present in the second, but
not in the second
"""
function second_αexclusive(D1::Determinant, D2::Determinant)

    x = D1.α ⊻ D2.α & D1.α
    y = 64 - leading_zeros(x&(~(x-1)))
    x = x &(~(2^(y-1)))
    w = 64 - leading_zeros(x&(~(x-1)))
end

"""
    Fermi.ConfigurationInteraction.DetOperations.βexclusive(D1::Determinant, D2::Determinant)

Compare two determinants to return the indexes of the beta electrons present in the first, but
not in the second
"""
function βexclusive(D1::Determinant, D2::Determinant)

    βexcl = D1.β ⊻ D2.β & D1.β

    out = []
    i = 1
    # Save betas exclusives, in crescent order
    while 1<<(i-1) ≤ βexcl
        if 1<<(i-1) & βexcl ≠ 0
            push!(out, i)
        end
        i += 1
    end
    return out
end

"""
    Fermi.ConfigurationInteraction.DetOperations.first_βexclusive(D1::Determinant, D2::Determinant)

Compare two determinants to return the index of the first beta electron present in the first, but
not in the second
"""
function first_βexclusive(D1::Determinant, D2::Determinant)

    x = D1.β ⊻ D2.β & D1.β
    return 64 - leading_zeros(x&(~(x-1)))

end

"""
    Fermi.ConfigurationInteraction.DetOperations.second_βexclusive(D1::Determinant, D2::Determinant)

Compare two determinants to return the index of the second beta electron present in the second, but
not in the second
"""
function second_βexclusive(D1::Determinant, D2::Determinant)

    x = D1.β ⊻ D2.β & D1.β
    y = 64 - leading_zeros(x&(~(x-1)))
    x = x &(~(2^(y-1)))
    w = 64 - leading_zeros(x&(~(x-1)))
end

"""
    Fermi.ConfigurationInteraction.DetOperations.exclusive(D1::Determinant, D2::Determinant)

Returns a list with tuples (orbital index, spin) of orbitals populated in the first, but
not in the second Determinant.
"""
function exclusive(D1::Determinant, D2::Determinant)

    αexcl = D1.α ⊻ D2.α & D1.α
    βexcl = D1.β ⊻ D2.β & D1.β

    out = []
    i = 1
    # Save alphas exclusives, in crescent order
    while 1<<((i-1)) ≤ αexcl
        if 1<<((i-1)) & αexcl ≠ 0
            push!(out, (i, 'α'))
        end
        i += 1
    end
    i = 1
    # Save betas exclusives, in crescent order
    while 1<<((i-1)) ≤ βexcl
        if 1<<((i-1)) & βexcl ≠ 0
            push!(out, (i, 'β'))
        end
        i += 1
    end
    return out
end

"""
    Fermi.ConfigurationInteraction.DetOperations.annihilate(D::Determinant, orb::Int, spin::Char)

Creates a copy of the determinat where an electron from a specified orbital (orb) and spin (spin) was deleted. Returns the phase
factor along with the new determinant
"""
function annihilate(D::Determinant, orb::Int, spin::Char)

    if spin == 'α'
        #if D.α & (1 << (orb-1)) == 0
        #    error("Annihilation error. Orbital $orb is not occupied")
        #end

        # Determine sign
        l = 0
        i = 1
        while i < (1 << ((orb-1)))
            l += D.α & i ≠ 0 ? 1 : 0
            i = i << 1
        end

        newα = D.α ⊻ (1 << ((orb-1)))

        return (-1)^l, Determinant(newα, D.β)

    elseif spin =='β'
        #if D.β & (1 << (orb-1)) == 0
        #    error("Annihilation error. Orbital $orb is not occupied")
        #end

        # Determine sign
        l = count(i->(i=='1'), bitstring(D.α))
        i = 1
        while i < (1 << (orb-1))
            l += D.β & i ≠ 0 ? 1 : 0
            i = i << 1
        end

        newβ = D.β ⊻ (1 << (orb-1))

        return (-1)^l, Determinant(D.α, newβ)
    end
end

"""
    Fermi.ConfigurationInteraction.DetOperations.create(D::Determinant, orb::Int, spin::Char)

Creates a copy of the determinat where an electron was added to a specified orbital (orb) with a given spin (spin).
Returns the phase factor along with the new determinant
"""
function create(D::Determinant, orb::Int, spin::Char)

    if spin == 'α'
        #if D.α & (1 << (orb-1)) ≠ 0
        #    error("Creation error. Orbital $orb is occupied")
        #end

        # Determine sign
        l = 0
        i = 1
        while i < (1 << (orb-1))
            l += D.α & i ≠ 0 ? 1 : 0
            i = i << 1
        end

        newα = D.α | (1 << (orb-1))

        return (-1)^l, Determinant(newα, D.β)

    elseif spin =='β'
        #if D.β & (1 << (orb-1)) ≠ 0
        #    error("Annihilation error. Orbital $orb is not occupied")
        #end

        # Determine sign
        l = count(i->(i=='1'), bitstring(D.α))
        i = 1
        while i < (1 << (orb-1))
            l += D.β & i ≠ 0 ? 1 : 0
            i = i << 1
        end

        newβ = D.β | (1 << (orb-1))

        return (-1)^l, Determinant(D.α, newβ)
    end
end

#function α1phase

"""
    Fermi.ConfigurationInteraction.DetOperations.phase(D1::Determinant, D2::Determinant)

Returns the phase (+1 or -1) associated with the transformation of one determinant
into the other through second quantization operators
"""
function phase(D1::Determinant, D2::Determinant)

    p = 1
    _det = Determinant(D1.α, D1.β)

    # For a string of excitation operators: abc...kji. We apply the annihilation operations such that k < j < i
    # Thus, the reverse function
    for (i,σ) in reverse(exclusive(D1, D2))
        f, _det = annihilate(_det, i, σ)
        p = f*p
    end
    # For the creation operations a > b > c. Thus, no reverse.
    for (a,σ) in exclusive(D2, D1)
        f, _det = create(_det, a, σ)
        p = f*p
    end

    return p
end

"""
    Fermi.ConfigurationInteraction.DetOperations.phase(D::Determinant, l::Int=0)

Prints alpha and beta strings of a Determinant with a given length(l)
"""
function showdet(D::Determinant, l::Int = 0)

    if l == 0
        l = length(bitstring(D.α))
    end
    println("α: "*reverse(bitstring(D.α))[1:l])
    println("β: "*reverse(bitstring(D.β))[1:l])

end

function αocc!(D::Determinant, R::UnitRange{Int64}, Out::Array{Int64,1})

    e = 1
    for i in R
        if 1<<(i-1) & D.α != 0
            Out[e] = i
            e += 1
        end
    end

end

function βocc!(D::Determinant, R::UnitRange{Int64}, Out::Array{Int64,1})

    e = 1
    for i in R
        if 1<<(i-1) & D.β != 0
            Out[e] = i
            e += 1
        end
    end

end

function αvir!(D::Determinant, R::UnitRange{Int64}, Out::Array{Int64,1})

    e = 1
    for i in R
        if 1<<(i-1) & D.α == 0
            Out[e] = i
            e += 1
        end
    end

end

function βvir!(D::Determinant, R::UnitRange{Int64}, Out::Array{Int64,1})

    e = 1
    for i in R
        if 1<<(i-1) & D.β == 0
            Out[e] = i
            e += 1
        end
    end
end

function detstring(D::Determinant,l::Int=64)

    αstr = ""
    βstr = ""
    str = ""

    i = 1
    for i = 1:l
        if 1<<(i-1) & (D.α & D.β) != 0
            αstr *= "1"
            βstr *= "1"
            str  *= "2"
        elseif 1<<(i-1) & D.α != 0
            αstr *= "1"
            βstr *= "0"
            str  *= "+"
        elseif 1<<(i-1) & D.β != 0
            αstr *= "0"
            βstr *= "1"
            str  *= "-"
        else
            αstr *= "0"
            βstr *= "0"
            str  *= "0"
        end
        i += 1
    end

    return "$str     $αstr     $βstr"
end


text:
# import Random
import Random: RandomDevice
import RandomNumbers: AbstractRNG

"""
```julia
gen_seed(T[, n])
```

Generate a tuple of `n` truly random numbers in type `T`. If `n` is missing, return only one number.
The "truly" random numbers are provided by the random device of system. See
[`Random.RandomDevice`](https://github.com/JuliaLang/julia/blob/master/base/random.jl#L29).

# Examples
```julia
julia> RandomNumbers.gen_seed(UInt64, 2)  # The output should probably be different on different computers.
(0x26aa3fe5e306f725,0x7b9dc3c227d8acc9)

julia> RandomNumbers.gen_seed(UInt32)
0x9ba60fdc
```
"""
gen_seed(::Type{T}) where {T<:Number} = rand(RandomDevice(), T)
gen_seed(::Type{T}, n) where {T<:Number} = tuple(rand(RandomDevice(), T, n)...)

"Get the original output type of a RNG."
@inline output_type(::AbstractRNG{T}) where {T} = T

# For 1.0
# Random.gentype(::Type{T}) where {T<:AbstractRNG} = output_type(T)

"Get the default seed type of a RNG."
@inline seed_type(r::AbstractRNG) = seed_type(typeof(r))

@inline split_uint(x::UInt128, ::Type{UInt64}=UInt64) = (x % UInt64, (x >> 64) % UInt64)
@inline split_uint(x::UInt128, ::Type{UInt32}) = (x % UInt32, (x >> 32) % UInt32,
    (x >> 64) % UInt32, (x >> 96) % UInt32)
@inline split_uint(x::UInt64, ::Type{UInt32}=UInt32) = (x % UInt32, (x >> 32) % UInt32)
@inline union_uint(x::NTuple{2, UInt32}) = unsafe_load(Ptr{UInt64}(pointer(collect(x))), 1)
@inline union_uint(x::NTuple{2, UInt64}) = unsafe_load(Ptr{UInt128}(pointer(collect(x))), 1)
@inline union_uint(x::NTuple{4, UInt32}) = unsafe_load(Ptr{UInt128}(pointer(collect(x))), 1)

@inline function unsafe_copyto!(r1::R, r2::R, ::Type{T}, len) where {R, T}
    arr1 = Ptr{T}(pointer_from_objref(r1))
    arr2 = Ptr{T}(pointer_from_objref(r2))
    for i = 1:len
        unsafe_store!(arr1, unsafe_load(arr2, i), i)
    end
    r1
end

@inline function unsafe_compare(r1::R, r2::R, ::Type{T}, len) where {R, T}
    arr1 = Ptr{T}(pointer_from_objref(r1))
    arr2 = Ptr{T}(pointer_from_objref(r2))
    all(unsafe_load(arr1, i) == unsafe_load(arr2, i) for i in 1:len)
end


text:
"""
A `QuadProjectionSampling` is an operator that maps a function to its inner products
with a projection basis, with those inner products computed by a quadrature strategy.
"""
struct QuadProjectionSampling <: BasisFunctions.SamplingOperator
    dict		::  Dictionary
	src_space	::	BasisFunctions.FunctionSpace
	quad		::	QuadratureStrategy
end

QuadProjectionSampling(dict::Dictionary, quad::QuadratureStrategy) = QuadProjectionSampling(dict, Span(dict), quad)

BasisFunctions.name(op::QuadProjectionSampling) = "Projection by quadrature operator"

BasisFunctions.dictionary(op::QuadProjectionSampling) = op.dict

BasisFunctions.dest(op::QuadProjectionSampling) = dictionary(op)

BasisFunctions.src_space(op::QuadProjectionSampling) = op.src_space

quad(op::QuadProjectionSampling) = op.quad

BasisFunctions.apply!(result, op::QuadProjectionSampling, f) =
	quadproject!(result, f, dictionary(op), quad(op))

function quadproject!(result, f, dict::Dictionary1d, quad::QuadratureStrategy)
	@assert size(result) == size(dict)
    for i in eachindex(dict)
		result[i] = projectionintegral(quad, f, dict, i)
	end
	result
end


text:
"""
    barsandstripes(nsamples, nvariables)
Generates a test data set. To see the structure in the data set, run e. g.
`reshape(barsandstripes(1, 16), 4, 4)` a few times.

Example from:
MacKay, D. (2003). Information Theory, Inference, and Learning Algorithms
"""
function barsandstripes(nsamples::Int, nvariables::Int)
   squareside = sqrt(nvariables)
   if squareside != floor(squareside)
      error("Number of variables must be a square number")
   end
   squareside = round(Int, squareside)
   x = zeros(nsamples, nvariables)
   for i in 1:nsamples
      sample = hcat([(rand() < 0.5 ? ones(squareside) : zeros(squareside))
            for j in 1:squareside]...)
      if rand() < 0.5
         sample = transpose(sample)
      end
      x[i,:] .= sample[:]
      fill!(sample, 0.0)
   end
   x
end

"""
    batchparallelized(f, n, op)
Distributes the work for executing the function `f` `n` times
on all the available workers and reduces the results with the operator `op`.
`f` is a function that gets a number (of tasks) to execute the tasks.

# Example:
    batchparallelized(n -> aislogimpweights(dbm; nparticles = n), 100, vcat)
"""
function batchparallelized(f::Function, n::Int, op::Function)
   batches = mostevenbatches(n)
   if length(batches) > 1
      return @sync @distributed (op) for batch in batches
         f(batch)
      end
   else
      return f(n)
   end
end

"""
    crossvalidation(x, monitoredfit; ...)
Performs k-fold cross-validation, given
* the data set `x` and
* `monitoredfit`: a function that fits and evaluates a model.
  As arguments it must accept:
  - a training data data set
  - a `DataDict` containing the evaluation data.

The return values of the calls to the `monitoredfit` function
are concatenated with `vcat`.
If the monitoredfit function returns `Monitor` objects,
`crossvalidation` returns a combined `Monitor` object that can be displayed
by creating a cross-validation plot via
`BoltzmannMachinesPlots.crossvalidationplot`.

# Optional named argument:
* `kfold`: specifies the `k` in "`k`-fold" (defaults to 10).

    crossvalidation(x, monitoredfit, pars; ...)
If additionaly a vector of parameters `pars` is given, `monitoredfit`
also expects an additional parameter from the parameter set.
"""
function crossvalidation(x::Matrix{Float64}, monitoredfit::Function, pars...;
      kfold::Int = 10)

   vcat(pmap(monitoredfit, crossvalidationargs(x, pars...)...)...)
end

"""
    crossvalidationargs(x, pars...; )

Returns a tuple of argument vectors containing the parameters for a function
such as the `monitoredfit` argument in `crossvalidation`.

Usage example:
    map(monitoredfit, crossvalidationargs(x)...)

# Optional named argument:
* `kfold`: see `crossvalidation`.
"""
function crossvalidationargs(x::Matrix{Float64}, pars ...; kfold::Int = 10)

   nsamples = size(x, 1)
   batchranges = ranges(mostevenbatches(nsamples, kfold))

   args_data = Vector{Tuple{Matrix{Float64}, DataDict}}(undef, kfold)
   for i in eachindex(batchranges)
      rng = batchranges[i]
      trainingdata = x[[!(j in rng) for j in 1:nsamples], :]
      evaluationdata = x[rng, :]
      datadict = DataDict(string(i) => evaluationdata)
      args_data[i] = (trainingdata, datadict)
   end

   # this implementation is inefficient and overly complicated

   if isempty(pars)
      args = map(arg -> (arg,), args_data)
   else
      args = (args_data, map(p -> collect(p), pars)...)
      args = vec(collect(Iterators.product(args...)))
   end

   # vector of tuples to tuples of vectors
   (
      map(arg -> arg[1][1], args), # x
      map(arg -> arg[1][2], args), # datadict
      [map(arg -> arg[i], args) for i in 2:length(args[1])]...
   )
end

"""
A function with no arguments doing nothing.
Usable as default argument for functions as arguments.W
"""
function emptyfunc
end

"""
    log1pexp(x)
Calculates log(1+exp(x)). For sufficiently large values of x, the approximation
log(1+exp(x)) ≈ x is used. This is useful to prevent overflow.
"""
function log1pexp(x::Float64)
   if x > 34
      return x
   else
      return log1p(exp(x))
   end
end

function logit(p::Float64)
   -log.(1.0 ./ p .- 1.0)
end

"""
    mostevenbatches(ntasks)
    mostevenbatches(ntasks, nbatches)
Splits a number of tasks `ntasks` into a number of batches `nbatches`.
The number of batches is by default `min(nworkers(), ntasks)`.
The returned result is a vector containing the numbers of tasks for each batch.
"""
function mostevenbatches(ntasks::Int, nbatches::Int = min(nworkers(), ntasks))

   minparticlesperbatch, nbatcheswithplus1 = divrem(ntasks, nbatches)
   batches = Vector{Int}(undef, nbatches)
   for i = 1:nbatches
      if i <= nbatcheswithplus1
         batches[i] = minparticlesperbatch + 1
      else
         batches[i] = minparticlesperbatch
      end
   end
   batches
end

"""
    curvebundles(...)
Generates an example dataset that can be visualized as bundles of trend curves
with added noise. Additional binary columns with labels may be added.

## Optional named arguments:
* `nbundles`: number of bundles
* `nperbundle`: number of sequences per bundle
* `nvariables`: number of variables in the sequences
* `noisesd`: standard deviation of the noise added on all sequences
* `addlabels`: add leading columns to the resulting dataset, specifying the
   membership to a bundle
* `pbreak`: probability that an intermediate point in a sequence is a
   breakpoint, defaults to 0.2.
* `breakval`: a function that expects no input and generates a single
  (random) value for a defining point of a piecewise linear sequence.
  Defaults to `rand`.

# Example:
To quickly grasp the idea, plot generated samples against the variable index:

    x = BMs.curvebundles(nvariables = 10, nbundles = 3,
                       nperbundle = 4, noisesd = 0.03,
                       addlabels = true)
    BoltzmannMachinesPlots.plotcurvebundles(x)
"""
function curvebundles(;
      nbundles::Int = 3,
      nperbundle::Int = 50,
      nvariables::Int = 10,
      noisesd::Float64 = 0.05,
      addlabels::Bool = false,
      pbreak::Float64 = 0.2,
      breakval::Function = rand)

   x = piecewiselinearsequences(nbundles, nvariables;
         pbreak = pbreak, breakval = breakval)
   x = repeat(x, nperbundle)
   nsamples = size(x, 1)
   x .+= noisesd * randn(nsamples, nvariables)

   if addlabels
      nlabelvars = Int(ceil(log(nbundles)/log(2)))
      labels = Matrix{Float64}(nbundles, nlabelvars)
      labels[1,:] .= 0.0
      for j = 2:nbundles
         labels[j,:] .= labels[j-1,:]
         next!(view(labels, j, :))
      end
      labels = repeat(labels, nperbundle)
      x = hcat(labels, x)
   end

   x = x[randperm(nsamples), :]
   x
end

"""
    piecewiselinearsequences(nsequences, nvariables; ...)
Generates a dataset consisting of samples with values that
are piecewise linear functions of the variable index.

Optional named arguments: `pbreak`, `breakval`,
see `piecewiselinearsequencebundles`.
"""
function piecewiselinearsequences(nsequences::Int, nvariables::Int;
      pbreak::Float64 = 0.2, breakval::Function = rand)

   inbetweenvariables = 2:(nvariables-1)
   x = Matrix{Float64}(undef, nsequences, nvariables)
   for i = 1:nsequences
      breakpointindexes = [1; randsubseq(inbetweenvariables, pbreak); nvariables]
      breakpointvalues = [breakval() for i = 1:length(breakpointindexes)]
      lastbreakpoint = 1
      x[i, breakpointindexes] .= breakpointvalues
      for breakpoint in breakpointindexes[2:end]
         valdiff = (x[i, breakpoint] - x[i, lastbreakpoint]) /
               (breakpoint - lastbreakpoint)
         for j = (lastbreakpoint + 1):(breakpoint - 1)
            x[i, j] = x[i, j - 1] + valdiff
         end
         lastbreakpoint = breakpoint
      end
   end
   x
end

"""
    splitdata(x, ratio)
Splits the data set `x` randomly in two data sets `x1` and `x2`, such that
the ratio `n2`/`n1` of the numbers of lines/samples in `x1` and `x2` is
approximately equal to the given `ratio`.
"""
function splitdata(x::Matrix{Float64}, ratio::Float64)
   nsamples = size(x,1)
   testidxs = randperm(nsamples)[1:(round(Int, nsamples*ratio))]
   xtest = x[testidxs,:]
   trainidxs = setdiff(1:nsamples, testidxs)
   xtraining = x[setdiff(1:nsamples, testidxs),:]
   xtraining, xtest, trainidxs, testidxs
end

"""
    DBindex(x, cl; p, q)
Computes the Davies-Bouldin Index for cluster validation.
"""
function DBindex(x::Array{Float64,2}, cl::Array{Int64,1}; p::Int = 1, q::Int = 2)

   n = size(cl,1)
   k = sort(unique(cl))

   centers = Array{Float64,2}(undef,(size(k,1),size(x,2)))
   for i = 1:size(k,1)
      for j = 1:size(x,2)
         centers[i,j] =  mean(x[cl .== k[i], j])
      end
   end

   # S is the within-cluster dispersion
   S = zeros(size(k,1))
   for i = 1:size(k,1)
      ind = (cl .== k[i])
      if sum(ind) > 1
         centerI = centers[i,:]
         centerI = reshape(repeat(centerI', sum(ind)), (sum(ind), size(x,2)))
         tmp = (x[ind,:] .- centerI).^2
         S[i] = mean(sqrt.(mapslices(sum, tmp, dims = 2)).^q).^(1/q)
      else
         S[i] = 0.0
      end
   end

   # R is the between-cluster similarity as a function of the two within-cluster dispersions
   # M is the distance between cluster centroids
   M = pairwise(Euclidean(), centers, dims=1)
   R = Array{Float64,2}(undef,(size(k,1),size(k,1)))
   r = zeros(size(k,1))
   for i = 1:size(k,1)
      for j = 1:size(k,1)
         R[i,j] = (S[i] + S[j]) / M[i,j]
      end
      r[i] = maximum(R[i,:][isfinite.(R[i,:])])
   end

   DB = mean(r[isfinite.(r)])
   DB
end

# Number of parameters for Negative Binomial DBM

function nparams(dbm::Array{AbstractXBernoulliRBM,1})

   params = zeros()
   for i in 1:size(dbm,1)
      if  typeof(dbm[i]) == NegativeBinomialBernoulliRBM
         tmp = (size(dbm[i].weights,1) * size(dbm[i].weights,2)) + size(dbm[i].visbias,1) + size(dbm[i].hidbias,1) + size(dbm[i].inversedispersion,1)
         params .+= tmp
      else
         tmp = (size(dbm[i].weights,1) * size(dbm[i].weights,2)) + size(dbm[i].visbias,1) + size(dbm[i].hidbias,1)
         params .+= tmp
      end
   end
   params[1]
end


text:
facts("Miscellaneous") do

    context("Results storage") do
        results_storage = Dict()
        results_storage[new_processing_key(results_storage, "FTest")] = 4
        results_storage[new_processing_key(results_storage, "Turtle")] = 5
        results_storage[new_processing_key(results_storage, "FTest")] = 49

        @fact new_processing_key(results_storage, "FTest") --> "FTest3"

        @fact find_keys_containing(results_storage, "FTest") --> [1, 3]
        @fact find_keys_containing(results_storage, "Mum") --> []

    end

    context("File parts") do

        a, b, c = fileparts("")
        @fact a --> ""
        @fact b --> ""
        @fact c --> ""

        a, b, c = fileparts("/Users/test/subdir/test-file.bdf")
        @fact a --> "/Users/test/subdir/"
        @fact b --> "test-file"
        @fact c --> "bdf"

        a, b, c = fileparts("/Users/test/subdir/test_file.bdf")
        @fact a --> "/Users/test/subdir/"
        @fact b --> "test_file"
        @fact c --> "bdf"

        a, b, c = fileparts("test-file.bdf")
        @fact a --> ""
        @fact b --> "test-file"
        @fact c --> "bdf"
    end

    context("Find closest number index") do
        @fact _find_closest_number_idx([1, 2, 2.7, 3.2, 4, 3.1, 7], 3) --> 6
    end

end


text:
function gcd(x,y)
    if x == 0
        return y
    else
        return gcd(y%x,x)
    end
end
function find(arr,n)
    r = arr[1]
    for i=2:n
        r = gcd(arr[i],r)
    end
    return r
end
v = [8 16 8 4]
len = length(v)
println("The GCD of the array is = ",find(v,len))

max_stars_repo_name: landrumb/julia

max_issues_repo_name: landrumb/julia

max_forks_repo_name: landrumb/julia


text:
@testset "Cglobal" begin
	lib = Clibrary()

	val = Cglobal{Ptr{Cvoid}}(lib, :jl_nothing)
	@test val.handle != C_NULL
	@test val[] == unsafe_load(cglobal(:jl_nothing, Ptr{Cvoid}))
end


text:
# Velocity distribution plot in a large PIC domain.
#
#
# Hongyang Zhou, hyzhou@umich.edu 01/30/2020

using Batsrus, PyPlot, Printf, LinearAlgebra, Statistics
# For precise colorbar control
using PyCall
axes_grid1 = pyimport("mpl_toolkits.axes_grid1")
inset_locator = pyimport("mpl_toolkits.axes_grid1.inset_locator")
inset_axes = inset_locator.inset_axes

## Parameters
using Batsrus.UnitfulBatsrus: me, mp
const cAlfven = 253       # average Alfven speed in G8, [km/s]
const mi = 14             # average ion mass [amu]

"""
    dist_select(fnameParticle; ParticleType='e', dir=".")

Select particle in box regions centered at `xC`, `yC`, and `zC` with lengths `xL`, `yL` and
`zL`. `ParticleType` can be ['e','i'].
"""
function dist_select(fnameParticle, xC=-1.90, yC=0.0, zC=-0.1, xL=0.005, yL=0.2, zL=0.07;
   ParticleType='e', dir=".")

   if ParticleType == 'i'
      !occursin("region0_2", fnameParticle) && @error "Check filename!"
   elseif ParticleType == 'e'
      !occursin("region0_1", fnameParticle) && @error "Check filename!"
   end

   fnameField = "3d_var_region0_0_"*fnameParticle[end-22:end]

   nBox = 9 # number of box regions

   # Classify particles based on locations
   region = Array{Float32,2}(undef,6,nBox)
   region[:,1] = [xC-xL*3/2, xC-xL/2,   yC-yL/2, yC+yL/2, zC+zL/2,   zC+zL*3/2]
   region[:,2] = [xC-xL/2,   xC+xL/2,   yC-yL/2, yC+yL/2, zC+zL/2,   zC+zL*3/2]
   region[:,3] = [xC+xL/2,   xC+xL*3/2, yC-yL/2, yC+yL/2, zC+zL/2,   zC+zL*3/2]
   region[:,4] = [xC-xL*3/2, xC-xL/2,   yC-yL/2, yC+yL/2, zC-zL/2,   zC+zL/2]
   region[:,5] = [xC-xL/2,   xC+xL/2,   yC-yL/2, yC+yL/2, zC-zL/2,   zC+zL/2]
   region[:,6] = [xC+xL/2,   xC+xL*3/2, yC-yL/2, yC+yL/2, zC-zL/2,   zC+zL/2]
   region[:,7] = [xC-xL*3/2, xC-xL/2,   yC-yL/2, yC+yL/2, zC-zL*3/2, zC-zL/2]
   region[:,8] = [xC-xL/2,   xC+xL/2,   yC-yL/2, yC+yL/2, zC-zL*3/2, zC-zL/2]
   region[:,9] = [xC+xL/2,   xC+xL*3/2, yC-yL/2, yC+yL/2, zC-zL*3/2, zC-zL/2]

   particle = [Array{Float32}(undef, 3, 0) for _ in 1:nBox]

   data = readdata(fnameParticle; dir)

   x = @view data.x[:,:,:,1]
   y = @view data.x[:,:,:,2]
   z = @view data.x[:,:,:,3]

   ux_ = findfirst(x->x=="ux", data.head.wnames)
   uy_ = findfirst(x->x=="uy", data.head.wnames)
   uz_ = findfirst(x->x=="uz", data.head.wnames)

   ux = @view data.w[:,:,:,ux_]
   uy = @view data.w[:,:,:,uy_]
   uz = @view data.w[:,:,:,uz_]

   for ip in eachindex(x)
      for iR = 1:nBox
         if region[1,iR] < x[ip] < region[2,iR] &&
            region[3,iR] < y[ip] < region[4,iR] &&
            region[5,iR] < z[ip] < region[6,iR]

            particle[iR] = hcat(particle[iR], [ux[ip]; uy[ip]; uz[ip]])
            break
         end
      end
   end
   region, particle
end

"""
    dist_scan(region, particle, ParticleType='i', plottype=1)

Velocity distribution plot in 9 box regions.
`plottype`: 1: uy-ux; 2: ux-uz; 3:uy-uz; 4:u⟂O-u⟂I; 5:u⟂I-u∥; 5:u⟂O-u∥
"""
function dist_scan(region, particle, ParticleType='i', plottype=1; dir=".",
   fnameField::AbstractString, nbin=60, fs=10)

   if ParticleType == 'i'
      binRange = [[-3.,3.], [-3.,3.]]
   elseif ParticleType == 'e'
      binRange = [[-10.,12.], [-10.,10.]]
   end

   nBox = 9
   figure(figsize=(11,6))
   for iB = 1:nBox
      if plottype ≤ 3
         ux = particle[iB][1,:] ./ cAlfven
         uy = particle[iB][2,:] ./ cAlfven
         uz = particle[iB][3,:] ./ cAlfven
      else
         dBx, dBy, dBz = GetMeanField(fnameField, region[:,iB]; dir)

         dPar = [dBx; dBy; dBz] # Parallel direction
         dPerpI = cross([0; -1; 0], dPar) # Perpendicular direction in-plane
         dPerpO = cross(dPar, dPerpI) # Perpendicular direction out-of-plane

         uPar = transpose(particle[iB][1:3,:])*dPar ./ cAlfven
         uPerpI = transpose(particle[iB][1:3,:])*dPerpI ./ cAlfven
         uPerpO = transpose(particle[iB][1:3,:])*dPerpO ./ cAlfven
      end

      ax = subplot(3,4,iB+ceil(iB/3))
      if plottype == 1
         h = hist2D(uy, ux, bins=nbin,
            norm=matplotlib.colors.LogNorm(),density=true, range=binRange)
      elseif plottype == 2
         h = hist2D(ux, uz, bins=nbin,
            norm=matplotlib.colors.LogNorm(),density=true, range=binRange)
      elseif plottype == 3
         h = hist2D(uy, uz, bins=nbin,
            norm=matplotlib.colors.LogNorm(),density=true, range=binRange)
      elseif plottype == 4
         h = hist2D(uPerpO, uPerpI, bins=nbin,
            norm=matplotlib.colors.LogNorm(),density=true, range=binRange)
      elseif plottype == 5
         h = hist2D(uPerpI, uPar, bins=nbin,
            norm=matplotlib.colors.LogNorm(),density=true, range=binRange)
      elseif plottype==6
         h = hist2D(uPerpO, uPar, bins=nbin,
            norm=matplotlib.colors.LogNorm(),density=true, range=binRange)
      else
         @error "Unknown PlotVType!"
      end
      grid(true)
      axis("equal")

      if iB in (1,4,7,8,9)
         if plottype == 1
            xlabel(L"u_y",fontsize=fs)
            ylabel(L"u_x",fontsize=fs)
         elseif plottype == 2
            xlabel(L"u_x",FontSize=fs)
            ylabel(L"u_z",FontSize=fs)
         elseif plottype == 3
            xlabel(L"u_y",FontSize=fs)
            ylabel(L"u_z",FontSize=fs)
         elseif plottype == 4
            xlabel(L"u_{\perp Out}",fontsize=fs)
            ylabel(L"u_{\perp In}",fontsize=fs)
         elseif plottype == 5
            xlabel(L"u_{\perp In}",FontSize=fs)
            ylabel(L"u_\parallel",FontSize=fs)
         elseif plottype == 6
            xlabel(L"u_{\perp Out}",FontSize=fs)
            ylabel(L"u_\parallel",FontSize=fs)
         end
      end
      title(@sprintf("%d, x[%3.3f,%3.3f], z[%3.3f,%3.3f]",iB,region[1,iB],
         region[2,iB],region[5,iB],region[6,iB]))
      colorbar()
      plt.set_cmap("jet")
      #clim(1e-2,10^0.3)

      if ParticleType == 'e'
         str = "electron"
      elseif ParticleType == 'i'
         str = "ion"
      end
      text(0.05,0.05,str, FontSize=fs, transform=ax.transAxes)
   end

end

"""
    GetMeanField(fnameField, limits; dir=".")

Get the average field direction in region `limits`.
"""
function GetMeanField(fnameField, limits; dir=".")

   # Get the average field direction in limited region
   data = readdata(fnameField; dir)

   x = data.x[:,:,:,1]
   y = data.x[:,:,:,2]
   z = data.x[:,:,:,3]

   bx_ = findfirst(x->x=="Bx", data.head.wnames)
   by_ = findfirst(x->x=="By", data.head.wnames)
   bz_ = findfirst(x->x=="Bz", data.head.wnames)

   Bx = @view data.w[:,:,:,bx_]
   By = @view data.w[:,:,:,by_]
   Bz = @view data.w[:,:,:,bz_]

   xnew, ynew, znew, BxNew, ByNew, BzNew = subvolume(x,y,z, Bx,By,Bz, limits)

   # Average over the selected volume
   B̄x, B̄y, B̄z = mean(BxNew), mean(ByNew), mean(BzNew)

   # Normalize vector
   Length = √(B̄x^2 + B̄y^2 + B̄z^2)
   dBx, dBy, dBz = B̄x/Length, B̄y/Length, B̄z/Length

   dBx, dBy, dBz
end

function plotExCut(fnameField::String, region, xC, yC, zC, xL, yL, zL;
   dir=".", fs=16, sequence=129)

   plotrange = [xC-xL*16, xC+xL*16, zC-zL*5, zC+zL*5]
   # Sample region plot over contour
   data = readdata(fnameField, dir=dir)

   bx_ = findfirst(x->x=="Bx", data.head.wnames)
   bz_ = findfirst(x->x=="Bz", data.head.wnames)

   Bx = @view data.w[:,:,:,bx_]
   Bz = @view data.w[:,:,:,bz_]

   subplot(3,4,(1,9))
   cutplot(data, "Ex"; cut="y", sequence, plotrange)
   colorbar()
   axis("scaled")
   plt.set_cmap("RdBu_r")
   clim(-9e4,9e4)

   streamslice(data, "Bx;Bz"; cut="y", sequence, plotrange, color="k", density=1.0)

   xlabel(L"x [R_G]", fontsize=fs)
   ylabel(L"z [R_G]", fontsize=fs)
   title(L"Ex [\mu V/m]")

   nBox = 9
   for iB = 1:nBox
      rect = matplotlib.patches.Rectangle( (region[1,iB], region[5,iB]),
      region[2,iB]-region[1,iB], region[6,iB]-region[5,iB],
      ec="r", lw=1.2, fill=false) # facecolor="none"
      ax = gca()
      ax.add_patch(rect)
   end

end

function dist_plot(pType='e')
   dir = "."
   nBox = 4
   nbin = 60
   fs = 10
   subplotlist = [1,2,3,4,1,2,3,4]
   plottype = (1,1,1,1,3,3,3,3)

   fnameE = "cut_particles0_region0_1_t00001640_n00020369.out"
   fnameI = "cut_particles1_region0_2_t00001640_n00020369.out"

   fnameField = "3d_var_region0_0_"*fnameE[end-22:end]

   binRangeI = [[-3.,3.], [-3.,3.]]
   binRangeE = [[-7.,13.], [-10.,10.]]

   # Classify particles based on locations
   region = Array{Float32,2}(undef,6,nBox)
   particle = [Array{Float32}(undef, 3, 0) for _ in 1:nBox]

   # Electron
   if pType == 'e'
      #=
      region[:,1] = [-1.950, -1.945, -0.08, 0.08, -0.10, -0.06]
      region[:,2] = [-1.915, -1.910, -0.08, 0.08, -0.10, -0.06]
      region[:,3] = [-1.905, -1.900, -0.08, 0.08, -0.10, -0.06]
      region[:,4] = [-1.895, -1.890, -0.08, 0.08, -0.10, -0.06]
      =#
      region[:,1] = [-1.930, -1.925, -0.08, 0.08, -0.30, -0.26]
      region[:,2] = [-1.930, -1.925, -0.08, 0.08, -0.20, -0.16]
      region[:,3] = [-1.930, -1.925, -0.08, 0.08, -0.10, -0.06]
      region[:,4] = [-1.930, -1.925, -0.08, 0.08,  0.00,  0.04]

      data = readdata(fnameE, dir=dir)

      x = @view data.x[:,:,:,1]
      y = @view data.x[:,:,:,2]
      z = @view data.x[:,:,:,3]

      ux_ = findfirst(x->x=="ux", data.head.wnames)
      uy_ = findfirst(x->x=="uy", data.head.wnames)
      uz_ = findfirst(x->x=="uz", data.head.wnames)

      ux = @view data.w[:,:,:,ux_]
      uy = @view data.w[:,:,:,uy_]
      uz = @view data.w[:,:,:,uz_]

      for ip in eachindex(x)
         for iR = 1:nBox
            if region[1,iR] < x[ip] < region[2,iR] &&
               region[3,iR] < y[ip] < region[4,iR] &&
               region[5,iR] < z[ip] < region[6,iR]

               particle[iR] = hcat(particle[iR], [ux[ip]; uy[ip]; uz[ip]])
               break
            end
         end
      end

      binRange = binRangeE
   elseif pType == 'i'
      #=
      region[:,1] = [-1.865, -1.855, -0.08, 0.08, -0.29, -0.25]
      region[:,2] = [-1.910, -1.900, -0.08, 0.08, 0.11, 0.15]
      region[:,3] = [-1.970, -1.960, -0.08, 0.08, -0.10, -0.06]
      region[:,4] = [-1.850, -1.840, -0.08, 0.08, -0.10, -0.06]
      =#
      region[:,1] = [-1.930, -1.920, -0.08, 0.08, -0.30, -0.26]
      region[:,2] = [-1.930, -1.920, -0.08, 0.08, -0.20, -0.16]
      region[:,3] = [-1.930, -1.920, -0.08, 0.08, -0.10, -0.06]
      region[:,4] = [-1.930, -1.920, -0.08, 0.08,  0.00,  0.04]

      data = readdata(fnameI, dir=dir)

      x = @view data.x[:,:,:,1]
      y = @view data.x[:,:,:,2]
      z = @view data.x[:,:,:,3]

      ux_ = findfirst(x->x=="ux", data.head.wnames)
      uy_ = findfirst(x->x=="uy", data.head.wnames)
      uz_ = findfirst(x->x=="uz", data.head.wnames)

      ux = @view data.w[:,:,:,ux_]
      uy = @view data.w[:,:,:,uy_]
      uz = @view data.w[:,:,:,uz_]

      for ip in eachindex(x)
         for iR = 1:nBox
            if region[1,iR] < x[ip] < region[2,iR] &&
               region[3,iR] < y[ip] < region[4,iR] &&
               region[5,iR] < z[ip] < region[6,iR]

               particle[iR] = hcat(particle[iR], [ux[ip]; uy[ip]; uz[ip]])
               break
            end
         end
      end

      binRange = binRangeI
   end

   # Normalized quantities
   fig, ax = plt.subplots(4,2,figsize=(4.2,8.0))
   if pType == 'e'
      fig.suptitle("electron", fontsize=14, fontweight="bold")
   elseif pType == 'i'
      fig.suptitle("ion", fontsize=14, fontweight="bold")
   end

   xlPos = (0.50, -0.10)
   ylPos = (-0.25, 0.50)
   c = Vector{PyObject}(undef,nBox)
   plt.set_cmap("jet")

   for i = 1:nBox*2
      ax[i].tick_params(which="both", direction="in", top=true, right=true)
      ax[i].minorticks_on()
      iB = subplotlist[i]
      if plottype[i] ≤ 3
         ux = particle[iB][1,:] ./ cAlfven
         uy = particle[iB][2,:] ./ cAlfven
         uz = particle[iB][3,:] ./ cAlfven
      else
         dBx, dBy, dBz = GetMeanField(fnameField, region[:,iB]; dir=dir)

         dPar   = [dBx; dBy; dBz]         # Parallel direction
         dPerpI = cross([0; -1; 0], dPar) # Perpendicular direction in-plane
         dPerpO = cross(dPar, dPerpI)     # Perpendicular direction out-of-plane

         uPar   = transpose(particle[iB][1:3,:])*dPar   ./ cAlfven
         uPerpI = transpose(particle[iB][1:3,:])*dPerpI ./ cAlfven
         uPerpO = transpose(particle[iB][1:3,:])*dPerpO ./ cAlfven
      end

      if plottype[i] == 1
         h = ax[i].hist2d(uy, ux, bins=nbin,
            norm=matplotlib.colors.LogNorm(),density=true, range=binRange)
      elseif plottype[i] == 2
         h = ax[i].hist2d(ux, uz, bins=nbin,
            norm=matplotlib.colors.LogNorm(),density=true, range=binRange)
      elseif plottype[i] == 3
         h = ax[i].hist2d(uy, uz, bins=nbin,
            norm=matplotlib.colors.LogNorm(),density=true, range=binRange)
      elseif plottype[i] == 4
         h = ax[i].hist2d(uPerpO, uPerpI, bins=nbin,
            norm=matplotlib.colors.LogNorm(),density=true, range=binRange)
      elseif plottype[i] == 5
         h = ax[i].hist2d(uPerpI, uPar, bins=nbin,
            norm=matplotlib.colors.LogNorm(),density=true, range=binRange)
      elseif plottype[i] == 6
         h = ax[i].hist2d(uPerpO, uPar, bins=nbin,
            norm=matplotlib.colors.LogNorm(),density=true, range=binRange)
      else
         @error "Unknown PlotVType!"
      end
      ax[i].set_aspect("equal", "box")
      #ax[i].grid(true)
      ax[i].axvline(x=0.0, ls="--", linewidth=0.5, color="k")
      ax[i].axhline(y=0.0, ls="--", linewidth=0.5, color="k")
      (i % 4 != 0) && ax[i].axes.xaxis.set_ticklabels([])

      if plottype[i] == 1
         if i in (4,8)
            ax[i].annotate(L"u_y", xy=(0.50, -0.20),xycoords="axes fraction",
               fontsize=fs)
         else
            ax[i].annotate(L"u_y", xy=xlPos,xycoords="axes fraction",
               fontsize=fs)
         end
         ax[i].annotate(L"u_x", xy=ylPos,xycoords="axes fraction",
            fontsize=fs)
      elseif plottype[i] == 2
         ax[i].annotate(L"u_x", xy=xlPos,xycoords="axes fraction",
            fontsize=fs)
         ax[i].annotate(L"u_z", xy=ylPos,xycoords="axes fraction",
            fontsize=fs)
      elseif plottype[i] == 3
         if i in (4,8)
            ax[i].annotate(L"u_y", xy=(0.50, -0.20),xycoords="axes fraction",
               fontsize=fs)
         else
            ax[i].annotate(L"u_y", xy=xlPos,xycoords="axes fraction",
               fontsize=fs)
         end
         ax[i].annotate(L"u_z", xy=ylPos,xycoords="axes fraction",
            fontsize=fs)
      elseif plottype[i] == 4
         ax[i].annotate(L"u_{\perp Out}", xy=xlPos,xycoords="axes fraction",
            fontsize=fs)
         ax[i].annotate(L"u_{\perp In}", xy=ylPos,xycoords="axes fraction",
            fontsize=fs)
      elseif plottype[i] == 5
         ax[i].annotate(L"u_{\perp In}", xy=xlPos,xycoords="axes fraction",
            fontsize=fs)
         ax[i].annotate(L"u_\parallel", xy=ylPos,xycoords="axes fraction",
            fontsize=fs)
      elseif plottype[i] == 6
         ax[i].annotate(L"u_{\perp Out}", xy=xlPos,xycoords="axes fraction",
            fontsize=fs)
         ax[i].annotate(L"u_\parallel", xy=ylPos,xycoords="axes fraction",
            fontsize=fs)
      end

      if pType == 'e'
         ax[i].text(-0.05,1.02, @sprintf("x[%3.3f,%3.3f], z[%3.3f,%3.3f]",
            region[1,iB],region[2,iB],region[5,iB],region[6,iB]), FontSize=7,
            transform=ax[i].transAxes)
         ax[i].text(0.05,0.05,string(iB), FontSize=fs, color="r",
            transform=ax[i].transAxes)
      elseif pType == 'i'
         ax[i].text(-0.05,1.02, @sprintf("x[%3.3f,%3.3f], z[%3.3f,%3.3f]",
            region[1,iB],region[2,iB],region[5,iB],region[6,iB]),
            FontSize=7, transform=ax[i].transAxes)
         ax[i].text(0.83,0.05,string(iB+4), FontSize=fs, color="c",
            transform=ax[i].transAxes)
      end
   end

   im = plt.gca().get_children()[1]
   cax = fig.add_axes([0.1,0.93,0.8,0.02])
   colorbar(im, ax=ax, cax=cax, orientation="horizontal")
   cax.tick_params(which="both", axis="x", direction="out",color="k")

   return ax
end

function show_box_region()
   #dir = "/Users/hyzhou/Documents/Computer/Julia/BATSRUS/VisAnaJulia"
   dir = "/Users/hyzhou/Documents/Computer/DeepBlue"
   fnameE = "cut_particles0_region0_1_t00001640_n00020369.out"
   fnameI = "cut_particles1_region0_2_t00001640_n00020369.out"

   fnameField = "3d_var_region0_0_"*fnameE[end-22:end]

   fs = 10
   plotrange = [-2.0, -1.8, -0.35, 0.35]
   sequence = 129 # cut plane index starting from -

   data = readdata(fnameField; dir)

   X, Z, Bx = cutdata(data, "Bx"; dir="y", sequence, plotrange)
   X, Z, Bz = cutdata(data, "Bz"; dir="y", sequence, plotrange)
   X, Z, Ex = cutdata(data, "Ex"; dir="y", sequence, plotrange)

   fig, ax = plt.subplots(1,1,figsize=(8.0,4.0))
   #plt.set_cmap("RdBu_r")
   plt.set_cmap("seismic")

   region = Array{Float32,2}(undef,6,8)
   #=
   region[:,1] = [-1.950, -1.945, -0.08, 0.08, -0.10, -0.06]
   region[:,2] = [-1.915, -1.910, -0.08, 0.08, -0.10, -0.06]
   region[:,3] = [-1.905, -1.900, -0.08, 0.08, -0.10, -0.06]
   region[:,4] = [-1.895, -1.890, -0.08, 0.08, -0.10, -0.06]
   =#
   region[:,1] = [-1.930, -1.925, -0.08, 0.08, -0.30, -0.26]
   region[:,2] = [-1.930, -1.925, -0.08, 0.08, -0.20, -0.16]
   region[:,3] = [-1.930, -1.925, -0.08, 0.08, -0.10, -0.06]
   region[:,4] = [-1.930, -1.925, -0.08, 0.08,  0.00,  0.04]

   #=
   region[:,5] = [-1.865, -1.855, -0.08, 0.08, -0.29, -0.25]
   region[:,6] = [-1.910, -1.900, -0.08, 0.08, 0.11, 0.15]
   region[:,7] = [-1.970, -1.960, -0.08, 0.08, -0.10, -0.06]
   region[:,8] = [-1.850, -1.840, -0.08, 0.08, -0.10, -0.06]
   =#
   region[:,5] = [-1.930, -1.920, -0.08, 0.08, -0.30, -0.26]
   region[:,6] = [-1.930, -1.920, -0.08, 0.08, -0.20, -0.16]
   region[:,7] = [-1.930, -1.920, -0.08, 0.08, -0.10, -0.06]
   region[:,8] = [-1.930, -1.920, -0.08, 0.08,  0.00,  0.04]

   c = ax.contourf(Z,X,Ex./1e3, 100, norm=matplotlib.colors.DivergingNorm(0),
      vmin=-12e1, vmax=12e1)

   # create an axes on the right side of ax. The width of cax will be 5%
   # of ax and the padding between cax and ax will be fixed at 0.05 inch.
   divider = axes_grid1.make_axes_locatable(ax)
   cax = divider.append_axes("right", size="5%", pad=0.05)
   colorbar(c, cax=cax)
   ax.invert_yaxis()
   ax.set_aspect("equal", "box")

   ax.set_xlabel(L"z [R_G]", fontsize=14)
   ax.set_ylabel(L"x [R_G]", fontsize=14)
   ax.set_title(L"Ex [mV/m]", fontsize=14)

   x  = @view X[:,1]
   z  = @view Z[1,:]

   seeds = select_seeds(x[10:end-10],z[10:end-10]; nSeed=1)
   xstart, zstart = seeds[1,:], seeds[2,:]
   append!(xstart, [-1.87, -1.97, -1.93, -1.93, -1.93, -1.95, -1.82, -1.93, -1.93, -1.93])
   append!(zstart, [-0.1, -0.1, -0.15, 0.0, 0.2, -0.25, -0.1, -0.08, -0.095, 0.116])

   xl = [Vector{Float32}(undef,0) for _ in 1:length(xstart)]
   zl = [Vector{Float32}(undef,0) for _ in 1:length(xstart)]
   for i in eachindex(xstart)
      xs, zs = xstart[i], zstart[i]
      xl[i], zl[i] = trace2d_rk4(Bx, Bz, xs, zs, x, z, ds=0.02, maxstep=20000,
         gridType="ndgrid")
   end
   [ax.plot(zl[j],xl[j],"-",color="k",lw=1.0) for j in eachindex(xstart)]

   ax.contour(Z,X,Bz,[0.],colors="k",linestyles="dotted",linewidths=1.)

   for iB = 1:4
      rect = matplotlib.patches.Rectangle( (region[5,iB], region[1,iB]),
         region[6,iB]-region[5,iB], region[2,iB]-region[1,iB],
         lw=1.0, fill=true, alpha=0.9, fc="r") # facecolor="none"
      ax.add_patch(rect)
   end

   for iB = 5:8
      rect = matplotlib.patches.Rectangle( (region[5,iB], region[1,iB]),
         region[6,iB]-region[5,iB], region[2,iB]-region[1,iB],
         lw=1.0, fill=true, alpha=0.9, fc="c") # facecolor="none"
      ax.add_patch(rect)
   end

   ax.annotate("1", (region[5,1]-2e-2,region[2,1]), fontsize=10, color="r")
   ax.annotate("2", (region[5,2]-2e-2,region[2,2]), fontsize=10, color="r")
   ax.annotate("3", (region[5,3]-2e-2,region[2,3]+0.5e-2),fontsize=10,color="r")
   ax.annotate("4", (region[5,4]-2e-2,region[2,4]+1e-2), fontsize=10,color="r")

   for i = 5:8
      ax.annotate(string(i), (region[5,i]-2e-2, region[2,i]+0.2e-2), fontsize=10,
         color="c")
   end

   return ax
end

function HF_velocity()
   dir = "/Users/hyzhou/"
   fnameE = "cut_particles0_region0_1_t00001640_n00020369.out"
   fnameI = "cut_particles1_region0_2_t00001640_n00020369.out"

   fnameField = "3d_var_region0_0_"*fnameE[end-22:end]

   data = readdata(fnameField; dir)

   x = data.x[:,:,:,1]
   y = data.x[:,:,:,2]
   z = data.x[:,:,:,3]

   ux_ = findfirst(x->x=="uxs1", lowercase.(data.head.wnames))
   uy_ = findfirst(x->x=="uys1", lowercase.(data.head.wnames))
   uz_ = findfirst(x->x=="uzs1", lowercase.(data.head.wnames))

   Ux = @view data.w[:,:,:,ux_]
   Uy = @view data.w[:,:,:,uy_]
   Uz = @view data.w[:,:,:,uz_]

   #region[:,5] = [-1.850, -1.840, -0.08, 0.08, -0.29, -0.25]
   #region[:,6] = [-1.910, -1.900, -0.08, 0.08, 0.11, 0.15]
   limits = [-1.910, -1.900, -0.08, 0.08, 0.11, 0.15]
   xnew, ynew, znew, UxNew, UyNew, UzNew = subvolume(x,y,z, Ux,Uy,Uz, limits)

   # Average over the selected volume
   Ūx, Ūy, Ūz = mean(UxNew), mean(UyNew), mean(UzNew)

   return Ūx, Ūy, Ūz
end

dir = "/Users/hyzhou/Documents/Computer/DeepBlue"
fnameField = "3d_var_region0_0_t00001640_n00020369.out"
PType = 'e'
plottype = 2

if PType == 'e'
   fnameParticle = "cut_particles0_region0_1_t00001640_n00020369.out"
elseif PType == 'i'
   fnameParticle = "cut_particles1_region0_2_t00001640_n00020369.out"
end

# Define regions
xC, yC, zC = -1.75, 0.0, -0.2
xL, yL, zL = 0.008, 0.2, 0.03 # box length in x,y,z

@time region, particle = dist_select(
   fnameParticle, xC, yC, zC, xL, yL, zL,
   dir=dir, ParticleType=PType)

@time dist_scan(region, particle, PType, plottype; dir, fnameField)

@time plotExCut(fnameField, region, xC,yC,zC,xL,yL,zL; dir)

#ax = dist_plot('i')
#ux, uy, uz = HF_velocity()
#show_box_region() # horizontal


text:
# ### m8.3t

using Pkg, DrWatson

@quickactivate "StatisticalRethinkingTuring"
using Turing
using StatisticalRethinking
Turing.setprogress!(false)

# Turing model
@model ppl8_3(y) = begin
    α ~ Normal(1, 10)
    σ ~ truncated(Cauchy(0, 1), 0, Inf)

    y .~ Normal(α, σ)
end

y = [-1,1]

# Sample

m8_3t = ppl8_3(y)
nchains = 4; sampler = NUTS(0.65); nsamples=2000
chns8_3t = mapreduce(c -> sample(m8_3t, sampler, nsamples), chainscat, 1:nchains)

# Results rethinking

m83rethinking = "
      mean   sd  5.5% 94.5% n_eff Rhat
alpha 0.09 1.63 -2.13  2.39   959    1
sigma 2.04 2.05  0.68  4.83  1090    1
";

# End of `08/m8.3t.jl`


text:
#####
##### `AbstractLayout`
#####

abstract type AbstractLayout end

struct CPUDevice end

struct Layout{L, S, D} <: AbstractLayout
    length::L
    size::S
    device::D
    ismutable::Bool
end

Base.length(layout::Layout) = layout.length

Base.size(layout::Layout) = layout.size

device(layout::Layout) = layout.device

ismutable(layout::Layout) = layout.ismutable

#####
##### `AbstractDomain`
#####

abstract type AbstractDomain end

struct RealDomain <: AbstractDomain end

struct ComplexDomain <: AbstractDomain end

#####
##### `AbstractArgument`
#####

abstract type AbstractArgument end

#== `Signature` ==#

struct Signature{I <: Tuple{Vararg{AbstractArgument}},
                 O <: Tuple{Vararg{AbstractArgument}}}
    input::I
    output::O
end

Signature(input, output) = Signature(markupify.(tuplify(input)),
                                     markupify.(tuplify(output)))

tuplify(x) = tuple(x)
tuplify(x::Tuple) = x

markupify(x::AbstractArgument) = x
markupify(x) = markup(x)

#== `Ignore` ==#

struct Ignore <: AbstractArgument end

Base.length(::Ignore) = 0

#== `AbstractVariable` ==#

abstract type AbstractVariable <: AbstractArgument end

struct Scalar{D <: AbstractDomain} <: AbstractVariable
    domain::D
end

const RealScalar = Scalar{RealDomain}

RealScalar() = Scalar(RealDomain())

const ComplexScalar = Scalar{ComplexDomain}

ComplexScalar() = Scalar(ComplexDomain())

struct Tensor{D <: AbstractDomain, L <: AbstractLayout} <: AbstractVariable
    domain::D
    layout::L
end

const RealTensor = Tensor{RealDomain}

RealTensor(layout) = Tensor(RealDomain(), layout)

const ComplexTensor = Tensor{ComplexDomain}

ComplexTensor(layout) = Tensor(ComplexDomain(), layout)

Base.length(::Scalar) = 1
Base.length(t::Tensor) = length(t.layout)

#####
##### `@sig`
#####

macro sig(expr)
    signature_type_from_expr(expr)
end

const MALFORMED_SIG_ERROR_MESSAGE = "Malformed expression given to `@sig`; see `@sig` docstring for proper format."

function signature_type_from_expr(expr)
    @assert(Meta.isexpr(expr, :call) && expr.args[1] === :→ && length(expr.args) === 3, MALFORMED_SIG_ERROR_MESSAGE)
    input_types = map(parse_into_markup_type, split_infix_args(expr.args[2], :⊕))
    output_types = map(parse_into_markup_type, split_infix_args(expr.args[3], :⊕))
    return :(Signature{<:Tuple{$(input_types...)}, <:Tuple{$(output_types...)}})
end

split_infix_args(invocation::Symbol, ::Symbol) = (invocation,)

function split_infix_args(invocation::Expr, op::Symbol)
    if Meta.isexpr(invocation, :call) && invocation.args[1] === op
        return (split_infix_args(invocation.args[2], op)..., invocation.args[3])
    end
    return (invocation,)
end

function parse_into_markup_type(x)
    if x === :R
        return :(Scalar{RealDomain})
    elseif x === :C
        return :(Scalar{ComplexDomain})
    elseif x === :_
        return :(Ignore)
    elseif isa(x, Expr) && length(x.args) === 1
        if Meta.isexpr(x, :vect)
            domain = x.args[1]
            if domain === :R
                return :(Tensor{RealDomain})
            elseif domain === :C
                return :(Tensor{ComplexDomain})
            end
        elseif Meta.isexpr(x, :braces)
            vararg_type = parse_into_markup_type(x.args[1])
            return :(Vararg{$vararg_type})
        end
    end
    error(string("Encountered unparseable signature element `", x, "`. ",
                 MALFORMED_SIG_ERROR_MESSAGE))
end


text:
export GitConfig, lookup, set!, global_config

type GitConfig
    ptr::Ptr{Void}

    function GitConfig(ptr::Ptr{Void})
        @assert ptr != C_NULL
        cfg = new(ptr)
        finalizer(cfg, free!)
        return cfg
    end
end

GitConfig(path::AbstractString) = begin
    cfg_ptr = Ptr{Void}[0]
    @check ccall((:git_config_open_ondisk, libgit2), Cint,
                 (Ptr{Ptr{Void}}, Ptr{UInt8}), cfg_ptr, path)
    return GitConfig(cfg_ptr[1])
end

GitConfig(r::GitRepo) = begin
    cfg_ptr = Ptr{Void}[0]
    @check ccall((:git_repository_config, libgit2), Cint,
                 (Ptr{Ptr{Void}}, Ptr{Void}), cfg_ptr, r)
    return GitConfig(cfg_ptr[1])
end

free!(cfg::GitConfig) = begin
    if cfg.ptr != C_NULL
        ccall((:git_config_free, libgit2), Void, (Ptr{Void},), cfg.ptr)
        cfg.ptr == C_NULL
    end
end

Base.convert(::Type{Ptr{Void}}, cfg::GitConfig) = cfg.ptr

typealias GitConfigType Union(Bool,Int32,Int64,AbstractString)

Base.getindex(c::GitConfig, key::AbstractString) = lookup(AbstractString, c, key)
Base.setindex!{T<:GitConfigType}(c::GitConfig, v::T, key::AbstractString) = set!(T, c, key, v)

Base.delete!(c::GitConfig, key::AbstractString) = begin
    err = ccall((:git_config_delete_entry, libgit2), Cint, (Ptr{Void}, Ptr{UInt8}), c, key)
    return err == GitErrorConst.ENOTFOUND ? false : true
end

function cb_each_key(entry_ptr::Ptr{ConfigEntryStruct}, o::Ptr{Void})
    entry = unsafe_load(entry_ptr)::ConfigEntryStruct
    n = bytestring(entry.name)
    produce(n)
    return GitErrorConst.GIT_OK
end

function cb_each_val(entry_ptr::Ptr{ConfigEntryStruct}, o::Ptr{Void})
    entry = unsafe_load(entry_ptr)
    v = bytestring(entry.value)
    produce(v)
    return GitErrorConst.GIT_OK
end

function cb_each_pair(entry_ptr::Ptr{ConfigEntryStruct}, o::Ptr{Void})
    entry = unsafe_load(entry_ptr)
    n = bytestring(entry.name)
    v = bytestring(entry.value)
    produce((n, v))
    return GitErrorConst.GIT_OK
end

const c_cb_each_key    = cfunction(cb_each_key,  Cint, (Ptr{ConfigEntryStruct}, Ptr{Void}))
const c_cb_each_val    = cfunction(cb_each_val,  Cint, (Ptr{ConfigEntryStruct}, Ptr{Void}))
const c_cb_each_kvpair = cfunction(cb_each_pair, Cint, (Ptr{ConfigEntryStruct}, Ptr{Void}))

Base.keys(c::GitConfig) = begin
    @task ccall((:git_config_foreach, libgit2), Cint,
                (Ptr{Void}, Ptr{Void}, Ptr{Void}), c, c_cb_each_key, C_NULL)
end

Base.values(c::GitConfig) = begin
    @task ccall((:git_config_foreach, libgit2), Cint,
                (Ptr{Void}, Ptr{Void}, Ptr{Void}), c, c_cb_each_val, C_NULL)
end

Base.start(c::GitConfig) = begin
    t = @task ccall((:git_config_foreach, libgit2), Cint,
                    (Ptr{Void}, Ptr{Void}, Ptr{Void}), c, c_cb_each_kvpair, C_NULL)
    (consume(t), t)
end

Base.done(c::GitConfig, state) = istaskdone(state[2])
Base.next(c::GitConfig, state) = begin
    v = consume(last(state))
    (state[1], (v, state[2]))
end

Base.Dict(c::GitConfig) = Dict{ByteString, ByteString}(c)

function global_config()
    cfg_ptr = Ptr{Void}[0]
    @check ccall((:git_config_open_default, libgit2), Cint, (Ptr{Ptr{Void}},), cfg_ptr)
    return GitConfig(cfg_ptr[1])
end

function lookup(::Type{Bool}, c::GitConfig, name::AbstractString)
    out = Cint[0]
    @check ccall((:git_config_get_bool, libgit2), Cint,
                 (Ptr{Cint}, Ptr{Void}, Ptr{UInt8}), out, c, name)
    if err == GitErrorConst.GIT_OK
        return bool(out[1])
    elseif err == GitErrorConst.ENOTFOUND
        return nothing
    else
        throw(LibGitError(err))
    end
end

function set!(::Type{Bool}, c::GitConfig, name::AbstractString, value::Bool)
    @check ccall((:git_config_set_bool, libgit2), Cint,
                 (Ptr{Void}, Ptr{UInt8}, Cint), c, name, value)
    return c
end

function lookup(::Type{Int32}, c::GitConfig, name::AbstractString)
    out = Int32[0]
    err = ccall((:git_config_get_int32, libgit2), Cint,
                (Ptr{Int32}, Ptr{Void}, Ptr{UInt8}), out, c, name)
    if err == GitErrorConst.GIT_OK
        return out[1]
    elseif err == ENOTFOUND
        return nothing
    else
        throw(LibGitError(err))
    end
end

function set!(::Type{Int32}, c::GitConfig, name::AbstractString, value::Int32)
    @check ccall((:git_config_set_int32, libgit2), Cint,
                 (Ptr{Void}, Ptr{UInt8}, Int32), c, name, value)
    return c
end

function lookup(::Type{Int64}, c::GitConfig, name::AbstractString)
    out = Int64[0]
    err = ccall((:git_config_get_int64, libgit2), Cint,
                (Ptr{Int64}, Ptr{Void}, Ptr{UInt8}), out, c, name)
    if err == GitErrorConst.GIT_OK
        return out[1]
    elseif err == GitErrorConst.ENOTFOUND
        return nothing
    else
        throw(GitError(err))
    end
end

function set!(::Type{Int64}, c::GitConfig, name::AbstractString, value::Int64)
    @check ccall((:git_config_set_int64, libgit2), Cint,
                 (Ptr{Void}, Ptr{UInt8}, Int64), c, name, value)
    return c
end

function lookup{T<:AbstractString}(::Type{T}, c::GitConfig, name::AbstractString)
    out = Ptr{UInt8}[0]
    err = ccall((:git_config_get_string, libgit2), Cint,
                (Ptr{Ptr{UInt8}}, Ptr{Void}, Ptr{UInt8}), out, c, name)
    if err == GitErrorConst.GIT_OK
        return bytestring(out[1])
    elseif err == GitErrorConst.ENOTFOUND
        return nothing
    else
        throw(GitError(err))
    end
end

function set!{T<:AbstractString}(::Type{T}, c::GitConfig, name::AbstractString, value::AbstractString)
    @check ccall((:git_config_set_string, libgit2), Cint,
                 (Ptr{Void}, Ptr{UInt8}, Ptr{UInt8}), c, name, value)
    return c
end


text:
## fixme ** should work
# or add_dirs in project
VW_DIR="/home/graham_s/VirtualWorlds/projects/"
push!(LOAD_PATH,joinpath(VW_DIR,"TBComponents.jl","src"))
push!(LOAD_PATH,joinpath(VW_DIR,"DDIMeta.jl","src"))

STB_DIR=pwd() # joinpath(VW_DIR,"stb.jl")
for dr in ["core","web","data_mapping","persist", "general"]
    push!(LOAD_PATH, joinpath(STB_DIR,"src",dr))
end
for dr in ["test","scratch","scripts"]
    push!(LOAD_PATH,joinpath(STB_DIR,dr))
end

# using Revise


text:
module TestExamples

using Test
using Transducers

using ..Utils: foreach_executor

module OnlineExtrema
include("../../../examples/online_extrema.jl")
end

module OnlineHistogram
include("../../../examples/online_histogram.jl")
end

function test_extrema(ex)
    @test OnlineExtrema.online_extrema(N = 100, ex = ex) isa Any
end

function test_extrema()
    foreach_executor(test_extrema)
end

function test_histogram(ex)
    @test OnlineHistogram.online_histogram(N = 100, ex = ex) isa Any
end

function test_histogram()
    foreach_executor(test_histogram)
end

end  # module


text:
"""
	GLO{T1<:Real, T2<:Integer} a.k.a. GPLinearODE

A structure that holds all of the relevant information for constructing the
GLOM model

# Arguments
- `kernel::Function`: kernel function for the latent Gaussian process (GP)
- `n_kern_hyper::T2`: amount of hyperparameters for the kernel function
- `n_dif::T2`: amount of times you are differenting the base kernel
- `n_out::T2`: amount of outputs you are jointly modelling
- `x_obs::Vector{T1}`: the observation times/phases
- `y_obs::Vector{T1}`: the flattened, observed data
- `noise::Vector{T1}`: the measurement noise at all observations (σ, not σ^2)
- `normals::Vector{T1}`: the normalization of each section of y_obs
- `a::Matrix{T1}`: the meta kernel coefficients
- `non_zero_hyper_inds::Vector{T2}`: the indices of the non-zero hyperparameters
- `coeff_orders::AbstractArray{T2,6}`: The powers that each `a` coefficient is taken to for each part of the covariance matrix construction. Used for constructing differentiated versions of the kernel
- `coeff_coeffs::AbstractArray{T2,4}`: The coefficients that each `a` coefficient is multiplied by for each part of the covariance matrix construction. Used for constructing differentiated versions of the kernel
- `covariance::Array{T1, 3}`: the measurement covariance at all observations
- `has_covariance::Bool`: Whether or not there are measurement covariances (as opposed to just having measurement noises)
- `kernel_changes_with_output::Bool`: A flag for if the kernel is not the same for each output (e.g. there is a time delay kernel hyperparameter that changes)
"""
struct GLO{T1<:Real, T2<:Integer}
	kernel::Function
	n_kern_hyper::T2
	n_dif::T2
	n_out::T2
	x_obs::Vector{T1}
	y_obs::Vector{T1}
	noise::Vector{T1}
	normals::Vector{T1}
	a::Matrix{T1}
	non_zero_hyper_inds::Vector{T2}
	coeff_orders::AbstractArray{T2,6}
	coeff_coeffs::AbstractArray{T2,4}
	covariance::Array{T1, 3}
	has_covariance::Bool
	kernel_changes_with_output::Bool

	function GLO(
		kernel::Function,
		n_kern_hyper::T2,
		n_dif::T2,
		n_out::T2,
		x_obs::Vector{T1},
		y_obs::Vector{T1};
		noise::Vector{T1} = zeros(length(y_obs)),
		normals::Vector{T1} = ones(n_out),
		a::Matrix{T1} = ones(n_out, n_dif),
		covariance::Array{T1, 3} = zeros(length(x_obs), n_out, n_out),
		kernel_changes_with_output::Bool=false
		) where {T1<:Real, T2<:Integer}

		@assert size(a) == (n_out, n_dif)
		non_zero_hyper_inds = append!(findall(!iszero, collect(Iterators.flatten(a))), collect(1:n_kern_hyper) .+ length(a))
		coeff_orders, coeff_coeffs = coefficient_orders(n_out, n_dif, a=a)
		has_covariance = (covariance != zeros(length(x_obs), n_out, n_out))

		return GLO(kernel, n_kern_hyper, n_dif, n_out, x_obs, y_obs, noise, normals, a, non_zero_hyper_inds, coeff_orders, coeff_coeffs, covariance, has_covariance, kernel_changes_with_output)
	end
	function GLO(
		kernel::Function,  # kernel function
		n_kern_hyper::T2,  # amount of hyperparameters for the kernel function
		n_dif::T2,  # amount of times you are differenting the base kernel
		n_out::T2,  # amount of scores you are jointly modelling
		x_obs::Vector{T1}, # the observation times/phases
		y_obs::Vector{T1},  # the flattened, observed data
		noise::Vector{T1},  # the measurement noise at all observations
		normals::Vector{T1},  # the normalization of each section of y_obs
		a::Matrix{T1},
		non_zero_hyper_inds::Vector{T2},
		coeff_orders::Array{T2,6},
		coeff_coeffs::Array{T2,4},
		covariance::Array{T1, 3},
		has_covariance::Bool,
		kernel_changes_with_output::Bool) where {T1<:Real, T2<:Integer}

		@assert isfinite(kernel(ones(n_kern_hyper), randn(), zeros(Int64, 2 + n_kern_hyper)))  # make sure the kernel is valid by testing a sample input
		@assert 0 < n_dif <= 3
		@assert 0 < n_out
		n_meas = length(x_obs)
		@assert (n_meas * n_out) == length(y_obs) == length(noise)
		@assert n_meas == size(covariance, 1)
		@assert length(normals) == n_out
		@assert size(a) == (n_out, n_dif)
		@assert size(coeff_orders) == (n_out, n_out, n_dif, n_dif, n_out, n_dif)  # maybe unnecessary due to the fact that we construct it
		@assert size(coeff_coeffs) == (n_out, n_out, n_dif, n_dif)  # maybe unnecessary due to the fact that we construct it
		@assert length(non_zero_hyper_inds) == length(findall(!iszero, collect(Iterators.flatten(a)))) + n_kern_hyper
		@assert n_out == size(covariance, 2) == size(covariance, 3)
		@assert (covariance != zeros(n_meas, n_out, n_out)) == has_covariance

		return new{typeof(x_obs[1]),typeof(n_kern_hyper)}(kernel, n_kern_hyper, n_dif, n_out, x_obs, y_obs, noise, normals, a, non_zero_hyper_inds, coeff_orders, coeff_coeffs, covariance, has_covariance, kernel_changes_with_output)
	end
end

"""
	normalize_GLO!(glo)

Standardize the outputs in `glo.y_obs` (and `glo.noise`) to have 0-mean and unit
variance and store the normalization in `glo.normals`
"""
function normalize_GLO!(glo::GLO)
	renorms = ones(glo.n_out)
	for i in 1:glo.n_out
		inds = (i:glo.n_out:length(glo.y_obs))
		glo.y_obs[inds] .-= mean(glo.y_obs[inds])
		renorms[i] = std(glo.y_obs[inds])
	end
	normalize_GLO!(glo, renorms)
end

"""
	normalize_GLO!(glo, renorms)

Divide the outputs in `glo.y_obs` and `glo.noise` and multiply `glo.normals` by
renorms
"""
function normalize_GLO!(glo::GLO, renorms::Vector)
	@assert length(renorms) == glo.n_out
	for i in 1:glo.n_out
		inds = (i:glo.n_out:length(glo.y_obs))
		glo.normals[i] *= renorms[i]
		glo.y_obs[inds] /= renorms[i]
		glo.noise[inds] /= renorms[i]
	end
	if glo.has_covariance
		renorm_mat = renorms .* transpose(renorms)
		for i in 1:length(glo.x_obs)
			glo.covariance[i, :, :] ./= renorm_mat
		end
	end
end

max_stars_repo_name: christiangil/GPLinearODEMaker

max_issues_repo_name: christiangil/GPLinearODEMaker

max_forks_repo_name: christiangil/GPLinearODEMaker


text:
include("smear_FD.jl")
include("calc_Focc.jl")
include("calc_entropy.jl")
include("sum_upto_E_fermi.jl")

function test_main(kT::Float64)
    Nstates = 8
    Nelectrons = 6
    evals = Array{Float64}(Nstates)
    evals = [-2.4, -1.0, -0.5, -0.2, -0.19, -0.10, -0.05, -0.05]
    println("\nkT = ", kT)
    spinpol = false
    println("\nspinpol = ", spinpol)
    Focc, E_fermi = calc_Focc(evals, Nelectrons, kT)
    for ist = 1:Nstates
        @printf("%18.10f %18.10f\n", evals[ist], Focc[ist])
    end
    @printf("E_fermi = %18.10f\n", E_fermi)
    integFocc = sum_upto_E_fermi( Focc, evals, E_fermi )
    @printf("integFocc = %18.10f\n", integFocc)
    @printf("sum(Focc) = %18.10f\n", sum(Focc))
    @printf("Entropy (-TS) = %18.10f\n", calc_entropy(Focc, kT, is_spinpol=spinpol))
end

test_main(0.01)


text:
function clear_comments(data::String)
    res = replace(data, r"(%.*\n)", "")
    res = replace(res, r"(comment [^\n]*\n)", "")

    return res
end

mutable struct Parser
    data::String
    token
    token_type
    _next_token
    hashtable::Dict
    mode
    records::Dict
    line::Int64

    # regular expressions
    white::Regex
    nl::Regex
    token_re::Regex
end

function Parser(data::String)
    # compile some regexes
    white = r"[\n|\s]+"
    nl = r"[\n]"
    token_re = r"([^\s\"#%'(){}@,=]+|\n|@|\"|{|}|=|,)"

    Parser(data,
           1,
           white,
           nl,
           token_re)
end

function tokenize(bib::Parser)
end

function parser(bib::Parser)
    for item in matchall(token_re, data)
    for item in self.token_re.finditer(self.data):
        i = item.group(0)
        if self.white.match(i) :
            if self.nl.match(i) :
                self.line += 1
            continue
        else :
            yield i
end


text:
using Documenter, SolarLEDs

makedocs(
    modules = [SolarLEDs],
    format = Documenter.HTML(),
    checkdocs = :exports,
    sitename = "SolarLEDs.jl",
    pages = Any["index.md"]
)

deploydocs(
    repo = "github.com/yakir12/SolarLEDs.jl.git",
)


text:
__precompile__()

module Screen
	include("./Types.jl")
	include("./GPU.jl")

	using PyCall
	using OpenCL
	using .Types
	using .GPU

	struct ScreenConfig
		sct::PyObject
		checkHeight::Int
		monitor::Int
		sectors::Array{Int}
		gpuContext::GPUContext
		gpuOpts::GPUOpts
		sectorCount::Int
		sectorSize::Array{Int}
	end

	function config_screen(monitor, height, sectors)::ScreenConfig
		sct = pyimport("mss").mss()
		monitorSize = sct.monitors[monitor]

		gpuOpts = GPUOpts(
			Int32(monitorSize["width"]), # Screen width
			Int32(monitorSize["height"]), # Screen height
			Int32(height), # Ignore min vertical
			Int32(monitorSize["height"] - height), # Ignore max vertical
			Int32(height), # Ignore min horizontal
			Int32(monitorSize["width"] - height), # Ignore max horizontal
			Int32(floor(monitorSize["width"] / sectors[1])), # Sector width top
			Int32(floor((monitorSize["height"] - 2 * height) / sectors[2])), # Sector height right
			Int32(floor(monitorSize["width"] / sectors[3])), # Sector width bottom
			Int32(floor((monitorSize["height"] - 2 * height) / sectors[4])), # Sector height left
			0, # Sector offset top
			Int32(sectors[1]), # Sector offset right
			Int32(sum(sectors[1:2])), # Sector offset bottom
			Int32(sum(sectors[1:3])), # Sector offset left
			Int32(sectors[3]), # Sector bottom count
			Int32(sectors[4]), # Sector left count
			Int32(sum(sectors))
		)

		sectorSize = [
			gpuOpts.sector_width_top * height,
			gpuOpts.sector_width_right * height,
			gpuOpts.sector_width_bottom * height,
			gpuOpts.sector_width_left * height
		]

		gpuContext = init_kernel()

		return ScreenConfig(
			sct,
			height,
			monitor,
			sectors,
			gpuContext,
			gpuOpts,
			sum(sectors),
			sectorSize
		)
	end

	function get_sector(config::ScreenConfig, index::Int)
		if index > config.gpuOpts.sector_offset_left
			return 4
		elseif index > config.gpuOpts.sector_offset_bottom
			return 3
		elseif index > config.gpuOpts.sector_offset_right
			return 2
		else
			return 1
		end
	end

	function hsv_to_rgb(h::Float64, s::Float64, v::Float64)::RGB
		c = v * s
		x = c * (1 - abs((h % 2) - 1))
		m = v - c
		col = [0, 0, 0]

		if h < 1
			col = (c, x, 0)
		elseif h < 2
			col = (x, c, 0)
		elseif h < 3
			col = (0, c, x)
		elseif h < 4
			col = (0, x, c)
		elseif h < 5
			col = (x, 0, c)
		else
			col = (c, 0, x)
		end

		return UInt8.(round.((col .+ m) .* 255))
	end

	function parse_data(config::ScreenConfig, raw::Array{Int32})
		avgColors = Array{RGB}(undef, config.sectorCount)

		@Threads.threads for i in 1:config.sectorCount
			offset = i * 4 - 3
			satWeight = raw[offset + 3]

			if satWeight == 0
				avgColors[i] = (0,0,0)
				continue
			end

			r = @fastmath(raw[offset] / satWeight / 255)
			g = @fastmath(raw[offset + 1] / satWeight / 255)
			b = @fastmath(raw[offset + 2] / satWeight / 255)

			cMin = min(r,g,b)
			cMax = max(r,g,b)
			delta = cMax - cMin
			sat = @fastmath((delta / cMax) ^ (1 / 2.5))

			hue = 0

			if cMax == r
				hue = @fastmath(((g - b) / delta) % 6)
			elseif cMax == g
				hue = @fastmath((b - r) / delta + 2)
			elseif cMax == b
				hue = @fastmath((r - g) / delta + 4)
			end

			avgColors[i] = hsv_to_rgb(abs(hue), sat, cMax)
		end

		return avgColors
	end

	function process_screen_gpu(config::ScreenConfig)
		screenshot = config.sct.grab(config.sct.monitors[config.monitor])
		rawData = process_data(config.gpuContext, config.gpuOpts, screenshot.raw)
		return parse_data(config, rawData)
	end

	export ScreenConfig, config_screen, process_screen_gpu
end


text:
####################
#
# Plot Shapefile
#
####################

# Shapefile Abstract
AbstractGeom = Shapefile.GeoInterface.AbstractGeometry

# Plot Shapefile - given array of shapes and a canvas
function plotshape(shparray::AbstractArray{T, 1},
                   canvas::Compose.Context;
                   convertcoords=lonlat_to_webmercator, img_width=12cm,
                   line_width=0.05mm, line_color="black", fill_color=nothing,
                   radius=[0.5cm]) where {T<:AbstractGeom}

    # Plot and return canvas
    draw_shp(shparray, canvas, convertcoords, line_width, line_color, fill_color, radius)

end

# Plot Shapefile - given array of shapes and an MBR
function plotshape(shparray::AbstractArray{T, 1},
                   MBR::Shapefile.Rect;
                   convertcoords=lonlat_to_webmercator, img_width=12cm,
                   options...) where {T<:AbstractGeom}

    # Create canvas
    canvas = create_canvas(MBR, convertcoords, img_width)

    # Plot
    plotshape(shparray, canvas,
              convertcoords=convertcoords, img_width=img_width; options...)

end

# Plot Shapefile - given an array of shapes
function plotshape(shparray::AbstractArray{T, 1};
                   options...) where {T<:AbstractGeom}

    # Create MBR
    MBR = minBR(shparray)

    # Plot
    plotshape(shparray, MBR; options...)
end

# Plot Shapefile - given the Shapefile Handle
function plotshape(shp::Shapefile.Handle; options...)
    plotshape(shp.shapes, shp.MBR; options...)
end

# Plot Shapefile - given the Shapefile Handle and an MBR
function plotshape(shp::Shapefile.Handle, MBR::Shapefile.Rect; options...)
    plotshape(shp.shapes, MBR; options...)
end

# Plot Shapefile - given the Shapefile Handle and a canvas
function plotshape(shp::Shapefile.Handle, canvas::Compose.Context; options...)
    plotshape(shp.shapes, canvas; options...)
end


text:
using VortexDistributions, Test

# @testset "Single vortex accuracy" begin

# simple test field
Nx = 400; Ny = 400
Lx = 200; Ly = Lx
x = LinRange(-Lx/2,Lx/2, Nx+1)[1:end-1]; y = LinRange(-Ly/2,Ly/2, Ny+1)[1:end-1];
psi0 = one.(x*y') |> complex

psi1 = Torus(copy(psi0),x,y)
psi2 = Torus(copy(psi0),x,y)
xp = 1.133
yp = 2.787
vp = PointVortex(xp,yp,1)
sp = ScalarVortex(vp)

vortex!(psi1,sp)
vfound1 = findvortices(psi1)
@test isapprox(vfound1[1].xv,xp,atol=0.01)
@test isapprox(vfound1[1].yv,yp,atol=0.01)

vortex!(psi2, vp)
vfound2 = findvortices(psi2)
@test isapprox(vfound2[1].xv,xp,atol=0.01)
@test isapprox(vfound2[1].yv,yp,atol=0.01)

@test (vfound1[1].xv === vfound2[1].xv && vfound1[1].yv === vfound2[1].yv && vfound1[1].qv === vfound2[1].qv)
# end


text:
function _bls_by_shape(bls::AbstractBL{SIDE},c::Circle) where SIDE
	if SIDE == 1
		stop = c.radius
		start = stop + bls.depth
	else
		throw("circle only has side=1, evaluated at $side")
	end
	return bls(start=start,stop=stop)
end

# Square is just special rectangle
_bls_by_shape(bls::AbstractBL,s::Square) = _bls_by_shape(bls,Rectangle(s))

function _bls_by_shape(bls::AbstractBL{SIDE},shape::Rectangle) where SIDE
	if SIDE==1
		stop = shape.origin.x - shape.a/2
		start = stop + bls.depth
	elseif SIDE==2
		stop = shape.origin.x + shape.a/2
		start = stop - bls.depth
	elseif SIDE==3
		stop = shape.origin.y - shape.b/2
		start = stop + bls.depth
	elseif SIDE==4
		stop = shape.origin.y + shape.b/2
		start = stop - bls.depth
	end
	return bls(start=start, stop=stop)
end

################################################################################
# Plotting

@recipe function f(bnd::Boundary{2})
	bnd.shape
end