text:
module BrickMesh

export brickmesh, centroidtocode, connectmesh, partition, mappings, commmapping

using MPI

"""
    linearpartition(n, p, np)

Partition the range `1:n` into `np` pieces and return the `p`th piece as a
range.

This will provide an equal partition when `n` is divisible by `np` and
otherwise the ranges will have lengths of either `floor(Int, n/np)` or
`ceil(Int, n/np)`.
"""
linearpartition(n, p, np) = range(div((p-1)*n, np) + 1, stop=div(p*n, np))

"""
    hilbertcode(Y::AbstractArray{T}; bits=8sizeof(T)) where T

Given an array of axes coordinates `Y` stored as `bits`-bit integers
the function returns the Hilbert integer `H`.

The encoding of the Hilbert integer is best described by example.
If 5-bits are used from each of 3 coordinates then the function performs

     X[2]|                       H[0] = A B C D E
         | /X[1]       ------->  H[1] = F G H I J
    axes |/                      H[2] = K L M N O
         0------ X[0]                   high low

where the 15-bit Hilbert integer = `A B C D E F G H I J K L M N O` is stored
in `H`

This function is based on public domain code from John Skilling which can be
found in <https://dx.doi.org/10.1063/1.1751381>.
"""
function hilbertcode(Y::AbstractArray{T}; bits=8sizeof(T)) where T
  # Below is Skilling's AxestoTranspose
  X = deepcopy(Y)
  n = length(X)
  M = one(T) << (bits-1)

  Q = M
  for j = 1:bits-1
    P = Q - one(T)
    for i = 1:n
      if X[i] & Q != zero(T)
        X[1] ⊻= P
      else
        t = (X[1] ⊻ X[i]) & P
        X[1] ⊻= t
        X[i] ⊻= t
      end
    end
    Q >>>= one(T)
  end

  for i = 2:n
    X[i] ⊻= X[i - 1]
  end

  t = zero(T)
  Q = M
  for j = 1:bits-1
    if X[n] & Q != zero(T)
      t ⊻= Q - one(T)
    end
    Q >>>= one(T)
  end

  for i = 1:n
    X[i] ⊻= t
  end

  # Below we transpose X and store it in H, i.e.:
  #
  #   X[0] = A D G J M               H[0] = A B C D E
  #   X[1] = B E H K N   <------->   H[1] = F G H I J
  #   X[2] = C F I L O               H[2] = K L M N O
  #
  # The 15-bit Hilbert integer is then = A B C D E F G H I J K L M N O
  H = zero(X)
  for i = 0:n-1, j = 0:bits-1
    k = i * bits + j
    bit = (X[n - mod(k,n)] >>> div(k,n)) & one(T)
    H[n - i] |= (bit << j)
  end

  return H
end

"""
    centroidtocode(comm::MPI.Comm, elemtocorner; coortocode, CT)

Returns a code for each element based on its centroid.

These element codes can be used to determine a linear ordering for the
partition function.

The communicator `comm` is used to calculate the bounding box for representing
the centroids in coordinates of type `CT`, defaulting to `CT=UInt64`.  These
integer coordinates are converted to a code using the function `coortocode`,
which defaults to `hilbertcode`.

The array containing the element corner coordinates, `elemtocorner`, is used
to compute the centroids.  `elemtocorner` is a dimension by number of corners
by number of elements array.
"""
function centroidtocode(comm::MPI.Comm, elemtocorner; coortocode=hilbertcode,
                        CT=UInt64)
  (d, nvert, nelem) = size(elemtocorner)

  centroids = sum(elemtocorner, dims=2) ./ nvert
  T = eltype(centroids)

  centroidmin = (nelem > 0) ? minimum(centroids, dims=3) : fill(typemax(T),d)
  centroidmax = (nelem > 0) ? maximum(centroids, dims=3) : fill(typemin(T),d)

  centroidmin = MPI.Allreduce(centroidmin, MPI.MIN, comm)
  centroidmax = MPI.Allreduce(centroidmax, MPI.MAX, comm)
  centroidsize = centroidmax - centroidmin

  # Fix centroidsize to be nonzero.  It can be zero for a couple of reasons.
  # For example, it will be zero if we have just one element.
  if iszero(centroidsize)
    centroidsize = ones(T, d)
  else
    for i = 1:d
      if iszero(centroidsize[i])
        centroidsize[i] = maximum(centroidsize)
      end
    end
  end

  code = Array{CT}(undef, d, nelem)
  for e = 1:nelem
    c = (centroids[:,1,e] .- centroidmin) ./ centroidsize
    X = CT.(floor.(typemax(CT).*BigFloat.(c, 16sizeof(CT))))
    code[:,e] = coortocode(X)
  end

  code
end

"""
    brickmesh(x, periodic; part=1, numparts=1; boundary)

Generate a brick mesh with coordinates given by the tuple `x` and the
periodic dimensions given by the `periodic` tuple.

The brick can optionally be partitioned into `numparts` and this returns
partition `part`.  This is a simple Cartesian partition and further
partitioning (e.g, based on a space-filling curve) should be done before the
mesh is used for computation.

By default boundary faces will be marked with a one and other faces with a
zero.  Specific boundary numbers can also be passed for each face of the brick
in `boundary`.  This will mark the nonperiodic brick faces with the given
boundary number.

# Examples

We can build a 3 by 2 element two-dimensional mesh that is periodic in the
\$x_2\$-direction with
```jldoctest brickmesh
julia> (elemtovert, elemtocoord, elemtobndy, faceconnections) =
        brickmesh((2:5,4:6), (false,true); boundary=((1,2), (3,4)));
```
This returns the mesh structure for

             x_2

              ^
              |
             6-  9----10----11----12
              |  |     |     |     |
              |  |  4  |  5  |  6  |
              |  |     |     |     |
             5-  5-----6-----7-----8
              |  |     |     |     |
              |  |  1  |  2  |  3  |
              |  |     |     |     |
             4-  1-----2-----3-----4
              |
              +--|-----|-----|-----|--> x_1
                 2     3     4     5

The (number of corners by number of elements) array `elemtovert` gives the
global vertex number for the corners of each element.
```jldoctest brickmesh
julia> elemtovert
4×6 Array{Int64,2}:
 1  2  3   5   6   7
 2  3  4   6   7   8
 5  6  7   9  10  11
 6  7  8  10  11  12
```
Note that the vertices are listed in Cartesian order.

The (dimension by number of corners by number of elements) array `elemtocoord`
gives the coordinates of the corners of each element.
```jldoctes brickmesh
julia> elemtocoord
2×4×6 Array{Int64,3}:
[:, :, 1] =
 2  3  2  3
 4  4  5  5

[:, :, 2] =
 3  4  3  4
 4  4  5  5

[:, :, 3] =
 4  5  4  5
 4  4  5  5

[:, :, 4] =
 2  3  2  3
 5  5  6  6

[:, :, 5] =
 3  4  3  4
 5  5  6  6

[:, :, 6] =
 4  5  4  5
 5  5  6  6
```

The (number of faces by number of elements) array `elemtobndy` gives the
boundary number for each face of each element.  A zero will be given for
connected faces.
```jldoctest brickmesh
julia> elemtobndy
4×6 Array{Int64,2}:
 1  0  0  1  0  0
 0  0  2  0  0  2
 0  0  0  0  0  0
 0  0  0  0  0  0
```
Note that the faces are listed in Cartesian order.

Finally, the periodic face connections are given in `faceconnections` which is a
list of arrays, one for each connection.
Each array in the list is given in the format `[e, f, vs...]` where
 - `e`  is the element number;
 - `f`  is the face number; and
 - `vs` is the global vertices that face associated with.
I the example
```jldoctest brickmesh
julia> faceconnections
3-element Array{Array{Int64,1},1}:
 [4, 4, 1, 2]
 [5, 4, 2, 3]
 [6, 4, 3, 4]
```
we see that face `4` of element `5` is associated with vertices `[2 3]` (the
vertices for face `1` of element `2`).
"""
function brickmesh(x, periodic; part=1, numparts=1,
                   boundary=ntuple(j->(1,1), length(x)))
  if boundary isa Matrix
    boundary = tuple(mapslices(x -> tuple(x...), boundary, dims=1)...)
  end

  @assert length(x) == length(periodic)
  @assert length(x) >= 1
  @assert 1 <= part <= numparts

  T = promote_type(eltype.(x)...)
  d = length(x)
  nvert = 2^d
  nface = 2d

  nelemdim = length.(x).-1
  elemlocal = linearpartition(prod(nelemdim), part, numparts)

  elemtovert = Array{Int}(undef, nvert, length(elemlocal))
  elemtocoord = Array{T}(undef, d, nvert, length(elemlocal))
  elemtobndy = zeros(Int, nface, length(elemlocal))
  faceconnections = Array{Array{Int, 1}}(undef, 0)

  verts = LinearIndices(ntuple(j->1:length(x[j]), d))
  elems = CartesianIndices(ntuple(j->1:length(x[j])-1, d))

  p = reshape(1:nvert, ntuple(j->2, d))
  fmask = hcat((p[ntuple(j->(j==div(f-1,2)+1) ? (mod(f-1,2)+1:mod(f-1,2)+1) :
                         (:), d)...,][:]
                for f=1:nface)...)

  for (e, ec) = enumerate(elems[elemlocal])
    corners = CartesianIndices(ntuple(j->ec[j]:ec[j]+1, d))
    for (v, vc) = enumerate(corners)
      elemtovert[v,e] = verts[vc]

      for j = 1:d
        elemtocoord[j,v,e] = x[j][vc[j]]
      end
    end

    for i=1:d
      if !periodic[i] && ec[i]==1
        elemtobndy[2(i-1)+1,e] = boundary[i][1]
      end
      if !periodic[i] && ec[i]==nelemdim[i]
        elemtobndy[2(i-1)+2,e] = boundary[i][2]
      end
    end

    for i=1:d
      if periodic[i] && ec[i]==nelemdim[i]
        neighcorners = CartesianIndices(ntuple(j->(i==j) ?
                                               (1:2) : (ec[j]:ec[j]+1), d))
        push!(faceconnections,
              vcat(e, 2i, verts[neighcorners[fmask[:,2i-1]]]))
      end
    end
  end

  (elemtovert, elemtocoord, elemtobndy, faceconnections)
end

"""
    parallelsortcolumns(comm::MPI.Comm, A;
                        alg::Base.Sort.Algorithm=Base.Sort.DEFAULT_UNSTABLE,
                        lt=isless,
                        by=identity,
                        rev::Union{Bool,Nothing}=nothing)

Sorts the columns of the distributed matrix `A`.

See the documentation of `sort!` for a description of the keyword arguments.

This function assumes `A` has the same number of rows on each MPI rank but can
have a different number of columns.
"""
function parallelsortcolumns(comm::MPI.Comm, A;
                             alg::Base.Sort.Algorithm=Base.Sort.DEFAULT_UNSTABLE,
                             lt=isless,
                             by=identity,
                             rev::Union{Bool,Nothing}=nothing)

  m, n = size(A)
  T = eltype(A)

  csize = MPI.Comm_size(comm)
  crank = MPI.Comm_rank(comm)
  croot = 0

  A = sortslices(A, dims=2, alg=alg, lt=lt, by=by, rev=rev)

  npivots = clamp(n, 0, csize)
  pivots = T[A[i, div(n*p,npivots)+1] for i=1:m, p=0:npivots-1]
  pivotcounts = MPI.Allgather(Cint(length(pivots)), comm)
  pivots = MPI.Allgatherv(pivots, pivotcounts, comm)
  pivots = reshape(pivots, m, div(length(pivots),m))
  pivots = sortslices(pivots, dims=2, alg=alg, lt=lt, by=by, rev=rev)

  # if we don't have any pivots then we must have zero columns
  if size(pivots) == (m, 0)
    return A
  end

  pivots =
    [pivots[i, div(div(length(pivots),m)*r,csize)+1] for i=1:m, r=0:csize-1]

  cols = map(i->view(A,:,i), 1:n)
  sendstarts = [(i<=csize) ? (searchsortedfirst(cols, pivots[:,i], lt=lt,
                                                by=by, rev=rev)-1)*m+1 : n*m+1
                for i=1:csize+1]
  sendcounts = [Cint(sendstarts[i+1]-sendstarts[i]) for i=1:csize]

  B = []
  for r = 0:csize-1
    counts = MPI.Allgather(sendcounts[r+1], comm)
    c = MPI.Gatherv(view(A, sendstarts[r+1]:sendstarts[r+2]-1), counts, r,
                    comm)
    if r == crank
      B = c
    end
  end
  B = reshape(B, m, div(length(B),m))

  sortslices(B, dims=2, alg=alg, lt=lt, by=by, rev=rev)
end

"""
    getpartition(comm::MPI.Comm, elemtocode)

Returns an equally weighted partition of a distributed set of elements by
sorting their codes given in `elemtocode`.

The codes for each element, `elemtocode`, are given as an array with a single
entry per local element or as a matrix with a column for each local element.

The partition is returned as a tuple three parts:

 - `partsendorder`: permutation of elements into sending order
 - `partsendstarts`: start entries in the send array for each rank
 - `partrecvstarts`: start entries in the receive array for each rank

Note that both `partsendstarts` and `partrecvstarts` are of length
`MPI.Comm_size(comm)+1` where the last entry has the total number of elements
to send or receive, respectively.
"""
getpartition(comm::MPI.Comm, elemtocode::AbstractVector) =
  getpartition(comm, reshape(elemtocode, 1, length(elemtocode)))

function getpartition(comm::MPI.Comm, elemtocode::AbstractMatrix)
  (ncode, nelem) = size(elemtocode)

  csize = MPI.Comm_size(comm)
  crank = MPI.Comm_rank(comm)

  CT = eltype(elemtocode)

  A = CT[elemtocode;                                 # code
         collect(CT, 1:nelem)';                      # original element number
         fill(CT(MPI.Comm_rank(comm)), (1, nelem));  # original rank
         fill(typemax(CT), (1, nelem))]              # new rank
  m, n = size(A)

  # sort by just code
  A = parallelsortcolumns(comm, A)

  # count the distribution of A
  counts = MPI.Allgather(last(size(A)), comm)
  starts = ones(Int, csize+1)
  for i=1:csize
    starts[i+1] = counts[i] + starts[i]
  end

  # loop to determine new rank
  j = range(starts[crank+1], stop=starts[crank+2]-1)
  for r = 0:csize-1
    k = linearpartition(starts[end]-1, r+1, csize)
    o = intersect(k,j) .- (starts[crank+1]-1)
    A[ncode+3,o] .= r
  end

  # sort by original rank and code
  A = sortslices(A, dims=2, by=x->x[[ncode+2,(1:ncode)...]])

  # count number of elements that are going to be sent
  sendcounts = zeros(Cint, csize)
  for i = 1:last(size(A))
    sendcounts[A[ncode+2,i]+1] += m
  end
  sendstarts = ones(Int, csize+1)
  for i=1:csize
    sendstarts[i+1] = sendcounts[i] + sendstarts[i]
  end

  # communicate columns of A to original rank
  B = []
  for r = 0:csize-1
    rcounts = MPI.Allgather(sendcounts[r+1], comm)
    c = MPI.Gatherv(view(A, sendstarts[r+1]:sendstarts[r+2]-1), rcounts, r,
                    comm)
    if r == crank
      B = c
    end
  end
  B = reshape(B, m, div(length(B),m))

  # check to make sure we didn't drop any elements
  @assert nelem == n == size(B)[2]

  partsendcounts = zeros(Cint, csize)
  for i = 1:last(size(B))
    partsendcounts[B[ncode+3,i]+1] += 1
  end
  partsendstarts = ones(Int, csize+1)
  for i=1:csize
    partsendstarts[i+1] = partsendcounts[i] + partsendstarts[i]
  end

  partsendorder = Int.(B[ncode+1,:])

  partrecvcounts = Cint[]
  for r = 0:csize-1
    c = MPI.Gather(partsendcounts[r+1], r, comm)
    if r == crank
      partrecvcounts = c
    end
  end

  partrecvstarts = ones(Int, csize+1)
  for i=1:csize
    partrecvstarts[i+1] = partrecvcounts[i] + partrecvstarts[i]
  end

  partsendorder, partsendstarts, partrecvstarts
end

"""
    partition(comm::MPI.Comm, elemtovert, elemtocoord, elemtobndy,
              faceconnections)

This function takes in a mesh (as returned for example by `brickmesh`) and
returns a Hilbert curve based partitioned mesh.
"""
function partition(comm::MPI.Comm, elemtovert, elemtocoord, elemtobndy,
                   faceconnections)
  (d, nvert, nelem) = size(elemtocoord)

  csize = MPI.Comm_size(comm)
  crank = MPI.Comm_rank(comm)

  nface = 2d
  nfacevert = 2^(d-1)

  # Here we expand the list of face connections into a structure that is easy
  # to partition.  The cost is extra memory transfer.  If this becomes a
  # bottleneck something more efficient may be implemented.
  #
  elemtofaceconnect = zeros(eltype(eltype(faceconnections)), nfacevert, nface,
                            nelem)
  for fc in faceconnections
    elemtofaceconnect[:,fc[2],fc[1]] = fc[3:end]
  end

  elemtocode = centroidtocode(comm, elemtocoord; CT=UInt64)
  sendorder, sendstarts, recvstarts = getpartition(comm, elemtocode)

  elemtovert = elemtovert[:,sendorder]
  elemtocoord = elemtocoord[:,:,sendorder]
  elemtobndy = elemtobndy[:,sendorder]
  elemtofaceconnect = elemtofaceconnect[:,:,sendorder]

  newelemtovert = []
  newelemtocoord = []
  newelemtobndy = []
  newelemtofaceconnect = []
  for r = 0:csize-1
    sendrange = sendstarts[r+1]:sendstarts[r+2]-1
    rcounts = MPI.Allgather(Cint(length(sendrange)), comm)

    netv = MPI.Gatherv(view(elemtovert, :, sendrange), rcounts.*Cint(nvert),
                       r, comm)

    netc = MPI.Gatherv(view(elemtocoord, :, :, sendrange),
                       rcounts.*Cint(d*nvert), r, comm)

    netb = MPI.Gatherv(view(elemtobndy, :, sendrange), rcounts.*Cint(nface),
                       r, comm)

    netfc = MPI.Gatherv(view(elemtofaceconnect, :, :, sendrange),
                        rcounts.*Cint(nfacevert*nface), r, comm)

    if r == crank
      newelemtovert = netv
      newelemtocoord = netc
      newelemtobndy = netb
      newelemtofaceconnect = netfc
    end
  end

  newnelem = recvstarts[end]-1
  newelemtovert = reshape(newelemtovert, nvert, newnelem)
  newelemtocoord = reshape(newelemtocoord, d, nvert, newnelem)
  newelemtobndy = reshape(newelemtobndy, nface, newnelem)
  newelemtofaceconnect = reshape(newelemtofaceconnect, nfacevert, nface,
                                 newnelem)

  # reorder local elements based on code of new elements
  A = UInt64[centroidtocode(comm, newelemtocoord; CT=UInt64);
             collect(1:newnelem)']
  A = sortslices(A, dims=2)
  newsortorder = view(A,d+1,:)
  newelemtovert = newelemtovert[:,newsortorder]
  newelemtocoord = newelemtocoord[:,:,newsortorder]
  newelemtobndy = newelemtobndy[:,newsortorder]
  newelemtofaceconnect = newelemtofaceconnect[:,:,newsortorder]

  newfaceconnections = similar(faceconnections, 0)
  for e = 1:newnelem, f = 1:nface
    if newelemtofaceconnect[1,f,e] > 0
      push!(newfaceconnections, vcat(e, f, newelemtofaceconnect[:,f,e]))
    end
  end

  (newelemtovert, newelemtocoord, newelemtobndy, newfaceconnections)
end

"""
    minmaxflip(x, y)

Returns `x, y` sorted lowest to highest and a bool that indicates if a swap
was needed.
"""
minmaxflip(x, y) = y < x ? (y, x, true) : (x, y, false)

"""
    vertsortandorder(a)

Returns `(a)` and an ordering `o==0`.
"""
vertsortandorder(a) = ((a,), 1)

"""
    vertsortandorder(a, b)

Returns sorted vertex numbers `(a,b)` and an ordering `o` depending on the
order needed to sort the elements.  This ordering is given below including the
vetex ordering for faces.

    o=    0      1

        (a,b)  (b,a)

          a      b
          |      |
          |      |
          b      a
"""
function vertsortandorder(a, b)
  a, b, s1 = minmaxflip(a, b)
  o = s1 ? 2 : 1
  ((a, b), o)
end

"""
    vertsortandorder(a, b, c)

Returns sorted vertex numbers `(a,b,c)` and an ordering `o` depending on the
order needed to sort the elements.  This ordering is given below including the
vetex ordering for faces.

    o=     1         2         3         4         5         6

        (a,b,c)   (c,a,b)   (b,c,a)   (b,a,c)   (c,b,a)   (a,c,b)

          /c\\      /b\\      /a\\      /c\\      /a\\      /b\\
         /   \\    /   \\    /   \\    /   \\    /   \\    /   \\
        /a___b\\  /c___a\\  /b___c\\  /b___a\\  /c___b\\  /a___c\\
"""
function vertsortandorder(a, b, c)
  # Use a (Bose-Nelson Algorithm based) sorting network from
  # <http://pages.ripco.net/~jgamble/nw.html> to sort the vertices.
  b, c, s1 = minmaxflip(b, c)
  a, c, s2 = minmaxflip(a, c)
  a, b, s3 = minmaxflip(a, b)

  if     !s1 && !s2 && !s3
    o = 1
  elseif !s1 &&  s2 &&  s3
    o = 2
  elseif  s1 && !s2 &&  s3
    o = 3
  elseif !s1 && !s2 &&  s3
    o = 4
  elseif  s1 &&  s2 &&  s3
    o = 5
  elseif  s1 && !s2 && !s3
    o = 6
  else
    error("Problem finding vertex ordering $((a,b,c)) with flips
          $((s1,s2,s3))")
  end

  ((a, b, c), o)
end

"""
    vertsortandorder(a, b, c, d)

Returns sorted vertex numbers `(a,b,c,d)` and an ordering `o` depending on the
order needed to sort the elements.  This ordering is given below including the
vetex ordering for faces.

    o=   1      2      3      4      5      6      7      8

       (a,b,  (a,c,  (b,a,  (b,d,  (c,a,  (c,d,  (d,b,  (d,c,
        c,d)   b,d)   c,d)   a,c)   d,b)   a,b)   c,a)   b,a)

       c---d  b---d  c---d  a---c  d---b  a---b  c---a  b---a
       |   |  |   |  |   |  |   |  |   |  |   |  |   |  |   |
       a---b  a---c  b---a  b---d  c---a  c---d  d---b  d---c
"""
function vertsortandorder(a, b, c, d)
  # Use a (Bose-Nelson Algorithm based) sorting network from
  # <http://pages.ripco.net/~jgamble/nw.html> to sort the vertices.
  a, b, s1 = minmaxflip(a, b)
  c, d, s2 = minmaxflip(c, d)
  a, c, s3 = minmaxflip(a, c)
  b, d, s4 = minmaxflip(b, d)
  b, c, s5 = minmaxflip(b, c)

 if     !s1 && !s2 && !s3 && !s4 && !s5
   o = 1
 elseif !s1 && !s2 && !s3 && !s4 &&  s5
   o = 2
 elseif  s1 && !s2 && !s3 && !s4 && !s5
   o = 3
 elseif !s1 && !s2 &&  s3 &&  s4 &&  s5
   o = 4
 elseif  s1 &&  s2 && !s3 && !s4 &&  s5
   o = 5
 elseif !s1 && !s2 &&  s3 &&  s4 && !s5
   o = 6
 elseif  s1 &&  s2 &&  s3 &&  s4 &&  s5
   o = 7
 elseif  s1 &&  s2 &&  s3 &&  s4 && !s5
   o = 8
 else
    # FIXME: some possible orientations are missing since there are a total of
    # 24. Missing orientations:
    #=
       d---c  d---c  b---c  a---d  c---b  d---a  a---b  b---a
       |   |  |   |  |   |  |   |  |   |  |   |  |   |  |   |
       a---b  b---a  a---d  b---c  d---a  c---b  d---c  c---d

       c---b  d---b  b---d  b---c  c---a  d---a  a---d  a---c
       |   |  |   |  |   |  |   |  |   |  |   |  |   |  |   |
       a---d  a---c  c---a  d---a  b---d  b---c  c---b  d---b
    =#
    error("Problem finding vertex ordering $((a,b,c,d))
            with flips $((s1,s2,s3,s4,s5))")
 end

  ((a, b, c, d), o)
end

"""
    connectmesh(comm::MPI.Comm, elemtovert, elemtocoord, elemtobndy,
                faceconnections)

This function takes in a mesh (as returned for example by `brickmesh`) and
returns a connected mesh.  This returns a `NamedTuple` of:

 - `elems` the range of element indices
 - `realelems` the range of real (aka nonghost) element indices
 - `ghostelems` the range of ghost element indices
 - `ghostfaces` ghost element to face is received;
   `ghostfaces[f,ge] == true` if face `f` of ghost element `ge` is received.
 - `sendelems` an array of send element indices
 - `sendfaces` send element to face is sent;
   `sendfaces[f,se] == true` if face `f` of send element `se` is sent.
 - `elemtocoord` element to vertex coordinates; `elemtocoord[d,i,e]` is the
    `d`th coordinate of corner `i` of element `e`
 - `elemtoelem` element to neighboring element; `elemtoelem[f,e]` is the
   number of the element neighboring element `e` across face `f`.  If there is
   no neighboring element then `elemtoelem[f,e] == e`.
 - `elemtoface` element to neighboring element face; `elemtoface[f,e]` is the
   face number of the element neighboring element `e` across face `f`.  If
   there is no neighboring element then `elemtoface[f,e] == f`.
 - `elemtoordr` element to neighboring element order; `elemtoordr[f,e]` is the
   ordering number of the element neighboring element `e` across face `f`.  If
   there is no neighboring element then `elemtoordr[f,e] == 1`.
 - `elemtobndy` element to bounday number; `elemtobndy[f,e]` is the
   boundary number of face `f` of element `e`.  If there is a neighboring
   element then `elemtobndy[f,e] == 0`.
 - `nabrtorank` a list of the MPI ranks for the neighboring processes
 - `nabrtorecv` a range in ghost elements to receive for each neighbor
 - `nabrtosend` a range in `sendelems` to send for each neighbor

"""
function connectmesh(comm::MPI.Comm, elemtovert, elemtocoord, elemtobndy,
                     faceconnections; dim = size(elemtocoord,1))
  d = dim
  (coorddim, nvert, nelem) = size(elemtocoord)
  nface, nfacevert = 2d, 2^(d-1)

  p = reshape(1:nvert, ntuple(j->2, d))
  fmask = hcat((p[ntuple(j->(j==div(f-1,2)+1) ? (mod(f-1,2)+1:mod(f-1,2)+1) :
                         (:), d)...][:]
                for f=1:nface)...)

  csize = MPI.Comm_size(comm)
  crank = MPI.Comm_rank(comm)

  VT = eltype(elemtovert)
  A = Array{VT}(undef, nfacevert+8, nface*nelem)

  MR, ME, MF, MO, NR, NE, NF, NO = nfacevert .+ (1:8)
  for e = 1:nelem
    v = reshape(elemtovert[:,e], ntuple(j->2, d))
    for f = 1:nface
      j = (e-1)*nface + f
      fv, o = vertsortandorder(v[fmask[:,f]]...)
      A[1:nfacevert, j] .= fv
      A[MR, j] = crank
      A[ME, j] = e
      A[MF, j] = f
      A[MO, j] = o
      A[NR, j] = typemax(VT)
      A[NE, j] = typemax(VT)
      A[NF, j] = typemax(VT)
      A[NO, j] = typemax(VT)
    end
  end

  # use neighboring vertices for connected faces
  for fc in faceconnections
    e = fc[1]
    f = fc[2]
    v = fc[3:end]
    j = (e-1)*nface + f
    fv, o = vertsortandorder(v...)
    A[1:nfacevert, j] .= fv
    A[MO, j] = o
  end

  A = parallelsortcolumns(comm, A, by=x->x[1:nfacevert])
  m, n = size(A)

  # match faces
  j = 1
  while j <= n
    if j+1 <= n && A[1:nfacevert,j]==A[1:nfacevert,j+1]
      # found connected face
      A[NR:NO, j  ] = A[MR:MO, j+1]
      A[NR:NO, j+1] = A[MR:MO, j  ]
      j += 2
    else
      # found unconnect face
      A[NR:NO, j] = A[MR:MO, j]
      j += 1
    end
  end

  A = sortslices(A, dims=2, by=x->(x[MR],x[NR],x[ME],x[MF]))

  # count number of elements that are going to be sent
  sendcounts = zeros(Cint, csize)
  for i = 1:last(size(A))
    sendcounts[A[MR,i]+1] += m
  end
  sendstarts = ones(Int, csize+1)
  for i=1:csize
    sendstarts[i+1] = sendcounts[i] + sendstarts[i]
  end

  # communicate columns of A to original rank
  B = []
  for r = 0:csize-1
    rcounts = MPI.Allgather(sendcounts[r+1], comm)
    c = MPI.Gatherv(view(A, sendstarts[r+1]:sendstarts[r+2]-1), rcounts, r,
                    comm)
    if r == crank
      B = c
    end
  end
  B = reshape(B, m, nface*nelem)

  # get element sending information
  B = sortslices(B, dims=2, by=x->(x[NR],x[ME]))
  sendelems = Int[]
  counts = zeros(Int, csize+1)
  counts[1] = (last(size(B)) > 0) ? 1 : 0
  sr, se = -1, 0
  for i = 1:last(size(B))
    r, e = B[NR,i], B[ME,i]
    # See if we need to send element `e` to rank `r` and make sure that we
    # didn't already mark it for sending.
    if r != crank && !(sr == r && se == e)
      counts[r+2] += 1
      append!(sendelems, e)
      sr, se = r, e
    end
  end

  # Mark which faces need to be sent
  sendfaces = BitArray(undef, nface, length(sendelems))
  sendfaces .= false
  sr, se, n = -1, 0, 0
  for i = 1:last(size(B))
    r, e, f = B[NR,i], B[ME,i], B[MF,i]
    if r != crank
      if !(sr == r && se == e)
        n += 1
        sr, se = r, e
      end
      sendfaces[f, n] = true
    end
  end

  sendstarts = cumsum(counts)
  nabrtosendrank = Int[r for r = 0:csize-1
                       if sendstarts[r+2]-sendstarts[r+1] > 0]
  nabrtosend = UnitRange{Int}[(sendstarts[r+1]:(sendstarts[r+2]-1))
                              for r = 0:csize-1
                              if sendstarts[r+2]-sendstarts[r+1] > 0]

  # get element receiving information
  B = sortslices(B, dims=2, by=x->(x[NR],x[NE]))
  counts = zeros(Int, csize+1)
  counts[1] = (last(size(B)) > 0) ? 1 : 0
  sr, se = -1, 0
  nghost = 0
  for i = 1:last(size(B))
    r, e = B[NR,i], B[NE,i]
    if r != crank
      # Check to make sure we have not already marked the element for
      # receiving since we could be connected to the receiving element across
      # multiple faces.
      if !(sr == r && se == e)
        nghost += 1
        counts[r+2] += 1
        sr, se = r, e
      end
      B[NE,i] = nelem + nghost
    end
  end

  # Mark which faces will be received
  ghostfaces = BitArray(undef, nface, nghost)
  ghostfaces .= false
  sr, se, ge = -1, 0, 0
  for i = 1:last(size(B))
    r, e, f = B[NR,i], B[NE,i], B[NF,i]
    if r != crank
      if !(sr == r && se == e)
        ge += 1
        sr, se = r, e
      end
      B[NR,i] = crank
      ghostfaces[f, ge] = true
    end
  end

  recvstarts = cumsum(counts)
  nabrtorecvrank = Int[r for r = 0:csize-1
                       if recvstarts[r+2]-recvstarts[r+1] > 0]
  nabrtorecv = UnitRange{Int}[(recvstarts[r+1]:(recvstarts[r+2]-1))
                              for r = 0:csize-1
                              if recvstarts[r+2]-recvstarts[r+1] > 0]

  @assert nabrtorecvrank == nabrtosendrank
  nabrtorank = nabrtorecvrank

  elemtoelem = repeat((1:nelem+nghost)', nface, 1)
  elemtoface = repeat(1:nface, 1, nelem+nghost)
  elemtoordr = ones(Int, nface, nelem+nghost)

  if d == 2
    for i = 1:last(size(B))
      me, mf, mo = B[ME,i], B[MF,i], B[MO,i]
      ne, nf, no = B[NE,i], B[NF,i], B[NO,i]

      elemtoelem[mf,me] = ne
      elemtoface[mf,me] = nf
      elemtoordr[mf,me] = (no == mo ? 1 : 2)
    end
  else
    for i = 1:last(size(B))
      me, mf, mo = B[ME,i], B[MF,i], B[MO,i]
      ne, nf, no = B[NE,i], B[NF,i], B[NO,i]

      elemtoelem[mf,me] = ne
      elemtoface[mf,me] = nf
      if no != 1 || mo != 1
        error("TODO add support for other orientations")
      end
      elemtoordr[mf,me] = 1
    end
  end

  # fill the ghost values in elemtocoord
  newelemtocoord = similar(elemtocoord, coorddim, nvert, nelem+nghost)
  newelemtobndy = similar(elemtobndy, nface, nelem+nghost)

  sendelemtocoord = elemtocoord[:,:,sendelems]
  sendelemtobndy = elemtobndy[:,sendelems]

  crreq = [MPI.Irecv!(view(newelemtocoord,:,:,nelem.+er), r, 666, comm)
           for (r, er) = zip(nabrtorank, nabrtorecv)]

  brreq = [MPI.Irecv!(view(newelemtobndy,:,nelem.+er), r, 666, comm)
           for (r, er) = zip(nabrtorank, nabrtorecv)]

  csreq = [MPI.Isend(view(sendelemtocoord,:,:,es), r, 666, comm)
          for (r, es) = zip(nabrtorank, nabrtosend)]

  bsreq = [MPI.Isend(view(sendelemtobndy,:,es), r, 666, comm)
          for (r, es) = zip(nabrtorank, nabrtosend)]

  newelemtocoord[:, :,1:nelem] .= elemtocoord
  newelemtobndy[:,1:nelem] .= elemtobndy

  MPI.Waitall!([csreq;crreq;bsreq;brreq])

  (elems=1:(nelem+nghost),       # range of          element indices
   realelems=1:nelem,            # range of real     element indices
   ghostelems=nelem.+(1:nghost), # range of ghost    element indices
   ghostfaces=ghostfaces,
   sendelems=sendelems,          # array of send     element indices
   sendfaces=sendfaces,
   elemtocoord=newelemtocoord,   # element to vertex coordinates
   elemtoelem=elemtoelem,        # element to neighboring element
   elemtoface=elemtoface,        # element to neighboring element face
   elemtoordr=elemtoordr,        # element to neighboring element order
   elemtobndy=newelemtobndy,     # element to boundary number
   nabrtorank=nabrtorank,        # list of neighboring processes MPI ranks
   nabrtorecv=nabrtorecv,        # neighbor receive ranges into `ghostelems`
   nabrtosend=nabrtosend)        # neighbor send ranges into `sendelems`
end

"""
    mappings(N, elemtoelem, elemtoface, elemtoordr)

This function takes in a polynomial order `N` and parts of a mesh (as returned
from `connectmesh`) and returns index mappings for the element surface flux
computation.  The returned `Tuple` contains:

 - `vmapM` an array of linear indices into the volume degrees of freedom where
   `vmapM[:,f,e]` are the degrees of freedom indices for face `f` of element
    `e`.

 - `vmapP` an array of linear indices into the volume degrees of freedom where
   `vmapP[:,f,e]` are the degrees of freedom indices for the face neighboring
   face `f` of element `e`.
"""
function mappings(N, elemtoelem, elemtoface, elemtoordr)
  nface, nelem = size(elemtoelem)

  d = div(nface, 2)
  Np, Nfp = (N+1)^d, (N+1)^(d-1)

  p = reshape(1:Np, ntuple(j->N+1, d))
  fd(f) =   div(f-1,2)+1
  fe(f) = N*mod(f-1,2)+1
  fmask = hcat((p[ntuple(j->(j==fd(f)) ? (fe(f):fe(f)) : (:), d)...][:]
                for f=1:nface)...)

  vmapM = similar(elemtoelem, Nfp, nface, nelem)
  vmapP = similar(elemtoelem, Nfp, nface, nelem)

  for e1 = 1:nelem, f1 = 1:nface
    e2 = elemtoelem[f1,e1]
    f2 = elemtoface[f1,e1]
    o2 = elemtoordr[f1,e1]

    # TODO support different orientations
    @assert o2 == 1

    vmapM[:,f1,e1] .= Np*(e1-1) .+ fmask[:,f1]
    vmapP[:,f1,e1] .= Np*(e2-1) .+ fmask[:,f2]
  end

  (vmapM, vmapP)
end

"""
   commmapping(N, commelems, commfaces, nabrtocomm)

This function takes in a polynomial order `N` and parts of a mesh (as returned
from `connectmesh` such as `sendelems`, `sendfaces`, and `nabrtosend`) and
returns index mappings for the element surface flux parallel communcation.
The returned `Tuple` contains:

 - `vmapC` an array of linear indices into the volume degrees of freedom to be
   communicated.

 - `nabrtovmapC` a range in `vmapC` to communicate with each neighbor.
"""
function commmapping(N, commelems, commfaces, nabrtocomm)
  nface, nelem = size(commfaces)

  @assert nelem == length(commelems)

  d = div(nface, 2)
  Nq = N+1
  Np = (N+1)^d

  vmapC = similar(commelems, nelem*Np)
  nabrtovmapC = similar(nabrtocomm)

  i = 1
  e = 1
  for neighbor in 1:length(nabrtocomm)
    rbegin = i
    for ne in nabrtocomm[neighbor]
      ce = commelems[ne]

      # Whole element sending
      # for n = 1:Np
      #   vmapC[i] = (ce-1)*Np + n
      #   i += 1
      # end

      CI = CartesianIndices(ntuple(_->1:Nq, d))
      for (ci, li) in zip(CI, LinearIndices(CI))
        addpoint = false
        for j = 1:d
          addpoint |= (commfaces[2*(j-1)+1, e] && ci[j] == 1 ) ||
                      (commfaces[2*(j-1)+2, e] && ci[j] == Nq)
        end

        if addpoint
          vmapC[i] = (ce-1)*Np + li
          i += 1
        end
      end

      e += 1
    end
    rend = i-1

    nabrtovmapC[neighbor] = rbegin:rend
  end

  resize!(vmapC, i-1)

  (vmapC, nabrtovmapC)
end

end # module

max_stars_repo_name: akshaysridhar/CLIMA

max_issues_repo_name: akshaysridhar/CLIMA

max_forks_repo_name: akshaysridhar/CLIMA


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder
name = "Cairo"
version = v"1.16.0"

sources = [
    ArchiveSource("https://www.cairographics.org/releases/cairo-$(version).tar.xz",
                  "5e7b29b3f113ef870d1e3ecf8adf21f923396401604bda16d44be45e66052331"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/cairo-*/

# Because `zlib` doesn't have a proper `.pc` file, configure fails to find.
export CPPFLAGS="-I${prefix}/include"

if [[ "${target}" == *-apple-* ]]; then
    BACKEND_OPTIONS="--enable-quartz --enable-quartz-image --disable-xcb --disable-xlib"
elif [[ "${target}" == *-mingw* ]]; then
    BACKEND_OPTIONS="--enable-win32 --disable-xcb --disable-xlib"
elif [[ "${target}" == *-linux-* ]] || [[ "${target}" == *freebsd* ]]; then
    BACKEND_OPTIONS="--enable-xlib --enable-xcb --enable-xlib-xcb"
fi

./configure --prefix=${prefix} --host=${target} \
    --disable-static \
    --enable-ft \
    --enable-tee \
    --enable-svg \
    --enable-ps \
    --enable-pdf \
    --enable-gobject \
    --disable-dependency-tracking \
    ${BACKEND_OPTIONS}
make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

# The products that we will ensure are always built
products = [
    LibraryProduct("libcairo-gobject", :libcairo_gobject),
    LibraryProduct("libcairo-script-interpreter", :libcairo_script_interpreter),
    LibraryProduct("libcairo", :libcairo),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    BuildDependency("Xorg_xorgproto_jll"),
    Dependency("Glib_jll"),
    Dependency("Pixman_jll"),
    Dependency("libpng_jll"),
    Dependency("Fontconfig_jll"),
    Dependency("FreeType2_jll"),
    Dependency("Bzip2_jll"),
    Dependency("Xorg_libXext_jll"),
    Dependency("Xorg_libXrender_jll"),
    Dependency("LZO_jll"),
    Dependency("Zlib_jll"),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; preferred_gcc_version=v"8")

ext: jl

max_stars_repo_name: blenessy/Yggdrasil

max_issues_repo_name: blenessy/Yggdrasil

max_forks_repo_name: blenessy/Yggdrasil


text:
# Copyright 2018 Erik Edin
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

using Behavior
using Behavior.Gherkin

# In case you want a more lenient parser, you can do something like this.
# For instance, the below options allows the Given/When/Then steps to be in any order
# exitcode = runspec(; parseoptions=ParseOptions(allow_any_step_order=true)) ? 0 : - 1

const use_experimental = length(ARGS) == 1 && ARGS[1] == "--experimental"
const options = Gherkin.ParseOptions(use_experimental=use_experimental)

exitcode = runspec(parseoptions=options) ? 0 : - 1
exit(exitcode)


text:
using Test

using MarketData
using TimeSeries

@testset "plotrecipes" begin

@testset "the Candlestick constructor works" begin
    candlestick = TimeSeries.Candlestick(ohlcv)
    @test length(candlestick.time) == length(ohlcv)

    candlestick = TimeSeries.Candlestick(ohlc)
    @test length(candlestick.time) == length(ohlc)

    @test_throws ArgumentError TimeSeries.Candlestick(cl)
    @test_throws ArgumentError TimeSeries.Candlestick(op)
end

end # @testset "plotrecipes


text:
# wrapper for arules.R

function init_arules()
    R"""
    test1p1 <- 1 + 1
    """
    if 2 == (@rget test1p1)
        R"""
        rm(test1p1)                     # remove to avoid memory usage
        """
    else
        stop("You must have R installed in order to use this package")
    end

    R"""
    if (\"arules\" %in% installed.packages() == FALSE) {
        install.packages(\"arules\", dependencies = TRUE, repos = \"http://cran.r-project.org\")
    }
    library(arules)

    all_factors <- function(dat) {
        p <- ncol(dat)
        out <- dat
        for (j in 1:p) {
            out[, j] <- factor(dat[, j])
        }
        return(out)
    }

    character_columns <- function(dat) {
        p <- ncol(dat)
        for (j in 1:p) {
            if (is.factor(dat[, j])) {
                dat[, j] <- as.character(dat[, j])
            }

        }
        return(dat)
    }
    """
    return nothing
end

# init_arules()

"""
    sanitize_input!(dat::DataFrame)
This is a simple function that assures we don't have strings
with a " => " as part of them. This is something R outputs
when it creates rules, and we'll uses it's presence to split
rules generated by R into left-hand side (lhs) and rhs.
"""
function sanitize_input!(dat)
    n, p = size(dat)
    for j = 1:p
        if eltype(dat[:, j]) <: Number
            continue
        end
        for i = 1:n
            dat[i, j] = replace(dat[i, j], " => ", " ")
        end
    end
end

function split_rule!(dat)
    n = size(dat, 1)
    dat[:lhs] = Array{String,1}(n)
    dat[:rhs] = Array{String,1}(n)
    for i = 1:n
        dat[i, :lhs], dat[i, :rhs] = split(dat[i, :rules], " => ")
    end
end

"""
    apriori(dat::DataFrame, supp::Float64, conf::Float64)
The apriori function implements the a-priori algorithm for
association rule mining. It returns a DataFrame with rules
with minimum support greater than or equal to `supp` and
confidence greater than or equal to `conf`.
"""
function apriori(dat::DataFrame, supp = 0.2, conf = 0.01)
    @rput dat

    rcode = """
        dat <- all_factors(dat)
        transacts <- as(dat, \"transactions\")

        rules1 <- apriori(transacts,
                          parameter = list(supp = $supp,
                                           conf = $conf,
                                           target = \"rules\"),
                          control = list(verbose = FALSE))

        rules1 <- if (length(rules1) == 0) data.frame() else rules1
        rules1 <- character_columns(as(rules1, \"data.frame\"))
    """
    reval(rcode);
    rules_df = @rget rules1;             # get dataframe from R

    R"""
    rm(dat, transacts, rules1)           # clean up R environment
    """
    split_rule!(rules_df);
    rules_df = rules_df[:, [:lhs, :rhs, :support, :confidence, :lift]]
    rules_df
end

# d = readtable("./data/adult.csv")
# sanitize_input!(d)
# a = apriori(d);


text:
module Bitcoin

using BitConverter, Secp256k1, MerkleTrees
import Secp256k1: serialize
using SHA, Ripemd, Base58
using Dates: unix2datetime, datetime2unix, now
using Sockets
import Base.show
export CompactSizeUInt
export Tx, TxIn, TxOut, Script, BlockHeader,
       VersionMessage, GetHeadersMessage, GetDataMessage,
       Node, BloomFilter
export point2address, wif, parse, serialize, id, fee, sig_hash,
       evaluate, fetch, verify, txsigninput,
       h160_2_address, script2address,
       iscoinbase, coinbase_height,
       blockparse, target, difficulty, check_pow, txoutparse
export get_tx, get_headers, get_blockhashbyheight

include("helper.jl")
include("CompactSizeUInt.jl")
include("constants.jl")
include("address.jl")
include("op.jl")
include("script.jl")
include("tx.jl")
include("rpc/rest.jl")
include("Block.jl")
include("BloomFilter.jl")
include("network.jl")
include("Node.jl")
include("murmur3.jl")

end # module


text:
module B

using A

export BType

struct BType
	x::AType
end

end


text:
"""Abstract histogram thresholding algorithm type"""
abstract type Thresholding end

## histogram separation type
"Cluster separation parameters"
struct Separation{T<:AbstractFloat}
    "Separation depth (depth between separated histogram modes)"
    depth::T
    "Separation discriminability (width between separated histogram modes)"
    discriminability::T
    "Distance threshold value"
    threshold::T
    "Global minimum as histogram bin index"
    globalmin::Int
    "Minimal distance"
    mindist::T
    "Maximal distance"
    maxdist::T
    "Number of bins in the histogram"
    bins::Int
end
Separation() = Separation{Float64}(0.0, 0.0, 0.0, 0, 0.0, 0.0, 0)

# properties
"Returns separation criteria value which is product of depth and discriminability."
criteria(sep::Separation) = sep.discriminability*sep.depth
"Returns distance threshold value for separation calculated on histogram of distances. It is used to determine which points belong to formed cluster."
threshold(sep::Separation) = sep.threshold
Base.extrema(sep::Separation) = (sep.mindist, sep.maxdist)

function Base.show(io::IO, S::Separation)
    print(io, "Separation($(criteria(S)), θ=$(threshold(S)))")
end

"Linear manifold cluster"
mutable struct Manifold{T<:AbstractFloat}
    "Dimension of the manifold"
    d::Int
    "Translation vector that contains coordinates of the linear subspace origin"
    μ::Vector{T}
    "Matrix of basis vectors that span the linear manifold"
    basis::Matrix{T}
    "Indexes of points associated with this cluster"
    points::Vector{Int}
    "Orthogonal subspace distance threshold"
    θ::T
    "Linear manifold subspace distance threshold"
    σ::T
end
Manifold(d::Int, μ::Vector{T}, basis::Matrix{T}, pnts::Vector{Int}) where {T<:AbstractFloat} =
    Manifold(d, μ, basis, pnts, zero(T), zero(T))
Manifold{T}(d::Int, pnts::Vector{Int}) where {T<:AbstractFloat} = Manifold(d, zeros(T,d), zeros(T,d,d), pnts)
Manifold{T}(d::Int) where {T<:AbstractFloat} = Manifold{T}(d, Int[])
Manifold() = Manifold{Float64}(0)

function Manifold(data::AbstractMatrix{T}, pnts::Vector{Int}) where {T<:AbstractFloat}
    X = @view data[:, pnts]
    R = fit(PCA, X; pratio=1.0)
    Manifold(outdim(R), mean(R), projection(R), pnts)
end

Base.copy(M::Manifold) = Manifold(outdim(M),mean(M),projection(M),points(M),threshold(M)...)

# properties
"Returns the dimension of the linear manifold cluster."
outdim(M::Manifold) = M.d
"Returns the dimension of the observation space."
indim(M::Manifold) = length(M.μ)
"Return an array of cluster assignments."
points(M::Manifold) = M.points
"Returns the cluster thresholds."
threshold(M::Manifold) = (M.θ, M.σ)
"Returns the matrix with columns corresponding to orthonormal vectors that span the linear manifold."
projection(M::Manifold) = M.basis[:,1:M.d]
"Returns the translation vector `μ` which contains coordinates of the linear manifold origin."
mean(M::Manifold) = M.μ
"Returns the cluster size - a number of points assosiated with the cluster."
Base.size(M::Manifold) = length(points(M))
"""Checks if the manifold `M` contains a full basis."""
hasfullbasis(M::Manifold) = size(M.basis,1) == size(M.basis,2)

function Base.show(io::IO, M::Manifold)
    print(io, "Manifold (dim = $(outdim(M)), size = $(size(M)), (θ,σ)=$(threshold(M)))")
end
function Base.dump(io::IO, M::Manifold)
    show(io, M)
    println(io)
    println(io, "thresholds (θ,σ): $(threshold(M)) ")
    println(io, "translation (μ): ")
    Base.showarray(io, mean(M)', header=false, repr=false)
    println(io)
    println(io, "basis: ")
    Base.showarray(io, projection(M), header=false, repr=false)
end

"""Algorithm specific exception"""
struct LMCLUSException <: Exception
    msg::String
end


text:
# using Pkg
# Pkg.activate(normpath(joinpath(@__DIR__, "..")))

using Oceananigans
using Oceananigans.Units
using Bioceananigans

# define the size and max depth of the simulation
const Nz = 48 # number of points in z
const H = 1000 # maximum depth

# create the grid of the model
grid = RectilinearGrid(GPU(),
    size=(Nz),
    z=(H * cos.(LinRange(π/2,0,Nz+1)) .- H)meters,
    topology=(Flat, Flat, Bounded)
)

#--------------- NP Model

# constants for the NP model
const μ₀ = 1/day   # surface growth rate
const m = 0.015/day # mortality rate due to virus and zooplankton grazing
const Kw = 0.059 # meter^-1
const kn = 0.75
const kr = 0.5

#  https://doi.org/10.1029/2017GB005850
const chl2c = 0.06 # average value for winter in North Atlantic

const α = 0.0538/day

const average = :growth
const shading = true

# create the mld field that will be updated at every timestep
h = Field{Center, Center, Nothing}(grid)
light_growth = Field{Center, Center, Center}(grid)

# time evolution of shortwave radiation (North Atlantic)
@inline Lₒ(t) = 116.5 * sin( 2π * ( t / days + 50 ) / 375.7 + 1.3 ) + 132.3
# evolution of the available light at the surface
@inline light_function(t, z) = 0.43 * Lₒ(t) * exp( z * Kw )
# light profile
@inline light_growth_function(light) = μ₀ * ( light * α ) / sqrt( μ₀^2 + ( light * α )^2 )

# nitrate and ammonium limiting
@inline N_lim(N, Nr) = (N/(N+kn)) * (kr/(Nr+kr))
@inline Nr_lim(Nr) =  (Nr/(Nr+kr))

# functions for the NP model
@inline P_forcing(light_growth, P, N, Nr)  =   light_growth * (N_lim(N, Nr) + Nr_lim(Nr)) * P - m * P^2
@inline N_forcing(light_growth, P, N, Nr)  = - light_growth * N_lim(N, Nr) * P
@inline Nr_forcing(light_growth, P, N, Nr) = - light_growth * Nr_lim(Nr) * P + m * P^2

# functions for the NP model
@inline P_forcing(i, j, k, grid, clock, fields, p)  = @inbounds P_forcing(p.light_growth[i, j, k], fields.P[i, j, k], fields.N[i, j, k], fields.Nr[i, j, k])
@inline N_forcing(i, j, k, grid, clock, fields, p)  = @inbounds N_forcing(p.light_growth[i, j, k], fields.P[i, j, k], fields.N[i, j, k], fields.Nr[i, j, k])
@inline Nr_forcing(i, j, k, grid, clock, fields, p) = @inbounds Nr_forcing(p.light_growth[i, j, k], fields.P[i, j, k], fields.N[i, j, k], fields.Nr[i, j, k])

# using the functions to determine the forcing
P_dynamics = Forcing(P_forcing, discrete_form=true, parameters=(; light_growth))
N_dynamics = Forcing(N_forcing, discrete_form=true, parameters=(; light_growth))
Nr_dynamics = Forcing(Nr_forcing, discrete_form=true, parameters=(; light_growth))

#--------------- Instantiate Model

# create the model
model = NonhydrostaticModel(grid = grid,
                            advection = WENO5(),
                            timestepper = :RungeKutta3,
                            closure=ScalarDiffusivity(ν=1e-5, κ=1e-5),
                            tracers = (:b, :P, :N, :Nr),
                            buoyancy = BuoyancyTracer(),
                            forcing = (P=P_dynamics, N=N_dynamics, Nr=Nr_dynamics))

#--------------- Initial Conditions

const cz = -200 # max. depth of phytoplankton
const g = 9.82 # gravity
const ρₒ = 1026 # reference density

# background density profile based on Argo data
@inline bg(z) = 0.25*tanh(0.0027*(-653.3-z))-6.8*z/1e5+1027.56

@inline B(x, y, z) = -(g/ρₒ) * bg(z)

# initial phytoplankton profile
@inline P(x, y, z) = ifelse(z>cz, 0.4, 0)

# setting the initial conditions
set!(model; b=B, P=P, N=13, Nr=0)

#--------------- Simulation

# # create a simulation
simulation = Simulation(model, Δt = 1minutes, stop_time = 5days)

wizard = TimeStepWizard(cfl=1.0, max_change=1.1, max_Δt=1hour)
simulation.callbacks[:wizard] = Callback(wizard, IterationInterval(10))

# buoyancy decrease criterium for determining the mixed-layer depth
const Δb=(g/ρₒ) * 0.03
compute_mixed_layer_depth!(simulation) = MixedLayerDepth!(h, simulation.model.tracers.b, Δb)
# add the function to the callbacks of the simulation
simulation.callbacks[:compute_mld] = Callback(compute_mixed_layer_depth!)

# compute_light_growth!(simulation) = LightGrowth!(light_growth, h, simulation.model.tracers.P, light_function, light_growth_function, time(simulation), average, shading, chl2c)
# # add the function to the callbacks of the simulation
# simulation.callbacks[:compute_light_growth] = Callback(compute_light_growth!)

# zeroing negative values
zero_P(sim) = parent(sim.model.tracers.P) .= max.(0, parent(sim.model.tracers.P))
simulation.callbacks[:zero_P] = Callback(zero_P)

zero_N(sim) = parent(sim.model.tracers.N) .= max.(0, parent(sim.model.tracers.N))
simulation.callbacks[:zero_N] = Callback(zero_N)

zero_Nr(sim) = parent(sim.model.tracers.Nr) .= max.(0, parent(sim.model.tracers.Nr))
simulation.callbacks[:zero_Nr] = Callback(zero_Nr)

# # merge light and h to the outputs
# outputs = merge(model.velocities, model.tracers, (; light_growth, h)) # make a NamedTuple with all outputs

# # writing the output
# simulation.output_writers[:fields] =
#     NetCDFOutputWriter(model, outputs, filepath = "data/output_average-$(average)_shading-$(shading).nc",
#                        schedule=TimeInterval(3hours))

using Printf

function print_progress(simulation)
    b, P, N, Nr = simulation.model.tracers

    # Print a progress message
    msg = @sprintf("i: %04d, t: %s, Δt: %s, P_max = %.1e, N_min = %.1e, Nr_max = %.1e, wall time: %s\n",
                   iteration(simulation),
                   prettytime(time(simulation)),
                   prettytime(simulation.Δt),
                   maximum(P), minimum(N), maximum(Nr),
                   prettytime(simulation.run_wall_time))

    @info msg

    return nothing
end

simulation.callbacks[:progress] = Callback(print_progress, TimeInterval(1hour))

# run the simulation
run!(simulation)

max_stars_repo_name: iuryt/biocenanigans

max_issues_repo_name: iuryt/biocenanigans


text:
mutable struct TrueSkillEnv
    mu::Float64
    sigma::Float64
    beta::Float64
    tau::Float64
    draw_probability::Float64
end

GLOBAL_ENV = TrueSkillEnv(25.0, 25.0/3, 25.0/6, 25.0/300, 0.10)
DELTA = 0.0001

function calc_draw_probability(draw_margin, size; env=GLOBAL_ENV)
    return 2 * cdf(draw_margin / (sqrt(size) * env.beta)) - 1.
end

function calc_draw_margin(draw_probability, size; env=GLOBAL_ENV)
    return ppf((draw_probability + 1) / 2.) * sqrt(size) * env.beta
end

function _team_sizes(rating_groups)
    szs = map(x -> length(x), rating_groups)
    cumsum(szs)
end

function _floating_point_error()
    error("Floating point error. Result out of range.")
end

struct Rating <: GaussianItem
    pi::Float64
    tau::Float64
end

function Rating(; mu=GLOBAL_ENV.mu, sigma=GLOBAL_ENV.sigma)
    pi = sigma^(-2)
    tau = pi * mu
    Rating(pi, tau)
end

function Rating(g::Gaussian)
    Rating(g.pi, g.tau)
end

@inline function mu(r::Rating)
    if r.pi != 0
        return r.tau / r.pi
    else
        return 0
    end
end

@inline function sigma(r::Rating)
    if r.pi != 0
        return(sqrt(1/r.pi))
    else
        return Inf
    end
end

@inline function isless(a::Rating, b::Rating)
    isless(mu(a), mu(b))
end

function show(io::IO, r::Rating)
    print(io, "Rating(mu=$(mu(r)), sigma=$(sigma(r)))")
end

function create_rating(;mu=GLOBAL_ENV.mu, sigma=GLOBAL_ENV.sigma)
    Rating(;mu=mu, sigma=sigma)
end

function v_win(diff::Real, draw_margin::Real)
    x = diff - draw_margin
    denom = cdf(x)
    if denom != 0.0
        return pdf(x) / denom
    else
        return -x
    end
end

function v_draw(diff::Real, draw_margin::Real)
    abs_diff = abs(diff)
    a, b = draw_margin - abs_diff, -draw_margin - abs_diff
    denom = cdf(a) - cdf(b)
    numer = pdf(b) - pdf(a)
    c = 0.0
    if denom != 0.0
        c = numer / denom
    else
        c = a
    end
    d = 0
    if diff < 0
        d = -1
    else
        d = 1
    end
    return c * d
end

function w_win(diff::Real, draw_margin::Real)
    x = diff - draw_margin
    v = v_win(diff, draw_margin)
    w = v * (v + x)
    if 0 < w < 1
        return w
    else
        _floating_point_error()
    end
end

function w_draw(diff::Real, draw_margin::Real)
    abs_diff = abs(diff)
    a, b = draw_margin - abs_diff, -draw_margin - abs_diff
    denom = cdf(a) - cdf(b)
    if denom == 0
        _floating_point_error()
    end
    v = v_draw(abs_diff, draw_margin)
    return (v ^ 2) + (a * pdf(a) - b * pdf(b)) / denom
end

function validate_rating_groups(rating_groups::Vector{Vector{Rating}})
    if length(rating_groups) < 2
        error("need multiple groups")
    end
    if any(map(length, rating_groups) .== 0)
        error("no empty groups allowed")
    end
    return rating_groups
end

function factor_graph_builders(rating_groups::Vector{Vector{Rating}}, ranks, weights; env=GLOBAL_ENV)
    # This is the main builder!!
    flatten_ratings = reduce(vcat, rating_groups)
    flatten_weights = reduce(vcat, weights)
    size = length(flatten_ratings)
    group_size = length(rating_groups)

    rating_vars = [Variable() for x in 1:size]
    perf_vars = [Variable() for x in 1:size]
    team_perf_vars = [Variable() for x in 1:group_size]
    team_diff_vars = [Variable() for x in 1:(group_size - 1)]
    team_sizes = _team_sizes(rating_groups)

    # layer builders
    function build_rating_layer()
        r = Vector{PriorFactor}()
        for (rating_var, rating) in zip(rating_vars, flatten_ratings)
            push!(r, PriorFactor(rating_var, rating, env.tau))
        end
        return r
    end

    function build_perf_layer()
        r = Vector{LikelihoodFactor}()
        for (rating_var, perf_var) in zip(rating_vars, perf_vars)
            push!(r, LikelihoodFactor(rating_var, perf_var, env.beta ^ 2))
        end
        return r
    end

    function build_team_perf_layer()
        r = Vector{SumFactor}()
        for (team, team_perf_var) in enumerate(team_perf_vars)
            if team > 1
                first = team_sizes[team - 1] + 1
            else
                first = 1
            end
            last = team_sizes[team]
            child_perf_vars = perf_vars[first:last]
            coeffs = flatten_weights[first:last]
            push!(r, SumFactor(team_perf_var, child_perf_vars, coeffs))
        end
        return r
    end

    function build_team_diff_layer()
        r = Vector{SumFactor}()
        for (team, team_diff_var) in enumerate(team_diff_vars)
            push!(r, SumFactor(team_diff_var,
                team_perf_vars[team:team+1], [+1, -1]))
        end
        return r
    end

    function build_trunc_layer()
        r = Vector{TruncateFactor}()
        for (x, team_diff_var) in enumerate(team_diff_vars)
            draw_probability = env.draw_probability
            size = sum(map(length, rating_groups[x:x+1]))
            draw_margin = calc_draw_margin(draw_probability, size; env=env)
            if ranks[x] == ranks[x + 1] # in case of tie
                v_func, w_func = v_draw, w_draw
            else
                v_func, w_func = v_win, w_win
            end
            push!(r, TruncateFactor(team_diff_var, v_func, w_func, draw_margin))
        end
        return r
    end
    return (build_rating_layer, build_perf_layer, build_team_perf_layer,
        build_team_diff_layer, build_trunc_layer)
end

function run_schedule(build_rating_layer, build_perf_layer,
    build_team_perf_layer, build_team_diff_layer, build_trunc_layer; min_delta=DELTA)
    # This is the main runner!!
    if min_delta <= 0
        error("min_delta must be greater than 0")
    end
    layers = Vector{Vector{Node}}()

    rating_layer = build_rating_layer()
    perf_layer = build_perf_layer()
    team_perf_layer = build_team_perf_layer()

    push!(layers, rating_layer)
    push!(layers, perf_layer)
    push!(layers, team_perf_layer)

    for f in reduce(vcat, layers)
        down(f)
    end

    team_diff_layer = build_team_diff_layer()
    trunc_layer = build_trunc_layer()

    push!(layers, team_diff_layer)
    push!(layers, trunc_layer)

    team_diff_len = length(team_diff_layer)

    for i in 1:10
        if team_diff_len == 1
            # two teams
            down(team_diff_layer[1])
            delta = up(trunc_layer[1])
        else
            # three or more teams
            delta = 0.0
            for j in 1:(team_diff_len - 1)
                down(team_diff_layer[j])
                delta = max(delta, up(trunc_layer[j]))
                up(team_diff_layer[j]; index=2) # up to right Variable
            end
            for j in team_diff_len:-1:2
                down(team_diff_layer[j])
                delta = max(delta, up(trunc_layer[j]))
                up(team_diff_layer[j]; index=1) # up to left Variable
            end
        end
        if delta <= min_delta
            break
        end
    end

    # up both ends
    up(team_diff_layer[1]; index=1)
    up(team_diff_layer[team_diff_len]; index=2)

    # up the remainder of the black arrows
    for f in team_perf_layer
        for i in 1:(length(f.vars) - 1) # number of members.
            up(f; index=i)
        end
    end
    for f in perf_layer
        up(f)
    end
    return layers
end

function rate(rating_groups::Vector{Vector{Rating}}; ranks=1:length(rating_groups), weights=[ones(length(x)) for x in rating_groups], min_delta=DELTA, env=GLOBAL_ENV)
    validate_rating_groups(rating_groups)
    group_size = length(rating_groups)
    if length(ranks) != group_size
        error("wrong ranks")
    end

    # sort rating groups by ranks
    by_rank = x -> x[2][2]
    sorting = sort(collect(enumerate(zip(rating_groups, ranks, weights))), by=by_rank)
    sorted_rating_groups = Vector{Vector{Rating}}()
    sorted_ranks = Vector{Float64}()
    sorted_weights = Vector{Vector{Float64}}()
    for (x, (g, r, w)) in sorting
        push!(sorted_rating_groups, g)
        push!(sorted_ranks, r)
        push!(sorted_weights, [max(min_delta, w_) for w_ in w])
    end
    b1, b2, b3, b4, b5 = factor_graph_builders(sorted_rating_groups, sorted_ranks, sorted_weights; env=env)
    layers = run_schedule(b1, b2, b3, b4, b5; min_delta=min_delta)

    rating_layer, team_sizes = layers[1], _team_sizes(sorted_rating_groups)
    transformed_groups = Vector{Vector{Rating}}()
    for (firstm1, last) in zip(vcat([0], team_sizes[1:end-1]), team_sizes)
        group = Vector{Rating}()
        for f in rating_layer[(firstm1+1):last]
            push!(group, Rating(; mu=mu(var(f)), sigma=sigma(var(f))))
        end
        push!(transformed_groups, group)
    end

    by_hint = x -> x[1]
    unsorting = sort(collect(zip((x for (x, _) in sorting), transformed_groups)), by=by_hint)
    return [g for (x, g) in unsorting]
end

function expose(rating::Rating; env::TrueSkillEnv=GLOBAL_ENV)
    k = env.mu / env.sigma
    return mu(rating) - k * sigma(rating)
end

function make_as_global(env::TrueSkillEnv)
    GLOBAL_ENV = env
end

function rate(rating1::Rating, rating2::Rating; drawn=false, min_delta=DELTA)
    if drawn
        ranks = [0, 0]
    else
        ranks = [0 ,1]
    end
    teams = rate([[rating1], [rating2]]; ranks=ranks,  min_delta = min_delta)
    return teams[1][1], teams[2][1]
end

function rate(ratings::Vector{Rating}; ranks=1:length(ratings), min_delta=DELTA)
    teams = rate(map(x -> [x], ratings); ranks=ranks, min_delta=min_delta)
    return map(x -> x[1], teams)
end


text:
mutable struct VLGRNSentence

    original_sentence::String
    sentence_actor_clause::String
    sentence_action_clause::String
    sentence_target_clause::String
    sentence_delimiter::Char

    function VLGRNSentence()
      this = new()
    end
end

mutable struct VLMetabolicReaction

    record::String
    reaction_name::String
    ec_number::String
    left_phrase::String
    right_phrase::String
    reversible::String

    function VLMetabolicReaction()
		  this = new()
	end
end

mutable struct VLMetabolicSpeciesReference

  species_symbol::String
  stoichiometry::String

  function VLMetabolicSpeciesReference()
    this = new()
  end
end

mutable struct VLSequenceRecord

  record::String
  operationType::Symbol
  sequence::String
  enzymeSymbol::String
  modelSpecies::String
  length::Int

  function VLSequenceRecord()
    this = new()
  end
end

mutable struct VLSpeciesSymbol

  species_type::Symbol
  species_symbol::String

  function VLSpeciesSymbol()
		this = new()
	end
end

mutable struct VLProgramComponent

  filename::String
  buffer::String
  vector::Array{Float64,1}
  matrix::Array{Float64,2}
  type::Symbol

  function VLProgramComponent()
    this = new()
  end

end

mutable struct VLGRNConnectionObject

  connection_symbol::String
  connection_actor_set::Array{VLSpeciesSymbol,1}
  connection_target_set::Array{VLSpeciesSymbol,1}
  connection_type::Symbol

  function VLGRNConnectionObject()
    this = new()
  end
end

max_stars_repo_name: varnerlab/CFMG

max_issues_repo_name: varnerlab/CFMG

max_forks_repo_name: varnerlab/CFMG


text:
# Licensed to the Apache Software Foundation (ASF) under one
# or more contributor license agreements.  See the NOTICE file
# distributed with this work for additional information
# regarding copyright ownership.  The ASF licenses this file
# to you under the Apache License, Version 2.0 (the
# "License"); you may not use this file except in compliance
# with the License.  You may obtain a copy of the License at
#
#   http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing,
# software distributed under the License is distributed on an
# "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied.  See the License for the
# specific language governing permissions and limitations
# under the License.

module TestNameManager
using MXNet
using Base.Test

function test_default()
  info("NameManager::default")

  name = :_____aaaaa_____
  @test get!(mx.DEFAULT_NAME_MANAGER, name, "") == name
  @test get!(mx.DEFAULT_NAME_MANAGER, string(name), "") == name

  hint = name
  @test get!(mx.DEFAULT_NAME_MANAGER, "", hint) == Symbol("$(hint)0")
  @test get!(mx.DEFAULT_NAME_MANAGER, "", string(hint)) == Symbol("$(hint)1")
end

function test_prefix()
  info("NameManager::prefix")

  name   = :_____bbbbb_____
  prefix = :_____foobar_____

  prefix_manager = mx.PrefixNameManager(prefix)
  @test get!(prefix_manager, name, "") == Symbol("$prefix$name")
  @test get!(prefix_manager, "", name) == Symbol("$prefix$(name)0")
end

@testset "Name Test" begin
  test_default()
  test_prefix()
end

end

max_stars_repo_name: Vikas-kum/incubator-mxnet

max_issues_repo_name: Vikas-kum/incubator-mxnet

max_forks_repo_name: Vikas-kum/incubator-mxnet


text:
"""
    The bare minimum algorithm of CG for forward error analysis.
    It's for the forward error analysis package.
"""
function cgs(
        A::AbstractMatrix,
        b::Abstractvector;
        x0=nothing,
        k::Int=nothing
    )
    if k === nothing
        k = length(b)
    end
    if x0 === nothing
        x0 = ones(b|>length)
    end
    r = b - A*x0
    x = x0
    p = r
    for II = 1:k
        Ap = A*p
        a = dot(r, r)/dot(p, Ap)
        b = 1/dot(r, r)
        x += x*p
        r -= -a*Ap
        b *= dot(r, r)
        p = r + b*p
    end
    return x
end


text:
@testset "1.0 always to 1.0" begin
    N = 10000

    for i = 1:N
        @test 1.0f0 == Float32(Float32sr(1.0))
        @test 1.0f0 == Float32(Float32_stochastic_round(1.0))
    end
end

@testset "Odd floats are never round away" begin
    # chances are 2^-28 = 4e-9
    # so technically one should use N > 100_000_000 but takes too long
    N = 100_000

    f_odd = prevfloat(one(Float32sr))
    f_odd_f64 = Float64(f_odd)
    for _ = 1:N
        @test f_odd == Float32sr(f_odd_f64)
        @test f_odd == Float32_stochastic_round(f_odd_f64)
    end
end

@testset "Even floats are never round away" begin
    # chances are 2^-28 = 4e-9
    # so technically one should use N > 100_000_000 but takes too long
    N = 100_000

    f_odd = prevfloat(prevfloat(one(Float32sr)))
    f_odd_f64 = Float64(f_odd)
    for _ = 1:N
        @test f_odd == Float32sr(f_odd_f64)
        @test f_odd == Float32_stochastic_round(f_odd_f64)
    end
end

@testset "Sign flip" begin
    @test one(Float32sr) == -(-(one(Float32sr)))
    @test zero(Float32sr) == -(zero(Float32sr))
end

@testset "Integer promotion" begin
    f = Float32sr(1)
    @test 2f == Float32sr(2)
    @test 0 == Float32sr(0)
end

@testset "Rounding" begin
    @test 1 == Int(round(Float32sr(1.2)))
    @test 1 == Int(floor(Float32sr(1.2)))
    @test 2 == Int(ceil(Float32sr(1.2)))

    @test -1 == Int(round(Float32sr(-1.2)))
    @test -2 == Int(floor(Float32sr(-1.2)))
    @test -1 == Int(ceil(Float32sr(-1.2)))
end

@testset "Nextfloat prevfloat" begin
    o = one(Float32sr)
    @test o == nextfloat(prevfloat(o))
    @test o == prevfloat(nextfloat(o))
end

@testset "Comparisons" begin
    @test Float32sr(1f0) <  Float32sr(2f0)
    @test Float32sr(1)   <  Float32sr(2)
    @test Float32sr(1.0) <  Float32sr(2.0)
    @test Float32sr(1)   <= Float32sr(2)
    @test Float32sr(1f0) <= Float32sr(2f0)
    @test Float32sr(1.0) <= Float32sr(2.0)
    @test Float32sr(2)   >  Float32sr(1)
    @test Float32sr(2f0) >  Float32sr(1f0)
    @test Float32sr(2)   >= Float32sr(1)
    @test Float32sr(2f0) >= Float32sr(1f0)
    @test Float32sr(2.0) >= Float32sr(1.0)
end

@testset "NaN and Inf" begin
    @test isnan(NaN32sr)
    @test ~isfinite(NaN32sr)
    @test ~isfinite(Inf32sr)
    @test isnan(Float32_stochastic_round(NaN))
    @test Inf32 == Float32(Float32_stochastic_round(Inf))
    @test -Inf32 == Float32(Float32_stochastic_round(-Inf))

    N = 1000
    for i in 1:N
        @test Inf32sr == Float32_stochastic_round(Inf)
        @test -Inf32sr == Float32_stochastic_round(-Inf)
        @test isnan(Float32_stochastic_round(NaN))
    end
end

@testset "No stochastic round to NaN" begin
    f1 = nextfloat(0.0)
    f2 = prevfloat(0.0)
    for i in 1:1000
        @test isfinite(Float32_stochastic_round(f1))
        @test isfinite(Float32_stochastic_round(f2))
    end
end

@testset "No stochastic round to NaN" begin
    f1 = Float64(floatmax(Float32))   # larger can map to Inf though!
    f2 = -f1
    for i in 1:1000
        @test isfinite(Float32_stochastic_round(f1))
        @test isfinite(Float32_stochastic_round(f2))
    end
end

function test_chances_round(f64::Float64;N::Int=100_000)
    p = Float32_chance_roundup(f64)

    f32_round = Float32sr(f64)
    if Float64(f32_round) <= f64
        f32_rounddown = f32_round
        f32_roundup = nextfloat(f32_round)
    else
        f32_roundup = f32_round
        f32_rounddown = prevfloat(f32_round)
    end

    Ndown = 0
    Nup = 0
    for _ in 1:N
        f32 = Float32_stochastic_round(f64)
        if f32 == f32_rounddown
            Ndown += 1
        elseif f32 == f32_roundup
            Nup += 1
        end
    end

    test1 = Ndown + Nup == N
    test2 = isapprox(Ndown/N,1-p,atol=1e-2)
    test3 = isapprox(Nup/N,p,atol=1e-2)

    return test1 && test2 && test3
end

@testset "Test for N(0,1)" begin
    for x in randn(10_000)
        @test test_chances_round(x)
    end
end

@testset "Test for U(0,1)" begin
    for x in rand(10_000)
        @test test_chances_round(x)
    end
end

@testset "Test for powers of two" begin
    for x in Float64[2,4,8,16,32,64,128,256]
        @test test_chances_round(x)
    end

    for x in -Float64[2,4,8,16,32,64,128,256]
        @test test_chances_round(x)
    end

    for x in Float64[1/2,1/4,1/8,1/16,1/32,1/64,1/128,1/256]
        @test test_chances_round(x)
    end

    for x in -Float64[1/2,1/4,1/8,1/16,1/32,1/64,1/128,1/256]
        @test test_chances_round(x)
    end
end

@testset "Test for subnormals" begin

    floatminF32 = reinterpret(UInt64,Float64(floatmin(Float32)))
    minposF32 = reinterpret(UInt64,Float64(nextfloat(zero(Float32))))

    N = 100
    subnormals = reinterpret.(Float64,rand(minposF32:floatminF32,N))
    for x in subnormals
        @test test_chances_round(x)
    end

    for x in -subnormals
        @test test_chances_round(x)
    end
end


text:
module EGraphs

include("../docstrings.jl")

using RuntimeGeneratedFunctions
RuntimeGeneratedFunctions.init(@__MODULE__)

# import ..Rule
# import ..getrhsfun

using TermInterface
using ..Util
using ..Patterns
using ..Rules
import ..@log

include("enode.jl")
export ENode
export EClassId

include("abstractanalysis.jl")
export AbstractAnalysis

include("eclass.jl")
export EClass
export hasdata
export getdata
export setdata!

include("intdisjointmap.jl")
export IntDisjointSet
export in_same_set

include("metadata_analysis.jl")
export MetadataAnalysis

include("egraph.jl")
export find
export lookup
export geteclass
export arity
export EGraph
export merge!
export in_same_class
export addexpr!
export rebuild!
export settermtype!
export gettermtype

include("analysis.jl")
export analyze!

include("subst.jl")
export Sub

include("ematch.jl")
include("Schedulers/Schedulers.jl")
export Schedulers
using .Schedulers

include("saturation/goal.jl")
export SaturationGoal
export EqualityGoal
export reached
include("saturation/reason.jl")
export ReportReasons
include("saturation/report.jl")
include("saturation/params.jl")
export SaturationParams
include("saturation/search.jl")
include("saturation/apply.jl")
include("saturation/saturation.jl")
export saturate!

include("equality.jl")
export areequal
export @areequal
export @areequalg

include("extraction.jl")
export extract!
export ExtractionAnalysis
export astsize
export astsize_inv
export @extract

end


text:
# This file contains functions for the Additive DEA model
"""
    AdditivelDEAModel
An data structure representing an additive DEA model.
"""
struct AdditiveDEAModel <: AbstractTechnicalDEAModel
    n::Int64
    m::Int64
    s::Int64
    weights::Symbol
    orient::Symbol
    rts::Symbol
    disposX::Symbol
    disposY::Symbol
    dmunames::Union{Vector{String},Nothing}
    eff::Vector
    slackX::Matrix
    slackY::Matrix
    lambda::SparseMatrixCSC{Float64, Int64}
    Xtarget::Matrix
    Ytarget::Matrix
end

"""
    deaadd(X, Y, model)
Compute related data envelopment analysis weighted additive models for inputs `X` and outputs `Y`.

Model specification:
- `:Ones`: standard additive DEA model.
- `:MIP`: Measure of Inefficiency Proportions. (Charnes et al., 1987; Cooper et al., 1999)
- `:Normalized`: Normalized weighted additive DEA model. (Lovell and Pastor, 1995)
- `:RAM`: Range Adjusted Measure. (Cooper et al., 1999)
- `:BAM`: Bounded Adjusted Measure. (Cooper et al, 2011)
- `:Custom`: User supplied weights.

# Optional Arguments
- `orient=:Graph`: choose between graph oriented `:Graph`, input oriented `:Input`, or output oriented model `:Output`.
- `rts=:VRS`: choose between constant returns to scale `:CRS` or variable returns to scale `:VRS`.
- `rhoX`: matrix of weights of inputs. Only if `model=:Custom`.
- `rhoY`: matrix of weights of outputs. Only if `model=:Custom`.
- `Xref=X`: Identifies the reference set of inputs against which the units are evaluated.
- `Yref=Y`: Identifies the reference set of outputs against which the units are evaluated.
- `disposX=:Strong`: chooses strong disposability of inputs. For weak disposability choose `:Weak`.
- `disposY=:Strong`: chooses strong disposability of outputs. For weak disposability choose `:Weak`.
- `names`: a vector of strings with the names of the decision making units.

# Examples
```jldoctest
julia> X = [5 13; 16 12; 16 26; 17 15; 18 14; 23 6; 25 10; 27 22; 37 14; 42 25; 5 17];

julia> Y = [12; 14; 25; 26; 8; 9; 27; 30; 31; 26; 12];

julia> deaadd(X, Y, :MIP)
Weighted Additive DEA Model
DMUs = 11; Inputs = 2; Outputs = 1
Orientation = Graph; Returns to Scale = VRS
Weights = MIP
─────────────────────────────────────────────────────
      efficiency       slackX1  slackX2       slackY1
─────────────────────────────────────────────────────
1    0.0           0.0              0.0   0.0
2    0.507519      0.0              0.0   7.10526
3    0.0           0.0              0.0   0.0
4   -4.72586e-17  -8.03397e-16      0.0   0.0
5    2.20395       0.0              0.0  17.6316
6    1.31279e-16   8.10382e-16      0.0   8.64407e-16
7    0.0           0.0              0.0   0.0
8    0.0           0.0              0.0   0.0
9    0.0           0.0              0.0   0.0
10   1.04322      17.0             15.0   1.0
11   0.235294      0.0              4.0   0.0
─────────────────────────────────────────────────────
```
"""
function deaadd(X::Union{Matrix,Vector}, Y::Union{Matrix,Vector},
    model::Symbol = :Default; orient::Symbol = :Graph,
    rts::Symbol = :VRS,
    rhoX::Union{Matrix,Vector,Nothing} = nothing, rhoY::Union{Matrix,Vector,Nothing} = nothing,
    Xref::Union{Matrix,Vector,Nothing} = nothing, Yref::Union{Matrix,Vector,Nothing} = nothing,
    disposX::Symbol = :Strong, disposY::Symbol = :Strong,
    names::Union{Vector{String},Nothing} = nothing,
    optimizer::Union{DEAOptimizer,Nothing} = nothing)::AdditiveDEAModel

    # Check parameters
    nx, m = size(X, 1), size(X, 2)
    ny, s = size(Y, 1), size(Y, 2)

    if Xref === nothing Xref = X end
    if Yref === nothing Yref = Y end

    nrefx, mref = size(Xref, 1), size(Xref, 2)
    nrefy, sref = size(Yref, 1), size(Yref, 2)

    if nx != ny
        throw(DimensionMismatch("number of rows in X and Y ($nx, $ny) are not equal"));
    end
    if nrefx != nrefy
        throw(DimensionMismatch("number of rows in Xref and Yref ($nrefx, $nrefy) are not equal"));
    end
    if m != mref
        throw(DimensionMismatch("number of columns in X and Xref ($m, $mref) are not equal"));
    end
    if s != sref
        throw(DimensionMismatch("number of columns in Y and Yref ($s, $sref) are not equal"));
    end

    if disposX != :Strong && disposX != :Weak
        throw(ArgumentError("`disposX` must be :Strong or :Weak"));
    end
    if disposY != :Strong && disposY != :Weak
        throw(ArgumentError("`disposY` must be :Strong or :Weak"));
    end

    if orient == :Input && disposX == :Weak
        throw(ArgumentError("`disposX` must be :Strong if `orient` is :Input"));
    end
    if orient == :Output && disposY == :Weak
        throw(ArgumentError("`disposY` must be :Strong if `orient` is :Output"));
    end
    if orient == :Graph && (disposX == :Weak || disposY == :Weak)
        throw(ArgumentError("`disposX` and `disposY` must be :Strong if `orient` is :Graph"));
    end

    # Default behaviour
    if model == :Default
        # If no weights are specified use :Ones
        if rhoX === nothing && rhoY === nothing
            model = :Ones
        else
            model = :Custom
        end
    end

    # Get weights based on the selected model
    if model != :Custom
        # Display error if both model and weights are specified
        if rhoX !== nothing || rhoY !== nothing
            throw(ArgumentError("`rhoX` and `rhoY` not allowed if model != :Custom"));
        end

        # Get weights for selected model
        rhoX, rhoY = deaaddweights(X, Y, model, orient = orient)
    end

    if size(rhoX) != size(X)
        throw(DimensionMismatch("size of rhoX and X ($(size(rhoX)), $(size(X))) are not equal"));
    end
    if size(rhoY) != size(Y)
        throw(DimensionMismatch("size of rhoY and Y ($(size(rhoY)), $(size(Y))) are not equal"));
    end

    # Parameters for additional condition in BAM model
    minXref = minimum(X, dims = 1)
    maxYref = maximum(Y, dims = 1)

    # Default optimizer
    if optimizer === nothing
        optimizer = DEAOptimizer(:LP)
    end

    # Compute efficiency for each DMU
    n = nx
    nref = nrefx

    effi = zeros(n)
    slackX = zeros(n, m)
    slackY = zeros(n, s)
    lambdaeff = spzeros(n, nref)

    for i=1:n
        # Value of inputs and outputs to evaluate
        x0 = X[i,:]
        y0 = Y[i,:]

        # Value of weights to evaluate
        rhoX0 = rhoX[i,:]
        rhoY0 = rhoY[i,:]

        # Set weights to zero if Weak disposability
        if disposX == :Weak
            rhoX0 = zeros(1, n)
        end

        if disposY == :Weak
            rhoY0 = zeros(1, n)
        end

        # Create the optimization model
        deamodel = newdeamodel(optimizer)

        @variable(deamodel, sX[1:m] >= 0)
        @variable(deamodel, sY[1:s] >= 0)
        @variable(deamodel, lambda[1:nref] >= 0)

        if orient == :Graph
            @objective(deamodel, Max, sum(rhoX0[j] * sX[j] for j in 1:m) + sum(rhoY0[j] * sY[j] for j in 1:s) )
        elseif orient == :Input
            @objective(deamodel, Max, sum(rhoX0[j] * sX[j] for j in 1:m)  )
        elseif orient == :Output
            @objective(deamodel, Max, sum(rhoY0[j] * sY[j] for j in 1:s) )
        else
            throw(ArgumentError("`orient` must be :Input, :Output or :Graph"));
        end

        @constraint(deamodel, [j in 1:m], sum(Xref[t,j] * lambda[t] for t in 1:nref) == x0[j] - sX[j])
        @constraint(deamodel, [j in 1:s], sum(Yref[t,j] * lambda[t] for t in 1:nref) == y0[j] + sY[j])

        # Add return to scale constraints
        if rts == :CRS
            # Add constraints for BAM CRS model
            if model == :BAM
                @constraint(deamodel, [j in 1:m], sum(Xref[t,j] * lambda[t] for t in 1:nref) >= minXref[j])
                @constraint(deamodel, [j in 1:s], sum(Yref[t,j] * lambda[t] for t in 1:nref) <= maxYref[j])
            end
        elseif rts == :VRS
            @constraint(deamodel, sum(lambda) == 1)
        else
            throw(ArgumentError("`rts` must be :CRS or :VRS"));
        end

        # Fix values of slacks when weight are zero
        for j = 1:m
            if rhoX0[j] == 0
                fix(sX[j], 0, force = true)
            end
        end

        for j = 1:s
            if rhoY0[j] == 0
                fix(sY[j], 0, force = true)
            end
        end

        # Optimize and return results
        JuMP.optimize!(deamodel)

        effi[i]  = JuMP.objective_value(deamodel)
        lambdaeff[i,:] = JuMP.value.(lambda)

        slackX[i,:] = JuMP.value.(sX)
        slackY[i,:] = JuMP.value.(sY)

        # Check termination status
        if (termination_status(deamodel) != MOI.OPTIMAL) && (termination_status(deamodel) != MOI.LOCALLY_SOLVED)
            @warn ("DMU $i termination status: $(termination_status(deamodel)). Primal status: $(primal_status(deamodel)). Dual status: $(dual_status(deamodel))")
        end

    end

    # Get X and Y targets
    Xtarget = X - slackX
    Ytarget = Y + slackY

    return AdditiveDEAModel(n, m, s, model, orient, rts, disposX, disposY, names, effi, slackX, slackY, lambdaeff, Xtarget, Ytarget)

end

function Base.show(io::IO, x::AdditiveDEAModel)
    compact = get(io, :compact, false)

    n = nobs(x)
    m = ninputs(x)
    s = noutputs(x)
    disposX = x.disposX
    disposY = x.disposY
    dmunames = names(x)

    eff = efficiency(x)
    slackX = slacks(x, :X)
    slackY = slacks(x, :Y)

    if !compact
        print(io, "Weighted Additive DEA Model \n")
        print(io, "DMUs = ", n)
        print(io, "; Inputs = ", m)
        print(io, "; Outputs = ", s)
        print(io, "\n")
        print(io, "Orientation = ", string(x.orient))
        print(io, "; Returns to Scale = ", string(x.rts))
        print(io, "\n")
        print(io, "Weights = ", string(x.weights))
        print(io, "\n")

        if disposX == :Weak print(io, "Weak disposability of inputs \n") end
        if disposY == :Weak print(io, "Weak disposability of outputs \n") end

        show(io, CoefTable(hcat(eff, slackX, slackY), ["efficiency"; ["slackX$i" for i in 1:m ]; ["slackY$i" for i in 1:s ]], dmunames))
    end

end

"""
    deaaddweights(X, Y, model)
Compute corresponding weights for related data envelopment analysis weighted additive models for inputs `X` and outputs `Y`.

Model specification:
- `:Ones`: standard additive DEA model.
- `:MIP`: Measure of Inefficiency Proportions. (Charnes et al., 1987; Cooper et al., 1999)
- `:Normalized`: Normalized weighted additive DEA model. (Lovell and Pastor, 1995)
- `:RAM`: Range Adjusted Measure. (Cooper et al., 1999)
- `:BAM`: Bounded Adjusted Measure. (Cooper et al, 2011)

# Optional Arguments
- `orient=:Graph`: choose between graph oriented `:Graph`, input oriented `:Input`, or output oriented model `:Output`.

"""
function deaaddweights(X::Union{Matrix,Vector}, Y::Union{Matrix,Vector},
    model::Symbol; orient::Symbol = :Graph)

    # Check orientation
    if orient != :Graph && orient != :Input && orient != :Output
        throw(ArgumentError("`orient` must be :Input, :Output or :Graph"));
    end

    # Compute specific weights based on the model
    if model == :Ones
        # Standard Additive DEA model
        if orient == :Graph || orient == :Input
            rhoX = ones(size(X))
        end
        if orient == :Graph || orient == :Output
            rhoY = ones(size(Y))
        end

    elseif model == :MIP
        # Measure of Inefficiency Proportions
        if orient == :Graph || orient == :Input
            rhoX = 1 ./ X
        end
        if orient == :Graph || orient == :Output
            rhoY = 1 ./ Y
        end

    elseif model == :Normalized
        # Normalized weighted additive DEA model
        if orient == :Graph || orient == :Input

            rhoX = zeros(size(X))
            m = size(X, 2)

            for i=1:m
                rhoX[:,i] .= 1 ./ std(X[:,i])
            end

            rhoX[isinf.(rhoX)] .= 0
        end
        if orient == :Graph || orient == :Output

            rhoY = zeros(size(Y))
            s = size(Y, 2)

            for i=1:s
                rhoY[:,i] .= 1 ./ std(Y[:,i])
            end

            rhoY[isinf.(rhoY)] .= 0
        end

    elseif model == :RAM
        # Range Adjusted Measure
        m = size(X, 2)
        s = size(Y, 2)

        normalization = 0
        if orient == :Graph
            normalization = m + s
        elseif orient == :Input
            normalization = m
        elseif orient == :Output
            normalization = s
        end

        if orient == :Graph || orient == :Input

            rhoX = zeros(size(X))

            for i=1:m
                rhoX[:,i] .= 1 ./ (normalization * (maximum(X[:,i])  - minimum(X[:,i])))
            end

            rhoX[isinf.(rhoX)] .= 0
        end
        if orient == :Graph || orient == :Output

            rhoY = zeros(size(Y))

            for i=1:s
                rhoY[:,i] .= 1 ./ (normalization * (maximum(Y[:,i])  - minimum(Y[:,i])))
            end

            rhoY[isinf.(rhoY)] .= 0
        end

    elseif model == :BAM
        # Bounded Adjusted Measure
        m = size(X, 2)
        s = size(Y, 2)

        normalization = 0
        if orient == :Graph
            normalization = m + s
        elseif orient == :Input
            normalization = m
        elseif orient == :Output
            normalization = s
        end

        if orient == :Graph || orient == :Input

            rhoX = zeros(size(X))
            minX = zeros(m)

            for i=1:m
                minX[i] = minimum(X[:,i])
                rhoX[:,i] = 1 ./ (normalization .* (X[:,i] .- minX[i] ))
            end

            rhoX[isinf.(rhoX)] .= 0
        end
        if orient == :Graph || orient == :Output

            rhoY = zeros(size(Y))
            maxY = zeros(s)

            for i=1:s
                maxY[i] = maximum(Y[:,i])
                rhoY[:,i] = 1 ./ (normalization .* (maxY[i] .- Y[:,i]))
            end

            rhoY[isinf.(rhoY)] .= 0
        end

    else
        throw(ArgumentError("Invalid model"));
    end

    if orient == :Input
        rhoY = ones(size(Y))
    elseif orient == :Output
        rhoX = ones(size(X))
    end

    return rhoX, rhoY

end