text:
"""
This part of the code primarily dispatches on the `BoundedWignerCache` type on the first
argument of the internal functions. It reduces the number of allocations in heavy workloads
by more than an order of magnitude with a number of different optimizations, compared to the
code in main code file. This makes it accelerate computations with parallelism like
threading without being limited by the memory bandwidth.

* The number n! has a prime factorization containing all prime numbers less than n. Then
  the total number of prime factors one needs to store is bounded by the prime counting
  function π(n). We thus preallocate prime factorization buffers of size ~ π(n).
* Instead of modifying array sizes with push/pop, one keeps track of the largest nonzero
  prime factor during additions. Prime factorizations and sums are preallocated with bounds
  given in Johansson and Forssén 2016.
* BigInt is basically an array allocated on the heap, and thus generates a lot of
  painful allocations. Each thread reuses the same BigInts for different parts of the
  calculation, when possible.
"""

"""
The nj symbols can be written in the form r ⋅ √s, and we store such numbers in the result
caches `Wigner3j` and `Wigner6j`. Previously this was a `Tuple{T,T}`, but this was
challenging to use as a Number type. The default result key is `RSPair{Rational{BigInt}}``.
"""
struct RSPair{T<:Number} <: Number
    r::T
    s::T
end

"""
Convert a number of the form r ⋅ √s to some kind of floating point number. This needs to
be done by way of BigFloat to avoid overflows.
"""
function Base.convert(T::Type{<:AbstractFloat}, rs::RSPair)
    convert(T,  # convert to BigFloat first for floating point types
        convert(BigFloat, rs.s)*convert(BigFloat, signedroot(rs.r)))
end

"""
Convert a number of the form r ⋅ √s to an arbitrary number type involves just performing
the conversion and then the actual computation.
"""
Base.convert(T::Type{<:Real}, rs::RSPair) = convert(T, rs.s)*convert(T, signedroot(rs.r))

"""
The arbitrary-precision types in Julia are still unpolished, and they are unlike normal
numerical types as they mutable. This BigNumber type allows us to dispatch based on
mutability and make a deepcopy if it's a mutable type.
"""
const BigNumber = Union{Rational{BigInt},BigInt,BigFloat}

"""
We're trying to reuse BigInts when possible, but this means that if one wants to keep the
exact solutions in the dictionary, one must deepcopy the number.
"""
store_result(rs::RSPair{T}) where {T <: BigNumber} = deepcopy(rs)
store_result(rs) = rs

"""
The default global thread-safe dictionaries used to cache results of calculations.
They map the symmetrized form to r,s pairs.
"""
const Wigner3j = ThreadSafeDict{Tuple{UInt,UInt,UInt,Int,Int},RSPair{Rational{BigInt}}}()
const Wigner6j = ThreadSafeDict{NTuple{6,UInt},RSPair{Rational{BigInt}}}()

"""
Convenience function for generating thread-safe dictionaries which store the results.
Returns the dictionaries for Wigner3j and Wigner6j in a tuple.
"""
function wigner_dicts(Tdict::Type{<:Real})
    return (ThreadSafeDict{Tuple{UInt,UInt,UInt,Int,Int},Tdict}(),  # 3j
            ThreadSafeDict{NTuple{6,UInt},Tdict}())                 # 6j
end

"""
Convenience function for generating an array of caches for each thread. The default
behavior essentially calls `get_thread_caches(Rational{BigInt}, MAX_J[])` for the default
caches.
"""
function get_thread_caches(Tdict::Type{<:Real}, maxj::Integer)
    d3, d6 = WignerSymbols.wigner_dicts(Tdict)  # shared by all threads
    caches = BoundedWignerCache[]
    resize!(empty!(caches), Threads.nthreads())
    Threads.@threads for i in 1:Threads.nthreads()
        cache = BoundedWignerCache(d3, d6, maxj)
        caches[i] = cache
    end
    return caches
end

# utilities for StaticPrimeFactorization
isnonzero(x::T) where {T<:Integer} = (x != zero(T))
function get_last_nonzero(vec::Vector{T}) where T
    last_nonzero_index = findlast(isnonzero, vec)
    if (isnothing(last_nonzero_index))
        return 1
    end
    return last_nonzero_index
end

"""
A custom `Integer` subtype to store an integer as its prime factorization.

Unlike PrimeFactorization, this type keeps track of the index of the last nonzero
index, and it is not expected for the `powers` vector to change size. If more powers are
required, we reallocate every buffer (see `grow!`). When one computes with this number,
all loops over the factors stop at the last nonzero index, which buys ~20% in performance.
"""
mutable struct StaticPrimeFactorization{U<:Unsigned} <: AbstractPrimeFactorization{U}
    powers::Vector{U}
    sign::Int8
    last_nonzero_index::Int
end
StaticPrimeFactorization(powers::Vector{U}) where {U<:Unsigned} =
    StaticPrimeFactorization{U}(powers, one(Int8), get_last_nonzero(powers))
StaticPrimeFactorization(powers::Vector{U}, sign::Integer) where {U<:Unsigned} =
    StaticPrimeFactorization{U}(powers, sign > 0 ? one(Int8) : -one(Int8),
        get_last_nonzero(powers))

"""
This structure keeps the buffers into which prime factorizations, big integers, and Wigner
symbols are buffered. By allocating these buffers, one avoids having to reallocate for
each combination of jᵢ.

We parametrize on the type of the stored result `Tdict`.
"""
struct BoundedWignerCache{Tdict} <: AbstractWignerCache

    bounded_multiplier::Int  # 3 for 3j, 4 for 6j. determines the largets factorial
    max_j::Base.RefValue{Int}  # largest jᵢ allowed by the cache

    # the caches for each thread should all point to the same ThreadSafeDict
    Wigner3j::ThreadSafeDict{Tuple{UInt,UInt,UInt,Int,Int},Tdict}
    Wigner6j::ThreadSafeDict{NTuple{6,UInt},Tdict}

    # these arrays replace the global caches in the previous WignerSymbols.jl.
    primetable::Array{Int64,1}
    factortable::Array{Array{UInt8,1},1}
    factorialtable::Array{Array{UInt32,1},1}
    bigprimetable::Array{Array{BigInt,1},1}
    bigone::Base.RefValue{BigInt}

    # numerators and denominator buffers in order to GCD the sums
    nums::Array{StaticPrimeFactorization{UInt32},1}
    dens::Array{StaticPrimeFactorization{UInt32},1}

    # prime factorization buffers
    numbuf::StaticPrimeFactorization{UInt32}
    denbuf::StaticPrimeFactorization{UInt32}

    # buffers for adding
    addbuf::Base.RefValue{BigInt}
    seriesnum::Base.RefValue{BigInt}
    seriesden::Base.RefValue{BigInt}

    # BigInts for numerators and denominators of s ⋅ √r
    snumint::Base.RefValue{BigInt}
    sdenint::Base.RefValue{BigInt}
    rnumint::Base.RefValue{BigInt}
    rdenint::Base.RefValue{BigInt}

    # buffers for r and s calculations. s1n and s1d use numbuf and denbuf.
    s2n::StaticPrimeFactorization{UInt32}  # s2n is reused for sumlist
    snum::StaticPrimeFactorization{UInt32}
    rnum::StaticPrimeFactorization{UInt32}
    sden::StaticPrimeFactorization{UInt32}
    rden::StaticPrimeFactorization{UInt32}
end

"""
Convenience function for generating a StaticPrimeFactorization of length `nfact`.
"""
factorbuffer(nfact) = StaticPrimeFactorization(zeros(UInt32, nfact), one(Int8),1)

"""
Bounds on the prime-counting function.

Rosser, J. Barkley; Schoenfeld, Lowell (1962). "Approximate formulas for some functions of
prime numbers". Illinois J. Math. 6: 64–94.
"""
πbound(n) = (n < 2) ? 0 : Int(ceil(1.25506n / log(n)))

"""
This constructs the BoundedWignerCache.

Bounds on the maximum factorial and the number of series terms are based on Table 3 of
Johansson and Forssén 2016.

Set nj=3 for 3j, nj=6 for 6j.
"""
function BoundedWignerCache(
        dict3j::ThreadSafeDict{Tuple{UInt,UInt,UInt,Int,Int},Tdict},
        dict6j::ThreadSafeDict{NTuple{6,UInt},Tdict},
        max_j::Integer, nj=6) where {Tdict}

    # compute the bounds
    multiplier = Int(ceil(2 + round(nj / 3)))
    maxfactorial = multiplier * max_j + 1
    maxt = max_j + 1
    max_prime_factors = πbound(maxfactorial)

    cache = BoundedWignerCache{Tdict}(
        multiplier,
        Ref{Int}(max_j),
        dict3j,
        dict6j,
        [2,3,5],
        [UInt8[], UInt8[1], UInt8[0,1], UInt8[2], UInt8[0,0,1]],
        [UInt32[], UInt32[], UInt32[1], UInt32[1,1], UInt32[3,1], UInt32[3,1,1]],
        [[big(2)], [big(3)], [big(5)]],
        Ref{BigInt}(big(1)),
        [factorbuffer(max_prime_factors) for i in 1:maxt],  # nums
        [factorbuffer(max_prime_factors) for i in 1:maxt],  # dens
        factorbuffer(max_prime_factors),  # numbuf
        factorbuffer(max_prime_factors),  # denbuf

        Ref{BigInt}(big(0)),
        Ref{BigInt}(big(0)),
        Ref{BigInt}(big(0)),

        Ref{BigInt}(big(0)),
        Ref{BigInt}(big(0)),
        Ref{BigInt}(big(0)),
        Ref{BigInt}(big(0)),

        factorbuffer(max_prime_factors),  # s2n
        factorbuffer(max_prime_factors),  # snum
        factorbuffer(max_prime_factors),  # rnum
        factorbuffer(max_prime_factors),  # sden
        factorbuffer(max_prime_factors)   # rden
    )
    return cache
end

"""
Expand a StaticPrimeFactorization to a new size.
"""
function grow_factorization!(spf::StaticPrimeFactorization{T}, new_size::Integer) where T
    old_size = length(spf.powers)
    resize!(spf.powers, new_size)
    @inbounds for i in old_size:new_size
        spf.powers[i] = zero(T)
    end
end

"""
Change the maximum jᵢ of a BoundedWignerCache. This requires reallocating everything, it's
a real mess. In `wigner3j` and `wigner6j` we grow the BoundedWignerCache by a factor of 2
as necessary, i.e. amortized constant time.
"""
function grow!(new_max_j::Integer, cache::BoundedWignerCache)
    cache.max_j[] = new_max_j
    maxfactorial = cache.bounded_multiplier * new_max_j + 1
    maxt = new_max_j + 1
    max_prime_factors = πbound(maxfactorial)

    resize!(cache.nums, maxt)
    resize!(cache.dens, maxt)
    for i in 1:maxt
        cache.nums[i] = factorbuffer(max_prime_factors)
        cache.dens[i] = factorbuffer(max_prime_factors)
    end
    grow_factorization!(cache.numbuf, max_prime_factors)
    grow_factorization!(cache.denbuf, max_prime_factors)
    grow_factorization!(cache.s2n, max_prime_factors)
    grow_factorization!(cache.snum, max_prime_factors)
    grow_factorization!(cache.rnum, max_prime_factors)
    grow_factorization!(cache.sden, max_prime_factors)
    grow_factorization!(cache.rden, max_prime_factors)
end

"""
Computes the prime factors of an integer, and stores that in the cache.
"""
function precompute_primefactorial!(cache::BoundedWignerCache, n::Integer)
    n < 0 && throw(DomainError(n))
    m = length(cache.factorialtable)-1
    @inbounds while m < n
        prevfactorial = cache.factorialtable[m+1]
        m += 1
        f = primefactor(cache, m).powers
        powers = copy(prevfactorial)
        if length(f) > length(powers) # can at most be 1 larger
            push!(powers, 0)
        end
        for k = 1:length(f)
            powers[k] += f[k]
        end
        push!(cache.factorialtable, powers)
    end
end

"""
Perform an in-place multiplication of a prime factorization of a factorial.
"""
function mul_primefactorial!(cache::BoundedWignerCache,
        fact_dest::StaticPrimeFactorization{T}, n::Integer) where T
    precompute_primefactorial!(cache, n)
    row = cache.factorialtable[n+1]
    for i in 1:length(row)
        fact_dest.powers[i] += row[i]
    end
    fact_dest.last_nonzero_index = max(length(row), fact_dest.last_nonzero_index)
    return
end

"""
Deposit the factorization of an integer n to a destination fact_dest.
"""
function primefactorial!(cache::BoundedWignerCache,
                         fact_dest::StaticPrimeFactorization{T}, n::Integer) where {T, N}
    one!(fact_dest, 1)  # reset factor destination
    mul_primefactorial!(cache, fact_dest, n)
end

"""
Convenience function for multiplying together a bunch of primefactorials in a buffer.
"""
function primefactorial!(cache::BoundedWignerCache,
                         fact_dest::AbstractPrimeFactorization{T},
                         ns::NTuple{N, Integer}) where {T, N}
    one!(fact_dest, 1)  # reset factor destination to +1
    for n in ns
        mul_primefactorial!(cache, fact_dest, n)
    end
    return
end

"""
Sets a prime factorization to 1 by setting all the nonzero powers to zero. One can specify
a sign for ±1.
"""
function one!(fact_dest::StaticPrimeFactorization{T}, sign::Integer) where T
    for i in 1:fact_dest.last_nonzero_index
        fact_dest.powers[i] = zero(T)
    end
    # fill!(fact_dest.powers, zero(T))
    fact_dest.sign = sign < 0 ? -one(Int8) : one(Int8)
    fact_dest.last_nonzero_index = 1
    return
end

"""
# Style Note
Our operations which manipulate StaticPrimeFactorization don't ever decrease the last
nonzero index. It's maybe possible that something cancelled, but tests show that such an
optimization has a negligible effect on performance. We do increase it as necessary.
"""

"""
Divide common factors of StaticPrimeFactorization numbers.
"""
function divgcd!(cache::BoundedWignerCache,
        a::StaticPrimeFactorization, b::StaticPrimeFactorization)

    af, bf = a.powers, b.powers
    @inbounds for k = 1:min(a.last_nonzero_index, b.last_nonzero_index)
        gk = min(af[k], bf[k])
        af[k] -= gk
        bf[k] -= gk
    end
    return
end

"""
Find the minimum of each power between two factorizations, and put it into the first
argument `a`.
"""
function _vmin!(a::StaticPrimeFactorization{U},
    b::StaticPrimeFactorization{U}) where {U<:Unsigned}
    af = a.powers
    bf = b.powers
    @inbounds for k = 1:a.last_nonzero_index
        af[k] = min(af[k], bf[k])
    end
    return
end

"""
Subtract the powers of two StaticPrimeFactorization numbers.

The resulting number has the max of the two arguments' last nonzero indices.
"""
function _vsub!(a::StaticPrimeFactorization{U},
    b::StaticPrimeFactorization{U}) where {U<:Unsigned}
    af = a.powers
    bf = b.powers
    @inbounds for k = 1:b.last_nonzero_index
        af[k] -= bf[k]
    end
    a.last_nonzero_index = max(a.last_nonzero_index, b.last_nonzero_index)
    return a
end

"""
Add the powers of two StaticPrimeFactorization numbers.
"""
function _vadd!(a::StaticPrimeFactorization{U},
    b::StaticPrimeFactorization{U}) where {U<:Unsigned}
    af = a.powers
    bf = b.powers
    @inbounds for k = 1:b.last_nonzero_index
        af[k] = +(af[k], bf[k])
    end
    a.last_nonzero_index = max(a.last_nonzero_index, b.last_nonzero_index)
    return a
end

"""
Given a list of numerators and denominators, compute the common denominator and
the rescaled numerator after putting all fractions at the same common denominator.

This deposits the common denominator in the `cache.denbuf`.
"""
function commondenominator!(cache::BoundedWignerCache, nums::Vector{P},
        dens::Vector{P}, nk::Integer) where {P<:AbstractPrimeFactorization}

    isempty(nums) && return one(P)
    iszero(nk) && return one(P)
    # accumulate lcm of denominator

    den = cache.denbuf
    den.powers .= dens[1].powers
    den.last_nonzero_index = dens[1].last_nonzero_index

    for i = 2:nk
        den.powers .= max.(den.powers, dens[i].powers)
        den.last_nonzero_index = max(den.last_nonzero_index, dens[i].last_nonzero_index)
    end

    # rescale numerators
    for i = 1:nk
        _vadd!(nums[i], den)
        _vsub!(nums[i], dens[i])
    end
    return
end

"""
Recover a BigInteger from the factorization. This sets the input `A` to the result.
"""
function _convert!(cache::BoundedWignerCache, A::BigInt, a::StaticPrimeFactorization)
    MPZ.set!(A, cache.bigone[])
    @inbounds for n in 1:a.last_nonzero_index
        e = a.powers[n]
        if !iszero(e)
            MPZ.mul!(A, bigprime(cache, n, e))
        end
    end
    return a.sign < 0 ? MPZ.neg!(A) : A
end

"""
Recover a BigInteger from the factorization. This creates a new BigInt.
"""
function _convert(cache::BoundedWignerCache, T::Type{BigInt}, a::StaticPrimeFactorization)
    A = one(BigInt)
    return _convert!(cache, A, a)
end

"""
Auxiliary function to compute sums of a list of PrimeFactorizations as quickly as possible.

This mutates the input vector of StaticPrimeFactorization and returns the sum. This
unfortunately incurs some allocations by operating recursively.
"""
function sumlist!(cache::BoundedWignerCache, list::Vector{StaticPrimeFactorization{T}},
                  ind = 1:length(list)) where {T}
    # first compute gcd to take out common factors
    g = cache.snum
    one!(g, 1)
    first_p = list[ind[1]]
    @inbounds for i in 1:first_p.last_nonzero_index
        g.powers[i] = first_p.powers[i]
    end
    g.last_nonzero_index = first_p.last_nonzero_index
    @inbounds for k in ind
        _vmin!(g, list[k])
    end
    @inbounds for k in ind
        _vsub!(list[k], g)
    end
    gint = _convert(cache, BigInt, g)  # TODO: make this non-allocating
    L = length(ind)
    if L > 32
        l = L >> 1
        s = sumlist!(cache, list, first(ind).+(0:l-1)) + sumlist!(
                cache, list, first(ind).+(l:L-1))
    else
        # do sum
        s = big(0)
        # summand = big(1)
        for k in ind
            MPZ.add!(s, _convert!(cache, cache.addbuf[], list[k]))
        end
    end
    return MPZ.mul!(s, gint)
end

"""
Squared triangle coefficient.

This puts the results in the arguments s1n and s1d.
"""
function Δ²!(cache::BoundedWignerCache, s1n::T, s1d::T,
             j₁, j₂, j₃) where T <: StaticPrimeFactorization
    # also checks the triangle conditions by converting to unsigned integer:
    primefactorial!(cache, s1n, (
        convert(UInt, + j₁ + j₂ - j₃),
        convert(UInt, + j₁ - j₂ + j₃),
        convert(UInt, - j₁ + j₂ + j₃)))

    primefactorial!(cache, s1d, convert(UInt, j₁ + j₂ + j₃ + 1) )
end

"""
Compute the sum appearing in the 3j symbol.

This makes heavy use of preallocated variables in the BoundedWignerCache.
"""
function compute3jseries(cache::BoundedWignerCache, β₁, β₂, β₃, α₁, α₂)
    krange = max(α₁, α₂, zero(α₁)):min(β₁, β₂, β₃)
    numk = length(krange)
    T_int = eltype(eltype(cache.factorialtable))
    T = PrimeFactorization{T_int}

    # get the preallocated variables from the cache
    nums = cache.nums
    dens = cache.dens

    for (i, k) in enumerate(krange)
        num = nums[i]
        den = dens[i]
        one!(num, iseven(k) ? 1 : -1)
        primefactorial!(cache, den, (k, k-α₁, k-α₂, β₁-k, β₂-k, β₃-k))
        # divgcd!(cache, num, den)
    end

     # write gcd to buffer
    commondenominator!(cache, nums, dens, numk)
    totalden = _convert!(cache, cache.seriesden[], cache.denbuf)
    totalnum = sumlist!(cache, nums, 1:numk)
    return totalnum//totalden
end

"""
Remove one factor from each odd exponent to move into square root.
"""
function splitsquare!(s::T, r::T, a::T) where {T <: StaticPrimeFactorization}
    one!(r, a.sign)
    for i in 1:a.last_nonzero_index
        r.powers[i] = convert(UInt8, isodd(a.powers[i]))
    end

    one!(s, 1)
    for i in 1:a.last_nonzero_index
        s.powers[i] = (a.powers[i])>>1
    end

    r.last_nonzero_index = a.last_nonzero_index
    s.last_nonzero_index = a.last_nonzero_index
end

"""
Compute the Wigner-3j symbol.

The cache parametric type `Tdict` sets the type that the cache stores results. The type
`T` is what the final result will be converted to. This enables the default behavior
of storing results as pairs of r and s where r, s are Rational{BigInt} but offering
arbitrary output.
"""
function wigner3j(cache::BoundedWignerCache{Tdict}, T::Type{<:Real},
                  j₁, j₂, j₃, m₁, m₂, m₃ = -m₁-m₂) where Tdict <: Number

    # check angular momenta
    for (jᵢ,mᵢ) in ((j₁, m₁), (j₂, m₂), (j₃, m₃))
        ϵ(jᵢ, mᵢ) || throw(DomainError((jᵢ, mᵢ), "invalid combination (jᵢ, mᵢ)"))
    end
    # check triangle condition and m₁+m₂+m₃ == 0
    if !δ(j₁, j₂, j₃) || !iszero(m₁+m₂+m₃)
        return zero(T)
    end

    new_max_j = Int(ceil(max(abs(j₁), abs(j₂), abs(j₃))))
    if(cache.max_j[] < new_max_j)
        grow!(2 * new_max_j, cache)
    end

    # we reorder such that j₁ >= j₂ >= j₃ and m₁ >= 0 or m₁ == 0 && m₂ >= 0
    j₁, j₂, j₃, m₁, m₂, m₃, sgn = reorder3j(j₁, j₂, j₃, m₁, m₂, m₃)
    # TODO: do we also want to use Regge symmetries?
    α₁ = convert(Int, j₂ - m₁ - j₃ ) # can be negative
    α₂ = convert(Int, j₁ + m₂ - j₃ ) # can be negative
    β₁ = convert(Int, j₁ + j₂ - j₃ )
    β₂ = convert(Int, j₁ - m₁ )
    β₃ = convert(Int, j₂ + m₂ )

    # extra sign in definition: α₁ - α₂ = j₁ + m₂ - j₂ + m₁ = j₁ - j₂ + m₃
    sgn = isodd(α₁ - α₂) ? -sgn : sgn

    # dictionary lookup or compute
    if haskey(cache.Wigner3j, (β₁, β₂, β₃, α₁, α₂))
        rs = cache.Wigner3j[(β₁, β₂, β₃, α₁, α₂)]
    else
        # get buffered variables
        s1n, s1d, s2n = cache.numbuf, cache.denbuf, cache.s2n
        snum, rnum, sden, rden = cache.snum, cache.rnum, cache.sden, cache.rden

        # mutating versions of the functions in the main WignerSymbols.jl file
        Δ²!(cache, s1n, s1d, j₁, j₂, j₃)
        splitsquare!(sden, rden, s1d)
        primefactorial!(cache, s2n, (β₂, β₁ - α₁, β₁ - α₂, β₃, β₃ - α₁, β₂ - α₂))
        s2n_mul_s1n = _vadd!(s1n, s2n)  # multiply s2n and s1n
        splitsquare!(snum, rnum, s2n_mul_s1n)  # split square and store in snum and rnum

        s = _convert!(cache, cache.snumint[], snum) // _convert!(cache, cache.sdenint[], sden)
        r = _convert!(cache, cache.rnumint[], rnum) // _convert!(cache, cache.rdenint[], rden)

        series = compute3jseries(cache, β₁, β₂, β₃, α₁, α₂)

        # multiply s by series
        Base.GMP.MPZ.mul!(s.num, series.num)
        Base.GMP.MPZ.mul!(s.den, series.den)

        # store the r and s pair into the dictionary.
        rs = convert(Tdict, RSPair(r, s))
        cache.Wigner3j[(β₁, β₂, β₃, α₁, α₂)] = store_result(rs)  # deepcopy if necessary
    end
    return convert(T, sgn)*convert(T, rs)
end

# 6j stuff ---

function wigner6j(cache::BoundedWignerCache{Tdict}, T::Type{<:Real},
                  j₁, j₂, j₃, j₄, j₅, j₆) where Tdict <: Number
    # check validity of `jᵢ`s
    for jᵢ in (j₁, j₂, j₃, j₄, j₅, j₆)
        (ishalfinteger(jᵢ) && jᵢ >= zero(jᵢ)) || throw(DomainError("invalid jᵢ", jᵢ))
    end

    α̂₁ = (j₁, j₂, j₃)
    α̂₂ = (j₁, j₆, j₅)
    α̂₃ = (j₂, j₄, j₆)
    α̂₄ = (j₃, j₄, j₅)

    # check triangle conditions
    if !(δ(α̂₁...) && δ(α̂₂...) && δ(α̂₃...) && δ(α̂₄...))
        return zero(T)
    end
    # reduce
    α₁ = convert(UInt, +(α̂₁...))
    α₂ = convert(UInt, +(α̂₂...))
    α₃ = convert(UInt, +(α̂₃...))
    α₄ = convert(UInt, +(α̂₄...))
    β₁ = convert(UInt, j₁+j₂+j₄+j₅)
    β₂ = convert(UInt, j₁+j₃+j₄+j₆)
    β₃ = convert(UInt, j₂+j₃+j₅+j₆)

    # we should have αᵢ < βⱼ, ∀ i, j and ∑ᵢ αᵢ = ∑ⱼ βⱼ
    # now order them as β₁ >= β₂ >= β₃ >= α₁ >= α₂ >= α₃ >= α₄
    (β₁, β₂, β₃, α₁, α₂, α₃, α₄) = reorder6j(β₁, β₂, β₃, α₁, α₂, α₃, α₄)

    new_max_j = Int(ceil(max(abs(j₁), abs(j₂), abs(j₃))))
    if(cache.max_j[] < new_max_j)
        grow!(2 * new_max_j, cache)
    end

    # dictionary lookup or compute
    if haskey(cache.Wigner6j, (β₁, β₂, β₃, α₁, α₂, α₃))
        rs = cache.Wigner6j[(β₁, β₂, β₃, α₁, α₂, α₃)]
    else
        # order irrelevant: product remains the same under action of reorder6j
        n₁, d₁ = Δ²(cache, α̂₁...)
        n₂, d₂ = Δ²(cache, α̂₂...)
        n₃, d₃ = Δ²(cache, α̂₃...)
        n₄, d₄ = Δ²(cache, α̂₄...)

        snum, rnum = splitsquare(n₁ * n₂ * n₃ * n₄)
        sden, rden = splitsquare(d₁ * d₂ * d₃ * d₄)

        snu, sden = divgcd!(snum, sden)
        rnu, rden = divgcd!(rnum, rden)

        s = _convert(cache, BigInt, snum) // _convert(cache, BigInt, sden)
        r = _convert(cache, BigInt, rnum) // _convert(cache, BigInt, rden)
        s *= compute6jseries(cache, β₁, β₂, β₃, α₁, α₂, α₃, α₄)

        rs = convert(Tdict, RSPair(r, s))
        cache.Wigner6j[(β₁, β₂, β₃, α₁, α₂, α₃)] = store_result(rs)
    end
    return convert(T, rs.s) * convert(T, signedroot(rs.r))
end

# compute the sum appearing in the 6j symbol
function compute6jseries(cache::BoundedWignerCache, β₁, β₂, β₃, α₁, α₂, α₃, α₄)
                         krange = max(α₁, α₂, α₃, α₄):min(β₁, β₂, β₃)
                         T = PrimeFactorization{eltype(eltype(cache.factorialtable))}
    numk = length(krange)
    nums = cache.nums
    dens = cache.dens

    for (i, k) in enumerate(krange)
        num = nums[i]
        den = dens[i]
        primefactorial!(cache, num, k+1)
        num.sign = Int8(iseven(k) ? 1 : -1)
        primefactorial!(cache, den, (k-α₁, k-α₂, k-α₃, k-α₄, β₁-k, β₂-k, β₃-k))
        divgcd!(cache, num, den)
    end

    commondenominator!(cache, nums, dens, numk)
    totalden = _convert!(cache, cache.seriesden[], cache.denbuf)
    totalnum = sumlist!(cache, nums, 1:numk)

    return totalnum//totalden
end


text:
using Flux
using NeuralStyle
using Random
using Test

@testset "Image Size" begin
# Make a set of randomly sized arrays with random kernels/strides
# and ensure that the input and output sizes are equal after convolutions
    @info "Testing image size preservation"
    equal_size = []  # Array to hold boolean values
    for i in 1:10
        # Make an array of random floats with random sizes
        input_tensor =  rand(rand(100:500), rand(100:500), 3, 1);
        size_in = size(input_tensor)[1:2]  # size of input array

        k = (rand(1:9), rand(1:9))  # kernel size
        s = (rand(1:9), rand(1:9))  # stride size
        d = (1, 1)
        epad1 = pad2d(size_in, k, s)  # Basic padding calculation
        layer2size = nstrides2d(size_in, k, epad1, s)  # Estimate the size of layer 2
        epad2 = pad2d(layer2size, k, s)

        autoencode = Chain(Conv(k, 3=>8, relu, pad=epad1, stride=s, dilation=d),
                           Conv(k, 8=>16, relu, pad=epad2, stride=s, dilation=d),
                           ConvTranspose(k, 16=>8,relu, pad=epad2, stride=s, dilation=d),
                           ConvTranspose(k, 8=>3, relu, pad=epad1, stride=s, dilation=d))

        #println("\nPrinting results")
        encoded = autoencode(input_tensor)
        push!(equal_size, size(input_tensor) == size(encoded))
    end
    @test sum(equal_size) == 10  # Should be 10, (true=1)
end

ext: jl


text:
## Copyright (c) 2013 Miguel Bazdresch
##
## This file is distributed under the 2-clause BSD License.

# append x,y,z coordinates and configuration to current figure
function addcoords(x::Coord,y::Coord,Z::Array,conf::CurveConf)
    global gnuplot_state

    # check that at least one figure has been setup
    if gnuplot_state.current == 0
        figure(1)
    end

    # check coordinates: dimensions, sizes and types
    nex = !isempty(x); ney = !isempty(y); neZ = !isempty(Z)
    # check types
    if nex
        @assert(eltype(x)<:Real,"Invalid coordinates")
    end
    if ney
        @assert(eltype(y)<:Real,"Invalid coordinates")
    end
    if neZ
        @assert(eltype(Z)<:Real,"Invalid coordinates")
        @assert(1 < ndims(Z) < 4,"Invalid coordinates")
    end
    # fill missing x, y coordinates when Z is not empty
    if neZ && !nex
        x = 1:size(Z,2)
    end
    if neZ && !ney
        y = 1:size(Z,1)
    end
    # if x,y are matrices, convert to vectors
    if isa(x,Matrix)
        s = size(x)
        if s[1] == 1 || s[2] == 1
            x = squeeze(x)
        else
            error("Invalid abscissa coordinates")
        end
    elseif isa(x,UnitRange) || isa(x,Range)
        x = collect(x)
    end
    if isa(y,Matrix)
        s = size(y)
        if s[1] == 1 || s[2] == 1
            y = squeeze(y)
        else
            error("Invalid abscissa coordinates")
        end
    elseif isa(y,UnitRange) || isa(y,Range)
        y = collect(y)
    end
    # check number of elements
    if neZ
        if conf.plotstyle == "image" || conf.plotstyle == "rgbimage"
            @assert(size(Z,2) == length(x), "Wrong number of columns in Z")
            @assert(size(Z,1) == length(y), "Wrong number of rows in Z")
        else
            @assert(size(Z,1) == length(x), "Wrong number of columns in Z")
            @assert(size(Z,2) == length(y), "Wrong number of rows in Z")
        end
    else
        @assert(length(x) == length(y),
        "Abscissa and ordinate must have the same number of elements")
    end

    # check curve configuration: property names
    # TODO: figure out how to check valid color names -- gnuplot supports
    #  112 different color names.
    # check valid values of marker
    @assert(validate_marker(conf.marker), "Invalid mark name specified")
    if nex && !neZ ## 2-d plot
        # check valid values of plotstyle
        @assert(validate_2d_plotstyle(conf.plotstyle),
            "Invalid plotstyle specified")
    elseif nex && neZ ## 3-d plot or image
        @assert(validate_3d_plotstyle(conf.plotstyle),
            "Invalid plotstyle specified")
    end

    conf = copy(conf)       # we need to dereference conf
    # append data to figure
    c = findfigure(gnuplot_state.current)
    fig = gnuplot_state.figs[c]
    if fig.isempty == true
        fig.curves[1] = CurveData(x,y,Z,conf)
        fig.isempty = false
    else
        push!(fig.curves, CurveData(x,y,Z,conf))
    end
    gnuplot_state.figs[c] = fig
end
addcoords(y) = addcoords(1:length(y),y,Any[],CurveConf())
addcoords(y,c::CurveConf) = addcoords(1:length(y),y,Any[],c)
addcoords(x,y) = addcoords(x,y,Any[],CurveConf())
addcoords(x,y,c::CurveConf) = addcoords(x,y,Any[],c)
addcoords(x,y,Z) = addcoords(x,y,Z,CurveConf())
# X, Y data in matrix columns
function addcoords(X::Matrix,Y::Matrix,conf::CurveConf)
    for i = 1:size(X,2)
        addcoords(X[:,i],Y[:,i],Any[],conf)
    end
end
function addcoords(Y::Matrix,conf::CurveConf)
    x = 1:size(Y,1)
    X = zeros(size(Y))
    for i = 1:size(Y,2)
        X[:,i] = x
    end
    addcoords(X,Y,conf)
end
addcoords(X::Matrix, Y::Matrix) = addcoords(X,Y,CurveConf())
addcoords(Y::Matrix) = addcoords(Y,CurveConf())

# append error data to current set of coordinates
function adderror(yl::Coord,yh::Coord)
    global gnuplot_state
    # check that at least one figure has been setup
    if gnuplot_state.current == 0
        figure(1)
    end
    c = findfigure(gnuplot_state.current)

    # check arguments and convert to vectors
    if isempty(yl)
        error("Invalid error data")
    else
        @assert(eltype(yl)<:Real,"Invalid error data")
        if isa(yl,Matrix)
            s = size(yl)
            if s[1] == 1 || s[2] == 1
                yl = squeeze(yl)
            else
                error("Invalid error data")
            end
        elseif isa(yl,UnitRange) || isa(yl,Range)
            yl = collect(yl)
        end
    end
    if !isempty(yh)
        @assert(eltype(yh)<:Real,"Invalid error data")
        if isa(yl,Matrix)
            s = size(yh)
            if s[1] == 1 || s[2] == 1
                yh = squeeze(yh)
            else
                error("Invalid error data")
            end
        elseif isa(yh,UnitRange) || isa(yh,Range)
            yh = collect(yh)
        end
    end
    # verify vector sizes -- this also implies that x,y coordinates must be
    # added to figure, before error data can be attached to it
    @assert(length(gnuplot_state.figs[c].curves[end].x) == length(yl),
        "Error data vector must be of same size as abscissa")
    if !isempty(yh)
        @assert(length(yh) == length(yl),
            "Error data vectors must be of same size")
    end

    # set fields in current curve
    gnuplot_state.figs[c].curves[end].ylow = yl
    gnuplot_state.figs[c].curves[end].yhigh = yh

end
adderror(ydelta) = adderror(ydelta,Any[])

# add financial data to current set of coordinates
# input order is: open, low, high, close
function addfinancial(args...)
	global gnuplot_state
	    # check that at least one figure has been setup
    if gnuplot_state.current == 0
        figure(1)
    end
    c = findfigure(gnuplot_state.current)

    # check arguments and convert to vectors
    if length(args) != 4
        error("Invalid financial data")
	end
	for a in args
		@assert(!isempty(a),"Invalid financial data")
        @assert(eltype(a)<:Real,"Invalid financial data")
        @assert(typeof(a)<:Array,"Invalid financial data")
		@assert(length(gnuplot_state.figs[c].curves[end].x) == length(a),
			"Financial data vector must be of same size as abscissa")
	end

    # set fields in current curve
    gnuplot_state.figs[c].curves[end].finance.open = args[1]
    gnuplot_state.figs[c].curves[end].finance.low = args[2]
    gnuplot_state.figs[c].curves[end].finance.high = args[3]
    gnuplot_state.figs[c].curves[end].finance.close = args[4]

end

# add axes configuration to current figure
function addconf(conf::AxesConf)
    global gnuplot_state
    # check that at least one figure has been setup
    if gnuplot_state.current == 0
        figure(1)
    end

    # Perform argument validation
    # TODO: find a way to validate the box argument
    # validate axis type
    @assert(validate_axis(conf.axis),"Invalid axis type specified")
    # validate ranges
    @assert(validate_range(conf.xrange),"Invalid xrange specified")
    @assert(validate_range(conf.yrange),"Invalid yrange specified")
    @assert(validate_range(conf.zrange),"Invalid zrange specified")

    conf = copy(conf)
    # select current plot
    c = findfigure(gnuplot_state.current)
    gnuplot_state.figs[c].conf = conf
end

# llplot() is our workhorse plotting function
function llplot()
    global gnuplot_state
    global gaston_config

    # select current plot
    c = findfigure(gnuplot_state.current)
    if c == 0
        println("No current figure")
        return
    end
    figs = gnuplot_state.figs
    config = figs[c].conf

    # Build terminal setup string and send it to gnuplot
    ts = termstring(gaston_config.terminal)
    gnuplot_send(ts)

    # if figure has no data, stop here
    if isempty(figs[c].curves[1].x)
        return
    end

    # datafile filename
    (filename,f) = mktemp(gnuplot_state.tmpdir)

    # Send appropriate coordinates and data to gnuplot, depending on
    # whether we are doing 2-d, 3-d or image plots.

    # 2-d plot: Z is empty or plostyle is {,rgb}image
    if isempty(figs[c].curves[1].Z) ||
        figs[c].curves[1].conf.plotstyle == "image" ||
        figs[c].curves[1].conf.plotstyle == "rgbimage"
        # create data file
        for i in figs[c].curves
            ps = i.conf.plotstyle
            if ps == "errorbars" || ps == "errorlines"
                if isempty(i.yhigh)
                    # ydelta (single error coordinate)
                    writedlm(f,[i.x i.y i.ylow],' ')
                else
                    # ylow, yhigh (double error coordinate)
                    writedlm(f,[i.x i.y i.ylow i.yhigh],' ')
                end
			elseif ps == "financebars"
            	# data is written to tmparr, which is then written to disk
            	tmparr = zeros(length(i.x),5)
                # output matrix
                for col = 1:length(i.x)
					tmparr[col,1] = i.x[col]
					tmparr[col,2] = i.finance.open[col]
					tmparr[col,3] = i.finance.low[col]
					tmparr[col,4] = i.finance.high[col]
					tmparr[col,5] = i.finance.close[col]
                end
				writedlm(f,tmparr,' ')
            elseif ps == "image"
            	# data is written to tmparr, which is then written to disk
            	tmparr = zeros(length(i.x)*length(i.y),3)
                tmparr_row_index = 1  # index into tmparr row
                # output matrix
                for row = 1:length(i.y)
                    x = length(i.x)
                    for col = 1:length(i.x)
                    	tmparr[tmparr_row_index,1] = i.x[col]
                    	tmparr[tmparr_row_index,2] = i.x[row]
                    	tmparr[tmparr_row_index,3] = i.Z[row,col]
                    	tmparr_row_index = tmparr_row_index+1
                        x = x-1
                    end
                end
				writedlm(f,tmparr,' ')
            elseif ps == "rgbimage"
            	# data is written to tmparr, which is then written to disk
            	tmparr = zeros(length(i.x)*length(i.y), 5)
            	tmparr_row_index = 1
                # output matrix
                for col = 1:length(i.x)
                    y = length(i.y)
                    for row = 1:length(i.y)
                    	tmparr[tmparr_row_index,1] = i.x[col]
                    	tmparr[tmparr_row_index,2] = i.y[row]
                    	tmparr[tmparr_row_index,3] = i.Z[y,col,1]
                    	tmparr[tmparr_row_index,4] = i.Z[y,col,2]
                    	tmparr[tmparr_row_index,5] = i.Z[y,col,3]
                    	tmparr_row_index = tmparr_row_index+1
                        y = y-1
                    end
                end
                writedlm(f,tmparr,' ')
            else
                writedlm(f,[i.x i.y],' ')
            end
            write(f,"\n\n")
        end
        close(f)
        # send figure configuration to gnuplot
        gnuplot_send_fig_config(config)
        # send plot command to gnuplot
        gnuplot_send(linestr(figs[c].curves, "plot", filename))

    # 3-d plot: Z is not empty and plotstyle is not {,rgb}image
    elseif !isempty(figs[c].curves[1].Z) &&
            figs[c].curves[1].conf.plotstyle != "image" &&
            figs[c].curves[1].conf.plotstyle != "rgbimage"
        # create data file
        for i in figs[c].curves
			# data is written to tmparr, which is then written to disk
			tmparr = zeros(1, 3)
			tmparr_row_index = 1
            for row in 1:length(i.x)
                for col in 1:length(i.y)
					tmparr[1,1] = i.x[row]
					tmparr[1,2] = i.y[col]
					tmparr[1,3] = i.Z[row,col]
					writedlm(f,tmparr,' ')
                end
                write(f,"\n")
            end
            write(f,"\n\n")
        end
        close(f)
        # send figure configuration to gnuplot
        gnuplot_send_fig_config(config)
        # send command to gnuplot
        gnuplot_send(linestr(figs[c].curves, "splot",filename))
    end
    # Wait until gnuplot is finished plotting before returning. To do this,
    # we make gnuplot output "X\n" in its stdout. Gnuplot will only get to
    # do this when the plot is finished. Otherwise, gnuplot will output
    # something in its stderr. In either case, we know that the plot is
    # finished and can carry on.
    gnuplot_send("set print \"-\"\n")
    gnuplot_send("print \"X\"\n")

	# Loop until we read data from either gnuplot's stdout or stdin
	while true
		yield()  # let async tasks run
		if !isempty(gnuplot_state.gp_stdout)
			# We got data from stdin, meaning the plot went well.
			gnuplot_state.gp_stdout = ""
			gnuplot_state.gp_error = false
			break
		end
		if !isempty(gnuplot_state.gp_stderr)
			# Gnuplot met trouble while plotting.
			gnuplot_state.gp_lasterror = gnuplot_state.gp_stderr
			gnuplot_state.gp_stderr = ""
			gnuplot_state.gp_error = true
			break
		end
	end

    # Reset gnuplot settable options.
    gnuplot_send("\nreset\n")

    # If the environment is IJulia, redisplay the figure.
    if displayable("image/png")
    	redisplay(figs[c])
	end
end


text:
struct Unicycle{T<:Real}
    ΔT::T
end

function Base.getproperty(system::Unicycle, sym::Symbol)
    if sym === :n_states
        4
    elseif sym === :n_controls
        2
    else
        getfield(system, sym)
    end
end

function DynamicsModelInterface.next_x(system::Unicycle, x_t, u_t)
    ΔT = system.ΔT
    @assert only(size(x_t)) == 4
    @assert only(size(u_t)) == 2
    px, py, v, θ = x_t
    Δv, Δθ = u_t
    [px + ΔT * v * cos(θ), py + ΔT * v * sin(θ), v + Δv, θ + Δθ]
end

# These constraints encode the dynamics of a unicycle with state layout x_t = [px, pyL, v, θ] and
# inputs u_t = [Δv, Δθ].
function DynamicsModelInterface.add_dynamics_constraints!(system::Unicycle, opt_model, x, u)
    ΔT = system.ΔT
    T = size(x, 2)

    # auxiliary variables for nonlinearities
    cosθ = @variable(opt_model, [1:T])
    @NLconstraint(opt_model, [t = 1:T], cosθ[t] == cos(x[4, t]))
    sinθ = @variable(opt_model, [1:T])
    @NLconstraint(opt_model, [t = 1:T], sinθ[t] == sin(x[4, t]))

    @constraint(
        opt_model,
        [t = 1:(T - 1)],
        x[:, t + 1] .== [
            x[1, t] + ΔT * x[3, t] * cosθ[t],
            x[2, t] + ΔT * x[3, t] * sinθ[t],
            x[3, t] + u[1, t],
            x[4, t] + u[2, t],
        ]
    )
end

function DynamicsModelInterface.add_dynamics_jacobians!(system::Unicycle, opt_model, x, u)
    ΔT = system.ΔT
    n_states, T = size(x)
    n_controls = size(u, 1)

    # auxiliary variables for nonlinearities
    cosθ = @variable(opt_model, [1:T])
    @NLconstraint(opt_model, [t = 1:T], cosθ[t] == cos(x[4, t]))
    sinθ = @variable(opt_model, [1:T])
    @NLconstraint(opt_model, [t = 1:T], sinθ[t] == sin(x[4, t]))

    # jacobians of the dynamics in x
    dfdx = @variable(opt_model, [1:n_states, 1:n_states, 1:T])
    @constraint(
        opt_model,
        [t = 1:T],
        dfdx[:, :, t] .== [
            1 0 ΔT*cosθ[t] -ΔT*x[3, t]*sinθ[t]
            0 1 ΔT*sinθ[t] +ΔT*x[3, t]*cosθ[t]
            0 0 1 0
            0 0 0 1
        ]
    )

    # jacobians of the dynamics in u
    dfdu = [
        0 0
        0 0
        1 0
        0 1
    ] .* reshape(ones(T), 1, 1, :)

    (; dx = dfdx, du = dfdu)
end

#========================================== Visualization ==========================================#

function TrajectoryVisualization.trajectory_data(::Unicycle, x, player = 1)
    [
        (; px = xi[1], py = xi[2], v = xi[3], θ = xi[4], player, t) for
        (t, xi) in enumerate(eachcol(x))
    ]
end


text:
import Base: +, -, *, ^, |, ~, dot, cross
+(a::Coord, b::Coord)     = Coord(a.x+b.x, a.y+b.y, a.z+b.z)
-(a::Coord, b::Coord)     = Coord(a.x-b.x, a.y-b.y, a.z-b.z)
*(a::Float64, b::Coord)   = Coord(a*b.x, a*b.y, a*b.z)
*(b::Coord, a::Float64,)  = Coord(a*b.x, a*b.y, a*b.z)
*(a::Int, b::Coord)       = Coord(a*b.x, a*b.y, a*b.z)
*(b::Coord, a::Int)       = Coord(a*b.x, a*b.y, a*b.z)
dot(a::Coord, b::Coord)   = (a.x*b.x + a.y*b.y + a.z*b.z)
"""
    magnitude(a::Coord)

A utility function to determine the magnitude of a `Coord` object. Typical use case is to subtract two Coord objects and check the resulting Coord object's magnitude to determine the distance between the two Coords.
"""
magnitude(a::Coord)       = sqrt(dot(a,a))
"""
    unitize(a::Coord)

A utility function to unitize a `Coord`
"""
unitize(a::Coord)         = (1. / magnitude(a) * a)
cross(a::Coord, b::Coord) = Coord(a.y*b.z - a.z*b.y, a.z*b.x - a.x*b.z, a.x*b.y - a.y*b.x)


text:
## SubOperator

immutable SubOperator{T,B,I,DI,BI} <: Operator{T}
    parent::B
    indexes::I
    dims::DI
    bandwidths::BI
end

SubOperator(A,inds,dims,lu) =
    SubOperator{eltype(A),typeof(A),typeof(inds),
                typeof(dims),typeof(lu)}(A,inds,map(length,inds),lu)

# cannot infer ranges
SubOperator(A,inds,dims) = SubOperator(A,inds,dims,(dims[1]-1,dims[2]-1))
SubOperator(A,inds) = SubOperator(A,inds,map(length,inds))

function view(A::Operator,kr::AbstractCount,jr::AbstractCount)
    @assert isbanded(A) && isinf(size(A,1)) && isinf(size(A,2))
    st=step(kr)
    @assert st==step(jr)  # Otherwise, its not a banded operator
    kr1=first(kr)
    jr1=first(jr)
    l,u=(bandinds(A,1)+kr1-jr1)÷st,(bandinds(A,2)+kr1-jr1)÷st
    SubOperator(A,(kr,jr),size(A),(-l,u))
end

function view(A::Operator,kr::Range,jr::Range)
    st=step(kr)
    if isbanded(A) && st == step(jr)
        kr1=first(kr)
        jr1=first(jr)
        l,u=(bandinds(A,1)+kr1-jr1)÷st,(bandinds(A,2)+kr1-jr1)÷st
        SubOperator(A,(kr,jr),(length(kr),length(jr)),(-l,u))
    else
        SubOperator(A,(kr,jr))
    end
end

function view(A::Operator,kr::UnitRange,jr::UnitRange)
    if isbanded(A)
        shft=first(kr)-first(jr)
        l,u=max(bandwidth(A,1)-shft,0),max(bandinds(A,2)+shft,0)
        SubOperator(A,(kr,jr),(length(kr),length(jr)),(l,u))
    else
        SubOperator(A,(kr,jr))
    end
end

view(A::Operator,::Colon,jr::Range) = view(A,1:size(A,1),jr)
view(A::Operator,kr::Range,::Colon) = view(A,kr,1:size(A,2))

view(A::SubOperator,kr::UnitRange,jr::UnitRange) =
    view(A.parent,A.indexes[1][kr],A.indexes[2][jr])

bandwidth(S::SubOperator,k::Int) = S.bandwidths[k]
bandinds(S::SubOperator) = (-bandwidth(S,1),bandwidth(S,2))
function colstop(S::SubOperator,j::Integer)
    cs=findfirst(parentindexes(S)[1],colstop(parent(S),parentindexes(S)[2][j]))
    cs==0?size(S,1):cs
end
colstart(S::SubOperator,j::Integer) =
    max(findfirst(parentindexes(S)[1],colstart(parent(S),parentindexes(S)[2][j])),1)
rowstart(S::SubOperator,j::Integer) =
    max(1,findfirst(parentindexes(S)[2],rowstart(parent(S),parentindexes(S)[1][j])))
function rowstop(S::SubOperator,j::Integer)
    rs=findfirst(parentindexes(S)[2],rowstop(parent(S),parentindexes(S)[1][j]))
    rs==0?size(S,2):rs
end

function bbbzeros(S::SubOperator)
    kr,jr=parentindexes(S)
    KO=parent(S)
    l,u=blockbandinds(KO)
    λ,μ=subblockbandinds(KO)

    rt=rangetensorizer(KO)
    dt=domaintensorizer(KO)
    bbbzeros(eltype(KO),-l,u,-λ,μ,
            blocklengthrange(rt,kr),
            blocklengthrange(dt,jr))
end

domainspace(S::SubOperator) = SubSpace(domainspace(parent(S)),parentindexes(S)[2])
rangespace(S::SubOperator) = SubSpace(rangespace(parent(S)),parentindexes(S)[1])

size(V::SubOperator) = V.dims
size(V::SubOperator,k::Int) = V.dims[k]
unsafe_getindex(V::SubOperator,k::Integer,j::Integer) = V.parent[V.indexes[1][k],V.indexes[2][j]]
getindex(V::SubOperator,k::Integer,j::Integer) = V.parent[V.indexes[1][k],V.indexes[2][j]]
getindex(V::SubOperator,k::Integer,j::Range) = V.parent[V.indexes[1][k],V.indexes[2][j]]
getindex(V::SubOperator,k::Range,j::Integer) = V.parent[V.indexes[1][k],V.indexes[2][j]]
getindex(V::SubOperator,k::Range,j::Range) = V.parent[V.indexes[1][k],V.indexes[2][j]]
Base.parent(S::SubOperator) = S.parent
Base.parentindexes(S::SubOperator) = S.indexes


text:
include("bounds.jl")

# Compute the negative log prior:

function nlog_prior_noecc(x::Array{T,1}) where {T <: Real}
# We will place a joint prior on eccentricity vector
# such that each planet has an eccentricity which lies between
# 0 and emax2, with a gradual decrease from emax1 to emax2:
# Eccentricity priors:
emax1 = 0.2; emax2 = 0.3
# Mass priors:
lgm1 = -8.0; lgm2 = -7.0
lgm3 = -3.0; lgm4 = -2.0
# Period priors:
P1 = [1.49,2.40,4.03,6.08,9.18,12.33,18.75]; P2 = [1.50,2.41,4.04,6.09,9.19,12.34,18.76]
P3 = [1.52,2.43,4.06,6.11,9.22,12.36,18.78]; P4 = [1.53,2.44,4.07,6.12,9.23,12.37,18.79]
# t0 priors:
t01 = [7257.53,7258.57,7257.05,7257.81,7257.05,7257.70,7249.58]
t02 = [7257.54,7258.58,7257.06,7257.82,7257.06,7257.71,7249.59]
t03 = [7257.56,7258.60,7257.08,7257.84,7257.08,7257.73,7249.61]
t04 = [7257.57,7258.61,7257.09,7257.85,7257.09,7257.74,7249.62]
# Define -log(prior):
nlprior = zero(T)
nparam = size(x)[1]
n = div(nparam-2,5)+1
# Loop over planets:
for i=2:n
  # Place prior on eccentricities:
  ecc = sqrt(x[(i-2)*5+4]^2+x[(i-2)*5+5]^2)
  nlprior_tmp,dpdx = log_bounds_upper(ecc,emax1,emax2)
  nlprior -= nlprior_tmp

  # We will place broad lower and upper limits on log_10(mass):
  log10mass = log10(x[(i-2)*5+1])
  nlprior_tmp,dpdx = log_bounds_lower(log10mass,lgm1,lgm2)
  nlprior -= nlprior_tmp
  nlprior_tmp,dpdx = log_bounds_upper(log10mass,lgm3,lgm4)
  nlprior -= nlprior_tmp
  # Broad priors on period:
  P = x[(i-2)*5+2]
  nlprior_tmp,dpdx = log_bounds_lower(P,P1[i-1],P2[i-1])
  nlprior -= nlprior_tmp
  nlprior_tmp,dpdx = log_bounds_upper(P,P3[i-1],P4[i-1])
  nlprior -= nlprior_tmp
  # Broad priors on t0:
  t0 = x[(i-2)*5+3]
  nlprior_tmp,dpdx = log_bounds_lower(t0,t01[i-1],t02[i-1])
  nlprior -= nlprior_tmp
  nlprior_tmp,dpdx = log_bounds_upper(t0,t03[i-1],t04[i-1])
  nlprior -= nlprior_tmp
end
# Add in prior on ndof & lnV1:
ndof1 = 0.5; ndof2 = 1.0
ndof3 = 50.0; ndof4 = 100.0
ndof = exp(x[nparam-1])
nlprior_tmp,dpdx = log_bounds_lower(ndof,ndof1,ndof2)
nlprior -= nlprior_tmp
nlprior_tmp,dpdx = log_bounds_upper(ndof,ndof3,ndof4)
nlprior -= nlprior_tmp
lnV11 = -2.0; lnV12 = -1.0
lnV13 = 5.0; lnV14 = 10.0
lnV1 = log(x[nparam]/exp(0.5/ndof))
nlprior_tmp,dpdx = log_bounds_lower(lnV1,lnV11,lnV12)
nlprior -= nlprior_tmp
nlprior_tmp,dpdx = log_bounds_upper(lnV1,lnV13,lnV14)
nlprior -= nlprior_tmp
# That's it!
return nlprior::T
end

# Take gradient and hessian:
grad_nlog_prior_noecc= x -> ForwardDiff.gradient(nlog_prior_noecc,x)
hess_nlog_prior_noecc= x -> ForwardDiff.hessian(nlog_prior_noecc,x)


text:
VERSION < v"0.7.0-beta2.199" && __precompile__()

module CMakeWrapper

using Parameters: @with_kw
using BinDeps
using BinDeps: BuildProcess,
               BuildStep,
               @dependent_steps,
               LibraryDependency,
               gethelper,
               builddir,
               stringarray,
               adjust_env
using Libdl
import BinDeps: lower,
                provider,
                generate_steps
export cmake_executable, CMakeBuild, CMakeProcess

import CMake
const cmake_executable = CMake.cmake

const dlext = Libdl.dlext

@with_kw struct CMakeBuild <: BuildStep
    srcdir::AbstractString = ""
    builddir::AbstractString = ""
    prefix::AbstractString = ""
    libtarget::Vector{AbstractString} = String[]
    installed_libpath::Vector{AbstractString} = String[]
    cmake_args::Vector{AbstractString} = String[]
    configname::AbstractString = ""
    targetname::AbstractString = "install"
    env::Dict{Any, Any} = Dict{Any, Any}()
end

function lower(s::CMakeBuild, collection)
    env = adjust_env(s.env)
    cmake_command = `$cmake_executable -DCMAKE_INSTALL_PREFIX=$(s.prefix)`
    for arg in s.cmake_args
        cmake_command = `$cmake_command $arg`
    end
    build_command = `$cmake_executable --build .`
    if !isempty(s.configname)
        build_command = `$build_command --config $(s.configname)`
    end
    if !isempty(s.targetname)
        build_command = `$build_command --target $(s.targetname)`
    end
    @dependent_steps begin
        CreateDirectory(s.builddir)
        begin
            ChangeDirectory(s.builddir)
            FileRule(joinpath(s.builddir, "CMakeCache.txt"),
                     setenv(`$cmake_command $(s.srcdir)`, env))
            FileRule(s.installed_libpath,
                     setenv(build_command, env))
        end
    end
end

mutable struct CMakeProcess <: BuildProcess
    source
    opts
end

CMakeProcess(; opts...) = CMakeProcess(nothing, Dict{Any, Any}(opts))

provider(::Type{CMakeProcess}, cm::CMakeProcess; opts...) = cm

function generate_steps(dep::LibraryDependency, h::CMakeProcess, provider_opts)
    # Shamelessly copied from BinDeps.jl, which is  distributed under
    # the MIT License, and Copyright (c) 2012: Keno Fischer and other
    # contributors. See LICENSE.md for license terms.
    if h.source === nothing
        h.source = gethelper(dep,Sources)
    end
    if isa(h.source,Sources)
        h.source = (h.source,Dict{Symbol,Any}())
    end
    h.source[1] === nothing && error("Could not obtain sources for dependency $(dep.name)")
    steps = lower(generate_steps(dep,h.source...))
    opts = Dict{Symbol, Any}(
        :srcdir => srcdir(dep,h.source...),
        :prefix => usrdir(dep),
        :builddir => joinpath(builddir(dep),dep.name),
        h.opts...
    )
    if haskey(opts,:installed_libname)
        !haskey(opts,:installed_libpath) || error("Can't specify both installed_libpath and installed_libname")
        opts[:installed_libpath] = String[joinpath(libdir(dep),opts[:installed_libname])]
        delete!(opts, :installed_libname)
    elseif !haskey(opts,:installed_libpath)
        opts[:installed_libpath] = String[joinpath(libdir(dep),x)*"."*dlext for x in stringarray(get(dep.properties,:aliases,String[]))]
    end
    if !haskey(opts,:libtarget) && haskey(dep.properties,:aliases)
        opts[:libtarget] = String[x*"."*dlext for x in stringarray(dep.properties[:aliases])]
    end
    env = Dict{String,String}()
    if Sys.isunix()
        env["PATH"] = bindir(dep)*":"*ENV["PATH"]
    elseif Sys.iswindows()
        env["PATH"] = bindir(dep)*";"*ENV["PATH"]
    end
    haskey(opts,:env) && merge!(env,opts[:env])
    opts[:env] = env
    if get(provider_opts, :force_rebuild, false)
        @warn "force_rebuild option is not supported for CMakeProcess. It will be ignored"
    end
    steps |= CMakeBuild(; opts...)
    steps
end

end # module


text:
using HWconstrained
using Test

@testset "HWconstrained.jl" begin

@testset "testing components" begin
d=data()
A=d["a"]
u(x)=-exp(-A*x)
grad=ones(4)

@testset "tests gradient of objective function" begin

end

@testset "tests gradient of constraint function" begin

end

end

@testset "testing result of both maximization methods" begin
truth = HWconstrained.DataFrame(a=[0.5;1.0;5.0],

c = [1.008;1.004;1.0008],

omega1=[-1.41237;-0.20618;0.758763],

omega2=[0.801455;0.400728;0.0801455],

omega3=[1.60291;0.801455;0.160291],

fval=[-1.20821;-0.732819;-0.013422])

tol2=1e-2

@testset "checking result of NLopt maximization" begin
t1 = table_JuMP()
for c in names(truth)
    @test all(maximum.(abs.(t1[c].-truth[c])) .< tol2)
end
end

@testset "checking result of NLopt maximization" begin
@test all(maximum.(abs.(t1[c].-truth[c])) .< tol2)
end
end
end
end


text:
root = RectanglePacker(Rect2D(0,0,1024,1024))
push!(root, Rect2D(0,0,20,20))
push!(root, Rect2D(0,0,20,20))
push!(root, [Rect2D(0,0, rand(5:50), rand(5:50)) for i=1:20])

# function get_rectangles(packer::Nothing, rectangles=IRect2D[])
#     return rectangles
# end
#
# function get_rectangles(packer, rectangles)
#     push!(rectangles, packer.area)
#     get_rectangles(packer.children.left, rectangles)
#     get_rectangles(packer.children.right, rectangles)
#     return rectangles
# end
# function get_rectangles(packer)
#     rectangles = IRect2D[]
#     get_rectangles(packer.children.left, rectangles)
#     get_rectangles(packer.children.right, rectangles)
#     return rectangles
# end
#
# rectangles = get_rectangles(root)
# using Makie
# linesegments(FRect(0, 0, 1024, 1024))
# poly!(rectangles, color = (:red, 0.1), strokewidth=1, strokecolor=:black)
#


text:
push!(LOAD_PATH, "../src/")

using Documenter, ExactPredicates

makedocs(sitename="ExactPredicates.jl",
         pages = ["index.md", "api.md"],
         repo = "https://github.com/lairez/ExactPredicates.jl/blob/master{path}#{line}")


text:
module TestJuliaREPL

include("preamble.jl")
using IPython: init_repl, afterreplinit, REPL

@static if VERSION >= v"0.7.0-"
    using REPL: TextTerminal
else
    using Base.Terminals: TextTerminal
end

@testset "afterreplinit" begin
    repl = afterreplinit(identity)
    if isinteractive()
        @test repl === Base.active_repl
    else
        @test repl === nothing
    end
end

struct DummyTerminal <: TextTerminal
end

function dummy_repl()
    if VERSION >= v"0.7.0-"
        hascolor = false
        repl = REPL.LineEditREPL(DummyTerminal(), hascolor)
    else
        repl = REPL.LineEditREPL(DummyTerminal())
    end
    repl.interface = REPL.setup_interface(repl)
    return repl
end

function dummy_initialized_repl()
    repl = dummy_repl()
    init_repl(repl)
    return repl
end

@testset "init_repl" begin
    repl = dummy_initialized_repl()
    keymap_dict = repl.interface.modes[1].keymap_dict
    @test haskey(keymap_dict, '.')
end

end  # module


text:
using Yao, YaoExtensions, Yao.AD
using Yao.AD: generator
using Test, Random

@testset "BP diff" begin
    Random.seed!(2)
    c = put(4, 3=>Rx(0.5))

    circuit = chain(4, repeat(4, H, 1:4), put(4, 3=>Rz(0.5)), control(4, 2, 1=>shift(0.4)), control(2, 1=>X), put(4, 4=>Ry(0.2)))
    op = heisenberg(4)
    θ = [0.9, 0.2, 0.3]
    dispatch!(circuit, θ)

    for ψ0 in [rand_state(4), rand_state(4, nbatch=10)]
        dispatch!(circuit, θ)
        ψδ, g1 = expect'(op, ψ0 => circuit)
        loss!(θ) = (dispatch!(circuit, θ); expect(op, copy(ψ0) |> circuit))

        g2 = Vector{Any}(undef, length(θ))
        η = 1e-5
        for i in 1:length(θ)
            θ1 = copy(θ)
            θ2 = copy(θ)
            θ1[i] -= 0.5η
            θ2[i] += 0.5η
            g2[i] = (loss!(θ2) - loss!(θ1))/η |> real
        end
        g3 = faithful_grad(op, ψ0 => circuit)
        g2_ = nbatch(ψ0) == 1 ? g2 : dropdims(sum(hcat(g2...); dims=1), dims=1)

        @test isapprox.(g1, g2_, atol=1e-4) |> all
        @test isapprox.(g2, g3, atol=1e-4) |> all
    end
end

@testset "generator" begin
    @test generator(put(4, 1=>Rx(0.1))) == put(4, 1=>X)
    @test generator(Rx(0.1)) == X
    @test generator(cphase(4, 2, 1, 0.1)) == control(4, 2, 1=>Z)
end

@testset "numdiff & opdiff" begin
    @test collect_blocks(XGate, chain([X, Y, Z])) == [X]

    c = chain(put(4, 1=>Rx(0.5)))
    nd = numdiff(c) do cc
        res = expect(put(4, 1=>Z), zero_state(4) |> cc) |> real
        return res
    end

    ed = faithful_grad(put(4, 1=>Z), zero_state(4) => c)
    @test isapprox(nd, ed, atol=1e-4)

    reg = rand_state(4)
    c = chain(put(4, 1=>Rx(0.5)), control(4, 1, 2=>Ry(0.5)), control(4, 1, 2=>shift(0.3)),  kron(4, 2=>Rz(0.3), 3=>Rx(0.7)))
    loss1z(c) = expect(kron(4, 1=>Z, 2=>X), copy(reg) |> c) |> real  # return loss please
    nd = numdiff(loss1z, c)
    ed = faithful_grad(kron(4, 1=>Z, 2=>X), copy(reg) => c)
    @test isapprox(nd, ed, atol=1e-4)

    # the batched version
    reg = rand_state(4, nbatch=10)
    ed2 = faithful_grad(kron(4, 1=>Z, 2=>X), copy(reg) => c)
    nd2 = numdiff(loss1z, c)
    @test isapprox(nd, ed, atol=1e-4)
end

@testset "random diff circuit" begin
    c = variational_circuit(4, 3, [1=>3, 2=>4, 2=>3, 4=>1])
    rots = collect_blocks(RotationGate, c)
    @test length(rots) == nparameters(c) == 40
    @test dispatch!(+, c, ones(40)*0.1) |> parameters == ones(40)*0.1
    @test dispatch!(+, c, :random) |> parameters != ones(40)*0.1

    nbit = 4
    c = variational_circuit(nbit, 1, pair_ring(nbit), mode=:Split)
    reg = rand_state(4)
    dispatch!(c, randn(nparameters(c)))

    dbs = collect_blocks(Diff, c)
    op = kron(4, 1=>Z, 2=>X)
    loss1z(c) = expect(op, copy(reg) |> c)  # return loss please

    # back propagation
    _, bd = expect'(op, reg=>c)

    # get num gradient
    nd = numdiff(loss1z, c)
    ed = faithful_grad(op, copy(reg)=>c)

    @test isapprox.(nd, ed, atol=1e-4) |> all
    @test isapprox.(nd, bd, atol=1e-4) |> all
end


text:
"""
Java的手法,Julia的实现
"""
mutable struct Node
    left::Node
    right::Node
    weight::Int
    # byte
    data::UInt8

    function Node(data::UInt8, weight::Int)
        Node.data = data
        Node.weight = weight
    end

    function toString()
        return "Node{'
            weight = $weight
            , data = $data
            '}"
    end
end

max_stars_repo_name: Colcactus/Huffman

max_issues_repo_name: Colcactus/Huffman

max_forks_repo_name: Colcactus/Huffman


text:
module ApproxFun
using Base, Reexport,
		AbstractFFTs, FFTW, DualNumbers, FastTransforms,
		LinearAlgebra, RecipesBase, DomainSets, SpecialFunctions #, Arpack

import Calculus

@reexport using ApproxFunBase
@reexport using ApproxFunFourier
@reexport using ApproxFunOrthogonalPolynomials
@reexport using ApproxFunSingularities

import ApproxFunBase: normalize!, flipsign, FiniteRange, Fun, MatrixFun, UnsetSpace, VFun, RowVector,
                    UnivariateSpace, AmbiguousSpace, SumSpace, SubSpace, WeightSpace, NoSpace, Space,
                    HeavisideSpace, PointSpace,
                    IntervalOrSegment, RaggedMatrix, AlmostBandedMatrix,
                    AnyDomain, ZeroSpace, ArraySpace, TrivialInterlacer, BlockInterlacer,
                    AbstractTransformPlan, TransformPlan, ITransformPlan,
                    ConcreteConversion, ConcreteMultiplication, ConcreteDerivative, ConcreteIntegral, CalculusOperator,
                    ConcreteVolterra, Volterra, VolterraWrapper,
                    MultiplicationWrapper, ConversionWrapper, DerivativeWrapper, Evaluation, EvaluationWrapper,
                    Conversion, defaultConversion, defaultcoefficients, default_Fun, Multiplication, Derivative, Integral, bandwidths,
                    ConcreteEvaluation, ConcreteDefiniteLineIntegral, ConcreteDefiniteIntegral, ConcreteIntegral,
                    DefiniteLineIntegral, DefiniteIntegral, ConcreteDefiniteIntegral, ConcreteDefiniteLineIntegral, IntegralWrapper,
                    ReverseOrientation, ReverseOrientationWrapper, ReverseWrapper, Reverse, NegateEven,
                    Dirichlet, ConcreteDirichlet, DirichletWrapper,
                    TridiagonalOperator, SubOperator, Space, @containsconstants, spacescompatible,
                    hasfasttransform, canonicalspace, domain, setdomain, prectype, domainscompatible,
                    plan_transform, plan_itransform, plan_transform!, plan_itransform!, transform, itransform, hasfasttransform,
                    CanonicalTransformPlan, ICanonicalTransformPlan,
                    Integral,
                    domainspace, rangespace, boundary,
                    union_rule, conversion_rule, maxspace_rule, conversion_type, maxspace, hasconversion, points,
                    rdirichlet, ldirichlet, lneumann, rneumann, ivp, bvp,
                    linesum, differentiate, integrate, linebilinearform, bilinearform,
                    UnsetNumber, coefficienttimes, subspace_coefficients, sumspacecoefficients, specialfunctionnormalizationpoint,
                    Segment, IntervalOrSegmentDomain, PiecewiseSegment, isambiguous, Vec, eps, isperiodic,
                    arclength, complexlength,
                    invfromcanonicalD, fromcanonical, tocanonical, fromcanonicalD, tocanonicalD, canonicaldomain, setcanonicaldomain, mappoint,
                    reverseorientation, checkpoints, evaluate, mul_coefficients, coefficients, coefficientmatrix, isconvertible,
                    clenshaw, ClenshawPlan, sineshaw,
                    toeplitz_getindex, toeplitz_axpy!, sym_toeplitz_axpy!, hankel_axpy!, ToeplitzOperator, SymToeplitzOperator, hankel_getindex,
                    SpaceOperator, ZeroOperator, InterlaceOperator,
                    interlace!, reverseeven!, negateeven!, cfstype, pad!, alternatesign!, mobius,
                    extremal_args, hesseneigvals, chebyshev_clenshaw, recA, recB, recC, roots,splitatroots,
                    chebmult_getindex, intpow, alternatingsum,
                    domaintype, diagindshift, rangetype, weight, isapproxinteger, default_Dirichlet, scal!, dotu,
                    components, promoterangespace, promotedomainspace, choosedomainspace,
                    block, blockstart, blockstop, blocklengths, isblockbanded, pointscompatible,
                    AbstractProductSpace, MultivariateFun, BivariateSpace,
                    @wrapperstructure, @wrapperspaces, @wrapper, @calculus_operator, resizedata!, slnorm,
                    sample, chop!, isbanded, colrange, bandwidth

import ApproxFunOrthogonalPolynomials: order

import AbstractFFTs: Plan, fft, ifft
import FFTW: plan_r2r!, fftwNumber, REDFT10, REDFT01, REDFT00, RODFT00, R2HC, HC2R,
                r2r!, r2r,  plan_fft, plan_ifft, plan_ifft!, plan_fft!

import Base: values, convert, getindex, setindex!, *, +, -, ==, <, <=, >, |, !, !=, eltype, iterate,
                >=, /, ^, \, ∪, transpose, size, tail, broadcast, broadcast!, copyto!, copy, to_index, (:),
                similar, map, vcat, hcat, hvcat, show, summary, stride, sum, cumsum, sign, imag, conj, inv,
                complex, reverse, exp, sqrt, abs, abs2, sign, issubset, values, in, first, last, rand, intersect, setdiff,
                isless, union, angle, join, isnan, isapprox, isempty, sort, merge, promote_rule,
                minimum, maximum, extrema, argmax, argmin, findmax, findmin, isfinite,
                zeros, zero, one, promote_rule, repeat, length, resize!, isinf,
                getproperty, findfirst, unsafe_getindex, fld, cld, div, real, imag,
                @_inline_meta, eachindex, firstindex, lastindex, keys, isreal, OneTo,
                Array, Vector, Matrix, view, ones, @propagate_inbounds, print_array,
                split

import Base.Broadcast: BroadcastStyle, Broadcasted, AbstractArrayStyle, broadcastable,
                        DefaultArrayStyle, broadcasted

import LinearAlgebra: BlasInt, BlasFloat, norm, ldiv!, mul!, det, eigvals, dot, cross,
                        qr, qr!, rank, isdiag, istril, istriu, issymmetric, ishermitian,
                        Tridiagonal, diagm, diagm_container, factorize, nullspace,
                        Hermitian, Symmetric, adjoint, transpose, char_uplo

# import Arpack: eigs

import FastTransforms: ChebyshevTransformPlan, IChebyshevTransformPlan, plan_chebyshevtransform,
                        plan_chebyshevtransform!, plan_ichebyshevtransform, plan_ichebyshevtransform!

"""
`Curve` Represents a domain defined by the image of a Fun.  Example
usage would be

```julia
x=Fun(1..2)
Curve(exp(im*x))  # represents an arc
```
"""
const Curve{S,T} = Union{IntervalCurve{S,T},PeriodicCurve{S,T}}
Curve(f::Fun{<:Space{<:PeriodicDomain}}) = PeriodicCurve(f)

#TODO: Make type stable
Curve(f::Fun{<:Space{<:ChebyshevInterval}}) = IntervalCurve(f)

export Curve

import Base: view

##Testing
export bisectioninv

## Further extra features

include("Extras/Extras.jl")
include("Plot/Plot.jl")
include("docs.jl")

#include("precompile.jl")
#_precompile_()

end #module


text:
using ModelingToolkit, OrdinaryDiffEq, Test

@parameters t σ ρ β
@variables x(t) y(t) z(t) a(t) u(t) F(t)
@derivatives D'~t

test_equal(a, b) = @test isequal(simplify(a, polynorm=true), simplify(b, polynorm=true))

eqs = [D(x) ~ σ*(y-x),
       D(y) ~ x*(ρ-z)-y,
       D(z) ~ a*y - β*z,
       0 ~ x - a]

lorenz1 = ODESystem(eqs,t,[x,y,z,a],[σ,ρ,β],name=:lorenz1)

lorenz1_aliased = alias_elimination(lorenz1)
@test length(equations(lorenz1_aliased)) == 3
@test length(states(lorenz1_aliased)) == 3

eqs = [D(x) ~ σ*(y-x),
       D(y) ~ x*(ρ-z)-y,
       D(z) ~ x*y - β*z]

@test lorenz1_aliased == ODESystem(eqs,t,[x,y,z],[σ,ρ,β],observed=[a ~ x],name=:lorenz1)

# Multi-System Reduction

eqs1 = [D(x) ~ σ*(y-x) + F,
       D(y) ~ x*(ρ-z)-u,
       D(z) ~ x*y - β*z]

aliases = [u ~ x + y - z]

lorenz1 = ODESystem(eqs1,pins=[F],observed=aliases,name=:lorenz1)

eqs2 = [D(x) ~ F,
       D(y) ~ x*(ρ-z)-x,
       D(z) ~ x*y - β*z]

aliases2 = [u ~ x - y - z]

lorenz2 = ODESystem(eqs2,pins=[F],observed=aliases2,name=:lorenz2)

connections = [lorenz1.F ~ lorenz2.u,
               lorenz2.F ~ lorenz1.u]

connected = ODESystem([0 ~ a + lorenz1.x - lorenz2.y],t,[a],[],observed=connections,systems=[lorenz1,lorenz2])

# Reduced Unflattened System
#=

connections2 = [lorenz1.F ~ lorenz2.u,
                lorenz2.F ~ lorenz1.u,
                a ~ -lorenz1.x + lorenz2.y]
connected = ODESystem(Equation[],t,[],[],observed=connections2,systems=[lorenz1,lorenz2])
=#

# Reduced Flattened System

flattened_system = ModelingToolkit.flatten(connected)

aliased_flattened_system = alias_elimination(flattened_system)

@test isequal(states(aliased_flattened_system), [
        lorenz1.x
        lorenz1.y
        lorenz1.z
        lorenz2.x
        lorenz2.y
        lorenz2.z
       ])

@test setdiff(parameters(aliased_flattened_system), [
        lorenz1.σ
        lorenz1.ρ
        lorenz1.β
        lorenz1.F
        lorenz2.F
        lorenz2.ρ
        lorenz2.β
       ]) |> isempty

test_equal.(equations(aliased_flattened_system), [
       D(lorenz1.x) ~ lorenz1.σ*(lorenz1.y-lorenz1.x) + lorenz2.x - lorenz2.y - lorenz2.z,
       D(lorenz1.y) ~ lorenz1.x*(lorenz1.ρ-lorenz1.z)-(lorenz1.x + lorenz1.y - lorenz1.z),
       D(lorenz1.z) ~ lorenz1.x*lorenz1.y - lorenz1.β*lorenz1.z,
       D(lorenz2.x) ~ lorenz1.x + lorenz1.y - lorenz1.z,
       D(lorenz2.y) ~ lorenz2.x*(lorenz2.ρ-lorenz2.z)-lorenz2.x,
       D(lorenz2.z) ~ lorenz2.x*lorenz2.y - lorenz2.β*lorenz2.z])

test_equal.(observed(aliased_flattened_system), [
       lorenz1.F ~ lorenz2.x + -1 * (lorenz2.y + lorenz2.z),
       lorenz1.u ~ lorenz1.x + lorenz1.y + -1 * lorenz1.z,
       lorenz2.F ~ lorenz1.x + lorenz1.y + -1 * lorenz1.z,
       a ~ lorenz2.y + -1 * lorenz1.x,
       lorenz2.u ~ lorenz2.x + -1 * (lorenz2.y + lorenz2.z),
])

# issue #578

let
    @variables t x(t) y(t) z(t);
    @derivatives D'~t;
    eqs = [
        D(x) ~ x + y
        x ~ y
    ];
    sys = ODESystem(eqs, t, [x], []);
    asys = alias_elimination(ModelingToolkit.flatten(sys))

    test_equal.(asys.eqs, [D(x) ~ 2x])
    test_equal.(asys.observed, [y ~ x])
end


text:
abstract type AbstractStarAlgebra{O,T} end

struct StarAlgebra{O,T,M<:MultiplicativeStructure,B<:AbstractBasis{T}} <:
       AbstractStarAlgebra{O,T}
    object::O
    mstructure::M
    basis::B

    function StarAlgebra(obj, basis::AbstractBasis, mstr::MultiplicativeStructure)
        O = typeof(obj)
        T = eltype(basis)
        M = typeof(mstr)
        B = typeof(basis)

        return new{O,T,M,B}(obj, mstr, basis)
    end

    function StarAlgebra(obj, mstr::MultiplicativeStructure)
        O = typeof(obj)
        T = Symbol
        M = typeof(mstr)
        B = Basis{T,Int}

        return new{O,T,M,B}(obj, mstr)
    end
end

# TrivialMStructure:
StarAlgebra(obj, basis::AbstractBasis) = StarAlgebra{false}(obj, basis)

function StarAlgebra{Tw}(obj, basis::AbstractBasis) where {Tw}
    mstr = TrivialMStructure{Tw}(basis)
    return StarAlgebra(obj, basis, mstr)
end

# CachedMStructure:
StarAlgebra(obj, basis::AbstractBasis, cache_size::Tuple{<:Integer,Integer}; precompute=false) =
    StarAlgebra{false}(obj, basis, cache_size, precompute=precompute)

function StarAlgebra{Tw}(
    obj,
    basis::AbstractBasis,
    cache_size::Tuple{<:Integer,Integer};
    precompute=false
) where {Tw}
    mstr = CachedMTable{Tw}(basis, table_size = cache_size)
    precompute && complete!(mstr)
    return StarAlgebra(obj, basis, mstr)
end

hasbasis(A::StarAlgebra) = isdefined(A, :basis)

basis(A::StarAlgebra) = A.basis
object(A::StarAlgebra) = A.object
# Base.eltype(A::StarAlgebra{O,B}) where {O,B} = eltype(B)

struct AlgebraElement{A,T,V<:AbstractVector{T}}
    coeffs::V
    parent::A

    function AlgebraElement(coeffs::AbstractVector, A::AbstractStarAlgebra)
        if hasbasis(A)
            @assert length(coeffs) == length(basis(A))
        end
        return new{typeof(A),eltype(coeffs),typeof(coeffs)}(coeffs, A)
    end
end

coeffs(a::AlgebraElement) = a.coeffs
Base.parent(a::AlgebraElement) = a.parent
Base.eltype(a::AlgebraElement) = eltype(coeffs(a))

### constructing elements

function Base.zero(A::AbstractStarAlgebra, T = Int)
    if hasbasis(A)
        I = SparseArrays.indtype(basis(A))
        return AlgebraElement(sparsevec(I[], T[], length(basis(A))), A)
    end
    throw(
        "Algebra without basis; to construct zero use the `AlgebraElement` constructor directly.",
    )
end

function Base.one(A::AbstractStarAlgebra, T = Int)
    hasbasis(A) && return A(one(object(A)), T)
    throw(
        "Algebra without basis; to construct one use the `AlgebraElement` constructor directly.",
    )
end

Base.zero(a::AlgebraElement) = (b = similar(a); return zero!(b))
Base.one(a::AlgebraElement) = one(parent(a))
Base.iszero(a::AlgebraElement) = iszero(coeffs(a))

function Base.isone(a::AlgebraElement)
    b = basis(parent(a))
    k = findfirst(!iszero, coeffs(a))
    k === nothing && return false
    isone(a[k]) || return false
    return isone(b[k])
end
let SA = @static VERSION < v"1.3.0" ? :StarAlgebra : :AbstractStarAlgebra
    @eval begin
        function (A::$SA{O,T})(elt::T, S = Int) where {O,T}
            if hasbasis(A)
                b = basis(A)
                i = b[elt]
                return AlgebraElement(sparsevec([i], [one(S)], length(b)), A)
            else
                throw("Algebra without basis: cannot coerce $elt.")
            end
        end

        function (A::$SA)(x::Number)
            g = one(object(A))
            res = A(g, typeof(x))
            res[g] *= x
            return res
        end
    end
end

Base.similar(X::AlgebraElement, ::Type{T} = eltype(X)) where {T} =
    AlgebraElement(similar(coeffs(X), T), parent(X))

function AlgebraElement{T}(X::AlgebraElement) where T
    v = coeffs(X)
    w = similar(v, T)
    w .= v
    return AlgebraElement(w, parent(X))
end


text:
using Pkg
Pkg.activate("..")
using TimerOutputs, Statistics, ArgParse, CUDAnative, Distributions
include("peps_util.jl")

s = ArgParseSettings()
@add_arg_table! s begin
    "--prefix"
        help = "Prefix string to store magnetization data"
        arg_type = String
    "--do_mag"
        help = "Measure magnetization"
        action = :store_true
    "--Nx"
        help = "Number of columns in the lattice"
        arg_type = Int
    "--Ny"
        help = "Number of rows in the lattice"
        arg_type = Int
    "--D"
        help = "PEPS bond dimension"
        arg_type = Int
        default = 3
    "--chi"
        help = "Environment bond dimension"
        arg_type = Int
        default = 3
    "--sweep_count"
        help = "Total number of sweeps to perform"
        arg_type = Int
        default = 50
    "--simple_update_cutoff"
        help = "Number of sweeps to perform simple updates"
        arg_type = Int
        default = 3
    "--device"
        help = "Which GPU to use"
        arg_type = Int
        default = 1
    "--model"
        help = "Which Hamiltonian to use"
        arg_type = String
        default = "XXZ"
    "--J"
        help = "J-parameter"
        arg_type = Float64
        default = 1.0
    "--hz"
        help = "Field in the z parameter, or random # bound"
        arg_type = Float64
        default = 0.0
    "--hx"
        help = "Field in the x parameter, or random # bound"
        arg_type = Float64
        default = 0.0
    "--random_x"
        help = "Use random fields in the x direction"
        action = :store_true
    "--random_z"
        help = "Use random fields in the z direction"
        action = :store_true
end

# get basic simulation parameters
parsed_args = parse_args(s)
# log file which keeps track of more detailed info about the simulation, not super exciting
io = open(parsed_args["prefix"]*".txt", "w+")
logger = SimpleLogger(io)
global_logger(logger)
Nx = parsed_args["Nx"]
Ny = parsed_args["Ny"]
χ  = parsed_args["chi"]
D  = parsed_args["D"]

parsed_args["device"] != 1 && device!(parsed_args["device"])

# Hamiltonian parameters
J     = parsed_args["J"]
sites = siteinds("S=1/2",Nx*Ny)

# disallow scalar indexing on GPU, which is very slow
CuArrays.allowscalar(false)
A = checkerboardPEPS(sites, Nx, Ny, mindim=D)
# to the user, these appear as normal ITensors, but they have on-device storage
# Julia can detect this at runtime and appropriately dispatch to CUTENSOR
cA = cuPEPS(A)
H  = nothing
if parsed_args["model"] == "XXZ"
    H  = makeCuH_XXZ(Nx, Ny, J)
elseif parsed_args["model"] == "Ising"
    hz = parsed_args["hz"]
    hx = parsed_args["hx"]
    random_z = parsed_args["random_z"]
    random_x = parsed_args["random_x"]
    if random_x || random_z
        seed = abs(rand(Int))
        @show seed
        Random.seed!(seed)
    end
    if random_x
        dx = Uniform(-hx, hx)
        hx = rand(dx, Ny, Nx)
    end
    if random_z
        dz = Uniform(-hz, hz)
        hz = rand(dz, Ny, Nx)
    end
    H  = makeCuH_Ising(Nx, Ny, J, hx, hz)
end
@info "Built cA and H"
# run heaviest functions one time to make Julia compile everything
@info "Built cA and H"
Ls = buildLs(cA, H; mindim=D, maxdim=D)
@info "Built first Ls"
#Rs = buildRs(cA, H; mindim=D, maxdim=D)
Rs = Vector{Environments}(undef, Nx)
@info "Built first Rs"

# actual profiling run
cA, tS, bytes, gctime, memallocs = @timed doSweeps(cA, Ls, Rs, H; mindim=D, maxdim=D, simple_update_cutoff=parsed_args["simple_update_cutoff"], sweep_count=parsed_args["sweep_count"], cutoff=0.0, env_maxdim=χ, do_mag=parsed_args["do_mag"], prefix=parsed_args["prefix"])
println("Done sweeping GPU $tS")
flush(stdout)
flush(io)


text:
export LCO, NMC

## LCO

function LCO(θ, funcs)
    ## parameters section
    # everything here can be modified without regenerating the sol/jacobian.

    # Solid diffusion coefficient [m/s²]
    θ[:D_sp] = 1e-14
    # Electrolyte diffusion coefficient [m/s²]
    θ[:D_p] = 7.5e-10
    # BV reaction rate constant [m^2.5/(mol^1/2⋅s)]
    θ[:k_p] = 2.334e-11
    # MHC reaction, reorganization energy [J] (only needed for MHC reaction)
    θ[:λ_MHC_p] = 6.26e-20
    # Stoichiometry coefficients, θ_min_p > θ_max_p [-]
    θ[:θ_min_p] = 0.99174
    θ[:θ_max_p] = 0.49550
    # Thickness of the electrode [m]
    θ[:l_p] = 80e-6
    # Conductivity [S/m]
    θ[:σ_p] = 100.0
    # Porosity
    θ[:ϵ_p] = 0.385
    # Filler fraction [note: (active material fraction) = (1 - (porosity) - (filler fraction))]
    θ[:ϵ_fp] = 0.025
    # Bruggeman exponent
    θ[:brugg_p] = 4.0
    # Maximum solid interpolate_electrolyte_concentration
    θ[:c_max_p] = 51554.0
    # Solid particle radius
    θ[:Rp_p] = 2e-6

    ## Temperature parameter
    # Thermal conductivity [W/(m⋅K)]
    θ[:λ_p] = 2.1
    # Density [kg/m³]
    θ[:ρ_p] = 2500.0
    # Specific heat capacity [J/(kg⋅K)]
    θ[:Cp_p] = 700.0
    # Activation energy of solid diffusion equation
    θ[:Ea_D_sp] = 5000.0
    # Activation energy of reaction rate equation
    θ[:Ea_k_p] = 5000.0

    ## Custon functions
    # Reaction rate equation
    funcs.rxn_p = rxn_BV
    # Open circuit voltage (OCV or OCP) equation
    funcs.OCV_p = OCV_LCO

    return LiC6, system_LCO_LiC6
end

function LiC6(θ, funcs)

    # Solid diffusion coefficient [m/s²]
    θ[:D_sn] = 3.9e-14
    # Electrolyte diffusion coefficient [m/s²]
    θ[:D_n] = 7.5e-10
    # BV reaction rate constant [m^2.5/(mol^1/2⋅s)]
    θ[:k_n] = 5.0310e-11
    # MHC reaction, reorganization energy [J] (only needed for MHC reaction)
    θ[:λ_MHC_n] = 6.26e-20
    # Stoichiometry coefficients, θ_max_n > θ_min_n [-]
    θ[:θ_max_n] = 0.85510
    θ[:θ_min_n] = 0.01429
    # Thickness of the electrode [m]
    θ[:l_n] = 88e-6
    # Conductivity [S/m]
    θ[:σ_n] = 100.0
    # Porosity
    θ[:ϵ_n] = 0.485
    # Filler fraction [note: (active material fraction) = (1 - (porosity) - (filler fraction))]
    θ[:ϵ_fn] = 0.0326
    # Bruggeman exponent
    θ[:brugg_n] = 4.0
    # Maximum solid interpolate_electrolyte_concentration
    θ[:c_max_n] = 30555.0
    # Solid particle radius
    θ[:Rp_n] = 2e-6

    ## Temperature parameters
    # Thermal conductivity [W/(m⋅K)]
    θ[:λ_n] = 1.7
    # Density [kg/m³]
    θ[:ρ_n] = 2500.0
    # Specific heat capacity [J/(kg⋅K)]
    θ[:Cp_n] = 700.0
    # Activation energy of solid diffusion equation
    θ[:Ea_D_sn] = 5000.0
    # Activation energy of reaction rate equation
    θ[:Ea_k_n] = 5000.0

    ## Aging parameters
    # Initial SEI resistance value [Ω⋅m²]
    θ[:R_SEI] = 0.01
    # Molar weight [kg/mol]
    θ[:M_n] = 7.3e-4
    # Admittance [S/m]
    θ[:k_n_aging] = 1.0
    # Side reaction current density [A/m²]
    θ[:i_0_jside] = 1.5e-6
    # Open circuit voltage for side reaction [V]
    θ[:Uref_s] = 0.4
    # Weigthing factor used in the aging dynamics
    θ[:w] = 2.0

    ## Custon functions
    # Reaction rate equation
    funcs.rxn_n = rxn_BV
    # Open circuit voltage (OCV or OCP) equation
    funcs.OCV_n = OCV_LiC6
end

function system_LCO_LiC6(θ, funcs, cathode, anode;
    # State-of-charge between 0 and 1
    SOC = 1.0,
    ### Cell discretizations, `N` ###
    # Volume discretizations per cathode
    N_p = 10,
    # Volume discretizations per separator
    N_s = 10,
    # Volume discretizations per anode
    N_n = 10,
    # Volume discretizations per positive current collector (temperature only)
    N_a = 10,
    # Volume discretizations per negative current collector (temperature only)
    N_z = 10,
    # Volume discretizations per cathode particle (Fickian diffusion only)
    N_r_p = 10,
    # Volume discretizations per anode particle (Fickian diffusion only)
    N_r_n = 10,
    ### Numerical options, `numerics` ###
    # 1D temperature, true or false
    temperature =  false,
    # (:Fickian) Fickian diffusion, (:quadratic) quadratic approx., (:polynomial) polynomial approx.
    solid_diffusion = :Fickian,
    # if solid_diffusion = :Fickian, then this can either be (:finite_difference) or (:spectral)
    Fickian_method = :finite_difference,
    # (false) off, (:SEI) SEI resistance
    aging =  false,
    # (:symbolic) symbolic Jacobian, (:AD) automatic differenation Jacobian
    # use symbolic when speed is crucial
    jacobian = :symbolic,
    ### User-defined functions in `numerics` ###
    # Effective solid diffusion coefficient function
    D_s_eff = D_s_eff,
    # Reaction rate function
    rxn_rate = rxn_rate,
    # Effective electrolyte diffusion coefficient function
    D_eff = D_eff_linear,
    # Effective electrolyte conductivity function
    K_eff = K_eff,
    # Thermodynamic factor, ∂ln(f)/∂ln(c_e)
    thermodynamic_factor = thermodynamic_factor_linear,
    # By default, this
    ## Custon functions
    # Reaction rate equation will use the reaction defined by the cathode
    rxn_p = funcs.rxn_p,
    # Open circuit voltage (OCV or OCP) equation
    # By default, this will use the OCV defined by the cathode
    OCV_p = funcs.OCV_p,
    # By default, this
    ## Custon functions
    # Reaction rate equation will use the reaction defined by the anode
    rxn_n = funcs.rxn_n,
    # Open circuit voltage (OCV or OCP) equation
    # By default, this will use the OCV defined by the anode
    OCV_n = funcs.OCV_n,
    )

    ## Physical parameters for the system

    # Electrolyte diffusion coefficient [m/s²]
    θ[:D_s] = 7.5e-10

    # Electrode thicknesses [m/s²]
    θ[:l_s] = 25e-6
    θ[:l_a] = 10e-6
    θ[:l_z] = 10e-6

    # Conductivities [S/m]
    θ[:σ_a] = 3.55e7
    θ[:σ_z] = 5.96e7

    # Porosity
    θ[:ϵ_s] = 0.724

    # Bruggeman exponent
    θ[:brugg_s] = 4.0

    # Transference number [-]
    θ[:t₊] = 0.364

    # Initial electrolyte concentration [mol/m³]
    θ[:c_e₀] = 1000.0
    # Initial temperature [K]
    θ[:T₀] = 25 + 273.15
    # Ambient temperature [K]
    θ[:T_amb] = 25 + 273.15

    ## Temperature
    # Thermal conductivities [W/(m⋅K)]
    θ[:λ_s] = 0.16
    θ[:λ_a] = 237.0
    θ[:λ_z] = 401.0

    # Densities [kg/m³]
    θ[:ρ_s] = 1100.0
    θ[:ρ_a] = 2700.0
    θ[:ρ_z] = 8940.0

    # Heat capacities [J/(kg⋅K)]
    θ[:Cp_s] = 700.0
    θ[:Cp_a] = 897.0
    θ[:Cp_z] = 385.0

    # Heat transfer coefficient [W/m²⋅K]
    θ[:h_cell] = 1.0

    ## Options section
    # everything here can be modified freely

    # `NaN` deactivates the bound
    bounds = boundary_stop_conditions()
    # Maximum permitted voltage [V]
    bounds.V_min = 2.5
    # Minimum permitted voltage [V]
    bounds.V_max = 4.3
    # Maximum permitted SOC [-]
    bounds.SOC_min = 0.0
    # Minimum permitted SOC [-]
    bounds.SOC_max = 1.0
    # Maximum permitted temperature [K]
    bounds.T_max = 55 + 273.15
    # Maximum permitted solid surface concentration in the anode [mol/m³]
    bounds.c_s_n_max = NaN
    # Maximum permitted current [C-rate]
    bounds.I_max = NaN
    # Minimum permitted current [C-rate]
    bounds.I_min = NaN
    # Minimum permitted plating overpotential at the separator-anode interface [V]
    bounds.η_plating_min = NaN
    # Minimum permitted electrolyte concentration [mol/m³]
    bounds.c_e_min = NaN

    opts = options_simulation()
    # Initial state of charge for a new simulation between 0 and 1
    opts.SOC = SOC # defined above
    # Saving sol states is expensive. What states do you want to keep? See the output of solution below for more info. Must be a tuple
    opts.outputs = (:t, :V)
    # Absolute tolerance of DAE solver
    opts.abstol = 1e-6
    # Relative tolerance of DAE solver
    opts.reltol = 1e-3
    # Maximum iterations for the DAE solver
    opts.maxiters = 10_000
    # Flag to check the bounds during simulation (SOC max/min, V max/min, etc.)
    opts.check_bounds = true
    # Get a new initial guess for DAE initialization
    opts.reinit = true
    # Show some outputs during simulation runtime
    opts.verbose = false
    # Interpolate the final results to match the exact simulation end point
    opts.interp_final = true
    # Times when the DAE solver explicitly stops
    opts.tstops = Float64[]
    # For input functions, times when there is a known discontinuity. Unknown discontinuities are handled automatically but less efficiently
    opts.tdiscon = Float64[]
    # :interpolate or :extrapolate when interpolating the solution
    opts.interp_bc = :interpolate

    #### DO NOT MODIFY BELOW ###
    N = discretizations_per_section(N_p, N_s, N_n, N_a, N_z, N_r_p, N_r_n)
    numerics = options_numerical(temperature, solid_diffusion, Fickian_method, aging, cathode, anode, rxn_p, rxn_n, OCV_p, OCV_n, D_s_eff, rxn_rate, D_eff, K_eff, thermodynamic_factor, jacobian)

    return θ, bounds, opts, N, numerics
end

## NMC

function NMC(θ, funcs)
    ## parameters section
    # everything here can be modified without regenerating the sol/jacobian.

    # Solid diffusion coefficient [m/s²]
    θ[:D_sp] = 2e-14
    # BV reaction rate constant [m^2.5/(mol^1/2⋅s)]
    θ[:k_p] = 6.3066e-10
    # Stoichiometry coefficients, θ_min_p > θ_max_p [-]
    θ[:θ_min_p] = 0.955473
    θ[:θ_max_p] = 0.359749
    # Thickness of the electrode [m]
    θ[:l_p] = 41.6e-6
    # Conductivity [S/m]
    θ[:σ_p] = 100
    # Porosity
    θ[:ϵ_p] = 0.3
    # Filler fraction [note: (active material fraction) = (1 - (porosity) - (filler fraction))]
    θ[:ϵ_fp] = 0.12
    # Bruggeman exponent
    θ[:brugg_p] = 1.5
    # Maximum solid interpolate_electrolyte_concentration
    θ[:c_max_p] = 51830.0
    # Solid particle radius
    θ[:Rp_p] = 7.5e-6
    # Activation energy of solid diffusion equation
    θ[:Ea_D_sp] = 2.5e4
    # Activation energy of reaction rate equation
    θ[:Ea_k_p] = 3e4

    ## Custon functions
    # Reaction rate equation
    funcs.rxn_p = rxn_BV
    # Open circuit voltage (OCV or OCP) equation
    funcs.OCV_p = OCV_NMC

    return LiC6_NMC, system_NMC_LiC6
end

function LiC6_NMC(θ, funcs)

    # Solid diffusion coefficient [m/s²]
    θ[:D_sn] = 1.5e-14
    # BV reaction rate constant [m^2.5/(mol^1/2⋅s)]
    θ[:k_n] = 6.3466e-10
    # Stoichiometry coefficients, θ_max_n > θ_min_n [-]
    θ[:θ_max_n] = 0.790813
    θ[:θ_min_n] = 0.001
    # Thickness of the electrode [m]
    θ[:l_n] = 48e-6
    # Conductivity [S/m]
    θ[:σ_n] = 100
    # Porosity
    θ[:ϵ_n] = 0.3
    # Filler fraction [note: (active material fraction) = (1 - (porosity) - (filler fraction))]
    θ[:ϵ_fn] = 0.038
    # Bruggeman exponent
    θ[:brugg_n] = 1.5
    # Maximum solid interpolate_electrolyte_concentration
    θ[:c_max_n] = 31080.0
    # Solid particle radius
    θ[:Rp_n] = 10e-6
    # Activation energy of solid diffusion equation
    θ[:Ea_D_sn] = 4e4
    # Activation energy of reaction rate equation
    θ[:Ea_k_n] = 3e4

    ## Custon functions
    # Reaction rate equation
    funcs.rxn_n = rxn_BV
    # Open circuit voltage (OCV or OCP) equation
    funcs.OCV_n = OCV_LiC6_with_NMC
end

function system_NMC_LiC6(θ, funcs, cathode, anode;
    # State-of-charge between 0 and 1
    SOC = 1.0,
    ### Cell discretizations, `N` ###
    # Volume discretizations per cathode
    N_p = 10,
    # Volume discretizations per separator
    N_s = 10,
    # Volume discretizations per anode
    N_n = 10,
    # Volume discretizations per positive current collector (temperature only)
    N_a = 10,
    # Volume discretizations per negative current collector (temperature only)
    N_z = 10,
    # Volume discretizations per cathode particle (Fickian diffusion only)
    N_r_p = 10,
    # Volume discretizations per anode particle (Fickian diffusion only)
    N_r_n = 10,
    ### Numerical options, `numerics` ###
    # 1D temperature, true or false
    temperature =  false,
    # (:Fickian) Fickian diffusion, (:quadratic) quadratic approx., (:polynomial) polynomial approx.
    solid_diffusion = :Fickian,
    # if solid_diffusion = :Fickian, then this can either be (:finite_difference) or (:spectral)
    Fickian_method = :finite_difference,
    # (false) off, (:SEI) SEI resistance
    aging =  false,
    # (:symbolic) symbolic Jacobian, (:AD) automatic differenation Jacobian
    # use symbolic when speed is crucial
    jacobian = :symbolic,
    ### User-defined functions in `numerics` ###
    # Effective solid diffusion coefficient function
    D_s_eff = D_s_eff,
    # Reaction rate function
    rxn_rate = rxn_rate,
    # Effective electrolyte diffusion coefficient function
    D_eff = D_eff,
    # Effective electrolyte conductivity function
    K_eff = K_eff,
    # Thermodynamic factor, ∂ln(f)/∂ln(c_e)
    thermodynamic_factor = thermodynamic_factor_linear,
    # By default, this
    ## Custon functions
    # Reaction rate equation will use the reaction defined by the cathode
    rxn_p = funcs.rxn_p,
    # Open circuit voltage (OCV or OCP) equation
    # By default, this will use the OCV defined by the cathode
    OCV_p = funcs.OCV_p,
    # By default, this
    ## Custon functions
    # Reaction rate equation will use the reaction defined by the anode
    rxn_n = funcs.rxn_n,
    # Open circuit voltage (OCV or OCP) equation
    # By default, this will use the OCV defined by the anode
    OCV_n = funcs.OCV_n,
    )

    # Electrode thicknesses [m/s²]## Physical parameters for the system
    θ[:l_s] = 25e-6

    # Porosity
    θ[:ϵ_s] = 0.4

    # Thermal conductivities [W/(m⋅K)]
    θ[:brugg_s] = 1.5

    # Transference number [-]
    θ[:t₊] = 0.38

    # Initial electrolyte concentration [mol/m³]
    θ[:c_e₀] = 1200
    # Initial temperature [K]
    θ[:T₀] = 25 + 273.15
    # Ambient temperature [K]
    θ[:T_amb] = 25 + 273.15

    ## Options section
    # everything here can be modified freely

    # `NaN` deactivates the bound
    bounds = boundary_stop_conditions()
    # Maximum permitted voltage [V]
    bounds.V_min = 2.8
    # Minimum permitted voltage [V]
    bounds.V_max = 4.2
    # Maximum permitted SOC [-]
    bounds.SOC_min = 0.0
    # Minimum permitted SOC [-]
    bounds.SOC_max = 1.0
    # Maximum permitted temperature [K]
    bounds.T_max = NaN
    # Maximum permitted solid surface concentration in the anode [mol/m³]
    bounds.c_s_n_max = NaN
    # Maximum permitted current [C-rate]
    bounds.I_max = NaN
    # Minimum permitted current [C-rate]
    bounds.I_min = NaN
    # Minimum permitted plating overpotential at the separator-anode interface [V]
    bounds.η_plating_min = NaN
    # Minimum permitted electrolyte concentration [mol/m³]
    bounds.c_e_min = NaN

    opts = options_simulation()
    # Initial state of charge for a new simulation between 0 and 1
    opts.SOC = SOC # defined above
    # Saving sol states is expensive. What states do you want to keep? See the output of sol below for more info. Must be a tuple
    opts.outputs = (:t, :V)
    # Absolute tolerance of DAE solver
    opts.abstol = 1e-6
    # Relative tolerance of DAE solver
    opts.reltol = 1e-3
    # Maximum iterations for the DAE solver
    opts.maxiters = 10_000
    # Flag to check the bounds during simulation (SOC max/min, V max/min, etc.)
    opts.check_bounds = true
    # Get a new initial guess for DAE initialization
    opts.reinit = true
    # Show some outputs during simulation runtime
    opts.verbose = false
    # Interpolate the final results to match the exact simulation end point
    opts.interp_final = true
    # Times when the DAE solver explicitly stops
    opts.tstops = Float64[]
    # For input functions, times when there is a known discontinuity. Unknown discontinuities are handled automatically but less efficiently
    opts.tdiscon = Float64[]
    # :interpolate or :extrapolate when interpolating the sol
    opts.interp_bc = :interpolate

    #### DO NOT MODIFY BELOW ###
    N = discretizations_per_section(N_p, N_s, N_n, N_a, N_z, N_r_p, N_r_n)
    numerics = options_numerical(temperature, solid_diffusion, Fickian_method, aging, cathode, anode, rxn_p, rxn_n, OCV_p, OCV_n, D_s_eff, rxn_rate, D_eff, K_eff, thermodynamic_factor, jacobian)

    return θ, bounds, opts, N, numerics
end


text:
module MortalityTables

include("Mortality.jl")

export MortalityTable,
    q, p, qx,
    TableMetaData
end # module