text:
isdefined(Base, :__precompile__) && __precompile__()

module SIUnits

    using Compat; import Compat.String

    import Base: ==, +, -, *, /, .+, .-, .*, ./, //, ^
    import Base: promote_rule, promote_type, convert, show, mod

    immutable SIQuantity{T<:Number,m,kg,s,A,K,mol,cd,rad,sr} <: Number
        val::T
    end

    typealias UnitQuantity{T} SIQuantity{T,0,0,0,0,0,0,0,0,0}

    SIQuantity{T<:Number}(x::T) = UnitQuantity{T}(x)

    immutable SIUnit{m,kg,s,A,K,mol,cd,rad,sr} <: Number
    end

    if !isdefined(Base, :UnitRange)
        const Range = Ranges # Deprecations introduced early in the 0.3 cycle
        const UnitRange = Range1
    end

    abstract SIRanges{T,m,kg,s,A,K,mol,cd,rad,sr} <: Range{SIQuantity{T,m,kg,s,A,K,mol,cd,rad,sr}}

    immutable SIRange{R<:Range,T<:Real,m,kg,s,A,K,mol,cd,rad,sr} <: SIRanges{T,m,kg,s,A,K,mol,cd,rad,sr}
        val::R
    end

    typealias UnitTuple NTuple{9,Int}

    unit{T,m,kg,s,A,K,mol,cd,rad,sr}(x::SIRanges{T,m,kg,s,A,K,mol,cd,rad,sr}) = SIUnit{m,kg,s,A,K,mol,cd,rad,sr}()
    quantity{T,m,kg,s,A,K,mol,cd,rad,sr}(x::SIRanges{T,m,kg,s,A,K,mol,cd,rad,sr}) = SIQuantity{T,m,kg,s,A,K,mol,cd,rad,sr}

    import Base: length, getindex, unsafe_getindex, next, float, show, start, step, last, done, first, eltype, one, zero

    one(x::SIQuantity) = one(x.val)
    one{T,m,kg,s,A,K,mol,cd,rad,sr}(::Type{SIQuantity{T,m,kg,s,A,K,mol,cd,rad,sr}}) = one(T)
    zero(x::SIQuantity) = zero(x.val) * unit(x)
    zero{T,m,kg,s,A,K,mol,cd,rad,sr}(::Type{SIQuantity{T,m,kg,s,A,K,mol,cd,rad,sr}}) = zero(T) * SIUnit{m,kg,s,A,K,mol,cd,rad,sr}()

    # This is all nessecary because SIQuanity{T<:Real} !<: Real
    show(io::IO, x::SIRanges) = (show(io, x.val); show(io,unit(x)))
    function show(io::IO, r::SIRange)
        if step(r) == zero(quantity(r))
            print(io, "SIRange(",first(r),",",step(r),",",length(r),")")
        else
            print(io, first(r),':',step(r),':',last(r))
        end
    end
    show{T<:UnitRange}(io::IO, r::SIRange{T}) = print(io, first(r),':',last(r))
    getindex(r::SIRanges,i::Integer) = (quantity(r)(getindex(r.val,i)))
    getindex{T<:SIRanges}(r::T,i::Range) = T(getindex(r.val,i))
    unsafe_getindex(r::SIRanges,i::Integer) = (quantity(r)(unsafe_getindex(r.val,i)))
    unsafe_getindex{T<:SIRanges}(r::T,i::Range) = T(unsafe_getindex(r.val,i))
    function next(r::SIRanges, i)
        v, j = next(r.val,i)
        to_q(quantity(r),v), j
    end
    length(r::SIRanges) = length(r.val)
    start(r::SIRanges) = start(r.val)
    done(r::SIRanges,i) = done(r.val,i)
    eltype(r::SIRanges) = quantity(r)

    for func in (:first,:step,:last)
        @eval $(func)(r::SIRanges) = to_q(quantity(r),$(func)(r.val))
    end
    # Forward some linear range transformations to the wrapped range
    rangequantity{R<:Range}(::Type{R},tup::UnitTuple) = SIRange{R,eltype(R),tup[1],tup[2],tup[3],tup[4],tup[5],tup[6],tup[7],tup[8],tup[9]}
    for func in (VERSION < v"0.3-" ? (:+, :-) : (:.+, :.-)) # version 0.3 swaps fallbacks
        @eval $(func){T,S,m,kg,s,A,K,mol,cd,rad,sr}(x::SIRanges{T,m,kg,s,A,K,mol,cd,rad,sr}, y::SIQuantity{S,m,kg,s,A,K,mol,cd,rad,sr}) = (val = $(func)(x.val, y.val); SIRange{typeof(val),eltype(val),m,kg,s,A,K,mol,cd,rad,sr}(val))
        @eval $(func){T,S,m,kg,s,A,K,mol,cd,rad,sr}(x::SIQuantity{T,m,kg,s,A,K,mol,cd,rad,sr}, y::SIRanges{S,m,kg,s,A,K,mol,cd,rad,sr}) = (val = $(func)(x.val, y.val); SIRange{typeof(val),eltype(val),m,kg,s,A,K,mol,cd,rad,sr}(val))
    end
    ./(x::SIRanges, y::SIQuantity) = (val = ./(x.val, y.val); rangequantity(typeof(val),tup(x)-tup(y))(val))
    .*(x::SIRanges, y::SIQuantity) = (val = .*(x.val, y.val); rangequantity(typeof(val),tup(x)+tup(y))(val))
    .*(x::SIQuantity, y::SIRanges) = (val = .*(x.val, y.val); rangequantity(typeof(val),tup(x)+tup(y))(val))
    # Version 0.2 assumes all Ranges have start and len fields in ==, and
    # the fallback in 0.3 needlessly iterates through all values
    ==(r::SIRanges, s::SIRanges) = r.val == s.val && tup(r) == tup(s)
    ==(s::SIRanges, r::Range) = s.val == r && tup(s) == (0,0,0,0,0,0,0,0,0)
    ==(r::Range, s::SIRanges) = r == s.val && tup(s) == (0,0,0,0,0,0,0,0,0)

    tup2u(tup) = SIUnit{tup[1],tup[2],tup[3],tup[4],tup[5],tup[6],tup[7],tup[8],tup[9]}
    quantity(T::Type,tup::UnitTuple) = quantity(T,tup2u(tup)())
    -(tup::UnitTuple) = (-tup[1],-tup[2],-tup[3],-tup[4],-tup[5],-tup[6],-tup[7],-tup[8],-tup[9])

    for op in (:-,:*,:+)
        @eval function $(op)(tup1::UnitTuple,tup2::UnitTuple)
            ($(op)(tup1[1],tup2[1]),$(op)(tup1[2],tup2[2]),$(op)(tup1[3],tup2[3]),$(op)(tup1[4],tup2[4]),$(op)(tup1[5],tup2[5]),
                $(op)(tup1[6],tup2[6]),$(op)(tup1[7],tup2[7]),$(op)(tup1[8],tup2[8]),$(op)(tup1[9],tup2[9]))
        end
    end

    export quantity, @quantity

    function quantity{S}(T,quant::SIQuantity{S})
        quant.val == one(S) || error("Quantity value must be unity!")
        quantity(T,unit(quant))
    end
    quantity{m,kg,s,A,K,mol,cd,rad,sr}(T::(@compat Union{Type,TypeVar}),unit::SIUnit{m,kg,s,A,K,mol,cd,rad,sr}) = SIQuantity{T,m,kg,s,A,K,mol,cd,rad,sr}

    tup{m,kg,s,A,K,mol,cd,rad,sr}(u::SIUnit{m,kg,s,A,K,mol,cd,rad,sr}) = (m,kg,s,A,K,mol,cd,rad,sr)
    tup{T,m,kg,s,A,K,mol,cd,rad,sr}(u::SIQuantity{T,m,kg,s,A,K,mol,cd,rad,sr}) = (m,kg,s,A,K,mol,cd,rad,sr)
    tup{T,m,kg,s,A,K,mol,cd,rad,sr}(u::SIRanges{T,m,kg,s,A,K,mol,cd,rad,sr}) = (m,kg,s,A,K,mol,cd,rad,sr)

    macro quantity(expr,unit)
        esc(:(SIUnits.SIQuantity{$expr,SIUnits.tup($unit)...}))
    end

    # Irrationals propagate through units. Fancy!!
    promote_rule{sym,m,kg,s,A,K,mol,cd,rad,sr}(x::Type{Irrational{sym}},y::Type{SIUnit{m,kg,s,A,K,mol,cd,rad,sr}}) =
        SIQuantity{Irrational{sym},m,kg,s,A,K,mol,cd,rad,sr}
    promote_rule{sym,T,m,kg,s,A,K,mol,cd,rad,sr}(x::Type{Irrational{sym}},y::Type{SIQuantity{T,m,kg,s,A,K,mol,cd,rad,sr}}) =
        SIQuantity{promote_type(Irrational{sym},T)}

    promote_rule{T,S,mS,kgS,sS,AS,KS,molS,cdS,radS,srS,mT,kgT,sT,AT,KT,molT,cdT,radT,srT}(
        A::Type{SIQuantity{T,mT,kgT,sT,AT,KT,molT,cdT,radT,srT}},B::Type{SIQuantity{S,mS,kgS,sS,AS,KS,molS,cdS,radS,srS}}) = SIQuantity{promote_type(T,S)}
    promote_rule{T,mS,kgS,sS,AS,KS,molS,cdS,radS,srS,mT,kgT,sT,AT,KT,molT,cdT,radT,srT}(
        A::Type{SIQuantity{T,mT,kgT,sT,AT,KT,molT,cdT,radT,srT}},B::Type{SIUnit{mS,kgS,sS,AS,KS,molS,cdS,radS,srS}}) = SIQuantity{T}
    promote_rule{S,m,kg,s,A,K,mol,cd,rad,sr}(x::Type{Bool},y::Type{SIQuantity{S,m,kg,s,A,K,mol,cd,rad,sr}}) = SIQuantity{promote_type(Bool,S)}
    promote_rule{m,kg,s,A,K,mol,cd,rad,sr}(x::Type{Bool},y::Type{SIUnit{m,kg,s,A,K,mol,cd,rad,sr}}) = SIQuantity{Bool}
    promote_rule{T<:Number,S,m,kg,s,A,K,mol,cd,rad,sr}(x::Type{T},y::Type{SIQuantity{S,m,kg,s,A,K,mol,cd,rad,sr}}) = SIQuantity{promote_type(T,S)}
    promote_rule{T<:Number,m,kg,s,A,K,mol,cd,rad,sr}(x::Type{T},y::Type{SIUnit{m,kg,s,A,K,mol,cd,rad,sr}}) = SIQuantity{T}

    # One unspecified, units, one concrete (unspecified occurs as the promotion result from the rules above)
    promote_rule{T,S,m,kg,s,A,K,mol,cd,rad,sr}(x::Type{SIQuantity{T}},y::Type{SIQuantity{S,m,kg,s,A,K,mol,cd,rad,sr}}) = SIQuantity{promote_type(T,S)}

    # Unlike most other types, the promotion of two identitical SIQuantities is
    # not that type itself. As such, the promote_type behavior itself must be
    # overridden. C.f. https://github.com/Keno/SIUnits.jl/issues/27
    promote_type{T,m,kg,s,A,K,mol,cd,rad,sr}(::Type{SIQuantity{T,m,kg,s,A,K,mol,cd,rad,sr}}, ::Type{SIQuantity{T,m,kg,s,A,K,mol,cd,rad,sr}}) = SIQuantity{T}

    if VERSION >= v"0.4-dev"
        eval(quote
            convert{T}(::Type{SIQuantity{T}},x::Dates.Period) = error("Conversion from Period to SIQuantity not defined")
        end)
    end
    convert{T<:Number,m,kg,s,A,K,mol,cd,rad,sr}(::Type{SIQuantity{T}},x::SIUnit{m,kg,s,A,K,mol,cd,rad,sr}) = SIQuantity{T,m,kg,s,A,K,mol,cd,rad,sr}(one(T))
    convert{T<:Number}(::Type{SIQuantity{T}},x::T) = UnitQuantity{T}(x)
    convert{T<:Number,S<:Number}(::Type{SIQuantity{T}},x::S) = convert(SIQuantity{T},convert(T,x))
    convert{T<:Number}(::Type{SIQuantity{T}},x::SIQuantity{T}) = x
    convert{T<:Number,S,m,kg,s,A,K,mol,cd,rad,sr}(::Type{SIQuantity{T}},x::SIQuantity{S,m,kg,s,A,K,mol,cd,rad,sr}) = SIQuantity{T,m,kg,s,A,K,mol,cd,rad,sr}(convert(T,x.val))

    to_q{T,m,kg,s,A,K,mol,cd,rad,sr}(::Type{SIQuantity{T,m,kg,s,A,K,mol,cd,rad,sr}},val::T) = (0 == m == kg == s == A == K == mol == cd == rad == sr) ? val : SIQuantity{T,m,kg,s,A,K,mol,cd,rad,sr}(val)
    convert{T,m,kg,s,A,K,mol,cd,rad,sr}(::Type{SIQuantity{T,m,kg,s,A,K,mol,cd,rad,sr}},val::Number) = (SIQuantity{T,m,kg,s,A,K,mol,cd,rad,sr}(convert(T,val)))
    function convert{T,S,mS,kgS,sS,AS,KS,molS,cdS,radS,srS,mT,kgT,sT,AT,KT,molT,cdT,radT,srT}(::Type{SIQuantity{T,mT,kgT,sT,AT,KT,molT,cdT,radT,srT}},val::SIQuantity{S,mS,kgS,sS,AS,KS,molS,cdS,radS,srS})
        if mS != mT || kgS != kgT || sS != sT || AS != AT || KS != KT || molS != molT || cdS != cdT || radS != radT || srS != srT
            error("Dimension mismatch in convert. Attempted to convert a ($(repr(SIUnit{mS,kgS,sS,AS,KS,molS,cdS,radS,srS}))) to ($(repr(SIUnit{mT,kgT,sT,AT,KT,molT,cdT,radT,srT})))")
        end
        SIQuantity{T,mT,kgT,sT,AT,KT,molT,cdT,radT,srT}(convert(T,val.val))
    end

    for op in (:/,://)

        @eval function ($op){T}(x::Number,y::SIQuantity{T})
            val = ($op)(x,y.val)
            to_q(quantity(typeof(val),-tup(y)),val)
        end

        @eval function ($op)(x::SIQuantity,y::SIQuantity)
            val = $(op)(x.val,y.val)
            to_q(quantity(typeof(val),tup(x)-tup(y)),val)
        end

        @eval $(op)(x::SIUnit,y::SIUnit) = tup2u(tup(x)-tup(y))()
        @eval $(op){T}(x::SIQuantity{T},y::SIUnit) = to_q(quantity(T,tup(unit(x))-tup(y)),x.val)
        @eval $(op){T}(x::SIUnit,y::SIQuantity{T}) = to_q(quantity(T,tup(x)-tup(unit(y))),($op)(1,y.val))

        @eval $(op)(x::Number,y::SIUnit) = x*tup2u(-tup(y))()
    end

    inv(y::SIUnit) = tup2u(-tup(y))()

    function +{T,S,m,kg,s,A,K,mol,cd,rad,sr}(
        x::SIQuantity{T,m,kg,s,A,K,mol,cd,rad,sr},y::SIQuantity{S,m,kg,s,A,K,mol,cd,rad,sr})
        val = x.val+y.val
        SIQuantity{typeof(val),m,kg,s,A,K,mol,cd,rad,sr}(val)
    end

    function -{T,S,m,kg,s,A,K,mol,cd,rad,sr}(
        x::SIQuantity{T,m,kg,s,A,K,mol,cd,rad,sr},y::SIQuantity{S,m,kg,s,A,K,mol,cd,rad,sr})
        val = x.val-y.val
        SIQuantity{typeof(val),m,kg,s,A,K,mol,cd,rad,sr}(val)
    end

    function -{T,S,mS,kgS,sS,AS,KS,molS,cdS,radS,srS,mT,kgT,sT,AT,KT,molT,cdT,radT,srT}(
        x::SIQuantity{T,mT,kgT,sT,AT,KT,molT,cdT,radT,srT},y::SIQuantity{S,mS,kgS,sS,AS,KS,molS,cdS,radS,srS})
        error("Unit mismatch. Got ($(repr(unit(x)))) - ($(repr(unit(y))))")
    end

    function +{T,S,mS,kgS,sS,AS,KS,molS,cdS,radS,srS,mT,kgT,sT,AT,KT,molT,cdT,radT,srT}(
        x::SIQuantity{T,mT,kgT,sT,AT,KT,molT,cdT,radT,srT},y::SIQuantity{S,mS,kgS,sS,AS,KS,molS,cdS,radS,srS})
        error("Unit mismatch. Got ($(repr(unit(x)))) + ($(repr(unit(y))))")
    end

    function -{T,m,kg,s,A,K,mol,cd,rad,sr}(x::SIQuantity{T,m,kg,s,A,K,mol,cd,rad,sr})
        val = -(x.val)
        SIQuantity{typeof(val),m,kg,s,A,K,mol,cd,rad,sr}(val)
    end

    function ^{T,m,kg,s,A,K,mol,cd,rad,sr}(
        x::SIQuantity{T,m,kg,s,A,K,mol,cd,rad,sr},i::Integer)
        if i == 0
            return one(T)
        end
        val = x.val^i
        SIQuantity{typeof(val),m*i,kg*i,s*i,A*i,K*i,mol*i,cd*i,rad*i,sr*i}(val)
    end

    function ^{T,m,kg,s,A,K,mol,cd,rad,sr}(
        x::SIQuantity{T,m,kg,s,A,K,mol,cd,rad,sr},r::Rational)
        if r == 0
            return one(T)
        end
        val = x.val^r
        SIQuantity{typeof(val),convert(Int,m*r),convert(Int,kg*r),convert(Int,s*r),convert(Int,A*r),
        convert(Int,K*r),convert(Int,mol*r),convert(Int,cd*r),convert(Int,rad*r),convert(Int,sr*r)}(val)
    end

    ^{T,m,kg,s,A,K,mol,cd,rad,sr}(x::SIQuantity{T,m,kg,s,A,K,mol,cd,rad,sr},r::AbstractFloat) = x^rationalize(r)

    function ^{T,S,mS,kgS,sS,AS,KS,molS,cdS,radS,srS,mT,kgT,sT,AT,KT,molT,cdT,radT,srT}(
        x::SIQuantity{T,mT,kgT,sT,AT,KT,molT,cdT,radT,srT},y::SIQuantity{S,mS,kgS,sS,AS,KS,molS,cdS,radS,srS})
        error("Can not raise a number to a unitful quantity. Got ($(repr(unit(x))))^($(repr(unit(y))))")
    end

    ^{T,S,m,kg,s,A,K,mol,cd,rad,sr}(x::SIQuantity{T,m,kg,s,A,K,mol,cd,rad,sr},y::SIQuantity{S,0,0,0,0,0,0,0,0}) = x.val^(y.val)

    ==(x::SIQuantity,y::SIQuantity) = (tup(x) == tup(y)) && (x.val == y.val)
    =={T}(x::SIQuantity{T},y::SIUnit) = (tup(x) == tup(y)) && (x.val == one(T))
    =={T}(x::SIUnit,y::SIQuantity{T}) = (tup(x) == tup(y)) && (one(T) == y.val)
    ==(x::SIUnit,y::SIUnit) = tup(x) == tup(y)

    import Base: sqrt, abs, abs2, colon, isless, isfinite, isreal, real, imag, isnan

    function colon{T,S,X,m,kg,s,A,K,mol,cd,rad,sr}(start::SIQuantity{T,m,kg,s,A,K,mol,cd,rad,sr},step::SIQuantity{S,m,kg,s,A,K,mol,cd,rad,sr},stop::SIQuantity{X,m,kg,s,A,K,mol,cd,rad,sr})
        val = colon(start.val,step.val,stop.val)
        SIRange{typeof(val),eltype(val),m,kg,s,A,K,mol,cd,rad,sr}(val)
    end

    function colon{T,S,m,kg,s,A,K,mol,cd,rad,sr}(start::SIQuantity{T,m,kg,s,A,K,mol,cd,rad,sr},stop::SIQuantity{S,m,kg,s,A,K,mol,cd,rad,sr})
        val = colon(start.val,stop.val)
        SIRange{typeof(val),eltype(val),m,kg,s,A,K,mol,cd,rad,sr}(val)
    end

    function sqrt{T,m,kg,s,A,K,mol,cd,rad,sr}(x::SIQuantity{T,m,kg,s,A,K,mol,cd,rad,sr})
        val = sqrt(x.val)
        SIQuantity{typeof(val),convert(Int,m/2),convert(Int,kg/2),convert(Int,s/2),convert(Int,A/2),
        convert(Int,K/2),convert(Int,mol/2),convert(Int,cd/2),convert(Int,rad/2),convert(Int,rad/2)}(val)
    end

    function abs{T,m,kg,s,A,K,mol,cd,rad,sr}(x::SIQuantity{T,m,kg,s,A,K,mol,cd,rad,sr})
        SIQuantity{T,m,kg,s,A,K,mol,cd,rad,sr}(abs(x.val))
    end

    function abs2{T,m,kg,s,A,K,mol,cd,rad,sr}(x::SIQuantity{T,m,kg,s,A,K,mol,cd,rad,sr})
        val = abs2(x.val)
        SIQuantity{typeof(val),2*m,2*kg,2*s,2*A,2*K,2*mol,2*cd,2*rad,2*sr}(val)
    end

    function isfinite{T,m,kg,s,A,K,mol,cd,rad,sr}(x::SIQuantity{T,m,kg,s,A,K,mol,cd,rad,sr})
        isfinite(x.val)
    end

    isnan(x::SIQuantity) = isnan(x.val)
    isreal(x::SIQuantity) = isreal(x.val)
    real(x::SIQuantity) = typeof(x)(real(x.val))
    imag(x::SIQuantity) = typeof(x)(imag(x.val))

    function isless{T}(x::SIQuantity{T,0,0,0,0,0,0,0,0,0}, y::SIQuantity{T,0,0,0,0,0,0,0,0,0})
        return isless(x.val,y.val)
    end
    function isless{T,S}(x::SIQuantity{T,0,0,0,0,0,0,0,0,0}, y::SIQuantity{S,0,0,0,0,0,0,0,0,0})
        return isless(x.val,y.val)
    end
    function isless{T}(x::SIQuantity{T,0,0,0,0,0,0,0,0,0}, y::Real)
        return isless(x.val,y)
    end
    function isless{T}(x::Real, y::SIQuantity{T,0,0,0,0,0,0,0,0,0})
        return isless(x,y.val)
    end
    function isless{T,S,mT,kgT,sT,AT,KT,molT,cdT,radT,srT}(
        x::SIQuantity{T,mT,kgT,sT,AT,KT,molT,cdT,radT,srT},y::SIQuantity{S,mT,kgT,sT,AT,KT,molT,cdT,radT,srT})
        return isless(x.val,y.val)
    end

    function mod{T,S,mS,kgS,sS,AS,KS,molS,cdS,radS,srS,mT,kgT,sT,AT,KT,molT,cdT,radT,srT}(
        x::SIQuantity{T,mT,kgT,sT,AT,KT,molT,cdT,radT,srT},y::SIQuantity{S,mS,kgS,sS,AS,KS,molS,cdS,radS,srS})
        error("Unit mismatch. Got mod($(repr(unit(x))),$(repr(unit(y))))")
    end

    function mod{T,S,m,kg,s,A,K,mol,cd,rad,sr}(x::SIQuantity{T,m,kg,s,A,K,mol,cd,rad,sr},y::SIQuantity{S,m,kg,s,A,K,mol,cd,rad,sr})
        val = mod(x.val,y.val)
        SIQuantity{typeof(val),m,kg,s,A,K,mol,cd,rad,sr}(val)
    end

    import Base: sin, cos, tan, cot, sec, csc
    for func in (:sin,:cos,:tan,:cot,:sec,:csc)
        @eval $func{T}(θ::SIQuantity{T,0,0,0,0,0,0,0,1,0}) = $func(θ.val)
    end

    # Forwarding methods that do not affect units
    import Base: conj
    conj(x::SIQuantity) = typeof(x)(conj(x.val))

    float(x::SIQuantity) = float(x.val)
    if VERSION >= v"0.4.0"
        Base.Float64(x::SIQuantity) = Float64(x.val)
        Base.Int(x::SIQuantity) = Int(x.val)
    else
        Base.float64(x::SIQuantity) = float64(x.val)
        Base.int(x::SIQuantity) = int(x.val)
    end

    *(x::SIUnit,y::SIUnit) = tup2u(tup(x)+tup(y))()
    *{T}(x::SIUnit,y::SIQuantity{T}) = to_q(quantity(T,tup(y)+tup(x)),y.val)
    *{T}(x::SIQuantity{T},y::SIUnit) = to_q(quantity(T,tup(y)+tup(x)),x.val)
    *(x::Irrational,y::SIUnit) = quantity(typeof(x),y)(x)
    function *(x::SIQuantity,y::SIQuantity)
        ret = x.val * y.val
        to_q(quantity(typeof(ret),tup(x)+tup(y)),ret)
    end

    function ^{m,kg,s,A,K,mol,cd,rad,sr}(
        x::SIUnit{m,kg,s,A,K,mol,cd,rad,sr},i::Integer)
        SIUnit{m*i,kg*i,s*i,A*i,K*i,mol*i,cd*i,rad*i,sr*i}()
    end

    unit{T,m,kg,s,A,K,mol,cd,rad,sr}(x::SIQuantity{T,m,kg,s,A,K,mol,cd,rad,sr}) = SIUnit{m,kg,s,A,K,mol,cd,rad,sr}()

    export SIPrefix, Meter, KiloGram, Second, Ampere, Kelvin, Mole, Candela, Radian, Steradian, Kilo, Mega, Giga,
        Tera, Peta, Exa, Zetta, Centi, Milli, Micro, Nano, Pico, Femto, Atto, Zepto, Yocto,
        Gram, Joule, Coulomb, Volt, Farad, Newton, Ohm, CentiMeter, Siemens, Hertz, Watt, Pascal

    const SIPrefix  = SIUnit{0,0,0,0,0,0,0,0,0}()
    const Meter     = SIUnit{1,0,0,0,0,0,0,0,0}()
    const KiloGram  = SIUnit{0,1,0,0,0,0,0,0,0}()
    const Second    = SIUnit{0,0,1,0,0,0,0,0,0}()
    const Ampere    = SIUnit{0,0,0,1,0,0,0,0,0}()
    const Kelvin    = SIUnit{0,0,0,0,1,0,0,0,0}()
    const Mole      = SIUnit{0,0,0,0,0,1,0,0,0}()
    const Candela   = SIUnit{0,0,0,0,0,0,1,0,0}()
    const Radian    = SIUnit{0,0,0,0,0,0,0,1,0}()
    const Steradian = SIUnit{0,0,0,0,0,0,0,0,1}()

    const Kilo       = (1000)SIPrefix
    const Mega       = (10^6)SIPrefix
    const Giga       = (10^9)SIPrefix
    const Tera       = (10^12)SIPrefix
    const Peta       = (10^15)SIPrefix
    const Exa        = (10^18)SIPrefix
    const Zetta      = (10^21)SIPrefix
    const Yotta      = (10^24)SIPrefix
    const Centi      = (1//100)SIPrefix
    const Milli      = (1//1000)SIPrefix
    const Micro      = (1//10^6)SIPrefix
    const Nano       = (1//10^9)SIPrefix
    const Pico       = (1//10^12)SIPrefix
    const Femto      = (1//10^15)SIPrefix
    const Atto       = (1//10^18)SIPrefix
    const Zepto      = (1//10^21)SIPrefix
    const Yocto      = (1//10^24)SIPrefix

    const Gram       = (1//1000)KiloGram
    const Joule      = KiloGram*Meter^2/Second^2
    const Coulomb    = Ampere*Second
    const Volt       = Joule/Coulomb
    const Farad      = Coulomb^2/Joule
    const Newton     = KiloGram*Meter/Second^2
    const Ohm        = Volt/Ampere
    const Hertz      = inv(Second)
    const Siemens    = inv(Ohm)
    const Watt       = Joule/Second
    const Pascal     = Newton/Meter^2

    const CentiMeter = Centi*Meter

# Pretty Printing - Text
    superscript(i) = map(repr(i)) do c
        c   ==  '-' ? '\u207b' :
        c   ==  '1' ? '\u00b9' :
        c   ==  '2' ? '\u00b2' :
        c   ==  '3' ? '\u00b3' :
        c   ==  '4' ? '\u2074' :
        c   ==  '5' ? '\u2075' :
        c   ==  '6' ? '\u2076' :
        c   ==  '7' ? '\u2077' :
        c   ==  '8' ? '\u2078' :
        c   ==  '9' ? '\u2079' :
        c   ==  '0' ? '\u2070' :
        error("Unexpected Chatacter")
    end

    function spacing(idx::Int, x::SIUnit)
        # Only print a space if there are nonzero units coming after this one
        tup(x)[idx+1:end] == ntuple((i)->0, 9-idx) ? "" : " "
    end
    function show{m,kg,s,A,K,mol,cd,rad,sr}(io::IO,x::SIUnit{m,kg,s,A,K,mol,cd,rad,sr})
        kg  != 0 && print(io, "kg",  (kg  == 1 ? spacing(1,x) : superscript(kg)))
        m   != 0 && print(io, "m",   (m   == 1 ? spacing(2,x) : superscript(m)))
        s   != 0 && print(io, "s",   (s   == 1 ? spacing(3,x) : superscript(s)))
        A   != 0 && print(io, "A",   (A   == 1 ? spacing(4,x) : superscript(A)))
        K   != 0 && print(io, "K",   (K   == 1 ? spacing(5,x) : superscript(K)))
        mol != 0 && print(io, "mol", (mol == 1 ? spacing(6,x) : superscript(mol)))
        cd  != 0 && print(io, "cd",  (cd  == 1 ? spacing(7,x) : superscript(cd)))
        rad != 0 && print(io, "rad", (rad == 1 ? spacing(8,x) : superscript(rad)))
        sr  != 0 && print(io, "sr",  (sr  == 1 ? ""           : superscript(sr)))
        nothing
    end

    function show{T,m,kg,s,A,K,mol,cd,rad,sr}(io::IO,x::SIQuantity{T,m,kg,s,A,K,mol,cd,rad,sr})
        show(io,x.val)
        print(io," ")
        show(io,unit(x))
    end

    function sidims{m,kg,s,A,K,mol,cd,rad,sr}(::SIUnit{m,kg,s,A,K,mol,cd,rad,sr})
        (m,kg,s,A,K,mol,cd,rad,sr)
    end

    function sidims{T,m,kg,s,A,K,mol,cd,rad,sr}(::SIQuantity{T,m,kg,s,A,K,mol,cd,rad,sr})
        (m,kg,s,A,K,mol,cd,rad,sr)
    end

    export @prettyshow

    macro prettyshow(unit,string)
        esc(quote function Base.show(io::IO,::SIUnits.SIUnit{SIUnits.sidims($(unit))...})
            print(io,$(string))
        end
        @compat function Base.show(io::IO,::MIME"text/mathtex+latex",::SIUnits.SIUnit{SIUnits.sidims($(unit))...})
            show(io,MIME("text/mathtex+latex"),$(string))
        end
        end)
    end

# Pretty Printing - LaTeX

    using TexExtensions

    @compat import Base: show

    macro l(x)
        esc(quote
            $x != 0 && push!($x>0?num:den,string("\\text{",$(string(x)),"\}",abs($x) == 1 ? " " : string("^{",abs($x),"}")))
        end)
    end

    @compat function show{m,kg,s,A,K,mol,cd,rad,sr}(io::IO,::MIME"text/mathtex+latex",x::SIUnit{m,kg,s,A,K,mol,cd,rad,sr})
        num = String[]
        den = String[]
        @l kg
        @l m
        @l s
        @l A
        @l K
        @l mol
        @l cd
        @l rad
        @l sr
        if !isempty(den)
            if isempty(num)
                write(io,"\\frac{1}{",join(den,"\\;"),"}")
            else
                write(io,"\\frac{",join(num,"\\;"),"}{",join(den,"\\;"),"}")
            end
        else
            write(io,join(num,"\\;"))
        end
    end

    @compat function show{T,m,kg,s,A,K,mol,cd,rad,sr}(io::IO,::MIME"text/mathtex+latex",x::SIQuantity{T,m,kg,s,A,K,mol,cd,rad,sr})
        show(io,MIME("text/mathtex+latex"),x.val)
        write(io,"\\;")
        show(io,MIME("text/mathtex+latex"),unit(x))
    end

# Non-SI Units
immutable NonSIUnit{BaseUnit<:SIUnit,Unit}; end
immutable NonSIQuantity{T<:Number,Unit<:NonSIUnit} <: Number
    val::T
end

# Non-SI promote rules
promote_rule(x::Type{Irrational},y::Type{NonSIUnit}) =
    NonSIQuantity{x,y}
promote_rule{sym,T<:Number,Unit}(x::Type{Irrational{sym}},y::Type{NonSIQuantity{T,Unit}}) =
    NonSIQuantity{promote_type(Irrational{sym},T),Unit}

promote_rule{T<:Number,S<:Number,U1,U2}(
    A::Type{NonSIQuantity{T,U1}},B::Type{SIQuantity{S,U2}}) = NonSIQuantity{promote_type(T,S)}
promote_rule{T<:Number,U1}(
    A::Type{NonSIQuantity{T,U1}},U2::Type{NonSIUnit}) = NonSIQuantity{T}
promote_rule{S<:Number,U}(x::Type{Bool},y::Type{NonSIQuantity{S,U}}) = NonSIQuantity{promote_type(Bool,S),U}
promote_rule(x::Type{Bool},U::Type{NonSIUnit}) = NonSIQuantity{Bool,U}
promote_rule{T<:Number,S<:Number,U}(x::Type{T},y::Type{NonSIQuantity{S,U}}) = NonSIQuantity{promote_type(T,S),U}
promote_rule{T<:Number}(x::Type{T},U::Type{NonSIUnit}) = NonSIQuantity{T,U}

# Interaction between SI and non-SI quantities
promote_rule{S<:Number,T<:Number,U,m,kg,s,A,K,mol,cd,rad,sr}(x::Type{NonSIQuantity{S,U}},y::Type{SIQuantity{T,m,kg,s,A,K,mol,cd,rad,sr}}) =
    SIQuantity{promote_type(S,T)}
promote_rule{S<:Number,T<:Number,U,m,kg,s,A,K,mol,cd,rad,sr}(x::Type{SIQuantity{T,m,kg,s,A,K,mol,cd,rad,sr}},y::Type{NonSIQuantity{S,U}}) =
    SIQuantity{promote_type(S,T)}

siquantity{B}(T,U::NonSIUnit{B}) = quantity(T,B())
siquantity{B}(T,U::Type{NonSIUnit{B}}) = quantity(T,B())
#convert{T,S,U}(::Type{SIQuantity{T}},x::NonSIQuantity{S,U}) = (siquantity(promote_type(T,S),U())(x.val))

*{T<:NonSIUnit}(x::Number,t::T) = NonSIQuantity{typeof(x),T}(x)

baseunit{BaseUnit}(x::NonSIUnit{BaseUnit}) = BaseUnit()
baseunit{T,Unit}(x::NonSIQuantity{T,Unit}) = baseunit(unit(x))
unit{T,Unit}(x::NonSIQuantity{T,Unit}) = Unit()
quantity(T::(@compat Union{Type,TypeVar}),x::NonSIUnit) = NonSIQuantity{T,typeof(x)}
quantity(T::(@compat Union{Type,TypeVar}),U::Type{NonSIUnit}) = NonSIQuantity{T,U}

/{T,U}(x::NonSIQuantity{T,U},y::SIQuantity) = convert(SIQuantity,x)/y
/(x::NonSIUnit,y::SIUnit) = convert(SIQuantity,x)/y
/(x::SIUnit,y::NonSIUnit) = x/convert(SIQuantity,y)

/(x::SIQuantity,y::NonSIUnit) = x/convert(SIQuantity,y)
/(x::NonSIUnit,y::SIQuantity) = convert(SIQuantity,x)/y
-{T,U}(x::NonSIQuantity{T,U}) = NonSIQuantity{T,U}(-x.val)

^(x::(@compat Union{NonSIQuantity,NonSIUnit}),i::Integer) = convert(SIQuantity,x)^i

show{BaseUnit,Unit}(io::IO,x::NonSIUnit{BaseUnit,Unit}) = write(io,string(Unit))
function show(io::IO,x::NonSIQuantity)
    show(io,x.val)
    print(io," ")
    show(io,unit(x))
end

@compat function show{BaseUnit,Unit}(io::IO,::MIME"text/mathtex+latex",x::NonSIUnit{BaseUnit,Unit})
    write(io,"\\text{",string(Unit),"}")
end

@compat function show(io::IO,::MIME"text/mathtex+latex",x::NonSIQuantity)
    show(io,MIME("text/mathtex+latex"),x.val)
    write(io,"\\;")
    show(io,MIME("text/mathtex+latex"),unit(x))
end

# No, these two are not the same. Sometimes we get SIQuantities that are not specialized
# on the type out of promotion, hence the first method, but we still need the second method
# to be more specific that the convert methods of plain SIUnits above.
convert(::Type{SIQuantity},x::NonSIQuantity) = x.val * convert(SIQuantity,unit(x))
convert{T<:Number}(::Type{SIQuantity{T}},x::NonSIQuantity) = x.val * convert(SIQuantity,unit(x))

convert{T}(::Type{NonSIQuantity{T}},U::NonSIUnit) = NonSIQuantity{T,U}(one(T))
convert{T<:Number,U}(::Type{NonSIQuantity{T,U}},x::T) = UnitQuantity{T}(x)
#convert{T}(::Type{NonSIQuantity{T}},x::T) = UnitQuantity{T}(x)
#convert{T,S}(::Type{NonSIQuantity{T}},x::S) = convert(NonSIQuantity{T},convert(T,x))
if VERSION >= v"0.4-dev"
    eval(quote
        convert{T,U}(::Type{NonSIQuantity{T,U}},x::Dates.Period) = error("Conversion from Period to NonSIQuantity not defined")
    end)
end
convert{T,U}(::Type{NonSIQuantity{T,U}},x::Number) = convert(NonSIQuantity{T,U},convert(T,x))
convert{T,U}(::Type{NonSIQuantity{T,U}},x::NonSIQuantity{T,U}) = x
convert{T<:Number,S,U1,U2}(::Type{NonSIQuantity{T,U1}},x::NonSIQuantity{S,U2}) = NonSIQuantity{T,U2}(convert(T,x.val))

export as

function as{U<:NonSIUnit}(x::SIQuantity,y::U)
    val = x/y
    @assert !(typeof(val)<:SIQuantity)
    NonSIQuantity{typeof(val),U}(val)
end

function as{U<:NonSIUnit,Q<:SIQuantity}(X::AbstractArray{Q},y::U)
    val = [x/y for x in X]
    @assert !(typeof(eltype(val))<:SIQuantity)
    NonSIQuantity{typeof(val),U}(val)
end

function as(x::NonSIQuantity,y::SIUnit)
    @assert baseunit(x) == y
    convert(SIQuantity,x)
end

as(x::NonSIQuantity,y::SIQuantity) = as(x,unit(y))

include("nonsiunits.jl")
include("shortunits.jl")

end # module


text:
# # Instrument Line Shapes

# Using packages:
using Plots
using Plots.PlotMeasures
# This needs to be installed from https://github.com/RadiativeTransfer/RadiativeTransfer.jl
using RadiativeTransfer.Absorption
using InstrumentOperator

# ## Load HITRAN data and CO2 cross sections
hitran_data = read_hitran(artifact("CO2"), mol=2, iso=1, ν_min=6000, ν_max=6400)
line_co2 = make_hitran_model(hitran_data, Voigt(), architecture=CPU())

# ## Compute a high resolution transmission spectrum
Δν = 0.0025
ν_min = 6320;
ν_max = 6355;
ν = ν_min:Δν:ν_max;
# CO₂ cross section at 800hPa and 296K:
σ_co2   = absorption_cross_section(line_co2, ν, 800.0, 296.0);
# Transmission through a typical atmosphere (8e21 molec/cm²)
T = exp.(-8e21*σ_co2);

# plot high resolution transmission
plot(ν, T, lw=2, label="High resolution",  top_margin = 10mm, bottom_margin = 10mm, left_margin = 5mm, right_margin = 5mm)

# Define the instrument kernel grid
x = -8:Δν:8;
# ## Create a kernel at a center wavenumber of 6300cm⁻¹
# Use FTS with 5cm MOPD, FOV of 7.9mrad, no assymetry:
FTS = FTSInstrument(5.0, 7.9e-3, 0.0)
FTSkernel = create_instrument_kernel(FTS, x,6300.0)
margin = 5.0
sampling = 0.01
FTS_instr = FixedKernelInstrument(FTSkernel, collect(ν_min+margin:sampling:ν_max-margin))

# ## Convolve with instrument kernel
T_conv = conv_spectra(FTS_instr, ν, T);

# overplot convolved transmission:
plot!(FTS_instr.ν_out, T_conv, label="Convolved")

ext: jl


text:
module TransitionParser

include("beamsearch.jl")

end


text:
import CUDA.CUDNN: cudnnNormalizationForwardAD
import AutoGrad: forw
using CUDA.CUDNN: handle, CU_NULL, cudnnNormalizationBackward, cudnnNormalizationForwardTraining, cudnnNormalizationForwardInference, cudnnGetNormalizationForwardTrainingWorkspaceSize, cudnnGetNormalizationTrainingReserveSpaceSize, cudnnTempSpace, scalingParameter
using AutoGrad: AutoGrad, @primitive1, forwargs, recording, Result, value

@primitive1((cudnnNormalizationForwardAD(x, scale, bias, z; training, mean, variance, y, mode, normOps, algo, alpha, beta, epsilon, groupCnt, exponentialAverageFactor, savedMean, savedInvVariance, activationDesc, xDesc, yDesc, zDesc, normScaleBiasDesc, normMeanVarDesc, workspace, reserveSpace, dx, dscale, dbias, dz, dready), dy, _y),
            (dready[] || cudnnNormalizationBack(dy, x, scale, bias, z; training, mean, variance, y, mode, normOps, algo, alpha, beta, epsilon, groupCnt, exponentialAverageFactor, savedMean, savedInvVariance, activationDesc, xDesc, yDesc, zDesc, normScaleBiasDesc, normMeanVarDesc, workspace, reserveSpace, dx, dscale, dbias, dz, dready); dx[]),
            (dready[] || cudnnNormalizationBack(dy, x, scale, bias, z; training, mean, variance, y, mode, normOps, algo, alpha, beta, epsilon, groupCnt, exponentialAverageFactor, savedMean, savedInvVariance, activationDesc, xDesc, yDesc, zDesc, normScaleBiasDesc, normMeanVarDesc, workspace, reserveSpace, dx, dscale, dbias, dz, dready); dscale[]),
            (dready[] || cudnnNormalizationBack(dy, x, scale, bias, z; training, mean, variance, y, mode, normOps, algo, alpha, beta, epsilon, groupCnt, exponentialAverageFactor, savedMean, savedInvVariance, activationDesc, xDesc, yDesc, zDesc, normScaleBiasDesc, normMeanVarDesc, workspace, reserveSpace, dx, dscale, dbias, dz, dready); dbias[]),
            (dready[] || cudnnNormalizationBack(dy, x, scale, bias, z; training, mean, variance, y, mode, normOps, algo, alpha, beta, epsilon, groupCnt, exponentialAverageFactor, savedMean, savedInvVariance, activationDesc, xDesc, yDesc, zDesc, normScaleBiasDesc, normMeanVarDesc, workspace, reserveSpace, dx, dscale, dbias, dz, dready); dz[]))

function cudnnNormalizationBack(dy, x, scale, bias, z; training, mean, variance, y, mode, normOps, algo, alpha, beta, epsilon, groupCnt, exponentialAverageFactor, savedMean, savedInvVariance, activationDesc, xDesc, yDesc, zDesc, normScaleBiasDesc, normMeanVarDesc, workspace, reserveSpace, dx, dscale, dbias, dz, dready)
    # Make sure backward gets called only once
    dready[] && return
    dready[] = true
    (x, scale, bias, z) = value.((x, scale, bias, z))
    # Allocate gradient buffers if necessary
    !isassigned(dx) || dx[]===nothing ? dx[] = similar(x) : @assert issimilar(x, dx[])
    !isassigned(dscale) || dscale[]===nothing ? dscale[] = similar(scale) : @assert issimilar(scale, dscale[])
    !isassigned(dbias) || dbias[]===nothing ? dbias[]  = similar(bias) : @assert issimilar(bias, dbias[])
    z === nothing ? dz[] = nothing : !isassigned(dz) || dz[]===nothing ? dz[] = zero(z) : (@assert issimilar(z, dz[]); dz[] .= 0) # z may not be used, dz may not be modified, should be zeroed
    s0,s1 = scalingParameter(eltype(x),0),scalingParameter(eltype(x),1)
    cudnnNormalizationBackward(handle(), mode, normOps, algo, alpha, s0, alpha, s0, xDesc, x, yDesc, y, yDesc, dy, something(zDesc,C_NULL), something(dz[],CU_NULL), xDesc, dx[], normScaleBiasDesc, scale, bias, dscale[], dbias[], epsilon, something(normMeanVarDesc,C_NULL), something(savedMean,CU_NULL), something(savedInvVariance,CU_NULL), something(activationDesc,C_NULL), something(workspace,CU_NULL), sizeof(workspace), something(reserveSpace,CU_NULL), sizeof(reserveSpace), groupCnt)
end

# CUDNN API controls training vs inference mode with a simple keyword argument.
# In Knet things are a bit more complicated: whether or not we are in @diff mode and we have Value args enter the picture.
# With Value args the AutoGrad.forw method is called regardless of the kwarg / @diff mode.
# The training keyword argument defaults to Knet.training() in Knet API, so by default @diff=training, but the user can override this.
# We have 8 cases to consider for Value, training, @diff:
# When there are no Value arguments, i.e. scale/bias not Params and x is not a Param or Result:
# cudnnNormalizationForwardAD is called, @diff=false, training=false: AD calls Inference: ok
# cudnnNormalizationForwardAD is called, @diff=false, training=true: AD calls Training: user is forcing training mode
# cudnnNormalizationForwardAD is called, @diff=true, training=false: AD calls Inference: user is forcing inference mode
# cudnnNormalizationForwardAD is called, @diff=true, training=true: AD calls Training: ok (but no Value args no back pass)
# When there are some Value arguments:
# forw is called, @diff=false, training=false: Standard test mode. forw->AD->Inference: ok
# forw is called, @diff=false, training=true: forw->AD->Training: We are not taking gradients but user is forcing training mode.
# forw is called, @diff=true, training=false: forw->AD->Inference: User wants to force inference mode?
# forw is called, @diff=true, training=true: Standard training mode. forw->AD->Training: ok


text:
author_login(pull_request::GitHub.PullRequest) = pull_request.user.login

base_repo(pull_request::GitHub.PullRequest) = pull_request.base.repo

body(c::GitHub.Comment) = c.body

function created_at(pull_request::GitHub.PullRequest)
    result = time_is_already_in_utc(pull_request.created_at)
    return result
end

function delete_comment!(api::GitHub.GitHubAPI,
                         repo::GitHub.Repo,
                         pr::GitHub.PullRequest,
                         comment_to_delete::GitHub.Comment;
                         auth::GitHub.Authorization)
    GitHub.delete_comment(api, repo,
                          comment_to_delete,
                          :pr;
                          auth = auth)
    return nothing
end

function delete_merged_branch!(api::GitHub.GitHubAPI, repo::GitHub.Repo, pr::GitHub.PullRequest; auth::GitHub.Authorization)
    updated_pr = _get_updated_pull_request(api, pr; auth=auth)
    if is_merged(updated_pr)
        try
            head_branch = pull_request_head_branch(updated_pr)
            repo = head_branch.repo
            ref = "heads/$(head_branch.ref)"
            GitHub.delete_reference(api, repo, ref; auth=auth)
        catch ex
            # showerror(stderr, ex)
            # Base.show_backtrace(stderr, catch_backtrace())
            # println(stderr)
        end
    end
    return nothing
end

function edit_comment!(api::GitHub.GitHubAPI,
                       repo::GitHub.Repo,
                       pr::GitHub.PullRequest,
                       comment::GitHub.Comment,
                       body::String;
                       auth::GitHub.Authorization)
    myparams = Dict("body" => body)
    GitHub.edit_comment(api, repo, comment, :pr; auth=auth, params = myparams)
    return nothing
end

full_name(repo::GitHub.Repo) = repo.full_name

function _get_updated_pull_request(api::GitHub.GitHubAPI, pull_request::GitHub.PullRequest; auth::GitHub.Authorization)
    pr_base_repo = base_repo(pull_request)
    pr_number = number(pull_request)
    updated_pr = GitHub.pull_request(api, pr_base_repo, pr_number; auth=auth)
    return updated_pr
end

function get_all_my_pull_request_comments(api::GitHub.GitHubAPI,
                                          repo::GitHub.Repo,
                                          pr::GitHub.PullRequest;
                                          auth::GitHub.Authorization,
                                          whoami)
    all_comments = get_all_pull_request_comments(api, repo,
                                                 pr;
                                                 auth = auth)
    my_comments = Vector{GitHub.Comment}(undef, 0)
    for c in all_comments
        if c.user.login == whoami
            push!(my_comments, c)
        end
    end
    unique!(my_comments)
    my_comments = my_comments[sortperm([x.created_at for x in my_comments])]
    return my_comments
end

function get_all_pull_request_comments(api::GitHub.GitHubAPI,
                                       repo::GitHub.Repo,
                                       pr::GitHub.PullRequest;
                                       auth::GitHub.Authorization)
    all_comments = Vector{GitHub.Comment}(undef, 0)
    myparams = Dict("per_page" => 100, "page" => 1)
    cs, page_data = GitHub.comments(api, repo, pr, :pr; auth=auth, params = myparams, page_limit = 100)
    append!(all_comments, cs)
    while haskey(page_data, "next")
        cs, page_data =  GitHub.comments(api, repo, pr, :pr; auth=auth, page_limit = 100, start_page = page_data["next"])
        append!(all_comments, cs)
    end
    unique!(all_comments)
    all_comments = all_comments[sortperm([x.created_at for x in all_comments])]
    return all_comments
end

function get_all_pull_requests(api::GitHub.GitHubAPI,
                               repo::GitHub.Repo,
                               state::String;
                               auth::GitHub.Authorization)
    all_pull_requests = Vector{GitHub.PullRequest}(undef, 0)
    myparams = Dict("state" => state, "per_page" => 100, "page" => 1)
    prs, page_data = GitHub.pull_requests(api, repo; auth=auth, params = myparams, page_limit = 100)
    append!(all_pull_requests, prs)
    while haskey(page_data, "next")
        prs, page_data = GitHub.pull_requests(api, repo; auth=auth, page_limit = 100, start_page = page_data["next"])
        append!(all_pull_requests, prs)
    end
    unique!(all_pull_requests)
    return all_pull_requests
end

function get_changed_filenames(api::GitHub.GitHubAPI, repo::GitHub.Repo, pull_request::GitHub.PullRequest; auth::GitHub.Authorization)
    files = GitHub.pull_request_files(api, repo, pull_request; auth=auth)
    n = length(files)
    filenames = Vector{String}(undef, n)
    for i = 1:n
        filenames[i] = files[i].filename
    end
    return filenames
end

is_merged(pull_request::GitHub.PullRequest) = pull_request.merged

function is_open(pull_request::GitHub.PullRequest)
    result = pr_state(pull_request) == "open"
    !result && @error("Pull request is not open")
    return result
end

function merge!(api::GitHub.GitHubAPI,
                registry_repo::GitHub.Repo,
                pr::GitHub.PullRequest,
                approved_pr_head_sha::AbstractString;
                auth::GitHub.Authorization)
    pr = wait_pr_compute_mergeability(api, registry_repo, pr; auth = auth)
    _approved_pr_head_sha = convert(String, strip(approved_pr_head_sha))::String
    pr_number = number(pr)
    @info("Attempting to squash-merge pull request #$(pr_number)")
    @debug("sha = $(_approved_pr_head_sha)")
    @debug("pr.mergeable = $(pr.mergeable)")
    params = Dict("sha" => _approved_pr_head_sha, "merge_method" => "squash")
    try
        GitHub.merge_pull_request(api, registry_repo,
                                  pr_number;
                                  auth=auth,
                                  params=params)
    catch ex
        showerror(stderr, ex)
        Base.show_backtrace(stderr, catch_backtrace())
        println(stderr)
    end
    try
        delete_merged_branch!(api, registry_repo, pr; auth=auth)
    catch
    end
    return nothing
end

function wait_pr_compute_mergeability(api::GitHub.GitHubAPI,
                                      repo::GitHub.Repo,
                                      pr::GitHub.PullRequest;
                                      auth::GitHub.Authorization)
    sleep(5)
    max_tries = 10
    num_tries = 0
    pr = GitHub.pull_request(api, repo, pr.number; auth = auth)
    while !(pr.mergeable isa Bool) && num_tries <= max_tries
        num_tries = num_tries + 1
        sleep(5)
        pr = GitHub.pull_request(api, repo, pr.number; auth = auth)
    end
    return pr
end

num_changed_files(pull_request::GitHub.PullRequest) = pull_request.changed_files

number(pull_request::GitHub.PullRequest) = pull_request.number

function post_comment!(api::GitHub.GitHubAPI,
                       repo::GitHub.Repo,
                       pr::GitHub.PullRequest,
                       body::String;
                       auth::GitHub.Authorization)
    myparams = Dict("body" => body)
    GitHub.create_comment(api, repo, pr, :pr; auth=auth, params = myparams)
    return nothing
end

pull_request_head_branch(pull_request::GitHub.PullRequest) = pull_request.head

pull_request_head_sha(pull_request::GitHub.PullRequest) = pull_request.head.sha

repo_url(repo::GitHub.Repo) = repo.html_url.uri

pr_state(pull_request::GitHub.PullRequest) = pull_request.state

function time_since_pr_creation(pull_request::GitHub.PullRequest)
    _pr_created_at_utc = created_at(pull_request)
    _now_utc = now_utc()
    result = _now_utc - _pr_created_at_utc
    return result
end

title(pull_request::GitHub.PullRequest) = pull_request.title

function username(api::GitHub.GitHubAPI, auth::GitHub.Authorization)
    user_information = GitHub.gh_get_json(api,
                                          "/user";
                                          auth = auth)
    return user_information["login"]::String
end


text:
"""
    reverse(A; dims) => AbstractDimArray
    reverse(dim::Dimension) => Dimension

Reverse the array order, and update the dim to match.
"""
Base.reverse(A::AbstractDimArray; dims=1) =
    _reverse(IndexOrder, _reverse(ArrayOrder, A, dims), dims)
Base.reverse(s::AbstractDimStack; dims=1) = map(A -> reverse(A; dims=dims), s)
Base.reverse(ot::Type{<:SubOrder}, x; dims) = _reverse(ot, x, dims)
_reverse(ot::Type{<:SubOrder}, x, lookup) = set(x, reverse(ot, dims(x, lookup)))
_reverse(ot::Type{<:ArrayOrder}, x, lookup) = begin
    newdims = reverse(ArrayOrder, dims(x, lookup))
    newdata = _reversedata(x, dimnum(x, lookup))
    setdims(rebuild(x, newdata), newdims)
end
# Dimension
Base.reverse(ot::Type{<:SubOrder}, dims::DimTuple) = map(d -> reverse(ot, d), dims)
Base.reverse(ot::Type{<:SubOrder}, dim::Dimension) = reverse(ot, mode(dim), dim)
Base.reverse(ot::Type{<:SubOrder}, mode::IndexMode, dim::Dimension) =
    rebuild(dim; mode=reverse(ot, mode))
# Reverse the index
Base.reverse(ot::Type{<:IndexOrder}, mode::IndexMode, dim::Dimension) =
    rebuild(dim, reverse(index(dim)), reverse(ot, mode))
Base.reverse(ot::Type{<:IndexOrder}, mode::IndexMode, dim::Dimension{<:Val{I}}) where I =
    rebuild(dim, Val(reverse(I)), reverse(ot, mode))
Base.reverse(ot::Type{<:IndexOrder}, mode::NoIndex, dim::Dimension) = dim
Base.reverse(dim::Dimension) = reverse(IndexOrder, dim)
# Mode
Base.reverse(ot::Type{<:SubOrder}, mode::NoIndex) = mode
Base.reverse(ot::Type{<:SubOrder}, mode::IndexMode) =
    rebuild(mode; order=reverse(ot, order(mode)))
Base.reverse(ot::Type{<:SubOrder}, mode::AbstractSampled) =
    rebuild(mode; order=reverse(ot, order(mode)), span=reverse(ot, span(mode)))
# Order
Base.reverse(::Type{<:IndexOrder}, o::Ordered) =
    Ordered(reverse(indexorder(o)), arrayorder(o), reverse(relation(o)))
Base.reverse(::Type{<:Union{ArrayOrder,Relation}}, o::Ordered) =
    Ordered(indexorder(o), reverse(arrayorder(o)), reverse(relation(o)))
Base.reverse(::Type{<:SubOrder}, o::Unordered) = Unordered(reverse(relation(o)))
# SubOrder
Base.reverse(::ReverseIndex) = ForwardIndex()
Base.reverse(::ForwardIndex) = ReverseIndex()
Base.reverse(::ReverseArray) = ForwardArray()
Base.reverse(::ForwardArray) = ReverseArray()
Base.reverse(::ReverseRelation) = ForwardRelation()
Base.reverse(::ForwardRelation) = ReverseRelation()
# Span
Base.reverse(::Type{<:SubOrder}, span::Span) = span
Base.reverse(::Type{<:IndexOrder}, span::Irregular) = span
Base.reverse(::Type{<:IndexOrder}, span::Span) = reverse(span)
Base.reverse(span::Regular) = Regular(-step(span))
Base.reverse(span::Explicit) = Explicit(reverse(val(span), dims=2))

_reversedata(A::AbstractDimArray, dimnum) = reverse(parent(A); dims=dimnum)
_reversedata(s::AbstractDimStack, dimnum) =
    map(a -> reverse(parent(a); dims=dimnum), data(s))

"""
    fliparray(Order, A, dims) => AbstractDimArray
    fliparray(dim::Dimension) => Dimension

`Flip` the array order without changing any data.
"""
function flip end

flip(ot::Type{<:SubOrder}, x; dims) = flip(ot, x, dims)
flip(ot::Type{<:SubOrder}, x, lookupdims) = set(x, flip(ot, dims(x, lookupdims)))
flip(ot::Type{<:SubOrder}, dims::DimTuple) = map(d -> flip(ot, d), dims)
flip(ot::Type{<:SubOrder}, dim::Dimension) = _set(dim, flip(ot, mode(dim)))
flip(ot::Type{<:SubOrder}, mode::IndexMode) = _set(mode, flip(ot, order(mode)))
flip(ot::Type{<:SubOrder}, mode::NoIndex) = mode
flip(ot::Type{<:SubOrder}, o::Order) = _set(o, reverse(ot, o))

"""
    reorder(::Order, A::AbstractDimArray) => AbstractDimArray
    reorder(::Order, A::Dimension) => AbstractDimArray
    reorder(A::AbstractDimArray, order::Union{Order,Dimension{<:Order},Tuple}) => AbstractDimArray
    reorder(A::AbstractDimArray, order::Pair{<:Dimension,<:SubOrder}...) => AbstractDimArray

Reorder every dims index/array/relation to `order`, or reorder index for
the the given dimension(s) to the `Order` they wrap.

Reordering `Relation` will reverse the array, not the dimension index.

`order` can be an [`Order`](@ref), a single [`Dimension`](@ref)
or a `Tuple` of `Dimension`.

If no axis reversal is required the same objects will be returned, without allocation.
"""
function reorder end

reorder(x, args...; kw...) = reorder(x, (args..., _kwdims(kw)...))
reorder(x, nt::NamedTuple) = reorder(x, _kwdims(nt))
reorder(x, p::Pair, ps::Vararg{<:Pair}) = reorder(x, (p, ps...))
reorder(x, ps::Tuple{Vararg{<:Pair}}) = reorder(x, _pairdims(ps...))
# Reorder specific dims.
reorder(x, dimwrappers::Tuple) = _reorder(x, dimwrappers)
# Reorder all dims.
reorder(x, ot::Union{SubOrder,Type{<:SubOrder}}) =
    _reorder(x, map(d -> basetypeof(d)(ot), dims(x)))
reorder(dim::Dimension, ot::Union{SubOrder,Type{<:SubOrder}}) =
    ot == basetypeof(order(ot, dim)) ? dim : reverse(ot, dim)

# Recursive reordering. x may be reversed here
_reorder(x, orderdims::DimTuple) = _reorder(reorder(x, orderdims[1]), tail(orderdims))
_reorder(x, orderdims::Tuple{}) = x

reorder(x, orderdim::Dimension) = _reorder(val(orderdim), x, dims(x, orderdim))

_reorder(neworder::SubOrder, x, orderdim::DimOrDimType) = _reorder(basetypeof(neworder), x, orderdim)
# Reverse the dimension index
_reorder(ot::Type{<:IndexOrder}, x, orderdim::DimOrDimType) =
    ot == basetypeof(order(ot, orderdim)) ? x : set(x, reverse(ot, orderdim))
# If either ArrayOrder or Relation are reversed, we reverse the array
_reorder(ot::Type{<:Union{ArrayOrder,Relation}}, x, dim::DimOrDimType) =
    ot == basetypeof(order(ot, dim)) ? x : reverse(ArrayOrder, x; dims=dim)

"""
    modify(f, A::AbstractDimArray) => AbstractDimArray
    modify(f, s::AbstractDimStack) => AbstractDimStack
    modify(f, dim::Dimension) => Dimension
    modify(f, x, lookupdim::Dimension) => typeof(x)

Modify the parent data, rebuilding the object wrapper without
change. `f` must return a `AbstractArray` of the same size as the original.

This method is mostly useful as a way of swapping the parent array type of
an object.

## Example

If we have a previously-defined `DimArray`, we can copy it to an Nvidia GPU with:

```julia
A = DimArray(rand(100, 100), (X, Y))
modify(CuArray, A)
```

This also works for all the data layers in a `DimStack`.
"""
function modify end
modify(f, s::AbstractDimStack) = map(a -> modify(f, a), s)
function modify(f, A::AbstractDimArray)
    newdata = f(parent(A))
    size(newdata) == size(A) || error("$f returns an array with a different size")
    rebuild(A, newdata)
end
modify(f, x, dim::DimOrDimType) = set(x, modify(f, dims(x, dim)))
function modify(f, dim::Dimension)
    newindex = f(index(dim))
    size(newindex) == size(dim) || error("$f returns a vector with a different size")
    rebuild(dim, newindex)
end
function modify(f, dim::Dimension{<:Val{Index}}) where Index
    newindex = f(Index)
    length(newindex) == length(dim) || error("$f returns a Tuple with a different size")
    rebuild(dim, Val(newindex))
end

"""
    dimwise(f, A::AbstractDimArray{T,N}, B::AbstractDimArray{T2,M}) => AbstractDimArray{T3,N}

Dimension-wise application of function `f` to `A` and `B`.

## Arguments

- `a`: `AbstractDimArray` to broacast from, along dimensions not in `b`.
- `b`: `AbstractDimArray` to broadcast from all dimensions. Dimensions must be a subset of a.

This is like broadcasting over every slice of `A` if it is
sliced by the dimensions of `B`.
"""
function dimwise(f, A::AbstractDimArray, B::AbstractDimArray)
    dimwise!(f, similar(A, promote_type(eltype(A), eltype(B))), A, B)
end

"""
    dimwise!(f, dest::AbstractDimArray{T1,N}, A::AbstractDimArray{T2,N}, B::AbstractDimArray) => dest

Dimension-wise application of function `f`.

## Arguments

- `dest`: `AbstractDimArray` to update
- `A`: `AbstractDimArray` to broacast from, along dimensions not in `b`.
- `B`: `AbstractDimArray` to broadcast from all dimensions. Dimensions must be a subset of a.

This is like broadcasting over every slice of `A` if it is
sliced by the dimensions of `B`, and storing the value in `dest`.
"""
function dimwise!(
    f, dest::AbstractDimArray{T,N}, a::AbstractDimArray{TA,N}, b::AbstractDimArray{TB,NB}
) where {T,TA,TB,N,NB}
    N >= NB || error("B-array cannot have more dimensions than A array")
    comparedims(dest, a)
    common = commondims(a, dims(b))
    od = otherdims(a, common)
    # Lazily permute B dims to match the order in A, if required
    if !dimsmatch(common, dims(b))
        b = PermutedDimsArray(b, common)
    end
    # Broadcast over b for each combination of dimensional indices D
    if length(od) == 0
        dest .= f.(a, b)
    else
        map(DimIndices(od)) do D
            dest[D...] .= f.(a[D...], b)
        end
    end
    return dest
end

"""
    basetypeof(x) => Type

Get the "base" type of an object - the minimum required to
define the object without it's fields. By default this is the full
`UnionAll` for the type. But custom `basetypeof` methods can be
defined for types with free type parameters.

In DimensionalData this is primariliy used for comparing `Dimension`s,
where `Dim{:x}` is different from `Dim{:y}`.
"""
@inline basetypeof(x::T) where T = basetypeof(T)
@generated function basetypeof(::Type{T}) where T
    if T isa Union
        T
    else
        getfield(parentmodule(T), nameof(T))
    end
end

# Left pipe operator for cleaning up brackets
@deprecate f <| x f(x)

# Unwrap Val
unwrap(::Val{X}) where X = X
unwrap(::Type{Val{X}}) where X = X
unwrap(x) = x

# Get a tuple of unique keys for DimArrays. If they have the same
# name we call them layerI.
function uniquekeys(das::Tuple{AbstractDimArray,Vararg{<:AbstractDimArray}})
    uniquekeys(Symbol.(map(name, das)))
end
function uniquekeys(keys::Tuple{Symbol,Vararg{<:Symbol}})
    ids = ntuple(x -> x, length(keys))
    map(keys, ids) do k, id
        count(k1 -> k == k1, keys) > 1 ? Symbol(:layer, id) : k
    end
end

# Produce a 2 * length(dim) matrix of interval bounds from a dim
function dim2boundsmatrix(dim::Dimension)
    samp = sampling(dim)
    samp isa Intervals || error("Cannot create a bounds matrix for $(nameof(typeof(samp)))")
    _dim2boundsmatrix(locus(dim), span(dim), dim)
end

_dim2boundsmatrix(::Locus, span::Explicit, dim) = val(span)
_dim2boundsmatrix(::Locus, span::Regular, dim) =
    vcat(permutedims(_shiftindexlocus(Start(), dim)), permutedims(_shiftindexlocus(End(), dim)))
@noinline _dim2boundsmatrix(::Center, span::Regular{Dates.TimeType}, dim) =
    error("Cannot convert a Center TimeType index to Explicit automatically: use a bounds matrix e.g. Explicit(bnds)")
@noinline _dim2boundsmatrix(::Start, span::Irregular, dim) =
    error("Cannot convert Irregular to Explicit automatically: use a bounds matrix e.g. Explicit(bnds)")

#=
Shift the index from the current locus to the new locus. We only actually
shift Regular Intervals, and do this my multiplying the offset of
-1, -0.5, 0, 0.5 or 1 by the absolute value of the span.
=#
function shiftlocus(locus::Locus, dim::Dimension)
    samp = sampling(dim)
    samp isa Intervals || error("Cannot shift locus of $(nameof(typeof(samp)))")
    rebuild(dim; val=_shiftindexlocus(locus, dim), mode=DD._set(mode(dim), locus))
end

_shiftindexlocus(locus::Locus, dim::Dimension) = _shiftindexlocus(locus::Locus, mode(dim), dim)
_shiftindexlocus(locus::Locus, mode::IndexMode, dim::Dimension) = index(dim)
_shiftindexlocus(locus::Locus, mode::AbstractSampled, dim::Dimension) =
    _shiftindexlocus(locus, span(mode), sampling(mode), dim)
_shiftindexlocus(locus::Locus, span::Span, sampling::Sampling, dim::Dimension) = index(dim)
_shiftindexlocus(destlocus::Locus, span::Regular, sampling::Intervals, dim::Dimension) =
    index(dim) .+ (abs(step(span)) * _offset(locus(sampling), destlocus))
_shiftindexlocus(::Start, span::Explicit, sampling::Intervals, dim::Dimension) = val(span)[1, :]
_shiftindexlocus(::End, span::Explicit, sampling::Intervals, dim::Dimension) = val(span)[2, :]
_shiftindexlocus(destlocus::Center, span::Explicit, sampling::Intervals, dim::Dimension) =
    _shiftindexlocus(destlocus, locus(dim), span, sampling, dim)
_shiftindexlocus(::Center, ::Center, span::Explicit, sampling::Intervals, dim::Dimension) = index(dim)
_shiftindexlocus(::Center, ::Locus, span::Explicit, sampling::Intervals, dim::Dimension) =
    view(val(span), 2, :)  .- view(val(span), 1, :)

_offset(::Start, ::Center) = 0.5
_offset(::Start, ::End) = 1
_offset(::Center, ::Start) = -0.5
_offset(::Center, ::End) = 0.5
_offset(::End, ::Start) = -1
_offset(::End, ::Center) = -0.5
_offset(::T, ::T) where T<:Locus = 0

maybeshiftlocus(locus::Locus, dim::Dimension) = _maybeshiftlocus(locus, sampling(dim), dim)

_maybeshiftlocus(locus::Locus, sampling::Intervals, dim::Dimension) = shiftlocus(locus, dim)
_maybeshiftlocus(locus::Locus, sampling::Sampling, dim::Dimension) = dim

_astuple(t::Tuple) = t
_astuple(x) = (x,)


text:
"""
Gene struct.

# Fields
````
id :: String
notes :: Dict{String, Vector{String}}
annotation :: Dict{String, Union{Vector{String}, String}}
````
"""
mutable struct Gene
    id::String
    notes::Notes
    annotations::Annotations

    Gene(id::String = ""; notes = Notes(), annotations = Annotations()) =
        new(id, notes, annotations)
end


text:
struct S a end
struct T b end
function f(x::T) x.a end

max_stars_repo_name: gabriel-doriath-dohler/pjulia

max_issues_repo_name: gabriel-doriath-dohler/pjulia

max_forks_repo_name: gabriel-doriath-dohler/pjulia


text:
"""
    nonatomicdeferredacceptance_iid(students, students_dist, capacities;
                      verbose, rev, tol)

Nonatomic (continuous) analogue of `deferredacceptance()`, simplified to return only the score cutoffs
associated with each school, after Azevedo and Leshno (2016).

Students are a continuum of profiles distributed over a fixed set of student preference
lists, and school capacities are fractions of the total student population. Returns
only the cutoffs; use `assignmentfrompreflists()` to get the match array or `nonatomicdeferredacceptance()`
for a wrapper function.
"""
function nonatomicdeferredacceptance_iid(students         ::Union{AbstractArray{Int, 2}, AbstractArray{UInt, 2}},
                           students_dist    ::AbstractArray{<:AbstractFloat, 1},
                           capacities       ::AbstractArray{<:AbstractFloat, 1};
                           verbose          ::Bool=false,
                           rev              ::Bool=false,
                           tol              ::AbstractFloat=1e-12,
                           )::AbstractArray{<:AbstractFloat, 1}
    (m, n) = size(students)

    @assert (m,) == size(capacities)     "Dim mismatch between students and capacities"
    @assert (n,) == size(students_dist)  "Dim mismatch between students and students_dist"

    nit = 0
    done = false

    if !rev
        cutoffs = zeros(m)

        while done==false
            nit += 1
            done = true

            verbose ? println("\nRound $nit") : nothing

            demands = demandfrompreflists(students, students_dist, cutoffs)

            for (c, d) in enumerate(demands)
                if d - capacities[c] > tol
                    verbose ? println("  Demand at school $c was $d > capacity $(capacities[c])") : nothing
                    verbose ? println("    Old cutoff:  ", cutoffs[c]) : nothing
                    done = false
                    cutoffs[c] += (1 - capacities[c] / demands[c]) * (1 - cutoffs[c])
                    verbose ? println("    New cutoff:  ", cutoffs[c]) : nothing
                end
            end
        end

    else        # reverse
        # Selfish cutoff: Each school assumes it's everyone's first choice.
        # Thus this is the highest possible cutoff.
        cutoffs = max.(0., 1 .- capacities ./ sum(students_dist))

        while done==false
            nit += 1
            done = true

            verbose ? println("\nRound $nit") : nothing

            demands = demandfrompreflists(students, students_dist, cutoffs)

            for (c, d) in enumerate(demands)
                if cutoffs[c] > 0 && capacities[c] - d > tol   # If school has lots of remaining capacity
                    verbose ? println("  Demand at school $c was $d < capacity $(capacities[c])") : nothing
                    verbose ? println("    Old cutoff:  ", cutoffs[c]) : nothing
                    done = false
                    cutoffs[c] = max(cutoffs[c] + (1 - capacities[c] / demands[c]) * (1 - cutoffs[c]), 0)
                    verbose ? println("    New cutoff:  ", cutoffs[c]) : nothing
                end
            end
        end
    end

    return cutoffs
end

"""
    nonatomicdeferredacceptance_iid(demand, capacities;
                      verbose, rev, tol, maxit)

Nonatomic (continuous) analogue of `deferredacceptance()`, simplified to return only the score cutoffs
associated with each school, where demand is given by an arbitrary function that
takes score cutoffs as inputs. This algorithm assumes that student preferability
is independent at each school and that the demand function satisfies weak gross
substitutability. To relax the former assumption, use `nonatomictatonnement()`.
"""
function nonatomicdeferredacceptance_iid(demand      ::Function,
                           capacities  ::AbstractArray{<:AbstractFloat, 1};
                           verbose     ::Bool=false,
                           rev         ::Bool=false,
                           tol         ::AbstractFloat=1e-12,
                           maxit       ::Int=500,
                           )::AbstractArray{<:AbstractFloat, 1}

    (m, ) = size(capacities)

    if !rev
        cutoffs = zeros(m)

        for nit in 1:maxit
            done = true

            verbose ? println("\nRound $nit") : nothing

            demands = demand(cutoffs)

            for (c, d) in enumerate(demands)
                if d - capacities[c] > tol
                    verbose ? println("  Demand at school $c was $d > capacity $(capacities[c])") : nothing
                    verbose ? println("    Old cutoff:  ", cutoffs[c]) : nothing
                    done = false
                    cutoffs[c] += (1 - capacities[c] / demands[c]) * (1 - cutoffs[c])
                    verbose ? println("    New cutoff:  ", cutoffs[c]) : nothing
                end
            end

            if done == true
                break
            elseif nit == maxit
                @warn "Exceeded maximum number of iterations; try tuning parameters"
            end
        end

    else      # reverse
        # Selfish cutoff: Each school assumes it's everyone's first choice.
        # Thus this is the highest possible cutoff.
        cutoffs = max.(0., 1 .- capacities)

        for nit in 1:maxit
            done = true
            verbose ? println("\nRound $nit") : nothing

            demands = demand(cutoffs)

            for (c, d) in enumerate(demands)
                if cutoffs[c] > 0 && capacities[c] - d > tol   # If school has lots of remaining capacity
                    verbose ? println("  Demand at school $c was $d < capacity $(capacities[c])") : nothing
                    verbose ? println("    Old cutoff:  ", cutoffs[c]) : nothing
                    done = false
                    cutoffs[c] = max(cutoffs[c] + (1 - capacities[c] / demands[c]) * (1 - cutoffs[c]), 0)
                    verbose ? println("    New cutoff:  ", cutoffs[c]) : nothing
                end
            end

            if done == true
                break
            elseif nit == maxit
                @warn "Exceeded maximum number of iterations; try tuning parameters"
            end
        end

    end

    return cutoffs
end

"""
    nonatomicdeferredacceptance(students, students_dist, schools, capacities;
                 verbose, rev, return_cutoffs, tol)

Nonatomic (continuous) analogue of `deferredacceptance()`. Students are a continuum of profiles distributed
over a fixed set of student preference lists, and school capacities are fractions of the total
student population. School preferences are optional. If you pass `schools=nothing`, it assumes
that student preferability is uncorrelated with student profile, and the schools simply accept
the best students from each profile. This is the more realistic case; the algorithm is not
incentive compatible if schools can favor students based on the students' preferences.

When school preferences are not `nothing`, `rev` is a placeholder and reverse mode
has not yet been implemented.

Returns a tuple expressing the assignment as an array, where the `[i, j]`` entry indicates the
measure of students having preference list `j` admitted to school `i` and `i + 1` indicates nonassignment,
and a vector giving the average utility in each preference list.

Set `return_cutoffs=true` to get the score cutoffs associated with the match, after Azevedo
and Leshno (2016). These cutoffs are the state space of the algorithm and therefore more
numerically accurate than the assignment array itself. To get only the cutoffs, use
`nonatomicdeferredacceptance_iid()` (which this function wraps wheen `schools===nothing`)
or `nonatomictatonnement()` (which is not a DA algorithm but can compute equilibrium cutoffs in a wider range of scenarios).
"""
function nonatomicdeferredacceptance(students          ::Union{AbstractArray{Int, 2}, AbstractArray{UInt, 2}},
                      students_dist     ::AbstractArray{<:AbstractFloat, 1},
                      schools           ::Union{AbstractArray{Int, 2}, AbstractArray{UInt, 2}, Nothing},
                      capacities_in     ::AbstractArray{<:AbstractFloat, 1};
                      verbose           ::Bool=false,
                      rev               ::Bool=false,
                      return_cutoffs    ::Bool=false,
                      tol               ::AbstractFloat=1e-12)

    m, n = size(students)
    @assert (m,) == size(capacities_in)         "Dim mismatch between students and capacities"
    @assert (n,) == size(students_dist)         "Dim mismatch between students and students_dist"
    @assert schools === nothing || rev == false  "Reverse unavailable when schools have prefs"

    done = false
    nit = 0

    if schools === nothing                    # Homogenous student preferability
        students_inv = mapslices(invperm, students, dims=1)

        cutoffs = nonatomicdeferredacceptance_iid(students, students_dist, capacities_in;
                                    verbose=verbose, rev=rev, tol=tol)
        curr_assn = assignmentfrompreflists(students_inv, students_dist, cutoffs)

        verbose ? println("Cutoffs: ", cutoffs) : nothing
        DA_rank_dist = sum([col[students_inv[:, i]] for (i, col) in enumerate(eachcol(curr_assn))])
        append!(DA_rank_dist, sum(curr_assn[m + 1, :]))

        if return_cutoffs
            return curr_assn, DA_rank_dist, cutoffs
        else
            return curr_assn, DA_rank_dist
        end

    else            # Schools have preference order over student types
        @assert (n, m) == size(schools) "Dim mismatch between schools and students"

        students_inv = mapslices(invperm, students, dims=1)
        schools_inv = mapslices(invperm, schools, dims=1)
        if rev == false
            capacities = vcat(capacities_in, sum(students_dist))  # For students who never get assigned

            # Each entry indicates the volume of students from type j assigned to school i
            curr_assn = zeros(Float64, m + 1, n)
            for (s, C) in enumerate(eachcol(students_inv))
                curr_assn[C[1], s] = students_dist[s]
            end

            while !done
                nit += 1
                verbose ? println("\n\nRound $nit") : nothing
                done = true
                demands = sum(curr_assn, dims = 2)
                for c in 1:m    # Reject bin (c = m + 1) is never overdemanded
                    if demands[c] - capacities[c] > tol
                        capacity_remaining = capacities[c]
                        verbose ? print("\n  Total demand for school $c was ", demands[c],
                                ", but capacity is ", capacities[c]) : nothing
                        for s in filter(i->curr_assn[c, i] > 0, schools_inv[:, c])
                            if capacity_remaining ≥ curr_assn[c, s]
                                capacity_remaining -= curr_assn[c, s]
                                verbose ? print("\n    Accepting demand of ", curr_assn[c, s],
                                                " from profile $s; remaining capacity ",
                                                capacity_remaining) : nothing
                            else
                                done = false
                                prop_to_reject = 1 - capacity_remaining / curr_assn[c, s]
                                verbose ? print("\n    Demand from profile $s was ", curr_assn[c, s],
                                        "; rejecting $prop_to_reject of it") : nothing
                                next_school_id = get(students_inv, CartesianIndex(students[c, s] + 1, s), m + 1)
                                curr_assn[next_school_id, s] += prop_to_reject * curr_assn[c, s]
                                curr_assn[c, s] *= 1 - prop_to_reject
                                capacity_remaining = 0
                            end
                        end
                    end
                end
            end

            verbose ? println("\nDA terminated in $nit iterations") : nothing
            DA_rank_dist = sum([col[students_inv[:, i]] for (i, col) in enumerate(eachcol(curr_assn))])
            append!(DA_rank_dist, sum(curr_assn[m + 1, :]))

            return curr_assn, DA_rank_dist

        else
            error("No reverse mode yet when schools have pref over students")
        end
    end
end

"""
    nonatomicsecant(demand, capacities; verbose, tol, maxit)

Search for a market-clearing cutoff vector using a modified secant method with
tâtonnement as a fallback when change in demand is small. Less theoretically
legit than `nonatomictatonnement()` but tends to make fewer calls to `demand()`.

`demand` is a school demand function that takes cutoffs as input.
"""
function nonatomicsecant(demand      ::Function,
                          capacities  ::AbstractArray{<:AbstractFloat, 1};
                          verbose     ::Bool=false,
                          tol         ::AbstractFloat=1e-12,
                          maxit       ::Int=500,
                          )::AbstractArray{<:AbstractFloat, 1}

    (m, ) = size(capacities)

    UB = max.(0., 1 .- capacities)
    old_cutoffs = rand(m) .* UB
    new_cutoffs = rand(m) .* UB
    # old_cutoffs = zeros(m)
    # new_cutoffs = copy(UB)
    new_excess_demand = demand(old_cutoffs) - capacities

    for nit in 1:maxit
        verbose ? println("Round $nit") : nothing
        old_excess_demand, new_excess_demand = new_excess_demand, demand(new_cutoffs) - capacities

        verbose ? println("  Old cutoff vector: ", round.(old_cutoffs, digits = 4)) : nothing
        verbose ? println("  New cutoff vector: ", round.(new_cutoffs, digits = 4)) : nothing
        verbose ? println("  Excess demand: ", round.(new_excess_demand, digits = 4)) : nothing

        for c in 1:m
            # Fall back to tatonnement when denominator of secant update would
            # be close to zero.
            if isapprox(old_excess_demand[c], new_excess_demand[c], atol=tol)
                verbose ? println("    Updating school $c by tâtonnement because its demand didn't change") : nothing
                old_cutoffs[c], new_cutoffs[c] =
                    new_cutoffs[c], max(0, min(UB[c], new_cutoffs[c] + new_excess_demand[c]))
            else
                old_cutoffs[c], new_cutoffs[c] =
                    new_cutoffs[c],
                    max(0, min(UB[c],
                                 begin
                                     new_cutoffs[c] -
                                     new_excess_demand[c] * (new_cutoffs[c] - old_cutoffs[c]) /
                                     (new_excess_demand[c] - old_excess_demand[c])
                                 end ))
            end
        end

        if nit == maxit
            @warn "Exceeded maximum number of iterations; try tuning parameters"
        end

        # Check if market clears.
        if isapprox(new_cutoffs' * new_excess_demand, 0, atol=tol) &&
           all(new_excess_demand .≤ 0 + tol)
            break
        end
    end

    return cutbox(new_cutoffs)
end

"""
    nonatomictatonnement(demand, capacities; verbose, tol, β, maxit)

Search for a market-clearing cutoff vector using a modified tâtonnement process.
Analogous to `nonatomicdeferredacceptance_iid()` but more robust; `nonatomictatonnement()`
converges for any demand function that satisfies weak gross substitutability.

`demand` is a school demand function that takes cutoffs as input.

`β` is a step size parameter, where the step size is `1 / nit ^ β`, where `nit`
is the iteration number.
"""
function nonatomictatonnement(demand      ::Function,
                               capacities  ::AbstractArray{<:AbstractFloat, 1};
                               verbose     ::Bool=false,
                               tol         ::AbstractFloat=1e-12,
                               β           ::AbstractFloat=1e-3,
                               maxit       ::Int=500,
                               )::AbstractArray{<:AbstractFloat, 1}

    (m, ) = size(capacities)

    UB = max.(0., 1 .- capacities)
    cutoffs = rand(m) .* UB

    for nit in 1:maxit
        verbose ? println("Round $nit") : nothing
        excess_demand = demand(cutoffs) - capacities

        α = 1 / nit ^ β

        verbose ? println("  Cutoff vector: ", round.(cutoffs, digits = 4)) : nothing
        verbose ? println("  Excess demand: ", round.(excess_demand, digits = 4)) : nothing
        verbose ? println("  Step size:     ", α) : nothing

        new_cutoffs = max.(0, min.(UB, cutoffs + α * excess_demand))

        if nit == maxit
            @warn "Exceeded maximum number of iterations; try tuning parameters"
        end

        if isapprox(new_cutoffs, cutoffs, atol=tol)
            cutoffs = new_cutoffs
            break
        else
            cutoffs = new_cutoffs
        end
    end

    return cutoffs
end

"""
    ismarketclearing(students, students_dist, capacities, cutoffs;
                     tol=1e-6)

Check if a set of cutoffs is market clearing with respect to the given nonatomic
market. Nonatomic analogue of `isstable()` by Lemma 1 of Azevedo and Leshno (2016).
"""
function ismarketclearing(students      ::Union{AbstractArray{Int, 2}, AbstractArray{UInt, 2}},
                          students_dist ::AbstractArray{<:AbstractFloat, 1},
                          capacities    ::AbstractArray{<:AbstractFloat, 1},
                          cutoffs       ::AbstractArray{<:AbstractFloat, 1};
                          verbose       ::Bool=false,
                          tol           ::AbstractFloat=1e-6,
                          )::Bool
    demands = demandfrompreflists(students, students_dist, cutoffs)

    crit = falses(2)

    crit[1] = isapprox(sum(demands), min(sum(students_dist), sum(capacities)), atol=tol)
    crit[2] = all(demands .≤ capacities .+ tol)

    res = all(crit)

    if verbose
        for (test, pass) in zip(["Market clearing ",
                                 "School feas.    "], crit)
            println("$test:  $pass")
        end
    end

    return res
end

"""
    ismarketclearing(qualities, capacities, cutoffs;
                     tol=1e-6)

Check if a set of cutoffs is market clearing with respect to the given nonatomic
market. Nonatomic analogue of `isstable()` by Lemma 1 of Azevedo and Leshno (2016).
"""
function ismarketclearing(qualities     ::AbstractArray{<:AbstractFloat, 1},
                          capacities    ::AbstractArray{<:AbstractFloat, 1},
                          cutoffs       ::AbstractArray{<:AbstractFloat, 1};
                          verbose       ::Bool=false,
                          tol           ::AbstractFloat=1e-6,
                          )::Bool

    demands = demandfromMNL_iid(qualities, cutoffs)

    crit = falses(2)

    crit[1] = isapprox(sum(demands), min(1., sum(capacities)), atol=tol)
    crit[2] = all(demands .≤ capacities .+ tol)

    res = all(crit)

    if verbose
        for (test, pass) in zip(["Market clearing ",
                                 "School feas.    "], crit)
            println("$test:  $pass")
        end
    end

    return res
end

"""
    ismarketclearing(qualities, capacities, cutoffs;
                     tol=1e-6)

Check if a set of cutoffs is market clearing with respect to the given nonatomic
market. Nonatomic analogue of `isstable()` by Lemma 1 of Azevedo and Leshno (2016).
"""
function ismarketclearing(demand        ::Function,
                          capacities    ::AbstractArray{<:AbstractFloat, 1},
                          cutoffs       ::AbstractArray{<:AbstractFloat, 1};
                          verbose       ::Bool=false,
                          tol           ::AbstractFloat=1e-6,
                          )::Bool

    demands = demand(cutoffs)

    crit = falses(2)

    crit[1] = isapprox(sum(demands), min(1., sum(capacities)), atol=tol)
    crit[2] = all(demands .≤ capacities .+ tol)

    res = all(crit)

    if verbose
        for (test, pass) in zip(["Market clearing ",
                                 "School feas.    "], crit)
            println("$test:  $pass")
        end
    end

    return res
end

max_stars_repo_name: maxkapur/DeferredAcceptance

max_issues_repo_name: maxkapur/DeferredAcceptance

max_forks_repo_name: maxkapur/DeferredAcceptance


text:
using MVTs
using Base.Test

include("Label.jl")
include("Row.jl")
include("Table.jl")


text:
import RDates: WeekendCalendar, @rd_str, setcachedcalendar!, calendar
import Dates: Date

setcachedcalendar!(CALENDARS, "US/SETTLEMENT", WeekendCalendar() + RuleBasedCalendar([
    # New Year's Day
    PeriodicCalendarRule(rd"1y", rd"1JAN@WEEKEND[NEAR]"; cal_mgr=CALENDARS),
    # New Year's Day (handle case where day falls back into this year)
    PeriodicCalendarRule(rd"1y", rd"1y+1JAN@WEEKEND[NEAR]"; cal_mgr=CALENDARS),
    # # Martin Luther King Jr Birthday
    WindowCalendarRule(PeriodicCalendarRule(rd"1y", rd"1JAN+3rd MON"); from=Date(1983,1,1)),
    # # Washingthon's Birthday
    PeriodicCalendarRule(rd"1y", rd"1FEB+3rd MON"),
    # # Memorial Day
    PeriodicCalendarRule(rd"1y", rd"1MAY+Last MON"),
    # # Juneteenth
    WindowCalendarRule(PeriodicCalendarRule(rd"1y", rd"19JUN@WEEKEND[NEAR]"; cal_mgr=CALENDARS); from=Date(2021,1,1)),
    # # Independence Day
    PeriodicCalendarRule(rd"1y", rd"4JUL@WEEKEND[NEAR]"; cal_mgr=CALENDARS),
    # # Labor Day
    PeriodicCalendarRule(rd"1y", rd"1SEP+1st MON"),
    # # Columbus Day
    PeriodicCalendarRule(rd"1y", rd"1OCT+2nd MON"),
    # # Veteran's Day
    PeriodicCalendarRule(rd"1y", rd"11NOV@WEEKEND[NEAR]"; cal_mgr=CALENDARS),
    # # Thanksgiving Day
    PeriodicCalendarRule(rd"1y", rd"1NOV+4th THU"),
    # # Christmas
    PeriodicCalendarRule(rd"1y", rd"25DEC@WEEKEND[NEAR]"; cal_mgr=CALENDARS),
]))

setcachedcalendar!(CALENDARS, "NEW YORK", calendar(CALENDARS, "US/SETTLEMENT"))


text:
export WrightSampler

#=
This document implements a WrightSampler mutable struct type.
This object is used to generate Gibbs-sampler output on Wright equilibria
Also included are methods of the current state to generate stochastic gradients
And fisher information estimates
 =#

# This lists the neccesary parameters
# for sampling from the Wright-equilibrium distribution. The partition function is the most costly to compute and so it cahched
# θ_mut, θ_wt, σ are forward and backward diversity parameters respectively. At the moment, there is only code for  σ>0.
mutable struct WrightSampler
    θ_wt::Float64
    θ_mut::Float64
    σ::Float64
    k::Int64   # auxiliary variable for gibbs sampling
    x::Float64 # The freqeuncy of wt
end

function WrightSampler(θ_wt, θ_mut, σ)
#    z = gamma(θ_wt) * HypergeometricFunctions.pFq([θ_wt],[θ_mut + θ_wt],σ) /gamma(θ_wt +θ_mut)
    ws = WrightSampler(θ_wt, θ_mut, σ, 0, .5)
    for ii in 1:10
    incrementk!(ws)
    incrementx!(ws)
    end
    return ws
end

# This sampler is robuset for small-exponent gamma distributions
# Code was copied from somewhere unknown (citation needed)
function randlogGamma(a;scale = 1)
    if a > .2
        return log.(rand(Gamma(a)))
    else
        L = 1 / a - 1
        w = a / MathConstants.e / (1 - a)
        ww = 1 / (1 + w)
        η = z -> ( z >= 0 ? exp(-z) : w * L * exp(L * z))
        h = z -> exp(-z - exp(-z / a))
        function rh(a)
            while true
                U = rand()
                z = (U <= ww) ? -log(U / ww) : log(rand()) / L
                h(z) / η(z) > rand() && return(z)
            end
        end
        return log(scale) - rh(a) / a
    end
end

function sigmoid(x) # numerically stable sigmoid
    if x >= 0
        1 / ( 1 + exp(-x) )
    else
        exp(x) / ( 1 + exp(x) )
    end
end

function incrementk!(ws::WrightSampler)
    if ws.σ >= 0
        ws.k = rand(Poisson(abs(ws.σ * ws.x)))
    else
        ws.k = rand(Poisson(abs(ws.σ * (1-ws.x))))
    end
end

function incrementx!(ws::WrightSampler)
    if ws.σ >= 0
        a = randlogGamma(ws.θ_wt + ws.k;scale = 1)
        b = randlogGamma(ws.θ_mut;scale = 1)
    else
        a = randlogGamma(ws.θ_wt;scale = 1)
        b = randlogGamma(ws.θ_mut + ws.k;scale = 1)
    end
    ws.x = sigmoid(a-b)
end

function sample!(ws::WrightSampler)
    for ii in 1:5
    incrementk!(ws)
    incrementx!(ws)
    end
    return ws.x
end

function rejection_sample_k(θ_wt, θ_mut,σ)
    k = 0
    while true
        k = rand(Poisson(σ))
        k == 0 ? break : nothing
        p = exp( loggamma(θ_wt + θ_mut) + loggamma(k + θ_wt) - loggamma(k + θ_wt + θ_mut) -loggamma(θ_wt))
        rand(Bernoulli(p)) ? break : nothing
    end
    return k
end

max_stars_repo_name: StatPhysBio/EscapeSimulator

max_issues_repo_name: StatPhysBio/EscapeSimulator

max_forks_repo_name: StatPhysBio/EscapeSimulator


text:
module RayTraceGR

using DifferentialEquations
using Images
using LinearAlgebra
using StaticArrays

export Dual
struct Dual{T,DT} <: Real
    val::T
    eps::DT
end

@inline function Dual{T,DT}(val::T) where {T,DT}
    Dual{T,DT}(val, zeros(DT))
end
@inline function Dual{T,DT}(val) where {T,DT}
    Dual{T,DT}(val, zeros(DT))
end

Base.promote_rule(::Type{Dual{T,DT}}, ::Type{T}) where {T,DT} = Dual{T,DT}
Base.promote_rule(::Type{Dual{T,DT}}, ::Type{<:Integer}) where {T,DT} =
    Dual{T,DT}
Base.promote_rule(::Type{Dual{T,DT}}, ::Type{Dual{U,DU}}) where {T,DT,U,DU} =
    Dual{promote{T,U}, promote{DT,DU}}

# Cannot have a return type annotation here, as this would lead to an
# infinite recursion
@inline function Base.convert(::Type{Dual{T,DT}}, x::Dual{U,DU}
                              ) where {T,DT,U,DU}
    Dual{T,DT}(convert(T, x.val), convert(DT, x.eps))
end
@inline function Base.convert(::Type{Dual{T,DT}}, val::Number
                              )::Dual{T,DT} where {T,DT}
    Dual{T,DT}(convert(T, val))
end

@inline function Base.eps(::Type{Dual{T,DT}})::Dual{T,DT} where {T,DT}
    Dual{T,DT}(eps(T))
end

@inline function Base.zero(::Type{Dual{T,DT}})::Dual{T,DT} where {T,DT}
    Dual{T,DT}(zero(T))
end
@inline function Base.one(::Type{Dual{T,DT}})::Dual{T,DT} where {T,DT}
    Dual{T,DT}(one(T))
end

@inline function Base.:+(x::Dual{T,DT})::Dual{T,DT} where {T,DT}
    Dual{T,DT}(+x.val, +x.eps)
end

@inline function Base.:-(x::Dual{T,DT})::Dual{T,DT} where {T,DT}
    Dual{T,DT}(-x.val, -x.eps)
end

@inline function Base.:+(x::Dual{T,DT}, y::Dual{T,DT})::Dual{T,DT} where {T,DT}
    Dual{T,DT}(x.val+y.val, x.eps+y.eps)
end
@inline function Base.:+(x::Dual{T,DT}, a::T)::Dual{T,DT} where {T,DT}
    Dual{T,DT}(x.val+ a, x.eps)
end
@inline function Base.:+(a::T, x::Dual{T,DT})::Dual{T,DT} where {T,DT}
    Dual{T,DT}(a+x.val, x.eps)
end
@inline function Base.:+(x::Dual{T,DT}, a::Integer)::Dual{T,DT} where {T,DT}
    Dual{T,DT}(x.val+a, x.eps)
end
@inline function Base.:+(a::Integer, x::Dual{T,DT})::Dual{T,DT} where {T,DT}
    Dual{T,DT}(a+x.val, x.eps)
end

@inline function Base.:-(x::Dual{T,DT}, y::Dual{T,DT})::Dual{T,DT} where {T,DT}
    Dual{T,DT}(x.val-y.val, x.eps-y.eps)
end
@inline function Base.:-(x::Dual{T,DT}, a::T)::Dual{T,DT} where {T,DT}
    Dual{T,DT}(x.val- a, x.eps)
end
@inline function Base.:-(a::T, x::Dual{T,DT})::Dual{T,DT} where {T,DT}
    Dual{T,DT}(a-x.val, x.eps)
end
@inline function Base.:-(x::Dual{T,DT}, a::Integer)::Dual{T,DT} where {T,DT}
    Dual{T,DT}(x.val-a, x.eps)
end
@inline function Base.:-(a::Integer, x::Dual{T,DT})::Dual{T,DT} where {T,DT}
    Dual{T,DT}(a-x.val, x.eps)
end

@inline function Base.:*(x::Dual{T,DT}, y::Dual{T,DT})::Dual{T,DT} where {T,DT}
    Dual{T,DT}(x.val*y.val, x.eps.*y.val + x.val.*y.eps)
end
@inline function Base.:*(x::Dual{T,DT}, a::T)::Dual{T,DT} where {T,DT}
    Dual{T,DT}(x.val*a, x.eps.*a)
end
@inline function Base.:*(a::T, x::Dual{T,DT})::Dual{T,DT} where {T,DT}
    Dual{T,DT}(a*x.val, a.*x.eps)
end
@inline function Base.:*(x::Dual{T,DT}, a::Integer)::Dual{T,DT} where {T,DT}
    Dual{T,DT}(x.val*a, x.eps.*a)
end
@inline function Base.:*(a::Integer, x::Dual{T,DT})::Dual{T,DT} where {T,DT}
    Dual{T,DT}(a*x.val, a.*x.eps)
end

# 1/(a+ϵb) = (a-ϵb)/a^2
@inline function Base.inv(x::Dual{T,DT})::Dual{T,DT} where {T,DT}
    Dual{T,DT}(inv(x.val), -inv(x.val)^2 .* x.eps)
end

@inline function Base.:/(x::Dual{T,DT}, y::Dual{T,DT})::Dual{T,DT} where {T,DT}
    Dual{T,DT}(x.val/y.val, (x.eps * y.val - x.val * y.eps) / y.val^2)
end
@inline function Base.:/(x::Dual{T,DT}, a::T)::Dual{T,DT} where {T,DT}
    Dual{T,DT}(x.val/a, x.eps./a)
end
@inline function Base.:/(x::Dual{T,DT}, a::Integer)::Dual{T,DT} where {T,DT}
    Dual{T,DT}(x.val/a, x.eps./a)
end

@inline function Base.:\(x::Dual{T,DT}, y::Dual{T,DT})::Dual{T,DT} where {T,DT}
    Dual{T,DT}(x.val\y.val, x.val^2 \ (x.eps .* y.val - x.val .* y.eps))
end
@inline function Base.:\(a::T, x::Dual{T,DT})::Dual{T,DT} where {T,DT}
    Dual{T,DT}(a\x.val, a.\x.eps)
end
@inline function Base.:\(a::Integer, x::Dual{T,DT})::Dual{T,DT} where {T,DT}
    Dual{T,DT}(a\x.val, a.\x.eps)
end

@inline Base.literal_pow(::typeof(^), x::Dual, ::Val{0}) = one(x)
@inline Base.literal_pow(::typeof(^), x::Dual, ::Val{1}) = x
@inline Base.literal_pow(::typeof(^), x::Dual, ::Val{2}) = x*x
@inline Base.literal_pow(::typeof(^), x::Dual, ::Val{3}) = x*x*x
@inline Base.literal_pow(::typeof(^), x::Dual, ::Val{4}) = (x^2)^2

@inline function Base.:^(x::Dual{T,DT}, n::Integer)::Dual{T,DT} where {T,DT}
    n == 0 && return one(x)
    Dual{T,DT}(x.val^n, n * x.val^(n-1) .* x.eps)
end
@inline function Base.:^(x::Dual{T,DT}, a::T)::Dual{T,DT} where {T,DT}
    Dual{T,DT}(x.val^a, a * x.val^(a-1) .* x.eps)
end
@inline function Base.:^(x::Dual{T,DT}, y::Dual{T,DT})::Dual{T,DT} where {T,DT}
    r = x.val ^ y.val
    Dual{T,DT}(r, y.val/x.val * r .* x.eps + r * log(x.val) .* y.eps)
end

@inline function Base.abs(x::Dual{T,DT})::Dual{T,DT} where {T,DT}
    Dual{T,DT}(abs(x.val), copysign(one(T), x.val) .* x.eps)
end

function Base.acos(x::Dual{T,DT})::Dual{T,DT} where {T,DT}
    Dual{T,DT}(acos(x.val), -1 / sqrt(1 - x.val^2) .* x.eps)
end

function Base.asin(x::Dual{T,DT})::Dual{T,DT} where {T,DT}
    Dual{T,DT}(asin(x.val), 1 / sqrt(1 - x.val^2) .* x.eps)
end

function Base.atan(x::Dual{T,DT})::Dual{T,DT} where {T,DT}
    Dual{T,DT}(atan(x.val), 1 / (1 + x.val^2) .* x.eps)
end
function Base.atan(y::Dual{T,DT}, x::Dual{T,DT}
                           )::Dual{T,DT} where {T,DT}
    ρ2 = x.val^2 + y.val^2
    Dual{T,DT}(atan(y.val, x.val), x.val .* y.eps - y.val / ρ2 .* x.eps)
end

@inline function Base.cbrt(x::Dual{T,DT})::Dual{T,DT} where {T,DT}
    r = cbrt(x.val)
    Dual{T,DT}(r, r/(3*x.val) .* x.eps)
end

function Base.cos(x::Dual{T,DT})::Dual{T,DT} where {T,DT}
    Dual{T,DT}(cos(x.val), -sin(x.val) .* x.eps)
end

function Base.exp(x::Dual{T,DT})::Dual{T,DT} where {T,DT}
    r = exp(x.val)
    Dual{T,DT}(r, r .* x.eps)
end

function Base.log(x::Dual{T,DT})::Dual{T,DT} where {T,DT}
    Dual{T,DT}(log(x.val), 1 / x.val .* x.eps)
end

function Base.sin(x::Dual{T,DT})::Dual{T,DT} where {T,DT}
    Dual{T,DT}(sin(x.val), cos(x.val) .* x.eps)
end

@inline function Base.sqrt(x::Dual{T,DT})::Dual{T,DT} where {T,DT}
    r = sqrt(x.val)
    Dual{T,DT}(r, 1/(2*r) .* x.eps)
end

@inline function Base.:(==)(x::Dual{T,DT}, y::Dual{T,DT})::Bool where {T,DT}
    x.val == y.val
end
@inline function Base.:(==)(x::Dual{T,DT}, a::T)::Bool where {T,DT}
    x.val == a
end
@inline function Base.:(==)(x::Dual{T,DT}, a::Integer)::Bool where {T,DT}
    x.val == a
end
@inline function Base.:(==)(a::T, x::Dual{T,DT})::Bool where {T,DT}
    a == x.val
end
@inline function Base.:(==)(a::Integer, x::Dual{T,DT})::Bool where {T,DT}
    a == x.val
end

@inline function Base.:(<)(x::Dual{T,DT}, y::Dual{T,DT})::Bool where {T,DT}
    x.val < y.val
end
@inline function Base.:(<)(x::Dual{T,DT}, a::T)::Bool where {T,DT}
    x.val < a
end
@inline function Base.:(<)(x::Dual{T,DT}, a::Integer)::Bool where {T,DT}
    x.val < a
end
@inline function Base.:(<)(a::T, x::Dual{T,DT})::Bool where {T,DT}
    a < x.val
end
@inline function Base.:(<)(a::Integer, x::Dual{T,DT})::Bool where {T,DT}
    a < x.val
end

@inline function Base.isinf(x::Dual{T,DT})::Bool where {T,DT}
    isinf(x.val)
end
@inline function Base.isnan(x::Dual{T,DT})::Bool where {T,DT}
    isnan(x.val) || any(isnan.(x.eps))
end

@inline function Base.hash(x::Dual{T,DT}, h::UInt)::UInt where {T,DT}
    hash(0xdccda268, hash(x.val, hash(x.eps, h)))
end
@inline function Base.isequal(x::Dual{T,DT}, y::Dual{T,DT})::Bool where {T,DT}
    isequal(x.val, y.val) && isequal(x.eps, y.eps)
end
@inline function Base.isless(x::Dual{T,DT}, y::Dual{T,DT})::Bool where {T,DT}
    isless(x.val, y.val) && return true
    isless(y.val, x.val) && return false
    isless(x.eps, y.eps)
end

################################################################################

export D
const D = 4

export minkowski
"""
Minkowski metric g_ab
"""
function minkowski(x::SVector{D,T})::SMatrix{D,D,T} where {T}
    @SMatrix T[a==b ? (a==1 ? -1 : 1) : 0 for a in 1:D, b in 1:D]
end

export kerr_schild
"""
Kerr-Schild metric g_ab

Living Reviews in Relativity, Greg Cook, 2000, section 3.3.1
"""
function kerr_schild(xx::SVector{D,T})::SMatrix{D,D,T} where {T}
    M = 1
    a = 0                       # T(0.8)

    t,x,y,z = xx
    @assert !any(isnan, (t, x, y, z))

    # <https://en.wikipedia.org/wiki/Kerr_metric>
    η = @SMatrix T[a==b ? (a==1 ? -1 : 1) : 0 for a in 1:D, b in 1:D]
    ρ = sqrt(x^2 + y^2 + z^2)
    r = sqrt(ρ^2 - a^2)/2 + sqrt(a^2*z^2 + ((ρ^2 - a^2)/2)^2)
    f = 2*M*r^3 / (r^4 + a^2*z^2)
    k = SVector{D,T}(1,
                     (r*x + a*y) / (r^2 + a^2),
                     (r*y - a*x) / (r^2 + a^2),
                     z / r)

    g = @SMatrix T[η[a,b] + f * k[a] * k[b] for a in 1:D, b in 1:D]

    g
end

export dmetric
"""
Derivative of four-metric g_ab,c
"""
function dmetric(metric::Metric, x::SVector{D,T}) where {Metric, T}
    DT = SVector{D,T}
    TDT = Dual{T,DT}
    xdx = SVector{D,TDT}(TDT(x[1], SVector{D,T}(1, 0, 0, 0)),
                         TDT(x[2], SVector{D,T}(0, 1, 0, 0)),
                         TDT(x[3], SVector{D,T}(0, 0, 1, 0)),
                         TDT(x[4], SVector{D,T}(0, 0, 0, 1)))
    gdg = metric(xdx)
    g = @SMatrix T[gdg[a,b].val for a in 1:D, b in 1:D]
    dg = @SArray T[gdg[a,b].eps[c] for a in 1:D, b in 1:D, c in 1:D]
    g, dg
end

export christoffel
"""
Christoffel symbol Γ^a_bc
"""
function christoffel(metric::Metric, x::SVector{D,T}) where {Metric, T}
    g,dg = dmetric(metric, x)
    gu = inv(g)
    Γl = @SArray T[(dg[a,b,c] + dg[a,c,b] - dg[b,c,a]) / 2
                   for a in 1:D, b in 1:D, c in 1:D]
    @SArray T[gu[a,1] * Γl[1,b,c] +
              gu[a,2] * Γl[2,b,c] +
              gu[a,3] * Γl[3,b,c] +
              gu[a,4] * Γl[4,b,c]
              for a in 1:D, b in 1:D, c in 1:D]
end

export Ray
"""
State of a ray: position x, velocity u
"""
struct Ray{T}
    x::SVector{D,T}
    u::SVector{D,T}
end

export r2s, s2r
function r2s(r::Ray{T})::SVector{2D,T} where {T}
    @SVector T[a<=D ? r.x[a] : r.u[a-D] for a in 1:2D]
end
function s2r(s::SVector{2D,T})::Ray{T} where {T}
    x = @SVector T[s[a] for a in 1:D]
    u = @SVector T[s[D+a] for a in 1:D]
    Ray{T}(x, u)
end

export geodesic
"""
RHS of geodesic equation for ray p
"""
function geodesic(r::Ray{T}, metric::Metric, λ::T)::Ray{T} where {Metric, T}
    Γ = christoffel(metric, r.x)
    xdot = r.u
    udot = @SVector T[- sum(@SMatrix T[Γ[a,x,y] * r.u[x] * r.u[y]
                                       for x in 1:D, y in 1:D])
                      for a in 1:D]
    Ray{T}(xdot, udot)
end

function geodesic(s::SVector{2D,T}, metric::Metric,
                  λ::T)::SVector{2D,T} where {Metric, T}
    r2s(geodesic(s2r(s), metric, λ))
end

export Object
abstract type Object{T} end

"""
distance between object and point

The distance does not need to be the geodesic distance; any distance
measure is fine, as long as it is zero on the surface, positive
outside, and negative inside.
"""
function distance(obj::Object, pos::SVector{D,T})::T where{T}
    @error "Called distance on abstract object"
end
function objcolor(obj::Object, pos::SVector{D,T})::SVector{3,T} where {T}
    @error "Called color on abstract object"
end

export Plane
struct Plane{T} <: Object{T}
    # TOOD: add normal
    time::T
end

function distance(pl::Plane{T}, pos::SVector{D,T})::T where{T}
    pos[1] - pl.time
end
function objcolor(pl::Plane{T}, pos::SVector{D,T})::SVector{3,T} where{T}
    SVector{3,T}(0, T(1)/2, 0)
end

export Sphere
struct Sphere{T} <: Object{T}
    pos::SVector{D,T}
    vel::SVector{D,T}
    radius::T
end

function distance(sph::Sphere{T}, pos::SVector{D,T})::T where{T}
    # TODO: Use metric?
    sign(sph.radius) *
        (sum(@SVector T[(pos[a] - sph.pos[a])^2 for a in 2:D]) - sph.radius^2)
end
function objcolor(sph::Sphere{T}, pos::SVector{D,T})::SVector{3,T} where{T}
    x = pos[2] - sph.pos[2]
    y = pos[3] - sph.pos[3]
    z = pos[4] - sph.pos[4]
    r = sqrt(x^2 + y^2 + z^2)
    θ = acos(z / r)
    ϕ = atan(y, x)
    SVector{3,T}(mod(12 * θ / π, 1), mod(12 * ϕ / π, 1), T(1))
end

export min_distance
function min_distance(objs::Vector{Object{T}}, s::SVector{2D,T})::T where {T}
    r = s2r(s)
    dmin = T(Inf)
    for obj in objs
        d = distance(obj, r.x)::T
        dmin = min(dmin, d)
    end
    dmin
end

export Pixel
struct Pixel{T}
    pos::SVector{D,T}         # x^a
    normal::SVector{D,T}      # n^a, past null normal to pixel surface
    rgb::SVector{3,T}
end

export Canvas
struct Canvas{T}
    pixels::Array{Pixel{T}, 2}
end

export make_canvas
function make_canvas(metric::Metric,
                     pos::SVector{D,T},
                     widthx::SVector{D,T}, widthy::SVector{D,T},
                     normal::SVector{D,T},
                     ni::Int, nj::Int) where {Metric, D, T}
    pixels = Array{Pixel{T}}(undef, ni, nj)
    for j in 1:nj, i in 1:ni
        dx = (i-T(1)/2) / ni - T(1)/2
        dy = (j-T(1)/2) / nj - T(1)/2
        x = pos + dx * widthx + dy * widthy
        n = normal + dx * widthx + dy * widthy
        g = metric(x)
        gu = inv(g)
        t = gu * SVector{D,T}(1, 0, 0, 0)
        t2 = t' * g * t
        n2 = n' * g * n
        u = (t / sqrt(-t2) + n / sqrt(n2)) / sqrt(T(2))
        pixels[i,j] = Pixel{T}(x, u, zeros(SVector{3,T}))
    end
    Canvas{T}(pixels)
end

export trace_rays
function trace_rays(metric, objs::Vector{Object{T}},
                    c::Canvas{T})::Canvas{T} where {T}
    tol = eps(T)^(T(3)/4)
    l = length(c.pixels)

    condition(s, λ, integrator) = min_distance(objs, s)
    affect!(integrator) = terminate!(integrator)
    cb = ContinuousCallback(condition, affect!)

    function input_func(i::Int)
        p = c.pixels[i]
        x,u = p.pos,p.normal
        r2s(Ray{T}(x, u))
    end
    λ0,λ1 = T(0),T(100)
    prob = ODEProblem(geodesic, input_func(1), (λ0, λ1), metric)
    function prob_func(prob,i,repeat)
        ODEProblem(geodesic, input_func(i), (λ0, λ1), metric)
    end
    function output_func(sol,i)
        λend = sol.t[end]
        s2r(sol[end]), false
    end

    probs = EnsembleProblem(prob,
                            prob_func=prob_func,
                            output_func=output_func)
    sols = solve(probs, Tsit5(), callback=cb, trajectories=length(c.pixels),
                 reltol=tol, abstol=tol)

    pixels = Array{Pixel{T}}(undef, size(c.pixels))
    for i in 1:length(pixels)
        p = c.pixels[i]
        r = sols.u[i]
        x,u = r.x,r.u
        omin = 0
        dmin = T(0.01)
        for (o,obj) in enumerate(objs)
            d = distance(obj, x)
            if d < dmin
                omin = o
                dmin = d
            end
        end
        if omin == 0
            col = SVector{3,T}(1, 0, 0) # did not hit any object
        else
            col = objcolor(objs[omin], x) .* (T(omin) / length(objs))
        end
        pixels[i] = Pixel{T}(p.pos, p.normal, col)
    end

    Canvas{T}(pixels)
end

const outdir = "scenes"

function example1()
    T = Float64

    metric = minkowski
    caelum = Sphere{T}(SVector{D,T}(0, 0, 0, 0), SVector{D,T}(1, 0, 0, 0), -10)
    frustum = Plane{T}(-20)
    sphere = Sphere{T}(SVector{D,T}(0, 0, 0, 0), SVector{D,T}(1, 0, 0, 0),
                       T(1)/2)
    objs = Object{T}[caelum, frustum, sphere]

    ni = 200
    nj = 200
    pos = SVector{D,T}(0, 0, -2, 0)
    widthx = SVector{D,T}(0, 1, 0, 0)
    widthy = SVector{D,T}(0, 0, 0, 1)
    normal = SVector{D,T}(0, 0, 1, 0)
    canvas = make_canvas(metric, pos, widthx, widthy, normal, ni, nj)

    canvas = trace_rays(metric, objs, canvas)

    # scene = colorview(RGB,
    #                   T[mod(p.rgb[1] + 1/2, 1) for p in canvas.pixels]',
    #                   T[mod(p.rgb[2] + 1/2, 1) for p in canvas.pixels]',
    #                   T[mod(p.rgb[3] + 1/2, 1) for p in canvas.pixels]')
    scene = colorview(RGB,
                      T[p.rgb[1] for p in canvas.pixels]',
                      T[p.rgb[2] for p in canvas.pixels]',
                      T[p.rgb[3] for p in canvas.pixels]')

    mkpath(outdir)
    file = joinpath(outdir, "sphere.png")
    rm(file, force=true)
    println("Output file is \"$file\"")
    save(file, scene)
end

function example2()
    T = Float64

    metric = kerr_schild
    caelum = Sphere{T}(SVector{D,T}(0, 0, 0, 0), SVector{D,T}(1, 0, 0, 0), -10)
    frustum = Plane{T}(-20)
    sphere = Sphere{T}(SVector{D,T}(0, 4, 0, 0), SVector{D,T}(1, 0, 0, 0),
                       T(1)/2)
    objs = Object{T}[caelum, frustum, sphere]

    ni = 200
    nj = 200
    pos = SVector{D,T}(0, 4, -2, 0)
    widthx = SVector{D,T}(0, 1, 0, 0)
    widthy = SVector{D,T}(0, 0, 0, 1)
    normal = SVector{D,T}(0, 0, 1, 0)
    canvas = make_canvas(metric, pos, widthx, widthy, normal, ni, nj)

    canvas = trace_rays(metric, objs, canvas)

    # scene = colorview(RGB,
    #                   T[mod(p.rgb[1] + 1/2, 1) for p in canvas.pixels]',
    #                   T[mod(p.rgb[2] + 1/2, 1) for p in canvas.pixels]',
    #                   T[mod(p.rgb[3] + 1/2, 1) for p in canvas.pixels]')
    scene = colorview(RGB,
                      T[p.rgb[1] for p in canvas.pixels]',
                      T[p.rgb[2] for p in canvas.pixels]',
                      T[p.rgb[3] for p in canvas.pixels]')

    mkpath(outdir)
    file = joinpath(outdir, "sphere2.png")
    rm(file, force=true)
    println("Output file is \"$file\"")
    save(file, scene)
end

end


text:
# those do not need to dispatch I guess
struct SemLasso{P, W} <: LossFunction
    penalty::P
    which::W
end

function (lasso::SemLasso)(par, model)
    F = lasso.penalty*sum(transpose(par)[lasso.which])
    return F
end

struct SemRidge{P, W} <: LossFunction
    penalty::P
    which::W
end

function (ridge::SemRidge)(par, model)
    F = ridge.penalty*sum(transpose(par)[ridge.which].^2)
    return F
end

# function (lasso::SemLasso)(par, model)
# end
#
# function (lasso::SemLasso)(par, model, E, G) where {G <: Nothing}
#     lasso.F .= lasso.penalty*sum(transpose(par)[lasso.which])
#     return lasso.F[1]
# end
#
# function (lasso::SemLasso)(par, model, E, G) where {E <: Nothing}
#     model.imply()
#     ForwardDiff.gradient!(G, lasso(par, model), par)
# end