text:
"""
    known_first(::Type{T})

If `first` of an instance of type `T` is known at compile time, return it.
Otherwise, return `nothing`.

@test isnothing(known_first(typeof(1:4)))
@test isone(known_first(typeof(Base.OneTo(4))))
"""
known_first(x) = known_first(typeof(x))
function known_first(::Type{T}) where {T}
    if parent_type(T) <: T
        return nothing
    else
        return known_first(parent_type(T))
    end
end
known_first(::Type{Base.OneTo{T}}) where {T} = one(T)

"""
    known_last(::Type{T})

If `last` of an instance of type `T` is known at compile time, return it.
Otherwise, return `nothing`.

@test isnothing(known_last(typeof(1:4)))
using StaticArrays
@test known_last(typeof(SOneTo(4))) == 4
"""
known_last(x) = known_last(typeof(x))
function known_last(::Type{T}) where {T}
    if parent_type(T) <: T
        return nothing
    else
        return known_last(parent_type(T))
    end
end

"""
    known_step(::Type{T})

If `step` of an instance of type `T` is known at compile time, return it.
Otherwise, return `nothing`.

@test isnothing(known_step(typeof(1:0.2:4)))
@test isone(known_step(typeof(1:4)))
"""
known_step(x) = known_step(typeof(x))
function known_step(::Type{T}) where {T}
    if parent_type(T) <: T
        return nothing
    else
        return known_step(parent_type(T))
    end
end
known_step(::Type{<:AbstractUnitRange{T}}) where {T} = one(T)

"""
    OptionallyStaticUnitRange(start, stop) <: AbstractUnitRange{Int}

This range permits diverse representations of arrays to communicate common information
about their indices. Each field may be an integer or `Val(<:Integer)` if it is known
at compile time. An `OptionallyStaticUnitRange` is intended to be constructed internally
from other valid indices. Therefore, users should not expect the same checks are used
to ensure construction of a valid `OptionallyStaticUnitRange` as a `UnitRange`.
"""
struct OptionallyStaticUnitRange{F<:Integer,L<:Integer} <: AbstractUnitRange{Int}
    start::F
    stop::L

    function OptionallyStaticUnitRange(start, stop)
        if eltype(start) <: Int
            if eltype(stop) <: Int
                return new{typeof(start),typeof(stop)}(start, stop)
            else
                return OptionallyStaticUnitRange(start, Int(stop))
            end
        else
            return OptionallyStaticUnitRange(Int(start), stop)
        end
    end

    function OptionallyStaticUnitRange{F,L}(x::AbstractRange) where {F,L}
        if step(x) == 1
            return OptionallyStaticUnitRange(static_first(x), static_last(x))
        else
            throw(ArgumentError("step must be 1, got $(step(x))"))
        end
    end

    function OptionallyStaticUnitRange(x::AbstractRange)
        if step(x) == 1
            return OptionallyStaticUnitRange(static_first(x), static_last(x))
        else
            throw(ArgumentError("step must be 1, got $(step(x))"))
        end
    end
end

function Base.first(r::OptionallyStaticUnitRange)::Int
    if known_first(r) === nothing
        return r.start
    else
        return known_first(r)
    end
end
function Base.last(r::OptionallyStaticUnitRange)::Int
    if known_last(r) === nothing
        return r.stop
    else
        return known_last(r)
    end
end

known_first(::Type{<:OptionallyStaticUnitRange{StaticInt{F}}}) where {F} = F
known_step(::Type{<:OptionallyStaticUnitRange}) = 1
known_last(::Type{<:OptionallyStaticUnitRange{<:Any,StaticInt{L}}}) where {L} = L

"""
    OptionallyStaticStepRange(start, step, stop) <: OrdinalRange{Int,Int}

Similarly to [`OptionallyStaticUnitRange`](@ref), `OptionallyStaticStepRange` permits
a combination of static and standard primitive `Int`s to construct a range. It
specifically enables the use of ranges without a step size of 1. It may be constructed
through the use of `OptionallyStaticStepRange` directly or using static integers with
the range operator (i.e., `:`).

```julia
julia> using ArrayInterface

julia> x = ArrayInterface.StaticInt(2);

julia> x:x:10
ArrayInterface.StaticInt{2}():ArrayInterface.StaticInt{2}():10

julia> ArrayInterface.OptionallyStaticStepRange(x, x, 10)
ArrayInterface.StaticInt{2}():ArrayInterface.StaticInt{2}():10
```
"""
struct OptionallyStaticStepRange{F<:Integer,S<:Integer,L<:Integer} <: OrdinalRange{Int,Int}
    start::F
    step::S
    stop::L

    function OptionallyStaticStepRange(start, step, stop)
        if eltype(start) <: Int
            if eltype(stop) <: Int
                lst = _steprange_last(start, step, stop)
                return new{typeof(start),typeof(step),typeof(lst)}(start, step, lst)
            else
                return OptionallyStaticStepRange(start, step, Int(stop))
            end
        else
            return OptionallyStaticStepRange(Int(start), step, stop)
        end
    end

    function OptionallyStaticStepRange(x::AbstractRange)
        return OptionallyStaticStepRange(static_first(x), static_step(x), static_last(x))
    end
end

# to make StepRange constructor inlineable, so optimizer can see `step` value
@inline function _steprange_last(start::StaticInt, step::StaticInt, stop::StaticInt)
    return StaticInt(_steprange_last(Int(start), Int(step), Int(stop)))
end
@inline function _steprange_last(start::Integer, step::StaticInt, stop::StaticInt)
    if step === one(step)
        # we don't need to check the `stop` if we know it acts like a unit range
        return stop
    else
        return _steprange_last(start, Int(step), Int(stop))
    end
end
@inline function _steprange_last(start::Integer, step::Integer, stop::Integer)
    z = zero(step)
    if step === z
        throw(ArgumentError("step cannot be zero"))
    else
        if stop == start
            return Int(stop)
        else
            if step > z
                if stop > start
                    return stop - Int(unsigned(stop - start) % step)
                else
                    return Int(start - one(start))
                end
            else
                if stop > start
                    return Int(start + one(start))
                else
                    return stop + Int(unsigned(start - stop) % -step)
                end
            end
        end
    end
end
function Base.first(r::OptionallyStaticStepRange)::Int
    if known_first(r) === nothing
        return r.start
    else
        return known_first(r)
    end
end
function Base.step(r::OptionallyStaticStepRange)::Int
    if known_step(r) === nothing
        return r.step
    else
        return known_step(r)
    end
end
function Base.last(r::OptionallyStaticStepRange)::Int
    if known_last(r) === nothing
        return r.stop
    else
        return known_last(r)
    end
end

known_first(::Type{<:OptionallyStaticStepRange{StaticInt{F}}}) where {F} = F
known_step(::Type{<:OptionallyStaticStepRange{<:Any,StaticInt{S}}}) where {S} = S
known_last(::Type{<:OptionallyStaticStepRange{<:Any,<:Any,StaticInt{L}}}) where {L} = L

Base.:(:)(L::Integer, ::StaticInt{U}) where {U} = OptionallyStaticUnitRange(L, StaticInt(U))
Base.:(:)(::StaticInt{L}, U::Integer) where {L} = OptionallyStaticUnitRange(StaticInt(L), U)
function Base.:(:)(::StaticInt{L}, ::StaticInt{U}) where {L,U}
    return OptionallyStaticUnitRange(StaticInt(L), StaticInt(U))
end
function Base.:(:)(::StaticInt{F}, ::StaticInt{S}, ::StaticInt{L}) where {F,S,L}
    return OptionallyStaticStepRange(StaticInt(F), StaticInt(S), StaticInt(L))
end
function Base.:(:)(start::Integer, ::StaticInt{S}, ::StaticInt{L}) where {S,L}
    return OptionallyStaticStepRange(start, StaticInt(S), StaticInt(L))
end
function Base.:(:)(::StaticInt{F}, ::StaticInt{S}, stop::Integer) where {F,S}
    return OptionallyStaticStepRange(StaticInt(F), StaticInt(S), stop)
end
function Base.:(:)(::StaticInt{F}, step::Integer, ::StaticInt{L}) where {F,L}
    return OptionallyStaticStepRange(StaticInt(F), step, StaticInt(L))
end
function Base.:(:)(start::Integer, step::Integer, ::StaticInt{L}) where {L}
    return OptionallyStaticStepRange(start, step, StaticInt(L))
end
function Base.:(:)(start::Integer, ::StaticInt{S}, stop::Integer) where {S}
    return OptionallyStaticStepRange(start, StaticInt(S), stop)
end
function Base.:(:)(::StaticInt{F}, step::Integer, stop::Integer) where {F}
    return OptionallyStaticStepRange(StaticInt(F), step, stop)
end
function Base.:(:)(::StaticInt{F}, ::StaticInt{1}, ::StaticInt{L}) where {F,L}
    return OptionallyStaticUnitRange(StaticInt(F), StaticInt(L))
end
function Base.:(:)(start::Integer, ::StaticInt{1}, ::StaticInt{L}) where {L}
    return OptionallyStaticUnitRange(start, StaticInt(L))
end
function Base.:(:)(::StaticInt{F}, ::StaticInt{1}, stop::Integer) where {F}
    return OptionallyStaticUnitRange(StaticInt(F), stop)
end
function Base.:(:)(start::Integer, ::StaticInt{1}, stop::Integer)
    return OptionallyStaticUnitRange(start, stop)
end

function Base.isempty(r::OptionallyStaticUnitRange)
    if known_first(r) === oneunit(eltype(r))
        return unsafe_isempty_one_to(last(r))
    else
        return unsafe_isempty_unit_range(first(r), last(r))
    end
end

function Base.isempty(r::OptionallyStaticStepRange)
    return (r.start != r.stop) & ((r.step > zero(r.step)) != (r.stop > r.start))
end

unsafe_isempty_one_to(lst) = lst <= zero(lst)
unsafe_isempty_unit_range(fst, lst) = fst > lst

unsafe_length_one_to(lst::Int) = lst
unsafe_length_one_to(::StaticInt{L}) where {L} = L

# TODO this should probably be renamed because the point is that it is safe
@inline function unsafe_length_step_range(start::Int, step::Int, stop::Int)
    if step > 1
        return Base.checked_add(Int(div(unsigned(stop - start), step)), 1)
    elseif step < -1
        return Base.checked_add(Int(div(unsigned(start - stop), -step)), 1)
    elseif step > 0
        return Base.checked_add(Int(div(Base.checked_sub(stop, start), step)), 1)
    else
        return Base.checked_add(Int(div(Base.checked_sub(start, stop), -step)), 1)
    end
end

@propagate_inbounds function Base.getindex(
    r::OptionallyStaticUnitRange,
    s::AbstractUnitRange{<:Integer},
)
    @boundscheck checkbounds(r, s)
    f = static_first(r)
    fnew = f - one(f)
    return (fnew+static_first(s)):(fnew+static_last(s))
end

@propagate_inbounds function Base.getindex(r::OptionallyStaticUnitRange, i::Integer)
    if known_first(r) === oneunit(eltype(r))
        return get_index_one_to(r, i)
    else
        return get_index_unit_range(r, i)
    end
end

@inline function get_index_one_to(r, i)
    @boundscheck if ((i < 1) || (i > last(r)))
        throw(BoundsError(r, i))
    end
    return convert(eltype(r), i)
end

@inline function get_index_unit_range(r, i)
    val = first(r) + (i - 1)
    @boundscheck if (i < 1) || val > last(r)
        throw(BoundsError(r, i))
    end
    return convert(eltype(r), val)
end

@inline _try_static(::StaticInt{N}, ::StaticInt{N}) where {N} = StaticInt{N}()
@inline function _try_static(::StaticInt{M}, ::StaticInt{N}) where {M,N}
    @assert false "Unequal Indices: StaticInt{$M}() != StaticInt{$N}()"
end
@noinline unequal_error(x,y) = @assert false "Unequal Indices: x == $x != $y == y"
@inline function check_equal(x, y)
    x == y || unequal_error(x,y)
end
@propagate_inbounds function _try_static(::StaticInt{N}, x) where {N}
    @boundscheck check_equal(StaticInt{N}(), x)
    return StaticInt{N}()
end
@propagate_inbounds function _try_static(x, ::StaticInt{N}) where {N}
    @boundscheck check_equal(x, StaticInt{N}())
    return StaticInt{N}()
end
@propagate_inbounds function _try_static(x, y)
    @boundscheck check_equal(x, y)
    return x
end

###
### length
###
@inline function known_length(::Type{T}) where {T<:OptionallyStaticUnitRange}
    fst = known_first(T)
    lst = known_last(T)
    if fst === nothing || lst === nothing
        return nothing
    else
        if fst === oneunit(eltype(T))
            return unsafe_length_one_to(lst)
        else
            return unsafe_length_unit_range(fst, lst)
        end
    end
end

@inline function known_length(::Type{T}) where {T<:OptionallyStaticStepRange}
    fst = known_first(T)
    stp = known_step(T)
    lst = known_last(T)
    if fst === nothing || stp === nothing || lst === nothing
        return nothing
    else
        if stp === 1
            if fst === oneunit(eltype(T))
                return unsafe_length_one_to(lst)
            else
                return unsafe_length_unit_range(fst, lst)
            end
        else
            return unsafe_length_step_range(fst, stp, lst)
        end
    end
end

function Base.length(r::OptionallyStaticUnitRange)
    if isempty(r)
        return 0
    else
        if known_first(r) === 1
            return unsafe_length_one_to(last(r))
        else
            return unsafe_length_unit_range(first(r), last(r))
        end
    end
end

function Base.length(r::OptionallyStaticStepRange)
    if isempty(r)
        return 0
    else
        if known_step(r) === 1
            if known_first(r) === 1
                return unsafe_length_one_to(last(r))
            else
                return unsafe_length_unit_range(first(r), last(r))
            end
        else
            return unsafe_length_step_range(Int(first(r)), Int(step(r)), Int(last(r)))
        end
    end
end

unsafe_length_unit_range(start::Integer, stop::Integer) = Int((stop - start) + 1)

function Base.AbstractUnitRange{T}(r::OptionallyStaticUnitRange) where {T}
    if T <: Int
        return r
    else
        if known_first(r) === 1 && T <: Integer
            return OneTo{T}(last(r))
        else
            return UnitRange{T}(first(r), last(r))
        end
    end
end

const OptionallyStaticRange = Union{<:OptionallyStaticUnitRange,<:OptionallyStaticStepRange}

Base.eachindex(r::OptionallyStaticRange) = r
@inline function Base.iterate(r::OptionallyStaticRange)
    isempty(r) && return nothing
    fi = Int(first(r));
    fi, fi
end

Base.to_shape(x::OptionallyStaticRange) = length(x)
Base.to_shape(x::Slice{T}) where {T<:OptionallyStaticRange} = length(x)

@inline function Base.axes(S::Slice{T}) where {T<:OptionallyStaticRange}
    if known_first(T) === 1 && known_step(T) === 1
        return (S.indices,)
    else
        return (Base.IdentityUnitRange(S.indices),)
    end
end

@inline function Base.axes1(S::Slice{T}) where {T<:OptionallyStaticRange}
    if known_first(T) === 1 && known_step(T) === 1
        return S.indices
    else
        return Base.IdentityUnitRange(S.indices)
    end
end

Base.:(-)(r::OptionallyStaticRange) = -static_first(r):-static_step(r):-static_last(r)

Base.reverse(r::OptionallyStaticUnitRange) = static_last(r):static(-1):static_first(r)
function Base.reverse(r::OptionallyStaticStepRange)
    return OptionallyStaticStepRange(static_last(r), -static_step(r), static_first(r))
end

function Base.show(io::IO, r::OptionallyStaticRange)
    print(io, first(r))
    if known_step(r) === 1
        print(io, ":")
    else
        print(io, ":")
        print(io, step(r))
        print(io, ":")
    end
    print(io, last(r))
end

"""
  reduce_tup(f::F, inds::Tuple{Vararg{Any,N}}) where {F,N}

An optimized `reduce` for tuples. `Base.reduce`'s `afoldl` will often not inline.
Additionally, `reduce_tup` attempts to order the reduction in an optimal manner.

```julia
julia> using StaticArrays, ArrayInterface, BenchmarkTools

julia> rsum(v::SVector) = ArrayInterface.reduce_tup(+, v.data)
rsum (generic function with 2 methods)

julia> for n ∈ 2:16
           @show n
           v = @SVector rand(n)
           s1 = @btime  sum(\$(Ref(v))[])
           s2 = @btime rsum(\$(Ref(v))[])
       end
n = 2
  0.863 ns (0 allocations: 0 bytes)
  0.863 ns (0 allocations: 0 bytes)
n = 3
  0.862 ns (0 allocations: 0 bytes)
  0.863 ns (0 allocations: 0 bytes)
n = 4
  0.862 ns (0 allocations: 0 bytes)
  0.862 ns (0 allocations: 0 bytes)
n = 5
  1.074 ns (0 allocations: 0 bytes)
  0.864 ns (0 allocations: 0 bytes)
n = 6
  0.864 ns (0 allocations: 0 bytes)
  0.862 ns (0 allocations: 0 bytes)
n = 7
  1.075 ns (0 allocations: 0 bytes)
  0.864 ns (0 allocations: 0 bytes)
n = 8
  1.077 ns (0 allocations: 0 bytes)
  0.865 ns (0 allocations: 0 bytes)
n = 9
  1.081 ns (0 allocations: 0 bytes)
  0.865 ns (0 allocations: 0 bytes)
n = 10
  1.195 ns (0 allocations: 0 bytes)
  0.867 ns (0 allocations: 0 bytes)
n = 11
  1.357 ns (0 allocations: 0 bytes)
  1.400 ns (0 allocations: 0 bytes)
n = 12
  1.543 ns (0 allocations: 0 bytes)
  1.074 ns (0 allocations: 0 bytes)
n = 13
  1.702 ns (0 allocations: 0 bytes)
  1.077 ns (0 allocations: 0 bytes)
n = 14
  1.913 ns (0 allocations: 0 bytes)
  0.867 ns (0 allocations: 0 bytes)
n = 15
  2.076 ns (0 allocations: 0 bytes)
  1.077 ns (0 allocations: 0 bytes)
n = 16
  2.273 ns (0 allocations: 0 bytes)
  1.078 ns (0 allocations: 0 bytes)
```

More importantly, `reduce_tup(_pick_range, inds)` often performs better than `reduce(_pick_range, inds)`.
```julia
julia> using ArrayInterface, BenchmarkTools

julia> inds = (Base.OneTo(100), 1:100, 1:ArrayInterface.StaticInt(100))
(Base.OneTo(100), 1:100, 1:Static(100))

julia> @btime reduce(ArrayInterface._pick_range, \$(Ref(inds))[])
  6.405 ns (0 allocations: 0 bytes)
Base.Slice(Static(1):Static(100))

julia> @btime ArrayInterface.reduce_tup(ArrayInterface._pick_range, \$(Ref(inds))[])
  2.570 ns (0 allocations: 0 bytes)
Base.Slice(Static(1):Static(100))

julia> inds = (Base.OneTo(100), 1:100, 1:UInt(100))
(Base.OneTo(100), 1:100, 0x0000000000000001:0x0000000000000064)

julia> @btime reduce(ArrayInterface._pick_range, \$(Ref(inds))[])
  6.411 ns (0 allocations: 0 bytes)
Base.Slice(Static(1):100)

julia> @btime ArrayInterface.reduce_tup(ArrayInterface._pick_range, \$(Ref(inds))[])
  2.592 ns (0 allocations: 0 bytes)
Base.Slice(Static(1):100)

julia> inds = (Base.OneTo(100), 1:100, 1:UInt(100), Int32(1):Int32(100))
(Base.OneTo(100), 1:100, 0x0000000000000001:0x0000000000000064, 1:100)

julia> @btime reduce(ArrayInterface._pick_range, \$(Ref(inds))[])
  9.048 ns (0 allocations: 0 bytes)
Base.Slice(Static(1):100)

julia> @btime ArrayInterface.reduce_tup(ArrayInterface._pick_range, \$(Ref(inds))[])
  2.569 ns (0 allocations: 0 bytes)
Base.Slice(Static(1):100)
```
"""
@generated function reduce_tup(f::F, inds::Tuple{Vararg{Any,N}}) where {F,N}
    q = Expr(:block, Expr(:meta, :inline, :propagate_inbounds))
    if N == 1
        push!(q.args, :(inds[1]))
        return q
    end
    syms = Vector{Symbol}(undef, N)
    i = 0
    for n ∈ 1:N
        syms[n] = iₙ = Symbol(:i_, (i += 1))
        push!(q.args, Expr(:(=), iₙ, Expr(:ref, :inds, n)))
    end
    W =  1 << (8sizeof(N) - 2 - leading_zeros(N))
    while W > 0
        _N = length(syms)
        for _ ∈ 2W:W:_N
            for w ∈ 1:W
                new_sym = Symbol(:i_, (i += 1))
                push!(q.args, Expr(:(=), new_sym, Expr(:call, :f, syms[w], syms[w+W])))
                syms[w] = new_sym
            end
            deleteat!(syms, 1+W:2W)
        end
        W >>>= 1
    end
    q
end

@propagate_inbounds function _pick_range(x, y)
    fst = _try_static(static_first(x), static_first(y))
    lst = _try_static(static_last(x), static_last(y))
    return Base.Slice(OptionallyStaticUnitRange(fst, lst))
end

"""
    indices(x, dim)

Given an array `x`, this returns the indices along dimension `dim`.
"""
@inline indices(x, d) = indices(axes(x, d))

"""
    indices(x) -> AbstractUnitRange

Returns valid indices for the entire length of `x`.
"""
@inline function indices(x)
    inds = eachindex(x)
    if inds isa AbstractUnitRange && eltype(inds) <: Integer
        return Base.Slice(OptionallyStaticUnitRange(inds))
    else
        return inds
    end
end
@inline indices(x::AbstractUnitRange{<:Integer}) = Base.Slice(OptionallyStaticUnitRange(x))

"""
    indices(x::Tuple) -> AbstractUnitRange

Returns valid indices for the entire length of each array in `x`.
"""
@propagate_inbounds function indices(x::Tuple)
    inds = map(eachindex, x)
    return reduce_tup(_pick_range, inds)
end

"""
    indices(x::Tuple, dim) -> AbstractUnitRange

Returns valid indices for each array in `x` along dimension `dim`
"""
@propagate_inbounds function indices(x::Tuple, dim)
    inds = map(x_i -> indices(x_i, dim), x)
    return reduce_tup(_pick_range, inds)
end

"""
    indices(x::Tuple, dim::Tuple) -> AbstractUnitRange

Returns valid indices given a tuple of arrays `x` and tuple of dimesions for each
respective array (`dim`).
"""
@propagate_inbounds function indices(x::Tuple, dim::Tuple)
    inds = map(indices, x, dim)
    return reduce_tup(_pick_range, inds)
end

"""
    indices(x, dim::Tuple) -> Tuple{Vararg{AbstractUnitRange}}

Returns valid indices for array `x` along each dimension specified in `dim`.
"""
@inline indices(x, dims::Tuple) = _indices(x, dims)
_indices(x, dims::Tuple) = (indices(x, first(dims)), _indices(x, tail(dims))...)
_indices(x, ::Tuple{}) = ()


text:
##################
# GradientResult #
##################

type GradientResult{V,G} <: ForwardDiffResult
    value::V
    gradient::G
end

GradientResult(x) = GradientResult(first(x), similar(x))

value(result::GradientResult) = result.value
gradient(result::GradientResult) = result.gradient

###############
# API methods #
###############

function gradient{N}(f, x, chunk::Chunk{N} = pickchunk(x);
                     multithread::Bool = false,
                     usecache::Bool = true)
    if N == length(x)
        return vector_mode_gradient(f, x, chunk, usecache)
    elseif multithread
        return multithread_chunk_mode_gradient(f, x, chunk, usecache)
    else
        return chunk_mode_gradient(f, x, chunk, usecache)
    end
end

function gradient!{N}(out, f, x, chunk::Chunk{N} = pickchunk(x);
                      multithread::Bool = false,
                      usecache::Bool = true)
    if N == length(x)
        vector_mode_gradient!(out, f, x, chunk, usecache)
    elseif multithread
        multithread_chunk_mode_gradient!(out, f, x, chunk, usecache)
    else
        chunk_mode_gradient!(out, f, x, chunk, usecache)
    end
    return out
end

#######################
# workhorse functions #
#######################

# result extraction #
#-------------------#

@inline load_gradient_value!(out, dual) = out

function load_gradient_value!(out::GradientResult, dual)
    out.value = value(dual)
    return out
end

function load_gradient!(out, dual)
    for i in eachindex(out)
        out[i] = partials(dual, i)
    end
    return out
end

@inline function load_gradient!(out::GradientResult, dual)
    load_gradient!(out.gradient, dual)
    return out
end

function load_gradient_chunk!(out, dual, index, chunksize)
    offset = index - 1
    for i in 1:chunksize
        out[i + offset] = partials(dual, i)
    end
    return out
end

@inline function load_gradient_chunk!(out::GradientResult, dual, index, chunksize)
    load_gradient_chunk!(out.gradient, dual, index, chunksize)
    return out
end

# vector mode #
#-------------#

function compute_vector_mode_gradient(f, x, chunk, usecache)
    cache = jacobian_cachefetch!(x, chunk, usecache)
    xdual = cache.duals
    seed!(xdual, x, cache.seeds, 1)
    return f(xdual)
end

function vector_mode_gradient(f, x, chunk, usecache)
    dual = compute_vector_mode_gradient(f, x, chunk, usecache)
    out = similar(x, numtype(dual))
    return load_gradient!(out, dual)
end

function vector_mode_gradient!(out, f, x, chunk, usecache)
    dual = compute_vector_mode_gradient(f, x, chunk, usecache)
    load_gradient_value!(out, dual)
    load_gradient!(out, dual)
    return out
end

# chunk mode #
#------------#

function chunk_mode_gradient_expr(out_definition::Expr)
    return quote
        @assert length(x) >= N "chunk size cannot be greater than length(x) ($(N) > $(length(x)))"

        # precalculate loop bounds
        xlen = length(x)
        remainder = xlen % N
        lastchunksize = ifelse(remainder == 0, N, remainder)
        lastchunkindex = xlen - lastchunksize + 1
        middlechunks = 2:div(xlen - lastchunksize, N)

        # fetch and seed work vectors
        cache = jacobian_cachefetch!(x, chunk, usecache)
        xdual = cache.duals
        seeds = cache.seeds
        zeroseed = zero(eltype(seeds))
        seedall!(xdual, x, zeroseed)

        # do first chunk manually to calculate output type
        seed!(xdual, x, seeds, 1)
        dual = f(xdual)
        seed!(xdual, x, zeroseed, 1)
        $(out_definition)
        load_gradient_chunk!(out, dual, 1, N)

        # do middle chunks
        for c in middlechunks
            i = ((c - 1) * N + 1)
            seed!(xdual, x, seeds, i)
            dual = f(xdual)
            seed!(xdual, x, zeroseed, i)
            load_gradient_chunk!(out, dual, i, N)
        end

        # do final chunk
        seed!(xdual, x, seeds, lastchunkindex, lastchunksize)
        dual = f(xdual)
        load_gradient_chunk!(out, dual, lastchunkindex, lastchunksize)

        # load value, this is a no-op unless out is a GradientResult
        load_gradient_value!(out, dual)

        return out
    end
end

@eval function chunk_mode_gradient{N}(f, x, chunk::Chunk{N}, usecache)
    $(chunk_mode_gradient_expr(:(out = similar(x, numtype(dual)))))
end

@eval function chunk_mode_gradient!{N}(out, f, x, chunk::Chunk{N}, usecache)
    $(chunk_mode_gradient_expr(:()))
end

# multithreaded chunk mode #
#--------------------------#

if IS_MULTITHREADED_JULIA
    function multithread_chunk_mode_expr(out_definition::Expr)
        return quote
            @assert length(x) >= N "chunk size cannot be greater than length(x) ($(N) > $(length(x)))"

            # precalculate loop bounds
            xlen = length(x)
            remainder = xlen % N
            lastchunksize = ifelse(remainder == 0, N, remainder)
            lastchunkindex = xlen - lastchunksize + 1
            middlechunks = 2:div(xlen - lastchunksize, N)

            # fetch and seed work vectors
            caches = multithread_jacobian_cachefetch!(x, chunk, usecache)
            current_cache = caches[compat_threadid()]
            current_xdual = current_cache.duals
            current_seeds = current_cache.seeds
            zeroseed = zero(eltype(current_seeds))

            Base.Threads.@threads for t in 1:NTHREADS
                seedall!(caches[t].duals, x, zeroseed)
            end

            # do first chunk manually to calculate output type
            seed!(current_xdual, x, current_seeds, 1)
            current_dual = f(current_xdual)
            seed!(current_xdual, x, zeroseed, 1)
            $(out_definition)
            load_gradient_chunk!(out, current_dual, 1, N)

            # do middle chunks
            Base.Threads.@threads for c in middlechunks
                # see https://github.com/JuliaLang/julia/issues/14948
                local chunk_cache = caches[compat_threadid()]
                local chunk_xdual = chunk_cache.duals
                local chunk_seeds = chunk_cache.seeds
                local chunk_index = ((c - 1) * N + 1)
                seed!(chunk_xdual, x, chunk_seeds, chunk_index)
                local chunk_dual = f(chunk_xdual)
                seed!(chunk_xdual, x, zeroseed, chunk_index)
                load_gradient_chunk!(out, chunk_dual, chunk_index, N)
            end

            # do final chunk
            seed!(current_xdual, x, current_seeds, lastchunkindex, lastchunksize)
            current_dual = f(current_xdual)
            load_gradient_chunk!(out, current_dual, lastchunkindex, lastchunksize)

            # load value, this is a no-op unless `out` is a GradientResult
            load_gradient_value!(out, current_dual)

            return out
        end
    end

    @eval function multithread_chunk_mode_gradient{N}(f, x, chunk::Chunk{N}, usecache)
        $(multithread_chunk_mode_expr(:(out = similar(x, numtype(current_dual)))))
    end

    @eval function multithread_chunk_mode_gradient!{N}(out, f, x, chunk::Chunk{N}, usecache)
        $(multithread_chunk_mode_expr(:()))
    end
else
    function multithread_chunk_mode_gradient(args...)
        error("Multithreading is not enabled for this Julia installation.")
    end
    multithread_chunk_mode_gradient!(args...) = multithread_chunk_mode_gradient()
end

ext: jl


text:
BurgersLiquid =  RheoModelClass(
        # Model name
        name="burgerliquid",
        # Model parameters,
        p = [:η₁ , :k₁, :η₂, :k₂],
        # Relaxation modulus
        G = quote
              p1 = η₁/k₁ + η₁/k₂ + η₂/k₂
              p2 = η₁ * η₂ / (k₁ * k₂)
              q1 = η₁
              q2 = η₁ * η₂ / k₂
              A = sqrt(p1^2 - 4*p2)
              r1 = (p1-A)/(2*p2)
              r2 = (p1+A)/(2*p2)
              ((q1-q2*r1)*exp(-r1*t) - (q1-q2*r2)* exp(-r2*t))/A
            end,
        # Creep modulus
        J = quote
              1/k₁ + t/ η₁ + 1/k₂ * (1- exp(-k₂*t/η₂))
            end,
        # Storage modulus
        Gp = quote
                p1 = η₁/k₁ + η₁/k₂ + η₂/k₂
                p2 = η₁ * η₂ / (k₁ * k₂)
                q1 = η₁
                q2 = η₁ * η₂ / k₂
                numerator = p1 * q1 * ω^2 - q2 * ω^2 * (1-p2*ω^2)
                denominator = p1^2 * ω^2 + (1-p2*ω^2)^2
                numerator/denominator
             end,
        # Loss modulus
        Gpp = quote
                p1 = η₁/k₁ + η₁/k₂ + η₂/k₂
                p2 = η₁ * η₂ / (k₁ * k₂)
                q1 = η₁
                q2 = η₁ * η₂ / k₂
                numerator = p1 * q2 * ω^3 + q1 * ω * (1-p2*ω^2)
                denominator = p1^2 * ω^2 + (1-p2*ω^2)^2
                numerator/denominator
              end,
        # Network
        info= "

                                                         ___
                                                 _________| |________
                 ___                            |        _|_| η₂     |
             _____| |________╱╲  ╱╲  ╱╲  _______|                    |____
                 _|_|          ╲╱  ╲╱  ╲╱       |                    |
                    η₁                  k₁      |____╱╲  ╱╲  ╱╲  ____|
                                                       ╲╱  ╲╱  ╲╱
                                                                 k₂

                "
        )


text:
module BenchCubical

using BenchmarkTools
using Ripserer

include(joinpath(@__DIR__, "utils.jl"))

benchmarks = [(file="bonsai64.dipha", dim_max=3), (file="lena2048.dipha", dim_max=2)]

suite = BenchmarkGroup()

for b in benchmarks
    name = b.file
    data = load_data(joinpath(@__DIR__, "data", b.file))
    suite[name] = @benchmarkable ripserer(Cubical($data); dim_max=$(b.dim_max)) samples = 2
end

end
BenchCubical.suite


text:
import NearestNeighbors.MinkowskiMetric
# This contains a bunch of random tests that should hopefully detect if
# some edge case has been missed in the real tests

@testset "metric" for metric in fullmetrics
    @testset "tree type" for TreeType in trees_with_brute
        @testset "type" for T in (Float32, Float64)
            @testset "knn monkey" begin
                # Checks that we find existing point in the tree
                # and that it is the closest
                if TreeType == KDTree && !isa(metric, MinkowskiMetric)
                    continue
                elseif TreeType == BallTree && isa(metric, Hamming)
                    continue
                end
                for i in 1:30
                    dim_data = rand(1:4)
                    size_data = rand(1000:1300)
                    data = rand(T, dim_data, size_data)
                    for j = 1:5
                        tree = TreeType(data, metric; leafsize = rand(1:15))
                        n = rand(1:size_data)
                        idx, dist = knn(tree, data[:,n], rand(1:30), true)
                        @test issorted(dist) == true
                        @test n == idx[1]
                    end
                end

                # Compares vs Brute Force
                for i in 1:30
                    dim_data = rand(1:5)
                    size_data = rand(100:151)
                    data = rand(T, dim_data, size_data)
                    tree = TreeType(data, metric; leafsize = rand(1:15))
                    btree = BruteTree(data, metric)
                    k = rand(1:12)
                    p = rand(dim_data)
                    idx, dist = knn(tree, p, k)
                    bidx, bdist = knn(tree, p, k)
                    @test idx == bidx
                    @test dist ≈ bdist
                end
            end

            @testset "inrange monkey" begin
                # Test against brute force
                for i in 1:30
                    dim_data = rand(1:6)
                    size_data = rand(20:250)
                    data = rand(T, dim_data, size_data)
                    tree = TreeType(data, metric; leafsize = rand(1:8))
                    btree = BruteTree(data, metric)
                    p = 0.5 * ones(dim_data)
                    r = 0.3

                    idxs = inrange(tree, p, r, true)
                    bidxs = inrange(btree, p, r, true)

                    @test idxs == bidxs
                end
            end

            @testset "coupled monkey" begin
                for i in 1:50
                    dim_data = rand(1:5)
                    size_data = rand(100:1000)
                    data = randn(T, dim_data, size_data)
                    tree = TreeType(data, metric; leafsize = rand(1:8))
                    point = randn(dim_data)
                    idxs_ball = Int[]
                    r = 0.1
                    while length(idxs_ball) < 10
                        r *= 2.0
                        idxs_ball = inrange(tree,  point, r, true)
                    end
                    idxs_knn, dists = knn(tree, point, length(idxs_ball))

                    @test sort(idxs_knn) == sort(idxs_ball)
                end
            end
        end
    end
end


text:
tojs(s::Symbol) = string(s)
tojs(v::Vector) = [tojs(v_i) for v_i in v]

function tojs(v::Dict)
    res = Dict()
    for (k, v) in v
        res[k] = tojs(v)
    end
    return res
end
tojs(x::Any) = x

function tojs(x::VegaStream)
    res = Dict()
    if x._type != nothing
        res["type"] = tojs(x._type)
    end
    if x.expr != nothing
        res["expr"] = tojs(x.expr)
    end
    if x.scale != nothing
        res["scale"] = tojs(x.scale)
    end
    return res
end

function tojs(x::VegaDataRef)
    res = Dict()
    if x.data != nothing
        res["data"] = tojs(x.data)
    end
    if x.field != nothing
        res["field"] = tojs(x.field)
    end
    if x.fields != nothing
        res["fields"] = tojs(x.fields)
    end
    if x.sort != nothing
        res["sort"] = tojs(x.sort)
    end
    return res
end

function tojs(x::VegaValueRef)
    res = Dict()
    if x.value != nothing
        res["value"] = tojs(x.value)
    end
    if x.field != nothing
        res["field"] = tojs(x.field)
    end
    if x.group != nothing
        res["group"] = tojs(x.group)
    end
    if x.scale != nothing
        res["scale"] = tojs(x.scale)
    end
    if x.mult != nothing
        res["mult"] = tojs(x.mult)
    end
    if x.offset != nothing
        res["offset"] = tojs(x.offset)
    end
    if x.band != nothing
        res["band"] = tojs(x.band)
    end
    if x.signal != nothing
        res["signal"] = tojs(x.signal)
    end
    return res
end

function tojs(x::VegaPadding)
    res = Dict()
    if x.top != nothing
        res["top"] = tojs(x.top)
    end
    if x.left != nothing
        res["left"] = tojs(x.left)
    end
    if x.bottom != nothing
        res["bottom"] = tojs(x.bottom)
    end
    if x.right != nothing
        res["right"] = tojs(x.right)
    end
    return res
end

function tojs(x::VegaMarkPropertySet)
    res = Dict()
    if x.x != nothing
        res["x"] = tojs(x.x)
    end
    if x.x2 != nothing
        res["x2"] = tojs(x.x2)
    end
    if x.xc != nothing
        res["xc"] = tojs(x.xc)
    end
    if x.width != nothing
        res["width"] = tojs(x.width)
    end
    if x.y != nothing
        res["y"] = tojs(x.y)
    end
    if x.y2 != nothing
        res["y2"] = tojs(x.y2)
    end
    if x.yc != nothing
        res["yc"] = tojs(x.yc)
    end
    if x.height != nothing
        res["height"] = tojs(x.height)
    end
    if x.opacity != nothing
        res["opacity"] = tojs(x.opacity)
    end
    if x.fill != nothing
        res["fill"] = tojs(x.fill)
    end
    if x.fillOpacity != nothing
        res["fillOpacity"] = tojs(x.fillOpacity)
    end
    if x.stroke != nothing
        res["stroke"] = tojs(x.stroke)
    end
    if x.strokeWidth != nothing
        res["strokeWidth"] = tojs(x.strokeWidth)
    end
    if x.strokeOpacity != nothing
        res["strokeOpacity"] = tojs(x.strokeOpacity)
    end
    if x.strokeDash != nothing
        res["strokeDash"] = tojs(x.strokeDash)
    end
    if x.strokeDashOffset != nothing
        res["strokeDashOffset"] = tojs(x.strokeDashOffset)
    end
    if x.cursor != nothing
        res["cursor"] = tojs(x.cursor)
    end
    if x.clip != nothing
        res["clip"] = tojs(x.clip)
    end
    if x.size != nothing
        res["size"] = tojs(x.size)
    end
    if x.shape != nothing
        res["shape"] = tojs(x.shape)
    end
    if x.path != nothing
        res["path"] = tojs(x.path)
    end
    if x.innerRadius != nothing
        res["innerRadius"] = tojs(x.innerRadius)
    end
    if x.outerRadius != nothing
        res["outerRadius"] = tojs(x.outerRadius)
    end
    if x.startAngle != nothing
        res["startAngle"] = tojs(x.startAngle)
    end
    if x.endAngle != nothing
        res["endAngle"] = tojs(x.endAngle)
    end
    if x.interpolate != nothing
        res["interpolate"] = tojs(x.interpolate)
    end
    if x.tension != nothing
        res["tension"] = tojs(x.tension)
    end
    if x.orient != nothing
        res["orient"] = tojs(x.orient)
    end
    if x.url != nothing
        res["url"] = tojs(x.url)
    end
    if x.align != nothing
        res["align"] = tojs(x.align)
    end
    if x.baseline != nothing
        res["baseline"] = tojs(x.baseline)
    end
    if x.text != nothing
        res["text"] = tojs(x.text)
    end
    if x.dx != nothing
        res["dx"] = tojs(x.dx)
    end
    if x.dy != nothing
        res["dy"] = tojs(x.dy)
    end
    if x.radius != nothing
        res["radius"] = tojs(x.radius)
    end
    if x.theta != nothing
        res["theta"] = tojs(x.theta)
    end
    if x.angle != nothing
        res["angle"] = tojs(x.angle)
    end
    if x.font != nothing
        res["font"] = tojs(x.font)
    end
    if x.fontSize != nothing
        res["fontSize"] = tojs(x.fontSize)
    end
    if x.fontWeight != nothing
        res["fontWeight"] = tojs(x.fontWeight)
    end
    if x.fontStyle != nothing
        res["fontStyle"] = tojs(x.fontStyle)
    end
    return res
end

function tojs(x::VegaAxis)
    res = Dict()
    if x._type != nothing
        res["type"] = tojs(x._type)
    end
    if x.scale != nothing
        res["scale"] = tojs(x.scale)
    end
    if x.orient != nothing
        res["orient"] = tojs(x.orient)
    end
    if x.title != nothing
        res["title"] = tojs(x.title)
    end
    if x.titleOffset != nothing
        res["titleOffset"] = tojs(x.titleOffset)
    end
    if x.format != nothing
        res["format"] = tojs(x.format)
    end
    if x.formatType != nothing
        res["formatType"] = tojs(x.formatType)
    end
    if x.ticks != nothing
        res["ticks"] = tojs(x.ticks)
    end
    if x.values != nothing
        res["values"] = tojs(x.values)
    end
    if x.subdivide != nothing
        res["subdivide"] = tojs(x.subdivide)
    end
    if x.tickPadding != nothing
        res["tickPadding"] = tojs(x.tickPadding)
    end
    if x.tickSize != nothing
        res["tickSize"] = tojs(x.tickSize)
    end
    if x.tickSizeMajor != nothing
        res["tickSizeMajor"] = tojs(x.tickSizeMajor)
    end
    if x.tickSizeMinor != nothing
        res["tickSizeMinor"] = tojs(x.tickSizeMinor)
    end
    if x.tickSizeEnd != nothing
        res["tickSizeEnd"] = tojs(x.tickSizeEnd)
    end
    if x.offset != nothing
        res["offset"] = tojs(x.offset)
    end
    if x.layer != nothing
        res["layer"] = tojs(x.layer)
    end
    if x.grid != nothing
        res["grid"] = tojs(x.grid)
    end
    if x.properties != nothing
        res["properties"] = tojs(x.properties)
    end
    return res
end

function tojs(x::VegaFormat)
    res = Dict()
    if x._type != nothing
        res["type"] = tojs(x._type)
    end
    if x.parse != nothing
        res["parse"] = tojs(x.parse)
    end
    if x.property != nothing
        res["property"] = tojs(x.property)
    end
    if x.feature != nothing
        res["feature"] = tojs(x.feature)
    end
    if x.mesh != nothing
        res["mesh"] = tojs(x.mesh)
    end
    if x.children != nothing
        res["children"] = tojs(x.children)
    end
    return res
end

function tojs(x::VegaData)
    res = Dict()
    if x.name != nothing
        res["name"] = tojs(x.name)
    end
    if x.format != nothing
        res["format"] = tojs(x.format)
    end
    if x.values != nothing
        res["values"] = tojs(x.values)
    end
    if x.source != nothing
        res["source"] = tojs(x.source)
    end
    if x.url != nothing
        res["url"] = tojs(x.url)
    end
    if x.transform != nothing
        res["transform"] = tojs(x.transform)
    end
    if x.modify != nothing
        res["modify"] = tojs(x.modify)
    end
    return res
end

function tojs(x::VegaMarkProperties)
    res = Dict()
    if x.enter != nothing
        res["enter"] = tojs(x.enter)
    end
    if x.exit != nothing
        res["exit"] = tojs(x.exit)
    end
    if x.update != nothing
        res["update"] = tojs(x.update)
    end
    if x.hover != nothing
        res["hover"] = tojs(x.hover)
    end
    return res
end

function tojs(x::VegaMarkFrom)
    res = Dict()
    if x.data != nothing
        res["data"] = tojs(x.data)
    end
    if x.mark != nothing
        res["mark"] = tojs(x.mark)
    end
    if x.transform != nothing
        res["transform"] = tojs(x.transform)
    end
    if x.value != nothing
        res["value"] = tojs(x.value)
    end
    return res
end

function tojs(x::VegaSignal)
    res = Dict()
    if x.name != nothing
        res["name"] = tojs(x.name)
    end
    if x.init != nothing
        res["init"] = tojs(x.init)
    end
    if x.verbose != nothing
        res["verbose"] = tojs(x.verbose)
    end
    if x.expr != nothing
        res["expr"] = tojs(x.expr)
    end
    if x.scale != nothing
        res["scale"] = tojs(x.scale)
    end
    if x.streams != nothing
        res["streams"] = tojs(x.streams)
    end
    return res
end

function tojs(x::VegaScale)
    res = Dict()
    if x.name != nothing
        res["name"] = tojs(x.name)
    end
    if x._type != nothing
        res["type"] = tojs(x._type)
    end
    if x.domain != nothing
        res["domain"] = tojs(x.domain)
    end
    if x.domainMin != nothing
        res["domainMin"] = tojs(x.domainMin)
    end
    if x.domainMax != nothing
        res["domainMax"] = tojs(x.domainMax)
    end
    if x.range != nothing
        res["range"] = tojs(x.range)
    end
    if x.rangeMin != nothing
        res["rangeMin"] = tojs(x.rangeMin)
    end
    if x.rangeMax != nothing
        res["rangeMax"] = tojs(x.rangeMax)
    end
    if x.reverse != nothing
        res["reverse"] = tojs(x.reverse)
    end
    if x.round != nothing
        res["round"] = tojs(x.round)
    end
    if x.points != nothing
        res["points"] = tojs(x.points)
    end
    if x.padding != nothing
        res["padding"] = tojs(x.padding)
    end
    if x.outerPadding != nothing
        res["outerPadding"] = tojs(x.outerPadding)
    end
    if x.bandWidth != nothing
        res["bandWidth"] = tojs(x.bandWidth)
    end
    if x.sort != nothing
        res["sort"] = tojs(x.sort)
    end
    if x.clamp != nothing
        res["clamp"] = tojs(x.clamp)
    end
    if x.nice != nothing
        res["nice"] = tojs(x.nice)
    end
    if x.exponent != nothing
        res["exponent"] = tojs(x.exponent)
    end
    if x.zero != nothing
        res["zero"] = tojs(x.zero)
    end
    return res
end

function tojs(x::VegaMark)
    res = Dict()
    if x.name != nothing
        res["name"] = tojs(x.name)
    end
    if x.key != nothing
        res["key"] = tojs(x.key)
    end
    if x._type != nothing
        res["type"] = tojs(x._type)
    end
    if x.from != nothing
        res["from"] = tojs(x.from)
    end
    if x.properties != nothing
        res["properties"] = tojs(x.properties)
    end
    if x.delay != nothing
        res["delay"] = tojs(x.delay)
    end
    if x.ease != nothing
        res["ease"] = tojs(x.ease)
    end
    if x.marks != nothing
        res["marks"] = tojs(x.marks)
    end
    if x.scales != nothing
        res["scales"] = tojs(x.scales)
    end
    if x.axes != nothing
        res["axes"] = tojs(x.axes)
    end
    if x.interactive != nothing
        res["interactive"] = tojs(x.interactive)
    end
    return res
end

function tojs(x::VegaLegend)
    res = Dict()
    if x.size != nothing
        res["size"] = tojs(x.size)
    end
    if x.shape != nothing
        res["shape"] = tojs(x.shape)
    end
    if x.fill != nothing
        res["fill"] = tojs(x.fill)
    end
    if x.stroke != nothing
        res["stroke"] = tojs(x.stroke)
    end
    if x.orient != nothing
        res["orient"] = tojs(x.orient)
    end
    if x.offset != nothing
        res["offset"] = tojs(x.offset)
    end
    if x.title != nothing
        res["title"] = tojs(x.title)
    end
    if x.format != nothing
        res["format"] = tojs(x.format)
    end
    if x.formatType != nothing
        res["formatType"] = tojs(x.formatType)
    end
    if x.values != nothing
        res["values"] = tojs(x.values)
    end
    if x.properties != nothing
        res["properties"] = tojs(x.properties)
    end
    return res
end

function tojs(x::VegaScene)
    res = Dict()
    if x.fill != nothing
        res["fill"] = tojs(x.fill)
    end
    if x.fillOpacity != nothing
        res["fillOpacity"] = tojs(x.fillOpacity)
    end
    if x.stroke != nothing
        res["stroke"] = tojs(x.stroke)
    end
    if x.strokeOpacity != nothing
        res["strokeOpacity"] = tojs(x.strokeOpacity)
    end
    if x.strokeWidth != nothing
        res["strokeWidth"] = tojs(x.strokeWidth)
    end
    if x.strokeDash != nothing
        res["strokeDash"] = tojs(x.strokeDash)
    end
    if x.strokeDashOffset != nothing
        res["strokeDashOffset"] = tojs(x.strokeDashOffset)
    end
    return res
end

function tojs(x::VegaContainer)
    res = Dict()
    if x.scene != nothing
        res["scene"] = tojs(x.scene)
    end
    if x.scales != nothing
        res["scales"] = tojs(x.scales)
    end
    if x.axes != nothing
        res["axes"] = tojs(x.axes)
    end
    if x.legends != nothing
        res["legends"] = tojs(x.legends)
    end
    if x.marks != nothing
        res["marks"] = tojs(x.marks)
    end
    return res
end

function tojs(x::VegaPredicate)
    res = Dict()
    if x.name != nothing
        res["name"] = tojs(x.name)
    end
    if x._type != nothing
        res["type"] = tojs(x._type)
    end
    if x.operands != nothing
        res["operands"] = tojs(x.operands)
    end
    if x.range != nothing
        res["range"] = tojs(x.range)
    end
    if x.scale != nothing
        res["scale"] = tojs(x.scale)
    end
    if x.item != nothing
        res["item"] = tojs(x.item)
    end
    return res
end

function tojs(x::VegaVisualization)
    res = Dict()
    if x.name != nothing
        res["name"] = tojs(x.name)
    end
    if x.width != nothing
        res["width"] = tojs(x.width)
    end
    if x.height != nothing
        res["height"] = tojs(x.height)
    end
    if x.padding != nothing
        res["padding"] = tojs(x.padding)
    end
    if x.viewport != nothing
        res["viewport"] = tojs(x.viewport)
    end
    if x.background != nothing
        res["background"] = tojs(x.background)
    end
    if x.data != nothing
        res["data"] = tojs(x.data)
    end
    if x.scales != nothing
        res["scales"] = tojs(x.scales)
    end
    if x.axes != nothing
        res["axes"] = tojs(x.axes)
    end
    if x.marks != nothing
        res["marks"] = tojs(x.marks)
    end
    if x.legends != nothing
        res["legends"] = tojs(x.legends)
    end
    if x.signals != nothing
        res["signals"] = tojs(x.signals)
    end
    if x.predicates != nothing
        res["predicates"] = tojs(x.predicates)
    end
    if x.modify != nothing
        res["modify"] = tojs(x.modify)
    end
    return res
end

max_stars_repo_name: jpthor/vega


text:
"""
    BOX{N, AM, S, D, R} <: AbstractContinuousPost

Implementation of a reachability method for linear systems using box approximations.

## Fields

- `δ`            -- step-size of the discretization
- `approx_model` -- (optional, default: `Forward`) approximation model;
                    see `Notes` below for possible options
- `static`
- `dim`          -- (optional default: `missing`) ambient dimension
- `recursive`    -- (optional default: `false`) if `true`, use the implementation that
                    recursively computes each reach-set; otherwise, use the implementation
                    that unwraps the sequence until the initial set

## Notes

The type fields are:

- `N`  -- number type of the step-size
- `AM` -- approximation model
- `S`  -- value type for the static option
- `D`  -- value type for the dimension
- `R`  -- value type for the recursive option

The default approximation model is:

```julia
Forward(sih=:concrete, exp=:base, setops=:lazy)
```

This algorithm solves the set-based recurrence equation ``X_{k+1} = ΦX_k ⊕ V_k``
by computing a tight hyperrectangular over-approximation of ``X_{k+1}`` at each
step ``k ∈ \\mathbb{N}``. The recursive implementation uses the previously computed
set ``X_k`` to compute ``X_{k+1}``. However, it is known that this method incurs
wrapping effects. The non-recursive implementation instead computes ``X_{k+1}``
by unwrapping the discrete recurrence until ``X_0 = Ω₀``, at the expense of computing
powers of the matrix ``Φ``. These ideas are discussed in [[BFFPSV18]](@ref).

### References

This algorithm is essentially a non-decomposed version of the method in [[BFFPSV18]](@ref),
using hyperrectangles as set representation. For a general introduction we refer
to the dissertation [[LG09]](@ref).

Regarding the approximation model, by default we use an adaptation of the method
presented in [[FRE11]](@ref).
"""
struct BOX{N, AM, S, D, R} <: AbstractContinuousPost
    δ::N
    approx_model::AM
    static::S
    dim::D
    recursive::R
end

# convenience constructor using symbols
function BOX(; δ::N,
               approx_model::AM=Forward(sih=:concrete, exp=:base, setops=:lazy),
               static::Bool=false,
               dim::Union{Int, Missing}=missing,
               recursive::Bool=false) where {N, AM}
    n = !ismissing(dim) ? Val(dim) : dim
    return BOX(δ, approx_model, Val(static), n, Val(recursive))
end

step_size(alg::BOX) = alg.δ
numtype(::BOX{N}) where {N} = N

function setrep(::BOX{N, AM, Val{false}, D, R}) where {N, AM, D, R}
    VT = Vector{N}
    ST = Hyperrectangle{N, VT, VT}
end

function rsetrep(::BOX{N, AM, Val{true}, Val{n}, R}) where {N, AM, n, R}
    VT = SVector{n, N}
    ST = Hyperrectangle{N, VT, VT}
end

function rsetrep(alg::BOX{N}) where {N}
    ST = setrep(alg)
    RT = ReachSet{N, ST}
end

include("post.jl")
include("reach_homog.jl")
include("reach_inhomog.jl")


text:
# Autogenerated wrapper script for GStreamer_jll for armv7l-linux-gnueabihf
export gst_inspect, gst_launch, gst_stats, gst_tester, gst_typefind, libgstbase, libgstcheck, libgstcontroller, libgstnet, libgstreamer

using Glib_jll
using LibUnwind_jll
using Elfutils_jll
using GMP_jll
using GSL_jll
using libcap_jll
JLLWrappers.@generate_wrapper_header("GStreamer")
JLLWrappers.@declare_executable_product(gst_inspect)
JLLWrappers.@declare_executable_product(gst_launch)
JLLWrappers.@declare_executable_product(gst_stats)
JLLWrappers.@declare_executable_product(gst_tester)
JLLWrappers.@declare_executable_product(gst_typefind)
JLLWrappers.@declare_library_product(libgstbase, "libgstbase-1.0.so.0")
JLLWrappers.@declare_library_product(libgstcheck, "libgstcheck-1.0.so.0")
JLLWrappers.@declare_library_product(libgstcontroller, "libgstcontroller-1.0.so.0")
JLLWrappers.@declare_library_product(libgstnet, "libgstnet-1.0.so.0")
JLLWrappers.@declare_library_product(libgstreamer, "libgstreamer-1.0.so.0")
function __init__()
    JLLWrappers.@generate_init_header(Glib_jll, LibUnwind_jll, Elfutils_jll, GMP_jll, GSL_jll, libcap_jll)
    JLLWrappers.@init_executable_product(
        gst_inspect,
        "bin/gst-inspect-1.0",
    )

    JLLWrappers.@init_executable_product(
        gst_launch,
        "bin/gst-launch-1.0",
    )

    JLLWrappers.@init_executable_product(
        gst_stats,
        "bin/gst-stats-1.0",
    )

    JLLWrappers.@init_executable_product(
        gst_tester,
        "bin/gst-tester-1.0",
    )

    JLLWrappers.@init_executable_product(
        gst_typefind,
        "bin/gst-typefind-1.0",
    )

    JLLWrappers.@init_library_product(
        libgstbase,
        "lib/libgstbase-1.0.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgstcheck,
        "lib/libgstcheck-1.0.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgstcontroller,
        "lib/libgstcontroller-1.0.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgstnet,
        "lib/libgstnet-1.0.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgstreamer,
        "lib/libgstreamer-1.0.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using ZXCalculus, Test
using YaoHIR, YaoLocations
using CompilerPluginTools

chain = Chain()
push_gate!(chain, Val{:Sdag}(), 1)
push_gate!(chain, Val{:H}(), 1)
push_gate!(chain, Val{:S}(), 1)
push_gate!(chain, Val{:S}(), 2)
push_gate!(chain, Val{:H}(), 4)
push_gate!(chain, Val{:CNOT}(), 3, 2)
push_gate!(chain, Val{:CZ}(), 4, 1)
push_gate!(chain, Val{:H}(), 2)
push_gate!(chain, Val{:T}(), 2)
push_gate!(chain, Val{:CNOT}(), 3, 2)
push_gate!(chain, Val{:Tdag}(), 2)
push_gate!(chain, Val{:CNOT}(), 1, 4)
push_gate!(chain, Val{:H}(), 1)
push_gate!(chain, Val{:T}(), 2)
push_gate!(chain, Val{:S}(), 3)
push_gate!(chain, Val{:H}(), 4)
push_gate!(chain, Val{:T}(), 1)
push_gate!(chain, Val{:H}(), 2)
push_gate!(chain, Val{:H}(), 3)
push_gate!(chain, Val{:Sdag}(), 4)
push_gate!(chain, Val{:S}(), 3)
push_gate!(chain, Val{:X}(), 4)
push_gate!(chain, Val{:CNOT}(), 3, 2)
push_gate!(chain, Val{:H}(), 1)
push_gate!(chain, Val{:S}(), 4)
push_gate!(chain, Val{:X}(), 4)

ir = @make_ircode begin
end
bir = BlockIR(ir, 4, chain)
zxd = convert_to_zxd(bir)
pt_zxd = phase_teleportation(zxd)
@test tcount(pt_zxd) <= tcount(zxd)
pt_chain = convert_to_chain(pt_zxd)
@test length(pt_chain) <= length(chain)

zxg = clifford_simplification(zxd)
cl_chain = circuit_extraction(zxg)

zxg = full_reduction(zxd)
fl_chain = circuit_extraction(zxg)

pt_bir = phase_teleportation(bir)
cl_bir = clifford_simplification(bir)
fl_bir = full_reduction(bir)

@test length(pt_chain) == length(pt_bir.circuit)
@test length(cl_chain) == length(cl_bir.circuit)
@test length(fl_chain) == length(fl_bir.circuit)


text:
"""
    CelestialBody

Struct representing a gravitating body with a NAIF ID code.
"""
struct CelestialBody
    gm::Float64
    naif_id::Int64
end

"""
    CelestialBody(gm)

Construct a new [`CelestialBody`](@ref) from its gravitational parameter.
"""
CelestialBody(gm) = CelestialBody(gm, -1)

"""
    CelestialBody(gm, naif_id, radius::Float64)

Get the gravitational parameter of the body, often denoted μ.
"""
gravitational_parameter(cb::CelestialBody) = cb.gm

"""
    mass(::CelestialBody)

Return the mass of the body calculated as μ/G.

!!! note

    The mass calculation is performed by dividing the gravitational parameter
    by the universal gravitational constant, G. The uncertainty in the
    gravitational constant is much greater than that in the body's gravitational
    parameter. Therefore, try to use the parameter when possible.
"""
mass(cb::CelestialBody) = cb.gm / GRAVITATIONAL_CONSTANT

"""
    naif_id(::CelestialBody)

Return the NAIF ID of the body
"""
naif_id(cb::CelestialBody) = cb.naif_id


text:
# This file contains the functions to compute the emissivities of bulks and multilayers

# Factors
"""
    bose_einstein(w,T)

Bose-Einstein factor using the definition:
```math
\\Theta(\\omega,T) = \\frac{\\hbar\\omega}{e^{\\frac{\\hbar\\omega}{k_B T}}-1}
```

"""
function bose_einstein(w,T)
    u = w*ħ/kb/T
    if T==0
        return 0.0
    end
    if w==0.0
        return kb*T
    end

    return u*kb*T/(exp(u)-1.0)
end

"""
    planck(w,T)

Planck's Distribution using the definition:

```math
q^{\\text{BB}}(\\omega,T) =  \\Theta(\\omega,T) \\frac{k_0^2}{4\\pi^2}
```
"""
function planck(w,T)
    return  bose_einstein(w,T)*(w/c0)^2/(2.0*pi)^2
end

"""
    planck_fraction(w1,w2,T)

fraction of the energy between pulsations w1 and w2

```math
q^{\\text{BB}}(\\omega_1,\\omega_2,T) =  \\frac{1}{\\sigma T^4}\\int_{\\omega_1}^{\\omega_2}q^{\\text{BB}}(\\omega)d\\omega
```
"""
function planck_fraction(w1,w2,T)
  pl(w) = planck(w,T)
  (val, err) = quadgk(pl, w1, w2 ;rtol=1e-8)
  return val/sigma/T^4
end

"""
    wien(T)

Wien frequency: return the frequency for which Planck distribution is maximum.
T is given in Kelvin and the result is given in radHz
"""
function wien(T)
    return   2.8214393721220787*kb*T/ħ
end

"""
    invwien(w)

Inverse of Wien law : return the temperature (K)
for a given frequency in radHz
"""
function invwien(w)
    return   w*ħ/(2.8214393721220787*kb)
end

" Unit conversion from rad/s to wavelength in m"
unitconv(w) = 2.0*pi*c0/w

" Wien's wavelength obtained from wien's law in m"
lambda_wien(T) = 2.8977729e-3/T

"Monocromatic directional emissivity"
function emissivity_kx_w(structure, kx, w)

    (Rte,Tte)=power_rt(structure, te(), kx ,w)
    (Rtm,Ttm)=power_rt(structure, tm(), kx ,w)

    if imag(permittivity(substrate(structure),w)) == 0.0
        integr1 = 1.0 - Rte - Tte
        integr2 = 1.0 - Rtm - Ttm
    else
        integr1 = 1.0 - Rte
        integr2 = 1.0 - Rtm
    end

    return integr1 + integr2
end

"Total directional emissivity"
function emissivity_kx(structure,T, kx, wi,wf; rtol=1e-8)

    e(u) = emissivity_kx_w(structure, kx, u*kb*T/ħ)*u^3/(exp(u)-1.0)
    val :: Float64  = 0.0
    err :: Float64  = 0.0
    (val,err) = quadgk(e, wi*ħ/kb/T , wf*ħ/kb/T ; rtol=rtol)

    return val*kb^4/ħ^3/c0^2/(2.0*pi)^2/sigma
end

" Monocromatic hemispherical emissivity"
function emissivity_w(structure, w; rtol = 1e-8)
    e_kx(kx) = kx*emissivity_kx_w(structure ,kx, w)

    val :: Float64  = 0.0
    err :: Float64  = 0.0
    (val,err) = quadgk(e_kx, 0.0, w/c0 ; rtol=rtol)

    return val/(w/c0)^2
end

" total hemispherical emissivity"
function emissivity(structure,T;rtol = 1e-8)
    e(u) = emissivity_w(structure, u*kb*T/ħ;rtol = rtol)*u^3/(exp(u)-1.0)
    e2(t) = e(t/(1.0-t))/(1.0-t)^2

    val :: Float64  = 0.0
    err :: Float64  = 0.0
    (val,err) = quadgk(e2, 0.0, 1.0 ; rtol=rtol)

    return val*kb^4/ħ^3/c0^2/(2.0*pi)^2/sigma
end

" total hemispherical emissivity with integration bounds for frequency"
function emissivity(structure,T,wi,wf;rtol = 1e-8)
    e(u) = emissivity_w(structure, u*kb*T/ħ; rtol = rtol)*u^3/(exp(u)-1.0)
    val :: Float64  = 0.0
    err :: Float64  = 0.0
    (val,err) = quadgk(e, wi*ħ/kb/T , wf*ħ/kb/T ; rtol=rtol)

    return val*kb^4/ħ^3/c0^2/(2.0*pi)^2/sigma
end

function emissivity_fraction(structure,T,wi,wf;rtol = 1e-8)
    return emissivity_kx(structure,T,0.0,wi,wf;rtol = rtol)/planck_fraction(wi,wf,T)
end


text:
# ndgrid-repeat.jl
# test lazy and non-lazy version

#using LazyGrids: ndgrid, ndgrid_array
using LazyGrids: btime
include("../src/ndgrid-repeat.jl") # ndgrid_array_r, ndgrid_lazy, ndgrid_repeat

using Test: @test, @testset, @test_throws, @inferred

@testset "ndgrid_array_r" begin

    # easy test with vectors of the same type
    xg, yg, zg = @inferred ndgrid_array_r(1:4, 5:7, 8:9)
    @test xg == repeat(1:4, 1, 3, 2)
    @test yg == repeat((5:7)', 4, 1, 2)

    # harder test with vectors of different types
    x = 1:4
    y = [1//2, 3//4]
    z = [:a, :b, :c]

    @test_throws BoundsError ndgrid_repeat(x, (4,2,3), 0)
    @test_throws DimensionMismatch ndgrid_repeat(x, (4,2,3), 2)

    @inferred ndgrid_repeat(x, (4,2,3), 1)
    @inferred ndgrid_repeat(y, (4,2,3), 2)
    @inferred ndgrid_repeat(z, (4,2,3), 3)
    xg, yg, zg = @inferred ndgrid_array_r(x, y, z)
    @test zg == repeat(reshape(z, (1,1,3)), 4, 2, 1)
end

@testset "ndgrid_lazy" begin
    x = 1:4
    y = [1//2, 3//4]
    z = [:a, :b, :c]
    a = ndgrid_array_r(x, y)
    b = @inferred ndgrid_lazy(x, y)
    @test a == b
    a = ndgrid_array_r(x, y, z)
#   b = @inferred ndgrid_lazy(x, y, z) # @inferred fails here
    b = ndgrid_lazy(x, y, z)
    @test a == b

    # indexing
    (xn, _, _) = ndgrid_array_r(x, y, z)
    (xl, _, _) = ndgrid_lazy(x, y, z)

    @test xn == xl
    @test xn[3] == xl[3]
    @test xn[3] == xl[3]
    @test xn[:] == xl[:]
    @test size(xn) == size(xl)
    @test ndims(xn) == ndims(xl)
    @test xl isa AbstractArray

    @test sizeof(xl) < 100 # small!
end

# copyto! test using LazyArrays version for comparison
using BenchmarkTools: @benchmark

dims = (2^7,2^8,2^9)
x,y,z = map(rand, dims)
(xl, _, _) = ndgrid_lazy(map(Base.OneTo, dims)...)
xa = Array(xl)
out = zeros(eltype(x), dims)
sizeof.((xl, xa))

ta = @benchmark copyto!(out, xa) # 15.9ms
@show btime(ta)

tl = @benchmark copyto!(out, xl) # 179.9ms with lots of memory
@show btime(tl)

#=
@which getindex(xn, 3) # calls Base
@which getindex(xl, 1, 1, 1) # calls specialized getindex
=#


text:
using ViewerGL
using Distributed
if nprocs() == 1
    addprocs(3)
	@info "adding 3 more CPUs"
end
# using Revise
using Test
using Logging
using SparseArrays
@everywhere using LarSurf

    block_size = [2, 2, 2]
    segmentation = LarSurf.data234()

    LarSurf.lsp_setup(block_size)

    # @debug "Setup done"
    larmodel = LarSurf.lsp_get_surface(segmentation)
	V, FV = larmodel
    # @test LarSurf.check_surface_euler(larmodel[2])
    FVtri = LarSurf.triangulate_quads(FV)
    # @test LarSurf.check_surface_euler(FVtri)

# ViewerGL.VIEW([
#     ViewerGL.GLGrid(V,FV,ViewerGL.Point4d(1,1,1,0.1))
# 	ViewerGL.GLAxis(ViewerGL.Point3d(-1,-1,-1),ViewerGL.Point3d(1,1,1))
# ])
ViewerGL.VIEW([
    ViewerGL.GLGrid(V,FVtri,ViewerGL.Point4d(1,1,1,0.1))
	ViewerGL.GLAxis(ViewerGL.Point3d(-1,-1,-1),ViewerGL.Point3d(1,1,1))
])


text:
using MUMPSjInv
using Test

include("getDivGrad.jl")
A = getDivGrad(32, 32, 16)
n = size(A, 1)

# REAL: test for single rhs
println("Test for real SPD matrix: one sparse random rhs")
rhs = sprandn(n, 1, 0.05)
x = solveMUMPS(A, rhs, 1)
err = norm(A * x - rhs) / norm(rhs)
@test err < 1e-14

# rhs as sparse vector instead of n X 1 sparse matrix
rhs = SparseVector(rhs)
x = solveMUMPS(A, rhs, 1)
err = norm(A * x - rhs) / norm(rhs)
@test err < 1e-14

# REAL: test for multiple rhs
println("Test for real SPD matrix: multiple sparse random rhs")
nrhs = 10
rhs = sprandn(n, nrhs, 0.03)
x = solveMUMPS(A, rhs, 1)
err = zeros(nrhs)
for i = 1:nrhs
    err[i] = norm(A * x[:, i] - rhs[:, i]) / norm(rhs[:, i])
end
@test eltype(x) == Float64
@test maximum(err) < 1e-14

println("Test for real SPD matrix: multiple sparse delta function rhs")
rhs = sparse(1.0I, n, n)
rhs = rhs[:, 1:nrhs]
x = solveMUMPS(A, rhs, 1)
err = zeros(nrhs)
for i = 1:nrhs
    err[i] = norm(A * x[:, i] - rhs[:, i]) / norm(rhs[:, i])
end
@test eltype(x) == Float64
@test maximum(err) < 1e-14

# COMPLEX: test for one rhs
println("Test for complex SPD matrix: one rhs")
r = rand(n)
A = A + im*spdiagm(0 => r)
rhs = sprandn(n, 1, 0.03) + im*sprandn(n, 1, 0.03)
x = solveMUMPS(A, rhs, 2)
err = norm(A * x - rhs) / norm(rhs)
@test eltype(x) == ComplexF64
@test err < 1e-14

# sparse vector rhs
rhs = SparseVector(sprandn(n, 1, 0.03) + im*sprandn(n, 1, 0.03))
x = solveMUMPS(A, rhs, 2)
err = norm(A * x - rhs) / norm(rhs)
@test eltype(x) == ComplexF64
@test err < 1e-14

# COMPLEX : test for multiple rhs
println("Test for complex SPD matrix: multiple sparse random rhs")
nrhs = 10
rhs = sprandn(n, nrhs, 0.03) + im * sprandn(n, nrhs, 0.03)
x = solveMUMPS(A, rhs, 2)
err = zeros(nrhs)
for i = 1:nrhs
     err[i] = norm(A * x[:, i] - rhs[:, i]) / norm(rhs[:, i])
end
@test eltype(x) == ComplexF64
@test maximum(err) < 1e-14


text:
using FilePathsBase
using LinearAlgebra
using Test

include("testpaths.jl")

@testset "FilePathsBase" begin
    include("mode.jl")
    include("path.jl")

    # Test that our weird registered path works
    path = p"test|;;my;weird;test;path"
    @test path.parts == ("test|", "", "my", "weird", "test", "path")
end


text:
function weights!(w::Vector{Float64},c::Vector{Float64},ϵ::Float64,IS::Vector{Float64},
    s::Vector{Float64})
    for i = 1:length(w)
        w[i] = c[i]/((ϵ + IS[i])^2);
    end
    s[1] = sum(w);
    for i = 1:length(w)
        w[i] /= s[1];
    end
end

max_stars_repo_name: dcanuto/leg-hemodynamics

max_issues_repo_name: dcanuto/leg-hemodynamics

max_forks_repo_name: dcanuto/leg-hemodynamics


text:
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.6.2
#     language: julia
#     name: julia-1.6
# ---

# %% [markdown]
# https://projecteuler.net/index.php?section=problems&id=30
#
# https://discourse.julialang.org/t/reduce-number-of-allocations/66990

# %%
using BenchmarkTools

# %%
function power_digit_sum(pow, n)
    return sum(c^pow for c in reverse(digits(n)))
end

function Problem30()
    ans = sum(i for i in 2:1_000_000 if i == power_digit_sum(5, i))
    return ans
end

@show Problem30();
@btime Problem30();

# %% [markdown]
# https://docs.julialang.org/en/v1/base/numbers/#Base.digits!

# %%
module Rev1

function power_digit_sum!(pow, n, ws)
    return sum(c^pow for c in digits!(ws, n))
end

function Problem30(N, m)
    ws = Vector{Int}(undef, N)
    ans = sum(i for i in 2:(10^N - 1) if i == power_digit_sum!(m, i, ws))
    return ans
end

end

@show Rev1.Problem30(5, 4);
@show Rev1.Problem30(6, 5);
@btime Rev1.Problem30(6, 5);

# %%
module DNF1

function powsum(pow, n)
    s = 0
    while n > 0
        (n, r) = divrem(n, 10)
        s += r^pow
    end
    return s
end

problem30() = sum(i for i in 2:1_000_000 if i == powsum(5, i))

end

@show DNF1.problem30();
@btime DNF1.problem30();

# %%
module DNF2

function pow5(x)
    y = x^2
    return y^2 * x
end

function pow5sum(n)
    s = 0
    while n >= 10
        (n, r) = divrem(n, 10)
        s += pow5(r)
    end
    return s + pow5(n)
end

problem30() = sum(i for i in 2:1_000_000 if i == pow5sum(i))

end

@show DNF2.problem30();
@btime DNF2.problem30();

# %% [markdown]
# https://docs.julialang.org/en/v1/devdocs/cartesian/

# %%
using Base.Cartesian

# %%
(@macroexpand @nexprs 4 d -> s += 10^(d-1)*i_d) |> Base.remove_linenums!

# %%
(@macroexpand @nloops 4 i d -> 0:9 begin
    @nexprs 4 d -> s += 10^(d-1)*i_d
end) |> Base.remove_linenums!

# %%
N = 4
s = 0
@eval @nloops $N i d -> 0:9 begin
    @nexprs $N d -> global s += 10^(d-1)*i_d
end
s

# %%
sum(0:9999)

# %%
begin
    begin
        a = Int[]
        s = 0
        x = 0
        @nloops 5 i (d -> 0:9) begin
            y = 0
            @nexprs 5 d -> y += i_d ^ 4
            if x > 1 && x == y
                push!(a, x)
                s += x
            end
            x += 1
        end
        a, s
    end
end

# %%
(@macroexpand begin
    begin
        p = @ntuple 10 i -> (i - 1) ^ 5
        a = Int[]
        s = 0
        x = 0
        @inbounds @nloops 6 i d -> 0:9 d -> y_d = p[i_d+1] begin
            y = 0
            @nexprs 6 d -> y += y_d
            if x > 1 && x == y
                push!(a, x)
                s += x
            end
            x += 1
        end
        a, s
    end
end) |> Base.remove_linenums!

# %%
@generated function _projecteuler30v(::Val{N}, ::Val{m}) where {N, m}
    quote
        p = @ntuple 10 i -> (i - 1) ^ $m
        a = Int[]
        s = 0
        x = 0
        @inbounds @nloops $N i d -> 0:9 d -> y_d = p[i_d+1] begin
            y = 0
            @nexprs $N d -> y += y_d
            if x > 1 && x == y
                push!(a, x)
                s += x
            end
            x += 1
        end
        a, s
    end
end
projecteuler30v(N, m) = _projecteuler30v(Val(N), Val(m))

@show projecteuler30v(5, 4);
@show projecteuler30v(6, 5);
@btime projecteuler30v(6, 5);

# %%
@generated function _projecteuler30(::Val{N}, ::Val{m}) where {N, m}
    quote
        p = @ntuple 10 i -> (i - 1) ^ $m
        s = 0
        x = 0
        @inbounds @nloops $N i d -> 0:9 begin
            y = 0
            @nexprs $N d -> y += p[i_d+1]
            s += (x == y) * x
            x += 1
        end
        s - 1
    end
end
projecteuler30(N, m) = _projecteuler30(Val(N), Val(m))

@show projecteuler30(5, 4);
@show projecteuler30(6, 5);
@btime projecteuler30(6, 5);

# %%
@generated function _projecteuler30x(::Val{N}, ::Val{m}) where {N, m}
    quote
        p = @ntuple 10 i -> (i - 1) ^ $m
        s = 0
        x = 0
        @inbounds @nloops $N i d -> 0:9 d -> begin
            y_d = p[i_d+1]
        end begin
            y = 0
            @nexprs $N d -> y += y_d
            s += (x == y) * x
            x += 1
        end
        s - 1
    end
end
projecteuler30x(N, m) = _projecteuler30x(Val(N), Val(m))

@show projecteuler30x(5, 4);
@show projecteuler30x(6, 5);
@btime projecteuler30x(6, 5);

# %%
function projecteuler30p()
    p = @ntuple 10 i -> (i - 1) ^ 5
    s = 0
    x = 0
    @inbounds @nloops 6 i d -> 0:9 begin
        y = 0
        @nexprs 6 d -> y += p[i_d+1]
        s += (x == y) * x
        x += 1
    end
    s - 1
end

@show projecteuler30p();
@btime projecteuler30p();

# %%
function projecteuler30()
    s = 0
    x = 0
    @nloops 6 i d -> 0:9 begin
        y = 0
        @nexprs 6 d -> y += (z = i_d^2; z^2*i_d)
        s += (x == y) * x
        x += 1
    end
    s - 1
end

@show projecteuler30();
@btime projecteuler30();

# %%
function projecteuler30x()
    s = 0
    x = 0
    @nloops 6 i d -> 0:9 d -> begin
        z = i_d^2
        y_d = z^2 * i_d
    end begin
        y = 0
        @nexprs 6 d -> y += y_d
        s += (x == y) * x
        x += 1
    end
    s - 1
end

@show projecteuler30x();
@btime projecteuler30x();

# %%
[(k, k * 9^4, log10(k * 9^4)+1) for k in 1:9]

# %%
[(k, k * 9^5, log10(k * 9^5)+1) for k in 1:9]

# %%
[(k, k * 9^10, log10(k * 9^10)+1) for k in 1:12]

# %%
[(m, findfirst(k -> k > log10(k*10.0^m)+1, 1:2m)-1) for m in 2:100]

# %%
?@nloops

# %%
(@macroexpand @nloops 3 i d -> 0:9 d -> y_d = i_d^2 begin
    y = 0
    @nexprs 3 d -> y += y_d
end)|> Base.remove_linenums!

# %%
@ntuple 10 i -> (i - 1) ^ 3

# %%
(@macroexpand @ntuple 10 i -> (i - 1) ^ 3) |> Base.remove_linenums!

# %%

max_stars_repo_path: 0017/Project Euler 30.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_path: 0017/Project Euler 30.jl

max_issues_repo_name: genkuroki/public

max_forks_repo_path: 0017/Project Euler 30.jl

max_forks_repo_name: genkuroki/public


text:
module Day7

export get_inputs, get_solution1, get_solution2

using Memoize

## Input getting
read_inputs(f_name) = make_the_rules(readlines(joinpath(@__DIR__, f_name)))

function get_inputs()
    test_input1 = read_inputs("test_input1.txt")
    test_input2 = read_inputs("test_input2.txt")
    test_output1 = 4
    test_output2 = 126
    data = read_inputs("input.txt")
    return (; test_input1, test_input2, test_output1, test_output2, data)
end

# Statement type for easy parsing of rules
struct Statement
    num::Int
    desc::String
end

const NO_OTHER = Statement(0, "no other")

function Statement(str)
    if str=="no other"
        return NO_OTHER
    else
        num, desc = split(str, " ", limit=2)
        num = parse(Int, num)
        return Statement(num, desc)
    end
end

num_bags(statement) = statement.num
description(statement) = statement.desc

## Helper functions
# Parse a single line of the rules
function parse_rule(line)
    head, body = split(line, r" bags contain ")
    body = split(body, r" (bag|bags)(, |\.)", keepempty=false)
    return String(head) => Statement.(body)
end

# Make the rule dict
make_the_rules(lines) = Dict(parse_rule.(lines))

## Solution functions
# Part 1
@memoize function contains_bag(rules, outer, inner)
    statements = rules[outer]
    if statements[1]==NO_OTHER
        return false
    elseif inner in description.(statements)
        return true
    else
        return any(contains_bag(rules, description(statement), inner) for statement in statements)
    end
end

get_solution1(rules; bag="shiny gold") = count(contains_bag(rules, head, bag) for head in keys(rules))

# Part 2
@memoize function count_bags(rules, bag)
    statements = rules[bag]
    if statements[1]==NO_OTHER
        return 0
    else
        return sum(num_bags(s) * (1 + count_bags(rules, description(s))) for s in statements)
    end
end

get_solution2(rules; bag="shiny gold") = count_bags(rules, bag)

## Bonus
# Pretty print bag structure
function show_bag(rules, bag="shiny gold", pre="", first=true)
    first && println(bag)
    statements = rules[bag]
    # statements = sort(rules[bag], by=x->count_bags(rules, x.desc), rev=true)
    if statements[1] == NO_OTHER
        return nothing
    end
    for i in 1:length(statements)-1
        statement = statements[i]
        println(pre * " ├╴ $(statement.desc) × $(statement.num)")
        show_bag(rules, statement.desc, pre * " │  ", false)
    end
    println(pre * " ╰╴ $(statements[end].desc) × $(statements[end].num)")
    show_bag(rules, statements[end].desc, pre * "    ", false)
end

end


text:
# sample sequences of ergms with different parameters' values from R package ergm
# and test the PseudoLikelihoodScoreDrivenERGM filter
using DrWatson
@quickactivate "ScoreDrivenExponentialRandomGraphs"
DrWatson.greet()

using ScoreDrivenERGM
import ScoreDrivenERGM:StaticNets,DynNets, Utilities, ErgmRcall
using LinearAlgebra
using Statistics
using StatsBase
using DataStructures
using JLD2
using StatsBase,CSV, RCall

model = DynNets.SdErgmPml(staticModel = StaticNets.NetModeErgmPml("edges +  gwesp(decay = 0.25,fixed = TRUE)", true), indTvPar = [false, true], scoreScalingType="FISH_D")

# load cleaned data
r_data_path = datadir("US_congr_covoting", "clean_data") * "\\Rollcall_VCERGM.RData"
R"load($r_data_path)"
R"obsMat_T = Rollcall$network "
@rget obsMat_T
A_T = [Int8.(m) for m in obsMat_T[2:end]]
obsT = DynNets.seq_of_obs_from_seq_of_mats(model, A_T)

#Single snapshot MLE estimation using MCMC in ergm R package can take a very long time. Load previously estimated data
# estParSS_T = [ErgmRcall.get_one_mle(m , model.staticModel.ergmTermsString) for m in obsMat_T]
JLD2.@load( datadir("US_congr_covoting", "estimates") * "\\SS_est_MCMC.jld2",estParSS_T)#, obsMat_T, changeStats_T, stats_T)

ENV["JULIA_DEBUG"] = ScoreDrivenERGM

N= [size(A)[1] for A in A_T]

res_est = DynNets.estimate_and_filter(model, N, obsT; show_trace = true)
fig, ax = DynNets.plot_filtered(model, N, res[3])

res_conf = DynNets.conf_bands_given_SD_estimates(model, N, obsT, DynNets.unrestrict_all_par(model, res.vEstSdResParAll), res.ftot_0, [[0.975, 0.025]]; indTvPar = model.indTvPar, offset=0, plotFlag=true, parUncMethod = "WHITE-MLE")

DynNets.plot_filtered(model, N, estParSS_T; lineType = ".", lineColor = "k", fig = res_conf.fig, ax=res_conf.ax, gridFlag=false, xval = dates[t0_train:tend_train])

estParStatic = ErgmRcall.get_one_mple(obsT, model.staticModel.ergmTermsString)

JLD2.@save( datadir("US_congr_covoting", "estimates")*"\\SD_est_$(DynNets.name(model))", res_conf, res_est, estParStatic)

max_stars_repo_name: domenicodigangi/SDERGM

max_issues_repo_name: domenicodigangi/SDERGM

max_forks_repo_name: domenicodigangi/SDERGM


text:
## Factorial of a number
function fac(n)
    if n == 1
        return 1
    else
        return (n * fac(n-1))
    end
end

for i in 1:10
    println(fac(i))
end


text:
#=
Find local bindings

Downstreams: completions.jl, goto.jl, datatip.jl, refactor.jl
=#

struct LocalBinding
    name::String
    verbatim::String
    span::UnitRange{Int64}
    line::Int
    expr::EXPR
end

struct LocalScope
    name::String
    verbatim::String
    span::UnitRange{Int64}
    line::Int
    children::Vector{Union{LocalBinding,LocalScope}}
    expr::EXPR
end

const LocalBS = Union{LocalBinding,LocalScope}

struct ActualLocalBinding
    name::String
    verbatim::String
    root::String
    line::Int
    locality::Float64
    expr::EXPR
end
function ActualLocalBinding(bs::LocalBS, root::String, line::Int, byteoffset::Int)
    locality = distance(line, byteoffset, bs.line, bs.span)
    return ActualLocalBinding(bs.name, bs.verbatim, root, bs.line, locality, bs.expr)
end

"""
    locals(text::String, line::Int, col::Int)::Vector{ActualLocalBinding}

Returns local bindings in `text`, while computing localities based on `line` and `col`.
"""
function locals(text::String, line::Int, col::Int)::Vector{ActualLocalBinding}
    expr = CSTParser.parse(text, true)
    traverse_expr!(expr)
    bindings = localbindings(expr, text)
    actual_localbindings(bindings, line, byteoffset(text, line, col))
end

function localbindings(expr, text, bindings = LocalBS[], pos = 1, line = 1)
    # binding
    bind = bindingof(expr)
    hs = hasscope(expr)
    if bind !== nothing && !hs
        verbatim = str_value_verbatim(bind, text, pos)
        range = pos:pos+expr.span
        push!(bindings, LocalBinding(bind.name, verbatim, range, line, expr))
    end

    if hs
        typof(expr) === CSTParser.Kw && return bindings

        # destructure multiple returns
        if ismultiplereturn(expr)
            for arg in expr
                # don't update `pos` & `line`, i.e.: treat all the multiple returns as same
                localbindings(arg, text, bindings, pos, line)
            end
        # properly detect the parameters of a method with where clause: https://github.com/JunoLab/Juno.jl/issues/404
        elseif iswhereclause(expr)
            for arg in expr
                localbindings(arg, text, bindings, pos, line)
                line += countlines(arg, text, pos)
                pos += arg.fullspan
            end
        else
            # find local binds in a scope
            # calculate fields for `LocalScope` first
            verbatim = str_value_verbatim(expr, text, pos)
            range = pos:pos+expr.span
            name = bind === nothing ? "" : bind.name

            children = LocalBS[]
            for arg in expr
                localbindings(arg, text, children, pos, line)
                line += countlines(arg, text, pos)
                pos += arg.fullspan
            end

            push!(bindings, LocalScope(name, verbatim, range, line, children, expr))
        end

        return bindings
    end

    # look for more local bindings if exists
    for arg in expr
        localbindings(arg, text, bindings, pos, line)
        line += countlines(arg, text, pos)
        pos += arg.fullspan
    end

    return bindings
end

function byteoffset(text, line, col)
    byteoffset = 1
    current_line = 1
    current_char = 0
    for c in text
        if line == current_line
            current_char += 1
            c === '\n' && break
        end
        current_char == col && break
        byteoffset += @static VERSION >= v"1.1" ? ncodeunits(c) : ncodeunits(string(c))
        c === '\n' && (current_line += 1)
    end
    byteoffset
end

function actual_localbindings(bindings, line, byteoffset)
    actual_bindings = _actual_localbindings(bindings, line, byteoffset)

    filter!(b -> !isempty(b.name), actual_bindings)
    sort!(actual_bindings, lt = (b1, b2) -> b1.locality < b2.locality)
    return @static VERSION ≥ v"1.1" ? unique!(b->b.name, actual_bindings) : unique(b->b.name, actual_bindings)
end
function _actual_localbindings(bindings, line, byteoffset, root = "", actual_bindings = ActualLocalBinding[])
    for bind in bindings
        push!(actual_bindings, ActualLocalBinding(bind, root, line, byteoffset))
        if bind isa LocalScope && byteoffset in bind.span
            _actual_localbindings(bind.children, line, byteoffset, bind.name, actual_bindings)
        end
    end

    return actual_bindings
end

function distance(line, byteoffset, defline, defspan)
    abslinediff = abs(line - defline)
    absbytediff = abs(byteoffset - defspan[1]) # tiebreaker for bindings on the same line
    diff = if byteoffset in defspan
        Inf
    elseif line < defline
        (defline - line)*10 # bindings defined *after* the current line have a lower priority
    else
        line - defline
    end
    diff + absbytediff*1e-6
end


text:
# the core loop for several CRC calculations.  it is used by CRC32,
# whose zlib implementation in C runs ~2.5 times faster than this
# code.

function generic_loop{D<:Unsigned, A<:Unsigned
                      }(::Type{D}, remainder::A, data, table::Vector{A}, word_size)
    for word::D in data
        word $= convert(D, remainder)
        remainder = (remainder >>> word_size) $ table[1 + word]
    end
    remainder
end

# try accessing the lowest byte directly to avoid the shift (slow)

immutable Bytes
    a::Uint8
    b::Uint8
    c::Uint8
    d::Uint8
end

function struct_loop(::Type{Uint8}, remainder::Uint64, data::Vector{Uint8}, table::Vector{Uint64}, ignored)
    r::Uint64 = remainder
    b = reinterpret(Bytes, r)
    for i in 1:length(data)
        @inbounds word::Uint8 = data[i]
        word $= b.a
        r = (r >>> 8) $ table[1 + word]
        b = reinterpret(Bytes, r)
    end
    return r
end

# try accessing via aliased arrays

function array_loop(::Type{Uint8}, remainder::Uint64, data::Vector{Uint8}, table::Vector{Uint64}, ignored)
    r = Array(Uint64, 1)
    @inbounds r[1] = remainder
    b = reinterpret(Uint8, r)
    for i in 1:length(data)
        @inbounds word::Uint8 = data[i]
        @inbounds word $= b[1]
        @inbounds r[1] = (r[1] >>> 8) $ table[1 + word]  # @inbounds no help here
    end
    r[1]
end

# according to my notes this was slightly (~15%) faster than
# generic_loop.  repeated here, however, it seems to be identical to
# the generic case.

function byte_loop{A<:Unsigned
                   }(::Type{Uint8}, remainder::A, data::Vector{Uint8}, table::Vector{A}, ignored)
    for i in 1:length(data)
        @inbounds word::Uint8 = data[i]
        word $= convert(Uint8, remainder)
        remainder = (remainder >>> 8) $ table[1 + word]  # @inbounds no help here
    end
    remainder
end

# timing with a standard data set (the table is "meaningless")

srand(0)  # repeatable results
TABLE = Uint64[rand(Uint64) for _ in 1:256]
DATA = Uint8[rand(Uint8) for _ in 1:10_000_000]
#KNOWN = 17880998858342525629  # 1 million
KNOWN = 3923249127266084647  # 10 million
#KNOWN = 10790464939549839568  # 100 million

function run_std(loop)
    loop(Uint8, zero(Uint64), DATA[1:10], TABLE, 8)  # warmup
    @printf("%15s: ", loop)
    @time (result = loop(Uint8, zero(Uint64), DATA, TABLE, 8))
    @assert result == KNOWN
end

for loop in (generic_loop, struct_loop, byte_loop, array_loop)
    run_std(loop)
end

#  generic_loop: elapsed time: 0.023491408 seconds (16 bytes allocated)
#   struct_loop: elapsed time: 2.027252231 seconds (1680000168 bytes allocated)
#     byte_loop: elapsed time: 0.02304929 seconds (16 bytes allocated)
#    array_loop: elapsed time: 0.045178428 seconds (168 bytes allocated)


text:
# using GLVisualize, GeometryTypes, Colors, GLAbstraction
# using GLAbstraction
# using FixedSizeDictionaries
#
# window = glscreen()
#
# s = """wooot"""
# dict = Dict()
# color = map(collect(s)) do char
#     c = get!(dict, char, RGBA{Float32}(rand(), rand(), rand(), 1))
#     c::RGBA{Float32}
# end
#
# # _view and visualize it!
# _view(visualize(s,
#     model=translationmatrix(Vec3f0(0,600,0)), # move this up, since the text starts at 0 and goes down from there
#     color=color,
#     rotation=fill(Vec3f0(0,0,1), length(s))
# ), window)
#
# @async renderloop(window)
#
#
# rt = GLVisualize.RichText(
#     GLVisualize.Text(renderlist(window)[1])
# )
#
# gap = 200
# for i=1:10
#     x = string(rand())
#     insert!(rt, x, rt.cursor, Dict(
#         :start_position=>Point2f0(0, (i-1)*gap),
#     ))
#     rt.cursor += length(x)
# end
# gap = 300
#
# for i=1:10
#     x = string(rand())
#     insert!(rt, x, rt.cursor, Dict(
#         :color => RGBA{Float32}(1,0,0,1),
#         :start_position=>Point2f0((i-1)*gap, 0),
#         :rotation => Vec3f0(0,0,1)
#     ))
#     rt.cursor += length(x)
# end
# gap


text:
"""
Modia module with block component models (inspired from Modelica Standard Library).

* Developer: Hilding Elmqvist, Mogram AB, Martin Otter, DLR
* Copyright (c) 2016-2021: Hilding Elmqvist, Martin Otter
* License: MIT (expat)

"""
#module Blocks

using Modia

#export Gain, FirstOrder, Feedback, PI, Step, Ramp # Sine, Switch, MIMO

# Single-Input-Single-Output continuous control block
SISO = Model(
    u = input,
    y = output
)

# Gain
Gain = SISO | Model(
    k = 1, # (info = "Gain")
    equations = :[
        y = k*u ]
)

# First-order transfer function block (= 1 pole)
FirstOrder = SISO | Model(
    k = 1.0,
    T = 1.0*u"s",
    x = Var(init=0.0),
    equations = :[
        der(x) = (k * u - x) / T
        y = x ]
)

# Output difference between commanded and feedback input
Feedback = Model(
    u1 = input | info"Input 1",
    u2 = input | info"Input 2",
    y = output | info"Output signal",
    equations = :[
        y = u1 - u2 ]
)

# Proportional-Integral controller
PI = SISO | Model(
    k = 1.0, # (info = "Gain")
    T = 1.0u"s", # (min = 1E-10, info = "Time Constant (T>0 required)")
    x = Var(init=0.0),
    equations = :[
        der(x) = u / T
        y = k * (x + u) ]
)

# Single-Output continuous control block
SO = Model(
    y = output
)

# Base class for a continuous signal source
SignalSource = SO | Model(
    offset = 0.0, # info = "Offset of output signal y"
    startTime = 0.0*u"s" # info = "Output y = offset for time < startTime")
)

# Step signal
Step = SignalSource | Model(
    height = 1.0,
    equations = :[
        y = offset + (time < startTime ? 0*height : height) ]  # 0*height is needed, if height has a unit
)

# Ramp signal
Ramp = SignalSource | Model(
    height   = 1.0,
    duration = 2.0u"s",
    equations = :[
        y = offset + (time < startTime ? 0.0*height :   # 0*height is needed, if height has a unit
                     (time < startTime + duration ? (time - startTime)*height/duration :
                                                    height)) ]
)

# Linear state space system
StateSpace = Model(
    A = parameter | fill(0.0,0,0),
    B = parameter | fill(0.0,0,0),
    C = parameter | fill(0.0,0,0),
    D = parameter | fill(0.0,0,0),
    u = input,
    y = output,
    x = Var(init = zeros(0)),
    equations = :[
        der(x) = A*x + B*u
             y = C*x + D*u
    ]
)

# -------------------------------------------------------

#=
# Sinusoidal signal
@model Sine begin
    amplitude = Parameter(1.0, info = "Amplitude of sine wave")
    freqHz = Parameter(1.0,info = "Frequency of sine wave")
    phase = Parameter(0.0, info = "Phase of sine wave")
    offset = Parameter(0.0, info = "Offset of output signal y")
    startTime = Parameter(0.0, info = "Output y = offset for time < startTime")
    SO()
equations
    y = offset + if time < startTime;  0 else amplitude*sin(2*pi*freqHz*(time - startTime) + phase) end
end

@model Sine2 begin
    # Generate sine signal
    amplitude = Parameter(1.0, info = "Amplitude of sine wave")
    freqHz = Parameter(1.0,info = "Frequency of sine wave")
    phase = Parameter(0.0, info = "Phase of sine wave")
    SignalSource()
equations
    y = offset + if time < startTime;  0 else amplitude * sin(2 * pi * freqHz * (time - startTime) + phase) end
end

# Switch
@model Switch begin
    sw = Input(Boolean(info = "Switch position (if `true`, use `u1`, else use `u2`)"))
    u1 = Input(info = "Input 1")
    u2 = Input(info = "Input 2")
    y = Output(info = "Output signal")
equations
    y = if sw; u1 else u2 end
end

# ABCD model
@model ABCD(A = -1, B = 1, C = 1, D = 0) begin
    u = Input(info = "Input signal"); y = Output(info = "Output signal")
    x = Float(start = 0)
equations
    der(x) = A * x + B * u
    y = C * x + D * u
end
=#
#end