text:
# test_functions.jl
# File containing functions that reduce the suze of the test files

@doc """
###FreeFormDeformation.getnWallFaces

Get the number of element faces that exist on a geometric face of a mesh
despite the name of the function

**Inputs**

* `mesh` : Pumi DG mesh
* `bc_nums` : Array of boundary condition numbers that define the surface over
              which the number of element faces needs to be computed

**Output**

* `nwall_faces` : Number of element faces(edges in 2D) on each geometric face
"""->

function getnWallFaces(mesh::AbstractDGMesh, bc_nums::AbstractArray{Int,1})

  # Prepare the wall coordinates array for mesh warping
  nwall_faces = zeros(Int,length(bc_nums))
  vtx_per_face = mesh.dim # only true for simplex elements
  for (idx, itr) in enumerate(bc_nums)
    start_index = mesh.bndry_offsets[itr]
    end_index = mesh.bndry_offsets[itr+1]
    idx_range = start_index:(end_index-1)
    bndry_facenums = sview(mesh.bndryfaces, idx_range) # faces on geometric edge i
    nwall_faces[idx] = length(bndry_facenums)
  end      # End for itr = 1:length(geomfaces)

  return nwall_faces
end

@doc """
###FreeFormDeformation.getUniqueVertexArray

Get a 2D array of coordinates of unique vertices for the LOCAL portion of the
entire mesh.

**Input**

* `mesh` : Pumi DG mesh

**Output**

* `volNodes` : 2D array of unique vertices that exist on an MPI rank.
               size = [3, n_local_vertices]

"""->

function getUniqueVertexArray(mesh::AbstractMesh{Tmsh}) where Tmsh

  volNodes = zeros(Tmsh, 3, mesh.numVert)
  for i = 1:mesh.numEl
    for j = 1:size(mesh.vert_coords,2)
      # Get the vertex numbering on the portion of mesh owned by the processor
      local_vertnum = mesh.element_vertnums[j,i]
      for k = 1:mesh.dim
        volNodes[k, local_vertnum] = mesh.vert_coords[k,j,i] # mesh.element_vertnums
      end
    end
  end

  return volNodes
end

@doc """
###FreeFormDeformation.getUniqueWallCoordsArray

Get a 2D array of coordinates of unique vertices for the local portion of a
geometric face.

NOTE: Use this ONLY to compute the unique wall vertices on a given MPI rank.
      This function does not account for non-uniqueness that arises from the
      same vertex existing on 2 MPI ranks

**Inputs**

* `mesh` : Pumi DG mesh
* `bc_nums` : Array of boundary condition numbers that define hte surface over
              which the number of element faces needs to be computed

**Outputs**

* `wallCoords` : 2D array of coordinates of wall vertices. size = [3, total_vertices]

"""->

function getUniqueWallCoordsArray(mesh::AbstractMesh{Tmsh},
                            bc_nums::AbstractArray{Int,1}) where Tmsh

  vtx_arr = getUniqueVertexArray(mesh) # Get Unique vertex Array
  wallCoords = getUniqueWallCoordsArray(mesh, vtx_arr, bc_nums)

  return wallCoords
end

function getUniqueWallCoordsArray(mesh::AbstractMesh{Tmsh},
                            vtx_arr::AbstractArray{Tmsh,2},
                            bc_nums::AbstractArray{Int,1}) where Tmsh

  nface_verts = getLocalNumFaceVerts_unique(mesh, bc_nums)
  wallCoords = zeros(Tmsh, 3, nface_verts)

  local_vertnum_history = Int[]
  ctr = 1
  for (idx, itr) in enumerate(bc_nums)
    start_index = mesh.bndry_offsets[itr]
    end_index = mesh.bndry_offsets[itr+1]
    idx_range = start_index:(end_index-1)
    bndry_facenums = sview(mesh.bndryfaces, idx_range) # faces on geometric edge i
    nfaces = length(bndry_facenums)
    for i = 1:nfaces
      bndry_i = bndry_facenums[i]
      # get the local index of the vertices
      face_vtx_arr = mesh.topo.face_verts[:,bndry_i.face]
      for j = 1:length(face_vtx_arr)
        local_vertnum = mesh.element_vertnums[face_vtx_arr[j],bndry_i.element]
        if findfirst(local_vertnum_history, local_vertnum) == 0
          for k = 1:mesh.dim
            wallCoords[k,ctr] = vtx_arr[k, local_vertnum]
          end
          ctr += 1
        end
        push!(local_vertnum_history, local_vertnum)
      end     # for j = 1:length(face_vtx_arr)
    end       # for i = 1:nfaces
  end # End for itr = 1:length(bc_nums)

  return wallCoords
end

@doc """
###FreeFormDeformation.getGlobalUniqueWallCorrdsArray

Get a 2D array of coordinates of unique wall vertices that is owned by the MPI
rank. By default PUMI Mesh does not have ownership considering the unstructured
nature of the mesh where a vertex exists on several ranks. A sense of ownership
is created where the lowest rank containing the vertex owns it. The Array is
still distributed across the MPI ranks however the shared vertices appear on
only the lowest rank.

**Input**

* `mesh` : Pumi DG mesh
* `bc_nums` : Array of geometric faces (edges in 2D) over which the number
                 of element faces needs to be computed

**Output**

* `wallCoords` : 2D array of wall coordinates owned by the rank.
                 size = [3, n_owned_vertices]

"""->

function getGlobalUniqueWallCorrdsArray(mesh::AbstractMesh{Tmsh},
                                  bc_nums::AbstractArray{Int,1}) where Tmsh

  # Objective is to identify which MPI ranks have the boundary faces that lie
  # on a geometric face

  comm = MPI.COMM_WORLD
  comm_world = MPI.MPI_COMM_WORLD
  comm_self = MPI.COMM_SELF
  my_rank = MPI.Comm_rank(comm)
  comm_size = MPI.Comm_size(comm)

  wallCoords = zeros(Tmsh, 3, 0)
  vertex_coordinate = zeros(Tmsh,3)
  vtx_arr = getUniqueVertexArray(mesh) # Get Unique vertex Array

  # Gather nface_verts from all ranks
  nlocal_face_verts = getLocalNumFaceVerts_unique(mesh, bc_nums)
  nlocal_face_verts_arr = MPI.Allgather(nlocal_face_verts, comm)
  ranks_for_bndry_faces = Int[]
  for i = 1:length(nlocal_face_verts_arr)
    if nlocal_face_verts_arr[i] > 0
      push!(ranks_for_bndry_faces, i-1)
    end # End if nface_verts_arr[i] > 0
  end
  # println("rank = $my_rank, ranks_for_bndry_faces = $(ranks_for_bndry_faces)")
  if nlocal_face_verts != 0
    local_vertnum_history = Int[]
    for (idx, itr) in enumerate(bc_nums)
      start_index = mesh.bndry_offsets[itr]
      end_index = mesh.bndry_offsets[itr+1]
      idx_range = start_index:(end_index-1)
      bndry_facenums = sview(mesh.bndryfaces, idx_range)
      nfaces = length(bndry_facenums)
      for i = 1:nfaces
        bndry_i = bndry_facenums[i]
        face_vtx_arr = mesh.topo.face_verts[:,bndry_i.face]
        for j = 1:length(face_vtx_arr)
          local_vertnum = mesh.element_vertnums[face_vtx_arr[j],bndry_i.element]
          if findfirst(local_vertnum_history, local_vertnum) == 0
            # Now Check if this local vertex exists on multiple ranks
            if haskey(mesh.vert_sharing.rev_mapping, local_vertnum)

              # Check which all MPI ranks have this vertex and its position in
              # that rank's mesh.vert_sharing.
              rank_arr = first(mesh.vert_sharing.rev_mapping[local_vertnum])
              localIdx_arr = last(mesh.vert_sharing.rev_mapping[local_vertnum])
              # The Idea is that the lowest rank which has boundary face on the
              # geometric face always has ownership rights
              # if my_rank < minimum(rank_arr)
              #   for k = 1:mesh.dim
              #     vertex_coordinate[k] = vtx_arr[k, local_vertnum]
              #   end # End for k = 1:mesh.dim
              #   # Append a colum of coordinates
              #   wallCoords = hcat(wallCoords, vertex_coordinate)
              # else
              intersect_arr = intersect(rank_arr, ranks_for_bndry_faces)
              if intersect_arr == []
                for k = 1:mesh.dim
                  vertex_coordinate[k] = vtx_arr[k, local_vertnum]
                end # End for k = 1:mesh.dim
                # Append a colum of coordinates
                wallCoords = hcat(wallCoords, vertex_coordinate)
                # println("my_rank = $my_rank, rank_arr = $rank_arr, ranks_for_bndry_faces = $(ranks_for_bndry_faces)")
              elseif my_rank < minimum(intersect_arr)
                for k = 1:mesh.dim
                  vertex_coordinate[k] = vtx_arr[k, local_vertnum]
                end # End for k = 1:mesh.dim
                # Append a colum of coordinates
                wallCoords = hcat(wallCoords, vertex_coordinate)
              end # End if intersect_arr == []
              # end
            else # The vertex exist only on one MPI rank

              for k = 1:mesh.dim
                vertex_coordinate[k] = vtx_arr[k, local_vertnum]
              end
              # Append a colum of coordinates
              wallCoords = hcat(wallCoords, vertex_coordinate)

            end   # End if haskey(mesh.vert_sharing.rev_mapping, local_vertnum)
          end     # End if findfirst(local_vertnum_history, local_vertnum) == 0
          push!(local_vertnum_history, local_vertnum)
        end     # for j = 1:length(face_vtx_arr)
      end       # for i = 1:nfaces
    end # End for itr = 1:length(bc_nums)
  end

  MPI.Barrier(comm)

  return wallCoords
end
@doc """
### NOT NEEDED!!! ###
Keep this function around for future reference only

"""->
function getNeighborRank(mesh::AbstractDGMesh, bndry_face)

  npeers = length(mesh.peer_parts)
  mpi_neighbor_local_elem_no = zero(Int)
  ctr = 1
  for i = 1:npeers
    shared_interface_arr = mesh.shared_interfaces[i]
    for j = 1:length(shared_interface_arr)
      elementL = shared_interface_arr[j].elementL
      faceL = shared_interface_arr[j].faceL
      if elementL == bndry_face.element # && faceL == bndry_face.face
        mpi_neighbor_rank = mesh.peer_parts[i]
        return mpi_neighbor_rank
      end # End if elementL == local_elem_no
    end # End for j = 1:length(shared_interface_arr)
  end # End for i = 1:npeers

end

@doc """
###FreeFormDeformation.getLocalNumFaceVerts_unique

Obtain the number of unique face vertices that exist on one MPI rank. This
function DOES NOT take into account non-uniquness of vertices that exist across
multiple MPI ranks.

**Inputs**

* `mesh` : Pumi DG mesh
* `bc_nums : Array of boundary condition numbers that define the surface over
             which the number of element faces needs to be computed

**Output** :

* `ctr` : Number of unique vertices on a particular MPI rank
"""->
function getLocalNumFaceVerts_unique(mesh::AbstractDGMesh, bc_nums::AbstractArray{Int,1})

  local_vertnum_history = Int[]
  ctr = 0
  for (idx, itr) in enumerate(bc_nums)
    start_index = mesh.bndry_offsets[itr]
    end_index = mesh.bndry_offsets[itr+1]
    idx_range = start_index:(end_index-1)
    bndry_facenums = sview(mesh.bndryfaces, idx_range) # faces on geometric edge i
    nfaces = length(bndry_facenums)
    for i = 1:nfaces
      bndry_i = bndry_facenums[i]
      # get the local index of the vertices
      face_vtx_arr = mesh.topo.face_verts[:,bndry_i.face]
      for j = 1:length(face_vtx_arr)
        local_vertnum = mesh.element_vertnums[face_vtx_arr[j],bndry_i.element]
        if findfirst(local_vertnum_history, local_vertnum) == 0
          ctr += 1
        end
        push!(local_vertnum_history, local_vertnum)
      end # End for j = 1:length(face_vtx_arr)
    end   # End for i = 1:nfaces
  end

  return ctr
end

@doc """
###FreeFormDeformation.getWallCoords

Get a 2D array of coordinates that exist on a wall element face. This function
does not take into account uniqueness across adjacent element faces that share
a common vertex or across MPI ranks that share a vertex.

**Inputs**

* `mesh` : Pumi DG mesh
* `bc_nums` : Array of boundary condition numbers that define the surface over
              which the number of element faces needs to be computed

**Output**

* `wallCoords` : 2D array of coordinates of vertices that exist on a wall.
                 size = [3, sum(nwall_faces) \* vtx_per_face]

"""->

function getWallCoords(mesh::PumiMeshDG2, bc_nums::AbstractArray{Int, 1})

  nwall_faces = getnWallFaces(mesh, bc_nums)

  # Populate wallCoords
  vtx_per_face = mesh.dim # only true for simplex elements
  nWallCoords = sum(nwall_faces)*vtx_per_face
  wallCoords = zeros(3, nWallCoords)
  ctr = 1 # Counter for wall coordinates
  for (idx, itr) in enumerate(bc_nums)
    start_index = mesh.bndry_offsets[itr]
    end_index = mesh.bndry_offsets[itr+1]
    idx_range = start_index:(end_index-1)
    bndry_facenums = sview(mesh.bndryfaces, idx_range) # faces on geometric edge i
    nfaces = length(bndry_facenums)
    for i = 1:nfaces
      bndry_i = bndry_facenums[i]
      # get the local index of the vertices
      vtx_arr = mesh.topo.face_verts[:,bndry_i.face]
      for j = 1:length(vtx_arr)
        wallCoords[1:2, ctr] = mesh.vert_coords[:,vtx_arr[j],bndry_i.element]
        ctr += 1
      end  # End for j = 1:length(vtx_arr)
    end    # End for i = 1:nfaces
  end

  return wallCoords
end # End function getWallCoords


text:
# Before running this, please make sure to activate and instantiate the environment
# corresponding to [this `Project.toml`](https://raw.githubusercontent.com/alan-turing-institute/DataScienceTutorials.jl/master/Project.toml) and [this `Manifest.toml`](https://raw.githubusercontent.com/alan-turing-institute/DataScienceTutorials.jl/master/Manifest.toml)
# so that you get an environment which matches the one used to generate the tutorials:
#
# ```julia
# cd("DataScienceTutorials") # cd to folder with the *.toml
# using Pkg; Pkg.activate("."); Pkg.instantiate()
# ```

using RDatasets
using MLJ

iris = dataset("datasets", "iris")

first(iris, 3) |> pretty

iris2 = coerce(iris, :PetalWidth => OrderedFactor)
first(iris2[:, [:PetalLength, :PetalWidth]], 1) |> pretty

y, X = unpack(iris, ==(:Species), colname -> true)
first(X, 1) |> pretty

y, X = unpack(iris, ==(:Species), !=(:PetalLength))
first(X, 1) |> pretty

X, y = @load_iris;

for m in models(matching(X, y))
    if m.prediction_type == :probabilistic
        println(rpad(m.name, 30), "($(m.package_name))")
    end
end

knc = @load KNeighborsClassifier

linreg = @load LinearRegressor pkg="GLM"

# This file was generated using Literate.jl, https://github.com/fredrikekre/Literate.jl

ext: jl


text:
# WoodburyMatrices.SymWoodbury does not work with Distributions.MvNormal, and
# PDMatsExtras.WoodburyPDMat does not support non-diagonal `D`, so we here generalize
# PDMatsExtras.WoodburyPDMat

"""
    WoodburyPDMat <: PDMats.AbstractPDMat

Lazily represents a real positive definite (PD) matrix as an update to a full-rank PD matrix.

    WoodburyPDMat(A, B, D)

Constructs the ``n \\times n`` PD matrix
```math
W = A + B D B^\\mathrm{T},
```
where ``A`` is an ``n \\times n`` full rank positive definite matrix, ``D`` is an
``m \\times m`` symmetric matrix, and ``B`` is an ``n \\times m`` matrix. Note that ``B``
and ``D`` must be chosen such that ``W`` is positive definite; this is only implicitly
checked.

Overloads for `WoodburyPDMat` make extensive use of the following decomposition.
Let ``L_A L_A^\\mathrm{T} = A`` be the Cholesky decomposition of ``A``, and
let ``Q R = L_A^{-1} B`` be a thin QR decomposition. Define ``C = I + RDR^\\mathrm{T}``,
with the Cholesky decomposition ``L_C L_C^\\mathrm{T} = C``. Then, ``W = T T^\\mathrm{T}``,
where
```math
T = L_A Q \\begin{pmatrix} L_C & 0 \\\\ 0 & I \\end{pmatrix}.
```

The positive definite requirement is equivalent to the requirement that both ``A`` and
``C`` are positive definite.

For a derivation of this decomposition for the special case of diagonal ``A``, see
appendix A of [^Zhang2021].

[^Zhang2021]: Lu Zhang, Bob Carpenter, Andrew Gelman, Aki Vehtari (2021).
              Pathfinder: Parallel quasi-Newton variational inference.
              arXiv: [2108.03782](https://arxiv.org/abs/2108.03782) [stat.ML]
"""
struct WoodburyPDMat{
    T<:Real,
    TA<:AbstractMatrix{T},
    TB<:AbstractMatrix{T},
    TD<:AbstractMatrix{T},
    TUA<:Union{Diagonal{T},UpperTriangular{T}},
    TQ<:AbstractMatrix{T},
    TUC<:Union{Diagonal{T},UpperTriangular{T}},
} <: PDMats.AbstractPDMat{T}
    A::TA
    B::TB
    D::TD
    UA::TUA
    Q::TQ
    UC::TUC
end

function WoodburyPDMat(
    A::AbstractMatrix{T}, B::AbstractMatrix{T}, D::AbstractMatrix{T}
) where {T}
    cholA = cholesky(A)
    UA = cholA.U
    Q, R = qr(UA' \ B)
    cholC = cholesky(Symmetric(muladd(R, D * R', I)))
    return WoodburyPDMat(A, B, D, UA, Q, cholC.U)
end
function WoodburyPDMat(A, B, D)
    T = Base.promote_eltype(A, B, D)
    return WoodburyPDMat(
        convert(AbstractMatrix{T}, A),
        convert(AbstractMatrix{T}, B),
        convert(AbstractMatrix{T}, D),
    )
end

Base.Matrix(W::WoodburyPDMat) = Matrix(Symmetric(muladd(W.B, W.D * W.B', W.A)))

function Base.AbstractMatrix{T}(W::WoodburyPDMat) where {T}
    return WoodburyPDMat(
        map(k -> convert(AbstractMatrix{T}, getfield(W, k)), fieldnames(typeof(W)))...
    )
end

Base.getindex(W::WoodburyPDMat, inds...) = getindex(Matrix(W), inds...)

Base.adjoint(W::WoodburyPDMat) = W

Base.transpose(W::WoodburyPDMat) = W

function Base.inv(W::WoodburyPDMat)
    invUA = inv(W.UA)
    Anew = invUA * invUA'
    Bnew = invUA * Matrix(W.Q)
    invUC = inv(W.UC)
    Dnew = muladd(invUC, invUC', -I)
    return WoodburyPDMat(Anew, Bnew, Dnew)
end

LinearAlgebra.det(W::WoodburyPDMat) = exp(logdet(W))
LinearAlgebra.logdet(W::WoodburyPDMat) = 2 * (logdet(W.UA) + logdet(W.UC))
function LinearAlgebra.logabsdet(W::WoodburyPDMat)
    l = logdet(W)
    return (l, one(l))
end

function LinearAlgebra.diag(W::WoodburyPDMat)
    D = W.D isa Diagonal ? W.D : Symmetric(W.D)
    return diag(W.A) + map(b -> dot(b, D, b), eachrow(W.B))
end

function LinearAlgebra.lmul!(W::WoodburyPDMat, x::StridedVecOrMat)
    UA = W.UA
    UC = W.UC
    Q = W.Q
    k = minimum(size(W.B))
    lmul!(UA, x)
    lmul!(Q', x)
    x1 = x isa AbstractVector ? view(x, 1:k) : view(x, 1:k, :)
    lmul!(UC, x1)
    lmul!(UC', x1)
    lmul!(Q, x)
    lmul!(UA', x)
    return x
end

function LinearAlgebra.mul!(y::AbstractVector, W::WoodburyPDMat, x::StridedVecOrMat)
    return lmul!(W, copyto!(y, x))
end
function LinearAlgebra.mul!(y::AbstractMatrix, W::WoodburyPDMat, x::StridedVecOrMat)
    return lmul!(W, copyto!(y, x))
end

function Base.:*(W::WoodburyPDMat, c::Real)
    c > 0 || return Matrix(W) * c
    return WoodburyPDMat(W.A * c, W.B * one(c), W.D * c)
end

PDMats.dim(W::WoodburyPDMat) = size(W.A, 1)

function PDMats.invquad(W::WoodburyPDMat, x::AbstractVector{T}) where {T}
    v = W.Q' * (W.UA' \ x)
    n, m = size(W.B)
    k = min(m, n)
    return @views sum(abs2, W.UC' \ v[1:k]) + sum(abs2, v[(k + 1):n])
end

function PDMats.invquad!(r::AbstractArray, W::WoodburyPDMat, x::StridedMatrix{T}) where {T}
    v = lmul!(W.Q', W.UA' \ x)
    k = minimum(size(W.B))
    @views ldiv!(W.UC', v[1:k, :])
    colwise_sumsq!(r, v)
    return r
end

function PDMats.unwhiten!(
    r::StridedVecOrMat{T}, W::WoodburyPDMat, x::StridedVecOrMat{T}
) where {T}
    k = minimum(size(W.B))
    copyto!(r, x)
    @views lmul!(W.UC', x isa AbstractVector ? r[1:k] : r[1:k, :])
    lmul!(W.Q, r)
    lmul!(W.UA', r)
    return r
end

# adapted from https://github.com/JuliaStats/PDMats.jl/blob/master/src/utils.jl
function colwise_sumsq!(r::AbstractArray, a::AbstractMatrix)
    n = length(r)
    @assert n == size(a, 2)
    @inbounds for j in 1:n
        v = zero(eltype(a))
        @simd for i in axes(a, 1)
            v += abs2(a[i, j])
        end
        r[j] = v
    end
    return r
end


text:
using LiquidDSP
using Base.Test
import DSP

const dsp = DSP

# (Th, Tx) = (Float32, Float32)
for (Th, Tx) in [(Float32, Float32), (Float32, Complex64), (Complex64, Complex64)]
    h    = dsp.digitalfilter(dsp.Lowpass(0.5), dsp.FIRWindow(dsp.hamming(21)))
    h    = convert(Vector{Th}, h)
    x    = rand(Tx, 101)
    ff_d = dsp.FIRFilter(h, 2//1)
    ff_l = FIRInterp(Tx, 2, h)

    y_d  = dsp.filt(ff_d, x)
    y_l  = execute(ff_l, x)
    @test_approx_eq y_l y_d
end

max_stars_repo_name: JayKickliter/LiquidDSP

max_issues_repo_name: JayKickliter/LiquidDSP

max_forks_repo_name: JayKickliter/LiquidDSP


text:
module TestCoalescedCUDAEx

using Adapt
using CUDA
using Folds.Testing: parse_tests, test_with_sequential
using FoldsCUDA
using Transducers

# Note: each line must ends with `; init = .*)`. See below.
rawdata = """
mapreduce(x -> x^2, +, 1:10; init = 0)
mapreduce(*, +, 1:10, 11:20; init = 0)
mapreduce(*, +, 1:10, 11:20, 21:30; init = 0)
sum(1:10; init = 0)
sum(1:2^30; init = 0)
sum(first, pairs(1:10); init = 0)
sum(last, pairs(1:10); init = 0)
sum(x^2 for x in 1:11; init = 0)
sum(x^2 for x in 1:11 if isodd(x); init = 0)
sum(y for x in 1:11 if isodd(x) for y in 1:x:x^2; init = 0)
"""

tests = parse_tests(rawdata, @__MODULE__)

lines = map(eachline(IOBuffer(rawdata))) do ln
    replace(ln, r"; init = .*?\)$" => ")")
end
tests = append!(tests, parse_tests(join(lines, "\n"), @__MODULE__))

upload(x) = adapt(CuArray, x)

test() = test_with_sequential(tests, [CoalescedCUDAEx()]; upload = upload)

end  # module


text:
abstract type Shape end

struct Rectangle <: Shape
    corner_upper_left  :: Tuple{Float64, Float64}
    corner_lower_right :: Tuple{Float64, Float64}
end
shape_get_corner_upper_left(s::Shape) = shape_get_corner_upper_left(s.parent)
shape_get_corner_upper_left(r::Rectangle) = r.corner_upper_left
shape_get_corner_lower_right(s::Shape) = shape_get_corner_lower_right(s.parent)
shape_get_corner_lower_right(r::Rectangle) = r.corner_lower_right

struct Circle <: Shape
    center      :: Tuple{Float64, Float64}
    radius      :: Float64
end
shape_get_center(s::Shape) = shape_get_center(s.parent)
shape_get_center(c::Circle) = c.center
shape_get_radius(s::Shape) = shape_get_radius(s.parent)
shape_get_radius(c::Circle) = c.radius

struct GrayShape{shapeType <: Shape} <: Shape
    parent      :: shapeType
    gray_value  :: Float64
end
shape_get_gray_value(s::Shape) = shape_get_gray_value(s.parent)
shape_get_gray_value(g::GrayShape) = g.gray_value

struct AlphaShape{shapeType <: Shape} <: Shape
    parent      :: shapeType
    alpha_value :: Float64
end

shape_get_alpha_value(s::Shape) = shape_get_alpha_value(s.parent)
shape_get_alpha_value(a::AlphaShape) = a.alpha_value

r1 = AlphaShape(GrayShape(Rectangle((0.0, 1.0), (1.0, 0.0)), 0.5), 0.9)
r2 = GrayShape(AlphaShape(Rectangle((0.0, 1.0), (1.0, 0.0)), 0.9), 0.5)

println("r1: gray_value: ", shape_get_gray_value(r1))
println("r2: gray_value: ", shape_get_gray_value(r2))

# vim:syn=julia:cc=79:fdm=marker:sw=4:ts=4:


text:
module GnuplotJulia

mutable struct Gnuplot
    io::Base.Process
    counter::Int
    padding::Int
    ext::String
    function Gnuplot()
        gp = new()
        gp.io = open(`gnuplot -p`, "w")
        print(gp.io, "set colorsequence podo;")
        gp.counter = 0 # File numbering. Useful when animating.
        gp.padding = 3 # File name padding, e.g., 009, 010.
        return gp
    end
end

# TODO: using variadic function causes overhead?
import Base: println
function println(gp::Gnuplot, x::T, y::T) where {T <:Real}
    println(gp.io, x, " ", y)
end
function println(gp::Gnuplot, x::T, y::T, z::T) where {T <:Real}
    println(gp.io, x, " ", y, " ", z)
end
function println(gp::Gnuplot, x::T, y::T, z::T, u::T) where {T <:Real}
    println(gp.io, x, " ", y, " ", z, " ", u)
end
function println(gp::Gnuplot,
                 x::T, y::T, z::T,
                 u::T, v::T) where {T <:Real}
    println(gp.io, x, " ", y, " ", z, " ", u, " ", v)
end
function println(gp::Gnuplot,
                 x::T, y::T, z::T,
                 u::T, v::T, w::T) where {T <:Real}
    println(gp.io, x, " ", y, " ", z, " ", u, " ", v, " ", w)
end

function tex(gp, x::Int, y::Int, opts::String = "")
    gp.ext = ".tex"
    println(gp.io, "
            reset;
            set terminal epslatex standalone size $(x)cm, $(y)cm\\
            color background 'white'\\
            header '\\usepackage[utf8]{inputenc} \\usepackage{amsmath}';
            set colorsequence podo;
            ", opts)
end

function png(gp, x::Int, y::Int, opts::String = "")
    gp.ext = ".png"
    println(gp.io, "
            reset;
            set terminal png enhanced nocrop font 'verdana, 8' size $x, $y;
            set colorsequence podo;
            ", opts)
end

macro plot(opts::String, ex::Expr)
    quote
        println(gp.io, "\$x$(gp.counter) << EOD")
        $(esc(ex))
        println(gp.io, "EOD")
        println(gp.io, "plot \$x$(gp.counter) ", $opts)
        println(gp.io, "e")
        gp.counter += 1
    end
end
macro replot(opts::String, ex::Expr)
    quote
        println(gp.io, "\$x$(gp.counter) << EOD")
        $(esc(ex))
        println(gp.io, "EOD")
        println(gp.io, "replot \$x$(gp.counter) ", $opts)
        println(gp.io, "e")
        gp.counter += 1
    end
end

#macro plot(opts::String, ex::Expr)
#    quote
#        println(gp.io, "set output 'tmp-",
#                lpad(gp.counter, gp.padding, '0'), gp.ext, "'")
#        println(gp.io, "\$tmp << EOD")
#        $(esc(ex))
#        println(gp.io, "EOD")
#        println(gp.io, "plot \$tmp ", opts)
#        println(gp.io, "e")
#        gp.counter += 1
#    end
#end
#
#macro plot!(opts::String, ex::Expr)
#    quote
#        println(gp.io, "set output 'tmp-",
#                lpad(gp.counter, gp.padding, '0'), gp.ext, "'")
#        println(gp.io, "replot ", opts)
#        println(gp.io, "$id << EOD")
#        $(esc(ex))
#        println(gp.io, "EOD")
#        println(gp.io, "e")
#        gp.counter += 1
#    end
#end

function tex2pdf()
    println(gp.io, "set output")
    run(pipeline(`latexmk`, stdout=devnull))
    for i = 0:gp.counter-1
        x::String = lpad(i, gp.padding, '0')
        run(pipeline(`dvips tmp-$x.dvi`, stdout=devnull))
        run(pipeline(`gs -dBATCH -dSAFER -dNOPAUSE
                     -sDEVICE=pdfwrite -sOutputFile=output-$x.pdf
                     tmp-$x.ps`, stdout=devnull))
    end
    run(`find -x . -name tmp-\* -delete`)
    #display("application/pdf", read("output.pdf"))
end

function tex2png()
    println(gp.io, "set output")
    run(pipeline(`latexmk`, stdout=devnull))
    for i = 0:gp.counter-1
        x::String = lpad(i, gp.padding, '0')
        run(pipeline(`dvips tmp-$x.dvi`, stdout=devnull))
        run(pipeline(`gs -dBATCH -dSAFER -dNOPAUSE
                     -sDEVICE=pngalpha -sOutputFile=output-$x.png
                     tmp-$x.ps`, stdout=devnull))
    end
    run(`find -x . -name tmp-\* -delete`)
    #display("image/png", read("output.png"))
end

function mpeg()
    run(pipeline(`ffmpeg -y -i output-%0$(gp.padding)d.png output.mpeg`,
                 stdout=devnull, stderr=devnull))
    #display("text/html", """
    #        <video autoplay controls>
    #        <source src="output.mpeg" type="video/mpeg">
    #        </video>
    #        """)
end

function cleanup()
    close(gp.io)
    wait(gp.io.closenotify)
    gp.counter = 0
end

gp = Gnuplot()
@plot "with circles" begin
    println(gp, 1,2,3)
end
@replot "with circles" begin
    println(gp, 1,2,2)
end

end #module

max_stars_repo_name: gomezcruz/GnuplotJulia

max_issues_repo_name: gomezcruz/GnuplotJulia

max_forks_repo_name: gomezcruz/GnuplotJulia


text:
# find expression of projections on cones and their derivatives here:
#   https://stanford.edu/~boyd/papers/pdf/cone_prog_refine.pdf

"""
    projection_on_set(::DefaultDistance, v::AbstractVector{T}, ::MOI.Zeros) where {T}

projection of vector `v` on zero cone i.e. K = {0}
"""
function projection_on_set(::DefaultDistance, v::AbstractVector{T}, ::MOI.Zeros) where {T}
    return FillArrays.Zeros{T}(size(v))
end

"""
    projection_on_set(::AbstractDistance, ::MOI.Reals, v::Array{T}) where {T}

projection of vector `v` on real cone i.e. K = R
"""
function projection_on_set(::DefaultDistance, v::AbstractVector{T}, ::MOI.Reals) where {T}
    return v
end

function projection_on_set(::DefaultDistance, v::T, set::MOI.EqualTo) where {T}
    return zero(T) + set.value
end

function projection_on_set(::DefaultDistance, v::T, set::MOI.LessThan) where {T}
    return min(v, MOI.constant(set))
end

function projection_on_set(::DefaultDistance, v::T, set::MOI.GreaterThan) where {T}
    return max(v, MOI.constant(set))
end

"""
    projection_on_set(::DefaultDistance, v::AbstractVector{T}, ::MOI.Nonnegatives) where {T}

projection of vector `v` on Nonnegative cone i.e. K = R^n+
"""
function projection_on_set(::DefaultDistance, v::AbstractVector{T}, ::MOI.Nonnegatives) where {T}
    return max.(v, zero(T))
end

"""
    projection_on_set(::DefaultDistance, v::AbstractVector{T}, ::MOI.Nonpositives) where {T}

projection of vector `v` on Nonpositive cone i.e. K = R^n-
"""
function projection_on_set(::DefaultDistance, v::AbstractVector{T}, ::MOI.Nonpositives) where {T}
    return min.(v, zero(T))
end

"""
    projection_on_set(::NormedEpigraphDistance{p}, v::AbstractVector{T}, ::MOI.SecondOrderCone) where {T}

projection of vector `v` on second order cone i.e. K = {(t, x) ∈ R+ × Rn |  ||x|| ≤ t }
"""
function projection_on_set(::NormedEpigraphDistance{p}, v::AbstractVector{T}, ::MOI.SecondOrderCone) where {p, T}
    t = v[1]
    x = v[2:length(v)]
    norm_x = LinearAlgebra.norm(x, p)
    if norm_x <= t
        return copy(v)
    elseif norm_x <= -t
        return zeros(T, size(v))
    end
    result = zeros(T, size(v))
    result[1] = one(T)
    result[2:length(v)] = x / norm_x
    result *= (norm_x + t) / 2
    return result
end

function projection_on_set(::DefaultDistance, v::AbstractVector{T}, cone::MOI.SecondOrderCone) where {T}
    return projection_on_set(NormedEpigraphDistance{2}(), v, cone)
end

"""
    projection_on_set(::DefaultDistance, v::AbstractVector{T}, ::MOI.PositiveSemidefiniteConeTriangle) where {T}

projection of vector `v` on positive semidefinite cone i.e. K = S^n⨥
"""
function projection_on_set(::DefaultDistance, v::AbstractVector{T}, ::MOI.PositiveSemidefiniteConeTriangle) where {T}
    dim = isqrt(2*length(v))
    X = unvec_symm(v, dim)
    λ, U = LinearAlgebra.eigen(X)
    D = LinearAlgebra.Diagonal(max.(λ, 0))
    return vec_symm(U * D * U')
end

"""
    unvec_symm(x, dim)

Returns a dim-by-dim symmetric matrix corresponding to `x`.

`x` is a vector of length dim*(dim + 1)/2, corresponding to a symmetric matrix
```
X = [ X11  X12 ... X1k
      X21  X22 ... X2k
      ...
      Xk1  Xk2 ... Xkk ],
```
where `vec(X) = (X11, X12, X22, X13, X23, X33, ..., Xkk)`.

### Note on inner products

Note that the scalar product for the symmetric matrix in its vectorized form is
the sum of the pairwise product of the diagonal entries plus twice the sum of
the pairwise product of the upper diagonal entries; see [p. 634, 1].
Therefore, this transformation breaks inner products:
```
dot(unvec_symm(x, dim), unvec_symm(y, dim)) != dot(x, y).
```

### References

[1] Boyd, S. and Vandenberghe, L.. *Convex optimization*. Cambridge university press, 2004.
"""
function unvec_symm(x, dim=isqrt(2length(x)))
    X = zeros(eltype(x), dim, dim)
    idx = 1
    for i in 1:dim
        for j in 1:i
            X[j,i] = X[i,j] = x[idx]
            idx += 1
        end
    end
    return X
end

"""
    vec_symm(X)

Returns a vectorized representation of a symmetric matrix `X`.
`vec(X) = (X11, X12, X22, X13, X23, X33, ..., Xkk)`

### Note on inner products

Note that the scalar product for the symmetric matrix in its vectorized form is
the sum of the pairwise product of the diagonal entries plus twice the sum of
the pairwise product of the upper diagonal entries; see [p. 634, 1].
Therefore, this transformation breaks inner products:
```
dot(vec_symm(X), vec_symm(Y)) != dot(X, Y).
```

### References

[1] Boyd, S. and Vandenberghe, L.. *Convex optimization*. Cambridge university press, 2004.

"""
function vec_symm(X)
    return X[LinearAlgebra.triu(trues(size(X)))]
end

"""
    projection_on_set(::DefaultDistance, v::AbstractVector{T}, ::MOI.ExponentialCone) where {T}

projection of vector `v` on closure of the exponential cone
i.e. `cl(Kexp) = {(x,y,z) | y e^(x/y) <= z, y>0 } U {(x,y,z)| x <= 0, y = 0, z >= 0}`.

References:
* [Proximal Algorithms, 6.3.4](https://web.stanford.edu/~boyd/papers/pdf/prox_algs.pdf)
by Neal Parikh and Stephen Boyd.
* [Projection, presolve in MOSEK: exponential, and power cones]
(https://docs.mosek.com/slides/2018/ismp2018/ismp-friberg.pdf) by Henrik Friberg
"""
function projection_on_set(::DefaultDistance, v::AbstractVector{T}, s::MOI.ExponentialCone) where {T}
    _check_dimension(v, s)

    if _in_exp_cone(v; dual=false)
        return v
    end
    if _in_exp_cone(-v; dual=true)
        # if in polar cone Ko = -K*
        return zeros(T, 3)
    end
    if v[1] <= 0 && v[2] <= 0
        return [v[1], 0, max(v[3],0)]
    end

    return _exp_cone_proj_case_4(v)
end

function _in_exp_cone(v::AbstractVector{T}; dual=false, tol=1e-8) where {T}
    if dual
        return (
            (isapprox(v[1], 0, atol=tol) && v[2] >= 0 && v[3] >= 0) ||
            (v[1] < 0 && v[1]*exp(v[2]/v[1]) + ℯ*v[3] >= tol)
        )
    else
        return (
            (v[1] <= 0 && isapprox(v[2], 0, atol=tol) && v[3] >= 0) ||
            (v[2] > 0 && v[2] * exp(v[1] / v[2]) - v[3] <= tol)
        )
    end
end

function _exp_cone_proj_case_4(v::AbstractVector{T}) where {T}
    # Ref: https://docs.mosek.com/slides/2018/ismp2018/ismp-friberg.pdf, p47-48
    # Thm: h(x) is smooth, strictly increasing, and changes sign on domain
    r, s, t = v[1], v[2], v[3]
    h(x) = (((x-1)*r + s) * exp(x) - (r - x*s)*exp(-x))/(x^2 - x + 1) - t

    # Note: won't both be Inf by case 3 of projection
    lb = r > 0 ? 1 - s/r : -Inf
    ub = s > 0 ? r/s : Inf

    # Deal with ±Inf bounds
    if isinf(lb)
        lb = min(ub-0.125, -0.125)
        for _ in 1:10
            h(lb) < 0 && break
            ub = lb
            lb *= 2
        end
    end
    if isinf(ub)
        ub = max(lb+0.125, 0.125)
        for _ in 1:10
            h(ub) > 0 && break
            lb = ub
            ub *= 2
        end
    end

    # Check bounds
    if !(h(lb) < 0 && h(ub) > 0)
        error("Failure to find bracketing interval for exp cone projection.")
    end

    x = _bisection(h, lb, ub)
    if x === nothing
        error("Failure in root-finding for exp cone projection with boundaries ($lb, $ub).")
    end

    return ((x - 1) * r + s)/(x^2 - x + 1) * [x, 1, exp(x)]
end

"""
    projection_on_set(::DefaultDistance, v::AbstractVector{T}, ::MOI.DualExponentialCone) where {T}

projection of vector `v` on the dual exponential cone
i.e. `Kexp^* = {(u,v,w) | u < 0, -u*exp(v/u) <= ew } U {(u,v,w)| u == 0, v >= 0, w >= 0}`.

References:
* [Proximal Algorithms, 6.3.4](https://web.stanford.edu/~boyd/papers/pdf/prox_algs.pdf)
by Neal Parikh and Stephen Boyd.
* [Projection, presolve in MOSEK: exponential, and power cones](https://docs.mosek.com/slides/2018/ismp2018/ismp-friberg.pdf)
by Henrik Friberg
"""
function projection_on_set(d::DefaultDistance, v::AbstractVector{T}, ::MOI.DualExponentialCone) where {T}
    return v + projection_on_set(d, -v, MOI.ExponentialCone())
end

"""
    projection_on_set(::DefaultDistance, v::AbstractVector{T}, sets::Array{<:MOI.AbstractSet})

Projection onto `sets`, a product of sets
"""
function projection_on_set(::DefaultDistance, v::AbstractVector{T}, sets::Array{<:MOI.AbstractSet}) where {T}
    length(v) == length(sets) || throw(DimensionMismatch("Mismatch between value and set"))
    return reduce(vcat, (projection_on_set(DefaultDistance(), v[i], sets[i]) for i in eachindex(sets)))
end

"""
    projection_gradient_on_set(::DefaultDistance, v::AbstractVector{T}, ::MOI.Zeros) where {T}

derivative of projection of vector `v` on zero cone i.e. K = {0}^n
"""
function projection_gradient_on_set(::DefaultDistance, v::AbstractVector{T}, ::MOI.Zeros) where {T}
    return FillArrays.Zeros(length(v), length(v))
end

"""
    projection_gradient_on_set(::DefaultDistance, v::AbstractVector{T}, ::MOI.Reals) where {T}

derivative of projection of vector `v` on real cone i.e. K = R^n
"""
function projection_gradient_on_set(::DefaultDistance, v::AbstractVector{T}, ::MOI.Reals) where {T}
    return LinearAlgebra.Diagonal(ones(length(v)))
end

"""
    projection_gradient_on_set(::DefaultDistance, v::T, ::MOI.EqualTo)
"""
function projection_gradient_on_set(::DefaultDistance, ::T, ::MOI.EqualTo) where {T}
    return zero(T)
end

function projection_gradient_on_set(::DefaultDistance, v::T, s::MOI.LessThan) where {T}
    return oneunit(T) * (v <= MOI.constant(s))
end

function projection_gradient_on_set(::DefaultDistance, v::T, s::MOI.GreaterThan) where {T}
    return oneunit(T) * (v >= MOI.constant(s))
end

"""
    projection_gradient_on_set(::DefaultDistance, v::AbstractVector{T}, ::MOI.Nonnegatives) where {T}

derivative of projection of vector `v` on Nonnegative cone i.e. K = R^n+
"""
function projection_gradient_on_set(::DefaultDistance, v::AbstractVector{T}, ::MOI.Nonnegatives) where {T}
    y = (sign.(v) .+ one(T))/2
    return LinearAlgebra.Diagonal(y)
end

"""
    projection_gradient_on_set(::DefaultDistance, v::AbstractVector{T}, ::MOI.Nonpositives) where {T}

derivative of projection of vector `v` on Nonpositives cone i.e. K = R^n-
"""
function projection_gradient_on_set(::DefaultDistance, v::AbstractVector{T}, ::MOI.Nonpositives) where {T}
    y = (-sign.(v) .+ one(T))/2
    return LinearAlgebra.Diagonal(y)
end

"""
    projection_gradient_on_set(::NormedEpigraphDistance{p}, v::AbstractVector{T}, ::MOI.SecondOrderCone) where {T}

derivative of projection of vector `v` on second order cone i.e. K = {(t, x) ∈ R+ × Rn |  ||x|| ≤ t }
"""
function projection_gradient_on_set(::NormedEpigraphDistance{p}, v::AbstractVector{T}, ::MOI.SecondOrderCone) where {p,T}
    n = length(v)
    t = v[1]
    x = v[2:n]
    norm_x = LinearAlgebra.norm(x, p)
    if norm_x <= t
        return Matrix{T}(LinearAlgebra.I,n,n)
    elseif norm_x <= -t
        return zeros(T, n, n)
    end
    result = [
        norm_x     x';
        x          (norm_x + t)*Matrix{T}(LinearAlgebra.I,n-1,n-1) - (t/(norm_x^2))*(x*x')
    ]
    result /= (2 * norm_x)
    return result
end

function projection_gradient_on_set(::DefaultDistance, v::AbstractVector{T}, cone::MOI.SecondOrderCone) where {T}
    return projection_gradient_on_set(NormedEpigraphDistance{2}(), v, cone)
end

"""
    projection_gradient_on_set(::DefaultDistance, v::AbstractVector{T}, cone::MOI.PositiveSemidefiniteConeTriangle) where {T}

derivative of projection of vector `v` on positive semidefinite cone i.e. K = S^n⨥
"""
function projection_gradient_on_set(::DefaultDistance, v::AbstractVector{T}, ::MOI.PositiveSemidefiniteConeTriangle) where {T}
    n = length(v)
    dim = isqrt(2n)
    X = unvec_symm(v, dim)
    λ, U = LinearAlgebra.eigen(X)
    Tp = promote_type(T, Float64)

    # if all the eigenvalues are >= 0
    if all(λi ≥ zero(λi) for λi in λ)
        return Matrix{Tp}(LinearAlgebra.I, n, n)
    end

    # k is the number of negative eigenvalues in X minus ONE
    k = count(λi < 1e-4 for λi in λ)

    y = zeros(Tp, n)
    D = zeros(Tp, n, n)

    for idx in 1:n
        # set eigenvector
        y[idx] = 1

        # defining matrix B
        X̃ = unvec_symm(y, dim)
        B = U' * X̃ * U

        for i in 1:size(B)[1] # do the hadamard product
            for j in 1:size(B)[2]
                if (i <= k && j <= k)
                    @inbounds B[i, j] = 0
                elseif (i > k && j <= k)
                    λpi = max(λ[i], zero(Tp))
                    λmj = -min(λ[j], zero(Tp))
                    @inbounds B[i, j] *= λpi / (λmj + λpi)
                elseif (i <= k && j > k)
                    λmi = -min(λ[i], zero(Tp))
                    λpj = max(λ[j], zero(Tp))
                    @inbounds B[i, j] *= λpj / (λmi + λpj)
                end
            end
        end
        @inbounds D[idx, :] = vec_symm(U * B * U')
        # reset eigenvector
        @inbounds y[idx] = 0
    end
    return D
end

"""
    projection_gradient_on_set(::DefaultDistance, v::AbstractVector{T}, ::MOI.ExponentialCone) where {T}

derivative of projection of vector `v` on closure of the exponential cone,
i.e. `cl(Kexp) = {(x,y,z) | y e^(x/y) <= z, y>0 } U {(x,y,z)| x <= 0, y = 0, z >= 0}`.

References:
* [Solution Refinement at Regular Points of Conic Problems](https://stanford.edu/~boyd/papers/cone_prog_refine.html)
by Enzo Busseti, Walaa M. Moursi, and Stephen Boyd
"""
function projection_gradient_on_set(::DefaultDistance, v::AbstractVector{T}, s::MOI.ExponentialCone) where {T}
    _check_dimension(v, s)

    if _in_exp_cone(v; dual=false)
        return Matrix{T}(I, 3, 3)
    end
    if _in_exp_cone(-v; dual=true)
        # if in polar cone Ko = -K*
        return zeros(T, 3, 3)
    end
    if v[1] <= 0 && v[2] <= 0
        return LinearAlgebra.diagm(0 => T[1, 0, v[3] >= 0])
    end

    z1, z2, z3 = _exp_cone_proj_case_4(v)
    nu = z3 - v[3]
    rs = z1/z2
    exp_rs = exp(rs)

    mat = inv([
        1+nu*exp_rs/z2     -nu*exp_rs*rs/z2       0     exp_rs;
        -nu*exp_rs*rs/z2   1+nu*exp_rs*rs^2/z2    0     (1-rs)*exp_rs;
        0                  0                      1     -1
        exp_rs             (1-rs)*exp_rs          -1    0
    ])
    return mat[1:3,1:3]
end

"""
    projection_gradient_on_set(::DefaultDistance, v::AbstractVector{T}, ::MOI.DualExponentialCone) where {T}

derivative of projection of vector `v` on the dual exponential cone,
i.e. `Kexp^* = {(u,v,w) | u < 0, -u*exp(v/u) <= ew } U {(u,v,w)| u == 0, v >= 0, w >= 0}`.

References:
* [Solution Refinement at Regular Points of Conic Problems]
(https://stanford.edu/~boyd/papers/cone_prog_refine.html)
by Enzo Busseti, Walaa M. Moursi, and Stephen Boyd
"""
function projection_gradient_on_set(::DefaultDistance, v::AbstractVector{T}, ::MOI.DualExponentialCone) where {T}
    # from Moreau decomposition: x = P_K(x) + P_-K*(x)
    return I - projection_gradient_on_set(DefaultDistance(), -v, MOI.ExponentialCone())
end

"""
    projection_gradient_on_set(::DefaultDistance, v::AbstractVector{T}, sets::Array{<:MOI.AbstractSet})

Derivative of the projection of vector `v` on product of `sets`
projection_gradient_on_set[i,j] = ∂projection_on_set[i] / ∂v[j] where `projection_on_set` denotes projection of `v` on `cone`

Find expression of projections on cones and their derivatives here: https://stanford.edu/~boyd/papers/pdf/cone_prog_refine.pdf
"""
function projection_gradient_on_set(::DefaultDistance, v::AbstractVector{T}, sets::Array{<:MOI.AbstractSet}) where {T}
    length(v) == length(sets) || throw(DimensionMismatch("Mismatch between value and set"))
    return BlockDiagonal([projection_gradient_on_set(DefaultDistance(), v[i], sets[i]) for i in eachindex(sets)])
end


text:
module Coconet

# Write your package code here.
include("tonnetz.jl")
include("batching.jl")
include("utils.jl")
include("training.jl")

export pitch_to_tonnetz, tonnetz_to_pitch,
        get_batch, mask,
        +, -, *, /,
        train_model
end

max_stars_repo_name: MacKenzieHnC/Coconet

max_issues_repo_name: MacKenzieHnC/Coconet

max_forks_repo_name: MacKenzieHnC/Coconet


text:
module Stats

import ..GreyNumbers: GreyNumber

function I(n::Int)::Array{GreyNumber{Float64}, 2}
    mat = zeros(GreyNumber, n, n)
    for i in 1:n
        mat[i, i] = GreyNumber{Float64}(1.0, 1.0)
    end
    return mat
end

function mean(x::Array{GreyNumber{T}, 1})::GreyNumber{T} where {T <: Real}
     return sum(x) / length(x)
end

function var(x::Array{GreyNumber{T}, 1})::GreyNumber{T} where {T <: Real}
    m = mean(x)
    xm = map(a -> a - m, x)
    return sum((xm .* xm) / (length(x) - 1))
end

function std(x::Array{GreyNumber{T}, 1})::GreyNumber{T} where {T <: Real}
    return var(x) |> sqrt
end

function geomean(x::Array{GreyNumber{T}, 1})::GreyNumber{T} where {T <: Real}
    return exp(sum(log.(x)) / length(x))
end

function cor(x::Array{GreyNumber{T}, 1}, y::Array{GreyNumber{T}, 1})::GreyNumber{T} where {T <: Real}
    mx = mean(x)
    my = mean(y)
    xmd= map(a -> a - mx, x)
    ymd =map(a -> a - my, y)
    return (sum(xmd .* ymd) / sqrt(var(x) * var(y))) / (length(x) - 1.0)
end

function cor(x::Array{GreyNumber{T}, 1})::GreyNumber{T} where {T <: Real}
    _, p = size(x)
    cormat = ones(p, p)
    for i in 1:p
        for j in i:p
            if i != j
                cr = cor(x[:,i], x[:,j])
                cormat[i, j] = cr
                cormat[j, i] = cr
            end
        end
    end
    return cormat
end

end # end of module Stats


text:
# High level API exports

export Devices, dB, gainrange

"""
    Devices()

Enumerates all detectable SDR devices on the system.
Indexing into the returned `Devices` object returns a list of
keywords used to create a `Device` struct.
"""
struct Devices
    kwargslist::KWArgsList
    Devices() = new(KWArgsList(SoapySDRDevice_enumerate()...))
end
Base.length(d::Devices) = length(d.kwargslist)

function Base.show(io::IO, d::Devices)
    if length(d) == 0
        println(io, "No devices available! Make sure a supported SDR module is included.")
    end
    for (i, dev) in enumerate(d.kwargslist)
        print(io, "[$i] ")
        join(io, dev, ", ")
        println(io)
    end
end

Base.getindex(d::Devices, i::Integer) = d.kwargslist[i]
Base.iterate(d::Devices, state=1) = state > length(d) ? nothing : (d[state], state+1)
Base.open(d::KWArgsListRef) = Device(SoapySDRDevice_make(ptr(d)))

####################################################################################################
#    Device
####################################################################################################

"""
    Device

A device is a collection of SDR channels, obtained from the `Devices()` list.

Fields:
- `info`
- `driver`
- `hardware`
- `tx`
- `rx`
- `sensors`
- `time_source`
- `time_sources`
- `frontendmapping_rx`
- `frontendmapping_tx`
"""
mutable struct Device
    ptr::Ptr{SoapySDRDevice}
end
SoapySDRDevice_unmake(d::Device) = SoapySDRDevice_unmake(d.ptr)
Base.cconvert(::Type{<:Ptr{SoapySDRDevice}}, d::Device) = d
Base.unsafe_convert(::Type{<:Ptr{SoapySDRDevice}}, d::Device) = d.ptr
Base.close(d::Device) = SoapySDRDevice_unmake(d)

function Base.show(io::IO, d::Device)
    println(io, "SoapySDR ", d.hardware, " device")
    println(io, "  driver: ", d.driver)
    println(io, "  number of TX channels:", length(d.tx))
    println(io, "  number of RX channels:", length(d.rx))
    println(io, "  sensors: ", d.sensors)
    println(io, "  time_source: ", d.time_source)
    println(io, "  time_sources:", d.time_sources)
    println(io, "  frontendmapping_rx: ", d.frontendmapping_rx)
    println(io, "  frontendmapping_tx: ", d.frontendmapping_tx)
end

function Base.getproperty(d::Device, s::Symbol)
    if s === :info
        OwnedKWArgs(SoapySDRDevice_getHardwareInfo(d))
    elseif s === :driver
        Symbol(unsafe_string(SoapySDRDevice_getDriverKey(d)))
    elseif s === :hardware
        Symbol(unsafe_string(SoapySDRDevice_getHardwareKey(d)))
    elseif s === :tx
        ChannelList(d, Tx)
    elseif s === :rx
        ChannelList(d, Rx)
    elseif s === :sensors
        ComponentList(SensorComponent, d)
    elseif s === :time_sources
        ComponentList(TimeSource, d)
    elseif s === :time_source
        TimeSource(Symbol(unsafe_string(SoapySDRDevice_getTimeSource(d.ptr))))
    elseif s === :frontendmapping_tx
        unsafe_string(SoapySDRDevice_getFrontendMapping(d, Tx))
    elseif s === :frontendmapping_rx
        unsafe_string(SoapySDRDevice_getFrontendMapping(d, Rx))
    else
        getfield(d, s)
    end
end

function Base.setproperty!(c::Device, s::Symbol, v)
    if s === :frontendmapping_tx
        SoapySDRDevice_setFrontendMapping(c.ptr, Tx, v)
    elseif s === :frontendmapping_rx
        SoapySDRDevice_setFrontendMapping(c.ptr, Rx, v)
    elseif s === :time_source
        SoapySDRDevice_setTimeSource(c.ptr, v)
    else
        return setfield!(c, s, v)
    end
end

function Base.propertynames(::Device)
    return (:ptr, :info, :driver, :hardware, :tx, :rx, :sensors, :time)
end

####################################################################################################
#    Channel
####################################################################################################
"""
    Channel

A channel on the given `Device`.

Note!!: A Channel can be created from a `Device` or extracted from a `ChannelList`. It should rarely
be necessary to create a Channel directly.

Has the following properties:
- `device::Device` - `Device` to which the `Channel` belongs
- `direction` - either `Tx` or `Rx`
- `idx` - channel index used by Soapy
- `info` - channel info consiting of `OwnedKWArgs`
- `antenna` - antenna name
- `gain_mode` - Automatic Gain control, `true`, `false`, or `missing`
- `gain_elements` - list of `GainElements` of the channel
- `gain` - effective gain, distributed amongst the `GainElements`
- `dc_offset_mode` - Automatic DC offset mode, `true`, `false` or `missing`
- `dc_offset` -  DC offset value
- `iq_balance_mode` -  Automatic IQ balance mode, `true`, `false` or `missing`
- `iq_balance` - IQ balance value
- `frequency_correction` - frequency correction value
- `sample_rate` - sample rate
- `bandwidth` - bandwidth
- `frequency` - center frequency
- `fullduplex` - full duplex mode with other (TX/RX) channels
- `native_stream_format` - native stream format
- `stream_formats` - supported stream formats (converted by Soapy)
- `fullscale` - full scale value
- `sensors` - sensor list

## Reading and writing to Components

gains, antennas, and sensors may consist of a chain or selectable subcomponets.
To set or read e.g. a sensors, one may use the following syntax:

dev = Devices()[1]
cr = dev.rx[1]

# read a sensor value
s1 = cr.sensors[1]
cr[s1]

# read and set the gain element
g1 = cr.gain_elements[1]
cr[g1]
cr[g1] = 4*u"dB"
"""
struct Channel
    device::Device
    direction::Direction
    idx::Int
end
Base.show(io::IO, c::Channel) =
    print(io, "Channel(", c.device.hardware, ", ", c.direction, ", ", c.idx, ")")

####################################################################################################
#    Unit Printing
####################################################################################################

# Express everything in (kHz, MHz, GHz)
function pick_freq_unit(val::Quantity)
    iszero(val.val) && return val
    abs(val) >= 1.0u"GHz" ? uconvert(u"GHz", val) :
    abs(val) >= 1.0u"MHz" ? uconvert(u"MHz", val) :
                            uconvert(u"kHz", val)
end

# Print to 3 digits of precision, no decimal point
print_3_digit(io::IO, val::Quantity) =  print(io, Base.Ryu.writeshortest(round(val.val, sigdigits=3),
    #= plus =# false,
    #= space =# false,
    #= hash =# false,
    ))

function print_unit(io::IO, val::Quantity)
    print_3_digit(io, val)
    print(io, " ", unit(val))
end

function print_unit_interval(io::IO, min, max)
    if unit(min) == unit(max) || iszero(min.val)
        print_3_digit(io, min)
        print(io, "..")
        print_3_digit(io, max)
        print(io, " ", unit(max))
    else
        print_unit(io, min)
        print(io, " .. ")
        print_unit(io, max)
    end
end

function print_unit_steprange(io::IO, min, max, step)
    print_unit(io, min)
    print(io, ":")
    print_unit(io, step)
    print(io, ":")
    print_unit(io, max)
end

print_unit_interval(io::IO, x::Interval{<:Any, Closed, Closed}) =
    print_unit_interval(io, minimum(x), maximum(x))

using Intervals: Closed
function print_hz_range(io::IO, x::Interval{<:Any, Closed, Closed})
    min, max = pick_freq_unit(minimum(x)), pick_freq_unit(maximum(x))
    print_unit_interval(io, min, max)
end

function print_hz_range(io::IO, x::AbstractRange)
    min, step, max = pick_freq_unit(first(x)), pick_freq_unit(Base.step(x)), pick_freq_unit(last(x))
    print_unit_steprange(io, min, max, step)
end

function Base.show(io::IO, ::MIME"text/plain", c::Channel)
    println(io, c.direction == Tx ? "TX" : "RX", " Channel #", c.idx + 1, " on ", c.device.hardware)
    if !get(io, :compact, false)
        println(io, "  antenna: ", c.antenna)
        println(io, "  antennas: ", c.antennas)
        print(io, "  bandwidth [ ")
            join(io, map(x->sprint(print_hz_range, x), bandwidth_ranges(c)), ", ")
            println(io, " ]: ", pick_freq_unit(c.bandwidth))
        print(io, "  frequency [ ")
                join(io, map(x->sprint(print_hz_range, x), frequency_ranges(c)), ", ")
                println(io, " ]: ", pick_freq_unit(c.frequency))
            for element in FrequencyComponentList(c)
                print(io, "    ", element, " [ ")
                join(io, map(x->sprint(print_hz_range, x), frequency_ranges(c, element)), ", ")
                println(io, " ]: ", pick_freq_unit(c[element]))
            end
        println(io, "  gain_mode (AGC=true/false/missing): ", c.gain_mode)
        println(io, "  gain: ", c.gain)
        println(io, "  gain_elements: ", c.gain_elements)
        println(io, "  fullduplex: ", c.fullduplex)
        println(io, "  stream_formats: ", c.stream_formats)
        println(io, "  native_stream_format: ", c.native_stream_format)
        print(io, "  sample_rate [ ", )
            join(io, map(x->sprint(print_hz_range, x), sample_rate_ranges(c)), ", ")
            println(io, " ]: ", pick_freq_unit(c.sample_rate))
        println(io, "  dc_offset_mode (true/false/missing): ", c.dc_offset_mode)
        println(io, "  dc_offset (if has dc_offset_mode): ", c.dc_offset)
        println(io, "  iq_balance_mode (true/false/missing): ", c.iq_balance_mode)
        println(io, "  iq_balance: ", c.iq_balance)
        println(io, "  frequency_correction: ", c.frequency_correction, " ppm")
    end
end

"""
    ChannelList

A grouping of channels on the Device.
Note: This should not be called directly, but rather through the Device.rx and Device.tx properties.
"""
struct ChannelList <: AbstractVector{Channel}
    device::Device
    direction::Direction
end

function Base.size(cl::ChannelList)
    (SoapySDRDevice_getNumChannels(cl.device, cl.direction),)
end

function Base.getindex(cl::ChannelList, i::Integer)
    checkbounds(cl, i)
    Channel(cl.device, cl.direction, i-1)
end

function Base.getproperty(c::Channel, s::Symbol)
    if s === :info
        return OwnedKWArgs(SoapySDRDevice_getChannelInfo(c.device.ptr, c.direction, c.idx))
    elseif s === :antenna
        return Antenna(Symbol(unsafe_string(SoapySDRDevice_getAntenna(c.device.ptr, c.direction, c.idx))))
    elseif s === :antennas
        return AntennaList(c)
    elseif s === :gain
        return SoapySDRDevice_getGain(c.device.ptr, c.direction, c.idx)*dB
    elseif s === :gain_elements
        return GainElementList(c)
    elseif s === :dc_offset_mode
        if !SoapySDRDevice_hasDCOffsetMode(c.device.ptr, c.direction, c.idx)
            return missing
        end
        return SoapySDRDevice_getDCOffsetMode(c.device.ptr, c.direction, c.idx)
    elseif s === :dc_offset
        if !SoapySDRDevice_hasDCOffset(c.device.ptr, c.direction, c.idx)
            return missing
        end
        return SoapySDRDevice_getDCOffset(c.device.ptr, c.direction, c.idx)
    elseif s === :iq_balance_mode
        if !SoapySDRDevice_hasIQBalanceMode(c.device.ptr, c.direction, c.idx)
            return missing
        end
        return SoapySDRDevice_getIQBalanceMode(c.device.ptr, c.direction, c.idx)
    elseif s === :iq_balance
        if !SoapySDRDevice_hasIQBalance(c.device.ptr, c.direction, c.idx)
            return missing
        end
        return SoapySDRDevice_getIQBalance(c.device.ptr, c.direction, c.idx)
    elseif s === :gain_mode
        if !SoapySDRDevice_hasGainMode(c.device.ptr, c.direction, c.idx)
            return missing
        end
        return SoapySDRDevice_getGainMode(c.device.ptr, c.direction, c.idx)
    elseif s === :frequency_correction
        if !SoapySDRDevice_hasFrequencyCorrection(c.device.ptr, c.direction, c.idx)
            return missing
        end
        # TODO: ppm unit?
        return SoapySDRDevice_getFrequencyCorrection(c.device.ptr, c.direction, c.idx)
    elseif s === :sample_rate
        return SoapySDRDevice_getSampleRate(c.device.ptr, c.direction, c.idx) * Hz
    elseif s === :bandwidth
        return SoapySDRDevice_getBandwidth(c.device.ptr, c.direction, c.idx) * Hz
    elseif s === :frequency
        return SoapySDRDevice_getFrequency(c.device.ptr, c.direction, c.idx) * Hz
    elseif s === :fullduplex
        return Bool(SoapySDRDevice_getFullDuplex(c.device.ptr, c.direction, c.idx))
    elseif s === :stream_formats
        slist = StringList(SoapySDRDevice_getStreamFormats(c.device.ptr, c.direction, c.idx)...)
        return map(_stream_map_soapy2jl, slist)
    elseif s === :native_stream_format
        fmt, _ = SoapySDRDevice_getNativeStreamFormat(c.device.ptr, c.direction, c.idx)
        return _stream_map_soapy2jl(unsafe_string(fmt))
    elseif s === :fullscale
        _, fullscale = SoapySDRDevice_getNativeStreamFormat(c.device.ptr, c.direction, c.idx)
        return fullscale
    else
        return getfield(c, s)
    end
end

function Base.propertynames(::SoapySDR.Channel)
    return (:device,
            :direction,
            :idx,
            :info,
            :antenna,
            :antennas,
            :gain_elements,
            :gain,
            :dc_offset_mode,
            :dc_offset,
            :iq_balance_mode,
            :iq_balance,
            :gain_mode,
            :frequency_correction,
            :sample_rate,
            :bandwidth,
            :frequency,
            :fullduplex,
            :native_stream_format,
            :stream_formats,
            :fullscale,
            :sensors)
end

function Base.setproperty!(c::Channel, s::Symbol, v)
    if s === :frequency
        if isa(v, Quantity)
            SoapySDRDevice_setFrequency(c.device.ptr, c.direction, c.idx, uconvert(u"Hz", v).val, C_NULL)
        elseif isa(v, FreqSpec)
            error("TODO")
        else
            error("Frequency must be specified as either a Quantity or a FreqSpec!")
        end
    elseif s === :bandwidth
        SoapySDRDevice_setBandwidth(c.device.ptr, c.direction, c.idx, uconvert(u"Hz", v).val)
    elseif s === :gain_mode
        SoapySDRDevice_setGainMode(c.device.ptr, c.direction, c.idx, v)
    elseif s === :gain
        SoapySDRDevice_setGain(c.device.ptr, c.direction, c.idx, uconvert(u"dB", v).val)
    elseif s === :sample_rate
        SoapySDRDevice_setSampleRate(c.device.ptr, c.direction, c.idx, uconvert(u"Hz", v).val)
    else
        return setfield!(c, s, v)
    end
end

####################################################################################################
#    Components
####################################################################################################

# Components is an internal mechaism to allow for dispatch and interface through the Julia API
# For example there may be several GainElements we list in a Channel. A Julian idiom for this is
# the set/getindex class of functions.

abstract type AbstractComponent; end
Base.print(io::IO, c::AbstractComponent) = print(io, c.name)
Base.convert(::Type{T}, s::Symbol) where {T <: AbstractComponent} = T(s)
Base.convert(::Type{T}, s::String) where {T <: AbstractComponent} = T(Symbol(s))
Base.convert(::Type{Cstring}, s::AbstractComponent) =  Cstring(unsafe_convert(Ptr{UInt8}, s.name))

for e in (:TimeSource, :GainElement, :Antenna, :FrequencyComponent, :SensorComponent)
    @eval begin
        struct $e <: AbstractComponent;
            name::Symbol
        end
    end
end

struct ComponentList{T<:AbstractComponent} <: AbstractVector{T}
    s::StringList
end
Base.size(l::ComponentList) = (length(l.s),)
Base.getindex(l::ComponentList{T}, i::Integer) where {T} = T(Symbol(l.s[i]))
function Base.show(io::IO, l::ComponentList{T}) where {T}
    print(io, T, "[")
    for c in l
        print(io, c, ",")
    end
    print(io, "]")
end

for (e, f) in zip((:Antenna, :GainElement, :FrequencyComponent, :SensorComponent),
               (:SoapySDRDevice_listAntennas, :SoapySDRDevice_listGains, :SoapySDRDevice_listFrequencies, :SoapySDRDevice_listSensors))
    @eval begin
        $f(channel::Channel) =
            $f(channel.device, channel.direction, channel.idx)
        $(Symbol(string(e, "List")))(channel::Channel) = ComponentList{$e}(StringList($f(channel)...))
    end
end

function ComponentList(::Type{T}, d::Device) where {T <: AbstractComponent}
    if T <: SensorComponent
        ComponentList{SensorComponent}(StringList(SoapySDRDevice_listSensors(d.ptr)...))
    elseif T <: TimeSource
        ComponentList{TimeSource}(StringList(SoapySDRDevice_listTimeSources(d.ptr)...))
    end
end

using Base: unsafe_convert
function Base.getindex(c::Channel, ge::GainElement)
    SoapySDRDevice_getGainElement(c.device, c.direction, c.idx, Cstring(unsafe_convert(Ptr{UInt8}, ge.name))) * dB
end

function Base.getindex(c::Channel, fe::FrequencyComponent)
    SoapySDRDevice_getFrequencyComponent(c.device, c.direction, c.idx, Cstring(unsafe_convert(Ptr{UInt8}, fe.name))) * Hz
end

function Base.getindex(c::Channel, se::SensorComponent)
    unsafe_string(SoapySDRDevice_readChannelSensor(c.device, c.direction, c.idx, Cstring(unsafe_convert(Ptr{UInt8}, se.name))))
end

function Base.getindex(d::Device, se::SensorComponent)
    unsafe_string(SoapySDRDevice_readSensor(d.ptr, Cstring(unsafe_convert(Ptr{UInt8}, se.name))))
end

function Base.setindex!(c::Channel, gain::typeof(1.0dB), ge::GainElement)
    SoapySDRDevice_setGainElement(c.device, c.direction, c.idx, Cstring(unsafe_convert(Ptr{UInt8}, ge.name)), gain.val)
    return gain
end

## GainElement

function _gainrange(soapyr::SoapySDRRange)
    if soapyr.step == 0.0
        # Represents an interval rather than a range
        return (soapyr.minimum*dB)..(soapyr.maximum*dB)
    end
    # TODO
    @warn "Step ranges are not supported for gain elements. Returning Interval instead. Step: $(soapyr.step*dB)"
    #return range(soapyr.minimum*dB; stop=soapyr.maximum*dB, step=soapyr.step*dB)
    return (soapyr.minimum*dB)..(soapyr.maximum*dB)
end

function gainrange(c::Channel)
    return _gainrange(SoapySDRDevice_getGainRange(c.device, c.direction, c.idx))
end

gainrange(c::Channel, ge::GainElement) =
    return _gainrange(SoapySDRDevice_getGainElementRange(c.device, c.direction, c.idx, Cstring(unsafe_convert(Ptr{UInt8}, ge.name))))

function _hzrange(soapyr::SoapySDRRange)
    if soapyr.step == 0.0
        # Represents an interval rather than a range
        return (soapyr.minimum*Hz)..(soapyr.maximum*Hz)
    end
    return range(soapyr.minimum*Hz; stop=soapyr.maximum*Hz, step=soapyr.step*Hz)
end

function bandwidth_ranges(c::Channel)
    (ptr, len) = SoapySDRDevice_getBandwidthRange(c.device.ptr, c.direction, c.idx)
    arr = map(_hzrange, unsafe_wrap(Array, Ptr{SoapySDRRange}(ptr), (len,)))
    SoapySDR_free(ptr)
    arr
end

function frequency_ranges(c::Channel)
    (ptr, len) = SoapySDRDevice_getFrequencyRange(c.device.ptr, c.direction, c.idx)
    arr = map(_hzrange, unsafe_wrap(Array, Ptr{SoapySDRRange}(ptr), (len,)))
    SoapySDR_free(ptr)
    arr
end

function frequency_ranges(c::Channel, fe::FrequencyComponent)
    (ptr, len) = SoapySDRDevice_getFrequencyRangeComponent(c.device.ptr, c.direction, c.idx, Cstring(unsafe_convert(Ptr{UInt8}, fe.name)))
    arr = map(_hzrange, unsafe_wrap(Array, Ptr{SoapySDRRange}(ptr), (len,)))
    SoapySDR_free(ptr)
    arr
end

function sample_rate_ranges(c::Channel)
    (ptr, len) = SoapySDRDevice_getSampleRateRange(c.device.ptr, c.direction, c.idx)
    arr = map(_hzrange, unsafe_wrap(Array, Ptr{SoapySDRRange}(ptr), (len,)))
    SoapySDR_free(ptr)
    arr
end

"""
    list_sample_rates(::Channel)

List the natively supported sample rates for a given channel.
"""
function list_sample_rates(c::Channel)
    (ptr, len) = SoapySDRDevice_listSampleRates(c.device.ptr, c.direction, c.idx)
    arr = unsafe_wrap(Array, Ptr{Float64}(ptr), (len,)) * Hz
    SoapySDR_free(ptr)
    arr
end

### Frequency Setting
struct FreqSpec{T}
    val::T
    kwargs::Dict{Any, String}
end

####################################################################################################
#    Stream
####################################################################################################

"""
    SoapySDR.Stream(channels)
    SoapySDR.Stream(::Type{T}, channels)

Constructs a `Stream{T}` where `T` is the stream type of the device. If unspecified,
the native format will be used.

Fields:
- nchannels - The
- mtu - The stream Maximum Transmission Unit

## Example
```
SoapySDR.Stream(Devices()[1].rx)
```
"""
mutable struct Stream{T}
    d::Device
    nchannels::Int
    ptr::Ptr{SoapySDRStream}
    function Stream{T}(d::Device, nchannels, ptr::Ptr{SoapySDRStream}) where {T}
        this = new{T}(d, Int(nchannels), ptr)
        finalizer(SoapySDRDevice_closeStream, this)
        return this
    end
end
Base.cconvert(::Type{<:Ptr{SoapySDRStream}}, s::Stream) = s
Base.unsafe_convert(::Type{<:Ptr{SoapySDRStream}}, s::Stream) = s.ptr
SoapySDRDevice_closeStream(s::Stream) = SoapySDRDevice_closeStream(s.d, s)
Base.close(s::Stream) = SoapySDRDevice_closeStream(s)

streamtype(::Stream{T}) where T = T

function Base.setproperty!(c::Stream, s::Symbol, v)
    return setfield!(c, s, v)
end

function Base.getproperty(stream::Stream, s::Symbol)
    if s === :mtu
        SoapySDRDevice_getStreamMTU(stream.d.ptr, stream.ptr)
    else
        return getfield(stream, s)
    end
end

function Base.show(io::IO, s::Stream)
    print(io, "Stream on ", s.d.hardware)
end

function Stream(format::Type, device::Device, direction::Direction; kwargs...)
    isempty(kwargs) || error("TODO")
    Stream{T}(device, 1, SoapySDRDevice_setupStream(device, direction, string(format), C_NULL, 0, C_NULL))
end

function Stream(format::Type, channels::AbstractVector{T}; kwargs...) where {T <: Channel}
    soapy_format = _stream_map_jl2soapy(format)
    isempty(kwargs) || error("TODO")
    isempty(channels) && error("Must specify at least one channel or use the device/direction constructor for automatic.")
    device = first(channels).device
    direction = first(channels).direction
    if !all(channels) do channel
                channel.device == device && channel.direction == direction
            end
        throw(ArgumentError("Channels must agree on device and direction"))
    end
    Stream{format}(device, length(channels), SoapySDRDevice_setupStream(device, direction, soapy_format, map(x->x.idx, channels), length(channels), C_NULL))
end

function Stream(channels::AbstractVector{T}; kwargs...) where {T <: Channel}
    native_format = promote_type(map(c -> c.native_stream_format, channels)...)
    if native_format <: AbstractComplexInteger
        @warn "$(string(native_format)) may be poorly supported, it is recommend to specify a different type with Stream(format::Type, channels)"
    end
    Stream(native_format, channels; kwargs...)
end

"""
    read!(s::SoapySDR.Stream{T}, buffer::NTuple{N, Vector{T}}; [timeout])

Read data from the device into the given buffer.
"""
function Base.read!(s::Stream{T}, buffer::NTuple{N, Vector{T}}; timeout=nothing) where {N, T}
    timeout === nothing && (timeout = 0.1u"s") # Default from SoapySDR upstream

    total_nread = 0
    to_read_ct = length(first(buffer)) # TODO assert all equal
    while total_nread < to_read_ct
        nread, flags, timens = SoapySDRDevice_readStream(s.d, s, Ref(map(b -> pointer(b, total_nread+1), buffer)), to_read_ct-total_nread, uconvert(u"μs", timeout).val)
        timens = timens * u"ns"
        total_nread += nread
        #@assert flags & SOAPY_SDR_MORE_FRAGMENTS == 0
    end
    buffer
end

"""
    read(s::SoapySDR.Stream{T}, nb::Integer; [timeout])

Read at most `nb` samples from s
"""
function Base.read(s::Stream{T}, n::Integer; timeout=nothing) where {T}
    bufs = ntuple(_->Vector{T}(undef, n), s.nchannels)
    read!(s, bufs; timeout)
    bufs
end

function activate!(s::Stream; flags = 0, timens = nothing, numElems = nothing)
    SoapySDRDevice_activateStream(s.d, s, flags, timens === nothing ? 0 : uconvert(u"ns", timens).val, numElems === nothing ? s.mtu : numElems)
    nothing
end

function deactivate!(s::Stream; flags = 0, timens = nothing)
    SoapySDRDevice_deactivateStream(s.d, s, flags, timens === nothing ? 0 : uconvert(u"ns", timens).val)
    nothing
end

"""
    write(s::SoapySDR.Stream{T}, buffer::NTuple{N, Vector{T}}; [timeout]) where {N, T}

Write data from the device into the given buffer.
"""
function Base.write(s::Stream{T}, buffer::NTuple{N, Vector{T}}; timeout = nothing) where {N, T}
    timeout === nothing && (timeout = 0.1u"s") # Default from SoapySDR upstream

    total_nwritten = 0
    to_write_ct = length(buffer[1])

    while total_nwritten < to_write_ct
        nelem, flags = SoapySDRDevice_writeStream(s.d, s, Ref(map(b -> pointer(b,total_nwritten+1), buffer)), to_write_ct-total_nwritten, 0, 0, uconvert(u"μs", timeout).val)
        total_nwritten += nelem
    end

    buffer
end

"""
    get_sensor_info(::Device, ::String)

Read the sensor extracted from `list_sensors`.
Returns: the value as a string.
Note: Appropriate conversions need to be done by the user.
"""
function get_sensor_info(d::Device, name)
    SoapySDRDevice_getSensorInfo(d.ptr, name)
end

## Time API

"""
    has_hardware_time(::Device, what::String)

Query if the Device has hardware time for the given source.
"""
function has_hardware_time(d::Device, what::String)
    SoapySDRDevice_hasHardwareTime(d.ptr, what)
end

"""
    get_hardware_time(::Device, what::String)

Get hardware time for the given source.
"""
function get_hardware_time(d::Device, what::String)
    SoapySDRDevice_getHardwareTime(d.ptr, what)
end

"""
    set_hardware_time(::Device, timeNs::Int64 what::String)

Set hardware time for the given source.
"""
function set_hardware_time(d::Device, timeNs::Int64, what::String)
    SoapySDRDevice_setHardwareTime(d.ptr, timeNs, what)
end


text:
# Ambegaokar, V., & Troyer, M. (2010). Estimating errors reliably in Monte
# Carlo simulations of the Ehrenfest model. American Journal of Physics, 78(2),
# 150-157.

function stderr_bin(xs::AbstractVector)
    stderr_full = std(xs) / sqrt(length(xs))
    stderr_max = stderr_full

    while length(xs) >= 60
        xs = [0.5 * (xs[i] + xs[i+1]) for i in (1+length(xs)%2):2:length(xs)]
        stderr_max = max(stderr_max, std(xs) / sqrt(length(xs)))
    end

    autocorrelation_steps = 0.5 * ((stderr_max / stderr_full)^2 - 1)

    stderr_max, autocorrelation_steps
end

function aggregate(xs::AbstractVector)
    mean(xs), stderr_bin(xs)...
end