text:
"""
    HashSet(N::Integer)

Construct a `HashSet` with at least `N` slots. A `HashSet` contains `UInt` hash
values and is optimized for quick inserts and membership test.

A `HashSet` is a completely unsafe type internal type, not to be used unless
you know what you are doing.
"""
mutable struct HashSet
    data::Vector{UInt}
    mask::UInt
    len::Int
    haszero::Bool # 0 is interpreted as absence of hash, so we need this

    function HashSet(N::Integer)
        # Min 8 otherwise a size 4 will never repopulate and get stuck
        len = max(8, nextpow(2, N))
        new(zeros(UInt, len), (len - 1) % UInt, 0, false)
    end
end

function Base.show(io::IO, s::HashSet)
    print(io, "HashSet($(s.len) / $(Int(s.mask) + 1))")
end

Base.length(s::HashSet) = s.len

# Hashes are already uniformly distributed, so we just use hash as index then search
# linearly until we find the hash or an empty slot (0 is empty)
function Base.in(h::UInt, s::HashSet)
    pos = (h & s.mask) + 1
    @inbounds v = s.data[pos]
    @inbounds while !iszero(v)
        v == h && return true
        pos = (pos & s.mask) + 1
        v = s.data[pos]
    end
    return iszero(h) & s.haszero
end

function Base.empty!(s::HashSet)
    fill!(s.data, 0)
    s.len = 0
    s.haszero = false
    return s
end

# This simple structure do not support deletion, so when it gets too full, we simply
# empty it and then repopulate it from a vector containing all the hashes we need
function repopulate!(s::HashSet, h::Vector{UInt})
    empty!(s)
    @inbounds for i in h
        unsafe_push!(s, i)
    end
    s.len = length(h)
end

function unsafe_push!(s::HashSet, h::UInt)
    s.haszero |= iszero(h)
    pos = (h & s.mask) + 1
    @inbounds v = s.data[pos]
    @inbounds while !iszero(v)
        pos = (pos & s.mask) + 1
        v = s.data[pos]
    end
    s.len += 1
    @inbounds s.data[pos] = h
end

# The HashSet doesn't actually check if the hash is already in the set or not
# it just adds it - so only push when item is not already in set.
function Base.push!(s::HashSet, h::UInt, heap::Vector{UInt})
    s.len > (s.mask - s.mask >>> 2) && repopulate!(s, heap)
    unsafe_push!(s, h)
end


text:
function hqrt(t, b)
#
#  transformation using the t from hessqr
#
   z = copy(b)
   k = min(size(t,2), size(z,1)-1)
   for j = 1:k
       T = [t[1,j] t[2,j]; -conj(t[2,j]) conj(t[1,j])]
       z[j:j+1,1:end] = T * z[j:j+1,1:end]
   end
   z
end

function hqrt1(t::AbstractArray{T,2}, n::Int=1) where T
#
#  transformation using the t from hessqr on unit vector e(n)
#
   k = size(t,2)
   z = zeros(T, k+1)
   z[n] = one(T)
   for j = n:k
	   zj = z[j]
	   z[j] *= t[1,j]
	   z[j+1] = t[2,j] * zj
   end
   z
end

ext: jl


text:
using Random, SparseArrays

@testset "SparseMatrixCSC with matrices" begin
    rng = MersenneTwister(123)
    A = sprand(rng, 10, 5, 0.2)
    C = zeros(10, 2)
    B = randn(5, 2)

    cols = 1:5
    colsmul!(C, A, B, cols)
    @test C ≈ A[:, cols]*B[cols, :]

    cols = 2:4
    colsmul!(C, A, B, cols)
    @test C ≈ A[:, cols]*B[cols, :]

    cols = [1, 5]
    colsmul!(C, A, B, cols)
    @test C ≈ A[:, cols]*B[cols, :]

    α, β = rand(), rand()
    Cc = copy(C)
    colsmul!(C, A, B, cols, α, β)
    @test C ≈ α.*A[:, cols]*B[cols, :] .+ β.*Cc
end

@testset "Transposed SparseMatrixCSC with matrices" begin
    rng = MersenneTwister(123)
    A = sprand(rng, 10, 5, 0.2)
    C = zeros(5, 2)
    B = randn(10, 2)

    cols = 1:5
    tcolsmul!(C, A, B, cols)
    @test C[cols, :] ≈ A[:, cols]'*B

    cols = 2:4
    tcolsmul!(C, A, B, cols)
    @test C[cols, :] ≈ A[:, cols]'*B

    cols = [1, 5]
    tcolsmul!(C, A, B, cols)
    @test C[cols, :] ≈ A[:, cols]'*B

    α, β = rand(), rand()
    Cc = copy(C)
    tcolsmul!(C, A, B, cols, α, β)
    @test C[cols, :] ≈ α.*A[:, cols]'*B .+ β.*Cc[cols, :]
end

@testset "SparseMatrixCSC with vectors" begin
    rng = MersenneTwister(123)
    A = sprand(rng, 10, 5, 0.2)
    c = zeros(10)
    b = randn(5)

    cols = 1:5
    colsmul!(c, A, b, cols)
    @test c ≈ A[:, cols]*b[cols]

    cols = 2:4
    colsmul!(c, A, b, cols)
    @test c ≈ A[:, cols]*b[cols]

    cols = [1, 5]
    colsmul!(c, A, b, cols)
    @test c ≈ A[:, cols]*b[cols]

    α, β = rand(), rand()
    cc = copy(c)
    colsmul!(c, A, b, cols, α, β)
    @test c ≈ α.*A[:, cols]*b[cols] .+ β.*cc
end

@testset "Matrix with matrices" begin
    A = randn(10, 5)
    C = zeros(10, 2)
    B = randn(5, 2)

    cols = 1:5
    colsmul!(C, A, B, cols)
    @test C ≈ A[:, cols]*B[cols, :]

    cols = 2:4
    colsmul!(C, A, B, cols)
    @test C ≈ A[:, cols]*B[cols, :]

    cols = [1, 5]
    colsmul!(C, A, B, cols)
    @test C ≈ A[:, cols]*B[cols, :]

    α, β = rand(), rand()
    Cc = copy(C)
    colsmul!(C, A, B, cols, α, β)
    @test C ≈ α.*A[:, cols]*B[cols, :] .+ β.*Cc
end

@testset "Transposed Matrix with matrices" begin
    A = randn(10, 5)
    C = zeros(5, 2)
    B = randn(10, 2)

    cols = 1:5
    tcolsmul!(C, A, B, cols)
    @test C[cols, :] ≈ A[:, cols]'*B

    cols = 2:4
    tcolsmul!(C, A, B, cols)
    @test C[cols, :] ≈ A[:, cols]'*B

    cols = [1, 5]
    tcolsmul!(C, A, B, cols)
    @test C[cols, :] ≈ A[:, cols]'*B

    α, β = rand(), rand()
    Cc = copy(C)
    tcolsmul!(C, A, B, cols, α, β)
    @test C[cols, :] ≈ α.*A[:, cols]'*B .+ β.*Cc[cols, :]
end


text:
#=
BFSParentsExample:
- Julia version: 1.6.2
- Author: samuel
- Date: 2021-09-11
=#

using SuiteSparseGraphBLAS
using SparseArrays
using LinearAlgebra
using Random

include("/Users/samuel/Documents/Jupyter/GraphBlasAlgorithmsFinal/src/BFSParents.jl")
dim = 10000
sourcdim = 3
matrix =  sprand(Bool, dim, dim, 0.1)
matrix = GBMatrix(matrix)
source =  1
@time bfs_parents(matrix, source, dim)


text:
const Required = Union{String,Tuple{String,String}}

"""
    const rules = Dict([
        "article"       => ["author", "journal", "title", "year"]
        "book"          => [("author", "editor"), "publisher", "title", "year"]
        "booklet"       => ["title"]
        "eprint"        => ["author", "eprint", "title", "year"]
        "inbook"        => [("author", "editor"), ("chapter", "pages"), "publisher", "title", "year"]
        "incollection"  => ["author", "booktitle", "publisher", "title", "year"]
        "inproceedings" => ["author", "booktitle", "title", "year"]
        "manual"        => ["title"]
        "mastersthesis" => ["author", "school", "title", "year"]
        "misc"          => []
        "phdthesis"     => ["author", "school", "title", "year"]
        "proceedings"   => ["title", "year"]
        "techreport"    => ["author", "institution", "title", "year"]
        "unpublished"   => ["author", "note", "title"]
    ])
List of BibTeX rules bases on the entry type. A field value as a singleton represents a required field. A pair of values represents mutually exclusive required fields.
"""
const rules = Dict{String,Vector{Required}}(
    [
        "article" => ["author", "journal", "title", "year"]
        "book" => [("author", "editor"), "publisher", "title", "year"]
        "booklet" => ["title"]
        "eprint" => ["author", "eprint", "title", "year"]
        "inbook" =>
            [("author", "editor"), ("chapter", "pages"), "publisher", "title", "year"]
        "incollection" => ["author", "booktitle", "publisher", "title", "year"]
        "inproceedings" => ["author", "booktitle", "title", "year"]
        "manual" => ["title"]
        "mastersthesis" => ["author", "school", "title", "year"]
        "misc" => []
        "phdthesis" => ["author", "school", "title", "year"]
        "proceedings" => ["title", "year"]
        "techreport" => ["author", "institution", "title", "year"]
        "unpublished" => ["author", "note", "title"]
    ],
)

"""
    check_entry(fields::Fields)
Check the validity of the fields of a BibTeX entry.
"""
function check_entry(fields, check, id)
    errors = Vector{String}()

    entry_type = get(fields, "_type", "misc")
    if entry_type ∉ keys(rules)
        if check ∈ [:error, :warn]
            @warn """KeyError: key "software" not found in BibTeX rules, parsed from the entry `$id` with""" fields
        end
        check == :error && throw(KeyError(entry_type))
    end

    for t_field in get(rules, entry_type, Vector{Required}())
        at_least_one = false
        if typeof(t_field) == Tuple{String,String}
            for field in t_field
                if get(fields, field, "") != ""
                    at_least_one = true
                    break
                end
            end
            if !at_least_one
                push!(errors, "{" * foldl((x, y) -> "$x≡$y", t_field; init="")[2:end] * "}")
            end
        else
            get(fields, t_field, "") == "" && push!(errors, t_field)
        end
    end

    return errors
end

"""
    make_bibtex_entry(id::String, fields::Fields)
Make an entry if the entry follows the BibTeX guidelines. Throw an error otherwise.
"""
function make_bibtex_entry(id, fields; check=:error)
    # @info id fields
    "eprint" ∈ keys(fields) && (fields["_type"] = "eprint")
    fields = Dict(lowercase(k) => v for (k, v) in fields) # lowercase tag names
    errors = check_entry(fields, check, id)
    if length(errors) > 0 && check ∈ [:error, :warn]
        message =
            "Entry $id is missing the " *
            foldl(((x, y) -> x * ", " * y), errors) *
            " field(s)."
        check == :error ? (@error message) : (@warn message)
    end
    return Entry(id, fields)
end


text:
export JacobiWeight, WeightedJacobi

"""
    JacobiWeight(β,α,s::Space)

weights a space `s` by a Jacobi weight, which on `-1..1`
is `(1+x)^β*(1-x)^α`.
For other domains, the weight is inferred by mapping to `-1..1`.
"""
struct JacobiWeight{S,DD,RR,T<:Real} <: WeightSpace{S,DD,RR}
    β::T
    α::T
    space::S
    function JacobiWeight{S,DD,RR,T}(β::T,α::T,space::S) where {S<:Space,DD,RR,T}
        if isa(space,JacobiWeight)
            new(β+space.β,α+space.α,space.space) else
            new(β,α,space)
        end
    end
end

const WeightedJacobi{D,R} = JacobiWeight{Jacobi{D,R},D,R,R}

JacobiWeight{S,DD,RR,T}(β,α,space::Space) where {S,DD,RR,T} =
    JacobiWeight{S,DD,RR,T}(convert(T,β)::T, convert(T,α)::T, convert(S,space)::S)

JacobiWeight(a::Number, b::Number, d::Space) =
    JacobiWeight{typeof(d),domaintype(d),rangetype(d),promote_type(eltype(a),eltype(b))}(a,b,d)
JacobiWeight(β::Number, α::Number, d::JacobiWeight) =  JacobiWeight(β+d.β,α+d.α,d.space)
JacobiWeight(a::Number, b::Number, d::IntervalOrSegment) = JacobiWeight(a,b,Space(d))
JacobiWeight(a::Number, b::Number, d) = JacobiWeight(a,b,Space(d))
JacobiWeight(a::Number, b::Number) = JacobiWeight(a,b,Chebyshev())

JacobiWeight(a::Number, b::Number,s::PiecewiseSpace) = PiecewiseSpace(JacobiWeight(a,b,components(s)))

WeightedJacobi(β,α,d::Domain) = JacobiWeight(β,α,Jacobi(β,α,d))
WeightedJacobi(β,α) = JacobiWeight(β,α,Jacobi(β,α))

Fun(::typeof(identity), S::JacobiWeight) =
    isapproxinteger(S.β) && isapproxinteger(S.α) ? Fun(x->x,S) : Fun(identity,domain(S))

order(S::JacobiWeight{Ultraspherical{Int,D,R},D,R}) where {D,R} = order(S.space)

spacescompatible(A::JacobiWeight,B::JacobiWeight) =
    A.β ≈ B.β && A.α ≈ B.α && spacescompatible(A.space,B.space)
spacescompatible(A::JacobiWeight,B::Space{DD,RR}) where {DD<:IntervalOrSegment,RR<:Real} =
    spacescompatible(A,JacobiWeight(0,0,B))
spacescompatible(B::Space{DD,RR},A::JacobiWeight) where {DD<:IntervalOrSegment,RR<:Real} =
    spacescompatible(A,JacobiWeight(0,0,B))

transformtimes(f::Fun{JW1},g::Fun{JW2}) where {JW1<:JacobiWeight,JW2<:JacobiWeight}=
            Fun(JacobiWeight(f.space.β+g.space.β,f.space.α+g.space.α,f.space.space),
                coefficients(transformtimes(Fun(f.space.space,f.coefficients),
                                            Fun(g.space.space,g.coefficients))))
transformtimes(f::Fun{JW},g::Fun) where {JW<:JacobiWeight} =
    Fun(f.space,coefficients(transformtimes(Fun(f.space.space,f.coefficients),g)))
transformtimes(f::Fun,g::Fun{JW}) where {JW<:JacobiWeight} =
    Fun(g.space,coefficients(transformtimes(Fun(g.space.space,g.coefficients),f)))

jacobiweight(β,α,x) = -1 ≤ x ≤ 1 ? (1+x)^β*(1-x)^α : zero(x)
jacobiweight(β,α,d::Domain) = Fun(JacobiWeight(β,α,ConstantSpace(d)),[1.])
jacobiweight(β,α) = jacobiweight(β,α,ChebyshevInterval())

weight(sp::JacobiWeight,x) = jacobiweight(sp.β,sp.α,real(tocanonical(sp,x)))
dimension(sp::JacobiWeight) = dimension(sp.space)

Base.first(f::Fun{JW}) where {JW<:JacobiWeight} = space(f).β>0 ? zero(cfstype(f)) : f(leftendpoint(domain(f)))
Base.last(f::Fun{JW}) where {JW<:JacobiWeight} = space(f).α>0 ? zero(cfstype(f)) : f(rightendpoint(domain(f)))

setdomain(sp::JacobiWeight,d::Domain)=JacobiWeight(sp.β,sp.α,setdomain(sp.space,d))

# we assume that points avoids singularities

##TODO: paradigm for same space
function coefficients(f::AbstractVector,sp1::JacobiWeight{SJ1,DD},sp2::JacobiWeight{SJ2,DD}) where {SJ1,SJ2,DD<:IntervalOrSegment}
    β,α=sp1.β,sp1.α
    c,d=sp2.β,sp2.α

    if isapprox(c,β) && isapprox(d,α)
        # remove wrapper spaces and then convert
        coefficients(f,sp1.space,sp2.space)
    else
        # go back to default
        defaultcoefficients(f,sp1,sp2)
    end
end
coefficients(f::AbstractVector,sp::JacobiWeight{SJ,DD},
             S2::SubSpace{S,IT,DD,RR}) where {SJ,S,IT,DD<:IntervalOrSegment,RR<:Real} = subspace_coefficients(f,sp,S2)
coefficients(f::AbstractVector,S2::SubSpace{S,IT,DD,RR},
             sp::JacobiWeight{SJ,DD}) where {SJ,S,IT,DD<:IntervalOrSegment,RR<:Real} = subspace_coefficients(f,S2,sp)
#TODO: it could be possible that we want to JacobiWeight a SumSpace....
coefficients(f::AbstractVector,sp::JacobiWeight{SJ,DD},S2::SumSpace{SV,DD,RR}) where {SJ,SV,DD<:IntervalOrSegment,RR<:Real} =
    sumspacecoefficients(f,sp,S2)
coefficients(f::AbstractVector,sp::JacobiWeight{SJ,Segment{Vec{2,TT}}},S2::TensorSpace{SV,TTT,DD}) where {SJ,TT,SV,TTT,DD<:EuclideanDomain{2}} =
    coefficients(f,sp,JacobiWeight(0,0,S2))

coefficients(f::AbstractVector,sp::JacobiWeight{SJ,DD},S2::Space{DD,RR}) where {SJ,DD<:IntervalOrSegment,RR<:Real} =
    coefficients(f,sp,JacobiWeight(0,0,S2))
coefficients(f::AbstractVector,sp::ConstantSpace{DD},ts::JacobiWeight{SJ,DD}) where {SJ,DD<:IntervalOrSegment} =
    f.coefficients[1]*ones(ts).coefficients
coefficients(f::AbstractVector,S2::Space{DD,RR},sp::JacobiWeight{SJ,DD}) where {SJ,DD<:IntervalOrSegment,RR<:Real} =
    coefficients(f,JacobiWeight(0,0,S2),sp)

"""
`increase_jacobi_parameter(f)` multiplies by `1-x^2` on the unit interval.
`increase_jacobi_parameter(-1,f)` multiplies by `1+x` on the unit interval.
`increase_jacobi_parameter(+1,f)` multiplies by `1-x` on the unit interval.
On other domains this is accomplished by mapping to the unit interval.
"""
increase_jacobi_parameter(f) = Fun(f,JacobiWeight(f.space.β+1,f.space.α+1,space(f).space))
increase_jacobi_parameter(s,f) = s==-1 ? Fun(f,JacobiWeight(f.space.β+1,f.space.α,space(f).space)) :
                                       Fun(f,JacobiWeight(f.space.β,f.space.α+1,space(f).space))

function canonicalspace(S::JacobiWeight)
    if isapprox(S.β,0) && isapprox(S.α,0)
        canonicalspace(S.space)
    else
        #TODO: promote singularities?
        JacobiWeight(S.β,S.α,canonicalspace(S.space))
    end
end

function union_rule(A::ConstantSpace,B::JacobiWeight{P}) where P<:PolynomialSpace
    # we can convert to a space that contains contants provided
    # that the parameters are integers
    # when the parameters are -1 we keep them
    if isapproxinteger(B.β) && isapproxinteger(B.α)
        JacobiWeight(min(B.β,0.),min(B.α,0.),B.space)
    else
        NoSpace()
    end
end

## Algebra

function /(c::Number,f::Fun{JW}) where JW<:JacobiWeight
    g=c/Fun(space(f).space,f.coefficients)
    Fun(JacobiWeight(-f.space.β,-f.space.α,space(g)),g.coefficients)
end

function ^(f::Fun{JW}, k::AbstractFloat) where JW<:JacobiWeight
    S=space(f)
    g=Fun(S.space,coefficients(f))^k
    Fun(JacobiWeight(k*S.β,k*S.α,space(g)),coefficients(g))
end

function *(f::Fun{JW1},g::Fun{JW2}) where {JW1<:JacobiWeight,JW2<:JacobiWeight}
    @assert domainscompatible(f,g)
    fβ,fα=f.space.β,f.space.α
    gβ,gα=g.space.β,g.space.α
    m=(Fun(space(f).space,f.coefficients).*Fun(space(g).space,g.coefficients))
    if isapprox(fβ+gβ,0)&&isapprox(fα+gα,0)
        m
    else
        Fun(JacobiWeight(fβ+gβ,fα+gα,space(m)),m.coefficients)
    end
end

/(f::Fun{JW1},g::Fun{JW2}) where {JW1<:JacobiWeight,JW2<:JacobiWeight}=f*(1/g)

# O(min(m,n)) Ultraspherical conjugated inner product

function conjugatedinnerproduct(sp::Ultraspherical,u::AbstractVector{S},v::AbstractVector{V}) where {S,V}
    λ=order(sp)
    if λ==1
        mn = min(length(u),length(v))
        if mn > 0
            return dotu(u[1:mn],v[1:mn])*π/2
        else
            return zero(promote_type(eltype(u),eltype(v)))
        end
    else
        T,mn = promote_type(S,V),min(length(u),length(v))
        if mn > 1
            wi = sqrt(convert(T,π))*gamma(λ+one(T)/2)/gamma(λ+one(T))
            ret = u[1]*wi*v[1]
            for i=2:mn
              wi *= (i-2one(T)+2λ)/(i-one(T)+λ)*(i-2one(T)+λ)/(i-one(T))
              ret += u[i]*wi*v[i]
            end
            return ret
        elseif mn > 0
            wi = sqrt(convert(T,π))*gamma(λ+one(T)/2)/gamma(λ+one(T))
            return u[1]*wi*v[1]
        else
            return zero(promote_type(eltype(u),eltype(v)))
        end
    end
end

function conjugatedinnerproduct(::Chebyshev,u::AbstractVector,v::AbstractVector)
    mn = min(length(u),length(v))
    if mn > 1
        return (2u[1]*v[1]+dotu(u[2:mn],v[2:mn]))*π/2
    elseif mn > 0
        return u[1]*v[1]*π
    else
        return zero(promote_type(eltype(u),eltype(v)))
    end
end

function bilinearform(f::Fun{JacobiWeight{Ultraspherical{LT,D,R},D,R,TT}},g::Fun{Ultraspherical{LT,D,R}}) where {LT,D,R,TT}
    d = domain(f)
    @assert d == domain(g)
    λ = order(space(f).space)
    if order(space(g)) == λ && f.space.β == f.space.α == λ-0.5
        return complexlength(d)/2*conjugatedinnerproduct(Ultraspherical(λ,d),f.coefficients,g.coefficients)
    else
        return defaultbilinearform(f,g)
    end
end

function bilinearform(f::Fun{Ultraspherical{LT,D,R}},
                    g::Fun{JacobiWeight{Ultraspherical{LT,D,R},D,R,TT}}) where {LT,D,R,TT}
    d = domain(f)
    @assert d == domain(g)
    λ = order(space(f))
    if order(space(g).space) == λ && g.space.β == g.space.α == λ-0.5
        return complexlength(d)/2*conjugatedinnerproduct(Ultraspherical(λ,d),f.coefficients,g.coefficients)
    else
        return defaultbilinearform(f,g)
    end
end

function bilinearform(f::Fun{JacobiWeight{Ultraspherical{LT,D,R},D,R,TT1}},
                    g::Fun{JacobiWeight{Ultraspherical{LT,D,R},D,R,TT2}}) where {LT,D,R,TT1,TT2}
    d = domain(f)
    @assert d == domain(g)
    λ = order(space(f).space)
    if order(space(g).space) == λ && f.space.β+g.space.β == f.space.α+g.space.α == λ-0.5
        return complexlength(domain(f))/2*conjugatedinnerproduct(Ultraspherical(λ,d),f.coefficients,g.coefficients)
    else
        return defaultbilinearform(f,g)
    end
end

function linebilinearform(f::Fun{JacobiWeight{Ultraspherical{LT,D,R},D,R,TT}},g::Fun{Ultraspherical{LT,D,R}}) where {LT,D,R,TT}
    d = domain(f)
    @assert d == domain(g)
    λ = order(space(f).space)
    if order(space(g)) == λ && f.space.β == f.space.α == λ-0.5
        return arclength(d)/2*conjugatedinnerproduct(Ultraspherical(λ,d),f.coefficients,g.coefficients)
    else
        return defaultlinebilinearform(f,g)
    end
end

function linebilinearform(f::Fun{Ultraspherical{LT,D,R}},g::Fun{JacobiWeight{Ultraspherical{LT,D,R},D,R,TT}}) where {LT,D,R,TT}
    d = domain(f)
    @assert d == domain(g)
    λ = order(space(f))
    if order(space(g).space) == λ &&  g.space.β == g.space.α == λ-0.5
        return arclength(d)/2*conjugatedinnerproduct(Ultraspherical(λ,d),f.coefficients,g.coefficients)
    else
        return defaultlinebilinearform(f,g)
    end
end

function linebilinearform(f::Fun{JacobiWeight{Ultraspherical{LT,D,R},D,R,TT1}},g::Fun{JacobiWeight{Ultraspherical{LT,D,R},D,R,TT2}}) where {LT,D,R,TT1,TT2}
    d = domain(f)
    @assert d == domain(g)
    λ = order(space(f).space)
    if order(space(g).space) == λ &&  f.space.β+g.space.β == f.space.α+g.space.α == λ-0.5
        return arclength(d)/2*conjugatedinnerproduct(Ultraspherical(λ,d),f.coefficients,g.coefficients)
    else
        return defaultlinebilinearform(f,g)
    end
end

function bilinearform(f::Fun{JacobiWeight{J,DD,RR,TT}},g::Fun{J}) where {J<:Jacobi,DD<:IntervalOrSegment,RR,TT}
    @assert domain(f) == domain(g)
    if f.space.β == f.space.space.a == g.space.a && f.space.α == f.space.space.b == g.space.b
        return complexlength(domain(f))/2*conjugatedinnerproduct(g.space,f.coefficients,g.coefficients)
    else
        return defaultbilinearform(f,g)
    end
end

function bilinearform(f::Fun{J},
                      g::Fun{JacobiWeight{J,DD,RR,TT}}) where {J<:Jacobi,DD<:IntervalOrSegment,RR,TT}
    @assert domain(f) == domain(g)
    if g.space.β == g.space.space.a == f.space.a && g.space.α == g.space.space.b == f.space.b
        return complexlength(domain(f))/2*conjugatedinnerproduct(f.space,f.coefficients,g.coefficients)
    else
        return defaultbilinearform(f,g)
    end
end

function bilinearform(f::Fun{JacobiWeight{J,DD,RR,TT}},
                      g::Fun{JacobiWeight{J,DD,RR,TT}}) where {J<:Jacobi,DD<:IntervalOrSegment,RR,TT}
    @assert domain(f) == domain(g)
    if f.space.β + g.space.β == f.space.space.a == g.space.space.a && f.space.α + g.space.α == f.space.space.b == g.space.space.b
        return complexlength(domain(f))/2*conjugatedinnerproduct(f.space.space,f.coefficients,g.coefficients)
    else
        return defaultbilinearform(f,g)
    end
end

function linebilinearform(f::Fun{JacobiWeight{J,DD,RR,TT}},g::Fun{J}) where {J<:Jacobi,DD<:IntervalOrSegment,RR,TT}
    @assert domain(f) == domain(g)
    if f.space.β == f.space.space.a == g.space.a && f.space.α == f.space.space.b == g.space.b
        return arclength(domain(f))/2*conjugatedinnerproduct(g.space,f.coefficients,g.coefficients)
    else
        return defaultlinebilinearform(f,g)
    end
end

function linebilinearform(f::Fun{J},g::Fun{JacobiWeight{J,DD,RR,TT}}) where {J<:Jacobi,DD<:IntervalOrSegment,RR,TT}
    @assert domain(f) == domain(g)
    if g.space.β == g.space.space.a == f.space.a && g.space.α == g.space.space.b == f.space.b
        return arclength(domain(f))/2*conjugatedinnerproduct(f.space,f.coefficients,g.coefficients)
    else
        return defaultlinebilinearform(f,g)
    end
end

function linebilinearform(f::Fun{JacobiWeight{J,DD,RR,TT}}, g::Fun{JacobiWeight{J,DD,RR,TT}}) where {J<:Jacobi,DD<:IntervalOrSegment,RR,TT}
    @assert domain(f) == domain(g)
    if f.space.β + g.space.β == f.space.space.a == g.space.space.a && f.space.α + g.space.α == f.space.space.b == g.space.space.b
        return arclength(domain(f))/2*conjugatedinnerproduct(f.space.space,f.coefficients,g.coefficients)
    else
        return defaultlinebilinearform(f,g)
    end
end

function Derivative(S::WeightedJacobi{DDD,RR}) where {DDD<:IntervalOrSegment,RR}
    if S.β>0 && S.β>0 && S.β==S.space.b && S.α==S.space.a
        ConcreteDerivative(S,1)
    else
        jacobiweightDerivative(S)
    end
end

bandwidths(D::ConcreteDerivative{WeightedJacobi{DDD,RR}}) where {DDD<:IntervalOrSegment,RR} = 1,0
rangespace(D::ConcreteDerivative{WeightedJacobi{DDD,RR}}) where {DDD<:IntervalOrSegment,RR} =
    WeightedJacobi(domainspace(D).β-1,domainspace(D).α-1,domain(D))

getindex(D::ConcreteDerivative{WeightedJacobi{DDD,RR}},k::Integer,j::Integer) where {DDD<:IntervalOrSegment,RR} =
    j==k-1 ? eltype(D)(-4(k-1)./complexlength(domain(D))) : zero(eltype(D))

for (Func,Len,Sum) in ((:DefiniteIntegral,:complexlength,:sum),(:DefiniteLineIntegral,:arclength,:linesum))
    ConcFunc = Meta.parse("Concrete"*string(Func))

    @eval begin
        function getindex(Σ::$ConcFunc{JacobiWeight{Jacobi{D,R},D,R,TT},T},k::Integer) where {D<:IntervalOrSegment,R,T,TT}
            dsp = domainspace(Σ)

            if dsp.β == dsp.space.b && dsp.α == dsp.space.a
                # TODO: copy and paste
                k == 1 ? convert(T,$Sum(Fun(dsp,[one(T)]))) : zero(T)
            else
                convert(T,$Sum(Fun(dsp,[zeros(T,k-1);1])))
            end
        end

        function bandwidths(Σ::$ConcFunc{JacobiWeight{Jacobi{D,R},D,R,TT}}) where {D<:IntervalOrSegment,R,TT}
            β,α = domainspace(Σ).β,domainspace(Σ).α
            if domainspace(Σ).β == domainspace(Σ).space.b && domainspace(Σ).α == domainspace(Σ).space.a
                0,0  # first entry
            else
                0,∞
            end
        end
    end
end

## <: IntervalOrSegment avoids a julia bug
function Multiplication(f::Fun{JacobiWeight{C,DD,RR,TT}}, S::Jacobi) where {C<:ConstantSpace,DD<:IntervalOrSegmentDomain,RR,TT}
    # this implements (1+x)*P and (1-x)*P special case
    # see DLMF (18.9.6)
    d=domain(f)
    if ((space(f).β==1 && space(f).α==0 && S.b >0) ||
                        (space(f).β==0 && space(f).α==1 && S.a >0))
        ConcreteMultiplication(f,S)
    elseif isapproxinteger(space(f).β) && space(f).β ≥ 1 && S.b >0
        # decrement β and multiply again
        M=Multiplication(f.coefficients[1]*jacobiweight(1.,0.,d),S)
        MultiplicationWrapper(f,Multiplication(jacobiweight(space(f).β-1,space(f).α,d),rangespace(M))*M)
    elseif isapproxinteger(space(f).α) && space(f).α ≥ 1 && S.a >0
        # decrement α and multiply again
        M=Multiplication(f.coefficients[1]*jacobiweight(0.,1.,d),S)
        MultiplicationWrapper(f,Multiplication(jacobiweight(space(f).β,space(f).α-1,d),rangespace(M))*M)
    else
# default JacobiWeight
        M=Multiplication(Fun(space(f).space,f.coefficients),S)
        rsp=JacobiWeight(space(f).β,space(f).α,rangespace(M))
        MultiplicationWrapper(f,SpaceOperator(M,S,rsp))
    end
end

Multiplication(f::Fun{JacobiWeight{C,DD,RR,TT}},
               S::Union{Ultraspherical,Chebyshev}) where {C<:ConstantSpace,DD<:IntervalOrSegmentDomain,RR,TT} =
    MultiplicationWrapper(f,Multiplication(f,Jacobi(S))*Conversion(S,Jacobi(S)))

function rangespace(M::ConcreteMultiplication{JacobiWeight{C,DD,RR,TT},J}) where {J<:Jacobi,C<:ConstantSpace,DD<:IntervalOrSegmentDomain,RR,TT}
    S=domainspace(M)
    if space(M.f).β==1
        # multiply by (1+x)
        Jacobi(S.b-1,S.a,domain(S))
    elseif space(M.f).α == 1
        # multiply by (1-x)
        Jacobi(S.b,S.a-1,domain(S))
    else
        error("Not implemented")
    end
end

bandwidths(::ConcreteMultiplication{JacobiWeight{C,DD,RR,TT},J}) where {J<:Jacobi,C<:ConstantSpace,DD<:IntervalOrSegmentDomain,RR,TT} = 1,0

function getindex(M::ConcreteMultiplication{JacobiWeight{C,DD,RR,TT},J},k::Integer,j::Integer) where {J<:Jacobi,C<:ConstantSpace,DD<:IntervalOrSegmentDomain,RR,TT}
    @assert ncoefficients(M.f)==1
    a,b=domainspace(M).a,domainspace(M).b
    c=M.f.coefficients[1]
    if space(M.f).β==1
        @assert space(M.f).α==0
        # multiply by (1+x)
        if j==k
            c*2(k+b-1)/(2k+a+b-1)
        elseif k > 1 && j==k-1
            c*(2k-2)/(2k+a+b-3)
        else
            zero(eltype(M))
        end
    elseif space(M.f).α == 1
        @assert space(M.f).β==0
        # multiply by (1-x)
        if j==k
            c*2(k+a-1)/(2k+a+b-1)
        elseif k > 1 && j==k-1
            -c*(2k-2)/(2k+a+b-3)
        else
            zero(eltype(M))
        end
    else
        error("Not implemented")
    end
end

# We can exploit the special multiplication to construct a Conversion

for FUNC in (:maxspace_rule,:union_rule,:hasconversion)
    @eval function $FUNC(A::WeightedJacobi{DD},B::Jacobi) where DD<:IntervalOrSegment
        if A.β==A.α+1 && A.space.b>0
            $FUNC(Jacobi(A.space.b-1,A.space.a,domain(A)),B)
        elseif A.α==A.β+1 && A.space.a>0
            $FUNC(Jacobi(A.space.b,A.space.a-1,domain(A)),B)
        else
            $FUNC(A,JacobiWeight(0.,0.,B))
        end
    end
end


text:
export stirling2, stirling2_partitions, stirling2_matrices
export permutation_matrices, n_choose_k_matrices
export base_n_val

"""
    stirling2( n :: Int64, k :: Int64  ) :: Int64

Counts the number of ways to partition `n` items into `k` unlabelled groups.
This quantity is known as Stirling's number of the 2nd kind:

```math
\\left\\{n \\atop k \\right\\} = \\frac{1}{k!}\\sum_{i=0}^k (-1)^i\\binom{k}{i}(k-i)^n
```

Throws a `DomainError` if inputs do not satisfy `n ≥ k ≥ 1`.
"""
function stirling2(n :: Int64, k :: Int64) :: Int64
    if !(n ≥ k ≥ 1)
        throw(DomainError((n,k), "Inputs (n,k) do not satisfy `n ≥ k ≥ 1`."))
    end

    sum(map( i -> (-1)^i*binomial(k, i)*(k-i)^n, 0:k) )/factorial(k)
end

"""
    stirling2_partitions( n :: Int64, k :: Int64 ) :: Vector{Vector{Vector{Int64}}}

Enumerates the unique partitions of `n` items into `k` unlabelled sets.
Each partition is a vector containing a  set of `k` vectors designating each group.

E.g.

```jldoctest
julia> stirling2_partitions( 4, 2 )
7-element Array{Array{Array{Int64,1},1},1}:
 [[1, 2, 3], [4]]
 [[3], [1, 2, 4]]
 [[1, 2], [3, 4]]
 [[1, 3], [2, 4]]
 [[2], [1, 3, 4]]
 [[2, 3], [1, 4]]
 [[1], [2, 3, 4]]
```

This recursive algorithm was inspired by [this blog](https://devblogs.microsoft.com/oldnewthing/20140324-00/?p=1413).
"""
function stirling2_partitions(n :: Int64, k :: Int64) :: Vector{Vector{Vector{Int64}}}
    if (n ==  0) || (k == 0) || (n < k)
        return Vector{Vector{Vector{Int64}}}(undef,0)
    elseif (k == 1)
        return [[[1:n...]]]
    else
        # partitions include all (n-1, k-1) paritions with [n] as k-th group
        partitions1 = _stirling2_add_partition(n, stirling2_partitions(n-1, k-1))

        # partitions include all (n-1,  k) partitions with [n] appended to each existing group
        partitions2 = _stirling2_extend_partitions(n, k, stirling2_partitions(n-1, k))

        return cat(partitions1, partitions2, dims=1)
    end
end

# helper function for adding a new group to a partition
function _stirling2_add_partition(n :: Int64, partitions :: Vector{Vector{Vector{Int64}}}) :: Vector{Vector{Vector{Int64}}}
    if length(partitions) == 0
        return [[[n]]] :: Vector{Vector{Vector{Int64}}}
    else
        return map(s -> cat(s, [[n]], dims=1), partitions) :: Vector{Vector{Vector{Int64}}}
    end
end

# helper function for extending existing groups with a new element
function _stirling2_extend_partitions(n :: Int64, k :: Int64, partitions :: Vector{Vector{Vector{Int64}}}) :: Vector{Vector{Vector{Int64}}}
    new_partitions = Vector{Vector{Vector{Int64}}}(undef, length(partitions)*k)
    id = 1

    _extend_partition_set(partition :: Vector{Vector{Int64}}) = begin
        for i in 1:k
            new_partitions[id] = cat(partition[1:k .!= i], [cat(partition[i], [n], dims=1)], dims=1)
            id += 1
        end
    end

    map(_extend_partition_set, partitions)

    new_partitions
end

"""
    stirling2_matrices( n :: Int64, k :: Int64 ) :: Vector{Matrix{Bool}}

Generates the set of matrices with `k` rows and `n` columns where rows correspond
to the groups and columns are the grouped elements. A non-zero element designates
that the column id is grouped into the corresponding row.

E.g.

```jldoctest
julia> stirling2_matrices( 4, 2 )
7-element Array{Array{Bool,2},1}:
 [1 1 1 0; 0 0 0 1]
 [0 0 1 0; 1 1 0 1]
 [1 1 0 0; 0 0 1 1]
 [1 0 1 0; 0 1 0 1]
 [0 1 0 0; 1 0 1 1]
 [0 1 1 0; 1 0 0 1]
 [1 0 0 0; 0 1 1 1]
```

A `DomainError` is thrown if `n ≥ k ≥ 1` is not satisfied.
"""
function stirling2_matrices(n :: Int64, k :: Int64) :: Vector{Matrix{Bool}}
    if !(n ≥ k ≥ 1)
        throw(DomainError((n,k), "Inputs (n,k) do not satisfy `n ≥ k ≥ 1`."))
    end

    _partition_to_matrix(partition :: Vector{Vector{Int64}}) :: Matrix{Bool} = begin
        m = fill(false, k, n)
        for row_id in 1:k
            col_ids = partition[row_id]
            m[row_id,col_ids] .= true
        end

        return m
    end

    map(_partition_to_matrix, stirling2_partitions(n, k))
end

"""
    permutation_matrices( dim :: Int64 ) :: Vector{Matrix{Bool}}

Generates the set of square permutation matrices of dimension `dim`.

E.g.

```jldoctest
julia> permutation_matrices( 3 )
6-element Array{Array{Bool,2},1}:
 [1 0 0; 0 1 0; 0 0 1]
 [1 0 0; 0 0 1; 0 1 0]
 [0 1 0; 1 0 0; 0 0 1]
 [0 0 1; 1 0 0; 0 1 0]
 [0 1 0; 0 0 1; 1 0 0]
 [0 0 1; 0 1 0; 1 0 0]
```
"""
function permutation_matrices(dim :: Int64) :: Vector{Matrix{Bool}}
    map(perm_ids -> Matrix{Bool}(I,dim,dim)[:,perm_ids], permutations(1:dim))
end

"""
    n_choose_k_matrices( n :: Int64, k :: Int64 ) :: Vector{Matrix{Bool}}

Generates a set of `n` by `k` matrices which represent all combinations of selecting
`k` columns from `n` rows.  Each column, contains a single non-zero element and
`k` rows contain a non-zero element.

E.g.

```jldoctest
julia> n_choose_k_matrices( 4, 2 )
6-element Array{Array{Bool,2},1}:
 [1 0; 0 1; 0 0; 0 0]
 [1 0; 0 0; 0 1; 0 0]
 [1 0; 0 0; 0 0; 0 1]
 [0 0; 1 0; 0 1; 0 0]
 [0 0; 1 0; 0 0; 0 1]
 [0 0; 0 0; 1 0; 0 1]
```

A `DomainError` is thrown if `n ≥ k ≥ 1` is not satisfied.
"""
function n_choose_k_matrices(n :: Int64, k :: Int64) :: Vector{Matrix{Bool}}
    if !(n ≥ k ≥ 1)
        throw(DomainError((n,k), "Inputs (n,k) must satisfy `n ≥ k ≥ 1`."))
    end
    _combo_to_matrix(combo :: Vector{Int64}) :: Matrix{Bool} = begin
        m = fill(false,n,k)
        for col_id in 1:k
            m[combo[col_id],col_id] = true
        end

        return m
    end

    map(_combo_to_matrix, combinations(1:n, k))
end

"""
    base_n_val(
        num_array :: Vector{Int64},
        base :: Int64;
        big_endian=true :: Bool
    ) :: Int64

Given an array representing a number in base-n returns the value of that
number in base-10.

Inputs:
* `num_array` - Vector containing semi-positive integers less than base.
* `base` - The base-n number represented by num_array.
* `big_endian` - `true` if most significant place is at index 1, else `false`.
"""
function base_n_val(num_array :: Vector{Int64}, base :: Int64; big_endian = true :: Bool) :: Int64
    if maximum(num_array) >= base
        throw(ArgumentError("max value of $num_array is not valid in $base"))
    elseif minimum(num_array) < 0
        throw(ArgumentError("min value of $num_array must be zero or greater"))
    end

    pad = length(num_array)

    place_vals = big_endian ? map(
        x -> base^(pad - x[1]), enumerate(num_array)
    ) : map(
        x -> base^(x[1] - 1), enumerate(num_array)
    )

    val = place_vals'*num_array

    val
end


text:
"""
    dmrg(H::MPO,psi0::MPS,sweeps::Sweeps;kwargs...)

Use the density matrix renormalization group (DMRG) algorithm
to optimize a matrix product state (MPS) such that it is the
eigenvector of lowest eigenvalue of a Hermitian matrix H,
represented as a matrix product operator (MPO).
The MPS `psi0` is used to initialize the MPS to be optimized,
and the `sweeps` object determines the parameters used to
control the DMRG algorithm.

Returns:
* `energy::Float64` - eigenvalue of the optimized MPS
* `psi::MPS` - optimized MPS
"""
function dmrg(H::MPO,
              psi0::MPS,
              sweeps::Sweeps;
              kwargs...)
  PH = ProjMPO(H)
  return dmrg(PH,psi0,sweeps;kwargs...)
end

"""
    dmrg(Hs::Vector{MPO},psi0::MPS,sweeps::Sweeps;kwargs...)

Use the density matrix renormalization group (DMRG) algorithm
to optimize a matrix product state (MPS) such that it is the
eigenvector of lowest eigenvalue of a Hermitian matrix H.
The MPS `psi0` is used to initialize the MPS to be optimized,
and the `sweeps` object determines the parameters used to
control the DMRG algorithm.

This version of `dmrg` accepts a representation of H as a
Vector of MPOs, Hs = [H1,H2,H3,...] such that H is defined
as H = H1+H2+H3+...
Note that this sum of MPOs is not actually computed; rather
the set of MPOs [H1,H2,H3,..] is efficiently looped over at
each step of the DMRG algorithm when optimizing the MPS.

Returns:
* `energy::Float64` - eigenvalue of the optimized MPS
* `psi::MPS` - optimized MPS
"""
function dmrg(Hs::Vector{MPO},
              psi0::MPS,
              sweeps::Sweeps;
              kwargs...)
  PHS = ProjMPOSum(Hs)
  return dmrg(PHS,psi0,sweeps;kwargs...)
end

"""
    dmrg(H::MPO,Ms::Vector{MPS},psi0::MPS,sweeps::Sweeps;kwargs...)

Use the density matrix renormalization group (DMRG) algorithm
to optimize a matrix product state (MPS) such that it is the
eigenvector of lowest eigenvalue of a Hermitian matrix H,
subject to the constraint that the MPS is orthogonal to each
of the MPS provided in the Vector `Ms`. The orthogonality
constraint is approximately enforced by adding to H terms of
the form w|M1><M1| + w|M2><M2| + ... where Ms=[M1,M2,...] and
w is the "weight" parameter, which can be adjusted through the
optional `weight` keyword argument.
The MPS `psi0` is used to initialize the MPS to be optimized,
and the `sweeps` object determines the parameters used to
control the DMRG algorithm.

Returns:
* `energy::Float64` - eigenvalue of the optimized MPS
* `psi::MPS` - optimized MPS
"""
function dmrg(H::MPO,
              Ms::Vector{MPS},
              psi0::MPS,
              sweeps::Sweeps;
              kwargs...)
  weight = get(kwargs,:weight,1.0)
  PMM = ProjMPO_MPS(H,Ms;weight=weight)
  return dmrg(PMM,psi0,sweeps;kwargs...)
end

function dmrg(PH,
              psi0::MPS,
              sweeps::Sweeps;
              kwargs...)
  # Debug level checks
  @debug begin
    checkflux(psi0)
    checkflux(PH)
  end

  which_decomp::Union{String, Nothing} = get(kwargs, :which_decomp, nothing)
  svd_alg::String = get(kwargs, :svd_alg, "recursive")
  obs = get(kwargs, :observer, NoObserver())
  outputlevel::Int = get(kwargs, :outputlevel, 1)

  # eigsolve kwargs
  eigsolve_tol::Float64   = get(kwargs, :eigsolve_tol, 1e-14)
  eigsolve_krylovdim::Int = get(kwargs, :eigsolve_krylovdim, 3)
  eigsolve_maxiter::Int   = get(kwargs, :eigsolve_maxiter, 1)
  eigsolve_verbosity::Int = get(kwargs, :eigsolve_verbosity, 0)

  # TODO: add support for non-Hermitian DMRG
  # get(kwargs, :ishermitian, true)
  ishermitian::Bool = true

  # TODO: add support for targeting other states with DMRG
  # (such as the state with the largest eigenvalue)
  # get(kwargs, :eigsolve_which_eigenvalue, :SR)
  eigsolve_which_eigenvalue::Symbol = :SR

  # Keyword argument deprecations
  if haskey(kwargs, :maxiter)
    error("""maxiter keyword has been replaced by eigsolve_krylovdim.
             Note: compared to the C++ version of ITensor,
             setting eigsolve_krylovdim 3 is the same as setting
             a maxiter of 2.""")
  end

  if haskey(kwargs, :errgoal)
    error("errgoal keyword has been replaced by eigsolve_tol.")
  end

  if haskey(kwargs, :quiet)
    error("quiet keyword has been replaced by outputlevel")
  end

  psi = copy(psi0)
  N = length(psi)

  position!(PH, psi0, 1)
  energy = 0.0

  for sw=1:nsweep(sweeps)
    sw_time = @elapsed begin

    for (b, ha) in sweepnext(N)

      @debug begin
        checkflux(psi)
        checkflux(PH)
      end

@timeit_debug GLOBAL_TIMER "position!" begin
      position!(PH, psi, b)
end

      @debug begin
        checkflux(psi)
        checkflux(PH)
      end

@timeit_debug GLOBAL_TIMER "psi[b]*psi[b+1]" begin
      phi = psi[b] * psi[b+1]
end

@timeit_debug GLOBAL_TIMER "eigsolve" begin
      vals, vecs = eigsolve(PH, phi, 1, eigsolve_which_eigenvalue;
                            ishermitian = ishermitian,
                            tol = eigsolve_tol,
                            krylovdim = eigsolve_krylovdim,
                            maxiter = eigsolve_maxiter)
end
      energy, phi = vals[1], vecs[1]

      ortho = ha == 1 ? "left" : "right"

      drho = nothing
      if noise(sweeps, sw) > 0.0
        # Use noise term when determining new MPS basis
        drho = noise(sweeps, sw) * noiseterm(PH,phi,ortho)
      end

      @debug begin
        checkflux(phi)
      end

@timeit_debug GLOBAL_TIMER "replacebond!" begin
      spec = replacebond!(psi, b, phi; maxdim = maxdim(sweeps, sw),
                                       mindim = mindim(sweeps, sw),
                                       cutoff = cutoff(sweeps, sw),
                                       eigen_perturbation = drho,
                                       ortho = ortho,
                                       normalize = true,
                                       which_decomp = which_decomp,
                                       svd_alg = svd_alg)
end

      @debug begin
        checkflux(psi)
        checkflux(PH)
      end

      if outputlevel >= 2
        @printf("Sweep %d, half %d, bond (%d,%d) energy=%.12f\n",sw,ha,b,b+1,energy)
        @printf("(Truncated using cutoff=%.1E maxdim=%d mindim=%d)\n",
                cutoff(sweeps, sw),maxdim(sweeps, sw),mindim(sweeps, sw))
        @printf("Trunc. err=%.1E, bond dimension %d\n\n",spec.truncerr,dim(linkind(psi,b)))
      end

      measure!(obs; energy = energy,
                    psi = psi,
                    bond = b,
                    sweep = sw,
                    half_sweep = ha,
                    spec = spec,
                    outputlevel = outputlevel)
    end
    end
    if outputlevel >= 1
      @printf("After sweep %d energy=%.12f maxlinkdim=%d time=%.3f\n",
              sw, energy, maxlinkdim(psi), sw_time)
    end
    isdone = checkdone!(obs;energy=energy,
                            psi=psi,
                            sweep=sw,
                            outputlevel=outputlevel)

    isdone && break
  end
  return (energy, psi)
end


text:
"""
    $(SIGNATURES)

Calculate the ΔG'⁰ of formation for all `mets`, as well as the covariance matrix of the
errors.

# Example
```
equilibrator = Equilibrator()

mets = bigg.(["atp", "adp", "amp", "pi", "h2o"])

standard_dgf_prime_mu, standard_dgf_cov = standard_dg_formation(eq, mets)
```
"""
function standard_dg_formation(equilibrator::Equilibrator, mets::Vector{String})

    pymets = [equilibrator.cc.get_compound(cname) for cname in mets]

    met_formations = [equilibrator.cc.standard_dg_formation(met) for met in pymets]
    dgs = [met_formation[1] for met_formation in met_formations]

    sigmas_fin = hcat([met_formation[2] for met_formation in met_formations]...)
    sigmas_inf = hcat([met_formation[3] for met_formation in met_formations]...)

    leg_trans = [
        pymet.transform(
            equilibrator.cc.p_h,
            equilibrator.cc.ionic_strength,
            equilibrator.cc.temperature,
            equilibrator.cc.p_mg,
        ).m_as(
            "kJ/mol",
        ) for pymet in pymets
    ]

    standard_dgf_prime_mu = dgs + leg_trans
    standard_dgf_cov = sigmas_fin' * sigmas_fin + 1e6 .* sigmas_inf' * sigmas_inf

    return (standard_dgf_prime_mu)u"kJ/mol", (standard_dgf_cov)u"kJ/mol"
end


text:
using Test: @testset

@testset "Tests" begin
    @testset "Single Player: Forward and Inverse Optimal Control" begin
        include("lq_control.jl")
        include("nonlq_control.jl")
    end

    @testset "Multi Player: Forward and Inverse Infinite Dynamnic Games" begin
        include("shared_unicycle_game.jl")
        include("collision_avoidance_game.jl")
    end
end


text:
"""
   IO functions for the HARPS-N spectrograph
    https://plone.unige.ch/HARPS-N/
Author: Alex Wise and collaborators
Created: April 2021
"""

"""
Create Dataframe containing filenames and key data for all files YYYY-MM-DDThh:mm:ss.mss (year-month-day T hour:minute:second:milisecond in directory
timeseries_data_fn is the filename of the harps-n timeseries RV data (e.g. harpn_sun_release_timeseries_2015-2018.csv)
"""
function make_manifest(data_path::String ; timeseries_data_fn::String = "", max_spectra_to_use::Int = 1000)
    df_filenames = EchelleInstruments.make_manifest(data_path,r"^20\d{2}-\d{2}-\d{2}T\d{2}:\d{2}:\d{2}.\d{3}$")

    df_files = DataFrame(read_metadata(joinpath(df_filenames.Filename[1],"s2d.fits")))
    keys = propertynames(df_files)
    allowmissing!(df_files, keys[map(k->k∉[:Filename, :bjd, :target, :airmass],keys)] )

    if length(df_filenames.Filename) >= 2
        map(fn->add_metadata_from_fits!(df_files,joinpath(fn,"s2d.fits")),df_filenames.Filename[2:end])
    end

    #load rv data file for harps-n solar data
    if timeseries_data_fn != ""
        if isfile(joinpath(data_path,timeseries_data_fn))
            RV_data = CSV.read(joinpath(data_path,timeseries_data_fn), DataFrame)
            # berv = RV_data[:, [:berv]]
            berv_bary_to_helio = RV_data[:, [:berv_bary_to_helio]] .* -1.0
            rv_diff_extinction = RV_data[:, [:rv_diff_extinction]] .* -1.0

            df_files = hcat(df_files, berv_bary_to_helio, rv_diff_extinction)
            append!(keys,[:ssb_rv, :diff_ext_rv])
            rename!(df_files,keys)
        else
            @info string("HARPSN.make_manifest did not find a file at",joinpath(data_path,timeseries_data_fn),".\n RV data cannot be loaded.")
        end
    end

    df_files

end

"""Create Dict containing filename and default metadata from file."""
function read_metadata(fn::Union{String,FITS})
    fields_to_save = metadata_symbols_default(HARPSN2D())
    fields_str_to_save = metadata_strings_default(HARPSN2D())
    dict = read_metadata_from_fits(fn,fields=fields_to_save,fields_str=fields_str_to_save)
    check_metadata_fields_expected_present(dict,fields_to_save)
    dict[:Filename] = fn
    if haskey(dict,:airmass) && typeof(dict[:airmass]) == String
       dict[:airmass] = parse(Float64,dict[:airmass])
    end
    return dict
end

function add_metadata_from_fits!(df::DataFrame, fn::String)
    println("# Reading metadata from ", fn)
    metadata_keep = read_metadata(fn)
    #reject file from df if it is a calibration (only keep spectra)
    #if metadata_keep[:target] != "Solar" return df end
    #calibration_target_substrings = ["Etalon","LFC","ThArS","LDLS","FlatBB"]
    #if any(map(ss->contains( metadata_keep[:target],ss),calibration_target_substrings))
    #   return df
    #end

    for k in propertynames(df)
        if typeof(metadata_keep[k]) == Missing continue end
        if typeof(metadata_keep[k]) == Nothing
            metadata_keep[k] = missing
            continue
        end
        if (eltype(df[!,k]) <: Union{Real,Union{<:Real,Missing}}) && !(typeof(metadata_keep[k]) <: Union{Real,Union{<:Real,Missing}})
            metadata_keep[k] = missing
        end
    end
    #println("metadata_keep = ", metadata_keep)
    if haskey(metadata_keep,:airmass) && typeof(metadata_keep[:airmass]) == String
       metadata_keep[:airmass] = parse(Float64,metadata_keep[:airmass])
    end

    push!(df, metadata_keep)
    return df
end

""" Read HARPS-N data from FITS file, and return in a Spectra2DBasic object."""
function read_data   end

function read_data(f::FITS, metadata::Dict{Symbol,Any} )
    @assert length(f) >= 2
    @assert all(map(i->typeof(f[i]) <: FITSIO.ImageHDU,2:length(f)))
    img_idx = Dict(map(i->first(read_key(f[i],"EXTNAME"))=>i,2:length(f)))
    λ, flux, var  = FITSIO.read(f[img_idx["WAVEDATA_AIR_BARY"]]), FITSIO.read(f[img_idx["SCIDATA"]]), FITSIO.read(f[img_idx["ERRDATA"]]).^2
    metadata[:normalization] = :raw
    spectrum = Spectra2DBasic(λ, flux, var, HARPSN2D(), metadata=metadata)
    apply_doppler_boost!(spectrum,metadata)
    return spectrum
end

function read_data(f::FITS, metadata::Dict{Symbol,Any}, orders_to_read::AR ) where  { AR<:AbstractRange }
    @assert length(f) >= 2
    @assert all(map(i->typeof(f[i]) <: FITSIO.ImageHDU,2:length(f)))
    img_idx = Dict(map(i->first(read_key(f[i],"EXTNAME"))=>i,2:length(f)))
    λ, flux, var  = FITSIO.read(f[img_idx["WAVEDATA_AIR_BARY"]]), FITSIO.read(f[img_idx["SCIDATA"]]), FITSIO.read(f[img_idx["ERRDATA"]]).^2
    metadata[:normalization] = :raw
    spectrum = Spectra2DBasic(λ, flux, var, HARPSN2D(), metadata=metadata)
    apply_doppler_boost!(spectrum,metadata)
    return spectrum
end

function read_data(fn::String, metadata::Dict{Symbol,Any} )
    f = FITS(fn)
    read_data(f, metadata)
end

function read_data(fn::String)
    f = FITS(fn)
    hdr = FITSIO.read_header(f[1])
    metadata = Dict(zip(map(k->Symbol(k),hdr.keys),hdr.values))
    read_data(f, metadata)
end

function read_data(dfr::DataFrameRow{DataFrame,DataFrames.Index})
    fn = dfr.Filename
    metadata = Dict(zip(keys(dfr),values(dfr)))
    read_data(fn,metadata)
end

function read_data(fn::String, orders_to_read::AR ) where  { AR<:AbstractRange }
    f = FITS(fn)
    metadata = read_metadata(f)
    read_data(f,metadata, orders_to_read)
end

function read_data(dfr::DataFrameRow{DataFrame,DataFrames.Index}, orders_to_read::AR ) where  { AR<:AbstractRange }
    fn = dfr.Filename
    metadata = Dict(zip(keys(dfr),values(dfr)))
    read_data(fn,metadata, orders_to_read)
end

function apply_barycentric_correction!(λ::AbstractArray{T,1}, z::Real) where { T<:Real }

end

""" Read normalization function for approximate continuum+blaze correction
TODO: Figure out mapping between these orders and what's stored in the FITS file before using. """
function read_normalization(fn::String = joinpath(pkgdir(EchelleInstruments),"data","neid","HARPSN_normalization_function.csv") )
    CSV.read(fn, DataFrame)
end

max_stars_repo_name: RvSpectML/EchelleInstruments

max_issues_repo_name: RvSpectML/EchelleInstruments

max_forks_repo_name: RvSpectML/EchelleInstruments


text:
using Documenter
using SPHKernels

makedocs(
    sitename = "SPHKernels",
    format = Documenter.HTML(),
    modules = [SPHKernels],
    pages = [
            "Table of Contents" => "index.md",
            "Install" => "install.md",
            "Kernels" => "kernels.md",
            "Extending" => "extending.md",
            "API reference" => "api.md"
            ]
)

# Documenter can also automatically deploy documentation to gh-pages.
# See "Hosting Documentation" and deploydocs() in the Documenter manual
# for more information.
deploydocs(
    repo = "github.com/LudwigBoess/SPHKernels.jl.git"
)


text:
# convergence stats
close("all")
transD_GP.getchi2forall(opt, fsize=8, alpha=0.5)
ax = gcf().axes;
r = vec(img.d[img.select] - img.f[img.select])
χ² = length(r)
figure(1)
ax[2].set_ylim(χ²/2 - 100, χ²/2 + 300)
ax[2].plot(xlim(), [χ²/2 , χ²/2], "--", color="gray")
figure(2)
ax = gcf().axes
ax[2].set_ylim(χ²/2 - 100, χ²/2 + 300)
ax[2].plot(xlim(), [χ²/2 , χ²/2], "--", color="gray")
## posterior marginals
rownum, colnum = 195, 85
transD_GP.plot_image_posterior(opt, img, burninfrac=0.5, rownum=rownum, colnum=colnum, nbins=100)
savefig("post_s.png", dpi=300)

max_stars_repo_name: GeoscienceAustralia/HiQGA

max_issues_repo_name: GeoscienceAustralia/HiQGA


text:
#!/usr/bin/env julia

function sam(genome_path, bam_paths...; alt_reads=3, alt_frac=0.1,
	min_mapq=10, region="", max_depth=1_000_000)
	lines = []
	genome_path = expanduser(genome_path)
	bam_paths = map(expanduser, bam_paths)

	samples = map(s -> replace(s, r".bam$", ""), bam_paths)
	S = length(samples)
	println("CHROM\tPOSITION\tREF\tALT\t", join(samples, '\t'))

	#spileup = expanduser("~/tools/pypette/compiled/spileup")
	spileup = expanduser("/home/gnanavel/Work/spileup")
	mpileup_options = ""
	if region != ""; mpileup_options = "-r$(region)"; end

	max_alleles = 100
	push!(lines, pipeline(`samtools mpileup -d $(max_depth) -A -x -R -sB -q0 $(mpileup_options) -f $(genome_path) $(bam_paths)`, `$(spileup) $(alt_reads) $(min_mapq)`))

	return lines
end

function format(lines)
	for line in lines
		cols = split(rstrip(line), '\t')
		assert(length(cols) >= 3)
		if cols[3] == "N"; continue; end

		total_reads = zeros(Int32, length(samples))
		allele_reads = zeros(Int32, length(samples), max_alleles)
		alleles = Array{ASCIIString}(0)

		for (s, pileup) in enumerate(cols[4:end])
			if pileup == ""; continue; end
			t = split(pileup, ' ')
			for a in 1:3:length(t)
				count = int(t[a+1])
				total_reads[s] += count
				if t[a] == "."; continue; end
				idx = findfirst(allele -> allele == t[a], alleles)
				if idx == 0 && length(alleles) >= max_alleles; continue; end
				if idx == 0; push!(alleles, t[a]); idx = length(alleles); end
				allele_reads[s, idx] = count
			end
		end

		for a in 1:length(alleles)
			frac = allele_reads[:, a] ./ total_reads
			alt_gt = (allele_reads[:, a] .>= alt_reads) & (frac .>= alt_frac)
			if !any(alt_gt); continue; end

			# Reformat indels in VCF4 format
			ref = cols[3]; alt = alleles[a]
			if length(alt) >= 2
				if alt[2] == '+'    # Insertion
					alt = string(alt[1] == '.' ? ref : alt[1], alt[3:end])
				elseif alt[2] == '-'  # Deletion
					ref = string(ref, alt[3:end])
					alt = alt[1] == '.' ? ref[1] : alt[1]
				end
			end

			print("$(cols[1])\t$(cols[2])\t$(uppercase(ref))\t$(uppercase(alt))")
			for s in 1:length(samples)
				print("\t$(allele_reads[s, a]):$(total_reads[s])")
				if alt_gt[s]; print(":*"); end
			end
			println()
		end
	end
end

function call(genome_path, bam_paths...; alt_reads=3, alt_frac=0.1,
	min_mapq=10, region="", max_depth=1_000_000)

	genome_path = expanduser(genome_path)
	bam_paths = map(expanduser, bam_paths)

	samples = map(s -> replace(s, r".bam$", ""), bam_paths)
	S = length(samples)
	println("CHROM\tPOSITION\tREF\tALT\t", join(samples, '\t'))

	#spileup = expanduser("~/tools/pypette/compiled/spileup")
	spileup = expanduser("/home/gnanavel/Work/spileup")
	mpileup_options = ""
	if region != ""; mpileup_options = "-r$(region)"; end

	max_alleles = 100
	for line in eachline(pipeline(`samtools mpileup -d $(max_depth) -A -x -R -sB -q0 $(mpileup_options) -f $(genome_path) $(bam_paths)`, `$(spileup) $(alt_reads) $(min_mapq)`))
		cols = split(rstrip(line), '\t')
		assert(length(cols) >= 3)
		if cols[3] == "N"; continue; end

		total_reads = zeros(Int32, length(samples))
		allele_reads = zeros(Int32, length(samples), max_alleles)
		alleles = Array{ASCIIString}(0)

		for (s, pileup) in enumerate(cols[4:end])
			if pileup == ""; continue; end
			t = split(pileup, ' ')
			for a in 1:3:length(t)
				count = int(t[a+1])
				total_reads[s] += count
				if t[a] == "."; continue; end
				idx = findfirst(allele -> allele == t[a], alleles)
				if idx == 0 && length(alleles) >= max_alleles; continue; end
				if idx == 0; push!(alleles, t[a]); idx = length(alleles); end
				allele_reads[s, idx] = count
			end
		end

		for a in 1:length(alleles)
			frac = allele_reads[:, a] ./ total_reads
			alt_gt = (allele_reads[:, a] .>= alt_reads) & (frac .>= alt_frac)
			if !any(alt_gt); continue; end

			# Reformat indels in VCF4 format
			ref = cols[3]; alt = alleles[a]
			if length(alt) >= 2
				if alt[2] == '+'    # Insertion
					alt = string(alt[1] == '.' ? ref : alt[1], alt[3:end])
				elseif alt[2] == '-'  # Deletion
					ref = string(ref, alt[3:end])
					alt = alt[1] == '.' ? ref[1] : alt[1]
				end
			end

			print("$(cols[1])\t$(cols[2])\t$(uppercase(ref))\t$(uppercase(alt))")
			for s in 1:length(samples)
				print("\t$(allele_reads[s, a]):$(total_reads[s])")
				if alt_gt[s]; print(":*"); end
			end
			println()
		end
	end
end


text:
function axis(gr::Graph; dir::String="xyz", stl::String="", opt::String="")
    cmgl.axis(gr.ptr, dir, stl, opt)
end

#=function axis(gr::Graph, dir::Symbol=:xyz; stl::String="", opt::String="")=#
    #=cmgl.axis(gr.ptr, string(dir), stl, opt)=#
#=end=#

function colorbar(gr::Graph; stl::String="")
    cmgl.colorbar(gr.ptr, stl)
end

function colorbar(gr::Graph, v::Array{cmgl.Float, 1}; stl::String="")
    cmgl.colorbar_val(gr.ptr, Data(v).ptr, stl)
end

function colorbar(gr::Graph, x::Real, y::Real; w::Real=1., h::Real=1., stl::String="")
    cmgl.colorbar_ext(gr.ptr, stl, x, y, w, h)
end

function colorbar(gr::Graph, v::Array{cmgl.Float, 1}, x::Real, y::Real; w::Real=1., h::Real=1., stl::String="")
    cmgl.colorbar_val_ext(gr.ptr, Data(v).ptr, stl, x, y, w, h)
end

function grid(gr::Graph, dir::String="xyz"; stl::String="h-", opt::String="")
    cmgl.axis_grid(gr.ptr, dir, stl, opt)
end

function subgrid(gr::Graph, dir::String="xyz"; stl::String="h=", opt::String="")
    cmgl.axis_grid(gr.ptr, "!"*dir, stl, opt)
end

function box(gr::Graph; stl::String="", ticks::Bool=true)
    cmgl.box_str(gr.ptr, stl, ticks)
end

#=function box(gr::Graph) <- causes overwrite-warning=#
    #=cmgl.box(gr.ptr)=#
#=end=#

function label(gr::Graph, dir::Char, text::String; pos::Real=0., opt::String="")
    cmgl.label(gr.ptr, dir, text, pos, opt)
end

#=function label(gr::Graph, dir::Char, text::UTF8String;  pos::Real=0., opt::String="")=#
    #=cmgl.labelw(gr.ptr, dir, text, pos, opt)=#
#=end=#

function xlabel(gr::Graph, text::String; pos::Real=0., opt::String="")
    cmgl.label(gr.ptr, 'x', text, pos, opt)
end

function ylabel(gr::Graph, text::String; pos::Real=0., opt::String="")
    cmgl.label(gr.ptr, 'y', text, pos, opt)
end

function zlabel(gr::Graph, text::String; pos::Real=0., opt::String="")
    cmgl.label(gr.ptr, 'z', text, pos, opt)
end

#=function xlabel(gr::Graph, text::UTF8String; pos::Real=0., opt::String="")=#
    #=cmgl.labelw(gr.ptr, 'x', text, pos, opt)=#
#=end=#

#=function ylabel(gr::Graph, text::UTF8String; pos::Real=0., opt::String="")=#
    #=cmgl.labelw(gr.ptr, 'y', text, pos, opt)=#
#=end=#

#=function zlabel(gr::Graph, text::UTF8String; pos::Real=0., opt::String="")=#
    #=cmgl.labelw(gr.ptr, 'z', text, pos, opt)=#
#=end=#

function legend(gr::Graph; pos::Int=3, fnt::String="", opt::String="")
    cmgl.legend(gr.ptr, pos, fnt, opt)
end

function legend(gr::Graph, x::Real, y::Real; fnt::String="", opt::String="")
    cmgl.legend_pos(gr.ptr, x, y, fnt, opt)
end

function addlegend(gr::Graph, text::String; stl::String="")
    cmgl.add_legend(gr.ptr, text, stl)
end

#=function addlegend(gr::Graph, text::UTF8String; stl::String="")=#
    #=cmgl.add_legendw(gr.ptr, text, stl)=#
#=end=#

function clearlegend(gr::Graph)
    cmgl.clear_legend(gr.ptr)
end

function legendamarks(gr::Graph, num::Int)
    cmgl.set_legend_marks(gr.ptr, num)
end


text:
#####
##### Gradient, forward finite-differences, periodic boundaries
#####

"""
    gradfp(
        u::AbstractArray{<:AbstractFloat, 3},
        h::NTuple{3, Real}
    ) -> NTuple{3, typeof(similar(u))}

First order forward difference gradient with periodic boundaries.

### Arguments
- `u::AbstractArray{<:AbstractFloat, 3}`: input array
- `h::NTuple{3, Real}`: grid spacing

### Returns
- `dx::typeof(similar(u))`: x-component of gradient of `u`
- `dy::typeof(similar(u))`: y-component of gradient of `u`
- `dz::typeof(similar(u))`: z-component of gradient of `u`
"""
gradfp(u::AbstractArray{<:AbstractFloat, 3}, h::NTuple{3, Real}) =
    gradfp!(similar(u), similar(u), similar(u), u, h)

"""
    gradfp!(
        dx::AbstractArray{<:AbstractFloat, 3},
        dy::AbstractArray{<:AbstractFloat, 3},
        dz::AbstractArray{<:AbstractFloat, 3},
        u::AbstractArray{<:AbstractFloat, 3},
        h::NTuple{3, Real},
    ) -> (dx, dy, dz)

First order forward difference gradient with periodic boundaries.

### Arguments
- `dx::AbstractArray{<:AbstractFloat, 3}`: x-component of gradient of `u`
- `dy::AbstractArray{<:AbstractFloat, 3}`: y-component of gradient of `u`
- `dz::AbstractArray{<:AbstractFloat, 3}`: z-component of gradient of `u`
- `u::AbstractArray{<:AbstractFloat, 3}`: input array
- `h::NTuple{3, Real}`: grid spacing

### Returns
- `dx`: x-component of gradient of `u`
- `dy`: y-component of gradient of `u`
- `dz`: z-component of gradient of `u`
"""
function gradfp!(dx, dy, dz, u, h)
    size(dx) == size(u) || throw(DimensionMismatch())
    size(dy) == size(u) || throw(DimensionMismatch())
    size(dz) == size(u) || throw(DimensionMismatch())
    return _gradfp!(dx, dy, dz, u, h)
end

function _gradfp!(
    dx::AbstractArray{Tdx, 3},
    dy::AbstractArray{Tdy, 3},
    dz::AbstractArray{Tdz, 3},
    u::AbstractArray{<:AbstractFloat, 3},
    h::NTuple{3, Real},
) where {Tdx<:AbstractFloat, Tdy<:AbstractFloat, Tdz<:AbstractFloat}
    nx, ny, nz = size(u)
    hx = convert(Tdx, inv(h[1]))
    hy = convert(Tdy, inv(h[2]))
    hz = convert(Tdz, inv(h[3]))

    @inbounds @batch for k in 1:nz-1
        for j in 1:ny-1
            for i in 1:nx-1
                dx[i,j,k] = hx*(u[i+1,j,k] - u[i,j,k])
                dy[i,j,k] = hy*(u[i,j+1,k] - u[i,j,k])
                dz[i,j,k] = hz*(u[i,j,k+1] - u[i,j,k])
            end

            dx[nx,j,k] = hx*(u[1,j,k] - u[nx,j,k])
            dy[nx,j,k] = hy*(u[nx,j+1,k] - u[nx,j,k])
            dz[nx,j,k] = hz*(u[nx,j,k+1] - u[nx,j,k])
        end

        for i in 1:nx-1
            dy[i,ny,k] = hy*(u[i,1,k] - u[i,ny,k])
            dx[i,ny,k] = hx*(u[i+1,ny,k] - u[i,ny,k])
            dz[i,ny,k] = hz*(u[i,ny,k+1] - u[i,ny,k])
        end

        dy[nx,ny,k] = hy*(u[nx,1,k] - u[nx,ny,k])
        dx[nx,ny,k] = hx*(u[1,ny,k] - u[nx,ny,k])
        dz[nx,ny,k] = hz*(u[nx,ny,k+1] - u[nx,ny,k])
    end

    @inbounds @batch for j in 1:ny-1
        for i in 1:nx-1
            dz[i,j,nz] = hz*(u[i,j,1] - u[i,j,nz])
            dx[i,j,nz] = hx*(u[i+1,j,nz] - u[i,j,nz])
            dy[i,j,nz] = hy*(u[i,j+1,nz] - u[i,j,nz])
        end

        dz[nx,j,nz] = hz*(u[nx,j,1] - u[nx,j,nz])
        dx[nx,j,nz] = hx*(u[1,j,nz] - u[nx,j,nz])
        dy[nx,j,nz] = hy*(u[nx,j+1,nz] - u[nx,j,nz])
    end

    @inbounds @batch for i in 1:nx-1
        dz[i,ny,nz] = hz*(u[i,ny,1] - u[i,ny,nz])
        dy[i,ny,nz] = hy*(u[i,1,nz] - u[i,ny,nz])
        dx[i,ny,nz] = hx*(u[i+1,ny,nz] - u[i,ny,nz])
    end

    @inbounds begin
        dz[nx,ny,nz] = hz*(u[nx,ny,1] - u[nx,ny,nz])
        dy[nx,ny,nz] = hy*(u[nx,1,nz] - u[nx,ny,nz])
        dx[nx,ny,nz] = hx*(u[1,ny,nz] - u[nx,ny,nz])
    end

    return dx, dy, dz
end

"""
    gradfp_adj(
        dx::AbstractArray{<:AbstractFloat, 3},
        dy::AbstractArray{<:AbstractFloat, 3},
        dz::AbstractArray{<:AbstractFloat, 3},
        h::NTuple{3, Real}
    ) -> typeof(similar(dx, promote_eltype(dx, dy, dz)))

Adjoint of first order forward difference gradient with periodic boundaries.

### Arguments
- `dx::AbstractArray{<:AbstractFloat, 3}`: x-component
- `dy::AbstractArray{<:AbstractFloat, 3}`: y-component
- `dz::AbstractArray{<:AbstractFloat, 3}`: z-component
- `h::NTuple{3, Real}`: grid spacing

### Returns
- `u::typeof(similar(dx, promote_eltype(dx, dy ,dz)))`: divergence of [dx, dy, dz]
"""
function gradfp_adj(
    dx::AbstractArray{Tdx, 3},
    dy::AbstractArray{Tdy, 3},
    dz::AbstractArray{Tdz, 3},
    h::NTuple{3, Real}
) where {Tdx, Tdy, Tdz}
    Td2u = foldl(promote_type, (Tdy, Tdz), init=Tdx)
    d2u = similar(dx, Td2u)
    gradfp_adj!(d2u, dx, dy, dz, h)
end

"""
    gradfp_adj!(
        u::AbstractArray{<:AbstractFloat, 3}
        dx::AbstractArray{<:AbstractFloat, 3},
        dy::AbstractArray{<:AbstractFloat, 3},
        dz::AbstractArray{<:AbstractFloat, 3},
        h::NTuple{3, Real}
    ) -> u

Adjoint of first order forward difference gradient with periodic boundaries.

### Arguments
- `u::AbstractArray{<:AbstractFloat, 3}`: divergence of [dx, dy, dz]
- `dx::AbstractArray{<:AbstractFloat, 3}`: x-component
- `dy::AbstractArray{<:AbstractFloat, 3}`: y-component
- `dz::AbstractArray{<:AbstractFloat, 3}`: z-component
- `h::NTuple{3, Real}`: grid spacing

### Returns
- `u`: divergence of [dx, dy, dz]
"""
function gradfp_adj!(d2u, dx, dy, dz, h)
    size(d2u) == size(dx) || throw(DimensionMismatch())
    size(d2u) == size(dy) || throw(DimensionMismatch())
    size(d2u) == size(dz) || throw(DimensionMismatch())
    return _gradfp_adj!(d2u, dx, dy, dz, h)
end

function _gradfp_adj!(
    d2u::AbstractArray{<:AbstractFloat, 3},
    dx::AbstractArray{Tdx, 3},
    dy::AbstractArray{Tdy, 3},
    dz::AbstractArray{Tdz, 3},
    h::NTuple{3, Real},
) where {Tdx, Tdy, Tdz}
    nx, ny, nz = size(d2u)
    hx = convert(Tdx, -inv(h[1]))
    hy = convert(Tdy, -inv(h[2]))
    hz = convert(Tdz, -inv(h[3]))

    @inbounds begin
        d2u[1,1,1] =
            hx*(dx[1,1,1] - dx[nx,1,1]) +
            hy*(dy[1,1,1] - dy[1,ny,1]) +
            hz*(dz[1,1,1] - dz[1,1,nz])
    end

    @inbounds @batch for i in 2:nx
        ux = dx[i,1,1] - dx[i-1,1,1]
        uy = dy[i,1,1] - dy[i,ny,1]
        uz = dz[i,1,1] - dz[i,1,nz]

        u = hx*ux
        u = muladd(hy, uy, u)
        u = muladd(hz, uz, u)
        d2u[i,1,1] = u
    end

    @inbounds @batch for j in 2:ny
        uy = dy[1,j,1] - dy[1,j-1,1]
        ux = dx[1,j,1] - dx[nx,j,1]
        uz = dz[1,j,1] - dz[1,j,nz]

        u = hy*uy
        u = muladd(hx, ux, u)
        u = muladd(hz, uz, u)
        d2u[1,j,1] = u

        for i in 2:nx
            uy = dy[i,j,1] - dy[i,j-1,1]
            ux = dx[i,j,1] - dx[i-1,j,1]
            uz = dz[i,j,1] - dz[i,j,nz]

            u = hy*uy
            u = muladd(hx, ux, u)
            u = muladd(hz, uz, u)
            d2u[i,j,1] = u
        end
    end

    @inbounds @batch for k in 2:nz
        uz = dz[1,1,k] - dz[1,1,k-1]
        ux = dx[1,1,k] - dx[nx,1,k]
        uy = dy[1,1,k] - dy[1,ny,k]

        u = hz*uz
        u = muladd(hx, ux, u)
        u = muladd(hy, uy, u)
        d2u[1,1,k] = u

        for i in 2:nx
            uz = dz[i,1,k] - dz[i,1,k-1]
            uy = dy[i,1,k] - dy[i,ny,k]
            ux = dx[i,1,k] - dx[i-1,1,k]

            u = hz*uz
            u = muladd(hy, uy, u)
            u = muladd(hx, ux, u)
            d2u[i,1,k] = u
        end

        for j in 2:ny
            uz = dz[1,j,k] - dz[1,j,k-1]
            uy = dy[1,j,k] - dy[1,j-1,k]
            ux = dx[1,j,k] - dx[nx,j,k]

            u = hz*uz
            u = muladd(hy, uy, u)
            u = muladd(hx, ux, u)
            d2u[1,j,k] = u

            for i in 2:nx
                uz = dz[i,j,k] - dz[i,j,k-1]
                uy = dy[i,j,k] - dy[i,j-1,k]
                ux = dx[i,j,k] - dx[i-1,j,k]

                u = hz*uz
                u = muladd(hy, uy, u)
                u = muladd(hx, ux, u)
                d2u[i,j,k] = u
            end
        end
    end

    return d2u
end

#####
##### Laplacian, central finite-differences
#####

"""
    lap(
        u::AbstractArray{<:AbstractFloat, 3},
        h::NTuple{3, Real}
    ) -> typeof(similar(u))

Second order central difference Laplacian.

### Arguments
- `u::AbstractArray{<:AbstractFloat, 3}`: input array
- `h::NTuple{3, Real}`: grid spacing

### Returns
- `d2u::typeof(similar(u))`: discrete Laplacian of `u`
"""
lap(u::AbstractArray{<:AbstractFloat}, h::NTuple{3, Real}) =
    lap!(tzero(u), u, h)

"""
    lap!(
        d2u::AbstractArray{<:AbstractFloat, 3},
        u::AbstractArray{<:AbstractFloat, 3},
        h::NTuple{3, Real}
    ) -> d2u

Second order central difference Laplacian.

### Arguments
- `d2u::AbstractArray{<:AbstractFloat, 3}`: discrete Laplacian of `u`
- `u::AbstractArray{<:AbstractFloat, 3}`: input array
- `h::NTuple{3, Real}`: grid spacing

### Returns
- `d2u`: discrete Laplacian of `u`
"""
function lap!(d2u, u, h)
    size(d2u) == size(u) || throw(DimensionMismatch())
    return _lap!(d2u, u, h)
end

function _lap!(
    d2u::AbstractArray{<:AbstractFloat, 3},
    u::AbstractArray{T, 3},
    h::NTuple{3, Real},
) where {T<:AbstractFloat}
    idx2 = convert(T, inv(h[1]*h[1]))
    idy2 = convert(T, inv(h[2]*h[2]))
    idz2 = convert(T, inv(h[3]*h[3]))
    D = -2*(idx2 + idy2 + idz2)

    ax = map(a -> first(a)+1:last(a)-1, axes(u))

    tsz = padded_tilesize(T, (2, 2, 2), 1)
    R = vec(collect(TileIterator(ax, tsz)))

    @inbounds @batch for (I, J, K) in R
        for k in K
            for j in J
                for i in I
                    z1 = u[i,j,k-1]
                    y1 = u[i,j-1,k]
                    x1 = u[i-1,j,k]
                    u0 = u[i,j,k]
                    x2 = u[i+1,j,k]
                    y2 = u[i,j+1,k]
                    z2 = u[i,j,k+1]

                    Δz = z1 + z2
                    Δy = y1 + y2
                    Δx = x1 + x2

                    Δ = D * u0
                    Δ = muladd(idz2, Δz, Δ)
                    Δ = muladd(idy2, Δy, Δ)
                    Δ = muladd(idx2, Δx, Δ)

                    d2u[i,j,k] = Δ
                end
            end
        end
    end

    return d2u
end

function _lap!(
    d2u::AbstractArray{<:AbstractFloat, 4},
    u::AbstractArray{T, 4},
    h::NTuple{3, Real},
) where {T<:AbstractFloat}
    idx2 = convert(T, inv(h[1]*h[1]))
    idy2 = convert(T, inv(h[2]*h[2]))
    idz2 = convert(T, inv(h[3]*h[3]))
    D = -2*(idx2 + idy2 + idz2)

    ax = map(a -> first(a)+1:last(a)-1, axes(u)[1:3])

    tsz = padded_tilesize(T, (2, 2, 2), 1)
    R = vec(collect(TileIterator(ax, tsz)))

    @inbounds for t in axes(u, 4)
        ut = @view(u[:,:,:,t])
        d2ut = @view(d2u[:,:,:,t])
        @batch for (I, J, K) in R
            for k in K
                for j in J
                    for i in I
                        z1 = ut[i,j,k-1]
                        y1 = ut[i,j-1,k]
                        x1 = ut[i-1,j,k]
                        u0 = ut[i,j,k]
                        x2 = ut[i+1,j,k]
                        y2 = ut[i,j+1,k]
                        z2 = ut[i,j,k+1]

                        Δz = z1 + z2
                        Δy = y1 + y2
                        Δx = x1 + x2

                        Δ = D * u0
                        Δ = muladd(idz2, Δz, Δ)
                        Δ = muladd(idy2, Δy, Δ)
                        Δ = muladd(idx2, Δx, Δ)

                        d2ut[i,j,k] = Δ
                    end
                end
            end
        end
    end

    return d2u
end


text:
export linsolve!,linsolve ,  LU , GMRES , BICGSTAB
export eigsolve, EigArpack

include("EigenSolver.jl")
include("LinearSolver.jl")


text:
function HansenTesseralLinear(arg0::jint, arg1::jint, arg2::jint, arg3::jint, arg4::jint)
    return HansenTesseralLinear((jint, jint, jint, jint, jint), arg0, arg1, arg2, arg3, arg4)
end

function compute_init_values(obj::HansenTesseralLinear, arg0::jdouble, arg1::jdouble, arg2::jdouble)
    return jcall(obj, "computeInitValues", void, (jdouble, jdouble, jdouble), arg0, arg1, arg2)
end

function get_derivative(obj::HansenTesseralLinear, arg0::jint, arg1::jdouble)
    return jcall(obj, "getDerivative", jdouble, (jint, jdouble), arg0, arg1)
end

function get_value(obj::HansenTesseralLinear, arg0::jint, arg1::jdouble)
    return jcall(obj, "getValue", jdouble, (jint, jdouble), arg0, arg1)
end


text:
#####################
# Generic evaluation
#####################

"Compute the evaluation matrix of the given dict on the given set of points."
function evaluation_matrix(Φ::Dictionary, pts, T = codomaintype(Φ))
    A = Array{T}(undef, length(pts), length(Φ))
    evaluation_matrix!(A, Φ, pts)
end

function evaluation_matrix!(A::AbstractMatrix, Φ::Dictionary, pts)
    @assert size(A,1) == length(pts)
    @assert size(A,2) == length(Φ)

    for (j,ϕ) in enumerate(Φ), (i,x) in enumerate(pts)
        A[i,j] = ϕ(x)
    end
    A
end

function dense_evaluation(::Type{T}, Φ::Dictionary, gb::GridBasis; options...) where {T}
    A = evaluation_matrix(Φ, grid(gb), T)
    ArrayOperator(A, Φ, gb)
end

evaluation(Φ::Dictionary, grid::AbstractGrid; options...) =
    evaluation(coefficienttype(Φ), Φ, grid; options...)

evaluation(::Type{T}, Φ::Dictionary, grid::AbstractGrid; options...) where {T} =
    evaluation(T, Φ, GridBasis{T}(grid); options...)

evaluation(::Type{T}, Φ::Dictionary, grid::SubGrid; options...) where {T} =
     restriction(T, supergrid(grid), grid; options...) * evaluation(T, Φ, supergrid(grid); options...)

function evaluation(::Type{T}, Φ::Dictionary, gb::GridBasis, grid; warnslow = true, options...) where {T}
    # warnslow && @debug "No fast evaluation available in $grid for dictionary $Φ, using dense evaluation matrix instead."
    dense_evaluation(T, Φ, gb; options...)
end

function resize_and_transform(::Type{T}, Φ::Dictionary, gb::GridBasis, grid; options...) where {T}
    if size(Φ) == size(grid)
        transform_to_grid(T, Φ, gb, grid; options...)
    elseif length(grid) > length(Φ)
        dlarge = resize(Φ, size(grid))
        transform_to_grid(T, dlarge, gb, grid; options...) * extension(T, Φ, dlarge; options...)
    else
        @debug "Resize and transform: dictionary evaluated in small grid"
        dense_evaluation(T, Φ, gb; options...)
    end
end


text:
using DashBootstrapComponents

nav_contents = [
    dbc_navitem(dbc_navlink("Active", href = "#", active = true)),
    dbc_navitem(dbc_navlink("A much longer link label", href = "#")),
    dbc_navitem(dbc_navlink("Link", href = "#")),
]

nav1 = dbc_nav(nav_contents, pills = true, fill = true);

nav2 = dbc_nav(nav_contents, pills = true, justified = true);

navs = html_div([nav1, html_hr(), nav2]);

max_issues_repo_name: sthagen/dash-bootstrap-components

max_forks_repo_name: sthagen/dash-bootstrap-components


text:
#   Debug.Graft:
# ================
# Debug instrumentation of code, and transformation of ASTs to act as if they
# were evaluated inside such code (grafting)

module Graft
using Debug.AST, Debug.Meta, Debug.Analysis, Debug.Runtime
using Compat
import Debug.Meta, Debug.AST.is_emittable
export instrument, graft, @localscope, @notrap

# ---- @localscope: returns the Scope instance for the local scope ------------

type GetLocalScope; end
macro localscope()
    code_analyzed_only(Node(GetLocalScope()),
        "@localscope can only be used within @debug or @debug_analyze")
end

# ---- @notrap: prevents generation of traps in the wrapped code --------------

type NoTrap; end
is_emittable(::Node{NoTrap}) = false
macro notrap(ex)
    # Avoid creating any line number nodes so that the resulting block has
    # exactly two arguments, with Node(NoTrap()) being the first,
    # as expected below
    esc(Expr(:block, Node(NoTrap()), ex))
end

# ---- instrument -------------------------------------------------------------
# Add Scope creation and debug traps to (analyzed) code

type Context
    trap_pred::Function
    trap_ex
    env::Env
    scope_ex
end
Context(c::Context,e::Env) = Context(c.trap_pred,c.trap_ex,e,nothing)

function get_scope_ex(c::Context)
    # Generate a gensym on demand; if c.scope_ex is never used
    # it remains nothing
    c.scope_ex === nothing ? c.scope_ex = gensym("scope") : c.scope_ex
end

function instrument(trap_pred::Function, trap_ex, ex)
    @gensym scope
    ex = instrument(Context(trap_pred,trap_ex,NoEnv(),scope), analyze(ex,true))
    quote
        $scope = $(quot(ModuleScope))(eval)
        $ex
    end
end

function code_getset(sym::Symbol)
    val = gensym(string(sym))
    :( ()->$sym, $val->($sym=$val) )
end
function code_scope(scopesym::Symbol, parent, env::Env, syms)
    pairs = [Expr(:(=>), quot(sym), code_getset(sym)) for sym in syms]
    :(local $scopesym = $(quot(LocalScope))(
        $parent,
        $(Expr(:typed_dict,
               :($(quot(Symbol))=>$(quot(@compat(Tuple{Function,Function})))), pairs...)),
        $(quot(env))
    ))
end

code_trap(c::Context, node) = Expr(:call, c.trap_ex, quot(node), get_scope_ex(c))
code_trap_if(c::Context,node) = c.trap_pred(node) ? code_trap(c,node) : nothing

function instrument(c::Context, node::Node)
    if isa(node.state, Rhs) && !is_in_type(node) && c.trap_pred(node)
        Expr(:block, code_trap(c, node), instrument_node(c, node))
    else
        instrument_node(c, node)
    end
end

if VERSION >= v"0.4.0-dev"
    code_enterleave(::Void, ex, ::Void) = ex
    code_enterleave(enter,  ex, ::Void) = quote; $enter; $ex; end
    code_enterleave(::Void, ex, leave) = :(try $ex; finally $leave; end)
else
    code_enterleave(::Nothing, ex, ::Nothing) = ex
    code_enterleave(enter,     ex, ::Nothing) = quote; $enter; $ex; end
    code_enterleave(::Nothing, ex, leave) = :(try $ex; finally $leave; end)
end
code_enterleave(enter,     ex, leave) = :(try $enter; $ex; finally $leave; end)

function instrument_node(c::Context, node::Node)
    if !is_emittable(node); return quot(nothing); end
    ex = instrument_args(c, node)
    if is_scope_node(node)
        enter, leave = code_trap_if(c,Enter(node)), code_trap_if(c,Leave(node))
        if is_function(node)
            @assert is_function(ex)
            Expr(ex.head, ex.args[1], code_enterleave(enter,ex.args[2],leave))
        else
            code_enterleave(enter, ex, leave)
        end
    else
        ex
    end
end

never_trap(ex) = false

instrument_args(c::Context, node::Node) = exof(node)
instrument_args(c::Context, ::Node{GetLocalScope}) = get_scope_ex(c)
function instrument_args(c::Context, node::ExNode)
    args = Any[]
    orig_c = c

    # Detect @notrap
    if isblocknode(node) && nargsof(node) == 2 && isa(argof(node, 1), Node{NoTrap})
        old_trap, c.trap_pred = c.trap_pred, never_trap
        result = instrument(c, argof(node, 2))
        c.trap_pred = old_trap
        return result
    end

    if isblocknode(node) && !is(envof(node), c.env)
        # node introduces a scope, see if we need to reify it
        c = Context(c, envof(node))
        node.introduces_scope = true

        for arg in argsof(node); push!(args, instrument(c, arg)); end

        if c.scope_ex !== nothing
            # create new Scope
            syms, e = Set{Symbol}(), envof(node)
            while !is(e, orig_c.env);  union!(syms, e.defined); e = e.parent  end

            unshift!(args, code_scope(c.scope_ex, get_scope_ex(orig_c), envof(node), syms))
        end
    else
        for arg in argsof(node); push!(args, instrument(c, arg)); end
    end

    Expr(headof(node), args...)
end

# ---- graft ------------------------------------------------------------------
# Rewrite an (analyzed) AST to work as if it were inside
# the given scope, when evaluated in global scope.
# Replaces reads and writes to variables from that scope
# with getter/setter calls.

graft(env::Env, scope::Scope, ex) = rawgraft(scope, analyze(env, ex, false))
graft(scope::Scope, ex) =           graft(child(NoEnv()), scope, ex)

rawgraft(s::LocalScope, ex)         = ex
rawgraft(s::LocalScope, node::Node) = exof(node)
function rawgraft(s::LocalScope, ex::SymNode)
    sym = exof(ex)
    (haskey(s,sym) && !haskey(envof(ex),sym)) ? Expr(:call,quot(getter(s,sym))) : sym
end
function rawgraft(s::LocalScope, ex::Ex)
    head, args = headof(ex), argsof(ex)
    if head == :(=)
        lhs, rhs = args
        if isa(lhs, SymNode)             # assignment to symbol
            rhs = rawgraft(s, rhs)
            sym = exof(lhs)
            if haskey(envof(lhs), sym) || !(sym in s.env.assigned); return :($sym = $rhs)
            elseif haskey(s, sym);   return Expr(:call, quot(setter(s,sym)), rhs)
            else; error("No setter in scope found for $(sym)!")
            end
        elseif is_expr(lhs, :tuple)  # assignment to tuple
            tup = Node(Plain(gensym("tuple"))) # don't recurse into tup
            return rawgraft(s, Expr(:block,
                 :($tup  = $rhs    ),
                [:($dest = $tup[$k]) for (k,dest)=enumerate(argsof(lhs))]...))
        elseif is_expr(lhs, [:ref, :.]) || isa(lhs, PLeaf)# need no lhs rewrite
        else error("graft: not implemented: $ex")
        end
    elseif haskey(Analysis.updating_ops, head) && isa(args[1], SymNode)
        # x+=y ==> x=x+y etc.
        op = Analysis.updating_ops[head]
        return rawgraft(s, :( $(args[1]) = ($op)($(args[1]), $(args[2])) ))
    end
    Expr(head, [rawgraft(s,arg) for arg in args]...)
end

end # module


text:
@testset "xlogx & xlogy" begin
    @test iszero(xlogx(0))
    @test xlogx(2) ≈ 2.0 * log(2.0)
    @test_throws DomainError xlogx(-1)
    @test isnan(xlogx(NaN))

    @test iszero(xlogy(0, 1))
    @test xlogy(2, 3) ≈ 2.0 * log(3.0)
    @test_throws DomainError xlogy(1, -1)
    @test isnan(xlogy(NaN, 2))
    @test isnan(xlogy(2, NaN))
    @test isnan(xlogy(0, NaN))

    # Since we allow complex/negative values, test for them. See comments in:
    # https://github.com/JuliaStats/StatsFuns.jl/pull/95

    @test xlogx(1 + im) == (1 + im) * log(1 + im)
    @test isnan(xlogx(NaN + im))
    @test isnan(xlogx(1 + NaN * im))

    @test xlogy(-2, 3) == -xlogy(2, 3)
    @test xlogy(1 + im, 3) == (1 + im) * log(3)
    @test xlogy(1 + im, 2 + im) == (1 + im) * log(2 + im)
    @test isnan(xlogy(1 + NaN * im, -1 + im))
    @test isnan(xlogy(0, -1 + NaN * im))
    @test isnan(xlogy(Inf + im * NaN, 1))
    @test isnan(xlogy(0 + im * 0, NaN))
    @test iszero(xlogy(0 + im * 0, 0 + im * Inf))
end

@testset "logistic & logit" begin
    @test logistic(2) ≈ 1.0 / (1.0 + exp(-2.0))
    @test logistic(-750.0) === 0.0
    @test logistic(-740.0) > 0.0
    @test logistic(+36.0) < 1.0
    @test logistic(+750.0) === 1.0
    @test iszero(logit(0.5))
    @test logit(logistic(2)) ≈ 2.0
end

@testset "log1psq" begin
    @test iszero(log1psq(0.0))
    @test log1psq(1.0) ≈ log1p(1.0)
    @test log1psq(2.0) ≈ log1p(4.0)
end

# log1pexp, log1mexp, log2mexp & logexpm1

@testset "log1pexp" begin
    @test log1pexp(2.0)    ≈ log(1.0 + exp(2.0))
    @test log1pexp(-2.0)   ≈ log(1.0 + exp(-2.0))
    @test log1pexp(10000)  ≈ 10000.0
    @test log1pexp(-10000) ≈ 0.0

    @test log1pexp(2f0)      ≈ log(1f0 + exp(2f0))
    @test log1pexp(-2f0)     ≈ log(1f0 + exp(-2f0))
    @test log1pexp(10000f0)  ≈ 10000f0
    @test log1pexp(-10000f0) ≈ 0f0
end

@testset "log1mexp" begin
    @test log1mexp(-1.0)  ≈ log1p(- exp(-1.0))
    @test log1mexp(-10.0) ≈ log1p(- exp(-10.0))
end

@testset "log2mexp" begin
    @test log2mexp(0.0)  ≈ 0.0
    @test log2mexp(-1.0) ≈ log(2.0 - exp(-1.0))
end

@testset "logexpm1" begin
    @test logexpm1(2.0)            ≈  log(exp(2.0) - 1.0)
    @test logexpm1(log1pexp(2.0))  ≈  2.0
    @test logexpm1(log1pexp(-2.0)) ≈ -2.0

    @test logexpm1(2f0)            ≈  log(exp(2f0) - 1f0)
    @test logexpm1(log1pexp(2f0))  ≈  2f0
    @test logexpm1(log1pexp(-2f0)) ≈ -2f0
end

@testset "log1pmx" begin
    @test iszero(log1pmx(0.0))
    @test log1pmx(1.0) ≈ log(2.0) - 1.0
    @test log1pmx(2.0) ≈ log(3.0) - 2.0
end

@testset "logmxp1" begin
    @test iszero(logmxp1(1.0))
    @test logmxp1(2.0) ≈ log(2.0) - 1.0
    @test logmxp1(3.0) ≈ log(3.0) - 2.0
end

@testset "logsumexp" begin
    @test logaddexp(2.0, 3.0)     ≈ log(exp(2.0) + exp(3.0))
    @test logaddexp(10002, 10003) ≈ 10000 + logaddexp(2.0, 3.0)

    @test @inferred(logsumexp([1.0])) == 1.0
    @test @inferred(logsumexp((x for x in [1.0]))) == 1.0
    @test @inferred(logsumexp([1.0, 2.0, 3.0])) ≈ 3.40760596444438
    @test @inferred(logsumexp((1.0, 2.0, 3.0))) ≈ 3.40760596444438
    @test logsumexp([1.0, 2.0, 3.0] .+ 1000.) ≈ 1003.40760596444438

    @test @inferred(logsumexp([[1.0, 2.0, 3.0] [1.0, 2.0, 3.0] .+ 1000.]; dims=1)) ≈ [3.40760596444438 1003.40760596444438]
    @test @inferred(logsumexp([[1.0 2.0 3.0]; [1.0 2.0 3.0] .+ 1000.]; dims=2)) ≈ [3.40760596444438, 1003.40760596444438]
    @test @inferred(logsumexp([[1.0, 2.0, 3.0] [1.0, 2.0, 3.0] .+ 1000.]; dims=[1,2])) ≈ [1003.4076059644444]

    # check underflow
    @test logsumexp([1e-20, log(1e-20)]) ≈ 2e-20

    let cases = [([-Inf, -Inf], -Inf),   # correct handling of all -Inf
                 ([-Inf, -Inf32], -Inf), # promotion
                 ([-Inf32, -Inf32], -Inf32), # Float32
                 ([-Inf, Inf], Inf),
                 ([-Inf, 9.0], 9.0),
                 ([Inf, 9.0], Inf),
                 ([0, 0], log(2.0))] # non-float arguments
        for (arguments, result) in cases
            @test logaddexp(arguments...) ≡ result
            @test logsumexp(arguments) ≡ result
        end
    end

    @test isnan(logsubexp(Inf, Inf))
    @test isnan(logsubexp(-Inf, -Inf))
    @test logsubexp(Inf, 9.0) ≡ Inf
    @test logsubexp(-Inf, 9.0) ≡ 9.0
    @test logsubexp(1f2, 1f2) ≡ -Inf32
    @test logsubexp(0, 0) ≡ -Inf
    @test logsubexp(3, 2) ≈ 2.541324854612918108978

    # NaN propagation
    @test isnan(logaddexp(NaN, 9.0))
    @test isnan(logaddexp(NaN, Inf))
    @test isnan(logaddexp(NaN, -Inf))

    @test isnan(logsubexp(NaN, 9.0))
    @test isnan(logsubexp(NaN, Inf))
    @test isnan(logsubexp(NaN, -Inf))

    @test isnan(logsumexp([NaN, 9.0]))
    @test isnan(logsumexp([NaN, Inf]))
    @test isnan(logsumexp([NaN, -Inf]))

    # logsumexp with general iterables (issue #63)
    xs = range(-500, stop = 10, length = 1000)
    @test @inferred(logsumexp(x for x in xs)) == logsumexp(xs)
end

@testset "softmax" begin
    x = [1.0, 2.0, 3.0]
    r = exp.(x) ./ sum(exp.(x))
    @test softmax(x) ≈ r
    softmax!(x)
    @test x ≈ r

    x = [1, 2, 3]
    r = exp.(x) ./ sum(exp.(x))
    @test softmax(x) ≈ r
    @test eltype(softmax(x)) == Float64

    x = [1//2, 2//3, 3//4]
    r = exp.(x) ./ sum(exp.(x))
    @test softmax(x) ≈ r
    @test eltype(softmax(x)) == Float64

    x = Float32[1, 2, 3]
    r = exp.(x) ./ sum(exp.(x))
    @test softmax(x) ≈ r
    @test eltype(softmax(x)) == Float32
end


text:
"""
    ProjectiveTransforms(sz; [augmentations, buffered]) <: Encoding

Encoding for spatial data that resizes blocks to a common size `sz` and applies
projective augmentations.

Encodes all spatial blocks, preserving the block type:
- `Image{N}` -> `Image{N}`
- `Mask{N}` -> `Mask{N}`
- `Keypoints{N}` -> `Keypoints{N}`

The behavior differs based on the `context` of encoding:

{.tight}
- [`DLPipelines.Training`](#):
    1. Resizes the data so the smallest side equals
    a side length in `sz` while keeping the aspect ratio.
    2. Applies `augmentations`.
    3. Crops a random `sz`-sized portion of the data
- [`DLPipelines.Validation`](#):
    1. Resizes the data so the smallest side equals
    a side length in `sz` while keeping the aspect ratio.
    2. Crops a `sz`-sized portion from the center.
- [`DLPipelines.Inference`](#):
    1. Resizes the data so the smallest side equals
    a side length in `sz` while keeping the aspect ratio.
        Note that in this context, the data does not have
        size `sz`, since no cropping happens and aspect ratio
        is preserved.

`ProjectiveTransforms` is not limited to 2D data, and works on 3D data as well.
Note, however, that some transformations in `augs_projection` (rotation, warping, flipping)
are 2D only so `augs_projection` cannot be used for 3D data.

## Keyword arguments

- `augmentations::`[`DataAugmentation.Transform`](#)` = Identity()`: Projective
    augmentation to apply during training. See [`augs_projection`](#).
- `buffered = true`: Whether to use inplace transformations. Reduces memory usage.
- `sharestate = true`: Whether to use the same random state and bounds for all blocks
    in a sample

"""
struct ProjectiveTransforms{N,T} <: StatefulEncoding
    sz::NTuple{N,Int}
    buffered::Bool
    augmentations::Any
    tfms::T
    sharestate::Bool
end

function ProjectiveTransforms(
    sz;
    augmentations = Identity(),
    inferencefactor = 8,
    buffered = true,
    sharestate = true,
)
    traintfm = ScaleKeepAspect(sz) |> augmentations |> RandomCrop(sz) |> PinOrigin()
    validtfm = CenterResizeCrop(sz)
    tfms = (;
        training = buffered ? BufferedThreadsafe(traintfm) : traintfm,
        validation = buffered ? BufferedThreadsafe(validtfm) : validtfm,
        inference = ResizePadDivisible(sz, inferencefactor),
    )

    return ProjectiveTransforms(sz, buffered, augmentations, tfms, sharestate)
end

function encodestate(enc::ProjectiveTransforms{N}, context, blocks, obss) where {N}
    bounds = getsamplebounds(blocks, obss, N)
    tfm = _gettfm(enc.tfms, context)
    randstate = DataAugmentation.getrandstate(tfm)
    return bounds, randstate
end

_gettfm(tfms, ::Training) = tfms.training
_gettfm(tfms, ::Validation) = tfms.validation
_gettfm(tfms, ::Inference) = tfms.inference

function encodedblock(enc::ProjectiveTransforms{N}, block::Block) where {N}
    return isnothing(blockitemtype(block, N)) ? nothing : Bounded(block, enc.sz)
end

function encode(
    enc::ProjectiveTransforms{N},
    context,
    block::Block,
    obs;
    state = nothing,
) where {N}
    ItemType = blockitemtype(block, N)
    isnothing(ItemType) && return obs
    # only init state if block is encoded
    bounds, randstate =
        (isnothing(state) || !enc.sharestate) ? encodestate(enc, context, block, obs) :
        state
    # don't encode if bounds have wrong dimensionality
    bounds isa DataAugmentation.Bounds{N} || return obs

    tfm = _gettfm(enc.tfms, context)
    item = ItemType(obs, bounds)
    tobs = apply(tfm, item; randstate = randstate) |> itemdata
    return copy(tobs)
end

# ProjectiveTransforms is not invertible, hence no `decode` method!

# Conversion of `Block` to `DataAugmentation.Item`

"""
    blockitemtype(block, N)

Return a constructor for a `DataAugmentation.Item` that can be projected.
Return `nothing` by default, indicating that `block` cannot be turned into a
projectable item for bounds with dimensionality `N`.
For example, we have

    blockitemtype(Image{2}(), 2) -> DataAugmentation.Image

but

    blockitemtype(Image{3}(), 2) -> nothing
"""
blockitemtype(block::Block, n) = nothing
blockitemtype(block::Image{N}, n::Int) where {N} = N == n ? DataAugmentation.Image : nothing
function blockitemtype(block::Mask{N}, n::Int) where {N}
    return if N == n
        (obs, bounds) -> DataAugmentation.MaskMulti(obs, block.classes, bounds)
    else
        nothing
    end
end
blockitemtype(block::Keypoints{N}, n::Int) where {N} =
    N == n ? DataAugmentation.Keypoints : nothing
blockitemtype(block::WrapperBlock, n::Int) = blockitemtype(wrapped(block), n)

"""
    grabbounds(blocks, obss, N)

Looks through `blocks` for block data that carries `N`-dimensional
bounds information needed for projective transformations.
"""
function grabbounds(blocks::Tuple, obss::Tuple, N::Int)
    for (block, obs) in zip(blocks, obss)
        bounds = grabbounds(block, obs, N)
        !isnothing(bounds) && return bounds
    end
end

grabbounds(block::Image{N}, a, n) where {N} =
    N == n ? DataAugmentation.Bounds(size(a)) : nothing
grabbounds(block::Mask{N}, a, n) where {N} =
    N == n ? DataAugmentation.Bounds(size(a)) : nothing
grabbounds(block::WrapperBlock, a, n) = grabbounds(wrapped(block), a, n)

function getsamplebounds(blocks, obss, N::Int)
    bounds = grabbounds(blocks, obss, N)
    isnothing(bounds) && error("Could not detect $N-dimensional bounds needed for projective
transformations from blocks $(blocks)! Bounds can be grabbed from arrays
like `Image{$N}` block data.")
    return bounds
end

# Augmentation helper

"""
    augs_projection([; kwargs...])

Helper to create a set of projective transformations for image, mask
and keypoint data. Similar to fastai's
[`aug_transforms`](https://github.com/fastai/fastai/blob/bdc58846753c6938c63344fcaebea7149585fd5c/fastai/vision/augment.py#L946).

## Keyword arguments

- `flipx = true`: Whether to perform a horizontal flip with probability `1/2`. See [`FlipX`](#).
- `flipy = false`: Whether to perform a vertical flip with probability `1/2`. See [`FlipY`](#).
- `max_zoom = 1.5`: Maximum factor by which to zoom. Set to `1.` to disable. See [`Zoom`](#).
- `max_rotate = 10`: Maximum absolute degree by which to rotate. Set to `0.` to disable. See [`Rotate`](#).
- `max_warp = 0.05`: Intensity of corner warp. Set to `0.` to disable. See [`WarpAffine`](#).
"""
function augs_projection(;
    flipx = true,
    flipy = false,
    max_zoom = 1.5,
    max_rotate = 10.0,
    max_warp = 0.05,
)
    tfms = []

    flipx && push!(tfms, Maybe(FlipX()))
    flipy && push!(tfms, Maybe(FlipY()))
    max_warp > 0 && push!(tfms, WarpAffine(max_warp))
    max_rotate > 0 && push!(tfms, Rotate(max_rotate))
    push!(tfms, Zoom((1.0, max_zoom)))
    return DataAugmentation.compose(tfms...)
end

# Pretty-printing

function Base.show(io::IO, p::ProjectiveTransforms)
    show(io, ShowTypeOf(p))
    fields = (
        sz = ShowLimit(p.sz, limit = 80),
        buffered = ShowLimit(p.buffered, limit = 80),
        augmentations = ShowLimit(p.augmentations, limit = 80),
    )
    show(io, ShowProps(fields, new_lines = true))
end

# ## Tests

@testset "ProjectiveTransforms [encoding]" begin
    @testset "image" begin
        encoding = ProjectiveTransforms((32, 32))
        image = rand(RGB, 64, 96)
        block = FastAI.Image{2}()

        ## We run `ProjectiveTransforms` in the different [`Context`]s:
        imagetrain = encode(encoding, Training(), block, image)
        @test size(imagetrain) == (32, 32)

        imagevalid = encode(encoding, Validation(), block, image)
        @test size(imagevalid) == (32, 32)

        imageinference = encode(encoding, Inference(), block, image)
        @test size(imageinference) == (32, 48)

        ## During inference, the aspect ratio should stay the same
        @test size(image, 1) / size(image, 2) ==
              size(imageinference, 1) / size(imageinference, 2)
    end

    @testset "keypoints" begin
        encoding = ProjectiveTransforms((32, 48))
        ks = [SVector(0.0, 0), SVector(64, 96)]
        block = FastAI.Keypoints{2}(10)
        bounds = DataAugmentation.Bounds((1:32, 1:48))
        r = DataAugmentation.getrandstate(encoding.tfms.training)
        kstrain = encode(encoding, Training(), block, ks; state = (bounds, r))
        ksvalid = encode(encoding, Validation(), block, ks; state = (bounds, r))
        ksinference = encode(encoding, Inference(), block, ks; state = (bounds, r))

        @test kstrain == ksvalid == ksinference
    end

    @testset "image and keypoints" begin
        encoding = ProjectiveTransforms((32, 32))
        image = rand(RGB, 64, 96)
        ks = [SVector(0.0, 0), SVector(64, 96)]
        blocks = (FastAI.Image{2}(), FastAI.Keypoints{2}(10))

        @test_nowarn encode(encoding, Training(), blocks, (image, ks))
        @test_nowarn encode(encoding, Validation(), blocks, (image, ks))
        @test_nowarn encode(encoding, Inference(), blocks, (image, ks))
    end

    @testset "custom tests" begin
        enc = ProjectiveTransforms((32, 32), buffered = false)
        block = Image{2}()
        image = rand(RGB, 100, 50)

        testencoding(enc, block, image)
        @testset "randstate is shared" begin
            im1, im2 = encode(enc, Training(), (block, block), (image, image))
            @test im1 ≈ im2
        end

        @testset "don't transform data that doesn't need to be resized" begin
            imagesmall = rand(RGB, 32, 32)
            @test imagesmall ≈ encode(enc, Validation(), block, imagesmall)
        end

        @testset "3D" begin

            testencoding(
                ProjectiveTransforms((16, 16, 16)),
                Image{3}(),
                rand(RGB{N0f8}, 32, 24, 24),
            )
        end
    end

    #= depends on buffered interface
    @testset ExtendedTestSet "`run!`" begin
        encoding = ProjectiveTransforms((32, 32))
        image1 = rand(RGB, 64, 96)
        image2 = rand(RGB, 64, 96)
        buf = encode(encoding, Validation(), image1)
        cbuf = copy(buf)
        encode!(buf, encoding, Validation(), image1)
        @test buf ≈ cbuf

        encode!(buf, encoding, Validation(), image2)
        # buf should be modified on different input
        @test !(buf ≈ cbuf)
    end
    =#
end