text:
function AbstractNFFTs.deconvolve!(p::NFFTPlan{T,D,R}, f::AbstractArray{U,D},
                      g::StridedArray{Complex{T},D}) where {T,D,R,U}
  dstart = Val(p.dims[1])
  dend = Val(p.dims[end])
  return _deconvolve!(p, f, g, dstart, dend)
end

@generated function _deconvolve!(p::NFFTPlan{T,D,R}, f::AbstractArray{U,D},
                                  g::StridedArray{Complex{T},D}, dstart::Val{DS},
                                  dend::Val{DE}) where {T,D,R,U,DS,DE}
    quote
        @nloops $D l f d->begin
            # preexpr
            if $DS ≤ d ≤ $DE
                offset_d = round(Int, p.n[d] - p.N[d]÷2) - 1
                gidx_d = rem(l_d+offset_d, p.n[d]) + 1
            else
                gidx_d = l_d
            end
        end begin
            # bodyexpr
            v = (@nref $D f l)
            @nexprs $(DE-DS+1) d -> v *= p.windowHatInvLUT[d][l_{d+$DS-1}]
            (@nref $D g gidx) = v
        end
    end
end

function AbstractNFFTs.deconvolve_transpose!(p::NFFTPlan{T,D,R},
                 g::AbstractArray{Complex{T},D}, f::StridedArray{U,D}) where {T,D,R,U}
  dstart = Val(p.dims[1])
  dend = Val(p.dims[end])
  return _deconvolve_transpose!(p, g, f, dstart, dend)
end

@generated function _deconvolve_transpose!(p::NFFTPlan{T,D,R},
           g::AbstractArray{Complex{T},D}, f::StridedArray{U,D}, dstart::Val{DS},
           dend::Val{DE}) where {T,D,R,U,DS,DE}
    quote
        @nloops $D l f d->begin
            # preexpr
            if $DS ≤ d ≤ $DE
                offset_d = round(Int, p.n[d] - p.N[d]÷2) - 1
                gidx_d = rem(l_d+offset_d, p.n[d]) + 1
            else
                gidx_d = l_d
            end
        end begin
            # bodyexpr
            v =  (@nref $D g gidx)
            @nexprs $(DE-DS+1) d -> v *= p.windowHatInvLUT[d][l_{d+$DS-1}]
            (@nref $D f l) = v
        end
    end
end

function AbstractNFFTs.convolve!(p::NFFTPlan{T,D,R}, g::AbstractArray{Complex{T},D},
      fHat::StridedArray{U,R}) where {T,D,R,U}
  l1 = Val(p.dims[1]-1)
  l2 = Val(length(p.dims))
  l3 = Val(D-p.dims[end])
  return _convolve!(p, g, fHat, l1, l2, l3)
end

@generated function _convolve!(p::NFFTPlan{T,D,R}, g::AbstractArray{Complex{T},D},
                   fHat::StridedArray{U,R}, l1::Val{L1}, l2::Val{L2}, l3::Val{L3}) where {T,D,R,U,L1,L2,L3}
  quote
    fill!(fHat, zero(T))
    scale = Int(p.params.LUTSize/(p.params.m))

    for k in 1:p.M
      ### The first L1 loops contain no NFFT
      @nloops_ $L1 d->l_{d} d->1:size(g,d) d->begin
         # preexpr
         gidx_d = l_d
         fidx_d = l_d
      end begin
          # bodyexpr
          @nexprs 1 d -> fidx_{$L1+1} = k
          @nexprs 1 d -> prodWin_{$L2} = one(T)
          @nexprs $L2 d -> xscale_d = p.x[d,k] * p.n[d+$L1]
          @nexprs $L2 d -> c_d = floor(Int, xscale_d)
          ### The next L2 loops contain the NFFT
          @nloops_ $L2 d->l_{d+$L1} d->(c_d-p.params.m+1):(c_d+p.params.m) d->begin
            # preexpr
            gidx_{d+$L1} = rem(l_{d+$L1}+p.n[d+$L1], p.n[d+$L1]) + 1
            idx = abs((xscale_d - l_{d+$L1} )*scale) + 1
            idxL = floor(idx)
            idxInt = Int(idxL)
            prodWin_{d-1} = prodWin_d * (p.windowLinInterp[idxInt] + ( idx-idxL ) * (p.windowLinInterp[idxInt+1] - p.windowLinInterp[idxInt]))
          end begin
             # bodyexpr
             ### The last L3 loops contain no NFFT
             @nloops_ $L3 d->l_{d+$L1+$L2} d->1:size(g,d+$L1+$L2) d->begin
                 # preexpr
                 gidx_{d+$L1+$L2} = l_{d+$L1+$L2}
                 fidx_{d+$L1+1} = l_{d+$L1+$L2}
              end begin
                (@nref $R fHat fidx) += (@nref $D g gidx) * prodWin_0
              end
          end
      end
    end
  end
end

function AbstractNFFTs.convolve_transpose!(p::NFFTPlan{T,D,R}, fHat::AbstractArray{U,R},
        g::StridedArray{Complex{T},D}) where {T,D,R,U}
  l1 = Val(p.dims[1]-1)
  l2 = Val(length(p.dims))
  l3 = Val(D-p.dims[end])
  return _convolve_transpose!(p, fHat, g, l1, l2, l3)
end

@generated function _convolve_transpose!(p::NFFTPlan{T,D,R}, fHat::AbstractArray{U,R},
            g::StridedArray{Complex{T},D}, l1::Val{L1}, l2::Val{L2}, l3::Val{L3}) where {T,D,R,U,L1,L2,L3}
  quote
    fill!(g, zero(T))
    scale = Int(p.params.LUTSize/(p.params.m))

    for k in 1:p.M
      ### The first L1 loops contain no NFFT
      @nloops_ $L1 d->l_{d} d->1:size(g,d) d->begin
         # preexpr
         gidx_d = l_d
         fidx_d = l_d
      end begin
          # bodyexpr
          @nexprs 1 d -> fidx_{$L1+1} = k
          @nexprs 1 d -> prodWin_{$L2} = one(T)
          @nexprs $L2 d -> xscale_d = p.x[d,k] * p.n[d+$L1]
          @nexprs $L2 d -> c_d = floor(Int, xscale_d)
          ### The next L2 loops contain the NFFT
          @nloops_ $L2 d->l_{d+$L1} d->(c_d-p.params.m+1):(c_d+p.params.m) d->begin
            # preexpr
            gidx_{d+$L1} = rem(l_{d+$L1}+p.n[d+$L1], p.n[d+$L1]) + 1
            idx = abs((xscale_d - l_{d+$L1})*scale) + 1
            idxL = floor(idx)
            idxInt = Int(idxL)
            prodWin_{d-1} = prodWin_d * (p.windowLinInterp[idxInt] + ( idx-idxL ) * (p.windowLinInterp[idxInt+1] - p.windowLinInterp[idxInt]))
          end begin
             # bodyexpr
             ### The last L3 loops contain no NFFT
             @nloops_ $L3 d->l_{d+$L1+$L2} d->1:size(g,d+$L1+$L2) d->begin
                 # preexpr
                 gidx_{d+$L1+$L2} = l_{d+$L1+$L2}
                 fidx_{d+$L1+1} = l_{d+$L1+$L2}
              end begin
                (@nref $D g gidx) += (@nref $R fHat fidx) * prodWin_0
              end
          end
      end
    end
  end
end


text:
# v0.6.0

using DataFrames

df = DataFrame(
EmployeeName=["Tyler Bennett", "John Rappl", "George Woltman", "Adam Smith",
"Claire Buckman", "David McClellan", "Rich Holcomb", "Nathan Adams",
"Richard Potter", "David Motsinger", "Tim Sampair", "Kim Arlich", "Timothy Grove"],
EmployeeID = ["E10297", "E21437", "E00127", "E63535", "E39876", "E04242",
"E01234", "E41298", "E43128", "E27002", "E03033", "E10001", "E16398"],
Salary = [32000, 47000, 53500, 18000, 27800, 41500, 49500, 21900, 15900, 19250,
27000, 57000, 29900],
Department = ["D101", "D050", "D101", "D202", "D202", "D101", "D202", "D050",
"D101", "D202", "D101", "D190", "D190"])

# To get only values
function firstnby(n::Int, y::Array, by::Array)
    # Check that each value belong to one and one only class
    if length(y) != length(by); error("y and by must have the same length"); end

    # Initialize resulting dictionary
    rst = Dict{eltype(by), Array{eltype(y)}}()

    # For each class...
    for cl in unique(by)
        # ...select the values of that class...
        i = find(x -> x == cl, by)
        # ...sort them and store them in result...
        rst[cl] = sort(y[i]; rev=true)
        # ...if length is greater than n select only first n elements
        if length(i) > n
            rst[cl] = rst[cl][1:n]
        end
    end
    return rst
end

for (cl, val) in firstnby(3, Array(df[:Salary]), Array(df[:Department]))
    println("$cl => $val")
end

# To get the full row...
function firstnby(n::Int, df::DataFrame, y::Symbol, by::Symbol)
    rst = Dict{eltype(df[by]), DataFrame}()

    for cl in unique(df[by])
        i = find(x -> x == cl, df[by])
        rst[cl] = sort(df[i, :]; cols=order(y; rev=true))
        if length(i) > n
            rst[cl] = rst[cl][1:n, :]
        end
    end
    return rst
end

for (cl, data) in firstnby(3, df, :Salary, :Department)
    println("\n$cl:\n$data")
end

ext: jl

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
mutable struct Rotational{T,N} <: Joint{T,N}
    V3::Adjoint{T,SVector{3,T}} # in body1's frame
    V12::SMatrix{2,3,T,6} # in body1's frame
    qoffset::UnitQuaternion{T} # in body1's frame

    Fτ::SVector{3,T}

    childid::Int64

    function Rotational{T,N}(body1::AbstractBody, body2::AbstractBody;
            axis::AbstractVector = zeros(3), qoffset::UnitQuaternion = one(UnitQuaternion)
        ) where {T,N}

        V1, V2, V3 = orthogonalrows(axis)
        V12 = [V1;V2]

        Fτ = zeros(T,3)

        childid = body2.id

        new{T,N}(V3, V12, qoffset, Fτ, childid), body1.id, body2.id
    end
end

Rotational0 = Rotational{T,0} where T
Rotational1 = Rotational{T,1} where T
Rotational2 = Rotational{T,2} where T
Rotational3 = Rotational{T,3} where T

function Base.show(io::IO, mime::MIME{Symbol("text/plain")}, constraint::Rotational{T,N}) where {T,N}
    summary(io, constraint)
    println(io,"")
    println(io, " V3:      "*string(constraint.V3))
    println(io, " V12:     "*string(constraint.V12))
    println(io, " qoffset: "*string(constraint.qoffset))
end

### Constraints and derivatives
## Position level constraints (for dynamics)
@inline function g(joint::Rotational, xa::AbstractVector, qa::UnitQuaternion, xb::AbstractVector, qb::UnitQuaternion)
    return Vmat(qa \ qb / joint.qoffset)
end
@inline function g(joint::Rotational, xb::AbstractVector, qb::UnitQuaternion)
    return Vmat(qb / joint.qoffset)
end

## Derivatives NOT accounting for quaternion specialness
@inline function ∂g∂posa(joint::Rotational{T}, xa::AbstractVector, qa::UnitQuaternion, xb::AbstractVector, qb::UnitQuaternion) where T
    X = szeros(T, 3, 3)
    Q = VRᵀmat(joint.qoffset) * Rmat(qb) * Tmat(T)

    return X, Q
end
@inline function ∂g∂posb(joint::Rotational{T}, xa::AbstractVector, qa::UnitQuaternion, xb::AbstractVector, qb::UnitQuaternion) where T
    X = szeros(T, 3, 3)
    Q = VRᵀmat(joint.qoffset) * Lᵀmat(qa)

    return X, Q
end
@inline function ∂g∂posb(joint::Rotational{T}, xb::AbstractVector, qb::UnitQuaternion) where T
    X = szeros(T, 3, 3)
    Q = VRᵀmat(joint.qoffset)

    return X, Q
end

## vec(G) Jacobian (also NOT accounting for quaternion specialness in the second derivative: ∂(∂ʳg∂posx)∂y)
@inline function ∂2g∂posaa(joint::Rotational{T}, xa::AbstractVector, qa::UnitQuaternion, xb::AbstractVector, qb::UnitQuaternion) where T
    XX = szeros(T, 9, 3)
    XQ = szeros(T, 9, 4)
    QX = szeros(T, 9, 3)
    QQ = kron(Vmat(T),VRᵀmat(joint.qoffset)*Rmat(qb)*Tmat(T))*∂L∂qsplit(T)

    return XX, XQ, QX, QQ
end
@inline function ∂2g∂posab(joint::Rotational{T}, xa::AbstractVector, qa::UnitQuaternion, xb::AbstractVector, qb::UnitQuaternion) where T
    XX = szeros(T, 9, 3)
    XQ = szeros(T, 9, 4)
    QX = szeros(T, 9, 3)
    QQ = kron(VLᵀmat(qa)*Tmat(T),VRᵀmat(joint.qoffset))*∂R∂qsplit(T)

    return XX, XQ, QX, QQ
end
@inline function ∂2g∂posba(joint::Rotational{T}, xa::AbstractVector, qa::UnitQuaternion, xb::AbstractVector, qb::UnitQuaternion) where T
    XX = szeros(T, 9, 3)
    XQ = szeros(T, 9, 4)
    QX = szeros(T, 9, 3)
    QQ = kron(VLᵀmat(qb),VRᵀmat(joint.qoffset))*∂Lᵀ∂qsplit(T)

    return XX, XQ, QX, QQ
end
@inline function ∂2g∂posbb(joint::Rotational{T}, xa::AbstractVector, qa::UnitQuaternion, xb::AbstractVector, qb::UnitQuaternion) where T
    XX = szeros(T, 9, 3)
    XQ = szeros(T, 9, 4)
    QX = szeros(T, 9, 3)
    QQ = kron(Vmat(T),VRᵀmat(joint.qoffset)*Lᵀmat(qa))*∂L∂qsplit(T)

    return XX, XQ, QX, QQ
end
@inline function ∂2g∂posbb(joint::Rotational{T}, xb::AbstractVector, qb::UnitQuaternion) where T
    XX = szeros(T, 9, 3)
    XQ = szeros(T, 9, 4)
    QX = szeros(T, 9, 3)
    QQ = kron(Vmat(T),VRᵀmat(joint.qoffset))*∂L∂qsplit(T)

    return XX, XQ, QX, QQ
end

### Forcing
## Application of joint forces (for dynamics)
@inline function applyFτ!(joint::Rotational{T}, statea::State, stateb::State, clear::Bool) where T
    τ = joint.Fτ
    xa, qa = posargsk(statea)
    xb, qb = posargsk(stateb)

    τa = vrotate(-τ, qa) # in world coordinates
    τb = -τa # in world coordinates

    τa = vrotate(τa,inv(qa)) # in local coordinates
    τb = vrotate(τb,inv(qb)) # in local coordinates

    statea.τk[end] += τa
    stateb.τk[end] += τb
    clear && (joint.Fτ = szeros(T,3))
    return
end
@inline function applyFτ!(joint::Rotational{T}, stateb::State, clear::Bool) where T
    τ = joint.Fτ
    xb, qb = posargsk(stateb)

    τa = -τ # in world coordinates
    τb = -τa # in world coordinates

    τb = vrotate(τb,inv(qb)) # in local coordinates

    stateb.τk[end] += τb
    clear && (joint.Fτ = szeros(T,3))
    return
end

## Forcing derivatives (for linearization)
# Control derivatives
@inline function ∂Fτ∂ua(joint::Rotational{T}, statea::State, stateb::State) where T
    BFa = (szeros(T, 3, 3))
    Bτa = -I

    return [BFa; Bτa]
end
@inline function ∂Fτ∂ub(joint::Rotational{T}, statea::State, stateb::State) where T
    xa, qa = posargsk(statea)
    xb, qb = posargsk(stateb)
    qbinvqa = qb\qa

    BFb = (szeros(T, 3, 3))
    Bτb = VLmat(qbinvqa) * RᵀVᵀmat(qbinvqa)

    return [BFb; Bτb]
end
@inline function ∂Fτ∂ub(joint::Rotational{T}, stateb::State) where T
    xb, qb = posargsk(stateb)

    BFb = (szeros(T, 3, 3))
    Bτb = VLᵀmat(qb) * RVᵀmat(qb)

    return [BFb; Bτb]
end

# Position derivatives
@inline function ∂Fτ∂posa(joint::Rotational{T}, statea::State, stateb::State) where T
    xa, qa = posargsk(statea)
    xb, qb = posargsk(stateb)
    τ = joint.Fτ

    FaXa = szeros(T,3,3)
    FaQa = szeros(T,3,3)
    τaXa = szeros(T,3,3)
    τaQa = szeros(T,3,3)
    FbXa = szeros(T,3,3)
    FbQa = szeros(T,3,3)
    τbXa = szeros(T,3,3)
    τbQa = 2*VLᵀmat(qb)*Rmat(qb)*Rᵀmat(qa)*Rmat(UnitQuaternion(τ))*LVᵀmat(qa)

    return FaXa, FaQa, τaXa, τaQa, FbXa, FbQa, τbXa, τbQa
end
@inline function ∂Fτ∂posb(joint::Rotational{T}, statea::State, stateb::State) where T
    xa, qa = posargsk(statea)
    xb, qb = posargsk(stateb)
    τ = joint.Fτ

    FaXb = szeros(T,3,3)
    FaQb = szeros(T,3,3)
    τaXb = szeros(T,3,3)
    τaQb = szeros(T,3,3)
    FbXb = szeros(T,3,3)
    FbQb = szeros(T,3,3)
    τbXb = szeros(T,3,3)
    τbQb = 2*VRᵀmat(qb)*Rᵀmat(qa)*Rmat(UnitQuaternion(τ))*Rmat(qa)*LVᵀmat(qb)

    return FaXb, FaQb, τaXb, τaQb, FbXb, FbQb, τbXb, τbQb
end
@inline function ∂Fτ∂posb(joint::Rotational{T}, stateb::State) where T
    xb, qb = posargsk(stateb)
    τ = joint.Fτ

    FaXb = szeros(T,3,3)
    FaQb = szeros(T,3,3)
    τaXb = szeros(T,3,3)
    τaQb = szeros(T,3,3)
    FbXb = szeros(T,3,3)
    FbQb = szeros(T,3,3)
    τbXb = szeros(T,3,3)
    τbQb = 2*VRᵀmat(qb)*Rmat(UnitQuaternion(τ))*LVᵀmat(qb)

    return FaXb, FaQb, τaXb, τaQb, FbXb, FbQb, τbXb, τbQb
end

### Minimal coordinates
## Position and velocity offsets
@inline function getPositionDelta(joint::Rotational, body1::AbstractBody, body2::Body, θ::SVector{N,T}) where {T,N}
    # axis angle representation
    θ = zerodimstaticadjoint(nullspacemat(joint)) * θ
    nθ = norm(θ)
    if nθ == 0
        q = one(UnitQuaternion{T})
    else
        q = UnitQuaternion(cos(nθ/2),(θ/nθ*sin(nθ/2))..., false)
    end

    Δq = q * joint.qoffset # in body1 frame
    return Δq
end
@inline function getVelocityDelta(joint::Rotational, body1::Body, body2::Body, ω::SVector)
    ω = zerodimstaticadjoint(nullspacemat(joint)) * ω
    Δω = vrotate(ω, inv(body2.state.qc)*body1.state.qc) # in body2 frame
    return Δω
end
@inline function getVelocityDelta(joint::Rotational, body1::Origin, body2::Body, ω::SVector)
    ω = zerodimstaticadjoint(nullspacemat(joint)) * ω
    Δω = vrotate(ω, inv(body2.state.qc)) # in body2 frame
    return Δω
end

## Minimal coordinate calculation
@inline function minimalCoordinates(joint::Rotational, body1::Body, body2::Body)
    statea = body1.state
    stateb = body2.state
    # q = g(joint, statea.xc, statea.qc, stateb.xc, stateb.qc)
    q = statea.qc \ stateb.qc / joint.qoffset
    return nullspacemat(joint) * rotation_vector(q)
end
@inline function minimalCoordinates(joint::Rotational, body1::Origin, body2::Body)
    stateb = body2.state
    # q = g(joint, stateb.xc, stateb.qc)
    q = stateb.qc / joint.qoffset
    return nullspacemat(joint) * rotation_vector(q)
end
@inline function minimalVelocities(joint::Rotational, body1::Body, body2::Body)
    statea = body1.state
    stateb = body2.state
    return nullspacemat(joint) * (vrotate(stateb.ωc,statea.qc\stateb.qc) - statea.ωc) # in body1's frame
end
@inline function minimalVelocities(joint::Rotational, body1::Origin, body2::Body)
    stateb = body2.state
    return nullspacemat(joint) * vrotate(stateb.ωc,stateb.qc) # in body1's frame
end


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

using Markdown
using InteractiveUtils

# ╔═╡ be0f3630-5028-11eb-0a61-a7a616623086
md"""
# Heaps : Ordered Trees

Elements are arranged in an increasing or decreasing order along the levels such that the root element is either the maximum (max heap) or minimum (min heap) element of the  tree. The parent must have higher

"""

# ╔═╡ c7e24be6-5029-11eb-1895-d5942b25dce0

# ╔═╡ Cell order:
# ╠═be0f3630-5028-11eb-0a61-a7a616623086
# ╠═c7e24be6-5029-11eb-1895-d5942b25dce0


text:
import Base: +, -, *, /
import Base.Broadcast: broadcasted

function +(x1::CuArray{T,N}, x2::CuArray{T,N}) where {T,N}
    @assert length(x1) == length(x2)
    y = fill!(similar(x1), 0)
    axpy!(T(1), x1, y)
    axpy!(T(1), x2, y)
    y
end

function -(x1::CuArray{T,N}, x2::CuArray{T,N}) where {T,N}
    @assert length(x1) == length(x2)
    y = fill!(similar(x1), 0)
    axpy!(T(1), x1, y)
    axpy!(T(-1), x2, y)
    y
end

function -(x::CuArray{T}) where T
    y = fill!(similar(x1), 0)
    axpy!(T(-1), x, y)
    y
end

*(A::CuMatrix{T}, x::CuVector{T}) where T = gemv('N', T(1), A, x)
*(A::CuMatrix{T}, B::CuMatrix{T}) where T = gemm('N', 'N', T(1), A, B)

function broadcasted(::typeof(+), x1::CuArray{T}, x2::CuArray{T}) where T
    x1, x2 = length(x1) >= length(x2) ? (x1,x2) : (x2,x1)
    y = copy(x1)
    CUDNN.add!(1, x2, 1, y)
    y
end

function broadcasted(::typeof(-), x1::CuArray{T}, x2::CuArray{T}) where T
    x1, x2 = length(x1) >= length(x2) ? (x1,x2) : (x2,x1)
    y = copy(x1)
    CUDNN.add!(-1, x2, 1, y)
    y
end

function broadcast_addto!(dest::CuArray{T}, src::CuArray{T}) where T
    CUDNN.add!(1, src, 1, dest)
end

function broadcasted(::typeof(*), x1::CuArray{T}, x2::CuArray{T}) where T
    elemtimes(x1, x2)
end

@generated function elemtimes(x1::CuArray{T,N}, x2::CuArray{T,N}) where {T,N}
    Ct = cstring(T)
    k = Kernel("""
    __global__ void elemtimes($Ct *y, Array<$Ct,$N> x1, Array<$Ct,$N> x2) {
        int idx = blockIdx.x * blockDim.x + threadIdx.x;
        if (idx >= x1.length()) return;
        int ndidxs[$N];
        x1.ndindex(ndidxs, idx);
        y[idx] = x1[idx] * x2(ndidxs);
    }
    """)
    quote
        length(x1) < length(x2) && return elemtimes(x2, x1)
        y = similar(x1)
        gdims, bdims = cudims(length(y))
        $k(gdims, bdims, pointer(y), x1, x2)
        y
    end
end


text:
abstract type AbstractAxis{T, BL, BR, I} <: AbstractVector{T} end

"""
    struct DiscreteAxis{T, BL, BR, I} <: AbstractAxis{T, BL, BR, I}

Axis with discrete ticks which is used to define a dimension of a [`Grid`](@ref).

## Parametric types
* `T`: Type of ticks
* `BL`: Boundary condition at the left endpoint.
* `BR`: Boundary condition at the right endpoint.
* `I`: IntervalSets.Interval (closed or open boundaries)

The boundary conditions of a `DiscreteAxis` can be
`BL, BR ∈ {:periodic, :reflecting, :infinite, :r0, :fixed}`.

## Fields
* `interval::I`: Interval that defines the range of the axis.
* `ticks::Vector{T}`: Array of values that correspond to the discrete ticks of the axis.

See also [`Grid`](@ref).
"""
struct DiscreteAxis{T, BL, BR, I} <: AbstractAxis{T, BL, BR, I}
    interval::I
    ticks::Vector{T}
end
@inline size(dx::DiscreteAxis{T, BL, BR}) where {T, BL, BR} = size(dx.ticks)
@inline IndexStyle(::Type{<:DiscreteAxis}) = IndexLinear()
@inline length(dx::DiscreteAxis{T, BL, BR}) where {T, BL, BR} = length(dx.ticks)
@inline getindex(dx::DiscreteAxis{T, BL, BR}, i::Int) where {T, BL, BR} = dx.ticks[i]
@inline setindex!(dx::DiscreteAxis{T, BL, BR}, v::T, i::Int) where {T, BL, BR} = setindex!(dx.ticks, v, i)
@inline axes(dx::DiscreteAxis{T, BL, BR}) where {T, BL, BR} = axes(dx.ticks)

function DiscreteAxis{T, BL, BR}(int::I, ticks::Vector{T})::DiscreteAxis{T, BL, BR, typeof(int)} where {T, BL, BR, I}
    return DiscreteAxis{T, BL, BR, typeof(int)}( int, ticks )
end

"""
    DiscreteAxis(left_endpoint::T, right_endpoint::T, BL::Symbol, BR::Symbol, L::Symbol, R::Symbol, ticks::AbstractVector{T}) where {T}

Constructor of a `DiscreteAxis`.

## Arguments
* `left_endpoint::T`: Left endpoint of the interval of the `DiscreteAxis`.
* `right_endpoint::T`: Right endpoint of the interval of the `DiscreteAxis`.
* `BL::Symbol`: Boundary condition at the left endpoint.
* `BR::Symbol`: Boundary condition at the right endpoint.
* `L::Symbol`: Boundary type of the left endpoint.
* `R::Symbol`: Boundary type of the right endpoint.
* `ticks::AbstractVector{T}`: Array of values that correspond to the discrete ticks of the axis.

The boundary conditions of a `DiscreteAxis` can be
`BL, BR ∈ {:periodic, :reflecting, :infinite, :r0, :fixed}`.

The boundary types of a `DiscreteAxis` can be `L, R ∈ {:closed, :open}`.

## Examples
    DiscreteAxis(-2.0, 2.0, :infinite, :infinite, :closed, :closed, collect(-2:0.1:2))
"""
function DiscreteAxis(left_endpoint::T, right_endpoint::T, BL::Symbol, BR::Symbol, L::Symbol, R::Symbol, ticks::AbstractVector{T}) where {T}
    int::Interval{L, R, T} = Interval{L, R, T}( left_endpoint, right_endpoint )
    return DiscreteAxis{T, BL, BR, typeof(int)}( int, ticks )
end

function sizeof(dx::DiscreteAxis{T, BL, BR}) where {T, BL, BR}
    return sizeof(dx.interval) + sizeof(dx.ticks)
end

function print(io::IO, dx::DiscreteAxis{T, BL, BR}) where {T, BL, BR}
    print(io, dx.interval, " - length = ", length(dx))
end
function println(io::IO, dx::DiscreteAxis{T, BL, BR}) where {T, BL, BR}
    println(io, dx.interval)
    println(io, "length = ", length(dx))
end
function show(io::IO, dx::DiscreteAxis{T, BL, BR}) where {T, BL, BR}
    print(io, dx)
end
function show(io::IO, ::MIME"text/plain", dx::DiscreteAxis{T, BL, BR}) where {T, BL, BR}
    show(io, dx)
end

function get_boundary_types(int::Interval{L, R})::Tuple{Symbol, Symbol} where {L, R}
    return L, R
end

function get_boundary_types(ax::DiscreteAxis{T,LB,RB})::NTuple{4, Symbol} where {T, LB, RB}
    return LB, RB, get_boundary_types(ax.interval)...
end

function get_extended_ticks( ax::DiscreteAxis{T, :reflecting, :reflecting} )::Vector{T} where {T}
    ticks_ext::Vector{T} = Array{T}(undef, length(ax.ticks) + 2)
    ticks_ext[2:end-1] = ax.ticks
    set_periodic_boundary_ticks!(ticks_ext, ax.interval)
    return ticks_ext
end
function get_extended_ticks( ax::DiscreteAxis{T, :fixed, :reflecting} )::Vector{T} where {T}
    ticks_ext::Vector{T} = Array{T}(undef, length(ax.ticks) + 2)
    ticks_ext[2:end-1] = ax.ticks
    set_periodic_boundary_ticks!(ticks_ext, ax.interval)
    return ticks_ext
end
function get_extended_ticks( ax::DiscreteAxis{T, :reflecting, :fixed} )::Vector{T} where {T}
    ticks_ext::Vector{T} = Array{T}(undef, length(ax.ticks) + 2)
    ticks_ext[2:end-1] = ax.ticks
    set_periodic_boundary_ticks!(ticks_ext, ax.interval)
    return ticks_ext
end
function get_extended_ticks( ax::DiscreteAxis{T, :periodic, :periodic} )::Vector{T} where {T}
    ticks_ext::Vector{T} = Array{T}(undef, length(ax.ticks) + 2)
    ticks_ext[2:end-1] = ax.ticks
    set_periodic_boundary_ticks!(ticks_ext, ax.interval)
    return ticks_ext
end
function get_extended_ticks( ax::DiscreteAxis{T, :infinite, :infinite} )::Vector{T} where {T}
    ticks_ext::Vector{T} = Array{T}(undef, length(ax.ticks) + 2)
    ticks_ext[2:end-1] = ax.ticks
    Δ::T = 1 * (ticks_ext[end-1] - ticks_ext[2])
    ticks_ext[1] = ticks_ext[2] - Δ
    ticks_ext[end] = ticks_ext[end - 1] + Δ
    return ticks_ext
end
function get_extended_ticks( ax::DiscreteAxis{T, :r0, :infinite} )::Vector{T} where {T}
    ticks_ext::Vector{T} = Array{T}(undef, length(ax.ticks) + 2)
    ticks_ext[2:end-1] = ax.ticks
    ticks_ext[1] = ticks_ext[2] - (ticks_ext[3] - ticks_ext[2])
    Δ::T = 1 * (ticks_ext[end-1] - ticks_ext[2])
    ticks_ext[end] = ticks_ext[end - 1] + Δ
    return ticks_ext
end
function get_extended_ticks( ax::DiscreteAxis{T, :r0, :fixed} )::Vector{T} where {T}
    ticks_ext::Vector{T} = Array{T}(undef, length(ax.ticks) + 2)
    ticks_ext[2:end-1] = ax.ticks
    ticks_ext[1] = ticks_ext[2] - (ticks_ext[3] - ticks_ext[2])
    Δ::T = 1 * (ticks_ext[end-1] - ticks_ext[2])
    ticks_ext[end] = ticks_ext[end - 1] + Δ
    return ticks_ext
end
function get_extended_ticks( ax::DiscreteAxis{T, :r0, :reflecting} )::Vector{T} where {T}
    ticks_ext::Vector{T} = Array{T}(undef, length(ax.ticks) + 2)
    ticks_ext[2:end-1] = ax.ticks
    ticks_ext[1] = ticks_ext[2] - (ticks_ext[3] - ticks_ext[2])
    Δ::T = ticks_ext[end-1] - ticks_ext[end - 2]
    ticks_ext[end] = ticks_ext[end - 1] + Δ
    return ticks_ext
end
function get_extended_ticks( ax::DiscreteAxis{T, :fixed, :fixed} )::Vector{T} where {T}
    # same as get_extended_ticks( ax::DiscreteAxis{T, :reflecting, :reflecting} )::Vector{T} where {T}
    ticks_ext::Vector{T} = Array{T}(undef, length(ax.ticks) + 2)
    ticks_ext[2:end-1] = ax.ticks
    set_periodic_boundary_ticks!(ticks_ext, ax.interval)
    return ticks_ext
end
function get_extended_ticks( ax::DiscreteAxis{T, :infinite, :fixed} )::Vector{T} where {T}
    ticks_ext::Vector{T} = Array{T}(undef, length(ax.ticks) + 2)
    ticks_ext[2:end-1] = ax.ticks
    set_periodic_boundary_ticks!(ticks_ext, ax.interval)
    Δ::T = 1 * (ticks_ext[end-1] - ticks_ext[2])
    ticks_ext[1] = ticks_ext[2] - Δ
    return ticks_ext
end
function get_extended_ticks( ax::DiscreteAxis{T, :infinite, :reflecting} )::Vector{T} where {T}
    ticks_ext::Vector{T} = Array{T}(undef, length(ax.ticks) + 2)
    ticks_ext[2:end-1] = ax.ticks
    set_periodic_boundary_ticks!(ticks_ext, ax.interval)
    Δ::T = 1 * (ticks_ext[end-1] - ticks_ext[2])
    ticks_ext[1] = ticks_ext[2] - Δ
    return ticks_ext
end
function get_extended_ticks( ax::DiscreteAxis{T, :fixed, :infinite} )::Vector{T} where {T}
    ticks_ext::Vector{T} = Array{T}(undef, length(ax.ticks) + 2)
    ticks_ext[2:end-1] = ax.ticks
    set_periodic_boundary_ticks!(ticks_ext, ax.interval)
    Δ::T = 1 * (ticks_ext[end-1] - ticks_ext[2])
    ticks_ext[end] = ticks_ext[end - 1] + Δ
    return ticks_ext
end
function get_extended_ticks( ax::DiscreteAxis{T, :reflecting, :infinite} )::Vector{T} where {T}
    ticks_ext::Vector{T} = Array{T}(undef, length(ax.ticks) + 2)
    ticks_ext[2:end-1] = ax.ticks
    set_periodic_boundary_ticks!(ticks_ext, ax.interval)
    Δ::T = 1 * (ticks_ext[end-1] - ticks_ext[2])
    ticks_ext[end] = ticks_ext[end - 1] + Δ
    return ticks_ext
end

function set_periodic_boundary_ticks!( ticks::Vector{T}, interval::Interval{:closed, :open, T})::Nothing where {T}
    ticks[1] = ticks[2] - (interval.right - ticks[end - 1])
    ticks[end] = interval.right
    nothing
end
function set_periodic_boundary_ticks!( ticks::Vector{T}, interval::Interval{:open, :closed, T})::Nothing where {T}
    ticks[1] = interval.left
    ticks[end] = ticks[end - 1] + (ticks[2] - interval.left)
    nothing
end
function set_periodic_boundary_ticks!( ticks::Vector{T}, interval::Interval{:open, :open, T})::Nothing where {T}
    ticks[1] = interval.left
    ticks[end] = interval.right
    nothing
end

function set_periodic_boundary_ticks!( ticks::Vector{T}, interval::Interval{:closed, :closed, T})::Nothing where {T, ispolaraxis}
    if length(ticks) == 3
        ticks[1] = ticks[2] - 2π
        ticks[end] = ticks[2] + 2π # -> Δmidpoint_φ = 2π -> area of circle is 2π * 0.5*r^2
    else
        ticks[1] = ticks[2] - (ticks[3] - ticks[2])
        ticks[end] = ticks[end - 1] + (ticks[end - 1] - ticks[end - 2])
    end
    nothing
end

function searchsortednearest(a::AbstractVector{T}, x::T)::Int where {T <: Real}
    idx::Int = searchsortedfirst(a, x)
    if (idx == 1) return idx end
    if (idx > length(a)) return length(a) end
    if (a[idx] == x) return idx end
    if (abs(a[idx] - x) < abs(a[idx - 1] - x))
        return idx
    else
        return idx - 1
    end
end
@inline function searchsortednearest(ax::DiscreteAxis{T}, x::T)::Int where {T <: Real}
    return searchsortednearest(ax.ticks, x)
end
function searchsortednearest(ax::DiscreteAxis{T, :periodic, :periodic}, x::T)::Int where {T <: Real}
    if x in ax.interval
        return searchsortednearest(ax.ticks, x)
    else
        period::T = width(ax.interval)
        v::T = x
        while v >= ax.interval.right
            v -= period
        end
        while v < ax.interval.left
            v += period
        end
        return searchsortednearest(ax.ticks, v)
    end

end

function DiscreteAxis(nt::NamedTuple; unit = u"m/m")
    T = typeof(ustrip(nt.knots[1]))
    knots::Vector{T} = convert(Vector{T}, ustrip.(uconvert.(unit, nt.knots)))
    lep::T = ustrip(uconvert.(unit, nt.interval.left_boundary.endpoint ))
    rep::T = ustrip(uconvert.(unit, nt.interval.right_boundary.endpoint))
    int = Interval{nt.interval.left_boundary.closedopen, nt.interval.right_boundary.closedopen}( lep, rep )
    return DiscreteAxis{T, nt.interval.left_boundary.boundaryhandling, nt.interval.right_boundary.boundaryhandling, typeof(int)}(
        int, knots
    )
end

Base.convert(T::Type{DiscreteAxis}, x::NamedTuple; unit = u"m/m") = T(x)

function NamedTuple(ax::DiscreteAxis{T, BL, BR}; unit = u"m/m") where {T, BL, BR}
    int::Interval = ax.interval
    int_types::Tuple{Symbol, Symbol} = get_boundary_types(int)
    return (
        knots = ax.ticks * unit,
        interval = (
            left_boundary = (
                endpoint = int.left * unit,
                closedopen = int_types[1],
                boundaryhandling = BL,
            ),
            right_boundary = (
                endpoint = int.right * unit,
                closedopen = int_types[2],
                boundaryhandling = BR,
            ),
        )
    )
end

Base.convert(T::Type{NamedTuple}, x::DiscreteAxis; unit = u"m/m") = T(x)

function merge_closest_ticks!(v::AbstractVector{T}, n::Int = length(v); min_diff::T = T(1e-6)) where {T}
    n == 1 && return n
    Δv = diff(v[1:n])
    Δ_min, Δv_min_indx = findmin(Δv)
    vFirst = v[1]
    vLast  = v[n]
    if Δ_min < min_diff
        v[Δv_min_indx] = (v[Δv_min_indx]+v[Δv_min_indx+1]) / 2
        v[Δv_min_indx+1:end-1] = v[Δv_min_indx+2:end]
        v[1] = vFirst
        n -= 1
        v[n] = vLast
        n
    else
        n
    end
end
function merge_close_ticks(v::AbstractVector{T}; min_diff::T = T(1e-6)) where {T}
    l = length(v)
    l <= 1 && return v
    n = merge_closest_ticks!(v, min_diff = min_diff)
    reduced = n < l
    l = n
    while reduced
        n = merge_closest_ticks!(v, n, min_diff = min_diff)
        reduced = n < l
        l = n
    end
    v[1:n]
end

"""
    merge_second_order_important_ticks(imp::Vector{T}, imp_second_order::Vector{T}; min_diff::T = T(1e-6)) where {T}

Merge all elements of the second vector, `imp_second_order`, into the first vector, `imp`,
if they are not too close (via `min_diff`) to elements of the first vector.
Returns the merged vector sorted.
"""
function merge_second_order_important_ticks(imp::Vector{T}, imp_second_order::Vector{T}; min_diff::T = T(1e-6)) where {T}
    sorted_imp = sort(imp)
    nearest_inds = map(x -> searchsortednearest(sorted_imp, x), imp_second_order)
    merge_inds = filter(i -> abs(sorted_imp[nearest_inds[i]] - imp_second_order[i]) >= min_diff, eachindex(imp_second_order))
    return sort!(vcat(imp, imp_second_order[merge_inds]))
end

function get_new_ticks_to_equalize_ratios_on_side(t::AbstractVector{T}; max_ratio = T(2)) where {T}
    @assert length(t) == 3
    Δt = diff(t)
    r = Δt[2] / Δt[1]
    min_ratio = inv(max_ratio)
    new_ticks = T[]
    if r > max_ratio
        x0 = t[2]
        Δnt = max_ratio * Δt[1]
        nt = x0 + Δnt # new tick
        while t[3] - nt > Δnt
            push!(new_ticks, nt)
            Δnt = max_ratio * Δnt
            nt = nt + Δnt # new tick
        end
        vcat(t[1:2], new_ticks, t[3:3])
    elseif r < min_ratio
        x0 = t[2]
        Δnt = max_ratio * Δt[2]
        nt = x0 - Δnt # new tick
        while nt - t[1] > Δnt
            push!(new_ticks, nt)
            Δnt = max_ratio * Δnt
            nt = nt - Δnt # new tick
        end
        sort!(new_ticks)
        vcat(t[1:1], new_ticks, t[2:3])
    else
        t
    end
end

function initialize_axis_ticks(t::AbstractVector{T}; max_ratio = T(2)) where {T}
    @assert max_ratio >= 1
    if length(t) <= 2 return t end
    # First: Left and right side intervals:
    new_ticks_left = get_new_ticks_to_equalize_ratios_on_side(t[1:3]; max_ratio)
    if length(t) == 3 return new_ticks_left end
    if length(t) == 4
        return vcat(new_ticks_left, t[4:4])
    end
    new_ticks_right = get_new_ticks_to_equalize_ratios_on_side(t[end-2:end]; max_ratio)
    if length(t) == 5
        return vcat(new_ticks_left, new_ticks_right[2:end])
    end
    ticks = vcat(new_ticks_left, t[4:end-3], new_ticks_right)
    # Second: Intervals in between
    iL = length(new_ticks_left) - 1
    iR = iL + 3
    n_sub_ints = length(t) - 5
    for i_sub_int in 1:n_sub_ints
        subticks = ticks[iL:iR]
        Δst = diff(subticks)
        r_l = Δst[2] / Δst[1]
        r_r = Δst[2] / Δst[3]
        r_lr = Δst[1] / Δst[3]
        new_ticks = T[]
        if r_lr < 1 # interval to the left is smaller than interval to the right
            if r_l > max_ratio
                x0 = subticks[2]
                Δnt = max_ratio * Δst[1]
                nt = x0 + Δnt # new tick
                while subticks[3] - nt > Δnt
                    push!(new_ticks, nt)
                    Δnt = max_ratio * Δnt
                    nt = nt + Δnt # new tick
                end
            end
        else # interval to the right is smaller than interval to the left
            if r_r > max_ratio
                x0 = subticks[3]
                Δnt = max_ratio * Δst[3]
                nt = x0 - Δnt # new tick
                while nt - subticks[2] > Δnt
                    push!(new_ticks, nt)
                    Δnt = max_ratio * Δnt
                    nt = nt - Δnt # new tick
                end
                sort!(new_ticks)
            end
        end

        ticks = vcat(ticks[1:iL+1], new_ticks, ticks[iR-1:end])

        iL += 1 + length(new_ticks)
        iR += 1 + length(new_ticks)
    end
    @assert issorted(ticks)
    @assert allunique(ticks)
    ticks
end

function fill_up_ticks(v::AbstractVector{T}, max_diff::T) where {T}
    if length(v) == 1 return v end
    Δv = diff(v)
    add_n_points = Int.(round.(Δv ./ max_diff, RoundUp)) .- 1
    r = Vector{T}(undef, length(v) + sum(add_n_points))
    r[:] .= 0
    i = 0
    for j in eachindex(Δv)
        x0 = v[j]
        n = add_n_points[j]
        Δ = Δv[j] / (n + 1)
        for l in 0:n
            i += 1
            r[i] = x0 + Δ * l
        end
    end
    r[end] = v[end]
    r
end

function even_tick_axis(ax::DiscreteAxis)
    if isodd(length(ax))
        int = ax.interval
        ticks = ax.ticks
        imax = findmax(diff(ticks))[2]
        push!(ticks, (ticks[imax] + ticks[imax+1]) / 2)
        sort!(ticks)
        typeof(ax)(int, ticks)
    else
        ax
    end
end

multiplicity(g::DiscreteAxis{T, :infinite, :infinite, I}, ::Type{Cartesian}) where {T, I} = one(T)
multiplicity(g::DiscreteAxis{T, :reflecting, :infinite, I}, ::Type{Cartesian}) where {T, I} = T(2)
multiplicity(g::DiscreteAxis{T, :infinite, :reflecting, I}, ::Type{Cartesian}) where {T, I} = T(2)
multiplicity(g::DiscreteAxis{T, :fixed, :fixed, I}, ::Type{Cartesian}) where {T, I} = one(T)
multiplicity(g::DiscreteAxis{T, :fixed, :reflecting, I}, ::Type{Cartesian}) where {T, I} = T(2)
multiplicity(g::DiscreteAxis{T, :reflecting, :fixed, I}, ::Type{Cartesian}) where {T, I} = T(2)
multiplicity(g::DiscreteAxis{T, :fixed, :infinite, I}, ::Type{Cartesian}) where {T, I} = one(T)
multiplicity(g::DiscreteAxis{T, :infinite, :fixed, I}, ::Type{Cartesian}) where {T, I} = one(T)

function multiplicity(g::DiscreteAxis{T, :reflecting, :reflecting, I}, ::Type{Cartesian}) where {T, I}
    @warn "Multiplicity of Cartesian axis, $(g) (:reflecting, :reflecting), would be infinite. It is set to 1 here."
    one(T)
end


text:
## https://developer.spotify.com/documentation/web-api/reference/#/operations/get-an-album

"""
    album_get_single(album_id; market="US")

**Summary**: Get Spotify catalog information for a single album.

# Arguments
- `album_id` : The Spotify ID for the album.

# Optional keywords
- `market` : An ISO 3166-1 alpha-2 country code. Default is set to "US".

# Example
```julia-repl
julia> Spotify.album_get_single("5XgEM5g3xWEwL4Zr6UjoLo")[1]
JSON3.Object{Base.CodeUnits{UInt8, String}, Vector{UInt64}} with 18 entries:
  :album_type             => "album"
  :artists                => JSON3.Object[{…
  :copyrights             => JSON3.Object[{…
  :external_ids           => {…
  :external_urls          => {…
  :genres                 => Union{}[]
```
"""
function album_get_single(album_id; market="US")
    return spotify_request("albums/$album_id?market=$market")
end

## https://developer.spotify.com/documentation/web-api/reference/#/operations/get-an-albums-tracks

"""
    album_get_tracks(album_id; limit=20, offset=0, market="US")

**Summary**: Get Spotify catalog information about an album's tracks. Optional parameters
             can be used to limit the number of tracks returned.

# Arguments
- `album_id` : The Spotify ID for the album

# Optional keywords
- `limit` : The maximum number of tracks to return. Default is set to 20.
- `offset` : The index of the first track to return. Default is 0.
- `market` : An ISO 3166-1 alpha-2 country code. Default is set to "US".

# Example
```julia-repl
julia> Spotify.album_get_tracks("5XgEM5g3xWEwL4Zr6UjoLo")[1]
JSON3.Object{Base.CodeUnits{UInt8, String}, Vector{UInt64}} with 7 entries:
  :href     => "https://api.spotify.com/v1/albums/5XgEM5g3xWEwL4Zr6UjoLo/tracks?offset=0&limit=20&market=US"
  :items    => JSON3.Object[{…
  :limit    => 20
  :next     => "https://api.spotify.com/v1/albums/5XgEM5g3xWEwL4Zr6UjoLo/tracks?offset=20&limit=20&market=US"
```
"""
function album_get_tracks(album_id; limit=20, offset=0, market="US")
    return spotify_request("albums/$album_id/tracks?limit=$limit&offset=$offset&market=$market")
end

## https://developer.spotify.com/documentation/web-api/reference/#/operations/get-multiple-albums

"""
    album_get_multiple(album_ids; market="US")

**Summary**: Get Spotify catalog information for multiple albums identified by their Spotify IDs.

# Arguments
- `album_ids` : A comma-separated list of the Spotify IDs for the albums. Maximum: 20 IDs.

# Optional keywords
- `market` : An ISO 3166-1 alpha-2 country code. Default is set to "US".

# Example
```julia-repl
julia> Spotify.album_get_multiple("5XgEM5g3xWEwL4Zr6UjoLo, 2rpT0freJsmUmmPluVWqg5")[1]
JSON3.Object{Base.CodeUnits{UInt8, String}, Vector{UInt64}} with 1 entry:
  :albums => JSON3.Object[{…
```
"""
function album_get_multiple(album_ids; market="US")
    return spotify_request("albums?ids=$album_ids&market=$market")
end

## https://developer.spotify.com/documentation/web-api/reference/#/operations/get-users-saved-albums

"""
    album_get_saved(;limit::Int64=20, market::String="US", offset::Int64=0)

**Summary**: Get a list of the albums saved in the current Spotify user's 'Your Music' library.

# Optional keywords
- `limit::Int64` : Maximum number of items to return, default is set to 20
- `market::String` : An ISO 3166-1 alpha-2 country code. If a country code is specified,
                     only episodes that are available in that market will be returned.
                     Default is set to "US".
- `offset::Int64` : Index of the first item to return, default is set to 0

# Example
```julia-repl
julia> Spotify.album_get_saved()[1]
[ Info: We try the request without checking if current grant includes scope user-library-read.
JSON3.Object{Base.CodeUnits{UInt8, String}, Vector{UInt64}} with 7 entries:
  :href     => "https://api.spotify.com/v1/me/albums?offset=0&limit=20&market=US"
  :items    => JSON3.Object[{…
  :limit    => 20
```
"""
function album_get_saved(;limit::Int64=20, market::String="US", offset::Int64=0)

    return spotify_request("me/albums?limit=$limit&market=$market&offset=$offset";
    scope = "user-library-read")

end

## https://developer.spotify.com/documentation/web-api/reference/#/operations/check-users-saved-albums

"""
    album_get_contains(ids)

**Summary**: Check if one or more albums is already saved in the current Spotify user's 'Your Music' library.

# Arguments
- `ids` : A comma-separated list of the Spotify album IDs

# Example
```julia-repl
julia> Spotify.album_get_contains("382ObEPsp2rxGrnsizN5TX,1A2GTWGtFfWp7KSQTwWOyo")[1]
[ Info: We try the request without checking if current grant includes scope user-library-read.
2-element JSON3.Array{Bool, Base.CodeUnits{UInt8, String}, Vector{UInt64}}:
 0
 0
```
"""
function album_get_contains(ids)
    return spotify_request("me/albums/contains?ids=$ids"; scope = "user-library-read")
end


text:
#= src/spfft2.jl
=#

abstract type SpFFTPlan2{T,K} end

# complex transforms

struct cSpFFTPlan2{T<:SpFFTComplex,K} <: SpFFTPlan2{T,K}
  P1::cSpFFTPlan1{T,K}
  P2::cSpFFTPlan1{T,K}
  k1n2::Matrix{T}
  n2k1::Matrix{T}
  k2k1::Matrix{T}
end

for (fcn, K) in ((:fft, FORWARD), (:bfft, BACKWARD))
  sf  = Symbol("sp", fcn)
  psf = Symbol("plan_", sf)
  f2s = Symbol(sf, "_f2s!")
  s2f = Symbol(sf, "_s2f!")
  @eval begin
    function $psf(::Type{T}, n1::Integer, n2::Integer,
                  idx1::AbstractVector{Ti}, idx2::AbstractVector{Ti};
                  args...) where {T<:SpFFTComplex,Ti<:Integer}
      P1 = $psf(T, n1, idx1; args...)
      P2 = $psf(T, n2, idx2; args...)
      k1 = length(idx1)
      k2 = length(idx2)
      k1n2 = Array{T}(undef, k1, n2)
      n2k1 = Array{T}(undef, n2, k1)
      k2k1 = Array{T}(undef, k2, k1)
      cSpFFTPlan2{T,$K}(P1, P2, k1n2, n2k1, k2k1)
    end

    function $f2s(Y::AbstractMatrix{Ty}, P::cSpFFTPlan2{T,K},
                  X::AbstractMatrix{Tx}; args...) where {T,K,Tx,Ty}
      l1, m1, k1 = spfft_size(P.P1)
      l2, m2, k2 = spfft_size(P.P2)
      n1 = l1*m1
      n2 = l2*m2
      size(X) == (n1, n2) || throw(DimensionMismatch)
      size(Y) == (k1, k2) || throw(DimensionMismatch)
      @inbounds for i = 1:n2
        $f2s(view(P.k1n2,:,i), P.P1, view(X,:,i); args...)
      end
      transpose!(P.n2k1, P.k1n2)
      @inbounds for i = 1:k1
        $f2s(view(P.k2k1,:,i), P.P2, view(P.n2k1,:,i); args...)
      end
      transpose_f!(z->spfft_rc(Ty,z), Y, P.k2k1)
    end

    function $s2f(Y::AbstractMatrix{Ty}, P::cSpFFTPlan2{T,K},
                  X::AbstractMatrix{Tx}; args...) where {T,K,Tx,Ty}
      l1, m1, k1 = spfft_size(P.P1)
      l2, m2, k2 = spfft_size(P.P2)
      n1 = l1*m1
      n2 = l2*m2
      size(X) == (k1, k2) || throw(DimensionMismatch)
      size(Y) == (n1, n2) || throw(DimensionMismatch)
      transpose!(P.k2k1, X)
      @inbounds for i = 1:k1
        $s2f(view(P.n2k1,:,i), P.P2, view(P.k2k1,:,i); args...)
      end
      transpose!(P.k1n2, P.n2k1)
      @inbounds for i = 1:n2
        $s2f(view(Y,:,i), P.P1, view(P.k1n2,:,i); args...)
      end
      Y
    end
  end
end

# real transforms

## (r2c, f2s)

struct rSpFFTPlan2{T<:SpFFTReal,Tc<:SpFFTComplex} <: SpFFTPlan2{T,FORWARD}
  P1::rSpFFTPlan1{T}
  P2::cSpFFTPlan1{Tc,FORWARD}
  k1n2::Matrix{Tc}
  n2k1::Matrix{Tc}
  k2k1::Matrix{Tc}
end

function plan_sprfft(::Type{T}, n1::Integer, n2::Integer,
                     idx1::AbstractVector{Ti}, idx2::AbstractVector{Ti};
                     args...) where {T<:SpFFTReal,Ti<:Integer}
  Tc = Complex{T}
  P1 = plan_sprfft(T, n1, idx1; args...)
  P2 = plan_spfft(Tc, n2, idx2; args...)
  k1 = length(idx1)
  k2 = length(idx2)
  k1n2 = Array{Tc}(undef, k1, n2)
  n2k1 = Array{Tc}(undef, n2, k1)
  k2k1 = Array{Tc}(undef, k2, k1)
  rSpFFTPlan2{T,Tc}(P1, P2, k1n2, n2k1, k2k1)
end

function sprfft_f2s!(
    Y::AbstractMatrix{Ty}, P::rSpFFTPlan2{T}, X::AbstractMatrix{Tx};
    args...) where {T,Tx<:Real,Ty<:Complex}
  l1, m1, k1 = spfft_size(P.P1)
  l2, m2, k2 = spfft_size(P.P2)
  n1 = l1*m1
  n2 = l2*m2
  size(X) == (n1, n2) || throw(DimensionMismatch)
  size(Y) == (k1, k2) || throw(DimensionMismatch)
  @inbounds for i = 1:n2
    sprfft_f2s!(view(P.k1n2,:,i), P.P1, view(X,:,i); args...)
  end
  transpose!(P.n2k1, P.k1n2)
  @inbounds for i = 1:k1
    spfft_f2s!(view(P.k2k1,:,i), P.P2, view(P.n2k1,:,i); args...)
  end
  transpose!(Y, P.k2k1)
end

## (c2r, s2f)
## - input must contain only nonredundant frequencies, i.e., up to index
##   n1/2 + 1 in the first dimension for a full size of n1

struct brSpFFTPlan2{T<:SpFFTReal,Tc<:SpFFTComplex} <: SpFFTPlan2{T,BACKWARD}
  P1::brSpFFTPlan1{T,Tc}
  P2::cSpFFTPlan1{Tc,BACKWARD}
  k1n2::Matrix{Tc}
  n2k1::Matrix{Tc}
  k2k1::Matrix{Tc}
end

function plan_spbrfft(::Type{T}, n1::Integer, n2::Integer,
                      idx1::AbstractVector{Ti}, idx2::AbstractVector{Ti};
                      args...) where {T<:SpFFTReal,Ti<:Integer}
  Tc = Complex{T}
  P1 = plan_spbrfft(T, n1, idx1; args...)
  P2 = plan_spbfft(Tc, n2, idx2; args...)
  k1 = length(idx1)
  k2 = length(idx2)
  k1n2 = Array{T }(undef, k1, n2)
  n2k1 = Array{Tc}(undef, n2, k1)
  k2k1 = Array{Tc}(undef, k2, k1)
  brSpFFTPlan2{T,Tc}(P1, P2, k1n2, n2k1, k2k1)
end

function spbrfft_s2f!(
    Y::AbstractMatrix{Ty}, P::brSpFFTPlan2{T}, X::AbstractMatrix{Tx};
    args...) where {T,Tx<:Complex,Ty<:Real}
  l1, m1, kc1, k1 = spbrfft_size(P.P1)
  l2, m2, k2 = spfft_size(P.P2)
  n1 = l1*m1
  n2 = l2*m2
  size(X) == (k1, k2) || throw(DimensionMismatch)
  size(Y) == (n1, n2) || throw(DimensionMismatch)
  transpose!(P.k2k1, X)
  @inbounds for i = 1:k1
    spbfft_s2f!(view(P.n2k1,:,i), P.P2, view(P.k2k1,:,i); args...)
  end
  transpose!(P.k1n2, P.n2k1)
  @inbounds for i = 1:n2
    spbrfft_s2f!(view(Y,:,i), P.P1, view(P.k1n2,:,i); args...)
  end
  Y
end


text:
############################################################
# joLinearFunction - operator  constructors ################
############################################################

# FFT operators: joDFT
include("joLinearFunctionConstructors/joDFT.jl")

# DCT operators: joDCT
include("joLinearFunctionConstructors/joDCT.jl")

# NFFT operators: joNFFT
include("joLinearFunctionConstructors/joNFFT.jl")

# DWT operators: joDWT
include("joLinearFunctionConstructors/joDWT.jl")

# SWT operators: joSWT
include("joLinearFunctionConstructors/joSWT.jl")

# Romberg operator
include("joLinearFunctionConstructors/joRomberg.jl")

# CurveLab operators: joCurvelet2D joCurvelet2DnoFFT
include("joLinearFunctionConstructors/joCurvelet2D.jl")
include("joLinearFunctionConstructors/joCurvelet2DnoFFT.jl")

# Restriction operator
include("joLinearFunctionConstructors/joRestriction.jl")

# Mask operator
include("joLinearFunctionConstructors/joMask.jl")

# Padding/extension operators: joExtend
include("joLinearFunctionConstructors/joExtend.jl")

# Permutation operator: joPermutation
include("joLinearFunctionConstructors/joPermutation.jl")

# Outer product operator: joOuterProd
include("joLinearFunctionConstructors/joOuterProd.jl")


text:
using FeldtLib: get_lat_lon, vgregion_lat_lon_bounds

@testset "geocoding" begin

@testset "test one known address" begin
    bgc, bgd = get_lat_lon("Bratteråsgatan, Sweden")
    # Manually taken from Nominatim 2020-04-23:
    @test bgc[1] == 57.7021997
    @test bgc[2] == 11.9131308
end

@testset "check that address of region Västra Götaland is within its bounding box" begin
    vgrcoord, vgrd = get_lat_lon("Västra Götaland, Sweden")
    lowleftlon, lowleftlat, uprightlon, uprightlat = vgregion_lat_lon_bounds()
    vgrlat, vgrlon = vgrcoord
    @test lowleftlat < vgrlat < uprightlat
    @test lowleftlon < vgrlon < uprightlon
end

end


text:
macro ArfFloat(x)
    convert(ArfFloat, string(:($x)))
end
macro ArfFloat(p,x)
    convert(ArfFloat{:($p)}, string(:($x)))
end

macro ArbFloat(x)
    convert(ArbFloat, string(:($x)))
end
macro ArbFloat(p,x)
    convert(ArbFloat{:($p)}, string(:($x)))
end

convert{T<:ArfFloat}(::Type{T}, x::T) = x
convert{T<:ArbFloat}(::Type{T}, x::T) = x
convert{P}(::Type{ArfFloat{P}}, x::ArfFloat{P}) = x
convert{P}(::Type{ArbFloat{P}}, x::ArbFloat{P}) = x

function convert{Q}(::Type{ArfFloat}, x::ArfFloat{Q})
   P = precision(ArfFloat)
   z = ArfFloat{P}()
   ccall(@libarb(arf_set_round), Void, (Ptr{ArfFloat{P}}, Ptr{ArfFloat{Q}}, Clong), &z, &x, Clong(P))
   return z
end
function convert{P,Q}(::Type{ArfFloat{P}}, x::ArfFloat{Q})
   z = ArfFloat{P}()
   ccall(@libarb(arf_set_round), Void, (Ptr{ArfFloat{P}}, Ptr{ArfFloat{Q}}, Clong), &z, &x, Clong(P))
   return z
end

function convert{Q}(::Type{ArbFloat}, x::ArbFloat{Q})
   P = precision(ArbFloat)
   z = ArbFloat{P}()
   ccall(@libarb(arb_set_round), Void, (Ptr{ArbFloat{P}}, Ptr{ArbFloat{Q}}, Clong), &z, &x, Clong(P))
   return z
end
function convert{P,Q}(::Type{ArbFloat{P}}, x::ArbFloat{Q})
   z = ArbFloat{P}()
   ccall(@libarb(arb_set_round), Void, (Ptr{ArbFloat{P}}, Ptr{ArbFloat{Q}}, Clong), &z, &x, Clong(P))
   return z
end

function convert{P}(::Type{ArbFloat{P}}, x::ArfFloat{P})
   z = ArbFloat{P}()
   ccall(@libarb(arb_set_arf), Void, (Ptr{ArbFloat{P}}, Ptr{ArfFloat{P}}), &z, &x)
   return z
end

function convert{P}(::Type{ArfFloat{P}}, x::ArbFloat{P})
    z = ArfFloat{P}()
    z.exponentOf2  = x.exponentOf2
    z.nwords_sign  = x.nwords_sign
    z.significand1 = x.significand1
    z.significand2 = x.significand2
    return z
end

function convert{P,Q}(::Type{ArbFloat{P}}, x::ArfFloat{Q})
    y = convert(ArfFloat{P}, x)
    z = convert(ArbFloat{P}, y)
    return z
end
function convert{Q}(::Type{ArbFloat}, x::ArfFloat{Q})
    P = precision(ArbFloat)
    y = convert(ArfFloat{P}, x)
    z = convert(ArbFloat{P}, y)
    return z
end

function convert{P,Q}(::Type{ArfFloat{P}}, x::ArbFloat{Q})
    y = convert(ArbFloat{P}, x)
    z = convert(ArfFloat{P}, y)
    return z
end
function convert{Q}(::Type{ArfFloat}, x::ArbFloat{Q})
    P = precision(ArfFloat)
    y = convert(ArbFloat{P}, x)
    z = convert(ArfFloat{P}, y)
    return z
end

# convert ArbFloat with other types

function convert{T<:ArbFloat}(::Type{T}, x::UInt64)
    z = T()
    ccall(@libarb(arb_set_ui), Void, (Ptr{T}, UInt64), &z, x)
    return z
end
function convert{T<:ArbFloat}(::Type{T}, x::Int64)
    z = T()
    ccall(@libarb(arb_set_si), Void, (Ptr{T}, Int64), &z, x)
    return z
end

function convert{T<:ArbFloat}(::Type{T}, x::Float64)
    z = T()
    ccall(@libarb(arb_set_d), Void, (Ptr{T}, Float64), &z, x)
    return z
end

function convert{T<:ArbFloat}(::Type{T}, x::String)
    P = precision(T)
    z = T()
    ccall(@libarb(arb_set_str), Void, (Ptr{T}, Ptr{UInt8}, Int), &z, x, P)
    return z
end

convert{T<:ArbFloat}(::Type{T}, x::UInt32)  = convert(T, x%UInt64)
convert{T<:ArbFloat}(::Type{T}, x::UInt16)  = convert(T, x%UInt64)
convert{T<:ArbFloat}(::Type{T}, x::UInt8)   = convert(T, x%UInt64)
convert{T<:ArbFloat}(::Type{T}, x::UInt128) = convert(T, string(x))

convert{T<:ArbFloat}(::Type{T}, x::Int32)  = convert(T, x%Int64)
convert{T<:ArbFloat}(::Type{T}, x::Int16)  = convert(T, x%Int64)
convert{T<:ArbFloat}(::Type{T}, x::Int8)   = convert(T, x%Int64)
convert{T<:ArbFloat}(::Type{T}, x::Int128) = convert(T, string(x))

convert{T<:ArbFloat}(::Type{T}, x::Float32) = convert(T, convert(Float64,x))
convert{T<:ArbFloat}(::Type{T}, x::Float16) = convert(T, convert(Float64,x))

function convert{T<:ArbFloat}(::Type{Float64}, x::T)
    ptr2mid = ptr_to_midpoint(x)
    fl = ccall(@libarb(arf_get_d), Float64, (Ptr{ArfFloat}, Int), ptr2mid, 4) # round nearest
    return fl
end

function convert{T<:ArbFloat}(::Type{Float32}, x::T)
    return convert(Float32, convert(Float64, x))
end
function convert{T<:ArbFloat}(::Type{Float16}, x::T)
    return convert(Float16, convert(Float64, x))
end

for I in (:UInt64, :UInt128)
  @eval begin
    function convert{T<:ArbFloat}(::Type{$I}, x::T)
        if isinteger(x)
           if notnegative(x)
               return convert($I, convert(BigInt,x))
           else
               throw( DomainError() )
           end
        else
           throw( InexactError() )
        end
    end
  end
end

convert{T<:ArbFloat}(::Type{UInt32}, x::T) = convert(UInt32, convert(UInt64,x))
convert{T<:ArbFloat}(::Type{UInt16}, x::T) = convert(UInt16, convert(UInt64,x))
convert{T<:ArbFloat}(::Type{UInt8}, x::T) = convert(UInt8, convert(UInt64,x))

for I in (:Int64, :Int128)
  @eval begin
    function convert{T<:ArbFloat}(::Type{$I}, x::T)
        if isinteger(x)
           return convert($I, convert(BigInt,x))
        else
           throw( InexactError() )
        end
    end
  end
end

convert{T<:ArbFloat}(::Type{Int32}, x::T) = convert(Int32, convert(Int64,x))
convert{T<:ArbFloat}(::Type{Int16}, x::T) = convert(Int16, convert(Int64,x))
convert{T<:ArbFloat}(::Type{Int8}, x::T) = convert(Int8, convert(Int64,x))

function parse{T<:ArbFloat}(::Type{T}, x::String)
    return T(x)
end

# =================

convert(::Type{BigInt}, x::String) = parse(BigInt,x)
convert(::Type{BigFloat}, x::String) = parse(BigFloat,x)

function convert{T<:ArbFloat}(::Type{T}, x::BigFloat)
     P = precision(T)+24
     x = round(x,P,2)
     s = string(x)
     z = T(s)
     return z
end

#=
function convert{T<:ArbFloat}(::Type{BigFloat}, x::T)
     s = string(midpoint(x))
     return parse(BigFloat, s)
end
=#
function convert{T<:ArbFloat}(::Type{BigFloat}, x::T)
    ptr2mid = ptr_to_midpoint(x)
    bf = zero(BigFloat)
    rounddir = ccall(@libarb(arf_get_mpfr), Int, (Ptr{BigFloat}, Ptr{ArfFloat}, Int), &bf, ptr2mid, 4) # round nearest
    return bf
end

function convert{I<:Integer,P}(::Type{Rational{I}}, x::ArbFloat{P})
    bf = convert(BigFloat, x)
    return convert(Rational{I}, bf)
end

for T in (:Integer, :Signed)
  @eval begin
    function convert{P}(::Type{$T}, x::ArbFloat{P})
        y = trunc(x)
        try
           return convert(Int64, x)
        catch
           try
              return convert(Int128, x)
           catch
              DomainError()
           end
        end
    end
  end
end

for F in (:BigInt, :Rational)
  @eval begin
    function convert{T<:ArbFloat}(::Type{T}, x::$F)
        P = precision(T)
        B = precision(BigFloat)
        if B < P+24
            return convert(ArbFloat{P}, string(x))
        else
            return convert(ArbFloat{P}, convert(BigFloat, x))
        end
    end
    function convert{P}(::Type{ArbFloat{P}}, x::$F)
        B = precision(BigFloat)
        if B < P+24
            return convert(ArbFloat{P}, string(x))
        else
            return convert(ArbFloat{P}, convert(BigFloat, x))
        end
    end
  end
end

function convert{T<:ArbFloat,S}(::Type{T}, x::Irrational{S})
    P = precision(T)
    a = ArbFloat{P}()
    setprecision(BigFloat, precision(T)+24) do
         a = convert(T, BigFloat(x))
    end
    return a
end

function convert{P}(::Type{BigInt}, x::ArbFloat{P})
   z = trunc(convert(BigFloat, x))
   return convert(BigInt, z)
end

# Promotion
for T in (:Int128, :Int64, :Int32, :Int16, :Float64, :Float32, :Float16,
          :(Rational{Int64}), :(Rational{Int32}), :(Rational{Int16}),
          :String)
  @eval promote_rule{P}(::Type{ArbFloat{P}}, ::Type{$T}) = ArbFloat{P}
end

float{P}(x::ArbFloat{P}) = x
float{T<:ArbFloat}(x::T) = x

promote_rule{P}(::Type{ArbFloat{P}}, ::Type{BigFloat}) = ArbFloat{P}
promote_rule{P}(::Type{ArbFloat{P}}, ::Type{BigInt}) = ArbFloat{P}
promote_rule{P}(::Type{ArbFloat{P}}, ::Type{Rational{BigInt}}) = Rational{BigInt}

promote_rule{P,Q}(::Type{ArbFloat{P}}, ::Type{ArbFloat{Q}}) =
    ifelse(P>Q, ArbFloat{P}, ArbFloat{Q})

@vectorize_1arg(BigFloat, ArbFloat)
@vectorize_1arg(Float64, ArbFloat)
@vectorize_1arg(Float32, ArbFloat)
@vectorize_1arg(BigInt, ArbFloat)
@vectorize_1arg(Int128, ArbFloat)
@vectorize_1arg(Int64, ArbFloat)
@vectorize_1arg(Int32, ArbFloat)
@vectorize_1arg(Rational, ArbFloat)

@vectorize_1arg(ArbFloat, Float64)
@vectorize_1arg(ArbFloat, Float32)
@vectorize_1arg(ArbFloat, BigInt)
@vectorize_1arg(ArbFloat, Int128)
@vectorize_1arg(ArbFloat, Int64)
@vectorize_1arg(ArbFloat, Int32)
@vectorize_1arg(ArbFloat, Rationa)


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    Cylinder(start, finish, radius)

A right circular cylinder with `start` and `finish` points,
and `radius` of revolution. See https://en.wikipedia.org/wiki/Cylinder.
"""
struct Cylinder{T} <: Primitive{3,T}
  start::Point{3,T}
  finish::Point{3,T}
  radius::T
end

Cylinder(start::Tuple, finish::Tuple, radius) =
  Cylinder(Point(start), Point(finish), radius)

paramdim(::Type{<:Cylinder}) = 3

isconvex(::Type{<:Cylinder}) = true

radius(c::Cylinder) = c.radius
height(c::Cylinder) = norm(c.finish - c.start)

measure(c::Cylinder) = π * radius(c)^2 * height(c)


text:
# Case G example
# Claudia December 2014

tree = "((((6:0.1,4:1.5)1:0.2,(7)11#H1)5:0.1,(11#H1,8)),10:0.1);" # Case G
#f = open("prueba_tree.txt","w")
#write(f,tree)
#close(f)
net = readTopologyLevel1(tree)

#printEdges(net)
#printNodes(net)


text:
function parsebinaryexp(ordering::Ordering,value::BinBoolType)
  return OBDD(ordering,value)
end

function parsebinaryexp(ordering::Ordering,var::Symbol)
  return OBDD(ordering,BDD(string(var),BDD(false),BDD(true)))
end

function parseordering(ordering::Expr)
  if ordering.head != :tuple
    throw(Meta.ParseError("expected a tuple, got $(ordering)"))
  end

  return ListOrdering(String[string(var) for var in ordering.args])
end

function parsebinaryexp(ordering::Ordering,binaryexp::Expr)
  if binaryexp.head == :&&
    return (parsebinaryexp(ordering,binaryexp.args[1])&
            parsebinaryexp(ordering,binaryexp.args[2]))
  end

  if binaryexp.head == :||
    return (parsebinaryexp(ordering,binaryexp.args[1])|
            parsebinaryexp(ordering,binaryexp.args[2]))
  end

  if binaryexp.head == :->
    return ((~parsebinaryexp(ordering,binaryexp.args[1]))|
            parsebinaryexp(ordering,binaryexp.args[2].args[2]))
  end

  if binaryexp.head != :call
    throw(Meta.ParseError("expected a bitwise binary operator, got $(binaryexp.head)"))
  end

  if binaryexp.args[1] == :~ || binaryexp.args[1] == :!
    return ~parsebinaryexp(ordering,binaryexp.args[2])
  end

  if binaryexp.args[1] == :& || binaryexp.args[1] == :*
    return (parsebinaryexp(ordering,binaryexp.args[2])&
            parsebinaryexp(ordering,binaryexp.args[3]))
  end

  if binaryexp.args[1] == :| || binaryexp.args[1] == :+
    return (parsebinaryexp(ordering,binaryexp.args[2])|
            parsebinaryexp(ordering,binaryexp.args[3]))
  end

  if binaryexp.args[1] == :$
    return (parsebinaryexp(ordering,binaryexp.args[2]) $
            parsebinaryexp(ordering,binaryexp.args[3]))
  end

  throw(Meta.ParseError("unknown bitwise binary operator $(binaryexp.args[1])"))
end

function parsebinaryfunct(binaryfunct::Expr)
  if binaryfunct.head == :->
    O=parseordering(binaryfunct.args[1])

    #if ordering!=nothing && ordering != O
    #  throw(ArgumentError("the parsed order $(O) and the passed order $(ordering) are not the same"))
    #end

    return parsebinaryexp(O,binaryfunct.args[2].args[2])
  end

  throw(Meta.ParseError("unknown bitwise binary function $(binaryfunct)"))
end

function OBDD(ordering::Ordering,binaryexpstr::String)
  binaryexp=Meta.parse(binaryexpstr)

  return parsebinaryexp(ordering,binaryexp)
end

function OBDD(ordering::Array{String,1},binaryexpstr::String)
  return OBDD(ListOrdering(ordering),binaryexpstr)
end

function OBDD(ordering::Array{Char,1},binaryexpstr::String)
  return OBDD(ListOrdering(ordering),binaryexpstr)
end

function OBDD(binaryfunctstr::String)
  binaryfunct=Meta.parse(binaryfunctstr)

  return parsebinaryfunct(binaryfunct)
end

function changeordering(A::OBDD,O::Ordering)
  return OBDD(O,string(A.root))
end

function changeordering(A::OBDD,O::Array{String,1})
  return changeordering(A,ListOrdering(O))
end

function changeordering(A::OBDD,O::Array{Char,1})
  return changeordering(A,ListOrdering(O))
end

function changeordering!(A::OBDD,O::Ordering)
  A.root=OBDD(O,string(A.root)).root
  A.ordering=O

  return A
end

function changeordering!(A::OBDD,O::Array{String,1})
  return changeordering!(A,ListOrdering(O))
end

function changeordering!(A::OBDD,O::Array{Char,1})
  return changeordering!(A,ListOrdering(O))
end


text:
include("LS-method-7s-init.jl");

begin
    @load "LS_pureIOU_2000_narrow_saved.jld2" LS_pureIOU_2000_narrow_saved;
    r_list=LS_pureIOU_2000_narrow_saved[1]
    a_list=LS_pureIOU_2000_narrow_saved[2]
    Λ1=LS_pureIOU_2000_narrow_saved[3]
    tbl2=LS_pureIOU_2000_narrow_saved[4]
end;

## Two important checkpoints!
sum(Λ1) # should print 1
[sum(Λ1[i,:]) for i=1:nos]-BSP # should print a list of zeros

begin
    A=hcat(r_list,a_list)
    Aplus=A[A[:,2].>0,:]
    Aminus=A[A[:,2].<=0,:]
    iloc1=findmin(Aplus[:,2])[2]
    iloc2=findmax(Aminus[:,2])[2]
end;

#the nearest interest rates
Aplus[iloc1,1],Aminus[iloc2,1]

#the respective aggregate demand
Aplus[iloc1,2],Aminus[iloc2,2]

#the gap in the demand a(r)
Aplus[iloc1,2]-Aminus[iloc2,2]

#the gap in the interest rate is
Aplus[iloc1,1]-Aminus[iloc2,1]

begin
    scatter(a_list,r_list,label="",color=:black,border=true)
    plot!((0.02:0.0001:0.032)*0,0.02:0.0001:0.032,color=:black,label="")
    xlabel!("average assets")
    ylabel!("interest rate")
end

begin
    r=r_list[end]
    w_bounds=(-1.6292720455478487, 7.18257249498523);
    Lg=2000 #the only place to change the grid size
    grda=[w_bounds[1]+i*(w_bounds[end]-w_bounds[1])/(Lg-1) for i=0:(Lg-1)]
    β=exp(-ρ)
end;

begin
    bond_price=ι/(1+r)
    cons=[(1+r)*grda[i]+wage*hours[k]-grda[Int64(tbl2[k,i])] for k=1:nos, i=1:Lg];
    aassets=[grda[Int64(tbl2[k,i])] for k=1:nos, i=1:Lg];
    shares=[grda[Int64(tbl2[k,i])]/bond_price for k=1:nos, i=1:Lg];
end;

begin
    SL=[sum(Λ1[state,:]) for state=1:nos];
    LF=[[Λ1[state,1]/SL[state]] for state=1:nos];
    for i=2:Lg
        for state=1:nos
            LF[state]=vcat(LF[state],LF[state][end]+Λ1[state,i]/SL[state])
        end
    end
end;

begin
    plot(aassets[1,2:920],LF[1][2:920],label="(LS)",color=:black)
    plot!(aassets[2,2:920],LF[2][2:920],label="",color=:black)
    plot!(aassets[3,2:920],LF[3][2:920],label="",color=:black)
    plot!(aassets[4,2:920],LF[4][2:920],label="",color=:black)
    plot!(aassets[5,2:920],LF[5][2:920],label="",color=:black)
    plot!(aassets[6,2:920],LF[6][2:920],label="",color=:black)
    plot!(aassets[7,2:920],LF[7][2:920],label="",color=:black)
end


text:
export BrickletIndustrialAnalogOutConfiguration
struct BrickletIndustrialAnalogOutConfiguration
    voltage_range::Integer
    current_range::Integer
end

export BrickletIndustrialAnalogOutIdentity
struct BrickletIndustrialAnalogOutIdentity
    uid::String
    connected_uid::String
    position::Char
    hardware_version::Vector{Integer}
    firmware_version::Vector{Integer}
    device_identifier::Integer
end

export BrickletIndustrialAnalogOut
"""
Generates configurable DC voltage and current, 0V to 10V and 4mA to 20mA
"""
mutable struct BrickletIndustrialAnalogOut <: TinkerforgeDevice
    ipcon::IPConnection
    deviceInternal::PyObject

    """
    Creates an object with the unique device ID *uid* and adds it to
    the IP Connection *ipcon*.
    """
    function BrickletIndustrialAnalogOut(uid::String, ipcon::IPConnection)
        package = pyimport("tinkerforge.bricklet_industrial_analog_out")
        deviceInternal = package.BrickletIndustrialAnalogOut(uid, ipcon.ipconInternal)

        return new(ipcon, deviceInternal)
    end
end

export enable
"""
    $(SIGNATURES)

Enables the output of voltage and current.

The default is disabled.
"""
function enable(device::BrickletIndustrialAnalogOut)
    device.deviceInternal.enable()
end

export disable
"""
    $(SIGNATURES)

Disables the output of voltage and current.

The default is disabled.
"""
function disable(device::BrickletIndustrialAnalogOut)
    device.deviceInternal.disable()
end

export is_enabled
"""
    $(SIGNATURES)

Returns *true* if output of voltage and current is enabled, *false* otherwise.
"""
function is_enabled(device::BrickletIndustrialAnalogOut)
    return device.deviceInternal.is_enabled()
end

export set_voltage
"""
    $(SIGNATURES)

Sets the output voltage.

The output voltage and output current are linked. Changing the output voltage
also changes the output current.
"""
function set_voltage(device::BrickletIndustrialAnalogOut, voltage)
    device.deviceInternal.set_voltage(voltage)
end

export get_voltage
"""
    $(SIGNATURES)

Returns the voltage as set by :func:`Set Voltage`.
"""
function get_voltage(device::BrickletIndustrialAnalogOut)
    return device.deviceInternal.get_voltage()
end

export set_current
"""
    $(SIGNATURES)

Sets the output current.

The output current and output voltage are linked. Changing the output current
also changes the output voltage.
"""
function set_current(device::BrickletIndustrialAnalogOut, current)
    device.deviceInternal.set_current(current)
end

export get_current
"""
    $(SIGNATURES)

Returns the current as set by :func:`Set Current`.
"""
function get_current(device::BrickletIndustrialAnalogOut)
    return device.deviceInternal.get_current()
end

export set_configuration
"""
    $(SIGNATURES)

Configures the voltage and current range.

Possible voltage ranges are:

* 0V to 5V
* 0V to 10V

Possible current ranges are:

* 4mA to 20mA
* 0mA to 20mA
* 0mA to 24mA

The resolution will always be 12 bit. This means, that the
precision is higher with a smaller range.
"""
function set_configuration(device::BrickletIndustrialAnalogOut, voltage_range, current_range)
    device.deviceInternal.set_configuration(voltage_range, current_range)
end

export get_configuration
"""
    $(SIGNATURES)

Returns the configuration as set by :func:`Set Configuration`.
"""
function get_configuration(device::BrickletIndustrialAnalogOut)
    return device.deviceInternal.get_configuration()
end

export get_identity
"""
    $(SIGNATURES)

Returns the UID, the UID where the Bricklet is connected to,
the position, the hardware and firmware version as well as the
device identifier.

The position can be 'a', 'b', 'c', 'd', 'e', 'f', 'g' or 'h' (Bricklet Port).
A Bricklet connected to an :ref:`Isolator Bricklet <isolator_bricklet>` is always at
position 'z'.

The device identifier numbers can be found :ref:`here <device_identifier>`.
|device_identifier_constant|
"""
function get_identity(device::BrickletIndustrialAnalogOut)
    return device.deviceInternal.get_identity()
end


text:
const BC = Base.Broadcast

struct CAStyle{InnerStyle<:BC.BroadcastStyle, Axes, N} <: BC.AbstractArrayStyle{N} end
CAStyle(::InnerStyle, ::Axes, N) where {InnerStyle, Axes} = CAStyle{InnerStyle, Axes, N}()
CAStyle(::InnerStyle, ::Type{<:Axes}, N) where {InnerStyle, Axes} = CAStyle{InnerStyle, Axes, N}()

function CAStyle(::InnerStyle, ax::Axes, ::Val{N}) where {InnerStyle, Axes, N}
    return CAStyle(InnerStyle(), ax, N)
end

function Base.BroadcastStyle(::Type{<:ComponentArray{T, N, A, Axes}}) where {T, A, N, Axes}
    return CAStyle(Base.BroadcastStyle(A), getaxes(Axes), ndims(A))
end
function Base.BroadcastStyle(AA::Type{<:Adjoint{T, <:ComponentArray{T, N, A, Axes}}}) where {T, N, A, Axes}
    return CAStyle(Base.BroadcastStyle(Adjoint{T,A}), getaxes(AA), ndims(AA))
end
function Base.BroadcastStyle(AA::Type{<:Transpose{T, <:ComponentArray{T, N, A, Axes}}}) where {T, N, A, Axes}
    return CAStyle(Base.BroadcastStyle(Transpose{T,A}), getaxes(AA), ndims(AA))
end

function Base.BroadcastStyle(::CAStyle{InnerStyle, Axes, N}, bc::BC.Broadcasted) where {InnerStyle, Axes, N}
    return CAStyle(Base.BroadcastStyle(InnerStyle(), bc), Axes, N)
end

function BC.BroadcastStyle(::CAStyle{<:In1, <:Ax1, <:N1}, ::CAStyle{<:In2, <:Ax2, <:N2}) where {In1, Ax1, N1, In2, Ax2, N2}
    ax, N = fill_flat(Ax1, Ax2, N1, N2)
    inner_style = BC.BroadcastStyle(In1(), In2())
    if inner_style isa BC.Unknown
        inner_style = BC.DefaultArrayStyle{N}()
    end
    return CAStyle(inner_style, ax, N)
end
function BC.BroadcastStyle(::CAStyle{In, Ax, N1}, ::Style) where Style<:BC.DefaultArrayStyle{N2} where {In, Ax, N1, N2}
    N = max(N1, N2)
    ax = fill_flat(Ax, max(N1, N2))
    inner_style = BC.BroadcastStyle(In(), Style())
    return CAStyle(inner_style, ax, N)
end
function BC.BroadcastStyle(CAS::CAStyle{In, Ax, N1}, ::BC.DefaultArrayStyle{0}) where {In, Ax, N1}
    return CAS
end
function BC.BroadcastStyle(CAS::CAStyle{In, Ax, N}, ::BC.DefaultArrayStyle{N}) where {In, Ax, N}
    return CAS
end
function BC.BroadcastStyle(::CAStyle{In, Ax, N1}, ::Style) where Style<:BC.AbstractArrayStyle{N2} where {In, Ax, N1, N2}
    N = max(N1, N2)
    ax = fill_flat(Ax, max(N1, N2))
    inner_style = BC.BroadcastStyle(In(), Style())
    return CAStyle(inner_style, ax, N)
end

Base.convert(::Type{<:BC.Broadcasted{Nothing}}, bc::BC.Broadcasted{<:CAStyle,Axes,F,Args}) where {Axes,F,Args} = getdata(bc)

getdata(bc::BC.Broadcasted{<:CAStyle}) = BC.broadcasted(bc.f, map(getdata, bc.args)...)

function Base.similar(bc::BC.Broadcasted{<:CAStyle{InnerStyle, Axes, N}}, args...) where {InnerStyle, Axes, N}
    return ComponentArray{Axes}(similar(BC.Broadcasted{InnerStyle}(bc.f, bc.args, bc.axes), args...))
end
function Base.similar(bc::BC.Broadcasted{<:CAStyle{InnerStyle, Axes, N}}, T::Type) where {InnerStyle, Axes, N}
    return ComponentArray{Axes}(similar(BC.Broadcasted{InnerStyle}(bc.f, bc.args, bc.axes), T))
end
function Base.similar(bc::BC.Broadcasted{<:CAStyle{<:BC.Unknown, Axes, N}}, T::Type) where {InnerStyle, Axes, N}
    return ComponentArray{Axes}(similar(BC.Broadcasted{BC.DefaultArrayStyle{N}}(bc.f, bc.args, bc.axes), T))
end

Base.Broadcast.broadcasted(f, x::ComponentArray) = ComponentArray(map(f, x), getaxes(x))

# function Base.copy(bc::BC.Broadcasted{<:CAStyle{InnerStyle, Axes, N}}) where {InnerStyle, Axes,  N}
#     return ComponentArray{Axes}(Base.copy(BC.broadcasted(bc.f, map(getdata, bc.args)...)))
# end
# function Base.copy(bc::BC.Broadcasted{<:CAStyle{InnerStyle, Axes, N}}) where {InnerStyle, Axes,  N}
#     return ComponentArray{Axes}(Base.copy(BC.Broadcasted(InnerStyle())))
# end

# From https://github.com/JuliaArrays/OffsetArrays.jl/blob/master/src/OffsetArrays.jl
Base.dataids(A::ComponentArray) = Base.dataids(parent(A))
Broadcast.broadcast_unalias(dest::ComponentArray, src) = getdata(dest) === getdata(src) ? src : Broadcast.unalias(dest, src)

# Helper for extruding axes
function fill_flat(Ax1, Ax2, N1, N2)
    if N1<N2
        N = N2
        ax1 = fill_flat(Ax1,N)
        ax2 = Ax2
    elseif N1>N2
        N = N1
        ax1 = Ax1
        ax2 = fill_flat(Ax2,N)
    else
        N = N1
        ax1, ax2 = Ax1, Ax2
    end
    Ax = promote.(getaxes(ax1), getaxes(ax2)) |> typeof
    return Ax, N
end
fill_flat(Ax::Type{<:VarAxes}, N) = fill_flat(getaxes(Ax), N) |> typeof
function fill_flat(Ax::VarAxes, N)
    axs = Ax
    n = length(axs)
    if N>n
        axs = (axs..., ntuple(x -> FlatAxis(), N-n)...)
    end
    return axs
end


text:
module Boltzmann

export AbstractRBM,
       RBM,
       BernoulliRBM,
       GRBM,
       ConditionalRBM,
       DBN,
       DAE,
       Bernoulli,
       Gaussian,
       Degenerate,
       fit,
       transform,
       generate,
       predict,
       coef,
       weights,
       hbias,
       vbias,
       unroll,
       save_params,
       load_params

include("core.jl")

end


text:
"""
    function addCheckBounds_!(p::Program_,platform::String)

Generate the functions related with boundary checking.
"""
function addCheckBounds_!(p::Program_,platform::String)

    if any([typeof(i) in [Periodic,Bounded] for i in p.agent.boundary.boundaries])
        push!(p.declareVar.args, #Checking that the particles are initially inside the simulation box
        Meta.parse(string("AgentBasedModels.checkInitiallyInBound",
                        typeof(p.agent.boundary),
                        "(Core.Array(simulationBox),Core.Array(localV),",
                        [typeof(i) for i in p.agent.boundary.boundaries],
                        ")"
                    )

                )
        )
    end

    if !isempty([i for i in keys(p.update["Local"])])
        ##Add bound code if bound symbols update
        code = returnBound_(p.agent.boundary,p) #This is defined in the boundary definitions

        if !emptyquote_(code)

            #Construct functions
            f = simpleFirstLoopWrapInFunction_(platform,:boundsCheck_!,code)
            f = vectorize_(p.agent,f,p)

            push!(p.declareF.args,
                f)

            push!(p.execInloop.args,
                    :(@platformAdapt boundsCheck_!(ARGS_))
                )
        end

    end

    return nothing
end


text:
@testset "EikonalInv" begin
	include("testGetData2D.jl")
	include("testGetData3D.jl")
	include("testSensMatVec.jl")
	include("testTravelTimeInversion.jl")
end