text:
"""
```
SeisDecimate(d ; <keyword arguments>)
```

Decimate a multidimensional array of data. Input and output have the same dimensions.

# Arguments
* `d`: Input data in the form of  2D, or 3D,4D,5D tensors. The first dimension is time.

# Keyword arguments
* `mode=random`: Data is decimated randomly. Otherwise it is regularly decimated.
* `perc=50`: percentage of traces to be decimated from original array
* `incx1=1,incx2=1,incx3=1,incx4=1`: traces to decimate in case of regular decimation

# Output
* `out`: Decimated data.

# Example
```
julia> d = SeisLinearEvents(d); SeisPlot(d)
julia> d_dec = SeisDecimate(d); SeisPlot(d_dec)
```
"""

function SeisDecimate(in;mode="random",perc=50,incx1=1,incx2=1,incx3=1,incx4=1)

    if (mode=="random") # decimate data randomly
        out = copy(in)
	mask = rand(1,size(in,2)*size(in,3)*size(in,4)*size(in,5));
	mask[find(mask .< perc/100)] = 0;
	mask[find(mask .>= perc/100)] = 1;
	for it = 1 : size(in,1)
	    out[it,:] = out[it:it,:].*mask;
	end
    else # decimate data regularly with respect to 4 spatial dimensions
        out = zeros(in)
	if (size(in),5) > 1
	    out[:,1:incx1:end,1:incx2:end,1:incx3:end,1:incx4:end] = in[:,1:incx1:end,1:incx2:end,1:incx3:end,1:incx4:end]
	elseif (size(in),4) > 1
            out[:,1:incx1:end,1:incx2:end,1:incx3:end] = in[:,1:incx1:end,1:incx2:end,1:incx3:end]
        elseif (size(in),3) > 1
            out[:,1:incx1:end,1:incx2:end] = in[:,1:incx1:end,1:incx2:end]
        elseif (size(in),2) > 1
            out[:,1:incx1:end] = in[:,1:incx1:end]
        end

    end

    return out
end


text:
const UF = CLIMAParameters.SurfaceFluxes.UniversalFunctions

UF.Pr_0_Businger(::AbstractEarthParameterSet) = 0.74
UF.a_m_Businger(::AbstractEarthParameterSet) = 4.7
UF.a_h_Businger(::AbstractEarthParameterSet) = 4.7

UF.Pr_0_Gryanik(::AbstractEarthParameterSet) = 0.98
UF.a_m_Gryanik(::AbstractEarthParameterSet) = 5.0
UF.a_h_Gryanik(::AbstractEarthParameterSet) = 5.0
UF.b_m_Gryanik(::AbstractEarthParameterSet) = 0.3
UF.b_h_Gryanik(::AbstractEarthParameterSet) = 0.4

UF.Pr_0_Grachev(::AbstractEarthParameterSet) = 0.98
UF.a_m_Grachev(::AbstractEarthParameterSet) = 5.0
UF.a_h_Grachev(::AbstractEarthParameterSet) = 5.0
UF.b_m_Grachev(ps::AbstractEarthParameterSet) = UF.a_m_Grachev(ps) / 6.5
UF.b_h_Grachev(::AbstractEarthParameterSet) = 5.0
UF.c_h_Grachev(::AbstractEarthParameterSet) = 3.0

ext: jl

max_stars_repo_name: climate-machine/CLIMAParameters

max_issues_repo_name: climate-machine/CLIMAParameters

max_forks_repo_name: climate-machine/CLIMAParameters


text:
function part1(input)
    n = parse(Int, readchomp(input))
    houses = zeros(Int, n)
    for i = 1:n
        houses[i:i:end] .+= 10 * i
    end
    return findfirst(>=(n), houses)
end

function part2(input)
    n = parse(Int, readchomp(input))
    houses = zeros(Int, n)
    for i = 1:n
        houses[i:i:min(end, (50 * i))] .+= 11 * i
    end
    return findfirst(>=(n), houses)
end


text:
"""
    $SIGNATURES

Return a random sample as a tangent vector from a belief represented by coordinates on a manifold at point p.

Notes

"""
function sampleTangent end

# Sampling MKD
function sampleTangent(M::AbstractGroupManifold, x::ManifoldKernelDensity, p=mean(x))
    # get legacy matrix of coordinates and selected labels
    coords, lbls = sample(x.belief,1)
    X = hat(x.manifold, p, coords)
    return X
end

function sampleTangent(x::ManifoldKernelDensity, p=mean(x))
  return sampleTangent(x.manifold, x, p)
end

# Sampling Distributions
function sampleTangent(M::AbstractManifold, z::Distribution, p, basis::AbstractBasis)
  return get_vector(M, p, rand(z), basis)
end

function sampleTangent(M::AbstractGroupManifold, z::Distribution, p=identity_element(M))
  return hat(M, p, rand(z,1)[:]) #TODO find something better than (z,1)[:]
end

"""
    $SIGNATURES

Return a random sample point on a manifold from a belief represented by coordinates at point p.

Notes

"""
function samplePoint(M::AbstractManifold, sbelief, p, basis::AbstractBasis, retraction_method::AbstractRetractionMethod=ExponentialRetraction())
  X = sampleTangent(M, sbelief, p, basis)
  return retract(M, p, X, retraction_method)
end
function samplePoint(M::AbstractGroupManifold, sbelief, p=identity_element(M), retraction_method::AbstractRetractionMethod=ExponentialRetraction())
  X = sampleTangent(M, sbelief, p)
  return retract(M, p, X, retraction_method)
end

function samplePoint(M::AbstractGroupManifold, sbelief::ManifoldKernelDensity, p=identity_element(M, mean(sbelief)), retraction_method::AbstractRetractionMethod=ExponentialRetraction())
  X = sampleTangent(M, sbelief, p)
  return retract(M, p, X, retraction_method)
end

function samplePoint(x::ManifoldKernelDensity, p=mean(x))
  return samplePoint(x.manifold, x, p)
end

# FIXME: rather use manifolds
function samplePoint(distr::SamplableBelief)
  Base.depwarn("samplePoint(distr::SamplableBelief) should be replaced by samplePoint(M<:AbstractManifold, distr::SamplableBelief, ...)", :samplePoint)
  rand(distr,1)
end

## default getSample
"""
    $SIGNATURES

Sample the factor in `CalcFactor`. A default `getSample` method is provided that should cover most use cases,
if more advanced sampling is required, the `getSample` function should be extended.

The default behavior for `getSample` is as follows:
- The `SamplableBelief`` shall be in the field `Z` and that shall be enough to fully define the factor, i.e. `Z<:SamplableBelief` should be the only field.
- Sampling on `<:AbstractManifoldMinimize` factors defined on Group Manifolds:
  - `getSample` normally returns a tangent vector at the identity element, however it should just match the custom factor definition.
- Sampling on prior (`<:AbstractPrior`) factors :
  - `getSample` must return a point on the manifold that matches the point representation of the variable.

Notes
- Users should overload this method should their factor not only use field `Z` for the `SamplableBelief`.
- See the Custom Factors section in the Caesar.jl documentation for more examples and details.
- Also see issue https://github.com/JuliaRobotics/IncrementalInference.jl/issues/1441

See also: [`getMeasurementParametric`](@ref)
"""
function getSample end

function getSample(cf::CalcFactor{<:AbstractPrior})
  M = getManifold(cf.factor)
  if hasfield(typeof(cf.factor), :Z)
    X = samplePoint(M, cf.factor.Z)
  else
    error("""Factor $(typeof(cf.factor)) does not have a field `Z`, to use the default `getSample` method, use `Z` for the measurement.
              Alternatively, provide a `getSample` method. See IIF issue #1441 and Custom Factors in the Caesar documentation.""")
  end
  return X
end

function getSample(cf::CalcFactor{<:AbstractRelative})
  M = getManifold(cf.factor)
  if hasfield(typeof(cf.factor), :Z)
    X = sampleTangent(M, cf.factor.Z)
  else
    error("""Factor $(typeof(cf.factor)) does not have a field `Z`, to use the default `getSample` method, use `Z` for the measurement.
              Alternatively, provide a `getSample` method. See IIF issue #1441 and Custom Factors in the Caesar documentation.""")
  end
  return X
end


text:
"""
    PrintSourceCode(PackageName::T; RootDir::T = pwd()) where T<:AbstractString

Convert `PackageName.jl` file in `src` and all `include()`ed files of package `PackageName` to a tex file: `PackageName.tex`
and store it in `RootDir`.

# Example
```jldoctest
PrintSourceCode("LinearAlgebra")
```
"""
function PrintSourceCode(
    PackageName::T;
    RootDir::T = pwd()
) where {T<:AbstractString}

    ## 1. Retrieve package directory
    if PackageName == "Base"
        PackagePath = normpath(joinpath(Base.Sys.BINDIR, "..", "share/julia/base/"))
    else
        PackagePath = dirname(Base.find_package(PackageName))
    end

    ## 2. Retrieve file name and file version
    if occursin("base", PackagePath)
        FileName = "Base"
        FileVersion = string(VERSION)
    elseif occursin("stdlib", PackagePath)
        FileName    = basename(dirname(PackagePath))
        FileVersion = basename(dirname(dirname(PackagePath)))
    else
        FileName    = basename(dirname(dirname(PackagePath)))
        FileVersion = basename(dirname(PackagePath))
    end

    ## 4. Retrieve Filelist from file "PackageName.jl"
    FileList = String[PackageName * ".jl"]
    for ln in readlines(joinpath(PackagePath, PackageName * ".jl"))
        if occursin(r"include\(\"(.*?)\"\)$", ln)
            file = replace(ln, r"include\(\"(.*?)\"\)$" => s"\1")
            push!(FileList, file)
        end
    end

    ## 5. Write the file list to a tex file
    pre = """
    \\documentclass[11pt]{article}
    \\usepackage{/Users/likanzhan/Documents/Meta/Technique/LaTex/Preamble_Article}
    \\begin{document}
    \\tableofcontents
    \\newpage
    """
    post = "\\end{document}"
    TexFile = joinpath(RootDir, FileName * ".tex")

    open(TexFile, "w") do io
        println(io, pre)
        for file in FileList
            fileRep = replace(file, "_" => "\\_")
            println(io, "\\section{", joinpath(fileRep), "}")
            println(io, "\\inputminted{julia}{", joinpath(PackagePath, file), "}")
            println(io, "\\newpage", "\n")
        end
        println(io, "\n", post)
        println(io, "%%%%%% PkgV: " * FileVersion)
    end

    @info """

    - Source Directory: `$PackagePath`
    - Tex File: `$TexFile`
    - FileName: `$FileName`;
    - FileVersion: `$FileVersion`

    """
end


text:
using ProgressMeter
using .Threads
using Lasso
using NPZ
using Statistics
using PyCall
using SparseArrays

println(ARGS[1])
println("Threads = $(nthreads())")

py"""
import scipy.sparse as sp

def spsave(x, name):
    return sp.save_npz(f"Combined_Results/combined_results_{name}.npz", x)
"""

PyObject(S::SparseMatrixCSC) =
    pyimport("scipy.sparse").csc_matrix((S.nzval, S.rowval .- 1, S.colptr .- 1), shape=size(S))

struct InvResult{T,S}
    i::Int
    λmin::T
    dn::T
    m::S
end

function fit_model(G,dv,i)
    dn = std(dv)
    dv ./= dn
    path = fit(LassoPath, G, dv, λminratio = 1e-3, standardize=false, intercept=false, maxncoef=size(G,2))
    λmin = path.λ[argmin(aicc(path))]
    m = coef(path; select=MinAICc())
    return InvResult(i, λmin, dn, m)
end

G = npzread("tmp/G.npy")
d = npzread("tmp/data.npy")
p = Progress(size(d,2), dt=0.5, barglyphs=BarGlyphs("[=> ]"), barlen=50, color=:yellow)

results = Array{Any,1}(undef, size(d,2))

@threads for i = 1:size(d,2)
    results[i] = fit_model(G,d[:,i],i)
    next!(p)
end

sres = reduce(hcat, [x.m.*x.dn for x in results])
py"spsave"(PyObject(sres),ARGS[1])

max_stars_repo_name: jbmuir/DAS-Reconstruction

max_issues_repo_name: jbmuir/DAS-Reconstruction

max_forks_repo_name: jbmuir/DAS-Reconstruction


text:
# These are methods to implement particle segregation - separating particles that touch for counting etc.

struct Watershed
    sz::Tuple
    steps:: LinRange
    bins::Vector{Vector{CartesianIndex}}

    function Watershed(data::Array, steps::LinRange)
        indexof(x,lr::LinRange)=min(lr.len,max(1,trunc(Int,(lr.lendiv*(x-lr.start))/(lr.stop-lr.start))+1))
        bins = [ CartesianIndex[] for i in 1:length(steps) ]
        foreach(ci->push!(bins[indexof(data[ci], steps)], ci), CartesianIndices(data))
        return new(size(data), steps, bins)
    end
end

bin(ws::Watershed, idx::Int) = ( ws.steps[idx-1], ws.steps[idx] )
NeXLSpectrum.depth(ws::Watershed) = length(ws.steps)

function masks(ws::Watershed)::BitArray
    prev = fill(true, ws.sz)
    res = BitArray(undef, (ws.sz..., length(ws.steps)))
    for (i, vci) in enumerate(ws.bins)
        foreach(ci->prev[ci]=false,vci)
        res[:,:,i]=prev
    end
    return res
end

function areas(ws::Watershed)::Vector{Int}
    res = [size(vci,1) for vci in reverse(ws.bins)]
    sum = 0
    for i in eachindex(res)
        res[i] = (sum += res[i])
    end
    return res
end

max_stars_repo_name: NicholasWMRitchie/NeXLParticle

max_issues_repo_name: NicholasWMRitchie/NeXLParticle

max_forks_repo_name: NicholasWMRitchie/NeXLParticle


text:
struct Derivative{T<:Union{Int,Tuple{Vararg{Int}}}}
    order :: T
    function Derivative{T}(order::T) where {T<:Int}
        order < 0 && return throw(DomainError(order, "Derivative is only defined for positive integers"))
        return new{T}(order)
    end
    function Derivative{T}(order::T) where {T<:Tuple{Vararg{Int}}}
        any(n -> n < 0, order) && return throw(DomainError(order, "Derivative is only defined for positive integers"))
        return new{T}(order)
    end
end
Derivative(order::T) where {T<:Union{Int,Tuple{Vararg{Int}}}} = Derivative{T}(order)

struct Integral{T<:Union{Int,Tuple{Vararg{Int}}}}
    order :: T
    function Integral{T}(order::T) where {T<:Int}
        order < 0 && return throw(DomainError(order, "Integral is only defined for positive integers"))
        return new{T}(order)
    end
    function Integral{T}(order::T) where {T<:Tuple{Vararg{Int}}}
        any(n -> n < 0, order) && return throw(DomainError(order, "Integral is only defined for positive integers"))
        return new{T}(order)
    end
end
Integral(order::T) where {T<:Union{Int,Tuple{Vararg{Int}}}} = Integral{T}(order)

for (F, f, f!) ∈ ((:Derivative, :differentiate, :differentiate!), (:Integral, :integrate, :integrate!))
    @eval begin
        Base.:*(ℱ₁::$F{Int}, ℱ₂::$F{Int}) = $F(ℱ₁.order + ℱ₂.order)
        Base.:*(ℱ₁::$F{NTuple{N,Int}}, ℱ₂::$F{NTuple{N,Int}}) where {N} = $F(map(+, ℱ₁.order, ℱ₂.order))

        Base.:^(ℱ::$F{Int}, n::Int) = $F(ℱ.order * n)
        Base.:^(ℱ::$F{NTuple{N,Int}}, n::Int) where {N} = $F(map(αᵢ -> *(αᵢ, n), ℱ.order))
        Base.:^(ℱ::$F{NTuple{N,Int}}, α::NTuple{N,Int}) where {N} = $F(map(*, ℱ.order, α))

        (ℱ::$F)(a::Sequence) = *(ℱ, a)
        Base.:*(ℱ::$F, a::Sequence) = $f(a, ℱ.order)

        function $f(a::Sequence, α=1)
            ℱ = $F(α)
            space_a = space(a)
            new_space = image(ℱ, space_a)
            CoefType = _coeftype(ℱ, space_a, eltype(a))
            c = Sequence(new_space, Vector{CoefType}(undef, dimension(new_space)))
            _apply!(c, ℱ, a)
            return c
        end

        function $f!(c::Sequence, a::Sequence, α=1)
            ℱ = $F(α)
            space(c) == image(ℱ, space(a)) || return throw(DomainError)
            _apply!(c, ℱ, a)
            return c
        end

        function project(ℱ::$F, domain::VectorSpace, codomain::VectorSpace, ::Type{T}) where {T}
            _iscompatible(domain, codomain) || return throw(DimensionMismatch)
            ind_domain = _findposition_nzind_domain(ℱ, domain, codomain)
            ind_codomain = _findposition_nzind_codomain(ℱ, domain, codomain)
            C = LinearOperator(domain, codomain, sparse(ind_codomain, ind_domain, zeros(T, length(ind_domain)), dimension(codomain), dimension(domain)))
            _project!(C, ℱ)
            return C
        end

        function project!(C::LinearOperator, ℱ::$F)
            _iscompatible(domain(C), codomain(C)) || return throw(DimensionMismatch)
            coefficients(C) .= zero(eltype(C))
            _project!(C, ℱ)
            return C
        end

        _findposition_nzind_domain(ℱ::$F, domain, codomain) =
            _findposition(_nzind_domain(ℱ, domain, codomain), domain)

        _findposition_nzind_codomain(ℱ::$F, domain, codomain) =
            _findposition(_nzind_codomain(ℱ, domain, codomain), codomain)
    end
end

# Sequence spaces

for F ∈ (:Derivative, :Integral)
    @eval begin
        image(ℱ::$F{NTuple{N,Int}}, s::TensorSpace{<:NTuple{N,BaseSpace}}) where {N} =
            TensorSpace(map((αᵢ, sᵢ) -> image($F(αᵢ), sᵢ), ℱ.order, spaces(s)))

        _coeftype(ℱ::$F, s::TensorSpace, ::Type{T}) where {T} =
            @inbounds promote_type(_coeftype($F(ℱ.order[1]), s[1], T), _coeftype($F(Base.tail(ℱ.order)), Base.tail(s), T))
        _coeftype(ℱ::$F, s::TensorSpace{<:Tuple{BaseSpace}}, ::Type{T}) where {T} =
            @inbounds _coeftype($F(ℱ.order[1]), s[1], T)

        function _apply!(c::Sequence{<:TensorSpace}, ℱ::$F, a)
            space_a = space(a)
            A = _no_alloc_reshape(coefficients(a), dimensions(space_a))
            C = _no_alloc_reshape(coefficients(c), dimensions(space(c)))
            _apply!(C, ℱ, space_a, A)
            return c
        end

        _apply!(C, ℱ::$F, space::TensorSpace, A) =
            @inbounds _apply!(C, $F(ℱ.order[1]), space[1], _apply($F(Base.tail(ℱ.order)), Base.tail(space), A))

        _apply!(C, ℱ::$F, space::TensorSpace{<:Tuple{BaseSpace}}, A) =
            @inbounds _apply!(C, $F(ℱ.order[1]), space[1], A)

        _apply(ℱ::$F, space::TensorSpace{<:NTuple{N₁,BaseSpace}}, A::AbstractArray{T,N₂}) where {N₁,T,N₂} =
            @inbounds _apply($F(ℱ.order[1]), space[1], Val(N₂-N₁+1), _apply($F(Base.tail(ℱ.order)), Base.tail(space), A))

        _apply(ℱ::$F, space::TensorSpace{<:Tuple{BaseSpace}}, A::AbstractArray{T,N}) where {T,N} =
            @inbounds _apply($F(ℱ.order[1]), space[1], Val(N), A)
    end
end

for F ∈ (:Derivative, :Integral)
    for (_f, __f) ∈ ((:_nzind_domain, :__nzind_domain), (:_nzind_codomain, :__nzind_codomain))
        @eval begin
            $_f(ℱ::$F{NTuple{N,Int}}, domain::TensorSpace{<:NTuple{N,BaseSpace}}, codomain::TensorSpace{<:NTuple{N,BaseSpace}}) where {N} =
                Base.Iterators.ProductIterator($__f(ℱ, domain, codomain))
            $__f(ℱ::$F, domain::TensorSpace, codomain) =
                @inbounds ($_f($F(ℱ.order[1]), domain[1], codomain[1]), $__f($F(Base.tail(ℱ.order)), Base.tail(domain), Base.tail(codomain))...)
            $__f(ℱ::$F, domain::TensorSpace{<:Tuple{BaseSpace}}, codomain) =
                @inbounds ($_f($F(ℱ.order[1]), domain[1], codomain[1]),)
        end
    end

    @eval function _project!(C::LinearOperator{<:SequenceSpace,<:SequenceSpace}, ℱ::$F)
        domain_C = domain(C)
        codomain_C = codomain(C)
        CoefType = eltype(C)
        @inbounds for (α, β) ∈ zip(_nzind_codomain(ℱ, domain_C, codomain_C), _nzind_domain(ℱ, domain_C, codomain_C))
            C[α,β] = _nzval(ℱ, domain_C, codomain_C, CoefType, α, β)
        end
        return C
    end
end

@generated function _nzval(𝒟::Derivative{NTuple{N,Int}}, domain::TensorSpace{<:NTuple{N,BaseSpace}}, codomain::TensorSpace{<:NTuple{N,BaseSpace}}, ::Type{T}, α, β) where {N,T}
    p = :(_nzval(Derivative(𝒟.order[1]), domain[1], codomain[1], T, α[1], β[1]))
    for i ∈ 2:N
        p = :(_nzval(Derivative(𝒟.order[$i]), domain[$i], codomain[$i], T, α[$i], β[$i]) * $p)
    end
    return p
end

@generated function _nzval(ℐ::Integral{NTuple{N,Int}}, domain::TensorSpace{<:NTuple{N,BaseSpace}}, codomain::TensorSpace{<:NTuple{N,BaseSpace}}, ::Type{T}, α, β) where {N,T}
    p = :(_nzval(Integral(ℐ.order[1]), domain[1], codomain[1], T, α[1], β[1]))
    for i ∈ 2:N
        p = :(_nzval(Integral(ℐ.order[$i]), domain[$i], codomain[$i], T, α[$i], β[$i]) * $p)
    end
    return p
end

# Taylor

image(𝒟::Derivative, s::Taylor) = Taylor(max(0, order(s)-𝒟.order))

_coeftype(::Derivative, ::Taylor, ::Type{T}) where {T} = typeof(zero(T)*0)

function _apply!(c::Sequence{Taylor}, 𝒟::Derivative, a)
    n = 𝒟.order
    if n == 0
        coefficients(c) .= coefficients(a)
    else
        order_a = order(a)
        if order_a < n
            @inbounds c[0] = zero(eltype(c))
        elseif n == 1
            @inbounds for i ∈ 1:order_a
                c[i-1] = i * a[i]
            end
        else
            space_a = space(a)
            CoefType_a = eltype(a)
            @inbounds for i ∈ n:order_a
                c[i-n] = _nzval(𝒟, space_a, space_a, CoefType_a, i-n, i) * a[i]
            end
        end
    end
    return c
end

function _apply!(C::AbstractArray{T}, 𝒟::Derivative, space::Taylor, A) where {T}
    n = 𝒟.order
    if n == 0
        C .= A
    else
        ord = order(space)
        if ord < n
            C .= zero(T)
        elseif n == 1
            @inbounds for i ∈ 1:ord
                selectdim(C, 1, i) .= i .* selectdim(A, 1, i+1)
            end
        else
            CoefType_A = eltype(A)
            @inbounds for i ∈ n:ord
                selectdim(C, 1, i-n+1) .= _nzval(𝒟, space, space, CoefType_A, i-n, i) .* selectdim(A, 1, i+1)
            end
        end
    end
    return C
end

function _apply(𝒟::Derivative, space::Taylor, ::Val{D}, A::AbstractArray{T,N}) where {D,T,N}
    n = 𝒟.order
    CoefType = _coeftype(𝒟, space, T)
    if n == 0
        return convert(Array{CoefType,N}, A)
    else
        ord = order(space)
        if ord < n
            return zeros(CoefType, ntuple(i -> ifelse(i == D, 1, size(A, i)), Val(N)))
        elseif n == 1
            C = Array{CoefType,N}(undef, ntuple(i -> ifelse(i == D, ord, size(A, i)), Val(N)))
            @inbounds for i ∈ 1:ord
                selectdim(C, D, i) .= i .* selectdim(A, D, i+1)
            end
            return C
        else
            C = Array{CoefType,N}(undef, ntuple(i -> ifelse(i == D, ord-n+1, size(A, i)), Val(N)))
            @inbounds for i ∈ n:ord
                selectdim(C, D, i-n+1) .= _nzval(𝒟, space, space, T, i-n, i) .* selectdim(A, D, i+1)
            end
            return C
        end
    end
end

_nzind_domain(𝒟::Derivative, domain::Taylor, codomain::Taylor) =
    𝒟.order:min(order(domain), order(codomain)+𝒟.order)

_nzind_codomain(𝒟::Derivative, domain::Taylor, codomain::Taylor) =
    0:min(order(domain)-𝒟.order, order(codomain))

function _nzval(𝒟::Derivative, ::Taylor, ::Taylor, ::Type{T}, i, j) where {T}
    n = 𝒟.order
    p = one(T)*1
    for k ∈ 1:n
        p *= i+k
    end
    return convert(T, p)
end

image(ℐ::Integral, s::Taylor) = Taylor(order(s)+ℐ.order)

_coeftype(::Integral, ::Taylor, ::Type{T}) where {T} = typeof(inv(one(T)*1)*zero(T))

function _apply!(c::Sequence{Taylor}, ℐ::Integral, a)
    n = ℐ.order
    if n == 0
        coefficients(c) .= coefficients(a)
    elseif n == 1
        @inbounds c[0] = zero(eltype(c))
        @inbounds for i ∈ 0:order(a)
            c[i+1] = a[i] / (i+1)
        end
    else
        space_a = space(a)
        CoefType_a = eltype(a)
        @inbounds view(c, 0:n-1) .= zero(eltype(c))
        @inbounds for i ∈ 0:order(a)
            c[i+n] = _nzval(ℐ, space_a, space_a, CoefType_a, i+n, i) * a[i]
        end
    end
    return c
end

function _apply!(C::AbstractArray{T}, ℐ::Integral, space::Taylor, A) where {T}
    n = ℐ.order
    if n == 0
        C .= A
    elseif n == 1
        ord = order(space)
        @inbounds selectdim(C, 1, 1) .= zero(T)
        @inbounds for i ∈ 0:ord
            selectdim(C, 1, i+2) .= selectdim(A, 1, i+1) ./ (i+1)
        end
    else
        CoefType_A = eltype(A)
        ord = order(space)
        @inbounds selectdim(C, 1, 1:n) .= zero(T)
        @inbounds for i ∈ 0:ord
            selectdim(C, 1, i+n+1) .= _nzval(ℐ, space, space, CoefType_A, i+n, i) .* selectdim(A, 1, i+1)
        end
    end
    return C
end

function _apply(ℐ::Integral, space::Taylor, ::Val{D}, A::AbstractArray{T,N}) where {D,T,N}
    n = ℐ.order
    CoefType = _coeftype(ℐ, space, T)
    if n == 0
        return convert(Array{CoefType,N}, A)
    elseif n == 1
        ord = order(space)
        C = Array{CoefType,N}(undef, ntuple(i -> ifelse(i == D, ord+2, size(A, i)), Val(N)))
        @inbounds selectdim(C, D, 1) .= zero(CoefType)
        @inbounds for i ∈ 0:ord
            selectdim(C, D, i+2) .= selectdim(A, D, i+1) ./ (i+1)
        end
        return C
    else
        ord = order(space)
        C = Array{CoefType,N}(undef, ntuple(i -> ifelse(i == D, ord+n+1, size(A, i)), Val(N)))
        @inbounds selectdim(C, D, 1:n) .= zero(CoefType)
        @inbounds for i ∈ 0:ord
            selectdim(C, D, i+n+1) .= _nzval(ℐ, space, space, T, i+n, i) .* selectdim(A, D, i+1)
        end
        return C
    end
end

_nzind_domain(ℐ::Integral, domain::Taylor, codomain::Taylor) =
    0:min(order(domain), order(codomain)-ℐ.order)

_nzind_codomain(ℐ::Integral, domain::Taylor, codomain::Taylor) =
    ℐ.order:min(order(domain)+ℐ.order, order(codomain))

_nzval(ℐ::Integral, s₁::Taylor, s₂::Taylor, ::Type{T}, i, j) where {T} =
    convert(T, inv(_nzval(Derivative(ℐ.order), s₁, s₂, T, j, i)))

# Fourier

image(::Derivative, s::Fourier) = s

_coeftype(::Derivative, ::Fourier{T}, ::Type{S}) where {T,S} = complex(typeof(zero(T)*0*zero(S)))

function _apply!(c::Sequence{<:Fourier}, 𝒟::Derivative, a)
    n = 𝒟.order
    if n == 0
        coefficients(c) .= coefficients(a)
    else
        ω = one(eltype(a))*frequency(a)
        @inbounds c[0] = zero(eltype(c))
        if n == 1
            iω = im*ω
            @inbounds for j ∈ 1:order(c)
                iωj = iω*j
                c[j] = iωj * a[j]
                c[-j] = -iωj * a[-j]
            end
        else
            if isodd(n)
                iⁿ = complex(0, ifelse(n%4 == 1, 1, -1))
                @inbounds for j ∈ 1:order(c)
                    iⁿωⁿjⁿ = iⁿ*(ω*j)^n
                    c[j] = iⁿωⁿjⁿ * a[j]
                    c[-j] = -iⁿωⁿjⁿ * a[-j]
                end
            else
                iⁿ_real = ifelse(n%4 == 0, 1, -1)
                @inbounds for j ∈ 1:order(c)
                    iⁿωⁿjⁿ_real = iⁿ_real*(ω*j)^n
                    c[j] = iⁿωⁿjⁿ_real * a[j]
                    c[-j] = iⁿωⁿjⁿ_real * a[-j]
                end
            end
        end
    end
    return c
end

function _apply!(C::AbstractArray{T}, 𝒟::Derivative, space::Fourier, A) where {T}
    n = 𝒟.order
    if n == 0
        C .= A
    else
        ord = order(space)
        ω = one(eltype(A))*frequency(space)
        @inbounds selectdim(C, 1, ord+1) .= zero(T)
        if n == 1
            iω = im*ω
            @inbounds for j ∈ 1:ord
                iωj = iω*j
                selectdim(C, 1, ord+1+j) .= iωj .* selectdim(A, 1, ord+1+j)
                selectdim(C, 1, ord+1-j) .= -iωj .* selectdim(A, 1, ord+1-j)
            end
        else
            if isodd(n)
                iⁿ = complex(0, ifelse(n%4 == 1, 1, -1))
                @inbounds for j ∈ 1:ord
                    iⁿωⁿjⁿ = iⁿ*(ω*j)^n
                    selectdim(C, 1, ord+1+j) .= iⁿωⁿjⁿ .* selectdim(A, 1, ord+1+j)
                    selectdim(C, 1, ord+1-j) .= -iⁿωⁿjⁿ .* selectdim(A, 1, ord+1-j)
                end
            else
                iⁿ_real = ifelse(n%4 == 0, 1, -1)
                @inbounds for j ∈ 1:ord
                    iⁿωⁿjⁿ_real = iⁿ_real*(ω*j)^n
                    selectdim(C, 1, ord+1+j) .= iⁿωⁿjⁿ_real .* selectdim(A, 1, ord+1+j)
                    selectdim(C, 1, ord+1-j) .= iⁿωⁿjⁿ_real .* selectdim(A, 1, ord+1-j)
                end
            end
        end
    end
    return C
end

function _apply(𝒟::Derivative, space::Fourier, ::Val{D}, A::AbstractArray{T,N}) where {D,T,N}
    n = 𝒟.order
    CoefType = _coeftype(𝒟, space, T)
    if n == 0
        return convert(Array{CoefType,N}, A)
    else
        C = Array{CoefType,N}(undef, size(A))
        ord = order(space)
        ω = one(T)*frequency(space)
        @inbounds selectdim(C, D, ord+1) .= zero(CoefType)
        if n == 1
            iω = im*ω
            @inbounds for j ∈ 1:ord
                iωj = iω*j
                selectdim(C, D, ord+1+j) .= iωj .* selectdim(A, D, ord+1+j)
                selectdim(C, D, ord+1-j) .= -iωj .* selectdim(A, D, ord+1-j)
            end
        else
            if isodd(n)
                iⁿ = complex(0, ifelse(n%4 == 1, 1, -1))
                @inbounds for j ∈ 1:ord
                    iⁿωⁿjⁿ = iⁿ*(ω*j)^n
                    selectdim(C, D, ord+1+j) .= iⁿωⁿjⁿ .* selectdim(A, D, ord+1+j)
                    selectdim(C, D, ord+1-j) .= -iⁿωⁿjⁿ .* selectdim(A, D, ord+1-j)
                end
            else
                iⁿ_real = ifelse(n%4 == 0, 1, -1)
                @inbounds for j ∈ 1:ord
                    iⁿωⁿjⁿ_real = iⁿ_real*(ω*j)^n
                    selectdim(C, D, ord+1+j) .= iⁿωⁿjⁿ_real .* selectdim(A, D, ord+1+j)
                    selectdim(C, D, ord+1-j) .= iⁿωⁿjⁿ_real .* selectdim(A, D, ord+1-j)
                end
            end
        end
        return C
    end
end

function _nzind_domain(::Derivative, domain::Fourier, codomain::Fourier)
    frequency(domain) == frequency(codomain) || return throw(DomainError)
    ord = min(order(domain), order(codomain))
    return -ord:ord
end

function _nzind_codomain(::Derivative, domain::Fourier, codomain::Fourier)
    frequency(domain) == frequency(codomain) || return throw(DomainError)
    ord = min(order(domain), order(codomain))
    return -ord:ord
end

function _nzval(𝒟::Derivative, domain::Fourier, ::Fourier, ::Type{T}, i, j) where {T}
    n = 𝒟.order
    if n == 0
        return one(T)
    else
        ωⁿjⁿ = (one(T)*frequency(domain)*j)^n
        r = n % 4
        if r == 0
            return convert(T, ωⁿjⁿ)
        elseif r == 1
            return convert(T, im*ωⁿjⁿ)
        elseif r == 2
            return convert(T, -ωⁿjⁿ)
        else
            return convert(T, -im*ωⁿjⁿ)
        end
    end
end

image(::Integral, s::Fourier) = s

_coeftype(::Integral, ::Fourier{T}, ::Type{S}) where {T,S} = complex(typeof(inv(one(S)*one(T)*1)*zero(S)))

function _apply!(c::Sequence{<:Fourier}, ℐ::Integral, a)
    n = ℐ.order
    if n == 0
        coefficients(c) .= coefficients(a)
    else
        @inbounds iszero(a[0]) || return throw(DomainError)
        ω = one(eltype(a))*frequency(a)
        @inbounds c[0] = zero(eltype(c))
        if n == 1
            @inbounds for j ∈ 1:order(c)
                iω⁻¹j⁻¹ = im*inv(ω*j)
                c[j] = -iω⁻¹j⁻¹ * a[j]
                c[-j] = iω⁻¹j⁻¹ * a[-j]
            end
        else
            if isodd(n)
                iⁿ = complex(0, ifelse(n%4 == 1, 1, -1))
                @inbounds for j ∈ 1:order(c)
                    iⁿω⁻ⁿj⁻ⁿ = iⁿ*inv(ω*j)^n
                    c[j] = -iⁿω⁻ⁿj⁻ⁿ * a[j]
                    c[-j] = iⁿω⁻ⁿj⁻ⁿ * a[-j]
                end
            else
                iⁿ_real = ifelse(n%4 == 0, 1, -1)
                @inbounds for j ∈ 1:order(c)
                    iⁿω⁻ⁿj⁻ⁿ_real = iⁿ_real*inv(ω*j)^n
                    c[j] = iⁿω⁻ⁿj⁻ⁿ_real * a[j]
                    c[-j] = iⁿω⁻ⁿj⁻ⁿ_real * a[-j]
                end
            end
        end
    end
    return c
end

function _apply!(C::AbstractArray{T}, ℐ::Integral, space::Fourier, A) where {T}
    n = ℐ.order
    if n == 0
        C .= A
    else
        ord = order(space)
        @inbounds iszero(selectdim(A, 1, ord+1)) || return throw(DomainError)
        ω = one(eltype(A))*frequency(space)
        @inbounds selectdim(C, 1, ord+1) .= zero(T)
        if n == 1
            @inbounds for j ∈ 1:ord
                iω⁻¹j⁻¹ = im*inv(ω*j)
                selectdim(C, 1, ord+1+j) .= -iω⁻¹j⁻¹ .* selectdim(A, 1, ord+1+j)
                selectdim(C, 1, ord+1-j) .= iω⁻¹j⁻¹ .* selectdim(A, 1, ord+1-j)
            end
        else
            if isodd(n)
                iⁿ = complex(0, ifelse(n%4 == 1, 1, -1))
                @inbounds for j ∈ 1:ord
                    iⁿω⁻ⁿj⁻ⁿ = iⁿ*inv(ω*j)^n
                    selectdim(C, 1, ord+1+j) .= -iⁿω⁻ⁿj⁻ⁿ .* selectdim(A, 1, ord+1+j)
                    selectdim(C, 1, ord+1-j) .= iⁿω⁻ⁿj⁻ⁿ .* selectdim(A, 1, ord+1-j)
                end
            else
                iⁿ_real = ifelse(n%4 == 0, 1, -1)
                @inbounds for j ∈ 1:ord
                    iⁿω⁻ⁿj⁻ⁿ_real = iⁿ_real*inv(ω*j)^n
                    selectdim(C, 1, ord+1+j) .= iⁿω⁻ⁿj⁻ⁿ_real .* selectdim(A, 1, ord+1+j)
                    selectdim(C, 1, ord+1-j) .= iⁿω⁻ⁿj⁻ⁿ_real .* selectdim(A, 1, ord+1-j)
                end
            end
        end
    end
    return C
end

function _apply(ℐ::Integral, space::Fourier, ::Val{D}, A::AbstractArray{T,N}) where {D,T,N}
    n = ℐ.order
    CoefType = _coeftype(ℐ, space, T)
    if n == 0
        return convert(Array{CoefType,N}, A)
    else
        ord = order(space)
        @inbounds iszero(selectdim(A, D, ord+1)) || return throw(DomainError)
        ω = one(T)*frequency(space)
        C = Array{CoefType,N}(undef, size(A))
        @inbounds selectdim(C, D, ord+1) .= zero(CoefType)
        if n == 1
            @inbounds for j ∈ 1:ord
                iω⁻¹j⁻¹ = im*inv(ω*j)
                selectdim(C, D, ord+1+j) .= -iω⁻¹j⁻¹ .* selectdim(A, D, ord+1+j)
                selectdim(C, D, ord+1-j) .= iω⁻¹j⁻¹ .* selectdim(A, D, ord+1-j)
            end
        else
            if isodd(n)
                iⁿ = complex(0, ifelse(n%4 == 1, 1, -1))
                @inbounds for j ∈ 1:ord
                    iⁿω⁻ⁿj⁻ⁿ = iⁿ*inv(ω*j)^n
                    selectdim(C, D, ord+1+j) .= -iⁿω⁻ⁿj⁻ⁿ .* selectdim(A, D, ord+1+j)
                    selectdim(C, D, ord+1-j) .= iⁿω⁻ⁿj⁻ⁿ .* selectdim(A, D, ord+1-j)
                end
            else
                iⁿ_real = ifelse(n%4 == 0, 1, -1)
                @inbounds for j ∈ 1:ord
                    iⁿω⁻ⁿj⁻ⁿ_real = iⁿ_real*inv(ω*j)^n
                    selectdim(C, D, ord+1+j) .= iⁿω⁻ⁿj⁻ⁿ_real .* selectdim(A, D, ord+1+j)
                    selectdim(C, D, ord+1-j) .= iⁿω⁻ⁿj⁻ⁿ_real .* selectdim(A, D, ord+1-j)
                end
            end
        end
        return C
    end
end

function _nzind_domain(::Integral, domain::Fourier, codomain::Fourier)
    frequency(domain) == frequency(codomain) || return throw(DomainError)
    ord = min(order(domain), order(codomain))
    return -ord:ord
end

function _nzind_codomain(::Integral, domain::Fourier, codomain::Fourier)
    frequency(domain) == frequency(codomain) || return throw(DomainError)
    ord = min(order(domain), order(codomain))
    return -ord:ord
end

function _nzval(ℐ::Integral, domain::Fourier, ::Fourier, ::Type{T}, i, j) where {T}
    n = ℐ.order
    if n == 0
        return one(T)
    else
        if j == 0
            return zero(T)
        else
            ω⁻ⁿj⁻ⁿ = inv(one(T)*frequency(domain)*j)^n
            r = n % 4
            if r == 0
                return convert(T, ω⁻ⁿj⁻ⁿ)
            elseif r == 1
                return convert(T, -im*ω⁻ⁿj⁻ⁿ)
            elseif r == 2
                return convert(T, -ω⁻ⁿj⁻ⁿ)
            else
                return convert(T, im*ω⁻ⁿj⁻ⁿ)
            end
        end
    end
end

# Chebyshev

image(𝒟::Derivative, s::Chebyshev) = Chebyshev(max(0, order(s)-𝒟.order))

_coeftype(::Derivative, ::Chebyshev, ::Type{T}) where {T} = typeof(zero(T)*0)

function _apply!(c::Sequence{Chebyshev}, 𝒟::Derivative, a)
    n = 𝒟.order
    if n == 0
        coefficients(c) .= coefficients(a)
    elseif n == 1
        CoefType = eltype(c)
        order_a = order(a)
        if order_a < n
            @inbounds c[0] = zero(CoefType)
        else
            @inbounds for i ∈ 0:order_a-1
                c[i] = zero(CoefType)
                @inbounds for j ∈ i+1:2:order_a
                    c[i] += j * a[j]
                end
                c[i] *= 2
            end
        end
    else # TODO: lift restriction
        return throw(DomainError)
    end
    return c
end

function _apply!(C::AbstractArray{T}, 𝒟::Derivative, space::Chebyshev, A) where {T}
    n = 𝒟.order
    if n == 0
        C .= A
    elseif n == 1
        ord = order(space)
        if ord < n
            C .= zero(T)
        else
            @inbounds for i ∈ 0:ord-1
                Cᵢ = selectdim(C, 1, i+1)
                Cᵢ .= zero(T)
                @inbounds for j ∈ i+1:2:ord
                    Cᵢ .+= (2j) .* selectdim(A, 1, j+1)
                end
            end
        end
    else # TODO: lift restriction
        return throw(DomainError)
    end
    return C
end

function _apply(𝒟::Derivative, space::Chebyshev, ::Val{D}, A::AbstractArray{T,N}) where {D,T,N}
    n = 𝒟.order
    CoefType = _coeftype(𝒟, space, T)
    if n == 0
        return convert(Array{CoefType,N}, A)
    elseif n == 1
        ord = order(space)
        if ord < n
            return zeros(CoefType, ntuple(i -> i == D ? 1 : size(A, i), Val(N)))
        else
            C = zeros(CoefType, ntuple(i -> i == D ? ord : size(A, i), Val(N)))
            @inbounds for i ∈ 0:ord-1
                Cᵢ = selectdim(C, D, i+1)
                @inbounds for j ∈ i+1:2:ord
                    Cᵢ .+= (2j) .* selectdim(A, D, j+1)
                end
            end
            return C
        end
    else # TODO: lift restriction
        return throw(DomainError)
    end
end

function _nzind_domain(𝒟::Derivative, domain::Chebyshev, codomain::Chebyshev)
    if 𝒟.order == 0
        return collect(0:min(order(domain), order(codomain)))
    elseif 𝒟.order == 1
        len = sum(j -> length((j-1)%2:2:min(j-1, order(codomain))), 1:order(domain))
        v = Vector{Int}(undef, len)
        l = 0
        @inbounds for j ∈ 1:order(domain)
            lnext = l+length((j-1)%2:2:min(j-1, order(codomain)))
            view(v, 1+l:lnext) .= j
            l = lnext
        end
        return v
    else # TODO: lift restriction
        return throw(DomainError)
    end
end

function _nzind_codomain(𝒟::Derivative, domain::Chebyshev, codomain::Chebyshev)
    if 𝒟.order == 0
        return collect(0:min(order(domain), order(codomain)))
    elseif 𝒟.order == 1
        len = sum(j -> length((j-1)%2:2:min(j-1, order(codomain))), 1:order(domain))
        v = Vector{Int}(undef, len)
        l = 0
        @inbounds for j ∈ 1:order(domain)
            r = (j-1)%2:2:min(j-1, order(codomain))
            lnext = l+length(r)
            view(v, 1+l:lnext) .= r
            l = lnext
        end
        return v
    else # TODO: lift restriction
        return throw(DomainError)
    end
end

function _nzval(𝒟::Derivative, ::Chebyshev, ::Chebyshev, ::Type{T}, i, j) where {T}
    n = 𝒟.order
    if n == 0
        return one(T)
    elseif n == 1
        return convert(T, 2j)
    else # TODO: lift restriction
        return throw(DomainError)
    end
end

image(ℐ::Integral, s::Chebyshev) = Chebyshev(order(s)+ℐ.order)

_coeftype(::Integral, ::Chebyshev, ::Type{T}) where {T} = typeof(zero(T)/1)

function _apply!(c::Sequence{Chebyshev}, ℐ::Integral, a)
    n = ℐ.order
    if n == 0
        coefficients(c) .= coefficients(a)
    elseif n == 1
        order_a = order(a)
        if order_a == 0
            @inbounds c[0] = a[0]
            @inbounds c[1] = a[1] / 2
        elseif order_a == 1
            @inbounds c[0] = a[0] - a[1] / 2
            @inbounds c[1] = a[0] / 2
            @inbounds c[2] = a[1] / 4
        else
            @inbounds c[0] = zero(eltype(c))
            @inbounds for i ∈ 2:2:order_a-1
                c[0] += a[i+1] / ((i+1)^2-1) - a[i] / (i^2-1)
            end
            if iseven(order_a)
                @inbounds c[0] -= a[order_a] / (order_a^2-1)
            end
            @inbounds c[0] = 2 * c[0] + a[0] - a[1] / 2
            @inbounds c[1] = (a[0] - a[2]) / 2
            @inbounds for i ∈ 2:order_a-1
                c[i] = (a[i-1] - a[i+1]) / (2i)
            end
            @inbounds c[order_a] = a[order_a-1] / (2order_a)
            @inbounds c[order_a+1] = a[order_a] / (2(order_a+1))
        end
    else # TODO: lift restriction
        return throw(DomainError)
    end
    return c
end

function _apply!(C::AbstractArray{T}, ℐ::Integral, space::Chebyshev, A) where {T}
    n = ℐ.order
    if n == 0
        C .= A
    elseif n == 1
        ord = order(space)
        @inbounds C₀ = selectdim(C, 1, 1)
        @inbounds C₁ = selectdim(C, 1, 2)
        @inbounds A₀ = selectdim(A, 1, 1)
        if ord == 0
            C₀ .= A₀
            C₁ .= A₀ ./ 2
        elseif ord == 1
            @inbounds A₁ = selectdim(A, 1, 2)
            C₀ .= A₀ .- A₁ ./ 2
            C₁ .= A₀ ./ 2
            @inbounds selectdim(C, 1, 3) .= A₁ ./ 4
        else
            C₀ .= zero(T)
            @inbounds for i ∈ 2:2:ord-1
                C₀ .+= selectdim(A, 1, i+2) ./ ((i+1)^2-1) .- selectdim(A, 1, i+1) ./ (i^2-1)
            end
            if iseven(ord)
                @inbounds C₀ .-= selectdim(A, 1, ord+1) ./ (ord^2-1)
            end
            @inbounds C₀ .= 2 .* C₀ .+ A₀ .- selectdim(A, 1, 2) ./ 2
            @inbounds C₁ .= (A₀ .- selectdim(A, 1, 3)) ./ 2
            @inbounds for i ∈ 2:ord-1
                selectdim(C, 1, i+1) .= (selectdim(A, 1, i) .- selectdim(A, 1, i+2)) ./ (2i)
            end
            @inbounds selectdim(C, 1, ord+1) .= selectdim(A, 1, ord) ./ (2ord)
            @inbounds selectdim(C, 1, ord+2) .= selectdim(A, 1, ord+1) ./ (2(ord+1))
        end
    else # TODO: lift restriction
        return throw(DomainError)
    end
    return C
end

function _apply(ℐ::Integral, space::Chebyshev, ::Val{D}, A::AbstractArray{T,N}) where {D,T,N}
    n = ℐ.order
    CoefType = _coeftype(ℐ, space, T)
    if n == 0
        return convert(Array{CoefType,N}, A)
    elseif n == 1
        ord = order(space)
        C = Array{CoefType,N}(undef, ntuple(i -> i == D ? ord+2 : size(A, i), Val(N)))
        @inbounds C₀ = selectdim(C, D, 1)
        @inbounds C₁ = selectdim(C, D, 2)
        @inbounds A₀ = selectdim(A, D, 1)
        if ord == 0
            C₀ .= A₀
            C₁ .= A₀ ./ 2
        elseif ord == 1
            @inbounds A₁ = selectdim(A, D, 2)
            C₀ .= A₀ .- A₁ ./ 2
            C₁ .= A₀ ./ 2
            @inbounds selectdim(C, D, 3) .= A₁ ./ 4
        else
            C₀ .= zero(CoefType)
            @inbounds for i ∈ 2:2:ord-1
                C₀ .+= selectdim(A, D, i+2) ./ ((i+1)^2-1) .- selectdim(A, D, i+1) ./ (i^2-1)
            end
            if iseven(ord)
                @inbounds C₀ .-= selectdim(A, D, ord+1) ./ (ord^2-1)
            end
            @inbounds C₀ .= 2 .* C₀ .+ A₀ .- selectdim(A, D, 2) ./ 2
            @inbounds C₁ .= (A₀ .- selectdim(A, D, 3)) ./ 2
            @inbounds for i ∈ 2:ord-1
                selectdim(C, D, i+1) .= (selectdim(A, D, i) .- selectdim(A, D, i+2)) ./ (2i)
            end
            @inbounds selectdim(C, D, ord+1) .= selectdim(A, D, ord) ./ (2ord)
            @inbounds selectdim(C, D, ord+2) .= selectdim(A, D, ord+1) ./ (2(ord+1))
        end
        return C
    else # TODO: lift restriction
        return throw(DomainError)
    end
end

function _nzind_domain(ℐ::Integral, domain::Chebyshev, codomain::Chebyshev)
    if ℐ.order == 0
        return collect(0:min(order(domain), order(codomain)))
    elseif ℐ.order == 1
        v = mapreduce(vcat, 0:order(domain)) do j
            if j < 2
                j+1 ≤ order(codomain) && return [j, j]
                return [j]
            else
                j+1 ≤ order(codomain) && return [j, j, j]
                j-1 ≤ order(codomain) && return [j, j]
                return [j]
            end
        end
        return v
    else # TODO: lift restriction
        return throw(DomainError)
    end
end

function _nzind_codomain(ℐ::Integral, domain::Chebyshev, codomain::Chebyshev)
    if ℐ.order == 0
        return collect(0:min(order(domain), order(codomain)))
    elseif ℐ.order == 1
        v = mapreduce(vcat, 0:order(domain)) do j
            if j < 2
                j+1 ≤ order(codomain) && return [0, j+1]
                return [0]
            else
                j+1 ≤ order(codomain) && return [0, j-1, j+1]
                j-1 ≤ order(codomain) && return [0, j-1]
                return [0]
            end
        end
        return v
    else # TODO: lift restriction
        return throw(DomainError)
    end
end

function _nzval(ℐ::Integral, ::Chebyshev, ::Chebyshev, ::Type{T}, i, j) where {T}
    n = ℐ.order
    if n == 0
        return one(T)
    elseif n == 1
        if i == 0
            if j == 0
                return one(T)
            elseif j == 1
                return convert(T, -one(T)/2)
            elseif iseven(j)
                return convert(T, 2one(T)/(1-j^2))
            else
                return convert(T, 2one(T)/(j^2-1))
            end
        elseif i == 1 && j == 0
            return convert(T, one(T)/2)
        elseif i == 2 && j == 1
            return convert(T, one(T)/4)
        else
            if i+1 == j
                return convert(T, -one(T)/(2i))
            else # i == j+1
                return convert(T, one(T)/(2i))
            end
        end
    else # TODO: lift restriction
        return throw(DomainError)
    end
end

# Cartesian spaces

for F ∈ (:Derivative, :Integral)
    @eval begin
        image(ℱ::$F, s::CartesianPower) =
            CartesianPower(image(ℱ, space(s)), nb_cartesian_product(s))

        image(ℱ::$F, s::CartesianProduct) =
            CartesianProduct(map(sᵢ -> image(ℱ, sᵢ), spaces(s)))

        _coeftype(ℱ::$F, s::CartesianPower, ::Type{T}) where {T} =
            _coeftype(ℱ, space(s), T)

        _coeftype(ℱ::$F, s::CartesianProduct, ::Type{T}) where {T} =
            @inbounds promote_type(_coeftype(ℱ, s[1], T), _coeftype(ℱ, Base.tail(s), T))
        _coeftype(ℱ::$F, s::CartesianProduct{<:Tuple{VectorSpace}}, ::Type{T}) where {T} =
            @inbounds _coeftype(ℱ, s[1], T)

        function _apply!(c::Sequence{<:CartesianPower}, ℱ::$F, a)
            @inbounds for i ∈ 1:nb_cartesian_product(space(c))
                _apply!(component(c, i), ℱ, component(a, i))
            end
            return c
        end
        function _apply!(c::Sequence{CartesianProduct{T}}, ℱ::$F, a) where {N,T<:NTuple{N,VectorSpace}}
            @inbounds _apply!(component(c, 1), ℱ, component(a, 1))
            @inbounds _apply!(component(c, 2:N), ℱ, component(a, 2:N))
            return c
        end
        function _apply!(c::Sequence{CartesianProduct{T}}, ℱ::$F, a) where {T<:Tuple{VectorSpace}}
            @inbounds _apply!(component(c, 1), ℱ, component(a, 1))
            return c
        end

        function _findposition_nzind_domain(ℱ::$F, domain::CartesianSpace, codomain::CartesianSpace)
            u = map((dom, codom) -> _findposition_nzind_domain(ℱ, dom, codom), spaces(domain), spaces(codomain))
            len = sum(length, u)
            v = Vector{Int}(undef, len)
            δ = δδ = 0
            @inbounds for (i, uᵢ) in enumerate(u)
                δ_ = δ
                δ += length(uᵢ)
                view(v, 1+δ_:δ) .= δδ .+ uᵢ
                δδ += dimension(domain[i])
            end
            return v
        end

        function _findposition_nzind_codomain(ℱ::$F, domain::CartesianSpace, codomain::CartesianSpace)
            u = map((dom, codom) -> _findposition_nzind_codomain(ℱ, dom, codom), spaces(domain), spaces(codomain))
            len = sum(length, u)
            v = Vector{Int}(undef, len)
            δ = δδ = 0
            @inbounds for (i, uᵢ) in enumerate(u)
                δ_ = δ
                δ += length(uᵢ)
                view(v, 1+δ_:δ) .= δδ .+ uᵢ
                δδ += dimension(codomain[i])
            end
            return v
        end

        function _project!(C::LinearOperator{<:CartesianSpace,<:CartesianSpace}, ℱ::$F)
            @inbounds for i ∈ 1:nb_cartesian_product(domain(C))
                _project!(component(C, i, i), ℱ)
            end
            return C
        end
    end
end


text:
module {PKGNAME}

end # module


text:
#
# Calculates the SVD of a data matrix containing missing entries.
#
# Uses the iterative algorithm of Hastie et al. 1999
#

# Impute a missing entries using current approximation.
function impute!(X::Matrix, missing_entries::Vector,
                 U::Matrix, D::Vector, V::AbstractMatrix,
                 k::Integer)
    approximation = U[:, 1:k] * diagm(D[1:k]) * V[1:k, :]
    for indices in missing_entries
        X[indices[1], indices[2]] = approximation[indices[1], indices[2]]
    end
end

# Should be done with a proper N-dimensional Int array.
function findmissing(dm::DataMatrix)
    indices = Any[]
    n, p = size(dm)
    for i = 1:n
        for j = 1:p
            if ismissing(dm[i, j])
                push!(indices, [i, j])
            end
        end
    end
    indices
end

function global_mean(dm::DataMatrix)
    mu = 0.0
    n = 0
    n, p = size(dm)
    for i = 1:n
        for j = 1:p
            if !ismissing(dm[i, j])
                mu += dm[i, j]
                n += 1
            end
        end
    end
    return mu / n
end

function na_safe_rowmeans(dm::DataMatrix)
    n, p = size(dm)
    mus = DataArray(Float64, n)
    for i = 1:n
        mu = 0.0
        n = 0
        for j = 1:p
            if !ismissing(dm[i, j])
                mu += dm[i, j]
                n += 1
            end
        end
        if n != 0
            mus[i] = mu / n
        end
    end
    return mus
end

# TODO: Default to failure in the face of missings
function Base.svd(D::DataMatrix, k::Int; tracing = false, tolerance = 10e-4)

    # Make a copy of the data that we can alter in place
    dm = copy(D)

    # Cache the dimensions of the matrix
    n, p = size(dm)

    # Estimate missingness and print a message.
    missing_entries = findmissing(dm)
    missingness = length(missing_entries) / (n * p)
    if tracing
        @printf "Matrix is missing %.2f%% of entries\n" missingness * 100
    end

    # Initial imputation uses global mean and row means
    global_mu = global_mean(dm)
    mu_i = na_safe_rowmeans(dm)
    for i = 1:n
        for j = 1:p
            if ismissing(dm[i, j])
                if ismissing(mu_i[i])
                    dm[i, j] = global_mu
                else
                    dm[i, j] = mu_i[i]
                end
            end
        end
    end

    # Convert dm to a Float array now that we've removed all missings
    dm = float(dm)

    # Count iterations of proper imputation method
    itr = 0

    # Keep track of approximate matrices
    previous_dm = copy(dm)
    current_dm = copy(dm)

    # Keep track of Frobenius norm of changes in imputed matrix
    change = Inf

    # Iterate until imputation stops changing up to chosen tolerance
    while change > tolerance
        if tracing
            @printf "Iteration %d\nChange %f\n" itr change
        end

        # Impute missing entries using current SVD
        previous_dm = copy(current_dm)
        U, D, V = svd(current_dm)
        impute!(current_dm, missing_entries, U, D, V', k)

        # Compute the change in the matrix across iterations
        change = norm(previous_dm - current_dm) / norm(dm)

        # Increment the iteration counter
        itr = itr + 1
    end

    # Tell the user how many iterations were required to impute matrix
    if tracing
        @printf "Tolerance achieved after %d iterations" itr
    end

    # Return the rank-k SVD entries
    U, D, V = svd(current_dm)

    # Only return the SVD entries, not the imputation
    return (U[:, 1:k], D[1:k], V[:, 1:k])
end
Base.svd(dm::DataMatrix) = svd(dm, minimum(size(dm)))

function Base.eig(dm::DataMatrix)
    U, D, V = svd(dm)
    return eig(U * diagm(D) * V')
end


text:
using ValueHistories

function ValueHistories.History(iters::Vector, vals::Vector)
    hist = History(eltype(vals), eltype(iters))
    hist.iterations = copy(iters)
    hist.values = copy(vals)
    hist.lastiter = last(iters)
    return hist
end

function Base.real(hist::History)
    it, val = get(hist)
    return History(it, real(val))
end

function Base.imag(hist::History)
    it, val = get(hist)
    return History(it, imag(val))
end


text:
using CUDA
using Adapt
using DataStructures

struct DummyApp{A,B}
    number_time_steps::Int32
    number_inputs::Int64
    number_outputs::Int64
    number_checkboxes::Int64
    number_gui_elements::Int64
    gui_elements_rectangles::A
    checkboxes_checked::B
end

function DummyApp(configuration::OrderedDict, number_individuals::Int)

    number_checkboxes = 8
    number_gui_elements = number_checkboxes + 1

    # Initialize Cuda Array for positions of GUI elements
    gui_elements_rectangles = CUDA.fill(0.0, (number_gui_elements, 4))

    checkboxes_size = 20
    checkboxes_grid_size = 50
    checkboxes_border = 30

    # Place all 8 checkboxes in 4 colums and 2 rows
    n = 1
    for i = 1:4
        for j = 1:2

            x = checkboxes_grid_size * i + checkboxes_border
            y = checkboxes_grid_size * j + checkboxes_border

            gui_elements_rectangles[n, :] = [x, y, checkboxes_size, checkboxes_size]
            n += 1
        end
    end

    # Place button (below checkboxes)
    gui_elements_rectangles[n, :] = [200, 200, 150, 50]

    DummyApp(
        convert(Int32, configuration["number_time_steps"]),
        number_checkboxes,
        2,
        number_checkboxes,
        number_gui_elements,
        gui_elements_rectangles,
        CUDA.fill(false, (number_checkboxes, number_individuals)),
    )
end

Adapt.@adapt_structure DummyApp

function kernel_eval_fitness(individuals, rewards, environment_seeds, number_rounds, brains, environments)

    threadID = threadIdx().x
    blockID = blockIdx().x

    offset_shared_memory = 0

    reward = @cuDynamicSharedMem(Float32, environments.number_gui_elements, offset_shared_memory)
    offset_shared_memory += sizeof(reward)

    sync_threads()

    observation = @cuDynamicSharedMem(Float32, environments.number_inputs, offset_shared_memory)
    offset_shared_memory += sizeof(observation)

    sync_threads()

    action = @cuDynamicSharedMem(Float32, environments.number_outputs, offset_shared_memory)
    offset_shared_memory += sizeof(action)

    sync_threads()

    click_position = @cuDynamicSharedMem(Int32, environments.number_outputs, offset_shared_memory)
    offset_shared_memory += sizeof(click_position)

    sync_threads()

    gui_elements_states = @cuDynamicSharedMem(Int32, environments.number_gui_elements, offset_shared_memory)
    offset_shared_memory += sizeof(gui_elements_states)

    sync_threads()

    # Initialize brains
    initialize(brains, individuals)

    sync_threads()

    # Initialize states of gui elements
    if threadID <= environments.number_gui_elements
        gui_elements_states[threadID] = 0
    end

    sync_threads()

    # Initialize rewards
    if threadID <= environments.number_gui_elements
        reward[threadID] = 0.0
    end

    sync_threads()

    # Iterate over given number of time steps
    for time_step = 1:environments.number_time_steps

        # Set observations
        if threadID <= environments.number_checkboxes
            observation[threadID] = gui_elements_states[threadID]
        end

        sync_threads()

        # Brain step
        step(brains, observation, action, offset_shared_memory)

        # Scale actions to click positions
        if threadID <= 2
            click_position[threadID] = trunc(0.5 * (action[threadID] + 1.0) * 400.0)
            #click_position[threadID] = rand(1:400)
        end

        sync_threads()

        #if blockID == 1 && threadID == 1
        #    @cuprintln("tx=", threadID,"   action_x=", action[1], "  action_y=", action[2])
        #end

        # Process mouse click
        process_click(environments, click_position, gui_elements_states, reward)

        sync_threads()

    end

    # Transfer rewards to output vector
    if threadID == 1
        rewards[blockID] = 0.0

        for i in 1:environments.number_gui_elements
            rewards[blockID] += reward[i]
        end
    end

    return
end

function get_memory_requirements(environments::DummyApp)
    return sizeof(Float32) * environments.number_gui_elements +                             # Reward
           sizeof(Float32) * (environments.number_inputs + environments.number_outputs) +   # Observation + Action
           sizeof(Int32) * environments.number_gui_elements                                 # States of gui elements

end

function get_number_inputs(environments::DummyApp)
    return environments.number_inputs
end

function get_number_outputs(environments::DummyApp)
    return environments.number_outputs
end

function process_click(environments, point, gui_elements_states, reward)

    threadID = threadIdx().x

    if threadID <= environments.number_gui_elements
        rect_x = environments.gui_elements_rectangles[threadID, 1]
        rect_y = environments.gui_elements_rectangles[threadID, 2]
        width = environments.gui_elements_rectangles[threadID, 3]
        height = environments.gui_elements_rectangles[threadID, 4]

        if is_point_in_rect(point[1], point[2], rect_x, rect_y, width, height)

            # Toggle checkboxes
            if threadID <= environments.number_checkboxes
                # Bitwise xor with 1 does toggling: https://stackoverflow.com/questions/11604409/how-to-toggle-a-boolean
                #gui_elements_states[threadID] ⊻= 1

                if gui_elements_states[threadID] == 0
                    reward[threadID] += 1.0
                    gui_elements_states[threadID] = 1
                end
            end
        end
    end

    sync_threads()

    #if blockID == 1 && threadID == 1
    #    @cuprintln("tx=", threadID,"   x=", point[1], "  y=", point[2])
    #end

    return
end

function is_point_in_rect(point_x, point_y, rect_x, rect_y, width, height)

    return rect_x <= point_x <= (rect_x + width) && rect_y <= point_y <= (rect_y + height)
end


text:
mutable struct Posterior
  MeanHistogram::DataFrame
  Parameters::DataFrame
  Probability::Float64
end

mutable struct Results
  ABCsetup
  ABCresults
  VAF::Array{Float64, 1}
  Posterior::Array{Posterior, 1}
  ModelProb::DataFrame
  SampleName::String
end


text:
using Distributions
import Base.size, Base.length, Base.rand, Base.mean, Base.convert, Base.eltype, Distributions.logpdf, Distributions.entropy, Distributions.params

################### Hidden Markov Model distribution #####################
immutable HMM{S <: Real} <: DiscreteMatrixDistribution
    # convention: Matrices are state x time
    ψ::Matrix{S}  # matrix of conditional probabilites of symbol emission
    π0::Vector{S}  # vector of initial state probabilities
    A::Matrix{S}  # transition matrix (columns sum to 1)
    ξ::Matrix{S}  # mean of distribution
    Ξ::Array{S, 3}  # two-slice marginal of distribution (state x state x time)
    logZ::S  # log normalization constant for distribution

    function HMM(ψ, π0, A)
        size(ψ, 1) == length(π0) == size(A, 1) == size(A, 2) || error("A and π0 have incompatible shapes")

        M, T = size(ψ)

        sum(π0) ≈ 1 || warn("Entries of π0 do not sum to 1.")

        sum(A, 1) ≈ ones(1, M) || warn("Columns of A do not sum to 1.")

        ξ, logZ, Ξ = forwardbackward(π0, A, ψ)

        new(ψ, π0, A, ξ, Ξ, logZ)
    end
end

HMM{S <: Real}(ψ::Matrix{S}, π0::Vector{S}, A::Matrix{S}) = HMM{S}(ψ, π0, A)

#### Conversions
convert{S <: Real, V <: Real}(::Type{HMM{S}}, ψ::Matrix{V}, π0::Vector{V}, A::Matrix{V}) = HMM(convert(Matrix{S}, ψ), convert(Vector{S}, π0), convert(Matrix{S}, A))
convert{S <: Real, V <: Real}(::Type{HMM{S}}, d::HMM{V}) = HMM(convert(Matrix{S}, d.ψ), convert(Vector{S}, d.π0), convert(Matrix{S}, d.A))

nstates(d::HMM) = length(d.π0)
size(d::HMM) = size(d.ψ)
length(d::HMM) = prod(size(d))

mean(d::HMM) = d.ξ
cov(d::HMM) = d.Ξ
params(d::HMM) = (d.ψ, d.π0, d.A)

function rand(d::HMM)
    M, T = size(d)
    z = zeros(M, T)
    init_state = rand(Categorical(d.ξ[:, 1]))
    z[init_state, 1] = 1
    for t in 2:T
        pvec = d.Ξ[:, findfirst(z[:, t - 1]), t - 1]
        pvec /= sum(pvec)  # normalize, since Ξ is joint, not conditional
        newstate = rand(Categorical(pvec))
        z[newstate, t] = 1
    end
    z
end

function logpdf{S <: Real}(d::HMM{S}, x::Matrix{S})
    size(x) == size(d) || error("Input matrix x has wrong size.")

    T = size(x, 2)

    emission_piece = sum(x .* log(d.ψ))
    initial_piece = sum(x[:, 1] .* log(d.π0))
    transition_piece = 0.
    for t in 1:(T - 1)
        transition_piece += (x[:, t + 1]' * log(d.A) * x[:, t])[1]
    end

    emission_piece + initial_piece + transition_piece - d.logZ
end

function entropy(d::HMM)
    _, T = size(d)

    emission_piece = sum(d.ξ .* log(d.ψ))
    initial_piece = sum(d.ξ[:, 1] .* log(d.π0))
    transition_piece = 0.
    for t in 1:(T - 1)
        transition_piece += sum(d.Ξ[:, :, t] .* log(d.A))
    end

    d.logZ -(emission_piece + initial_piece + transition_piece)
end

function naturals(d::HMM)
    (log(d.ψ), log(d.π0), log(d.A))
end

function naturals_to_params{D <: HMM}(η, ::Type{D})
    map(exp, η)
end

################### Markov Chain distribution #####################
immutable MarkovChain{S <: Real} <: DiscreteMatrixDistribution
    # convention: Matrices are state x time
    π0::Vector{S}  # vector of initial state probabilities
    A::Matrix{S}  # transition matrix (columns sum to 1)
    T::Integer

    function MarkovChain(π0, A, T)
        isa(T, Int) || error("Supplied chain length is not an integer.")
        length(π0) == size(A, 1) == size(A, 2) || error("A and π0 have incompatible shapes")

        sum(π0) ≈ 1 || warn("Entries of π0 do not sum to 1.")

        sum(A, 1) ≈ ones(1, size(A, 1)) || warn("Columns of A do not sum to 1.")

        new(π0, A, T)
    end
end

MarkovChain{S <: Real}(π0::Vector{S}, A::Matrix{S}, T) = MarkovChain{S}(π0, A, T)
function MarkovChain{S <: Real, V <: Real}(π0::Vector{S}, A::Matrix{V}, T)
    W = promote_type(eltype(π0), eltype(A))
    MarkovChain(convert(Vector{W}, π0), convert(Matrix{W}, A), T)
end

#### Conversions
convert{S <: Real, V <: Real}(::Type{MarkovChain{S}}, π0::Vector{V}, A::Matrix{V}, T) = MarkovChain(convert(Vector{S}, π0), convert(Matrix{S}, A), T)
convert{S <: Real, V <: Real}(::Type{MarkovChain{S}}, d::MarkovChain{V}) = MarkovChain(convert(Vector{S}, d.π0), convert(Matrix{S}, d.A), d.T)

nstates(d::MarkovChain) = length(d.π0)
size(d::MarkovChain) = (nstates(d), d.T)
length(d::MarkovChain) = prod(size(d))

function mean(d::MarkovChain)
    M, T = size(d)
    m = zeros(M, T)
    m[:, 1] = d.π0
    for i in 2:T
        m[:, i] = d.A * m[:, i - 1]
    end
    m
end

params(d::MarkovChain) = (d.π0, d.A)

function rand(d::MarkovChain)
    M, T = size(d)
    z = zeros(M, T)
    init_state = rand(Categorical(d.π0))
    z[init_state, 1] = 1
    for t in 2:T
        pvec = d.A * z[:, t - 1]
        pvec /= sum(pvec)  # normalize, since Ξ is joint, not conditional
        newstate = rand(Categorical(pvec))
        z[newstate, t] = 1
    end
    z
end

function logpdf{S <: Real}(d::MarkovChain{S}, x::Matrix{S})
    size(x) == size(d) || error("Input matrix x has wrong size.")

    T = size(x, 2)

    initial_piece = sum(x[:, 1] .* log(d.π0))
    transition_piece = 0.
    for t in 1:(T - 1)
        transition_piece += (x[:, t + 1]' * log(d.A) * x[:, t])[1]
    end

    initial_piece + transition_piece
end

function entropy(d::MarkovChain)
    _, T = size(d)

    initial_piece = sum(d.π0 .* log(d.π0))
    transition_piece = 0.
    for t in 1:(T - 1)
        transition_piece += sum(d.A .* log(d.A))
    end

    -(initial_piece + transition_piece)
end

function naturals(d::MarkovChain)
    (log(d.π0), log(d.A))
end

function naturals_to_params{D <: MarkovChain}(η, ::Type{D})
    map(exp, η)
end
################### Markov Matrix distribution #####################
immutable MarkovMatrix{T <: Real} <: ContinuousMatrixDistribution
    cols::Vector{Dirichlet{T}}  # each column is a Dirichlet distribution

    function MarkovMatrix(cols)
        length(cols) == length(cols[1]) || error("Input is not a square matrix.")

        new(cols)
    end
end

MarkovMatrix{T <: Real}(cols::Vector{Dirichlet{T}}) = MarkovMatrix{T}(cols)

function MarkovMatrix{T <: Real}(A::Matrix{T})
    r, c = size(A)
    r == c || error("Input matrix must be square.")

    MarkovMatrix([Dirichlet(A[:, i]) for i in 1:c])
end

function MarkovMatrix{T <: Real}(cols::Vector{Vector{T}})
    p = length(cols)
    all(x -> length(x) == p, cols) || error("Probability vector lengths do not match.")

    MarkovMatrix([Dirichlet(c) for c in cols])
end

#### Conversions
convert{T <: Real, S <: Real}(::Type{MarkovMatrix{T}}, cols::Vector{Dirichlet{S}}) = MarkovMatrix([convert(Dirichlet{T}, c) for c in cols])
convert{T <: Real, S <: Real}(::Type{MarkovMatrix{T}}, A::Matrix{S}) = MarkovMatrix(convert(Matrix{T}, A))
convert{T <: Real, S <: Real}(::Type{MarkovMatrix{T}}, cols::Vector{Vector{S}}) = MarkovMatrix([convert(Vector{T}, c) for c in cols])
convert{T <: Real, S <: Real}(::Type{MarkovMatrix{T}}, d::MarkovMatrix{S}) = MarkovMatrix([convert(Dirichlet{T}, c) for c in d.cols])

eltype{T <: Real}(d::MarkovMatrix{T}) = T
nstates(d::MarkovMatrix) = length(d.cols)
size(d::MarkovMatrix) = (length(d.cols), length(d.cols))
length(d::MarkovMatrix) = length(d.cols)^2

function mean(d::MarkovMatrix)
    p = nstates(d)
    m = Array{Float64}(p, p)
    for i in 1:p
        m[:, i] = mean(d.cols[i])
    end
    m
end

function Elog(d::MarkovMatrix)
    p = nstates(d)
    m = Array{Float64}(p, p)
    for i in 1:p
        m[:, i] = Elog(d.cols[i])
    end
    m
end

function rand(d::MarkovMatrix)
    p = nstates(d)
    x = Array{Float64}(p, p)
    for i in 1:p
        x[:, i] = rand(d.cols[i])
    end
    x
end

function logpdf{T <: Real}(d::MarkovMatrix, x::Matrix{T})
    size(x) == size(d) || error("Input matrix x has wrong size.")
    s = 0.
    for i in 1:nstates(d)
        s += logpdf(d.cols[i], x[:, i])
    end
    s
end

function entropy(d::MarkovMatrix)
    s = 0.
    for c in d.cols
        s += entropy(c)
    end
    s
end

function naturals(d::MarkovMatrix)
    p = nstates(d)
    nats = Array{eltype(d)}(p, p)
    for i in 1:p
        nats[:, i] = naturals(d.cols[i])[1]  # naturals returns a tuple
    end
    (nats,)
end

function params(d::MarkovMatrix)
    p = nstates(d)
    nats = Array{eltype(d)}(p, p)
    for i in 1:p
        nats[:, i] = params(d.cols[i])[1]  # params returns a tuple
    end
    (nats,)
end

function naturals_to_params{D <: MarkovMatrix}(η, ::Type{D})
    nats = η[1]
    p, _ = size(nats)
    pars = similar(nats)
    for i in 1:p
        pars[:, i] = naturals_to_params((nats[:, i],), Dirichlet)[1]  # n2par returns a tuple
    end
    (pars,)
end

"""
Implement the forward-backward inference algorithm.
A is a matrix of transition probabilities that acts to the right:
new_state = A * old_state, so that columns of A sum to 1
ψ is the vector of evidence: p(y_t|z_t); it does not need to be
normalized, but the lack of normalization will be reflected in logZ
such that the end result using the given ψ will be properly normalized
when using the returned value of logZ.
"""
function forwardbackward(π0, A, ψ)
    if any(π0 .> 1) throw(ArgumentError()) end
    if any(A .> 1) throw(ArgumentError()) end

    # get shape parameters
    M, T = size(ψ)

    # allocate empty variables
    α = Array{Float64}(M, T)
    β = Array{Float64}(M, T)
    γ = Array{Float64}(M, T)
    Ξ = Array{Float64}(M, M, T - 1)

    # initialize
    a = ψ[:, 1] .* π0
    α[:, 1] = a / sum(a)
    logZ = log(sum(a))
    β[:, T] = 1

    # forward pass
    for t in 2:T
        asum = 0.
        for i in 1:M
            a[i] = 0.
            for j in 1:M
                a[i] += ψ[i, t] * A[i, j] * α[j, t - 1]
            end
            asum += a[i]
        end

        for i in 1:M
            α[i, t] = a[i] / asum
        end

        logZ += log(asum)
    end

    # backward pass
    for t in T:-1:2
        asum = 0.
        for i in 1:M
            a[i] = 0.
            for j in 1:M
                a[i] += β[j, t] * ψ[j, t] * A[j, i]
            end
            asum += a[i]
        end

        for i in 1:M
            β[i, t - 1] = a[i] / asum
        end
    end

    # calculate posterior
    for t in 1:T
        γsum = 0.
        for i in 1:M
            γ[i, t] = α[i, t] * β[i, t]
            γsum += γ[i, t]
        end

        for i in 1:M
            γ[i, t] /= γsum
        end
    end

    # calculate two-slice marginals
    for t in 1:(T - 1)
        xsum = 0.
        for i in 1:M
            for j in 1:M
                Ξ[i, j, t] = β[i, t + 1] * ψ[i, t + 1]
                Ξ[i, j, t] *= α[j, t] * A[i, j]
                xsum += Ξ[i, j, t]
            end
        end

        #normalize
        for i in 1:M
            for j in 1:M
                Ξ[i, j, t] /= xsum
            end
        end
    end

    γ, logZ, Ξ

end


text:
function bfs(G, source)
    visited = Set(source);
    q = Deque{typeof(source)}()
    push!(q, source);
    while ! isempty(q)
        v = popfirst!(q);
        for n in neighbors(G, v)
            if !(n in visited)
                push!(q, n);
                push!(visited, n);
            end
        end
    end
    return visited
end

function connected_components(G)
    seen = []
    components = []
    for v in vertices(G)
        if !(v in seen)
            c = bfs(G, v)
            [push!(seen, elem) for elem in c]
            push!(components, c)
        end
    end
    return components
end


text:
import REPL

using REPL.TerminalMenus
using CSV
using DataFrames, Query

export rebalance

function list_csv(path::String)
  files = [file for file in readdir(path) if isfile(abspath(joinpath(path, file)))]
  files = [m.match for m in match.(r".*\.csv", files) if m != nothing]
  return (path * "/").*(files)
end

function choose(prompt, options)
  menu = RadioMenu(options)
  choice = request(prompt, menu)
  return options[choice]
end

function request_csv(path::String, prompt::String)
  options = list_csv(path)
  return choose(prompt, options)
end

function open_fidelity_csv(path::String=homedir() * "/Downloads")
  filename = request_csv(path, "Open Fidelity CSV:")
  df = DataFrame(CSV.File(filename; footerskip=6, normalizenames=true, stringtype=String))
  return df[:,2:8]
end

function open_strategy_csv(path::String=homedir() * "/Documents")
  filename = request_csv(path, "Open Targets CSV:")
  return DataFrame(CSV.File(filename; normalizenames=true, stringtype=String))
end

function by_account(accounts::DataFrame, name::String)
  @from a in accounts begin
    @where a.Account_Name == name
    @select a
    @collect DataFrame
  end
end

function get_account(accounts::DataFrame)
  account_names = @from a in accounts begin
    @group a by a.Account_Name into n
    @select {Account_Name=key(n)}
    @collect DataFrame
  end
  account_name = choose("Choose account:", account_names[:, 1])
  return by_account(accounts, account_name)
end

function target_symbol(taxable::Bool, target)
  if taxable
    return target.Symbol_Tax
  else
    return target.Symbol_Adv
  end
end

function make_target(strategy::DataFrame)
  taxable = ("Taxable" == choose("Account Type:", ["Taxable", "Tax-Advantaged"]))
  levels = ["Conservative", "Moderate", "Balanced", "Growth/Income", "Growth", "Aggressive"]
  target = choose("Target:", levels)
  sort!(strategy, :Target, rev=true)
  # starts at 30%
  equity_percent = (findfirst(t -> t == target, levels) + 2) * 10
  equities = @from t in strategy begin
    @where !Bool(t.Income)
    @select {Symbol=target_symbol(taxable, t), t.Target}
    @collect DataFrame
  end
  equities.Target = allocate(equity_percent, equities.Target)
  income = @from t in strategy begin
    @where Bool(t.Income)
    @select {Symbol=target_symbol(taxable, t), t.Target}
    @collect DataFrame
  end
  income_percent = 100 - equity_percent
  income.Target = allocate(income_percent, income.Target)
  return sort(vcat(equities, income), :Target, rev=true)
end

function unshift(trades::DataFrame)
  return trades[1, :], trades[2:nrow(trades), :]
end

function is_core(symbol)
  return endswith(symbol, "**") || symbol == "Pending Activity"
end

function Base.get(a::String, b)
  return a
end

function sum_dtoi(c::Vector{String})
  return sum(dtoi.(c))
end

function sum_dtoi(c::Vector{Union{Missing, String}})
  result = 0
  for d in c
    result += dtoi(get(d, "\$0.00"))
  end
  return result
end

function get_core(account::DataFrame)
  core = @from asset in account begin
    @select {asset.Symbol, Current_Value=get(asset.Current_Value, asset.Last_Price_Change)}
    @collect DataFrame
  end
  return filter(asset -> is_core(get(asset.Symbol, "?")), core)
end

function get_non_core(account::DataFrame)
  nonCore = @from asset in account begin
    @select {asset.Symbol, asset.Description, asset.Quantity, asset.Current_Value}
    @collect DataFrame
  end
  return filter(asset -> !is_core(get(asset.Symbol, "?")), nonCore)
end

function get_holdings(target::DataFrame, nonCore::DataFrame)
  return @from asset in leftjoin(target, nonCore, on = :Symbol) begin
    @orderby descending(asset.Target), descending(dtoi(get(asset.Current_Value, "\$0.00")))
    @select {asset.Symbol, asset.Description, Quantity=get(asset.Quantity, 0.0), Current_Value=get(asset.Current_Value, "\$0.00"), asset.Target}
    @collect DataFrame
  end
end

function get_exiting(target::DataFrame, nonCore::DataFrame)
  exits = @from asset in nonCore begin
    @orderby descending(dtoi(get(asset.Current_Value, "\$0.00")))
    @select {asset.Symbol, asset.Description, asset.Quantity, asset.Current_Value, Target=0, Trade_Type="SELL", Trade_Value=itod(-dtoi(get(asset.Current_Value, "\$0.00"))), Drift=NaN, DriftPct=NaN}
    @collect DataFrame
  end
  return filter(e -> !(e.Symbol in target.Symbol), exits)
end

function generate_trades(account::DataFrame, target::DataFrame, deposit::Int=0)
  # core cash and pending
  core = get_core(account)
  println("\n---- HOLDINGS ----")
  if nrow(core) > 0
    printframe(core)
  end
  # assets other than core position
  nonCore = get_non_core(account)
  # account totals
  core_value = sum_dtoi(core.Current_Value)
  non_core_value = sum_dtoi(nonCore.Current_Value)
  cash = core_value + deposit
  accountTotal = core_value + non_core_value
  newTotal = non_core_value + cash
  # assets not in target (sell them)
  exiting = get_exiting(target, nonCore)
  # target assets (trade them)
  holdings = get_holdings(target, nonCore)
  # calculate drift from target
  holdings[!, :Drift] = map(h -> round(((dtoi(h.Current_Value)*100/accountTotal) - h.Target)*10, RoundNearestTiesAway)/10, eachrow(holdings))
  holdings[!, :DriftPct] = map(h -> round(h.Drift*100/h.Target*10)/10, eachrow(holdings))
  printframe(holdings)
  println("CASH | ", itod(core_value))
  # trade value in dollars
  tradeAmount = allocate(newTotal, holdings.Target[:,1]) - map(h -> dtoi(h), holdings.Current_Value[:,1])
  tradeType = map(t -> t < 0 ? "SELL" : "BUY", tradeAmount)
  holdings[!, :Trade_Type] = tradeType
  holdings[!, :Trade_Value] = itod.(tradeAmount)
  # combine all tradable holdings
  trades = vcat(holdings, exiting)
  # create exchange trades, sort
  first, second = optimize_trades(trades, cash)
  println("\n---- TRADES ----")
  settle = print_trades(first, core_value)
  if nrow(second) > 0
    settle = print_trades(second, settle)
  end
end

function optimize_trades(trades::DataFrame, cash::Int)
  # split trades into buys and sells
  sells = @from trade in trades begin
    @where trade.Trade_Type == "SELL"
    @orderby dtoi(trade.Trade_Value)
    @select {trade.Symbol, trade.Description, trade.Quantity, trade.Current_Value, trade.Trade_Type, trade.Trade_Value, Trade_For="CASH", Trade_Shares=NaN}
    @collect DataFrame
  end
  buys = @from trade in trades begin
    @where trade.Trade_Type == "BUY"
    @orderby descending(dtoi(trade.Trade_Value))
    @select {trade.Symbol, trade.Description, trade.Quantity, trade.Current_Value, trade.Trade_Type, trade.Trade_Value, Trade_For="", Trade_Shares=NaN}
    @collect DataFrame
  end
  # find mutual fund exchanges
  trades = similar(sells, 0)
  while true
    if nrow(sells) == 0
      # out of sells, append remaining buys
      append!(trades, buys)
      break
    end
    sell, sells = unshift(sells)
    if !endswith(get(sell.Symbol, "?"), "X")
      # not exchangable sell, push it
      push!(trades, sell)
      continue
    end
    if nrow(buys) == 0
      # out of buys, append remaining sells
      push!(trades, sell)
      append!(trades, sells)
      break
    end
    buy, buys = unshift(buys)
    if !endswith(get(buy.Symbol, "?"), "X")
      # not exchangable buy, push it
      push!(trades, buy)
      push!(sells, sell)
      sells = @from trade in sells begin
        @orderby dtoi(trade.Trade_Value)
        @select trade
        @collect DataFrame
      end
      continue
    end
    # exchange mutual funds
    if abs(dtoi(sell.Trade_Value)) >= dtoi(buy.Trade_Value)
      # sell larger than buy - reduce sell and turn buy into exchange
      sell.Trade_Value = itod(dtoi(sell.Trade_Value) + dtoi(buy.Trade_Value))
      buy.Trade_For = buy.Symbol
      buy.Symbol = sell.Symbol
      buy.Description = sell.Description
      buy.Quantity = sell.Quantity
      buy.Current_Value = sell.Current_Value
      buy.Trade_Type = "EXCHANGE"
      buy.Trade_Value = itod(-dtoi(buy.Trade_Value))
      push!(sells, sell)
      push!(trades, buy)
      sells = @from trade in sells begin
        @orderby dtoi(trade.Trade_Value)
        @select trade
        @collect DataFrame
      end
    else
      # buy larger than sell, reduce buy and turn sell into exchange
      buy.Trade_Value = itod(dtoi(sell.Trade_Value) + dtoi(buy.Trade_Value))
      sell.Trade_Type = "EXCHANGE"
      sell.Trade_For = buy.Symbol
      push!(buys, buy)
      push!(trades, sell)
      buys = @from trade in buys begin
        @orderby descending(dtoi(trade.Trade_Value))
        @select trade
        @collect DataFrame
      end
    end
  end
  # filter-out zero-value trades
  trades = @from trade in trades begin
    @where dtoi(trade.Trade_Value) != 0
    @orderby trade.Symbol
    @select trade
    @collect DataFrame
  end
  # calculate number of shares for sells and exchanges
  symbol = ""
  remaining_quantity = 0
  remaining_value = 0
  for trade in eachrow(trades)
    if trade.Trade_Type == "BUY"
      # skip trade-quantity for buys
      trade.Trade_Shares = NaN
      trade.Trade_For = ""
      continue
    end
    if trade.Symbol != symbol
      # reset on next symbol
      symbol = trade.Symbol
      remaining_quantity = qtoi(trade.Quantity)
      remaining_value = dtoi(trade.Current_Value)
    end
    # value for this sell/exchange, and remainder
    trade_value = abs(dtoi(trade.Trade_Value))
    remaining_value -= trade_value
    if (remaining_value == 0)
      trade_quantity = remaining_quantity
    else
      # allocate shares in thousandths
      trade_quantity, remaining_quantity = allocate(remaining_quantity, [trade_value, remaining_value])
    end
    trade.Trade_Shares = itoq(trade_quantity)
  end
  # split trades into sells-first, and buys-second
  first = @from trade in trades begin
    @where trade.Trade_Type != "BUY"
    @orderby endswith(get(trade.Symbol, "?"), "XX"), descending(trade.Trade_Type), trade.Symbol, descending(abs(dtoi(trade.Trade_Value)))
    @select trade
    @collect DataFrame
  end
  second = @from trade in trades begin
    @where trade.Trade_Type == "BUY"
    # least-expensive buys first
    @orderby dtoi(trade.Trade_Value)
    @select trade
    @collect DataFrame
  end
  # non-mutual fund sales - for funding trades
  non_mutual = @from trade in first begin
    @where !endswith(get(trade.Symbol, "?"), "X")
    @select trade
    @collect DataFrame
  end
  # amount of funding for first trades
  available = max(0, cash) - sum_dtoi(non_mutual.Trade_Value)
  # move buys from second to first, if we can afford them
  while true
    if nrow(second) == 0
      break
    end
    # least-expensive buy...
    buy, second = unshift(second)
    if dtoi(buy.Trade_Value) <= available
      available -= dtoi(buy.Trade_Value)
      push!(first, buy)
    else
      push!(second, buy)
      break
    end
  end
  # final sort
  first = @from t in first begin
    @orderby endswith(get(t.Trade_For, "?"), "XX"), endswith(get(t.Symbol, "?"), "XX"), endswith(get(t.Symbol, "?"), "X"), descending(t.Trade_Type), t.Symbol, descending(abs(dtoi(t.Trade_Value)))
    @select {t.Symbol, t.Trade_Type, t.Trade_Value, t.Trade_Shares, t.Trade_For}
    @collect DataFrame
  end
  second = @from t in second begin
    @orderby endswith(get(t.Symbol, "?"), "XX"), endswith(get(t.Symbol, "?"), "X"), descending(t.Trade_Type), t.Symbol, descending(abs(dtoi(t.Trade_Value)))
    @select {t.Symbol, t.Trade_Type, t.Trade_Value, t.Trade_Shares, t.Trade_For}
    @collect DataFrame
  end
  return first, second
end

function printframe(df::DataFrame)
  println(replace(sprint(show, df, context=:compact=>false), r".*DataFrame" => ""))
end

function print_trades(trades::DataFrame, core::Int)
  printframe(trades)
  nonex = @from trade in trades begin
    @where trade.Trade_Type != "EXCHANGE"
    @select trade
    @collect DataFrame
  end
  settle = max(0, core) - sum_dtoi(nonex.Trade_Value)
  println("CASH | ", itod(settle))
  return settle
end

function rebalance(deposit::Int=0)
  accounts = open_fidelity_csv()
  strategy = open_strategy_csv()
  println()
  account = get_account(accounts)
  target = make_target(strategy)
  generate_trades(account, target, deposit)
end


text:
# Helper functions
# The implementation is based on the `particle` Python package from
# the SciKit-HEP group: https://github.com/scikit-hep/particle with the Git
# commit 6fe8211181d6a0b21933b2d1042eeb4e87d77b3c from 2020-11-27

"""
    isvalid(p::PDGID)

Returns true if the PDG ID is valid.
"""
function Base.isvalid(p::PDGID)
    fundamentalid(p) != 0 && return true
    ismeson(p) && return true
    isbaryon(p) && return true
    isgaugebosonorhiggs(p) && return true
    ispentaquark(p) && return true
    isSUSY(p) && return true
    isRhadron(p) && return true
    isdyon(p) && return true
    isdiquark(p) && return true
    ispentaquark(p) && return true
    isgeneratorspecific(p) && return true
    istechnicolor(p) && return true
    iscompositequarkorlepton(p) && return true
    !isstandard(p) && return (isQball(p) || isnucleus(p))
    false
end

"""
    isstandard(p::Union{Particle, PDGID, Integer})

Returns true if the PDG ID of the particle follows the standard numbering scheme.
"""
function isstandard(p)
    p = pdgid(p)
    p.N8 == 0 && p.N9 == 0 && p.N10 == 0
end

"""
    isfundamental(p::Union{Particle, PDGID, Integer})

"""
function isfundamental(p)
    p = pdgid(p)
    !isstandard(p) && return false
    p.Nq2 == 0 && p.Nq1 == 0 && return true
    abs(p.value) <= 100 && return true
    false
end

"""
    fundamentalid(p::Union{Particle, PDGID, Integer})

"""
function fundamentalid(p)
    p = pdgid(p)
    !isstandard(p) && return 0
    abspdgid = abs(p.value)
    p.Nq2 == 0 && p.Nq1 == 0 && return abspdgid % 10000
    abs(p.value) <= 100 && return abspgdid
    0
end

"""
    isquark(p::Union{Particle, PDGID, Integer})

"""
isquark(p) = 1 <= abs(pdgid(p).value) <= 8

"""
    islepton(p::Union{Particle, PDGID, Integer})

"""
function islepton(p)
    p = pdgid(p)
    !isstandard(p) && return false
    11 <= fundamentalid(p) <= 18 && return true
    false
end

"""
    ismeson(p::Union{Particle, PDGID, Integer})

"""
function ismeson(p)
    p = pdgid(p)
    !isstandard(p) && return false

    abspdgid = abs(p.value)

    abspdgid <= 100 && return false
    0 < fundamentalid(p) <= 100 && return false
    # Special IDs - K(L)0, ???, K(S)0
    abspdgid ∈ [130, 210, 310] && return true
    # Special IDs - B(L)0, B(sL)0, B(H)0, B(sH)0
    abspdgid ∈ [150, 350, 510, 530] && return true
    # Special particles - reggeon, pomeron, odderon
    p.value ∈ [110, 990, 9990] && return true

    if p.Nj > 0 && p.Nq3 > 0 && p.Nq2 > 0 && p.Nq1 == 0
        # check for illegal antiparticles
        p.Nq3 == p.Nq2 && p.value < 0 && return false
        return true
    end
    false
end

"""
    isbaryon(p::Union{Particle, PDGID, Integer})

"""
function isbaryon(p)
    p = pdgid(p)
    abspdgid = abs(p.value)

    abspdgid <= 100 && return false

    # Special case of proton and neutron:
    # needs to be checked first since isstandard() is false for nuclei
    abspdgid ∈ [1000000010, 1000010010] && return true

    !isstandard(p) && return false
    0 < fundamentalid(p) <= 100 && return false

    # Old codes for diffractive p and n (MC usage)
    abspdgid ∈ [2110, 2210] && return true

    p.Nj > 0 && p.Nq3 > 0 && p.Nq2 > 0 && p.Nq1 > 0 && return true

    (isRhadron(p) || ispentaquark(p)) && return false

    false
end

"""
    ishadron(p::Union{Particle, PDGID, Integer})

"""
function ishadron(p)
    p = pdgid(p)
    abs(p.value) ∈ [1000000010, 1000010010] && return true
    !isstandard(p) && return false
    ismeson(p) && return true
    isbaryon(p) && return true
    ispentaquark(p) && return true
    isRhadron(p) && return true
    false
end

"""
    isRhadron(p::Union{Particle, PDGID, Integer})

An R-hadron is of the form 10abcdj, 100abcj, or 1000abj,
where j = 2J + 1 gives the spin; b, c, and d are quarks or gluons;
and a (the digit following the zero's) is a SUSY particle.

"""
function isRhadron(p)
    p = pdgid(p)
    !isstandard(p) && return false
    (p.N != 1 || p.Nr != 0) && return false
    isSUSY(p) && return false
    # All R-hadrons have at least 3 core digits
    (p.Nq2 == 0 || p.Nq3 == 0 || p.Nj == 0) && return false
    true
end

"""
    isSUSY(p::Union{Particle, PDGID, Integer})

Fundamental SUSY particles have N = 1 or 2.

"""
function isSUSY(p)
    p = pdgid(p)
    !isstandard(p) && return false
    p.N != 1 && p.N != 2 && return false
    p.Nr != 0 && return false
    fundamentalid(p) == 0 && return false
    true
end

"""
    ispentaquark(p::Union{Particle, PDGID, Integer})

Pentaquark IDs are of the form +/- 9 Nr Nl Nq1 Nq2 Nq3 Nj,
where Nj = 2J + 1 gives the spin and Nr Nl Nq1 Nq2 Nq3 denote the quark numbers
in order Nr >= Nl >= Nq1 >= Nq2 and Nq3 gives the antiquark number.

"""
function ispentaquark(p)
    p = pdgid(p)
    !isstandard(p) && return false
    p.N != 9 && return false
    (p.Nr == 9 || p.Nr == 0) && return false
    (p.Nj == 9 || p.Nl == 0) && return false
    p.Nq1 == 0 && return false
    p.Nq2 == 0 && return false
    p.Nq3 == 0 && return false
    p.Nj == 0 && return false
    p.Nq2 > p.Nq1 && return false
    p.Nq1 > p.Nl && return false
    p.Nl > p.Nr && return false
    true
end

"""
    isgaugebosonorhiggs(p::Union{Particle, PDGID, Integer})

"""
function isgaugebosonorhiggs(p)
    p = pdgid(p)
    21 <= abs(p.value) <= 40
end

"""
    issmgaugebosonorhiggs(p::Union{Particle, PDGID, Integer})
"""
function issmgaugebosonorhiggs(p)
    p = pdgid(p)
    abs(p.value) == 24 || 21 <= p.value <= 25
end

"""
    istechnicolor(p::Union{Particle, PDGID, Integer})
"""
function istechnicolor(p)
    p = pdgid(p)
    !isstandard(p) && return false
    p.N == 3
end

"""
    iscompositequarkorlepton(p::Union{Particle, PDGID, Integer})

Excited (composite) quarks and leptons have N = 4 and Nr = 0.
"""
function iscompositequarkorlepton(p)
    p = pdgid(p)
    !isstandard(p) && return false
    fundamentalid(p) == 0 && return false
    !(p.N == 4 && p.Nr == 0) && return false
    true
end

"""
    isdyon(p::Union{Particle, PDGID, Integer})

Magnetic monopoles and Dyons are assumed to have one unit of Dirac monopole
charge and a variable integer number xyz units of electric charge, where xyz
stands for Nq1 Nq2 Nq3.
Codes 411xyz0 are used when the magnetic and electrical charge sign agree and
412xyz0 when they disagree, with the overall sign of the particle set by the
magnetic charge. For now, no spin information is provided.
"""
function isdyon(p)
    p = pdgid(p)
    !isstandard(p) && return false
    p.N != 4 && return false
    p.Nr != 1 && return false
    !(p.Nl ∈ [1, 2]) && return false
    p.Nq3 == 0 && return false
    p.Nj != 0 && return false
    true
end

"""
    isdiquark(p::Union{Particle, PDGID, Integer})
"""
function isdiquark(p)
    p = pdgid(p)
    !isstandard(p) && return false
    abs(p.value) <= 100 && return false
    0 < fundamentalid(p) <= 100 && return false
    p.Nj > 0 && p.Nq3 == 0 && p.Nq2 > 0 && p.Nq1 > 0 && return true
    false
end

"""
    isgeneratorspecific(p::Union{Particle, PDGID, Integer})

Codes 81-100 are reserved for generator-specific pseudoparticles and concepts.
Codes 901-930, 1901-1930, 2901-2930, and 3901-3930 are for additional components
of Standard Model parton distribution functions, where the latter three ranges
are intended to distinguish left/right/longitudinal components.
Codes 998 and 999 are reserved for GEANT tracking purposes.
"""
function isgeneratorspecific(p)
    p = pdgid(p)
    abspdgid = abs(p.value)
    81 <= abspdgid <= 100 && return true
    901 <= abspdgid <= 930 && return true
    1901 <= abspdgid <= 1930 &&  return true
    2901 <= abspdgid <= 2930 && return true
    3901 <= abspdgid <= 3930 && return true
    abspdgid ∈ [998, 999] && return true
    abspdgid ∈ [20022, 480000000] && return true  # Special cases of opticalphoton and geantino
    false
end

"""
    isspecial(p::Union{Particle, PDGID, Integer})

Special particle in the sense of the classification in the PDG MC particle
numbering scheme document, hence the graviton, the DM (S = 0, 1/2, 1) particles,
the reggeons (reggeon, pomeron and odderon), and all generator-specific
pseudo-particles and concepts, see `isgeneratorspecific`.
"""
function isspecial(p)
    p = pdgid(p)
    p.value ∈ [39, 41, 42, 51, 52, 53, 110, 990, 9990] || isgeneratorspecific(p)
end

"""
    isQball(p::Union{Particle, PDGID, Integer})

Does this PDG ID correspond to a Q-ball or any exotic particle with electric
charge beyond the qqq scheme?
Ad-hoc numbering for such particles is +/- 100XXXY0, where XXX.Y is the charge.
"""
function isQball(p)
    p = pdgid(p)
    p.N8 != 1 && return false
    p.N != 0 && return false
    p.Nr != 0 && return false
    (abs(p.value) ÷ 10) % 10000 == 0 && return false
    p.Nj != 0 && return false
    true
end

"""
    hasfundamentalanti(p::Union{Particle, PDGID, Integer})

If this is a fundamental particle, does it have a valid antiparticle?

Notes
-----
Based on the current list of defined particles/concepts
in the PDG Monte Carlo Particle Numbering Scheme document.

"""
function hasfundamentalanti(p)
    p = pdgid(p)
    fid = fundamentalid(p)  # always a positive integer

    # Check generator-specific PDGIDs
    81 <= fid <= 100 && return fid ∈ [82, 84, 85, 86, 87]

    # Check PDGIDs from 1 to 79
    cp_conjugates = [21, 22, 23, 25, 32, 33, 35, 36, 39, 40, 43]
    unassigned = vcat([9, 10, 19, 20, 26], 26:31, 45:79)  # not in conversion.csv

    if (1 <= fid <= 79) && !(fid ∈ cp_conjugates)
        fid ∈ unassigned && return false
        return true
    end
    false
end

"""
    isnucleus(p::Union{Particle, PDGID, Integer})

Ion numbers are +/- 10LZZZAAAI.
AAA is A - total baryon number
ZZZ is Z - total charge
L is the total number of strange quarks.
I is the isomer number, with I=0 corresponding to the ground state.
"""
function isnucleus(p)
    p = pdgid(p)
    # A proton can be a Hydrogen nucleus
    # A neutron can be considered as a nucleus when given the PDG ID 1000000010,
    # hence consistency demands that isnucleus(neutron) is True
    abs(p.value) ∈ [2112, 2212] && return true

    if p.N10 == 1 && p.N9 == 0
        # Charge should always be less than or equal to the baryon number
        A_ = A(p)
        Z_ = Z(p)

        (isnothing(A_) || isnothing(Z_)) && return false

        A_ >= abs(Z_) && return true
    end
    false
end

"""
    A(p::Union{Particle, PDGID, Integer})

Returns the atomic number A if the PDG ID corresponds to a nucleus, else it
returns `nothing`.
"""
function A(p)
    p = pdgid(p)
    abspdgid = abs(p.value)
    abspdgid ∈ [2112, 2212] && return 1
    (p.N10 != 1 || p.N9 != 0) && return nothing
    (abspdgid ÷ 10) % 1000
end

"""
    Z(p::Union{Particle, PDGID, Integer})

Returns the charge Z if the PDG ID corresponds to a nucleus, else it returns
`nothing`.
"""
function Z(p)
    p = pdgid(p)
    abspdgid = abs(p.value)
    abspdgid == 2212 && return sign(p.value)
    abspdgid == 2112 && return 0
    (p.N10 != 1 || p.N9 != 0) && return nothing
    ((abspdgid ÷ 10000) % 1000) * sign(p.value)
end

"""
    jspin(p::Union{Particle, PDGID, Integer})

Returns the total spin as 2J+1.
"""
function jspin(p)
    p = pdgid(p)
    !isvalid(p) && return nothing
    fid = fundamentalid(p)
    if fid > 0
        0 < fid < 7 && return 2  # 4th generation quarks not dealt with!
        fid == 9 && return 3
        10 < fid < 17 && return 2  # 4th generation quarks not deal with!
        20 < fid < 25 && return 3  # 4th generation quarks not deal with!
        return nothing
    end
    abspdgid = abs(p.value)
    abspdgid ∈ [1000000010, 1000010010]  && return 2  # neutrion, proton
    !isstandard(p) && return nothing
    p.value ∈ [130, 310] && return 1
    return abspdgid % 10
end

"""
    J(p::Union{Particle, PDGID, Integer})

Returns the total spin J.
"""
function J(p)
    jspin_ = jspin(p)
    isnothing(jspin_) && return nothing
    (jspin_ - 1) / 2
end

"""
    S(p::Union{Particle, PDGID, Integer})

Returns the spin S.

This is valid for mesons only. `nothing` is returned otherwise.
Mesons with PDGIDs of the kind 9XXXXXX (N=9) are not experimentally well-known
particles and `nothing` is returned too.
"""
function S(p)
    p = pdgid(p)
    !ismeson(p) && return nothing
    !isvalid(p) && return nothing

    abspdgid = abs(p.value)
    (abspdgid ÷ 1000000) % 10 == 9 && return nothing  # no knowledge so far

    nl = ((abspdgid) ÷ 10000) % 10
    js = abspdgid % 10

    !(js == 1 || js >= 3) && return 0

    if nl == 0
        js == 1 && return 0
        return 1
    end
    if nl == 1
        js == 1 && return 1
        return 0
    end
    if nl ∈ [2, 3]
        js >= 3 && return 1
        return 0
    end
    0
end

"""
    sspin(p::Union{Particle, PDGID, Integer})

Returns the spin S as 2S+1.

This is valid for mesons only. `nothing` is returned otherwise.
Mesons with PDGIDs of the kind 9XXXXXX (N=9) are not experimentally well-known
particles and `nothing` is returned too.
"""
function sspin(p)
    s = S(p)
    isnothing(s) && return nothing
    2s + 1
end

"""
    L(p::Union{Particle, PDGID, Integer})

Returns the orbital angular momentum L.

This is valid for mesons only. `nothing` is returned otherwise.
Mesons with PDGIDs of the kind 9XXXXXX (N=9) are not experimentally well-known
particles and `nothing` is returned too.
"""
function L(p)
    p = pdgid(p)
    !ismeson(p) && return nothing
    !isvalid(p) && return nothing

    abspdgid = abs(p.value)
    (abspdgid ÷ 1000000) % 10 == 9 && return nothing  # no knowledge so far

    nl = (abspdgid ÷ 10000) % 10
    js = abspdgid % 10

    if nl == 0
        js == 1 && return 0
        js == 3 && return 0
        js == 5 && return 1
        js == 7 && return 2
        js == 9 && return 3
    elseif nl == 1
        js == 1 && return 1
        js == 3 && return 1
        js == 5 && return 2
        js == 7 && return 3
        js == 9 && return 4
    elseif nl == 2
        js == 3 && return 1
        js == 5 && return 2
        js == 7 && return 3
        js == 9 && return 4
    elseif nl == 3
        js == 3 && return 2
        js == 5 && return 3
        js == 7 && return 4
        js == 9 && return 5
    end

    0
end

"""
    lspin(p::Union{Particle, PDGID, Integer})

Returns the orbital angular momentum L as 2S+1.

This is valid for mesons only. `nothing` is returned otherwise.
Mesons with PDGIDs of the kind 9XXXXXX (N=9) are not experimentally well-known
particles and `nothing` is returned too.
"""
function lspin(p)
    l = L(p)
    isnothing(l) && return nothing
    2l + 1
end

"""
    threecharge(p::Union{Particle, PDGID, Integer})

Returns 3 times the EM charge.
"""
function threecharge(p)
    p = pdgid(p)
    !isvalid(p) && return nothing
    abspdgid = abs(p.value)
    ch100 = [-1, 2, -1, 2, -1, 2, -1, 2, 0, 0, -3, 0, -3, 0, -3, 0, -3, 0, 0, 0,
             0, 0, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 3, 0, 0, 3, 0, 0, 0, 0, -1,
             0, 0, 0, 0, 0, 0, 0, 0, 0, 6, 3, 6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
             0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
             0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]

    fid = fundamentalid(p)

    charge = nothing

    if !isstandard(p)
        if isnucleus(p)  # ion
            Z_ = Z(p)
            isnothing(Z_) && return nothing
            return 3Z_
        elseif isQball(p)
            charge = 3 * ((abspdgid ÷ 10) % 10000)
        else
            # this should never be reached in the present numbering scheme
            # since extra bits exist only for Q-balls and nuclei
            return nothing
        end
    elseif isdyon(p)
        charge = 3 * ((abspdgid ÷ 10) % 1000)
        # this is half right
        # the charge sign will be changed below if pid < 0
        if p.Nl == 2
            charge = -charge
        end
    elseif 0 < fid <= 100  # use table from the SciKit-HEP particle project
        charge = ch100[fid]
        if abspdgid ∈ [1000017, 1000018, 1000034, 1000052, 1000053, 1000054]
            charge = 0
        end
        if abspdgid == 5100061 && abspdgid == 5100062
            charge = 6
        end
    elseif p.Nj == 0  # KL, KS, or undefined
        return 0
    elseif p.Nq1 == 0 || (isRhadron(p) && p.Nq1 == 9)  # mesons
        if p.Nq2 == 3 || p.Nq2 == 5
            charge = ch100[p.Nq3] - ch100[p.Nq2]
        else
            charge = ch100[p.Nq2] - ch100[p.Nq3]
        end
    elseif p.Nq3 == 0  # diquarks
        charge = ch100[p.Nq2] + ch100[p.Nq1]
    elseif isbaryon(p) || (isRhadron(p) && p.Nl == 9)  # baryons
        charge = ch100[p.Nq3] + ch100[p.Nq2] + ch100[p.Nq1]
    end

    if !isnothing(charge) && p.value < 0
        charge = -charge
    end
    charge
end

"""
    charge(p::Union{Particle, PDGID, Integer})

Returns the EM charge.
"""
function charge(p)
    c = threecharge(p)
    !isnothing(c) && return c / 3
    return nothing
end

"""
    _hasquark(p:::Union{Particle, PDGID, Integer}, q::Integer)

Note that q is always positive, so [1, 6] for Standard Model quarks
and [7, 8] for fourth-generation quarks.
"""
function _hasquark(p, q::Integer)
    p = pdgid(p)
    # Nuclei can also contain strange quarks,
    # cf. the definition of a nucleus PDG ID in isnucleus.
    # This check needs to be done first since isstandard is false for nuclei
    if isnucleus(p)
        q ∈ [1, 2]  && return true # Nuclei by construction contain up and down quarks
        if q == 3 && !(p.value ∈ [2112, 2212])
            p.N8 > 0 && return true
            return false
        end
    end

    !isstandard(p) && return false
    fundamentalid(p) > 0 && return false
    isdyon(p) && return false

    if isRhadron(p)
        _digits = digits(abs(p.value), pad=10)
        iz = 7
        for loc ∈ range(6; step=-1, stop=2)
            if _digits[loc] == 0
                iz = loc
            elseif loc == iz - 1
                # ignore squark or gluino
            else
                _digits[loc] == q && return true
            end
        end
        return false
    end

    (p.Nq3 == q || p.Nq2 == q || p.Nq1 == q) && return true

    ispentaquark(p) && (p.Nl == q || p.Nr == q) && return true
    false
end

"""
    hasdown(p::Union{Particle, PDGID, Integer})
"""
hasdown(p) = _hasquark(pdgid(p), 1)

"""
    hasdown(p::Union{Particle, PDGID, Integer})
"""
hasup(p) = _hasquark(pdgid(p), 2)

"""
    hasstrange(p::Union{Particle, PDGID, Integer})
"""
hasstrange(p) = _hasquark(pdgid(p), 3)

"""
    hascharm(p::Union{Particle, PDGID, Integer})
"""
hascharm(p) = _hasquark(pdgid(p), 4)

"""
    hasbottom(p::Union{Particle, PDGID, Integer})
"""
hasbottom(p) = _hasquark(pdgid(p), 5)

"""
    hastop(p::Union{Particle, PDGID, Integer})
"""
hastop(p) = _hasquark(pdgid(p), 6)