text:
# we will want to overload these operators and functions for our custom types
import Base: ==, ≈, *, +, -

export normal_form, comm

Base.length(A::BaseOpProduct) = length(A.v)
Base.length(A::Union{ExpVal,Corr}) = length(A.ops)
Base.length(A::QuTerm) = length(A.bares) + mapreduce(length,+,A.expvals;init=0) + mapreduce(length,+,A.corrs;init=0)

Base.zero(::Type{QuExpr}) = QuExpr()
Base.zero(::QuExpr) = zero(QuExpr)
Base.one(::Type{QuExpr}) = QuExpr(QuTerm())
Base.one(::Type{QuTerm}) = QuExpr(QuTerm())
Base.one(::T) where T<:Union{QuTerm,QuExpr} = one(T)

function Base.:^(A::Union{QuTerm,QuExpr},n::Int)
    if n > 0
        prod(A for _=1:n)
    elseif n == 0
        one(A)
    else
        throw(ArgumentError("Only positive exponents supported for QuantumAlgebra objects."))
    end
end
Base.literal_pow(::typeof(^),A::Union{QuTerm,QuExpr},::Val{n}) where n = A^n

noaliascopy(A::Param) = Param(A.name,A.state,copy(A.inds))
# δ is isbitstype
noaliascopy(A::δ) = A
noaliascopy(A::BaseOperator) = BaseOperator(A.t,A.name,copy(A.inds))
noaliascopy(A::T) where T<:Union{ExpVal,Corr} = T(noaliascopy(A.ops))
noaliascopy(A::BaseOpProduct) = BaseOpProduct(noaliascopy(A.v))
noaliascopy(v::Vector{T}) where T<:Union{δ,Param,BaseOperator,ExpVal,Corr} = noaliascopy.(v)
noaliascopy(A::QuTerm) = QuTerm(A.nsuminds, noaliascopy(A.δs), noaliascopy(A.params),
                                noaliascopy(A.expvals), noaliascopy(A.corrs), noaliascopy(A.bares))
# IMPT: create a dictionary directly here, do not go through _add_sum_term
noaliascopy(A::QuExpr) = QuExpr(Dict{QuTerm,Number}(noaliascopy(t) => s for (t,s) in A.terms))

==(A::BaseOperator,B::BaseOperator) = A.t == B.t && A.name == B.name && A.inds == B.inds
==(A::BaseOpProduct,B::BaseOpProduct) = A.v == B.v
==(A::Param,B::Param) = A.name == B.name && A.state == B.state && A.inds == B.inds
==(A::T,B::T) where T<:Union{ExpVal,Corr} = A.ops == B.ops
==(A::QuTerm,B::QuTerm) = (A.nsuminds == B.nsuminds && A.δs == B.δs && A.params == B.params &&
                           A.expvals == B.expvals && A.corrs == B.corrs && A.bares == B.bares)
==(A::QuExpr,B::QuExpr) = A.terms == B.terms

function ≈(A::QuExpr,B::QuExpr)
    length(A.terms) != length(B.terms) && return false
    for (tA,sA) in A.terms
        sB = get(B.terms,tA) do
            return false
        end
        sA ≈ sB || return false
    end
    return true
end

function Base.hash(v::Vector{T},h::UInt) where T <: QuantumObject
    h = hash(length(v),h)
    for x in v
        h = hash(x,h)
    end
    h
end
function Base.hash(ind::QuIndex,h::UInt)
    h = hash(ind.sym,h)
    h = hash(ind.num,h)
    h
end
function Base.hash(name::QuOpName,h::UInt)
    h = hash(name.i,h)
    h
end
function Base.hash(A::BaseOperator,h::UInt)
    h = hash(UInt(A.t),h)
    h = hash(A.name,h)
    h = hash(A.inds,h)
    h
end
function Base.hash(A::Param,h::UInt)
    h = hash(A.name,h)
    h = hash(A.state,h)
    h = hash(A.inds,h)
    h
end
function Base.hash(A::BaseOpProduct,h::UInt)
    h = hash(A.v,h)
    h
end
function Base.hash(A::T,h::UInt) where T<:Union{ExpVal,Corr}
    h = hash(T,h)
    h = hash(A.ops,h)
    h
end
function Base.hash(A::QuTerm,h::UInt)
    h = hash(A.nsuminds,h)
    h = hash(A.δs,h)
    h = hash(A.expvals,h)
    h = hash(A.corrs,h)
    h = hash(A.bares,h)
    h
end
function Base.hash(A::QuExpr,h::UInt)
    h = hash(A.terms,h)
    h
end

# Base.cmp(::AbstractArray,::AbstractArray) uses isequal and isless, so doesn't shortcut if a,b are themselves arrays
# https://github.com/JuliaLang/julia/blob/539f3ce943f59dec8aff3f2238b083f1b27f41e5/base/abstractarray.jl
function recursive_cmp(A, B)
    for (a, b) in zip(A, B)
        cc = recursive_cmp(a,b)
        cc == 0 || return cc
    end
    return cmp(length(A), length(B))
end

@inline recursive_cmp(A::BaseOperator,B::BaseOperator) = cmp(A,B)
@inline recursive_cmp(A::Param,B::Param) = cmp(A,B)
@inline recursive_cmp(A::δ,B::δ) = cmp(A,B)
@inline recursive_cmp(A::BaseOpProduct,B::BaseOpProduct) = (cc = cmp(length(A),length(B)); cc == 0 ? recursive_cmp(A.v,B.v) : cc)
@inline recursive_cmp(A::T,B::T) where T<:Union{ExpVal,Corr}  = recursive_cmp(A.ops,B.ops)

macro _cmpAB(member,rec=true)
    fun = rec ? :recursive_cmp : :cmp
    esc(quote
        cc = $fun(A.$member,B.$member)
        cc == 0 || return cc
    end)
end

@inline function Base.cmp(A::BaseOperator,B::BaseOperator)
    @_cmpAB t false
    @_cmpAB name false
    @_cmpAB inds false
    # they are equal
    return 0
end

@inline function Base.cmp(A::Param,B::Param)
    @_cmpAB name false
    @_cmpAB state false
    @_cmpAB inds false
    # they are equal
    return 0
end

@inline function Base.cmp(A::δ,B::δ)
    @_cmpAB iA false
    @_cmpAB iB false
    # they are equal
    return 0
end

function Base.cmp(A::QuTerm,B::QuTerm)
    # only evaluate each part that we need for each step of the comparison to avoid unnecessary work
    # order first by number of operators (also within expectation values)
    cc = cmp(length(A), length(B))
    cc == 0 || return cc
    # then order by operators, expectation values, correlations, params, deltas, nsuminds
    @_cmpAB bares
    @_cmpAB expvals
    @_cmpAB corrs
    @_cmpAB params
    @_cmpAB δs
    @_cmpAB nsuminds false
    # they are equal
    return 0
end

@inline Base.isless(A::BaseOperator,B::BaseOperator) = cmp(A,B) < 0
@inline Base.isless(A::Param,B::Param) = cmp(A,B) < 0
@inline Base.isless(A::δ,B::δ) = cmp(A,B) < 0
@inline Base.isless(A::BaseOpProduct,B::BaseOpProduct) = recursive_cmp(A,B) < 0
@inline Base.isless(A::T,B::T) where T<:Union{ExpVal,Corr} = recursive_cmp(A,B) < 0
@inline Base.isless(A::QuTerm,B::QuTerm) = cmp(A,B) < 0

comm(A,B) = A*B - B*A

Base.adjoint(A::BaseOperator) = BaseOperator(BaseOpType_adj[Int(A.t)],A.name,A.inds)
Base.adjoint(A::BaseOpProduct) = BaseOpProduct(adjoint.(view(A.v, lastindex(A.v):-1:1)))
Base.adjoint(A::T) where {T<:Union{ExpVal,Corr}} = T(adjoint(A.ops))
Base.adjoint(A::Param) = A.state=='r' ? A : Param(A.name,A.state=='n' ? 'c' : 'n',A.inds)
function Base.adjoint(A::QuTerm)
    B = QuTerm(A.nsuminds,A.δs,adjoint.(A.params),adjoint.(A.expvals),adjoint.(A.corrs),adjoint(A.bares))
    B.nsuminds>0 ? reorder_suminds()(B) : B
end
Base.adjoint(A::QuExpr) = QuExpr((adjoint(t),adjoint(s)) for (t,s) in A.terms)

is_normal_form(A::Union{ExpVal,Corr}) = is_normal_form(A.ops)
is_normal_form(ops::BaseOpProduct) = is_normal_form(ops.v)
function is_normal_form(v::Vector{BaseOperator})
    issorted(v) || return false
    # now check that no contractions occur
    for k in 1:length(v)-1
        O = v[k]
        dosimplify, fac, op = _contract(O,v[k+1])
        dosimplify && return false
        if O.t in (TLSx_,TLSy_,TLSz_)
            # lookahead while we can commute through
            for kp = k+2:length(v)
                v[kp].t in (TLSx_,TLSy_,TLSz_) || break
                _exchange(v[kp-1],v[kp]) == (1,nothing) || break
                dosimplify, fac, op = _contract(O,v[kp])
                dosimplify && return false
            end
        end
    end
    return true
end
function is_normal_form(t::QuTerm)
    (issorted(t.params) && issorted(t.expvals) && issorted(t.corrs)) || return false
    is_normal_form(t.bares) || return false
    for EV in t.expvals
        is_normal_form(EV) || return false
    end
    for CO in t.corrs
        is_normal_form(CO) || return false
    end

    # ensure that sum indices are ordered
    changeable_indices = indices(t,false)
    last_sumind = sumindex(0).num
    for ind in Base.Iterators.filter(issumindex,changeable_indices)
        # sumindex have to be increasing without jumps
        ind.num ≤ last_sumind + oneunit(last_sumind) || return false
        last_sumind = max(ind.num,last_sumind)
    end

    isempty(t.δs) && return true

    # all indices that are not in δs
    changeable_indices_set = Set{QuIndex}(changeable_indices)
    # accumulate indices that δs want to change (to compare with later δs)
    replace_inds = Set{QuIndex}()
    for (iδ,dd) in enumerate(t.δs)
        iA,iB = dd.iA, dd.iB
        # δ is not ordered (or disappears if iA == iB)
        iA ≥ iB && return false
        # the whole term disappears in cleanup
        isintindex(iA) && isintindex(iB) && return false
        (issumindex(iA) || issumindex(iB)) && return false
        (iA in replace_inds || iB in replace_inds) && return false
        if iδ > 1
            (dd > @inbounds t.δs[iδ-1]) || return false
        end
        iB in changeable_indices_set && return false
        push!(replace_inds, iB)
    end
    return true
end
function is_normal_form(A::QuExpr)
    for t in keys(A.terms)
        is_normal_form(t) || return false
    end
    return true
end

function _normalize_without_commutation(A::QuTerm)::Union{QuTerm,Nothing}
    # first, clean up the δs
    if isempty(A.δs)
        # we always want _normalize_without_commutation to return a new object so we can modify it later
        A = noaliascopy(A)
    else
        δs = sort(A.δs)
        #println("in _normalize_without_commutation for A = $A, starting with δs = $δs")
        replacements = Dict{QuIndex,QuIndex}()
        delsuminds = IndexInt[]
        iwrite = 1
        for dd in δs
            iA = get(replacements,dd.iA,dd.iA)
            iB = get(replacements,dd.iB,dd.iB)
            if iA == iB
                # this delta just gives one
                continue
            elseif isintindex(iA) && isintindex(iB)
                # the term is zero
                return nothing
            elseif iA > iB
                iA, iB = iB, iA
            end
            if issumindex(iB) # one sum disappears, and the delta does as well
                replacements[iB] = iA
                # we will later need to shift all larger sumindices down by one
                push!(delsuminds,iB.num)
            elseif issumindex(iA) # one sum disappears, and the delta does as well
                replacements[iA] = iB
                # we will later need to shift all larger sumindices down by one
                push!(delsuminds,iA.num)
            else
                # otherwise, replace the larger index by the smaller one
                replacements[iB] = iA
                # we keep the delta, but include any possible replacements made
                δs[iwrite] = δ(iA,iB)
                iwrite += 1
            end
        end
        resize!(δs,iwrite-1)
        # after replacements, δs can be out of order
        sort!(δs)
        if !isempty(delsuminds)
            sumrepls = Dict{QuIndex,QuIndex}()
            sort!(delsuminds)
            for (shift,startind) in enumerate(delsuminds)
                endind = shift == length(delsuminds) ? A.nsuminds : delsuminds[shift+1]-1
                for num = startind+1:endind
                    sumrepls[sumindex(num)] = sumindex(num-shift)
                end
            end
            # we want to first apply replacements and then sumrepls
            # this is the same as replacing results from replacements with sumrepls,
            # and then applying sumrepls on indices that have not been replaced yet
            for (iold,inew) in replacements
                x = get(sumrepls,inew,inew)
                x != inew && (replacements[iold] = x)
            end
            # sumrepl only for indices not in replacements
            merge!(sumrepls,replacements)
            replacements = sumrepls
        end
        nsuminds = A.nsuminds-length(delsuminds)
        f = replace_inds(replacements)
        #println("in _normalize_without_commutation, deleting sum indices $delsuminds and doing replacements $replacements.")
        A = QuTerm(nsuminds,δs,f.(A.params),f.(A.expvals),f.(A.corrs),f(A.bares))
    end
    # also sort all the commuting terms
    sort!(A.params)
    sort!(A.expvals)
    sort!(A.corrs)
    A
end

# levicivita_lut[a,b] contains the Levi-Cevita symbol ϵ_abc
# for c=6-a-b, i.e, when a,b,c is a permutation of 1,2,3
const levicivita_lut = [0 1 -1; -1 0 1; 1 -1 0]
function ϵ_ab(A::BaseOperator,B::BaseOperator)
    # a+b+c == 6 (since a,b,c is a permutation of 1,2,3)
    a = Int(A.t) - Int(TLSx_) + 1
    b = Int(B.t) - Int(TLSx_) + 1
    c = BaseOpType(Int(TLSx_) - 1 + (6 - a - b))
    s = @inbounds levicivita_lut[a,b]
    c, s
end

struct ExchangeResult
    pref::Complex{Rational{Int}}
    δs::Vector{δ}
    op::Union{Nothing,BaseOperator}
end

# rewrite B A as x A B + y, with A < B in our ordering
# returns (x::Int,y::Union{ExchangeResult,Nothing})
function _exchange(A::BaseOperator,B::BaseOperator)::Tuple{Int,Union{ExchangeResult,Nothing}}
    if A.t == B.t
        # these operators always commute with the same type
        A.t in (BosonDestroy_,BosonCreate_,TLSCreate_,TLSDestroy_,TLSx_,TLSy_,TLSz_) && return (1,nothing)
        # these operators anticommute if they refer to the same species (name), and commute otherwise
        A.t in (FermionDestroy_,FermionCreate_) && return (A.name == B.name ? -1 : 1,nothing)
    end

    # different types of operators commute
    if A.t in (BosonDestroy_,BosonCreate_) && B.t in (FermionDestroy_,FermionCreate_,TLSx_,TLSy_,TLSz_,TLSDestroy_,TLSCreate_)
        return (1,nothing)
    elseif A.t in (FermionDestroy_,FermionCreate_) && B.t in (BosonDestroy_,BosonCreate_,TLSx_,TLSy_,TLSz_,TLSDestroy_,TLSCreate_)
        return (1,nothing)
    elseif A.t in (TLSx_,TLSy_,TLSz_,TLSDestroy_,TLSCreate_) && B.t in (BosonDestroy_,BosonCreate_,FermionDestroy_,FermionCreate_)
        return (1,nothing)
    end

    # a(i) adag(j) = adag(j) a(i) + δij
    if A.t == BosonCreate_ && B.t == BosonDestroy_
        if A.name == B.name && (dd = δ(A.inds,B.inds)) !== nothing
            return (1, ExchangeResult(1,dd,nothing))
        else
            return (1, nothing)
        end
    end

    # f(i) fdag(j) = -fdag(j) f(i) + δij
    if A.t == FermionCreate_ && B.t == FermionDestroy_
        if A.name == B.name
            dd = δ(A.inds,B.inds)
            return (-1, dd === nothing ? nothing : ExchangeResult(1,dd,nothing))
        else
            return (1, nothing)
        end
    end

    if A.t == TLSCreate_ && B.t == TLSDestroy_
        if A.name != B.name || (dd = δ(A.inds,B.inds)) === nothing
            return (1,nothing)
        elseif isempty(dd)
            # indices were all the same, σ- σ+ = 1 - σ+ σ-
            return (-1, ExchangeResult(1, dd, nothing))
        else
            # σ-_i σ+_j = σ+_j σ-_i - δij σz_i = σ+_j σ-_i + δij (1 - 2 σ+_i σ-_i)
            # note that we return a TLSz to "fit" in ExchangeResult, but need to undo that later on
            # also note that we pass "+TLSz" instead of "-TLSz" so we do not have to flip the sign of the "1" term in the sort algorithm
            return (1, ExchangeResult(1, dd, TLSz(A.name,A.inds)))
        end
    end

    # B A as x A B + y, with A < B
    if A.t in (TLSx_,TLSy_,TLSz_) && B.t in (TLSx_,TLSy_,TLSz_)
        if A.name != B.name || (dd = δ(A.inds,B.inds)) === nothing
            return (1,nothing)
        else
            # we need ϵbac
            c, s = ϵ_ab(B,A)
            if isempty(dd)
                # indices were all the same,
                # σb σa = δab + i ϵbac σc
                # since A < B, we know that a != b
                return (0, ExchangeResult(im*s,dd,BaseOperator(c,A.name,A.inds)))
            else
                # [σb_i,σa_j] = 2i εbac δij σc_i =>
                # σb_j σa_i = σa_i σb_j + 2i εbac δij σc_i
                return (1, ExchangeResult(2im*s,dd,BaseOperator(c,A.name,A.inds)))
            end
        end
    end

    if (A.t == TLSCreate_ && B.t in (TLSx_,TLSy_,TLSz_)) || (A.t in (TLSx_,TLSy_,TLSz_) && B.t == TLSDestroy_)
        throw(ArgumentError("QuantumAlgebra currently does not support mixing 'normal' (x,y,z) Pauli matrices with jump operators (+,-)."))
    end

    error("_exchange should never reach this! A=$A, B=$B.")
end

const ComplexInt = Complex{Int}

function _contract(A::BaseOperator,B::BaseOperator)::Tuple{Bool,ComplexInt,Union{BaseOperator,Nothing}}
    if A.t in (TLSCreate_,TLSDestroy_,FermionDestroy_,FermionCreate_) && A == B
        return (true,zero(ComplexInt),nothing)
    elseif A.t in (TLSx_,TLSy_,TLSz_) && B.t in (TLSx_,TLSy_,TLSz_) && A.name == B.name && A.inds == B.inds
        # σa σb = δab + i ϵabc σc
        if B.t == A.t
            return (true,one(ComplexInt),nothing)
        else
            c, s = ϵ_ab(A,B)
            return (true, ComplexInt(0,s), BaseOperator(c,A.name,A.inds))
        end
    else
        return (false,zero(ComplexInt),nothing)
    end
end

function normal_order!(ops::BaseOpProduct,term_collector)
    # do an insertion sort to get to normal ordering
    # reference: https://en.wikipedia.org/wiki/Insertion_sort
    A = ops.v
    prefactor::ComplexInt = one(ComplexInt)
    for i = 2:length(A)
        j = i
        while j>1 && A[j]<A[j-1]
            # need to commute A[j-1] and A[j]
            pp, exc_res = _exchange(A[j],A[j-1])
            #println("exchanging A[$j] = $(A[j]) and A[$kk] = $(A[kk]) gave result: $pp, $exc_res.")
            if exc_res !== nothing
                if exc_res.op === nothing
                    onew = BaseOpProduct([A[1:j-2]; A[j+1:end]])
                elseif A[j].t == TLSCreate_
                    @assert exc_res.op.t == TLSz_
                    # we got σz_i, have to replace it by (1 - 2 σ+_i σ-_i) (which is really -σz, see explanation in _exchange)
                    onew = BaseOpProduct([A[1:j-2]; TLSCreate(exc_res.op.name,exc_res.op.inds); TLSDestroy(exc_res.op.name,exc_res.op.inds); A[j+1:end]])
                    t = QuTerm(exc_res.δs, onew)
                    _add_sum_term!(term_collector,t,-2exc_res.pref*prefactor)
                    # this one gives the "1" term that is used below
                    onew = BaseOpProduct([A[1:j-2]; A[j+1:end]])
                else
                    onew = BaseOpProduct([A[1:j-2]; exc_res.op; A[j+1:end]])
                end
                t = QuTerm(exc_res.δs, onew)
                _add_sum_term!(term_collector,t,exc_res.pref*prefactor)
                #println("adding term = $(exc_res.pref*prefactor) * $t, term_collector = $(term_collector)")
            end
            # only modify prefactor after the exchange
            prefactor *= pp
            iszero(prefactor) && return prefactor

            # now finally exchange the two
            A[j-1], A[j] = A[j], A[j-1]
            j -= 1
        end
    end
    # check if we have any products that simplify
    did_contractions = false
    k = 1
    while k < length(A)
        dosimplify, fac, op = _contract(A[k],A[k+1])
        if dosimplify
            did_contractions = true
            prefactor *= fac
            iszero(prefactor) && return prefactor
            if op === nothing
                deleteat!(A,(k,k+1))
                #println("deleted indices $(k-1) and $k, A is now: $A, prefactor is now: $prefactor")
            else
                A[k] = op
                deleteat!(A,k+1)
                #println("deleted index $k and replaced $(k-1) by $op. A is now: $A, prefactor is now: $prefactor")
            end
            # we replaced [...,A[k-1],A[k],A[k+1],A[k+2],...]
            # with        [...,A[k-1],A[k+2],...]
            # or with     [...,A[k-1],Anew,A[k+2]...]
            # so reduce k by one to compare the new operator at k with A[k-1] as well
            # (with the ordering we have, it's probably impossible for A[k-1] to contract
            # with Anew if it didn't contract with A[k], but let's be safe)
            k > 1 && (k -= 1)
        else
            if A[k].t in (TLSx_,TLSy_,TLSz_)
                # lookahead while we can commute through
                # we already checked with k+1, so start at k+2
                for kp = k+2:length(A)
                    A[kp].t in (TLSx_,TLSy_,TLSz_) || break
                    _exchange(A[kp-1],A[kp]) == (1,nothing) || break
                    dosimplify, fac, op = _contract(A[k],A[kp])
                    if dosimplify
                        did_contractions = true
                        prefactor *= fac
                        iszero(prefactor) && return prefactor
                        if op === nothing
                            deleteat!(A,(k,kp))
                            #println("deleted indices $(k-1) and $k, A is now: $A, prefactor is now: $prefactor")
                        else
                            A[k] = op
                            deleteat!(A,kp)
                            #println("deleted index $k and replaced $(k-1) by $op. A is now: $A, prefactor is now: $prefactor")
                        end
                        # k -= 2 here because we have k += 1 just below
                        k > 1 && (k -= 2)
                        break
                    end
                end
            end
            k += 1
        end
    end
    if did_contractions
        # do another round of normal ordering since this could have changed after contractions
        prefactor *= normal_order!(ops,term_collector)
    end
    prefactor
end

normal_order!(A::Union{Corr,ExpVal},term_collector) = normal_order!(A.ops,term_collector)
function normal_order!(v::Vector{T},commterms,pref=1) where {T<:Union{Corr,ExpVal}}
    for (ii,EC) in enumerate(v)
        newterms = QuExpr()
        pp = normal_order!(EC,newterms)
        for (t,s) in newterms.terms
            #@assert iszero(t.nsuminds) && isempty(t.params) && isempty(t.expvals) && isempty(t.corrs)
            nv = [isempty(t.bares) ? T[] : T(t.bares); v[1:ii-1]; v[ii+1:end]]
            nt = QuTerm(t.δs, noaliascopy(nv))
            # include current prefactor here
            _add_sum_term!(commterms,nt,pref*s)
        end
        # only modify prefactor after the exchange
        pref *= pp
    end
    filter!(EC -> !isempty(EC.ops), v)
    sort!(v)
    pref
end

# iterate over a single term t with prefactor s and a sum as if they were in the same sum
termsumiter(t::QuTerm,s,sum::QuExpr) = Base.Iterators.flatten((((t,s),), sum.terms))
termsumiter(t,s,sum::QuExpr) = termsumiter(QuTerm(t),s,sum)

function _add_with_normal_order!(A::QuExpr,t::QuTerm,s)
    # no need to do anything since everything is in normal order
    # in particular, no copy necessary!
    if is_normal_form(t)
       _add_sum_term!(A,t,s)
       return
    end

    t = _normalize_without_commutation(t)
    t === nothing && return
    newbareterms = QuExpr()
    newexpvterms = QuExpr()
    newcorrterms = QuExpr()
    prefbare = normal_order!(t.bares,   newbareterms)
    prefexpv = normal_order!(t.expvals, newexpvterms)
    prefcorr = normal_order!(t.corrs,   newcorrterms)

    # normal_form(t) = t.prefs * (prefexpv*tN.ev + ev_new) * (prefcorr*tN.co + co_new) * (prefbare*tN.bare + bare_new)
    # t has already been transformed to tN in-place
    pref = prefbare*prefexpv*prefcorr
    if !iszero(pref)
        if iszero(t.nsuminds)
            # should be guaranteed to be in normal order
            _add_sum_term!(A,t,s*pref)
        else
            # might need to reorder sum indices, which can "break" our normal order
            tn = reorder_suminds()(t)
            # this will check again if it is in normal form to make sure
            _add_with_normal_order!(A,tn,s*pref)
        end
    end

    isempty(newexpvterms) && isempty(newcorrterms) && isempty(newbareterms) && return

    firstterm = true
    for (tev,sev) in termsumiter(t.expvals, prefexpv, newexpvterms)
        #@assert iszero(tev.nsuminds) && isempty(tev.params) && isempty(tev.corrs) && isempty(tev.bares)
        for (tco,sco) in termsumiter(t.corrs, prefcorr, newcorrterms)
            #@assert iszero(tco.nsuminds) && isempty(tco.params) && isempty(tco.expvals) && isempty(tco.bares)
            for (tba,sba) in termsumiter(t.bares, prefbare, newbareterms)
                #@assert iszero(tba.nsuminds) && isempty(tba.params) && isempty(tba.expvals) && isempty(tba.corrs)

                if firstterm
                    # the first term has been added above in _add_sum_term!
                    firstterm = false
                    continue
                end

                pref = sev*sco*sba
                iszero(pref) && continue
                # IMPTE: nt is cleaned afterwards with _normalize_without_commutation,
                # which creates copies of everything, so we can reuse arrays here
                nt = QuTerm(t.nsuminds,[t.δs;tev.δs;tco.δs;tba.δs],t.params,tev.expvals,tco.corrs,tba.bares)
                _add_with_normal_order!(A,nt,s*pref)
            end
        end
    end
end

function normal_form(A::QuExpr)
    An = QuExpr()
    for (t,s) in A.terms
        _add_with_normal_order!(An,t,s)
    end
    An
end

function *(A::QuTerm,B::QuTerm)
    nsuminds = A.nsuminds + B.nsuminds
    if A.nsuminds > 0 && B.nsuminds > 0
        if A.nsuminds > B.nsuminds
            A = shift_sumind(B.nsuminds)(A)
        else
            B = shift_sumind(A.nsuminds)(B)
        end
    end
    if nsuminds > 0
        f = reorder_suminds()
        apflat = (args...) -> f.(Iterators.flatten(args))
    else
        apflat = (args...) -> vcat(args...)
    end
    δs = apflat(A.δs, B.δs)
    params = apflat(A.params, B.params)
    expvals = apflat(A.expvals, B.expvals)
    corrs = apflat(A.corrs, B.corrs)
    barevec = apflat(A.bares.v, B.bares.v)
    QuTerm(nsuminds,δs,params,expvals,corrs,BaseOpProduct(barevec))
end

*(A::Union{QuExpr,QuTerm}) = A
function *(A::QuExpr,B::QuExpr)
    QuExpr((tA*tB,sA*sB) for ((tA,sA),(tB,sB)) in Iterators.product(A.terms,B.terms))
end
*(A::Number,B::QuExpr) = QuExpr((tB,A*sB) for (tB,sB) in B.terms)
*(B::QuExpr,A::Number) = A*B
Base.:/(A::QuExpr,B::Number) = inv(B)*A

+(A::QuExpr) = A
function +(A::QuExpr,B::QuExpr)
    S = copy(A)
    for (t,s) in B.terms
        _add_with_auto_order!(S,t,s)
    end
    S
end
function +(A::QuExpr,B::Number)
    S = copy(A)
    _add_with_auto_order!(S,QuTerm(),B)
    S
end
+(B::Number,A::QuExpr) = A+B

-(A::QuExpr) = -1 * A
function -(A::QuExpr,B::QuExpr)
    S = copy(A)
    for (t,s) in B.terms
        _add_with_auto_order!(S,t,-s)
    end
    S
end
-(B::Number,A::QuExpr) = B + (-A)
-(A::QuExpr,B::Number) = A + (-B)


text:
function next_prime(n)

  p = n + 1

  if p == 2
    return 2
  end

  if mod(p, 2) == 0
    p = p + 1
  end

  while length(factors(p)) > 2
    p = p + 2
  end

  return p

end

ext: jl


text:
function makeneighbourhood(dims::Int)::Array{CartesianIndex}
    filter(c -> !all(==(0), c[i] for i in 1:length(c)), CartesianIndices(tuple(repeat([-1:1], dims)...)))
end

mutable struct CellularAutomaton
    grid::Set{CartesianIndex}
end

function rules(grid::Set{CartesianIndex}, location::CartesianIndex, neighbourhood=makeneighbourhood(3))::Bool
    neighbourcount = count(neighbour in grid for neighbour in [location] .+ neighbourhood)
    neighbourcount == 3 || neighbourcount == 2 && location in grid
end

function CellularAutomaton(initstate::String, active::Char; rowseparator='\n', size=3)
    grid = Set{CartesianIndex}()

    for (j, row) in Iterators.enumerate(split(initstate, rowseparator))
        for (i, cell) in Iterators.enumerate(row)
            if cell == active
                push!(grid, CartesianIndex((i, j, repeat([0], size - 2)...)))
            end
        end
    end

    CellularAutomaton(grid)
end

import Base.iterate
function iterate(cellularautomaton::CellularAutomaton, state=cellularautomaton; neighbourhood=makeneighbourhood(3))::Tuple{Set{CartesianIndex},CellularAutomaton}
    grid = Set{CartesianIndex}()
    unitvector = CartesianIndex(tuple(repeat([1], length(first(cellularautomaton.grid)))...))
    for location in (minimum(state.grid) - unitvector:maximum(state.grid) + unitvector)
        if rules(state.grid, location, neighbourhood)
            push!(grid, location)
        end
    end
    state.grid = grid
    (grid, state)
end

function solve(init::String, neighbourhood::Array{CartesianIndex})::Int64
    auto = CellularAutomaton(init, '#', size=length(first(neighbourhood)))
    for _ in 1:6
        iterate(auto; neighbourhood=neighbourhood)
    end
    length(auto.grid)
end

using Test

exampleinput = ".#.\n..#\n###"
input = string(strip(read("input.txt", String)))

moore3 = makeneighbourhood(3)
moore4 = makeneighbourhood(4)

println("Part one")
@test solve(exampleinput, moore3) == 112
@show solve(input, moore3)

println("Part two")
@test solve(exampleinput, moore4) == 848
@show solve(input, moore4)

max_stars_repo_name: stellartux/advent-of-code

max_issues_repo_name: stellartux/advent-of-code

max_forks_repo_name: stellartux/advent-of-code


text:
import LinearAlgebra: Matrix, Cholesky, cholesky, mul!, ldiv!, Hermitian, Factorization
import Base: size, getindex, enumerate, iterate, copy, fill!
import SparseArrays: SparseVector, SparseMatrixCSC, nonzeros, rowvals, getcolptr, sparse, sparsevec, nnz
# Defining an alias for matrices that are realized as resized contiguous view into a buffer
const ContiguousBufferedMatrix{T} = Base.ReshapedArray{T,2,SubArray{Float64,1,Array{T,1},Tuple{UnitRange{Int64}},true},Tuple{}}

# Empty constructor used for initializing arrays
function ContiguousBufferedMatrix{T}() where T
    return reshape(view(zeros(T, 0), 1:0), 0, 0)
end

abstract type  AbstractSupernodalArray{Tv, N} end
abstract type  AbstractSupernodalSparseArray{Tv, N} <: AbstractSupernodalArray{Tv,N} end

function data(in::AbstractSupernodalArray)
    return in.data
end

function nonzeros(in::AbstractSupernodalSparseArray)
    return nonzeros(data(in))
end

function rowvals(in::AbstractSupernodalSparseArray)
    return rowvals(data(in))
end

function getcolptr(in::AbstractSupernodalSparseArray)
    return getcolptr(data(in))
end

function getindex(in::AbstractSupernodalSparseArray, I...)
    return getindex(data(in), I...)
end

function enumerate(in::AbstractSupernodalArray)
    enumerate(data(in))
end

# A supernodal column that contains
struct SupernodalVector{RT} <: AbstractSupernodalArray{RT,1}
    data::Vector{ContiguousBufferedMatrix{RT}}
    buffer::Vector{RT}
    row_supernodes::Vector{Vector{Int}}
end

function copy(in::SupernodalVector{RT}) where RT
    return SupernodalVector(Matrix(in), copy(in.row_supernodes))
end

function fill!(in::SupernodalVector, a)
    fill!(in.buffer, a)
end

# Constructs a supernodal column from a dense matrix and supernodal ids
function SupernodalVector(𝐌, row_supernodes::AbstractVector{<:AbstractVector{Int}})
    buffer = zeros(eltype(𝐌), length(𝐌))
    data = Vector{ContiguousBufferedMatrix{eltype(𝐌)}}(undef, length(row_supernodes))
    N = size(𝐌, 2)
    start_index = 1
    for (k, node) in enumerate(row_supernodes)
        data[k] = reshape(view(buffer, start_index : (start_index + length(node) * N - 1)), length(node), N)
        start_index += length(node) * N
    end
    for (k, block) in enumerate(data)
        # Iterate over all indices in the block matrix
        for c in CartesianIndices(block)
            # Transform the i-index from the local index c[1] to the global index, as read of from the corresponding row_supernode
            block[c] = 𝐌[row_supernodes[k][c[1]], c[2]]
        end
    end
    return SupernodalVector{eltype(𝐌)}(data, buffer, row_supernodes)
end

function size(in::SupernodalVector)
    return (sum(length.(in.row_supernodes)), size(first(in.data), 2))
end

function size(in::SupernodalVector, dim)
    return size(in)[dim]
end

# A sparse supernodal column, which will be used for storing the basis functions in supernodal form
struct SupernodalSparseVector{RT} <: AbstractSupernodalSparseArray{RT,1}
    data::SparseVector{ContiguousBufferedMatrix{RT}, Int}
    buffer::Vector{RT}
    row_supernodes::Vector{Vector{Int}}
end

function size(in::SupernodalSparseVector)
    return (sum(length.(in.row_supernodes)), size(first(in.data.nzval), 2))
end

function size(in::SupernodalSparseVector, dim)
    return (sum(length.(in.row_supernodes)), size(first(in.data.nzval), 2))[dim]
end

# Construct a supernodal sparse vector from a sparse matrix
function SupernodalSparseVector(mat::SparseMatrixCSC, row_supernodes::Vector{Vector{Int}})
    M, N = size(mat)
    @assert sum(length.(row_supernodes)) == M
    @assert sort(vcat(row_supernodes...)) == 1 : M
    # A vector such that row_supernodes[supernode_dict[k][1]][supernode_dict[k][2]] = k
    supernode_dict = Vector{Tuple{Int,Int}}(undef, M)
    for (i, node) in enumerate(row_supernodes)
        for (j, k) in enumerate(node)
            supernode_dict[k] = (i, j)
        end
    end
    I_el, J_el, S_el = findnz(mat)
    # Re-expressing the dofs in terms of the supernodes
    I_super = supernode_dict[I_el]
    # keeping track of those supernodes that appear in the sparsity pattern
    # We preallocate the array that will hold the matrix elements of the sparse vector
    S_super = Vector{ContiguousBufferedMatrix{eltype(mat)}}(undef, length(unique(getindex.(I_super, 1))))
    # We sum up the lengths of the supernodes in the sparsity pattern and multiply them with the number of entries per row.
    buffer = Vector{eltype(mat)}(undef, sum(length.(row_supernodes[unique(getindex.(I_super, 1))])) * N)
    # we are sorting the entries according by the supernode that they form part of
    sp = sortperm(I_super)
    # reordering the sparsity entries
    I_super = I_super[sp]
    I_el = I_el[sp]
    J_el = J_el[sp]
    S_el = S_el[sp]

    # We now extract those supernodes that appear in the sparsity pattern and recreate the dict, with the first index indicating the position of the supernode among the used supernodes
    used_supernode_indices = unique(getindex.(I_super, 1))
    lookup = Dict(used_supernode_indices .=> 1 : length(used_supernode_indices))
    used_row_supernodes = row_supernodes[used_supernode_indices]

    # Setting up the buffer and nonzeros
    offset = 0
    for (k, node) in enumerate(used_row_supernodes)
        S_super[k] = reshape(view(buffer, offset .+ (1 : (length(node)* N))), length(node), N)
        # filling the new array with zeros
        S_super[k] .= 0
        # updating the offset
        offset += length(node) * N
    end

    for (i_super, j, s) in zip(I_super, J_el, S_el)
        # in the supernode \# i_super[1], we add s to the entry (i_super[2], j)
        S_super[lookup[i_super[1]]][i_super[2], j] += s
    end
    return SupernodalSparseVector{eltype(mat)}(sparsevec(used_supernode_indices, S_super), buffer, row_supernodes)
end

function SupernodalSparseVector(node::SuperNodeBasis, row_supernodes::AbstractVector{<:AbstractVector{Int}})
    return SupernodalSparseVector(hcat(coefficients.(basis_functions(node))...), row_supernodes)
end

# inverse of the constructor of the sparsesupernodalvector
function SparseMatrixCSC(in::SupernodalSparseVector)
    out_I = Vector{Int}(undef, 0)
    out_J = Vector{Int}(undef, 0)
    out_S = Vector{eltype(in.buffer)}(undef, 0)
    for (i, mat) in zip(findnz(data(in))...)
        node = in.row_supernodes[i]
        for (i, j) in Tuple.(CartesianIndices(mat))
            push!(out_I, node[i])
            push!(out_J, j)
            push!(out_S, mat[i, j])
        end
    end
    return sparse(out_I, out_J, out_S, sum(length.(in.row_supernodes)), size(first(data(in).nzval), 2))
end

# Constructs a supernodal column from a dense matrix and a list of domain supernodes
function SupernodalVector(𝐌, row_supernodes::AbstractVector{<:SuperNodeDomain})
    return SupernodalVector(𝐌, id.(domains(row_supernodes)))
end

function Matrix(col::SupernodalVector)
    M = sum(length.(col.row_supernodes))
    N = Int(length(col.buffer) / M)
    out = zeros(eltype(col.buffer), M, N)
    # Iterate over all block matrices in data
    for (k, block) in enumerate(col)
        # Iterate over all indices in the block matrix
        for c in CartesianIndices(block)
            # Transform the i-index from the local index c[1] to the global index, as read of from the corresponding row_supernode
            out[col.row_supernodes[k][c[1]], c[2]] = block[c]
        end
    end
    return out
end

# We store the factorization as an LL' factorization in column-major ordering
struct SupernodalFactorization{RT}<:AbstractSupernodalSparseArray{RT,2}
    # This is the underlying buffer storing all nonzeros of the factor
    # The nonzeros appear column-by-column and, within each column, supernode by supernode

    # The (supernodal) sparse matrix structure
    data::SparseMatrixCSC{ContiguousBufferedMatrix{RT},Int}

    # The member lists of the underlying supernodes
    row_supernodes::Vector{Vector{Int}}

    # column_supernodes, stored as a vector of vectors of sparse supernodal vectors representing the hcated basis functions in terms of the row_supernodes
    # The outermost vector loops over the different colors of the multicolor ordering.
    column_supernodes::Vector{Vector{SupernodalSparseVector{RT}}}

    # The buffer underlying nnz
    buffer::Vector{RT}
end

# A function that takes in a  supernodal factorization and
# integers m, n that give the number of (block) rows and columns,
# and returns a supernodal column that can be used as scratch space when multiplying the factorizatoin with a supernodal column.
function create_scratch_space(𝐅::SupernodalFactorization{RT}, m, n) where RT
    column_lengths = size.(vcat(𝐅.column_supernodes...)[1 : m], 2)
    delims = [1; cumsum(column_lengths) .+ 1]
    column_supernodes = [Vector{Int}(1 : (delims[end] - 1))[delims[k] : (delims[k + 1] - 1)] for k = 1 : m]
    @assert all(length.(column_supernodes) == column_lengths)
return SupernodalVector(zeros(RT, delims[end] - 1, n), column_supernodes)
end

function SparseMatrixCSC(𝐅::SupernodalFactorization)
    # vertically concatenates the vectors corresponding to the different colors,
    # then horizontally concatenates the resulting matrices
    row_supernodes = 𝐅.row_supernodes
    column_lengths = size.(vcat(𝐅.column_supernodes...), 2)
    delims = [1; cumsum(column_lengths) .+ 1]
    column_supernodes = [Vector{Int}(1 : (delims[end] - 1))[delims[k] : (delims[k + 1] - 1)] for k = 1 : size(𝐅.data, 2)]

    I = Int[]
    J = Int[]
    V = eltype(eltype(𝐅.data))[]
    for (i, j, mat) in zip(findnz(𝐅.data)...)
        for cart in CartesianIndices(mat)
            (i_local, j_local) = Tuple(cart)
            push!(I, row_supernodes[i][i_local])
            push!(J, column_supernodes[j][j_local])
            push!(V, mat[i_local, j_local])
        end
    end
    return sparse(I, J, V)
end

# TODO: provide way to extract sparse matrix representative from
# supernodal factorization
# function sparse

# Constructor for the supernodal factorization
function SupernodalFactorization{RT}(I::AbstractVector{Int}, J::AbstractVector{Int}, row_supernodes::AbstractVector{<:AbstractVector{Int}}, column_supernodes::AbstractVector{<:AbstractVector{SupernodalSparseVector{RT}}}) where {RT<:Real}
    lengths_column_supernodes = size.(vcat(vcat(column_supernodes...)...), 2)
    @assert length(I) == length(J)
    M = sum(length.(row_supernodes))
    N = sum(lengths_column_supernodes)

    @assert M == N
    # Check that row and column supernode indices are valid
    @assert sort(vcat(row_supernodes...)) == 1 : M
    m = length(row_supernodes)
    n = length(lengths_column_supernodes)
    # using the existing sparse functionality to construct rowval and colptr
    data = sparse(I, J, fill(ContiguousBufferedMatrix{RT}(), length(I)))
    # Length of buffer is equal to the sum of the product of the length of the supernodes involved in each entry
    # After removind duplicates
    I, J, ~ = findnz(data)
    row_lengths = length.(row_supernodes)
    cum_sum_product_lengths = cumsum(row_lengths[I] .* lengths_column_supernodes[J])
    buffer = Vector{RT}(undef, cum_sum_product_lengths[end])
    cum_sum_product_lengths = 1 .+ vcat([0], cum_sum_product_lengths)
    for k = 1 : (length(cum_sum_product_lengths) - 1)
        data.nzval[k] = reshape(view(buffer, cum_sum_product_lengths[k] : (cum_sum_product_lengths[k + 1]) - 1), row_lengths[I[k]], lengths_column_supernodes[J[k]])
    end
    return SupernodalFactorization{RT}(data, row_supernodes, column_supernodes, buffer)
end

# Constructs a supernodal factorization from a multicolor ordering
function SupernodalFactorization(multicolor_ordering::AbstractVector{<:AbstractVector{SuperNodeBasis{PT,RT}}}, domain_supernodes::AbstractVector{<:SuperNodeDomain}, tree_function) where {PT,RT<:Real}
    row_supernodes = [id.(domains(domain_supernodes[k])) for k = 1 : length(domain_supernodes)]
    # Constructing the column supernodes
    # The ordering implied by the row_supernodes
    column_supernodes = Vector{Vector{SupernodalSparseVector{RT}}}(undef, length(multicolor_ordering))
    for k = 1 : length(column_supernodes)
        column_supernodes[k] = Vector{SupernodalSparseVector{RT}}(undef, length(multicolor_ordering[k]))
        for (l, node) in enumerate(multicolor_ordering[k])
            column_supernodes[k][l] = SupernodalSparseVector(node, row_supernodes)
        end
    end
    # Gather the lengths of the domain supernodes
    # obtain the (indices of) the row_supernodes
    I = Vector{Int}(undef, 0)
    J = Vector{Int}(undef, 0)
    # the offsets used to transform position within a given color to the position in the full sparse array
    offsets = [0; cumsum(length.(multicolor_ordering))[ 1 : (end - 1)]]
    # iterate over all colors
    for (offset, color) in zip(offsets, multicolor_ordering)
        # for each domain supernode, compute its nearest neighbor within the present color-- the basis_supernode that it will be assigned to.
        neighbors = nn(tree_function(center.(color)), center.(domain_supernodes))[1]
        # Each entry of neighbors contains the id (with respect to the array color) of the column that it is added to
        for (i, j_color) in enumerate(neighbors)
            push!(I, i)
            push!(J, offset + j_color)
        end
    end
    SupernodalFactorization{RT}(I, J, row_supernodes, column_supernodes)
end

function supernodal_size(𝐋::SupernodalFactorization)
    return (length(𝐋.row_supernodes), length(vcat(𝐋.column_supernodes...)))
end

function supernodal_size(𝐋::SupernodalFactorization, dim)
    return supernodal_size(𝐋)[dim]
end

function inner_lengths(𝐋::SupernodalFactorization)
    return length.(𝐋.row_supernodes)
end

# Multiply the matrix implied by a supernodal factorization with a supernodal vector,
function partial_multiply!(out::SupernodalVector, 𝐋::SupernodalFactorization, in::SupernodalVector; max_k=supernodal_size(𝐋, 2), scratch=create_scratch_space(𝐋, max_k, size(in, 2)))

    scratch.buffer .= 0.0
    # computing scratch = 𝐋' * in
    for k = 1 : max_k
        for i_index = getcolptr(𝐋)[k] : (getcolptr(𝐋)[k + 1] - 1)
            i = rowvals(𝐋)[i_index]
            mul!(scratch.data[k], nonzeros(𝐋)[i_index]', in.data[i], 1, 1)
        end
    end

    out.buffer .= 0.0
    # Computing out = L * scratch
    for k = 1 : max_k
        for i_index = getcolptr(𝐋)[k] : (getcolptr(𝐋)[k + 1] - 1)
            i = rowvals(𝐋)[i_index]
            mul!(out.data[i], nonzeros(𝐋)[i_index], scratch.data[k], 1, 1)
        end
    end
end

# A function that computes the Cholesky factorization of the k-"diagonal" block of 𝐅 and
# divides the corresponding column by it
function normalize_column!(𝐅::SupernodalFactorization, k)
    # The constructor of the sparse matrix does not need the buffer, so we set it to undef
    column = SparseMatrixCSC(SupernodalSparseVector(𝐅.data[:, k], Vector{eltype(𝐅.buffer)}(undef, 0), 𝐅.row_supernodes))
    D = (Matrix(SparseMatrixCSC(vcat(𝐅.column_supernodes...)[k])' * column))
    LD = cholesky(Hermitian(D)).L'
    for s in findnz(𝐅.data[:, k])[2]
        s .= s / LD
    end
end

# scatters the Supernodal Column 𝐌 into the columns of 𝐋 given by color
function scatter_column!(𝐋::SupernodalFactorization, 𝐌::SupernodalVector,color::AbstractVector{Int})
    # We check that each children supernode appears at most once in the sparsity pattern of the columns of the color
    @assert (rowvals(𝐋)[getcolptr(𝐋)[color[1]]: getcolptr(𝐋)[color[end]]]) == unique(rowvals(𝐋)[getcolptr(𝐋)[color[1]]: getcolptr(𝐋)[color[end]]])
    for index = getcolptr(𝐋)[color[1]] : (getcolptr(𝐋)[color[end] + 1] - 1)
        # Setting the nzval
        nonzeros(𝐋)[index] .= view(𝐌.data[rowvals(𝐋)[index]], :, 1 : size(nonzeros(𝐋)[index], 2))
    end
end

# Finish function that creates the measurement matrix
function create_measurement_matrix(multicolor_ordering::AbstractVector{<:AbstractVector{SuperNodeBasis{PT,RT}}}, row_supernodes) where {PT<:AbstractArray{<:Real}, RT<:Real}
    # obtaining M by searching for the dimensionality of the coefficients of the first basis function of the first supernode of the first color
    M = size(coefficients(first(basis_functions(first(first(multicolor_ordering))))), 1)
    out = SupernodalVector{RT}[]
    for color in multicolor_ordering
        # We initialize the arrays that will store the data about the sparse matrix of entries
        I = Int[]; J = Int[]; S = RT[]
        # we iterate through all supernodes of the present color
        for node in color
            # we extract the supernodes of the resulting matrix
            I_loc, J_loc, S_loc = findnz(hcat(coefficients.(basis_functions(node))...))
            append!(I, I_loc)
            append!(J, J_loc)
            append!(S, S_loc)
        end
        # we check that the sparsity patterns of the different coefficient vectors are indeed disjoint in the sense that no coordinate occurs more than once
        @assert length([(I[k], J[k]) for k = 1 : length(I)]) == length(unique([(I[k], J[k]) for k = 1 : length(I)]))
        #We have now assembled all nonzero values associated to a given color
        push!(out, SupernodalVector(Matrix(sparse(I, J, S, M, maximum(J))), row_supernodes))
    end
    return out
end

# computes a vector of supernodal vectors 𝐎 from a linear operator Θ and a measurement matrix 𝐌
function measure(Θ, 𝐌::AbstractVector{<:SupernodalVector}, row_supernodes)
    𝐎 = Vector{eltype(𝐌)}(undef, length(𝐌))
    for k = 1 : length(𝐌)
        𝐎[k] = SupernodalVector(Θ * Matrix(𝐌[k]), row_supernodes)
    end
    return 𝐎
end

# computes a vector of supernodal vectors 𝐎 from a linear operator Θ and a measurement matrix 𝐌
function measure(A::Factorization, 𝐌::AbstractVector{<:SupernodalVector}, row_supernodes)
    𝐎 = Vector{eltype(𝐌)}(undef, length(𝐌))
    for k = 1 : length(𝐌)
        𝐎[k] = SupernodalVector(A \ Matrix(𝐌[k]), row_supernodes)
    end
    return 𝐎
end

# function that uses an existing supernodal factorization and a vector of measurements to reconstruct the solver from which the measurements arose.
function reconstruct!(𝐅::SupernodalFactorization{RT}, 𝐎::Vector{<:SupernodalVector{RT}}, 𝐌::Vector{<:SupernodalVector{RT}}, multicolor_ordering) where RT<:AbstractFloat
    # TODO: replace with a better assertion? This one is not true, in general
    # @assert length(𝐎) == supernodal_size(𝐅, 2)

    # we are forming an
    colors = cumsum(length.(multicolor_ordering))
    prepend!(colors, [0])
    colors .+= 1
    colors = [colors[k] : (colors[k + 1] - 1) for k = 1 : (length(colors) - 1)]

    for k = 1 : length(𝐎)
        # create copy of vector to be used in the partial multiply
        temp = copy(𝐎[k])
        fill!(temp.buffer, 0.0)

        # doing a partial multiply up to the last color
        partial_multiply!(temp, 𝐅, 𝐌[k]; max_k=colors[k][1] - 1)
        𝐎[k].buffer .= 𝐎[k].buffer .- temp.buffer
        scatter_column!(𝐅, 𝐎[k], colors[k])
        # normalizing by diagonal square root
        # Temporarily removing the
        for l in colors[k]
            normalize_column!(𝐅, l)
        end
    end
end


text:
export Tree, PairSet, Tree!, link!,
    neighbors, children, neighbor_edges, child_edges,
    edges, vertices, root_to_leaves, leaves_to_root, both_ways

# Tree data structure

type Tree
    neighbors::Vector{Tree}
end
Tree() = Tree([])

show(io::IO, v::Tree) = print(io, "Vertex ", UInt(pointer_from_objref(v)))

# Pair sets (used for edges)

immutable PairSet{T} t::Tuple{T,T} end
PairSet{T}(t1::T, t2::T) = PairSet((t1,t2))
getindex(s::PairSet, i) = s.t[i]
=={T}(s::PairSet{T}, t::PairSet{T}) = (s[1] == t[1] && s[2] == t[2]) || (s[1] == t[2] && s[2] == t[1])
hash(s::PairSet, h::UInt) = hash(s[1],h) $ hash(s[2],h)
show(io::IO, s::PairSet) = print(io, "PairSet{",eltype(s),"}(",s[1],",",s[2],")")

# Allows to write (a,b) = PairSet(a,b)
start(s::PairSet) = start(s.t)
next(s::PairSet, state) = next(s.t, state)
done(s::PairSet, state) = done(s.t, state)
eltype{T}(s::PairSet{T}) = T
length(s::PairSet) = 2

# Tree construction

function Tree!(neighbors::Tree...)
    v = Tree([neighbors...])
    for u in neighbors
        push!(u.neighbors,v)
    end
    return v
end
Tree(neighbors::Tree...) = Tree!(deepcopy(neighbors)...)

function link!(v::Tree,u::Tree)
    push!(v.neighbors,u)
    push!(u.neighbors,v)
    return v
end

# Property functions

neighbors(v::Tree) = v.neighbors
children(v::Tree, p::Tree) = filter(u -> u != p, neighbors(v))
neighbor_edges(v::Tree) = [PairSet(u,v) for u in neighbors(v)]
child_edges(v::Tree, p::Tree) = [PairSet(u,v) for u in children(v,p)]

# Tree traversal

abstract TraversalOrder
immutable RootToLeaves <: TraversalOrder end
const root_to_leaves = RootToLeaves()
immutable LeavesToRoot <: TraversalOrder end
const leaves_to_root = LeavesToRoot()
immutable BothWays <: TraversalOrder end
const both_ways = BothWays()

function edges_rec(v,p, order::TraversalOrder)
    for u in children(v,p)
        if isa(order, RootToLeaves) produce(PairSet(u,v)); end
        if isa(order, BothWays) produce(PairSet(v,u)); end
        edges_rec(u,v, order)
        if isa(order, LeavesToRoot) produce(PairSet(u,v)); end
        if isa(order, BothWays) produce(PairSet(u,v)); end
    end
end
edges_with_root(v::Tree, p::Tree, order) = @task begin
    if isa(order, RootToLeaves) produce(PairSet(v,p)) end
    if isa(order, BothWays) produce(PairSet(p,v)); end
    edges_rec(v,p, order)
    if isa(order, LeavesToRoot) produce(PairSet(v,p)) end
    if isa(order, BothWays) produce(PairSet(v,p)); end
end
edges(v::Tree, p::Tree, order) = @task edges_rec(v,p, order)

function vertices_rec(v,p, order::TraversalOrder)
    if isa(order, RootToLeaves) produce(v); end
    if isa(order, BothWays) produce(v); end
    for u in children(v,p)
        vertices_rec(u,v, order)
    end
    if isa(order, LeavesToRoot) produce(v); end
    if isa(order, BothWays) produce(v); end
end
vertices(v::Tree, p::Tree, order) = @task vertices_rec(v,p, order)

for f in (:edges, :edges_with_root, :vertices)
    @eval begin
        $f(v::Tree, order) = $f(v,v,order)
        $f(x, order) = $f(tree(x),order)

        $f(v::Tree, p::Tree) = $f(v,p,root_to_leaves)
        $f(v::Tree) = $f(v,v)
        $f(x) = $f(tree(x))
    end
end


text:
function Gaussian(Q)
    cholesky(Q).L*randn(size(Q, 1))
end


text:
#
println("Ea) Three-step cascade after 1s-3p photo-excitation of Si^+: Configuration model.")
#
setDefaults("print summary: open", "zzz-Cascade.sum")
setDefaults("method: continuum, asymptotic Coulomb")    ## setDefaults("method: continuum, Galerkin")
setDefaults("method: normalization, pure sine")         ## setDefaults("method: normalization, pure Coulomb")    setDefaults("method: normalization, pure sine")

wa = Cascade.Computation("Si+ (1s^-1) cascade", Nuclear.Model(18.),
                         JAC.Radial.Grid("grid: by given parameters"; rnt = 2.0e-5, h = 5.0e-2, hp = 2.0e-2, NoPoints = 600), AsfSettings(),
                         JAC.Cascade.AverageSCA(), [JAC.Radiative],          #  JAC.Cascade.AverageSCA()
                         [Configuration("1s^1 2s^2 2p^6 3s^2 3p^2")], [(1, 1.)],  0, 0, Shell[], Shell[], JAC.Cascade.Step[])
wb = perform(wa; output=true)
setDefaults("print summary: close", "")

## using JLD
## save("zzz-Cascade.jld", wb)


text:
## Read metadata of all food webs archived on mangal.io (generated from 01_import_mangal_metadata.jl)
mangal_foodwebs = DataFrame(CSV.File(joinpath("data", "raw", "mangal_foodwebs.csv")))
N = network.(mangal_foodwebs.id)
N = convert.(UnipartiteNetwork, N)

## Define flexible links model and infer model parameters

# S: number of species
# R: number of flexible links that are realized (links above the minimum)

@model FL(S,R) = begin
  # number of trials
  N = length(S)
  F = S.*S .- (S.-1)
  # parameters (priors)
  ϕ ~ Normal(3.0, 0.5)
  μ ~ Beta(3.0, 7.0)
  # flexible links model
  for i in 1:N
    R[i] ~ BetaBinomial(F[i], μ*exp(ϕ), (1-μ)*exp(ϕ))
  end
  return μ, ϕ
end

# Observations

# number of species
S = vec(mangal_foodwebs[:,:S])
# number of links
L = vec(mangal_foodwebs[:,:L])
# number of flexible links that are realized
R = L .- (S.-1)

# Hamiltonian Monte Carlo sampler with static trajectory
chain = sample(FL(S,R), HMC(0.01,10), 3000)

# Diagnostic plot
plot(chain[200:end,:,:])

## Simulate counterfactuals

# Simulate the number of links for a range of species richness

# numbers of species
sp = collect(5:1000)
# number of simulations
n = 1000

predicted_links = zeros(Int64, n, length(sp))

p = Progress(length(sp))
Threads.@threads for s in sp
    i = findall(sp .== s)
    predicted_links[:, i] = predict_links(s, n, chain)
    next!(p)
end

save(joinpath("data", "sim", "predicted_links.jld"), "data", predicted_links)


text:
import LinearAlgebra: Givens, givensAlgorithm
import LinearAlgebra.ldiv!

struct FastHessenberg{T<:AbstractMatrix}
    H::T # H is assumed to be Hessenberg of size (m + 1) × m
end

@inline Base.size(H::FastHessenberg, args...) = size(H.H, args...)

"""
Solve Hy = rhs for a non-square Hessenberg matrix.
Note that `H` is also modified as is it converted
to an upper triangular matrix via Given's rotations
"""
function ldiv!(H::FastHessenberg, rhs)
    # Implicitly computes H = QR via Given's rotations
    # and then computes the least-squares solution y to
    # |Hy - rhs| = |QRy - rhs| = |Ry - Q'rhs|

    width = size(H, 2)

    # Hessenberg -> UpperTriangular; also apply to r.h.s.
    @inbounds for i = 1 : width
        c, s, _ = givensAlgorithm(H.H[i, i], H.H[i + 1, i])

        # Skip the first sub-diagonal since it'll be zero by design.
        H.H[i, i] = c * H.H[i, i] + s * H.H[i + 1, i]

        # Remaining columns
        @inbounds for j = i + 1 : width
            tmp = -conj(s) * H.H[i, j] + c * H.H[i + 1, j]
            H.H[i, j] = c * H.H[i, j] + s * H.H[i + 1, j]
            H.H[i + 1, j] = tmp
        end

        # Right hand side
        tmp = -conj(s) * rhs[i] + c * rhs[i + 1]
        rhs[i] = c * rhs[i] + s * rhs[i + 1]
        rhs[i + 1] = tmp
    end

    # Solve the upper triangular problem.
    U = UpperTriangular(view(H.H, 1 : width, 1 : width))
    ldiv!(U, view(rhs, 1 : width))
    nothing
end


text:
# This file is a part of AstroLib.jl. License is MIT "Expat".
# Copyright (C) 2016 Mosè Giordano.

# Test adstring
@test adstring((30.4, -1.23), truncate=true) ==
    adstring([30.4, -1.23], truncate=true) == " 02 01 35.9  -01 13 48"
@test adstring(19.19321, truncate=true) == "+19 11 35.5"
@test adstring(ten(36,24,15.015), -ten(8,24,36.0428), precision=3) ==
    " 02 25 37.0010  -08 24 36.043"
@test adstring([30.4, -15.63], [-1.23, 48.41], precision=2) ==
    [" 02 01 36.000  -01 13 48.00", " 22 57 28.800  +48 24 36.00"]
@test adstring([(-58, 724)]) == [" 20 08 00.0  +724 00 00"]

# Test airtovac
@test_approx_eq airtovac([1234 6056.125]) [1234.0 6057.801930991426]
@test_approx_eq airtovac(2100) 2100.666421596007

# Test aitoff
@test aitoff([227.23, 130], [-8.890, -35]) ==
    ([-137.92196683723276, 115.17541338020645], [-11.772527357473054, -44.491889962090085])
@test aitoff([375], [2.437]) ==
    ([16.63760711611838],[2.712427279646118])
@test aitoff((227.23, -8.890)) ==
    (-137.92196683723276,-11.772527357473054)

# Test altaz2hadec
@test altaz2hadec(59.086111, 133.30806, 43.07833) ==
    altaz2hadec((59.086111, 133.30806), 43.07833) ==
    (336.68286017949157, 19.182449588316555)
@test altaz2hadec([15, 25, 35],
                  [25.12, 45.32, -20.3],
                  [-23.44, 45.0, 52.5]) ==
                      ([324.9881067314537, 256.7468302330436, 132.4919217875949],
                       [44.38225395397647, 48.542947077386664, 67.33061196497327])

# Test bprecess
let
    local ra, dec
    ra, dec = bprecess([ten(13, 42, 12.74)*15], [ten(8, 23, 17.69)],
                       reshape(100*[-15*0.0257, -0.090], 2, 1))
    @test_approx_eq ra  204.93552515632123
    @test_approx_eq dec 8.641287183886163
    ra, dec = bprecess(82, 19)
    @test_approx_eq ra  81.26467916346334
    @test_approx_eq dec 18.959495700195394
    ra, dec = bprecess([57], [23], 2024)
    @test_approx_eq ra  56.26105898810067
    @test_approx_eq dec 22.84693298145991
    ra, dec = bprecess([57], [23], reshape([9, 86], 2, 1), parallax=[1], radvel=[4])
    @test_approx_eq ra  56.25988479854577
    @test_approx_eq dec 22.83493370392355
    ra, dec = bprecess((-57, -23), 2024)
    @test_approx_eq ra   302.2593299643789
    @test_approx_eq dec -23.150089972802036
    ra, dec = bprecess((-57, -23), [9, 86], parallax=1, radvel=4)
    @test_approx_eq ra   302.2580376402947
    @test_approx_eq dec -23.16208183899836
end

# Test calz_unred
@test_approx_eq calz_unred(reshape(900:1000:9900, 10, 1), ones(Float64, 10), -0.1) [1.0,                0.43189326452379095,
                                                                                0.5203675483533704, 0.594996469192435,
                                                                                0.6569506252451913, 0.7080829505773865,
                                                                                0.7502392743978797, 0.7861262388745882,
                                                                                0.8151258710444882, 0.8390325371659836]

#test ct2lst
@test_approx_eq ct2lst(-76.72, -4, [DateTime(2008, 7, 30, 15, 53)]) [11.356505172312609]
@test_approx_eq ct2lst(9, [jdcnv(2015, 11, 24, 12, 21)]) [17.159574059885927]

# Test daycnv with Gregorian Calendar in force.
@test daycnv(2440000.0) == DateTime(1968, 05, 23, 12)
# Test daycnv with Julian Calendar in force (same result as IDL AstroLib's
# daycnv).
@test daycnv(2000000.0) == DateTime(763, 09, 18, 12)
@test daycnv(0.0) == DateTime(-4713, 11, 24, 12)

# Test deredd
let
    local by0, m0, c0, ub0
    by0, m0, c0, ub0 = deredd([0.5, -0.5], [0.2, 0.5], [1, 1], [1, 1], [0.1, 0.3])
    @test_approx_eq by0 [-0.3,0.5]
    @test_approx_eq m0  [1.165,1.0]
    @test_approx_eq c0  [0.905,1.0]
    @test_approx_eq ub0 [-0.665,0.3]
end

# Test eci2geo
let
    local lat, long, alt, jd
    lat, long, alt = eci2geo([0], [0], [0], [2452343])
    @test_approx_eq lat  [0]
    @test_approx_eq long [12.992783145436988]
    @test_approx_eq alt  [-6378.137]
    lat, long, alt = eci2geo((6978.137, 0, 0), jdcnv("2015-06-30T14:03:12.857"))
    @test_approx_eq lat  0
    @test_approx_eq long 230.87301833205856
    @test_approx_eq alt  600
    # Test `eci2geo' is the inverse of `geo2eci'
    jd = get_juldate()
    lat, long, alt = eci2geo(geo2eci(10, 10, 10, jd), jd)
    @test_approx_eq lat  10
    @test_approx_eq long 10
    @test_approx_eq alt  10
end

# Test eqpole
let
    local x, y
    x, y = eqpole([100], [35], southpole=true)
    @test_approx_eq x [-111.18287262822456]
    @test_approx_eq y [ -19.604540237028665]
    x, y = eqpole([80], [19])
    @test_approx_eq x [72.78853915267848]
    @test_approx_eq y [12.83458333897169]
end

# Test flux2mag
@test_approx_eq flux2mag([1.5e-12, 8.7e-15, 4.4e-10]) [8.459771852360795,
                                                       14.051201868453454,
                                                       2.291368308784527]
@test_approx_eq flux2mag(1) -21.1
@test_approx_eq flux2mag(5.2e-15) 14.609991640913002
@test_approx_eq flux2mag(5.2e-15, 15) 20.709991640913003
@test_approx_eq flux2mag(5.2e-15, ABwave=15) 27.423535345634598

# Test gal_uvw
let
    local u, v, w
    u, v, w = gal_uvw([ten(1,9,42.3)*15], [ten(61,32,49.5)], [627.89], [77.84],
                      [-321.4], [1e3/129], lsr=true)
    @test_approx_eq u [118.2110474553902]
    @test_approx_eq v [-466.4828898385057]
    @test_approx_eq w [88.16573278565097]
    u, v, w = gal_uvw(1, 2, 3, 4, 5, 6)
    @test_approx_eq u  4.0228405867158745
    @test_approx_eq v  3.7912174342038227
    @test_approx_eq w -3.1700191400725464
end

# Test geo2eci
let
    local x, y, z
    x, y, z = geo2eci([0], [0], [0], [2452343])
    @test_approx_eq x [6214.846433007192]
    @test_approx_eq y [-1433.9858454345972]
    @test_approx_eq z [0.0]
    x, y, z = geo2eci((0,0,0), jdcnv("2015-06-30T14:03:12.857"))
    @test_approx_eq x -4024.8671780315185
    @test_approx_eq y 4947.835465127513
    @test_approx_eq z 0.0
end

# Test geo2geodetic
let
    local lat, long, alt
    lat, long, alt = geo2geodetic([90], [0], [0], "Jupiter")
    @test_approx_eq lat [90]
    @test_approx_eq long [0]
    @test_approx_eq alt 4638
    lat, long, alt = geo2geodetic((90, 0, 0))
    @test_approx_eq lat 90
    @test_approx_eq long 0
    @test_approx_eq alt 21.38499999999931
    lat, long, alt = geo2geodetic((43.16, -24.32, 3.87), 8724.32, 8619.19)
    @test_approx_eq lat 43.849399515234516
    @test_approx_eq long -24.32
    @test_approx_eq alt 53.53354478670836
    lat, long, alt = geo2geodetic([43.16], [-24.32], [3.87], 8724.32, 8619.19)
    @test_approx_eq lat [43.849399515234516]
    @test_approx_eq long [-24.32]
    @test_approx_eq alt [53.53354478670836]
end

# Test geo2mag
let
    local lat, long
    lat, long = geo2mag(ten(35,0,42), ten(135,46,6), 2016)
    @test_approx_eq lat 36.86579228937769
    @test_approx_eq long -60.184060536651614
    lat, long = geo2mag([15], [24], 2016)
    @test_approx_eq lat [11.452100529696096]
    @test_approx_eq long [-169.86030510727102]
end

# Test geodetic2geo
let
    local lat, long, alt
    lat, long, alt = geodetic2geo([90], [0], [0], "Jupiter")
    @test_approx_eq lat [90]
    @test long == [0]
    @test_approx_eq alt -4638
    lat, long, alt = geodetic2geo((90, 0, 0))
    @test_approx_eq lat 90
    @test long == 0
    @test_approx_eq alt -21.38499999999931
    lat, long, alt = geodetic2geo((43.16, -24.32, 3.87), 8724.32, 8619.19)
    @test_approx_eq lat 42.46772711708433
    @test long == -24.32
    @test_approx_eq alt -44.52902080669082
    lat, long, alt = geodetic2geo([43.16], [-24.32], [3.87], 8724.32, 8619.19)
    @test_approx_eq lat [42.46772711708433]
    @test long == [-24.32]
    @test_approx_eq alt [-44.52902080669082]
    # Test geodetic2geo is the inverse of geo2geodetic, within a certain
    # tolerance.
    lat, long, alt = geodetic2geo(geo2geodetic(67.2,13.4,1.2))
    @test_approx_eq_eps lat 67.2 1e-8
    @test long == 13.4
    @test_approx_eq_eps alt 1.2 1e-9
end

# Test get_date with mixed keywords.
@test get_date(DateTime(2001,09,25,14,56,14), old=true,timetag=true) ==
    get_date(2001,09,25,14,56,14, old=true,timetag=true) ==
    get_date("2001-09-25T14:56:14", old=true,timetag=true) ==
    "25/09/2001T14:56:14"
@test get_date(DateTime(2001,09,25,14,56,14)) ==
    get_date(2001,09,25,14,56,14) == get_date("2001-09-25T14:56:14") ==
    "2001-09-25"
@test get_date([DateTime(2024), Date(2016, 3, 14)]) ==
    get_date([Date(2024), "2016-03-14"]) ==
    get_date(["2024-01", DateTime(2016, 3, 14)]) == ["2024-01-01", "2016-03-14"]

# Test gcirc.
@test_approx_eq gcirc(0, [0,1,2], [1,2,3], [2,3,4], [3,4,5]) [1.222450611061632,
                                                              2.500353926443337,
                                                              1.5892569925227757]
@test_approx_eq gcirc(0,  120, -43,   175, +22)  1.590442261600714
@test_approx_eq gcirc(1, (120, -43),  175, +22)  415908.56615322345
@test_approx_eq gcirc(2,  120, -43,  (175, +22)) 296389.3666794745
@test_approx_eq gcirc(0, (120, -43), (175, +22)) 1.590442261600714
@test_approx_eq gcirc(1, [120], [-43],  175, +22)  [415908.56615322345]
@test_approx_eq gcirc(2,  120, -43,  [175], [+22]) [296389.3666794745]
@test_throws ErrorException gcirc(3, 0, 0, 0, 0)

# Test hadec2altaz
let
    local alt1, az1, alt2, az2
    alt1, az1 = hadec2altaz([0], [11.978165], [ten(43,4,42)])
    @test_approx_eq alt1 [58.89983166666667]
    @test_approx_eq az1 [180.0]
    @test_approx_eq hadec2altaz((0, 11.978165), ten(43,4,42), ws=true)[2] 0.0
    alt1, az1 = 50, 20
    alt2, az2 = hadec2altaz(altaz2hadec(alt1, az1, 40), 40)
    @test_approx_eq alt1 alt2
    @test_approx_eq az1  az2
end

# Test helio_jd
@test_approx_eq helio_jd([juldate(2016, 6, 15, 11, 40)], ten(20, 9, 7.8)*15, ten(37, 9, 7)) 57554.98808289718
@test_approx_eq helio_jd(1000, 23, 67, B1950=true) 999.9997659545342
@test_approx_eq helio_jd(2000, 12, 88, diff=true) -167.24845957792076

# Test helio_rv
@test_approx_eq helio_rv(helio_jd(juldate(94, 10, 25, 17, 30),
                                  ten(04, 38, 16)*15, ten(20, 41, 05)),
                         46487.5303, 2.0563056, -6, 59.3) -62.965570109145034
@test_approx_eq helio_rv([0.1, 0.9], 0, 1, 0, 100, 0.6, 45) [-45.64994926111004,
                                                             89.7820347358485]

# Test jdcnv.
@test_approx_eq jdcnv(-4713, 11, 24, 12) 0.0
@test jdcnv(763, 09, 18, 12) == jdcnv("763-09-18T12") == 2000000.0
@test (jd=1234567.89; jdcnv(daycnv(jd)) == jd)
@test jdcnv([DateTime(2016, 07, 31), "1969-07-20"]) ==
    jdcnv([Date(2016, 07, 31), DateTime(1969, 07, 20)]) ==
    jdcnv(["2016-07-31", Date(1969, 07, 20)])

# Test jprecess
let
    local ra, dec
    ra, dec = jprecess([ten(13, 39, 44.526)*15], [ten(8, 38, 28.63)],
                       reshape(100*[-15*0.0259, -0.093], 2, 1))
    @test_approx_eq ra  205.5530845731372
    @test_approx_eq dec 8.388247441904628
    ra, dec = jprecess(82, 19)
    @test_approx_eq ra  82.73568745151148
    @test_approx_eq dec 19.036972917272056
    ra, dec = jprecess([57], [23], 2024)
    @test_approx_eq ra  57.74049975335702
    @test_approx_eq dec 23.150053754297726
    ra, dec = jprecess([57], [23], reshape([9, 86], 2, 1), parallax=[1], radvel=[4])
    @test_approx_eq ra  57.74180294549785
    @test_approx_eq dec 23.16200582079095
    ra, dec = jprecess((-57, -23), 2024)
    @test_approx_eq ra  303.73910971499015
    @test_approx_eq dec -22.846895476784482
    ra, dec = jprecess((-57, -23), [9, 86], parallax=1, radvel=4)
    @test_approx_eq ra  303.7402950607101
    @test_approx_eq dec -22.834931625610313
end

# Test juldate with Gregorian Calendar in force.  This also makes sure precision
# of the result is high enough.  Note that "juldate(dt::DateTime) =
# Dates.datetime2julian(dt)-2.4e6" would not be precise.
@test_approx_eq juldate([DateTime(2016, 1, 1, 8)]) (57388.5 + 1.0/3.0)

# Test juldate with Julian Calendar in force, for different centuries.  This
# also makes sure precision of the result is high enough.
@test_approx_eq juldate(1582, 10, 1, 20) (-100843.0 + 1.0/3.0)
@test_approx_eq juldate(["1000-01-01T20"]) (-313692.0 + 1.0/3.0)
@test_approx_eq juldate("100-10-25T20") (-642119.0 + 1.0/3.0)
@test_approx_eq juldate(-4713, 1, 1, 12) -2.4e6
# Test daycnv and juldate together, with Gregorian Calendar in force.  Note that
# they are not expected to be one the inverse of the other during Julian
# Calendar.
@test (dt=DateTime(2016, 1, 1, 20);
       daycnv(juldate(dt) + 2.4e6) == dt)

# Test kepler_solver
@test_approx_eq kepler_solver(8pi/3, 0.7) 2.5085279492864223
@test_approx_eq kepler_solver([pi/4, pi/6, 8pi/3], 0) [pi/4, pi/6, 2pi/3]
@test_approx_eq kepler_solver(3pi/2, 0.8) -2.2119306096084457
@test_approx_eq kepler_solver(0, 1) 0.0
@test_throws ErrorException kepler_solver(pi, -0.5)
@test_throws ErrorException kepler_solver(pi,  1.5)

# Test mag2flux
let
    local vel, lsf
    vel, lsf = lsf_rotate(3, 90)
    @test length(vel) == length(lsf) == 61
    vel, lsf = lsf_rotate(5, 10)
    @test_approx_eq vel collect(-10.0:5.0:10.0)
    @test_approx_eq lsf [0.0, 0.556914447710896, 0.6933098861837907, 0.556914447710896, 0.0]
end

# Test mag2flux
@test_approx_eq mag2flux(4.83, 21.12) 4.1686938347033296e-11
@test_approx_eq mag2flux([4.83], 21.12) 4.1686938347033296e-11
@test_approx_eq flux2mag(mag2flux(15, ABwave=12.), ABwave=12) 15.0
@test_approx_eq mag2flux(8.3) 1.7378008287493692e-12
@test_approx_eq mag2flux(8.3, 12) 7.58577575029182e-9
@test_approx_eq mag2flux(8.3, ABwave=12) 3.6244115683017193e-7

# Test mag2geo
let
    local lat, long
    lat, long = mag2geo(90, 0, 2016)
    @test_approx_eq lat 86.395
    @test_approx_eq long -166.29000000000002
    lat, long = mag2geo([15], [24], 2016)
    @test_approx_eq lat [11.702066965890157]
    @test_approx_eq long [-142.6357492442842]
    # Test geo2mag is approximately the inverse of mag2geo
    lat, long = geo2mag(mag2geo(12.34, 56.78, 2016)..., 2016)
    @test_approx_eq lat 12.34
    @test_approx_eq long 56.78
end

# Test month_conv
@test month_cnv([" januavv  ", "SEPPES ", " aUgUsT", "la"]) == [1, 9, 8, -1]
@test month_cnv([2, 12, 6], short=true, low=true) == ["feb", "dec", "jun"]
@test month_cnv(5, up=true) == "MAY"
@test (list=[1, 2, 3]; month_cnv(month_cnv(list)) == list)
@test (list=["July", "March", "November"]; month_cnv(month_cnv(list)) == list)

# Test moonpos
let
    local ra, dec, dis, lng, lat
    ra, dec, dis, lng, lat =
        moonpos(jdcnv(1992, 4, 12))
    @test_approx_eq ra  134.68846854844108
    @test_approx_eq dec 13.768366630560255
    @test_approx_eq dis 368409.68481612665
    @test_approx_eq lng 133.16726428105378
    @test_approx_eq lat -3.2291264192144356
    ra, dec, dis, lng, lat =
        moonpos([2457521], radians=true)
    @test_approx_eq ra  2.2587950290926178
    @test_approx_eq dec 0.26183388011392217
    @test_approx_eq dis 385634.68772395694
    @test_approx_eq lng 2.232459255739553
    @test_approx_eq lat -0.059294466326164315
end

# Test mphase
@test_approx_eq mphase([2457520, 2457530, 2457650]) [0.2781695910737857,
                                                     0.9969808583803166,
                                                     0.9580708477591693]

# Test nutate
let
    local long, obl
    long, obl = nutate(jdcnv(1987, 4, 10))
    @test_approx_eq long -3.787931077110755
    @test_approx_eq obl   9.442520698644401
    long, obl = nutate(2457521)
    @test_approx_eq long -4.401443629818089
    @test_approx_eq obl  -9.26823431959121
    long, obl = nutate([2457000, 2458000])
    @test_approx_eq long[1]  4.327189321653877
    @test_approx_eq obl[1]  -9.507794266102866
    @test_approx_eq long[2] -9.686089990639474
    @test_approx_eq obl[2]  -6.970768250588256
end

# Test paczynski
@test_approx_eq paczynski(-1e-10)  -1e10
@test_approx_eq paczynski(1e-1)    10.037461005722337
@test_approx_eq paczynski(-1)      -1.3416407864998738
@test_approx_eq paczynski(10)       1.0001922892047386
@test_approx_eq paczynski(-1e10)   -1

# Test planck_freq
@test_approx_eq planck_freq([2000], [5000]) [6.1447146126144004e-30]

# Test planck_wave
@test_approx_eq planck_wave([2000], [5000]) [8.127064833530511e-24]

# Test polrec
let
    local x=zeros(Float64, 3), y=zeros(Float64, 3)
    x, y = polrec([1, 2, 3], [pi, pi/2.0, pi/4.0])
    @test_approx_eq x [-1.0, 0.0, 1.5*sqrt(2.0)]
    @test_approx_eq y [ 0.0, 2.0, 1.5*sqrt(2.0)]
    x, y = polrec((2, 135), degrees=true)
    @test_approx_eq x -sqrt(2)
    @test_approx_eq y  sqrt(2)
end

# Test posang.
@test_approx_eq posang(1, ten(13, 25, 13.5), ten(54, 59, 17),
                       ten(13, 23, 55.5), ten(54, 55, 31)) -108.46011246802047
@test_approx_eq posang(0, [0,1,2], [1,2,3], [2,3,4], [3,4,5]) [1.27896824717634,
                                                               1.6840484573313608,
                                                               0.2609280020139511]
@test_approx_eq posang(0,  120, -43,   175, +22)  -1.5842896165356724
@test_approx_eq posang(1, (120, -43),  175, +22)  82.97831348792039
@test_approx_eq posang(2,  120, -43,  (175, +22)) 50.02816530382374
@test_approx_eq posang(0, (120, -43), (175, +22)) -1.5842896165356724
@test_approx_eq posang(1, [120], [-43],  175, +22)  [82.97831348792039]
@test_approx_eq posang(2,  120, -43,  [175], [+22]) [50.02816530382374]
@test_throws ErrorException posang(3, 0, 0, 0, 0)

# Test precess
let
    local ra1, dec1, ra2, dec2
    ra1, dec1 = precess((ten(2,31,46.3)*15, ten(89,15,50.6)), 2000, 1985)
    @test_approx_eq ra1  34.09470328718033
    @test_approx_eq dec1 89.19647174928589
    ra2, dec2 = precess([ten(21, 59, 33.053)*15], [ten(-56, 59, 33.053)], 1950, 1975, FK4=true)
    @test_approx_eq ra2  330.3144305418865
    @test_approx_eq dec2 -56.87186126487889
end

# Test precess_xyz
let
    local x1 ,y1, z1, x2, y2, z2
    x1, y1, z1 = precess_xyz((1.2, 2.3, 1.7), 2000, 2050)
    @test_approx_eq x1 1.165933061423247
    @test_approx_eq y1 2.313228746401996
    @test_approx_eq z1 1.7057470102860104
    x2, y2, z2 = precess_xyz([0.7, -2.4], [3.3, 6.6], [0, 4], 2000, 2016)
    @test_approx_eq x2 [0.688187142071843,   -2.429815562246262]
    @test_approx_eq y2 [3.3024835038223532,   6.591359330834213]
    @test_approx_eq z2 [0.001079105285993004, 3.9962455511755794]
end

# Test premat
@test_approx_eq premat(1967, 1982, FK4=true) [0.9999933170034135    -0.0033529069683496567 -0.0014573823699636742
                                              0.00335290696825777    0.9999943789886484    -2.443304965138481e-6
                                              0.0014573823701750721 -2.4431788671274868e-6  0.9999989380147651   ]

@test_approx_eq premat(1995, 2003) [ 0.9999980977132219    -0.0017889257711428855  -0.0007773766929507687
                                     0.0017889257711354528  0.9999983998707707     -6.953448226403318e-7
                                     0.0007773766929678732 -6.953257000046125e-7    0.9999996978424512  ]

# Test radec
@test radec(15.90, -0.85) == (1.0, 3.0, 36.0, -0.0, 51.0, 0.0)
@test radec(-0.85,15.9) == (23.0,56.0,36.0,15.0,54.0,0.0)
@test radec(-20,4,hours=true) == (4.0,0.0,0.0,4.0,0.0,0.0)
@test radec([15.90, -0.85], [-0.85,15.9]) ==
    ([1.0, 23.0], [3.0, 56.0], [36.0, 36.0],
     [-0.0, 15.0], [51.0, 54.0], [0.0, 0.0])

# Test recpol
let
    local r = a = zeros(Float64, 3)
    r, a = recpol([0, sqrt(2.0), 2.0*sqrt(3.0)],
                  [0, sqrt(2.0), 2.0])
    @test_approx_eq r [0.0,  2.0,  4.0]
    @test_approx_eq a [0.0, pi/4.0, pi/6.0]
    r, a = recpol(1, 1)
    @test_approx_eq r sqrt(2.0)
    @test_approx_eq a pi/4.0
    # Test polrec is the inverse of recpol
    local xi = 6.3, yi = -2.7, x = y = 0.0
    x, y = polrec(recpol((xi, yi), degrees=true), degrees=true)
    @test_approx_eq x xi
    @test_approx_eq y yi
end

# Test rhotheta
let
    local ρ, θ
    ρ, θ = rhotheta(41.623, 1934.008, 0.2763, 0.907, 59.025, 23.717, 219.907, 1980)
    @test_approx_eq ρ 0.41101776646245836
    @test_approx_eq θ 318.4242564860495
end

# Test "sixty".  Test also it's the reverse of ten.
@test_approx_eq sixty(-51.36) [-51.0, 21.0, 36.0]
@test_approx_eq ten(sixty(-0.10934835545824395)) -0.10934835545824395
@test_approx_eq sixty(1) [1.0, 0.0, 0.0]

# Test sphdist.
@test_approx_eq sphdist([0,1,2], [1,2,3], [2,3,4], [3,4,5]) [1.222450611061632,
                                                             2.500353926443337,
                                                             1.5892569925227762]
@test_approx_eq sphdist(120, -43, 175, +22) 1.5904422616007134
@test_approx_eq sphdist([120], [-43], 175, +22) 1.5904422616007134
@test_approx_eq sphdist(120, -43, [175], [+22]) 1.5904422616007134

# Test sunpos
let
    local ra, dec, lon, obl
    ra, dec, lon, obl = sunpos(jdcnv(1982, 5, 1))
    @test_approx_eq ra  37.88589057369026
    @test_approx_eq dec 14.909699471099517
    @test_approx_eq lon 2309.6312912217463
    @test_approx_eq obl 1343.0612563977593
    ra, dec, lon, obl = sunpos(jdcnv([DateTime(2016, 5, 10)]), radians=true)
    @test_approx_eq ra  0.8259691339090751
    @test_approx_eq dec 0.3085047454107549
    @test_approx_eq lon 49.77773359512005
    @test_approx_eq obl 23.434647165246304
    ra, dec, lon, obl = sunpos([2457531])
    @test_approx_eq ra  59.71655864208797
    @test_approx_eq dec 20.52127006818727
    @test_approx_eq lon 3542.279299068626
    @test_approx_eq obl 1342.7064622183311
end

# Test "ten" and "tenv".  Always make sure string and numerical inputs are
# consistent (IDL implementation of "ten" is not).
@test ten(0, -23, 34) == ten((0, -23, 34)) == ten([0, -23, 34]) ==
    ten(" : 0 -23 :: 34") == -0.37388888888888894
@test ten(-0.0, 60) == ten((-0.0, 60)) == ten([-0.0, 60]) ==
    ten("-0.0 60") == -1.0
@test ten(-5, -60, -3600) == ten((-5, -60, -3600)) ==
    ten([-5, -60, -3600]) == ten("  -5: :-60: -3600") == -3.0
@test ten("") == 0.0
@test tenv([0, -0.0, -5], [-23, 60, -60], [34, 0, -3600]) ==
    tenv([(0, -23,34), ":-0.0:60", (-5, -60, -3600)]) ==
    tenv(["0   -23 :: 34", (-0.0, 60), " -5:-60: -3600"]) ==
    [-0.37388888888888894, -1.0, -3.0]
@test tenv([12.0, -0.0], [24, 30]) == tenv([" 12::24", " -0:30: "]) == [12.4, -0.5]

# Test kepler_solver
@test_approx_eq trueanom(8pi/3, 0.7) 2.6657104039293764
@test_approx_eq trueanom([pi/4, pi/6, 8pi/3], 0) [pi/4, pi/6, 2pi/3]
@test_approx_eq trueanom(3pi/2, 0.8) -2.498091544796509
@test_approx_eq trueanom(0.1, 1) pi
@test_throws ErrorException trueanom(pi, -0.5)
@test_throws DomainError trueanom(pi,  1.5)

# Test vactoair and that airtovac is its inverse (it isn't true only around
# 2000, just avoid those values)
@test_approx_eq vactoair([2000]) [1999.3526230448367]
@test_approx_eq airtovac(vactoair(collect(1000:300:4000))) collect(1000:300:4000)

# Test xyz
let
    local x, y, z, vx, vy, vz
    x, y, z, vx, vy, vz = xyz([51200.5 + 64./86400.], 2000)
    @test_approx_eq x  [0.5145687092402946]
    @test_approx_eq y  [-0.7696326261820777]
    @test_approx_eq z  [-0.33376880143026394]
    @test_approx_eq vx [0.014947267514081075]
    @test_approx_eq vy [0.008314838205475709]
    @test_approx_eq vz [0.003606857607574784]
end

# Test ydn2md.
@test ydn2md(2016, [60, 234]) == [Date(2016, 02, 29), Date(2016, 08, 21)]
@test ymd2dn(ydn2md(2016, 60)) == 60

# Test ymd2dn
@test ymd2dn([Date(2015,3,5), Date(2016,3,5)]) == [64, 65]
@test ydn2md(2016, ymd2dn(Date(2016, 09, 16))) == Date(2016, 09, 16)


text:
export @parameterize

using Base.Meta

macro parameterize(params::Union{Expr, Symbol}, test::Expr)
  const name = test.args[1].args[1]
  const args = test.args[1].args[2:end]
  const body = test.args[2]
  const func_expr = Expr(:function, Expr(:tuple, args...), body)

  quote
    const $(esc(name)) = Parameterize($(esc(func_expr)), $(esc(params)))
  end
end

immutable Parameterize
  fn::Function
  params::Any
end

function push_test!(tests::Vector{Tuple{AbstractString, Function}}, name::AbstractString, test::Parameterize)
  for params in test.params
    push_test!(tests, "$(name)[$params]", ()->test.fn(params...))
  end
end