text:
import DataValues: DataValue

_is_subtype(::Type{S}, ::Type{T}) where {S, T} = promote_type(S, T) == T

Base.@pure function dataarrayof(T)
    if T<:DataValue
        DataValueArray{T.parameters[1],1}
    else
        Vector{T}
    end
end

"""
`collect_columns(itr)`

Collect an iterable as a `Columns` object if it iterates `Tuples` or `NamedTuples`, as a normal
`Array` otherwise.

## Examples

```jldoctest collect
julia> s = [(1,2), (3,4)];

julia> collect_columns(s)
2-element Columns{Tuple{Int64,Int64}}:
 (1, 2)
 (3, 4)

 julia> s = Iterators.filter(isodd, 1:8);

 julia> collect_columns(s)
 4-element Array{Int64,1}:
  1
  3
  5
  7
```
"""
collect_columns(itr) = collect_columns(itr, Base.iteratorsize(itr))

function collect_empty_columns(itr::T) where {T}
    S = Core.Inference.return_type(first, Tuple{T})
    similar(arrayof(S), 0)
end

function collect_columns(itr, ::Union{Base.HasShape, Base.HasLength})
    st = start(itr)
    done(itr, st) && return collect_empty_columns(itr)
    el, st = next(itr, st)
    dest = similar(arrayof(typeof(el)), length(itr))
    dest[1] = el
    collect_to_columns!(dest, itr, 2, st)
end

function collect_to_columns!(dest::AbstractArray{T}, itr, offs, st = start(itr)) where {T}
    # collect to dest array, checking the type of each result. if a result does not
    # match, widen the result type and re-dispatch.
    i = offs
    while !done(itr, st)
        el, st = next(itr, st)
        if fieldwise_isa(el, T)
            @inbounds dest[i] = el
            i += 1
        else
            new = widencolumns(dest, i, el, T)
            @inbounds new[i] = el
            return collect_to_columns!(new, itr, i+1, st)
        end
    end
    return dest
end

function collect_columns(itr, ::Base.SizeUnknown)
    st = start(itr)
    done(itr, st) && return collect_empty_columns(itr)
    el, st = next(itr, st)
    dest = similar(arrayof(typeof(el)), 1)
    dest[1] = el
    grow_to_columns!(dest, itr, st)
end

function collect_columns_flattened(itr)
    st = start(itr)
    el, st = next(itr, st)
    collect_columns_flattened(itr, el, st)
end

function collect_columns_flattened(itr, el, st)
    while isempty(el)
        done(itr, st) && return collect_empty_columns(el)
        el, st = next(itr, st)
    end
    dest = collect_columns(el)
    collect_columns_flattened!(dest, itr, el, st)
end

function collect_columns_flattened!(dest, itr, el, st)
    while !done(itr, st)
        el, st = next(itr, st)
        dest = grow_to_columns!(dest, el)
    end
    return dest
end

function collect_columns_flattened(itr, el::Pair, st)
    while isempty(el.second)
        done(itr, st) && return collect_empty_columns(el.first => i for i in el.second)
        el, st = next(itr, st)
    end
    dest_data = collect_columns(el.second)
    dest_key = collect_columns(el.first for i in dest_data)
    collect_columns_flattened!(Columns(dest_key => dest_data), itr, el, st)
end

function collect_columns_flattened!(dest::Columns{<:Pair}, itr, el::Pair, st)
    dest_key, dest_data = dest.columns
    while !done(itr, st)
        el, st = next(itr, st)
        n = length(dest_data)
        dest_data = grow_to_columns!(dest_data, el.second)
        dest_key = grow_to_columns!(dest_key, el.first for i in (n+1):length(dest_data))
    end
    return Columns(dest_key => dest_data)
end

function grow_to_columns!(dest::AbstractArray{T}, itr, st = start(itr)) where {T}
    # collect to dest array, checking the type of each result. if a result does not
    # match, widen the result type and re-dispatch.
    i = length(dest)+1
    while !done(itr, st)
        el, st = next(itr, st)
        if fieldwise_isa(el, T)
            push!(dest, el)
            i += 1
        else
            new = widencolumns(dest, i, el, T)
            push!(new, el)
            return grow_to_columns!(new, itr, st)
        end
    end
    return dest
end

# extra methods if we have widened to Vector{Tuple} or Vector{NamedTuple}
# better to not generate as this is the case where the user is sending heterogenoeus data
fieldwise_isa(el::S, ::Type{Tuple}) where {S<:Tup} = _is_subtype(S, Tuple)
fieldwise_isa(el::S, ::Type{NamedTuple}) where {S<:Tup} = _is_subtype(S, NamedTuple)

@generated function fieldwise_isa(el::S, ::Type{T}) where {S<:Tup, T<:Tup}
    if (fieldnames(S) == fieldnames(T)) && all(_is_subtype(s, t) for (s, t) in zip(S.parameters, T.parameters))
        return :(true)
    else
        return :(false)
    end
end

@generated function fieldwise_isa(el::S, ::Type{T}) where {S, T}
    if _is_subtype(S, T)
        return :(true)
    else
        return :(false)
    end
end

fieldwise_isa(el::Pair, ::Type{Pair{T1, T2}}) where {T1, T2}  =
    fieldwise_isa(el.first, T1) && fieldwise_isa(el.second, T2)

function widencolumns(dest, i, el::S, ::Type{T}) where{S <: Tup, T<:Tup}
    if fieldnames(S) != fieldnames(T) || T == Tuple || T == NamedTuple
        R = (S <: Tuple) && (T <: Tuple) ? Tuple :  (S <: NamedTuple) && (T <: NamedTuple) ? NamedTuple : Any
        new = Array{R}(length(dest))
        copy!(new, 1, dest, 1, i-1)
    else
        sp, tp = S.parameters, T.parameters
        idx = find(!(s <: t) for (s, t) in zip(sp, tp))
        new = dest
        for l in idx
            newcol = dataarrayof(promote_type(sp[l], tp[l]))(length(dest))
            copy!(newcol, 1, column(dest, l), 1, i-1)
            new = setcol(new, l, newcol)
        end
    end
    new
end

function widencolumns(dest, i, el::S, ::Type{T}) where{S, T}
    new = dataarrayof(promote_type(S, T))(length(dest))
    copy!(new, 1, dest, 1, i-1)
    new
end

function widencolumns(dest::Columns{<:Pair}, i, el::Pair, ::Type{Pair{T1, T2}}) where{T1, T2}
    dest1 = fieldwise_isa(el.first, T1) ? dest.columns.first : widencolumns(dest.columns.first, i, el.first, T1)
    dest2 = fieldwise_isa(el.second, T2) ? dest.columns.second : widencolumns(dest.columns.second, i, el.second, T2)
    Columns(dest1 => dest2)
end


text:
using Documenter, TightBinding

makedocs(;
    modules=[TightBinding],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/cometscome/TightBinding.jl/blob/{commit}{path}#L{line}",
    sitename="TightBinding.jl",
    authors="Yuki Nagai",
    assets=[],
)

deploydocs(;
    repo="github.com/cometscome/TightBinding.jl",
)

ext: jl


text:
#######################################################################################################################################################################################################
#
# Changes made to this function
# General
#     2021-Oct-22: add function to compute leaf level PAR and APAR
#     2022-Jan-13: use LeafBiophysics directly in the function rather than Leaf
#     2022-Feb-02: fix documentation
#     2022-Feb-02: unpack CONSTANTS only
#
#######################################################################################################################################################################################################
"""

    leaf_PAR(bio::LeafBiophysics{FT}, wls::WaveLengthSet{FT}, rad::HyperspectralRadiation{FT}; APAR_car::Bool = true) where {FT<:AbstractFloat}

Return leaf level PAR and APAR, given
- `bio` `ClimaCache.LeafBiophysics` type struct that contains leaf biophysical parameters
- `wls` `ClimaCache.WaveLengthSet` type struct that contains wave length bins
- `rad` `ClimaCache.HyperspectralRadiation` type struct that contains incoming radiation information
- `APAR_car` If true (default), account carotenoid absorption as APAR; otherwise, APAR is only by chlorophyll

---
# Examples
```julia
wls = WaveLengthSet{Float64}();
bio = LeafBiophysics{Float64}(wls);
rad = HyperspectralRadiation{Float64}();
par,apar = leaf_PAR(bio, wls, rad);
par,apar = leaf_PAR(bio, wls, rad; APAR_car=false);
```
"""
function leaf_PAR(bio::LeafBiophysics{FT}, wls::WaveLengthSet{FT}, rad::HyperspectralRadiation{FT}; APAR_car::Bool = true) where {FT<:AbstractFloat}
    @unpack IΛ_PAR, ΔΛ_PAR, Λ_PAR = wls;

    # APAR absorption feature
    _α = (APAR_car ? view(bio.α_cabcar, IΛ_PAR) : view(bio.α_cab, IΛ_PAR));

    # PAR energy from direct  and diffuse light
    _e_par_dir  = view(rad.e_direct , IΛ_PAR) .* view(bio.α_SW, IΛ_PAR);
    _e_par_diff = view(rad.e_diffuse, IΛ_PAR) .* view(bio.α_SW, IΛ_PAR);
    _par_dir    = photon.(Λ_PAR, _e_par_dir);
    _par_diff   = photon.(Λ_PAR, _e_par_diff);

    # absorbed PAR energy from direct and diffuse light
    _apar_dir  = _α .* _par_dir;
    _apar_diff = _α .* _par_diff;

    # total PAR and APAR in μmol photons m⁻² s⁻¹
    _∑par_dir   = numerical∫(_par_dir, ΔΛ_PAR);
    _∑par_diff  = numerical∫(_par_diff, ΔΛ_PAR);
    _∑apar_dir  = numerical∫(_apar_dir, ΔΛ_PAR);
    _∑apar_diff = numerical∫(_apar_diff, ΔΛ_PAR);

    return 1000 * (_∑par_dir + _∑par_diff), 1000 * (_∑apar_dir + _∑apar_diff)
end


text:
@testset "LLS" begin
  # least-square problem without constraints
  m, n = 20, 10
  A = rand(m, n)
  b = rand(m)
  nls = LLSModel(A, b)
  statsnls = ripqp(nls, itol = InputTol(ϵ_rb = sqrt(eps()), ϵ_rc = sqrt(eps())), display = false)
  x, r = statsnls.solution, statsnls.solver_specific[:r]
  @test norm(x - A \ b) ≤ norm(b) * sqrt(eps())
  @test norm(A * x - b - r) ≤ norm(b) * sqrt(eps())

  # least-square problem with constraints
  A = rand(m, n)
  b = rand(m)
  lcon, ucon = zeros(m), fill(Inf, m)
  C = ones(m, n)
  lvar, uvar = fill(-10.0, n), fill(200.0, n)
  nls = LLSModel(A, b, lvar = lvar, uvar = uvar, C = C, lcon = lcon, ucon = ucon)
  statsnls = ripqp(nls, display = false)
  x, r = statsnls.solution, statsnls.solver_specific[:r]
  @test length(r) == m
  @test norm(A * x - b - r) ≤ norm(b) * sqrt(eps())
end

@testset "matrixwrite" begin
  stats1 = ripqp(
    QuadraticModel(qps1),
    display = false,
    iconf = InputConfig(w = SystemWrite(write = true, name = "test_", kfirst = 4, kgap = 1000)),
  )
  @test isapprox(stats1.objective, -1.59078179, atol = 1e-2)
  @test stats1.status == :acceptable

  K = MatrixMarket.mmread("test_K_iter4.mtx")
  rhs_aff = readdlm("test_rhs_iter4_aff.rhs", Float64)[:]
  rhs_cc = readdlm("test_rhs_iter4_cc.rhs", Float64)[:]
  @test size(K, 1) == size(K, 2) == length(rhs_aff) == length(rhs_cc)
end

function QuadraticModelMaximize(qps, x0 = zeros(qps.nvar))
  QuadraticModel(
    -qps.c,
    qps.qrows,
    qps.qcols,
    -qps.qvals,
    Arows = qps.arows,
    Acols = qps.acols,
    Avals = qps.avals,
    lcon = qps.lcon,
    ucon = qps.ucon,
    lvar = qps.lvar,
    uvar = qps.uvar,
    c0 = -qps.c0,
    x0 = x0,
    minimize = false,
  )
end

@testset "maximize" begin
  stats1 = ripqp(QuadraticModelMaximize(qps1), display = false)
  @test isapprox(stats1.objective, 1.59078179, atol = 1e-2)
  @test stats1.status == :acceptable

  stats2 = ripqp(QuadraticModelMaximize(qps2), display = false)
  @test isapprox(stats2.objective, 9.99599999e1, atol = 1e-2)
  @test stats2.status == :acceptable

  stats3 = ripqp(QuadraticModelMaximize(qps3), display = false)
  @test isapprox(stats3.objective, -5.32664756, atol = 1e-2)
  @test stats3.status == :acceptable
end

@testset "presolve" begin
  stats5 = ripqp(QuadraticModel(qps5), display = false)
  @test isapprox(stats5.objective, 0.250000001, atol = 1e-2)
  @test stats5.status == :acceptable

  stats5 = ripqp(QuadraticModel(qps5), iconf = InputConfig(sp = K2KrylovParams()), display = false)
  @test isapprox(stats5.objective, 0.250000001, atol = 1e-2)
  @test stats5.status == :acceptable
end

qp_linop = QuadraticModel(
  c,
  LinearOperator(Q),
  A = LinearOperator(A),
  lcon = b,
  ucon = b,
  lvar = l,
  uvar = u,
  c0 = 0.0,
  name = "QM_LINOP",
)
qp_dense = QuadraticModel(
  c,
  tril(Q),
  A = A,
  lcon = b,
  ucon = b,
  lvar = l,
  uvar = u,
  c0 = 0.0,
  name = "QM_LINOP",
)

@testset "Dense and LinearOperator QPs" begin
  stats_linop = ripqp(
    qp_linop,
    iconf = InputConfig(sp = K2KrylovParams(), presolve = false, scaling = false),
    display = false,
  )
  @test isapprox(stats_linop.objective, 1.1249999990782493, atol = 1e-2)
  @test stats_linop.status == :acceptable

  stats_dense = ripqp(qp_dense, iconf = InputConfig(sp = K2KrylovParams(uplo = :U)))
  @test isapprox(stats_dense.objective, 1.1249999990782493, atol = 1e-2)
  @test stats_dense.status == :acceptable

  for fact_alg in [:bunchkaufman, :ldl]
    stats_dense = ripqp(
      qp_dense,
      iconf = InputConfig(
        sp = K2LDLDenseParams(fact_alg = fact_alg, ρ0 = 0.0, δ0 = 0.0),
        presolve = false,
        scaling = false,
      ),
      display = false,
    )
    @test isapprox(stats_dense.objective, 1.1249999990782493, atol = 1e-2)
    @test stats_dense.status == :acceptable
  end

  # test conversion
  stats_dense = ripqp(qp_dense)
  @test isapprox(stats_dense.objective, 1.1249999990782493, atol = 1e-2)
  @test stats_dense.status == :acceptable
end


text:
module GraphGLRM

using LowRankModels
using LightGraphs
using DataArrays, DataFrames
import LowRankModels: prox, prox!,
      evaluate, ObsArray, sort_observations, observations,
      fit!
import Base.BLAS: axpy!, gemm!
import Base.Threads

export #impute_means, impute_zeros, standardize, standardize!, #Simple imputation for pre-preprocessing
      #matrixRegFact, quadgraphRegFact, init_qqreg!,#Closed-form factorizations with regularization
      IndexGraph, #The data structure for easily initializing a GraphQuadReg

      #The regularizer
      AbstractGraphReg, GraphQuadReg, NonNegGraphReg, matrix, prox, prox!, evaluate,

      #The constructors for the GGLRM itself
      GGLRM, add_offset!, equilibrate_variance!,

      #The alternating minimization and objective calculation
      fit!, whole_objective, loss_objective, fit_multithread!, fit_sparse!,
      reconstruct_obs!, reconstruct_obs

# package code goes here
include("indexgraph.jl")
include("graphquadreg.jl")
include("nonneggraphreg.jl")
include("regularizers.jl")
include("gglrm.jl")
include("offsetscale.jl")
#include("initialize.jl")
include("fit.jl")
include("fit_multithread.jl")
include("fit_sparse.jl")
#include("constantstepsize.jl")

end # module


text:
using BigRationals
using Test
using Aqua

Aqua.test_all(BigRationals)

ONE = BigRational(1,1)
ZERO = BigRational(0,1)

@testset "Construction" begin
    @test ZERO == zero(BigRational) == 0
    @test iszero(ZERO)
    @test ONE == one(BigRational) == 1
    @test isone(ONE)
    @test BigRational(2,2) == UInt(1)
    @test BigRational(1,-3) == BigRational(-1,3)
    @test BigRational(-4,-6) == BigRational(2//3) == 2//3

    @test BigRational(-4) == BigRational(Int8(-4)) == BigRational(Int128(-4))
    @test BigRational(UInt(6)) == BigRational(BigInt(6)) == BigRational(6)
    @test BigRational(UInt16(1))  == BigRational(true) == BigRational(ONE) == 0x1

    @test BigRational(typemin(Int)+1,UInt(2)) == BigRational(typemax(Int),-2) == -BigRational(typemin(Int)+1,-2)
    @test BigRational(-4,Int8(-4)) == BigRational(typemax(UInt32), typemax(UInt32)) == Int16(1)
    @test BigRational(Int16(3),Int8(3)) == BigRational(UInt(1),1) == BigRational(Int16(4),0x4) == ONE
    @test BigRational(-3,UInt(5)) == BigRational(UInt16(6),-10) == -BigRational(-3,-5)
    @test BigRational(UInt(14),-7) == BigRational(typemax(UInt)-1,-typemax(Int)) == BigRational(0x6,-3) == BigRational(-2,1)
    @test BigRational(0x0007,0x7) == BigRational(0x9,9) == BigRational(BigInt(1),true) == one(BigInt)
    @test BigRational(typemin(Int16),-2) == cld(typemax(Int16),2)
    @test BigRational(-3,typemin(Int16)) == BigRational(3,typemax(Int16)+1)
    @test BigRational(typemin(Int32), UInt(3)) == BigRational(typemin(Int32)%Int64, UInt(3)) == BigRational(typemin(Int32), UInt16(3))
    @test BigRational(typemin(Int64), UInt(3)) == BigRational(typemin(Int64)%Int64, UInt(3)) == BigRational(typemin(Int64), UInt16(3))
    @test BigRational(typemin(Int), typemin(Int)) == -BigRational(-3,3) == -BigRational(7,-7) == ONE
    @test BigRational(3,false) == BigRational(true,false) == BigRational(true, 0) == BigRational(Int128(1),Int128(0))
    @test BigRational(true, 12) == 1//12
    @test iszero(BigRational(false, 0x1))

    @test BigRational(1.5) == BigRational(3,2) == 1.5
    @test BigRational(Float16(ℯ)) == Float16(ℯ)
    @test BigRational(BigFloat(13.25)) == BigRational(53,4) == 13.25
end

@testset "Rational functions" begin
    unsigned_types_test = (UInt8, UInt32, UInt64)
    for i in 1:3, j in i:3
        y1 = unsigned_types_test[i](9); y2 = unsigned_types_test[j](6)
        z = BigRational(y1, y2)
        @test numerator(z) == 3
        @test denominator(z) == 2
        if i != j
            z = BigRational(y2, y1)
            @test numerator(z) == 2
            @test denominator(z) == 3

        end
    end
    signed_types_test = (Int8, Int32, Int64, BigInt)
    for i in 1:4, j in i:4
        x1 = signed_types_test[i](12); x2 = signed_types_test[j](-4)
        z = BigRational(x1, x2)
        @test numerator(z) == -3
        @test denominator(z) == 1
        if i != j
            z = BigRational(x2, x1)
            @test numerator(z) == -1
            @test denominator(z) == 3
            for k in 1:3, t in i:3
                y1 = unsigned_types_test[k](3); y2 = unsigned_types_test[t](2)
                z = BigRational(x1, y1)
                @test numerator(z) == 4
                @test denominator(z) == 1
                z = BigRational(y1,x1)
                @test numerator(z) == 1
                @test denominator(z) == 4
                z = BigRational(x2, y1)
                @test numerator(z) == -4
                @test denominator(z) == 3
                z = BigRational(y1,x2)
                @test numerator(z) == -3
                @test denominator(z) == 4
                if k != t
                    z = BigRational(x1, y2)
                    @test numerator(z) == 6
                    @test denominator(z) == 1
                    z = BigRational(y2,x1)
                    @test numerator(z) == 1
                    @test denominator(z) == 6
                    z = BigRational(x2, y2)
                    @test numerator(z) == -2
                    @test denominator(z) == 1
                    z = BigRational(y2,x2)
                    @test numerator(z) == -1
                    @test denominator(z) == 2
                end
            end
        end
    end

    @test Rational(BigRational(3,4)) == Rational{BigInt}(BigRational(3,4)) == 3//4
    @test Rational{Int16}(BigRational(typemax(Int32),typemax(Int32))) === one(Rational{Int16})
    @test isinteger(BigRational(5))
    @test !isinteger(BigRational(3,4))
end

@testset "Conversions and promotions" begin
    @test BigFloat(BigRational(3,4)) == float(BigRational(3,4)) == Float16(BigRational(3,4)) == 0.75
    @test Bool(one(BigRational)) == !Bool(zero(BigRational)) == Int8(one(BigRational)) == true
    @test_throws InexactError Bool(BigRational(3,4))
    @test_throws InexactError Int(BigRational(3,4))
    @test promote(3//4, BigRational(2,1)) == promote(BigRational(3,4), 2) == (BigRational(3,4), BigRational(2))

    x1 = promote(BigRational(1,2), 0.75)
    x2 = promote(Float16(0.5), BigRational(3,4))
    x3 = (BigFloat(0.5), BigFloat(0.75))
    @test (x1[1] ≈ x2[1] ≈ x3[1]) && (x1[2] ≈ x2[2] ≈ x3[2])
end

@testset "IO operations" begin
    @test string(BigRational(6,3)) == "BigRational(2,1)"
    @test string(BigRational[BigRational(-2,3), BigRational(4,5)]) == "BigRational[-2//3, 4//5]"
    io = IOBuffer()
    show(io, BigRational(2,-4))
    @test String(take!(io)) == "BigRational(-1,2)"

    io = IOBuffer()
    # The following tests are broken because there is no read/write for BigInt in Base
    @test_broken write(io, BigRational(2,3)) == 32
    @test_broken read(io, BigRational) == BigRational(2,3)
end

@testset "Comparisons" begin
    @test BigRational(3,4) == BigRational(6,8) == Rational{Int}(3,4) == Rational{UInt}(3,4)
    @test BigRational(2,3) != BigRational(4,7) != 2//5
    @test 0x3 == BigRational(3) == 3
    @test Int16(5) != BigRational(3) != UInt(7)
    @test big(-4) == BigRational(-4) == Int128(-4)
    @test Int128(12) != BigRational(-4) != big(4)
    @test zero(BigRational) == false
    @test iszero(zero(BigRational))
    @test !iszero(typemin(BigRational))
    @test true == one(BigRational)
    @test isone(one(BigRational))
    @test !isone(typemax(BigRational))
    @test Rational{Int128}(7,9) == BigRational(7,9) != Rational{Bool}(true, true)

    @test false < BigRational(3,4) < BigRational(4,2) < 3
    @test BigRational(2,3) >= 2//3 >= BigRational(2,3)
    @test Rational{UInt}(3,5) <= BigRational(3,5) <= Rational{BigInt}(3,5)
    @test BigInt(1) <= BigRational(2,2) < 0x2
    @test Int128(3) <= BigRational(23,7) < Rational{Int128}(25,7)
    @test Inf >= BigRational(4,5) > 0.3
    @test !(BigRational(3,5) > NaN) && !(BigRational(4,7) < NaN) && !(BigRational(5,6) == NaN)

    @test isfinite(BigRational(-4,3))
    @test !isfinite(BigRational(2,0))
    @test !isfinite(BigRational(-3,0))
end

@testset "Arithmetic operations" begin
    @test numerator(BigRational(typemax(Int),typemin(Int))) == -typemax(Int)
    @test numerator(-BigRational(typemax(Int),typemin(Int))) == typemax(Int)
    for op in (:+, :-, :*, :/, ://, :div, :rem, :mod, :fld, :cld)
        ret1 = @eval $op(47//7, 13//5)
        ret2 = @eval $op(15//2, 26//11)
        for T in (:(BigRational), :(Rational{Int16}), :(Rational{UInt}), :(Rational{Int128}), :(Rational{BigInt}))
            @test @eval $op(BigRational(47,7), $T(13,5)) == $ret1
            @test @eval $op($T(15,2), BigRational(26,11)) == $ret2
        end
    end
    @test rem(BigRational(100,3), 2) == BigRational(100,3) - 2 * div(BigRational(100,3), 2)
    @test rem(5, BigRational(2,7)) == 5 - BigRational(2,7) * div(5, BigRational(2,7))
    @test fld(BigRational(100,3), 2) == div(BigRational(100,3), 2) == cld(BigRational(100,3), 2) - 1
    @test cld(-5, BigRational(2,7)) == div(-5, BigRational(2,7)) == fld(-5, BigRational(2,7)) + 1
    @test 4^(BigRational(3,2)) == 2.0^(BigRational(3,1)) == 8.0
    @test im^(BigRational(3,2)) ≈ -sqrt(2)/2 + sqrt(2)*im/2
    @test (BigRational(2,3)*im)^true == (2im)//3 == (BigRational(-3,2)*im)^(-1)
    @test (BigRational(1,1)*im)^false == BigRational(2,3)^0
    @test BigRational(4,9)^3 == BigRational(64,729)^1 == BigRational(8,27)^2
    @test BigRational(16,25)^(-1) == BigRational(25,16) == BigRational(4,5)^(-2)
    x = 2 # used to avoid ^ to Val lowering
    @test BigRational(16,25)^(x-3) == BigRational(25,16) == BigRational(4,5)^(-x)
    @test BigRational(2,3)^8 == BigRational(4,9)^4 == BigRational(16,81)^2
    @test ℯ^BigRational(3,4) ≈ 2.117
    for op in (:-, :inv)
        @test @eval $op(BigRational(12,7)) == $op(12//7)
        @test @eval $op(BigRational(-2,3)) == $op(-2//3)
    end
    for op in (:trunc, :floor, :ceil, :round)
        @test @eval $op(UInt, BigRational(14,5)) == $op(UInt, 14//5)
        for T in (:(Rational{Int8}), :Int, :(Union{BigInt,Missing}), :Float16)
            @test @eval $op($T, BigRational(45,13)) == $op($T, 45//13)
            @test @eval $op($T, BigRational(-7,9)) == $op($T, -7//9)
        end
    end
    @test round(Float64, typemax(BigRational)) === -round(Float64, typemin(BigRational)) === Inf
    @test round(BigRational(5,2), RoundUp) == BigRational(3)
    @test round(Union{Int,Missing}, BigRational(3,2), RoundDown) == round(Union{Int,Missing}, BigRational(3,2), RoundToZero)
    @test round(Union{Float64,Missing}, BigRational(1,4), RoundUp) == 1.0
    @test_throws DivideError round(Int, typemax(BigRational))
    if VERSION >= v"1.4.0-rc1"
        @test gcd(BigRational(3,4), BigRational(5,6)) == gcd(3//4, 5//6)
        @test lcm(BigRational(-2,7), BigRational(3,8)) == lcm(-2//7, 3//8)
        @test gcdx(BigRational(1,9), BigRational(5,-2)) == gcdx(1//9, -5//2)
        @test gcdx(BigRational(2,3), zero(BigRational)) == gcdx(2//3, 0//1)
        @test gcdx(BigRational(1,0), BigRational(5,3)) == BigRational.(gcdx(1//0, 5//3))
    end
    @test sign(BigRational(-3,4)) == -sign(BigRational(2,9))
    @test signbit(BigRational(5,-8)) && !signbit(BigRational(3,4))
    @test copysign(BigRational(2,3), -4.3) == copysign(BigRational(-4,6), BigRational(-1,3)) == -2//3
    @test copysign(BigRational(-4,5), BigInt(2)) == copysign(BigRational(4,5), 2//3) == BigRational(4,5)
    @test flipsign(BigRational(1,3), -3.0)  == flipsign(BigRational(2,6), -1) == flipsign(BigRational(-1,3), true) == -1//3
    @test flipsign(BigRational(2,5), 0x2) == copysign(BigRational(-2,5), 0x0003) == BigRational(2,5)
    @test abs(BigRational(typemin(Int)+2,typemax(Int))) == abs(BigRational(typemax(Int)-1,typemax(Int)))
end

@testset "Other functions" begin
    x = BigRational(3,4)
    @test widen(BigRational) == BigRational
    @test widen(x) isa BigRational
    @test big(BigRational) == BigRational
    @test big(x) === x
    @test Base.hastypemax(BigRational)
    @test denominator(typemax(BigRational)) == denominator(typemin(BigRational)) == 0
    @test numerator(typemax(BigRational)) == -numerator(typemin(BigRational)) == 1
    @test fma(x, x, x) == x*(x+1)
    @test fma(BigRational(2,3), BigRational(6,5), BigRational(12,7)) == 2//3 * 6//5 + 12//7
    y = deepcopy(x)
    BigRationals.MPQ.mul_2exp!(y, 3)
    @test x == 3//4
    @test y == x * (1<<3)
    @test hash(x) == hash(3//4)
    @test hash(y, UInt(123)) == hash(x*(1<<3), UInt(123))
    z = [BigRational(7,5)]
    push!(z, BigRationals.MPQ.div_2exp(first(z), 2))
    push!(z, first(z))
    t = deepcopy(z)
    @test first(t) === last(t)
    @test first(t) != t[2]
    @test hash(t) == hash(z)
    @test t == z == BigRational[7//5, 7//20, 7//5]
    BigRationals.MPQ.swap!(z[1], z[3])
    @test hash(t) == hash(z)
    @test t == z == BigRational[7//5, 7//20, 7//5]
    @test first(t) === last(t)
    BigRationals.MPQ.swap!(z[1], z[2])
    @test hash(t) != hash(z)
    @test t != z
    @test z == BigRational[7//20, 7//5, 7//20]
    @test first(t) === last(t)
end


text:
using Printf

struct Integrator
    step
    dt
    T
    const_dt
    kwargs
end
Integrator(step::Function,dt,T,const_dt=false;kwargs...) =
    Integrator(step,dt,T,const_dt,kwargs)

function integrate(HS,step::Function,dt,T;kwargs...)
    t=[.0]
    Q=[HS.q0]
    P=[HS.p0]
    N = Int(ceil(T/dt/100))
    progress = Progress(N, desc="Integrating using $(typeof(step).name.mt.name) dt = $dt...")
    lastupdate=0
    while(t[end]<T)
        DT,q,p=step(Q[end],P[end],dt,HS;kwargs...)
        push!(Q,q)
        push!(P,p)
        push!(t,t[end]+DT)
        if t[end]-lastupdate>T/N
            lastupdate=t[end]
            ProgressMeter.next!(progress)
        end
    end
    ProgressMeter.finish!(progress)
    return (t,hcat(Q...),hcat(P...))
end

function integrate_constant_dt(HS,step::Function,dt,T;kwargs...)
    N = Int(ceil(T/dt))
    Q=zeros(length(HS.q0),N+1)
    P=zeros(length(HS.q0),N+1)
    t=[dt*i for i in 0:N]
    Q[:,1]=HS.q0
    P[:,1]=HS.p0
    progress = Progress(N, desc="Integrating using $(typeof(step).name.mt.name) dt = $dt...")
    lastupdate=0
    for i in 1:N
        DT,q,p=step(Q[:,i],P[:,i],dt,HS;kwargs...)
        Q[:,i+1]=q
        P[:,i+1]=p
        ProgressMeter.next!(progress)
    end
    ProgressMeter.finish!(progress)
    return (t,Q,P)
end

integrate(HS,I::Integrator) = I.const_dt ?
    integrate_constant_dt(HS,I.step,I.dt,I.T;I.kwargs...) :
    integrate(HS,I.step,I.dt,I.T;I.kwargs...)

function integrate(HS,Is)
    ts = []
    Qs = []
    Ps = []
    for I in Is
        t,Q,P = integrate(HS,I)
        push!(ts,t)
        push!(Qs,Q)
        push!(Ps,P)
    end
    return (ts,Qs,Ps)
end

function getNames(Is)
     return ["$(typeof(I.step).name.mt.name) with dt=$(@sprintf("%.2e", I.dt))$( length(I.kwargs)>0 ? "\n$(join(["$k: $(I.kwargs[k])" for k in keys(I.kwargs)],", "))" : "")" for I in Is]
end

println("Finished loading HamiltonIntegrator")

max_stars_repo_name: BoergeSt/HamiltonODE

max_issues_repo_name: BoergeSt/HamiltonODE

max_forks_repo_name: BoergeSt/HamiltonODE


text:
__precompile__()

module IntervalTrees

export
    IntervalTree,
    IntervalMap,
    AbstractInterval,
    Interval,
    IntervalValue,
    depth,
    hasintersection,
    from,
    value

include("slice.jl")

"""
An `AbstractInterval{T}` must have a `first` and `last` function each returning
a value of type `T`, and `first(i) <= last(i)` must always be true.
"""
abstract AbstractInterval{T}

function Base.isless(u::AbstractInterval, v::AbstractInterval)
    return first(u) < first(v) || (first(u) == first(v) && last(u) < last(v))
end

"""
A basic interval.
"""
immutable Interval{T} <: AbstractInterval{T}
    first::T
    last::T
end

Base.first{T}(i::Interval{T}) = i.first
Base.last{T}(i::Interval{T}) = i.last

"""
An interval with some associated data.
"""
immutable IntervalValue{K, V} <: AbstractInterval{K}
    first::K
    last::K
    value::V
end

Base.first{K, V}(i::IntervalValue{K, V}) = i.first
Base.last{K, V}(i::IntervalValue{K, V}) = i.last
value{K, V}(i::IntervalValue{K, V}) = i.value

Base.print(io::IO, x::Interval) = print(io, "\n($(first(x)),$(last(x)))")
function Base.show(io::IO, x::Interval)
    show(io, typeof(x))
    print(io, x)
end

Base.print(io::IO, x::IntervalValue) = print(io, "\n($(first(x)),$(last(x))) => $(value(x))")
function Base.show(io::IO, x::IntervalValue)
    show(io, typeof(x))
    print(io, x)
end

# Each of these types is indexes by K, V, B, where
#   K : Interval type. Intervals are represented as (K, K) tuples.
#   V : Interval value type. This an interval type that may have associated
#        data. It must have `first`, `last`, and `isless` methods.
#   B : Integer giving the B-tree order.

abstract Node{K, V, B}

type InternalNode{K, V, B} <: Node{K, V, B}
    # Internal nodes are keyed by the minimum interval in the right subtree.  We
    # need internal node keys to be intervals themselves, since ordering by
    # start value alone only works if start values are unique. We don't
    # force that restriction, so we must break ties in order to split nodes that
    # are full of intervals with the same start value.
    keys::Slice{Interval{K}, B}

    # Maximum interval end-point in this subtree.
    maxend::K

    children::Slice{Node{K, V, B}, B}

    parent::Nullable{InternalNode{K, V, B}}

    # Sibling/cousin pointers.
    left::Nullable{InternalNode{K, V, B}}
    right::Nullable{InternalNode{K, V, B}}

    function InternalNode()
        t = new(Slice{Interval{K}, B}(), zero(K), Slice{Node{K, V, B}, B}(),
                Nullable{InternalNode{K, V, B}}(),
                Nullable{InternalNode{K, V, B}}(),
                Nullable{InternalNode{K, V, B}}())
        return t
    end
end

function minstart(t::InternalNode)
    return t.keys[1].first
end

type LeafNode{K, V, B} <: Node{K, V, B}
    entries::Slice{V, B}

    # maximum interval end-point in this leaf node.
    maxend::K

    parent::Nullable{InternalNode{K, V, B}}

    # Sibling/cousin pointers.
    left::Nullable{LeafNode{K, V, B}}
    right::Nullable{LeafNode{K, V, B}}

    function LeafNode()
        return new(
            Slice{V, B}(),
            zero(K),
            Nullable{InternalNode{K, V, B}}(),
            Nullable{LeafNode{K, V, B}}(),
            Nullable{LeafNode{K, V, B}}())
    end
end

function minstart(t::LeafNode)
    return first(t.entries[1])
end

type IntervalBTree{K, V, B}
    root::Node{K, V, B}
    n::Int # Number of entries

    function IntervalBTree()
        return new(LeafNode{K, V, B}(), 0)
    end

    # Construct an interval tree from a sorted array of intervals.
    #
    # This is generally much more efficient than constructing the same tree by
    # inserting intervals one by one.
    #
    # Args:
    #   entries: Interval entry values in sorted order.
    #
    function IntervalBTree(entries::AbstractVector{V})
        if !issorted(entries)
            error("Intervals must be sorted to construct an IntervalTree")
        end

        # Here's the plan: figure out how many leaf nodes we need, allocate
        # them all up front. Copy in the keys and values, then work up towards
        # the root.

        n = length(entries)

        if n == 0
            return new(LeafNode{K, V, B}(), 0)
        end

        numleaves = cld(n, B - 2)
        leaves = [LeafNode{K, V, B}() for _ in 1:numleaves]

        maxends = Array(K, numleaves)
        minkeys = Array(Interval{K}, numleaves)

        # divy up the keys and values among the leaves
        keys_per_leaf = cld(n, numleaves)
        for i in 1:numleaves
            u = (i - 1) * keys_per_leaf + 1
            v = min(n, i * keys_per_leaf)
            minkeys[i] = Interval{K}(first(entries[u]), last(entries[u]))
            maxends[i] = last(entries[u])
            for j in u:v
                push!(leaves[i].entries, entries[j])
                maxends[i] = max(maxends[i], last(entries[j]))
            end
            leaves[i].maxend = maxends[i]
        end

        children = leaves
        while length(children) > 1
            # sibling pointers
            for i in 1:length(children)-1
                children[i].right = children[i+1]
                children[i+1].left = children[i]
            end

            # make parents
            numparents = cld(length(children), B - 2)
            parents = [InternalNode{K, V, B}() for _ in 1:numparents]

            # divy up children among parents
            children_per_parent = cld(length(children), numparents)
            for i in 1:numparents
                u = (i - 1) * keys_per_leaf + 1
                v = min(length(children), i * keys_per_leaf)
                maxend = maxends[u]
                for j in u:v
                    push!(parents[i].children, children[j])
                    children[j].parent = parents[i]
                    maxend = max(maxend, maxends[j])
                    if j > u
                        push!(parents[i].keys, minkeys[j])
                    end
                end
                minkeys[i] = minkeys[u]
                parents[i].maxend = maxend
                maxends[i] = maxend
            end

            children = parents
        end
        @assert length(children) == 1

        return new(children[1], n)
    end
end

# Default B-tree order
typealias IntervalTree{K, V} IntervalBTree{K, V, 64}

# Show

function Base.show(io::IO, it::IntervalTree)
    show(io, typeof(it))
    n = length(it)
    if length(it) > 6
        # Hacky random access ...
        for (i, x) in enumerate(it)
            i < 4 && print(x)
        end
        print("\n\u22EE") # Vertical ellipsis
        for (i, x) in enumerate(it)
            i > (n-3) && print(x)
        end
    else
        for x in it
            print(x)
        end
    end
end

# Length
# ------

function Base.length(t::IntervalBTree)
    return t.n
end

function depth(t::IntervalBTree)
    k = 1
    node = t.root
    while !isa(node, LeafNode)
        k += 1
        node = node.children[1]
    end
    return k
end

function Base.length(t::InternalNode)
    return length(t.children)
end

function Base.length(t::LeafNode)
    return length(t.entries)
end

function Base.isempty(t::IntervalBTree)
    return t.n == 0
end

function Base.isempty(t::InternalNode)
    return isempty(t.children)
end

function Base.isempty(t::LeafNode)
    return isempty(t.entries)
end

# Iterating
# ---------

immutable IntervalBTreeIteratorState{K, V, B}
    leaf::Nullable{LeafNode{K, V, B}}
    i::Int
end

function Base.start{K, V, B}(t::IntervalBTree{K, V, B})
    # traverse to the first leaf node
    node = t.root
    while !isa(node, LeafNode{K, V, B})
        node = node.children[1]
    end
    return IntervalBTreeIteratorState(Nullable(node), 1)
end

function Base.next{K, V, B}(t::IntervalBTree{K, V, B},
                            state::IntervalBTreeIteratorState{K, V, B})
    leaf = get(state.leaf)
    entry = leaf.entries[state.i]
    if state.i < length(leaf)
        state = IntervalBTreeIteratorState{K, V, B}(leaf, state.i + 1)
    else
        state = IntervalBTreeIteratorState{K, V, B}(leaf.right, 1)
    end
    return entry, state
end

function Base.done{K, V, B}(t::IntervalBTree{K, V, B},
                            state::IntervalBTreeIteratorState{K, V, B})
    return isnull(state.leaf) || isempty(get(state.leaf))
end

# Iterate from a given starting from the first interval that intersects a given
# point.

immutable IntervalFromIterator{K, V, B}
    t::IntervalBTree{K, V, B}
    p::K
end

function from{K, V, B}(t::IntervalBTree{K, V, B}, p)
    return IntervalFromIterator{K, V, B}(t, convert(K, p))
end

function Base.start{K, V, B}(it::IntervalFromIterator{K, V, B})
    node, i = firstfrom(it.t, it.p)
    if i == 0
        return IntervalBTreeIteratorState{K, V, B}(Nullable{LeafNode{K, V, B}}(), i)
    else
        return IntervalBTreeIteratorState{K, V, B}(node, i)
    end
end

function Base.next{K, V, B}(it::IntervalFromIterator{K, V, B},
                            state::IntervalBTreeIteratorState{K, V, B})
    return next(it.t, state)
end

function Base.done{K, V, B}(it::IntervalFromIterator{K, V, B},
                            state::IntervalBTreeIteratorState{K, V, B})
    return done(it.t, state)
end

# Inserting
# ---------

# Split a leaf into two, returning (leftnode, rightnode)
function split!{K, V, B}(left::LeafNode{K, V, B})
    right = LeafNode{K, V, B}()
    right.right = left.right
    right.left = left
    left.right = right
    if !isnull(right.right)
        get(right.right).left = right
    end
    right.parent = left.parent

    m = length(left)
    resize!(right.entries, m - div(m, 2))
    copy!(right.entries, 1, left.entries, div(m, 2) + 1, length(right.entries))
    resize!(left.entries, div(m, 2))

    left.maxend = last(left.entries[1])
    for entry in left.entries
        left.maxend = max(left.maxend, last(entry))
    end

    right.maxend = last(right.entries[1])
    for entry in right.entries
        right.maxend = max(right.maxend, last(entry))
    end

    return left, right
end

# Split an internal node in two, returning (leftnode, rightnode)
function split!{K, V, B}(left::InternalNode{K, V, B})
    right = InternalNode{K, V, B}()
    right.right = left.right
    right.left = left
    left.right = right
    if !isnull(right.right)
        get(right.right).left = right
    end
    right.parent = left.parent

    m = length(left)

    resize!(right.children, m - div(m, 2))
    resize!(right.keys, m - div(m, 2) - 1)

    copy!(right.children, 1, left.children, div(m, 2)+1, length(right.children))
    copy!(right.keys, 1, left.keys, div(m, 2)+1, length(right.keys))

    resize!(left.children, div(m, 2))
    resize!(left.keys, div(m, 2) - 1)

    left.maxend = left.children[1].maxend
    for child in left.children
        left.maxend = max(left.maxend, child.maxend)
    end

    right.maxend = right.children[1].maxend
    for child in right.children
        right.maxend = max(right.maxend, child.maxend)
    end

    for child in right.children
        child.parent = right
    end

    return left, right
end

# Find the maximum interval end point is a subtree
function nodemaxend{K, V, B}(t::InternalNode{K, V, B})
    maxend = zero(K)
    for i in 1:length(t.children)
        maxend = max(maxend, t.children[i].maxend)
    end
    return maxend
end

function nodemaxend{K, V, B}(t::LeafNode{K, V, B})
    maxend = zero(K)
    for i in 1:length(t.entries)
        maxend = max(maxend, last(t.entries[i]))
    end
    return maxend
end

function sameparent{K, V, B}(u::Node{K, V, B}, v::Node{K, V, B})
    return (isnull(u.parent) && isnull(v.parent)) ||
           (!isnull(u.parent) && !isnull(v.parent) && get(u.parent) == get(v.parent))
end

# Find the first leaf node in the tree
function firstleaf(t::IntervalBTree)
    return firstleaf(t.root)
end

function firstleaf(t::InternalNode)
    return firstleaf(t.children[1])
end

function firstleaf(t::LeafNode)
    return t
end

# TODO: Having to construct Nullables here is killng us.
# How can we avoid that?
#
function nextleafkey{K, V, B}(t::LeafNode{K, V, B}, i::Int)
    if i < length(t)
        return (Nullable{LeafNode{K, V, B}}(t), i + 1)
    else
        return (t.right, 1)
    end
end

# This version of nextleafkey is used in nextintersection at aggressively avoid
# allocation, and dereferecing nullables more than once.
macro nextleafkey(leaf, nleaf, i)
    quote
        if $(esc(i)) < length($(esc(leaf)))
            $(esc(i)) += 1
        else
            $(esc(i)) = 1
            $(esc(nleaf)) = $(esc(leaf)).right
        end
    end
end

# Find the minimum interval in a subtree
function minkey(t::InternalNode)
    return minkey(t.children[1])
end

function minkey{K, V, B}(t::LeafNode{K, V, B})
    return Interval{K}(first(t.entries[1]), last(t.entries[1]))
end

function Base.push!{K, V, B}(t::IntervalBTree{K, V, B}, entry::V,
                             unique_key::Bool=false, update::Bool=true)
    return _push!(t, t.root, entry, unique_key, update)
end

function _push!{K, V, B}(t::IntervalBTree{K, V, B},
                         node::InternalNode{K, V, B},
                         entry::V, unique_key::Bool, update::Bool)
    i = findidx(node, entry) # key index
    j = i <= length(node) - 1 && entry >= node.keys[i] ? i + 1 : i # child index
    return _push!(t, node.children[j], entry, unique_key, update)
end

function _push!{K, V, B}(t::IntervalBTree{K, V, B},
                         node::LeafNode{K, V, B}, entry::V,
                         unique_key::Bool, update::Bool)
    i = max(1, findidx(node, entry))

    # key exists and we are replacing it
    if i <= length(node) && unique_key &&
        first(node.entries[i]) == first(entry) && last(node.entries[i]) == last(entry)
        if update
            node.entries[i] = entry
            return entry
        else
            return node.entries[i]
        end
    end

    #unsafe_insert!(node.entries, i, entry)
    insert!(node.entries, i, entry)
    if length(node) == 1 || last(entry) > node.maxend
        node.maxend = last(entry)
    end
    t.n += 1

    # split when full
    if length(node) == B
        maxend = node.maxend
        leftleaf, rightleaf = split!(node)
        median = Interval{K}(first(rightleaf.entries[1]),
                             last(rightleaf.entries[1]))

        # travel back up the tree setting maxend values and splitting as needed
        parent = node.parent
        leftnode = leftleaf
        rightnode = rightleaf
        hassplit = true
        while !isnull(parent)
            p = get(parent)
            if hassplit
                i = findidx(p, entry) # key index
                j = i <= length(p) - 1 && entry >= p.keys[i] ? i + 1 : i # child index
                insert!(p.children, j + 1, rightnode)
                p.maxend = max(p.maxend, maxend)
                insert!(p.keys, j, median)

                # split when full
                if length(p) == B
                    leftnode, rightnode = split!(p)
                    median = minkey(rightnode)
                    maxend = max(leftnode.maxend, rightnode.maxend)
                else
                    hassplit = false
                end
            else
                p.maxend = max(p.maxend, last(entry))
            end
            parent = p.parent
        end

        if hassplit
            t.root = InternalNode{K, V, B}()
            push!(t.root.keys, median)
            push!(t.root.children, leftnode)
            push!(t.root.children, rightnode)
            t.root.maxend = maxend
            leftnode.parent = t.root
            rightnode.parent = t.root
        end
    else
        # travel back up the tree setting maxend values
        parent = node.parent
        while !isnull(parent)
            p = get(parent)
            p.maxend = max(p.maxend, last(entry))
            parent = p.parent
        end
    end

    return entry
end

# Deleting
# --------

function deletefirst!{K, V, B}(t::IntervalBTree{K, V, B}, first::K, last::K)
    return deletefirst!(t, Interval{K}(first, last))
end

function deletefirst!{K, V, B}(t::IntervalBTree{K, V, B}, key::Tuple{K, K})
    return deletefirst!(t, Interval{K}(key[1], key[2]))
end

function deletefirst!{K, V, B}(t::IntervalBTree{K, V, B}, key::Interval{K})
    result = _deletefirst!(t.root, key)
    if result.keyfound
        t.n -= 1
    end

    # if the root has only one child, promote the child
    if isa(t.root, InternalNode) && length(t.root) == 1
        t.root = t.root.children[1]
        t.root.parent = Nullable{InternalNode{K, V, B}}()
    end

    return t
end

# Indicate what steps are need to account for an updated child.
immutable KeyFate
    value::UInt8
end

const KEYFATE_NONE         = KeyFate(0) # no changes
const KEYFATE_DELETE       = KeyFate(1) # delete the child node
const KEYFATE_DELETE_RIGHT = KeyFate(2) # delete the child's right sibling
const KEYFATE_UPDATE_LEFT  = KeyFate(3) # update the key separating the node
                                        # from its left sibling
const KEYFATE_UPDATE_RIGHT = KeyFate(4) # update the key separating the node
                                        # from its right sibling

# Information returned from a call to _delete
immutable DeleteResult
    keyfound::Bool # true if the key was found

    # Indicate what steps are need to account for an updated child. One of:
    fate::KeyFate
end

# Delete key from the subtree is present. Return a DeleteResult object.
function _deletefirst!{K, V, B}(t::InternalNode{K, V, B}, key::Interval{K})
    i = findidx(t, key)
    j = i <= length(t) - 1 && key >= t.keys[i] ? i + 1 : i # child index

    ans = _deletefirst!(t.children[j], key)

    if !ans.keyfound
        return DeleteResult(false, KEYFATE_NONE)
    end

    if ans.fate == KEYFATE_UPDATE_LEFT
        t.maxend = nodemaxend(t)
        if j > 1
            t.keys[j - 1] = minkey(t.children[j]);
            return DeleteResult(true, KEYFATE_NONE)
        else
            return DeleteResult(true, KEYFATE_UPDATE_LEFT)
        end
    elseif ans.fate == KEYFATE_UPDATE_RIGHT
        t.maxend = nodemaxend(t)
        if j < length(t)
            t.keys[j] = minkey(t.children[j+1])
        end

        if j > 1
            t.keys[j - 1] = minkey(t.children[j]);
            return DeleteResult(true, KEYFATE_NONE)
        else
            return DeleteResult(true, KEYFATE_UPDATE_LEFT)
        end
    elseif ans.fate == KEYFATE_DELETE || ans.fate == KEYFATE_DELETE_RIGHT
        deleteidx = ans.fate == KEYFATE_DELETE ? j : j + 1

        # deleteidx == 1 would only happen if we had only one child, but if
        # that were the case it would have been promoted.
        @assert deleteidx > 1

        splice!(t.keys, deleteidx - 1)
        splice!(t.children, deleteidx)
        t.maxend = nodemaxend(t)

        # not underfull
        minsize = div(B, 2)
        if length(t) >= minsize
            return DeleteResult(true, deleteidx == 2 ? KEYFATE_UPDATE_LEFT : KEYFATE_NONE)
        end

        # borrow right
        if !isnull(t.right)
            right = get(t.right)
            if sameparent(right, t) && length(right) > minsize
                splice!(right.keys, 1)
                push!(t.children, splice!(right.children, 1))
                push!(t.keys, minkey(t.children[end]))
                t.maxend = max(t.maxend, nodemaxend(t.children[end]))
                right.maxend = nodemaxend(right)
                t.children[end].parent = t

                return DeleteResult(true, KEYFATE_UPDATE_RIGHT)
            end
        end

        # borrow left
        if !isnull(t.left)
            left = get(t.left)
            if sameparent(left, t) && length(left) > minsize
                insert!(t.children, 1, pop!(left.children))
                pop!(left.keys)
                insert!(t.keys, 1, minkey(t.children[2]))
                t.maxend = max(t.maxend, nodemaxend(t.children[1]))
                left.maxend = nodemaxend(left)
                t.children[1].parent = t

                return DeleteResult(true, KEYFATE_UPDATE_LEFT)
            end
        end

        # merge with left
        if !isnull(t.left)
            left = get(t.left)
            if sameparent(left, t)
                merge!(left, t)
                left.right = t.right
                if !isnull(t.right)
                    get(t.right).left = left
                end

                return DeleteResult(true, KEYFATE_DELETE)
            end
        end

        # merge with right
        if !isnull(t.right)
            right = get(t.right)
            if sameparent(right, t)
                merge!(t, right)
                if !isnull(right.right)
                    get(right.right).left = t
                end
                t.right = right.right

                return DeleteResult(true, KEYFATE_DELETE_RIGHT)
            end
        end

        # Allow the root node to be underfull
        return DeleteResult(true, KEYFATE_NONE)
    else
        t.maxend = nodemaxend(t)
        return DeleteResult(true, KEYFATE_NONE)
    end
end

function _deletefirst!{K, V, B}(t::LeafNode{K, V, B}, key::Interval{K})
    i = findidx(t, key)

    # do nothing if the key isn't present
    if i < 1 || i > length(t) ||
        first(t.entries[i]) != key.first || last(t.entries[i]) != key.last
        return DeleteResult(false, KEYFATE_NONE)
    end

    splice!(t.entries, i)

    # This is the root node. Allow it to be empty.
    if isempty(t)
        return DeleteResult(true, KEYFATE_NONE)
    end

    minsize = div(B, 2)

    t.maxend = nodemaxend(t)

    # not underfull
    if length(t) >= minsize
        return DeleteResult(true, i == 1 ? KEYFATE_UPDATE_LEFT : KEYFATE_NONE)
    end

    # borrow right
    if !isnull(t.right)
        right = get(t.right)
        if sameparent(right, t) && length(right) > minsize
            push!(t.entries, splice!(right.entries, 1))
            t.maxend = max(t.maxend, last(t.entries[end]))
            right.maxend = nodemaxend(right)

            return DeleteResult(true, KEYFATE_UPDATE_RIGHT)
        end
    end

    # borrow left
    if !isnull(t.left)
        left = get(t.left)
        if sameparent(left, t) && length(left) > minsize
            insert!(t.entries, 1, pop!(left.entries))
            t.maxend = max(t.maxend, last(t.entries[1]))
            left.maxend = nodemaxend(left)

            return DeleteResult(true, KEYFATE_UPDATE_LEFT)
        end
    end

    # merge with left
    if !isnull(t.left)
        left = get(t.left)
        if sameparent(left, t)
            merge!(left, t)
            left.right = t.right
            if !isnull(t.right)
                get(t.right).left = left
            end
            return DeleteResult(true, KEYFATE_DELETE)
        end
    end

    # merge with right
    if !isnull(t.right)
        right = get(t.right)
        if sameparent(right, t)
            merge!(t, right)
            if !isnull(right.right)
                get(right.right).left = t
            end
            t.right = right.right
            return DeleteResult(true, KEYFATE_DELETE_RIGHT)
        end
    end

    # This must be the root node. Allow it to be underfull.
    return DeleteResult(true, KEYFATE_NONE)
end

# Join two leaf nodes into one.
function merge!{K, V, B}(left::LeafNode{K, V, B}, right::LeafNode{K, V, B})
    leftlen, rightlen = length(left), length(right)
    @assert leftlen + rightlen <= B
    resize!(left.entries, leftlen + rightlen)
    copy!(left.entries, leftlen+1, right.entries, 1, length(right.entries))
    left.maxend = nodemaxend(left)
    return
end

# Join two internal nodes into one
function merge!{K, V, B}(left::InternalNode{K, V, B}, right::InternalNode{K, V, B})
    leftlen, rightlen = length(left), length(right)
    @assert length(left) + length(right) <= B
    resize!(left.keys, leftlen + rightlen - 1)
    resize!(left.children, leftlen + rightlen)
    copy!(left.children, leftlen+1, right.children, 1, length(right.children))
    for i in leftlen+1:length(left.children)
        left.children[i].parent = left
    end
    left.keys[leftlen] = minkey(left.children[leftlen+1])
    copy!(left.keys, leftlen+1, right.keys, 1, length(right.keys))
    left.maxend = nodemaxend(left)
    return
end

# Searching
# ---------

function findidx{K, V, B}(t::LeafNode{K, V, B}, key::AbstractInterval{K})
    if isempty(t.entries)
        return 0
    end
    return searchsortedfirst(t.entries, key)
end

function findidx{K, V, B}(t::InternalNode{K, V, B}, key::AbstractInterval{K})
    if isempty(t.keys)
        return 0
    end
    i = searchsortedfirst(t.keys, key)
    return min(i, length(t.keys))
end

function Base.haskey{K, V, B}(t::LeafNode{K, V, B}, key::AbstractInterval{K})
    i = findidx(t, key)
    return 1 <= i <= length(t) && first(t.entries[i]) == key.first &&
           last(t.entries[i]) == key.last
end

function Base.haskey{K, V, B}(t::InternalNode{K, V, B}, key::AbstractInterval{K})
    i = findidx(t, key)
    if i <= length(t) - 1 && key >= t.keys[i]
        return haskey(t.children[i+1], key)
    else
        return haskey(t.children[i], key)
    end
end

function Base.haskey{K, V, B}(t::IntervalBTree{K, V, B}, key0::Tuple{Any, Any})
    key = Interval{K}(key0[1], key0[2])
    return haskey(t.root, key)
end

# Intersection
# ------------

# There are two ways we can tackle intersection: traversal (search the tree for
# intervals that intersect) and iteration (iterate through the tree's intervals
# in sorted order).
#
# Traversal is roughly O(log(n)) while iteration is roughly O(n). Obviously
# if we are doing a single intersection, we would choose traversal, but when
# intersecting two trees it becomes more interesting.
#
# If the second tree is of size m, then intersection by traversal is
# O(min(m * log(n), n * log(m))) while interesction by iteration is O(m + n).
# When m and n are both large, it becomes more efficient to intersect by
# iteration, but deciding when requires some heuristics.
#

# Return true iff two key1 and key2 intersect.
@inline function intersects{K}(key1::AbstractInterval{K}, key2::AbstractInterval{K})
    return first(key1) <= last(key2) && first(key2) <= last(key1)
end

# Return true if the tree has an intersection with the given position
function hasintersection{K, V, B}(t::IntervalBTree{K, V, B}, query)
    return hasintersection(t.root, convert(K, query))
end

function hasintersection{K, V, B}(t::InternalNode{K, V, B}, query::K)
    if isempty(t) || t.maxend < query
        return false
    end

    for (i, child) in enumerate(t.children)
        if child.maxend >= query && (i == 1 || t.keys[i-1].first <= query)
            if hasintersection(child, query)
                return true
            end
        elseif minstart(child) > query
            break
        end
    end

    return false
end

function hasintersection{K, V, B}(t::LeafNode{K, V, B}, query::K)
    if isempty(t) || t.maxend < query
        return false
    end

    for i in 1:length(t)
        if first(t.entries[i]) <= query <= last(t.entries[i])
            return true
        elseif query < first(t.entries[i])
            break
        end
    end

    return false
end

# Represent an intersection in an IntervalTree by pointing to a LeafNode
# and an index within that leaf node. No intersection is represented with
# index == 0.
type Intersection{K, V, B}
    index::Int
    node::LeafNode{K, V, B}

    Intersection(index, node) = new(index, node)
    Intersection() = new(0)
end

# Find the first interval in the tree that intersects the query and return
# as a (leafnode, index) pair, indicating that leafnode.keys[index] intersects.
# If no intersection is found, index is 0 and leafnode is the last node
# searched.
function firstintersection!{K, V, B}(t::IntervalBTree{K, V, B},
                                     query::AbstractInterval{K},
                                     out::Intersection{K, V, B})
    return firstintersection!(t.root, query, out)
end

function firstintersection!{K, V, B}(t::InternalNode{K, V, B},
                                     query::AbstractInterval{K},
                                     out::Intersection{K, V, B})
    if isempty(t) || t.maxend < first(query)
        out.index = 0
        return
    end

    for (i, child) in enumerate(t.children)
        if child.maxend >= first(query) && (i == 1 || t.keys[i-1].first <= last(query))
            firstintersection!(child, query, out)
            if out.index > 0
                return
            end
        elseif minstart(child) > last(query)
            break
        end
    end

    out.index = 0
    return
end

function firstintersection!{K, V, B}(t::LeafNode{K, V, B},
                                     query::AbstractInterval{K},
                                     out::Intersection{K, V, B})
    if isempty(t) || t.maxend < first(query)
        out.index = 0
        return
    end

    for i in 1:length(t)
        if intersects(t.entries[i], query)
            out.index = i
            out.node = t
            return
        elseif query.last < first(t.entries[i])
            break
        end
    end

    out.index = 0
    return
end

# Find the first key with an end point >= query
function firstfrom{K, V, B}(t::IntervalBTree{K, V, B}, query::K)
    return firstfrom(t.root, query)
end

function firstfrom{K, V, B}(t::InternalNode{K, V, B}, query::K)
    if isempty(t)
        return (t, 0)
    end

    for (i, child) in enumerate(t.children)
        if child.maxend >= query
            return firstfrom(child, query)
        end
    end

    return (t, 0)
end

function firstfrom{K, V, B}(t::LeafNode{K, V, B}, query::K)
    if isempty(t)
        return (t, 0)
    end

    for i in 1:length(t)
        if last(t.entries[i]) >= query
            return (t, i)
        end
    end

    return (t, 0)
end

# If query intersects node.keys[i], return the next intersecting key as
# a (leafnode, index) pair.
function nextintersection!{K, V, B}(t::LeafNode{K, V, B}, i::Integer,
                                    query::AbstractInterval{K},
                                    out::Intersection{K, V, B})
    j = i + 1
    while true
        while j <= length(t)
            if intersects(unsafe_getindex(t.entries, j), query)
                out.index = j
                out.node = t
                return
            end
            j += 1
        end
        j = 1
        if isnull(t.right)
            break
        end
        t = get(t.right)
        if minstart(t) > last(query)
            break
        end
    end

    out.index = 0
    return
end

type IntervalIntersectionIterator{K, V, B}
    intersection::Intersection{K, V, B}
    t::IntervalBTree{K, V, B}
    query::AbstractInterval{K}

    function IntervalIntersectionIterator(intersection, t, query)
        return new(intersection, t, query)
    end

    function IntervalIntersectionIterator()
        return new(Intersection{K, V, B}())
    end
end

# Intersect an interval tree t with a single interval, returning an iterator
# over the intersecting (key, value) pairs in t.
function Base.intersect{K, V, B}(t::IntervalBTree{K, V, B}, query0::Tuple{Any, Any})
    query = Interval{K}(query0[1], query0[2])
    return intersect(t, query)
end

function Base.intersect{K, V, B}(t::IntervalBTree{K, V, B}, first::K, last::K)
    query = Interval{K}(first, last)
    return intersect(t, query)
end

function Base.intersect{K, V, B}(t::IntervalBTree{K, V, B}, query::AbstractInterval{K})
    return IntervalIntersectionIterator{K, V, B}(Intersection{K, V, B}(), t, query)
end

function Base.start{K, V, B}(it::IntervalIntersectionIterator{K, V, B})
    return firstintersection!(it.t, it.query, it.intersection)
end

function Base.next{K, V, B}(it::IntervalIntersectionIterator{K, V, B}, ::Void)
    intersection = it.intersection
    entry = intersection.node.entries[intersection.index]
    nextintersection!(intersection.node, intersection.index,
                      it.query, intersection)
    return entry, nothing
end

function Base.done{K, V, B}(it::IntervalIntersectionIterator{K, V, B}, ::Void)
    return it.intersection.index == 0
end

type IntersectionIterator{K, V1, B1, V2, B2}
    t1::IntervalBTree{K, V1, B1}
    t2::IntervalBTree{K, V2, B2}

    # if true, use successive intersection, if false iterative.
    successive::Bool

    # true if done
    isdone::Bool

    # successive intersection state
    intersection::Intersection{K, V2, B2}
    t1_state::IntervalBTreeIteratorState{K, V1, B1}
    t1_value::V1

    # iterative intersection state
    u::Nullable{LeafNode{K, V1, B1}}
    v::Nullable{LeafNode{K, V2, B2}}
    w::Nullable{LeafNode{K, V2, B2}}
    i::Int
    j::Int
    k::Int

    function IntersectionIterator(t1, t2, successive::Bool)
        return new(t1, t2, successive)
    end
end

function Base.start{K, V1, B1, V2, B2}(it::IntersectionIterator{K, V1, B1, V2, B2})
    it.isdone = true
    if it.successive
        # Successive Query Intersection: intersect by repeatedly performing
        # queries against on tree ot the other.
        it.intersection = Intersection{K, V2, B2}()
        t1_state = start(it.t1)
        while !done(it.t1, t1_state)
            t1_value, t1_state = next(it.t1, t1_state)
            firstintersection!(it.t2, t1_value, it.intersection)
            if it.intersection.index != 0
                it.t1_state = t1_state
                it.t1_value = t1_value
                it.isdone = false
                return nothing
            end
        end
    else
        # Iterative Intersection: Intersect by iterating through the two
        # collections in unison.
        #
        # Notation notes. We proceed by finding all the the intervals
        # in t2 that intersect u.keys[i], before proceeding to the next
        # key in t1.
        #
        # We use v to keep track of the first leaf node in t2 that might contain an
        # intersecting interval, # and with (w, j) the current intersecting entry in
        # t2. We need to hang onto v since we will need to jump back when the
        # entry in t1 gets incremented.
        #
        # The thing to keep in mind is that this is just like the merge operation in
        # mergesort, except that some backtracking is needed since intersection
        # isn't as simple as ordering.

        it.u = isempty(it.t1) ? Nullable{LeafNode{K, V1, B1}}() : firstleaf(it.t1)
        it.v = isempty(it.t2) ? Nullable{LeafNode{K, V2, B2}}() : firstleaf(it.t2)
        it.w = it.v
        it.i, it.j, it.k = 1, 1, 1
        nextintersection!(it)
    end
    return nothing
end

function nextintersection!{K, V1, B1, V2, B2}(it::IntersectionIterator{K, V1, B1, V2, B2})
    u, v, w, i, j, k = it.u, it.v, it.w, it.i, it.j, it.k
    it.isdone = true
    while !isnull(u) && !isnull(v)
        unode = get(u)
        vnode = get(v)
        if isnull(w)
            @nextleafkey(unode, u, i)
            w, k = v, j
            continue
        end
        wnode = get(w)

        uentry = unsafe_getindex(unode.entries, i)
        ventry = unsafe_getindex(vnode.entries, j)
        wentry = unsafe_getindex(wnode.entries, k)

        if last(uentry) < first(ventry)
            @nextleafkey(unode, u, i)
        elseif last(ventry) < first(uentry)
            @nextleafkey(vnode, v, j)
        elseif first(wentry) <= last(uentry)
            if intersects(uentry, wentry)
                it.isdone = false
                break
            end
            @nextleafkey(wnode, w, k)
        else
            @nextleafkey(unode, u, i)
            w, k = v, j
        end
    end

    it.u = u
    it.v = v
    it.w = w
    it.i = i
    it.j = j
    it.k = k
    return
end

function Base.next{K, V1, B1, V2, B2}(it::IntersectionIterator{K, V1, B1, V2, B2}, state)
    if it.successive
        intersection = it.intersection
        entry = unsafe_getindex(intersection.node.entries, intersection.index)
        value = (it.t1_value, entry)
        nextintersection!(intersection.node, intersection.index, it.t1_value,
                          intersection)

        t1 = it.t1
        t2 = it.t2
        t1_state = it.t1_state
        t1_value = it.t1_value
        while intersection.index == 0 && !done(t1, t1_state)
            t1_value, t1_state = next(t1, t1_state)
            firstintersection!(t2, t1_value, intersection)
        end
        it.isdone = intersection.index == 0
        it.t1_state = t1_state
        it.t1_value = t1_value
    else
        u = get(it.u)
        w = get(it.w)
        value = (unsafe_getindex(u.entries, it.i),
                 unsafe_getindex(w.entries, it.k))
        @nextleafkey(w, it.w, it.k)
        nextintersection!(it)
    end
    return (value, nothing)
end

function Base.done{K, V1, B1, V2, B2}(it::IntersectionIterator{K, V1, B1, V2, B2}, state)
    return it.isdone
end

# Intersect two interval trees, returning an iterator yielding values of the
# form:
#   ((key1, value1), (key2, value2))
#
# Where key1 is from the first tree and key2 from the second, and they
# intersect.
function Base.intersect{K, V1, B1, V2, B2}(t1::IntervalBTree{K, V1, B1},
                                           t2::IntervalBTree{K, V2, B2};
                                           method=:auto)
    # We decide heuristically which intersection algorithm to use.
    n = length(t1)
    m = length(t2)

    iterative_cost  = n + m
    successive_cost = n * log(1 + m)

    if method == :auto
        iterative_cost  = n + m
        successive_cost = n * log(1 + m)
        if iterative_cost < successive_cost
            return IntersectionIterator{K, V1, B1, V2, B2}(t1, t2, false)
        else
            return IntersectionIterator{K, V1, B1, V2, B2}(t1, t2, true)
        end
    elseif method == :successive
        return IntersectionIterator{K, V1, B1, V2, B2}(t1, t2, true)
    elseif method == :iterative
        return IntersectionIterator{K, V1, B1, V2, B2}(t1, t2, false)
    else
        error("No such intersection method: $method")
    end
end

# Diagnostics
# -----------

# Dumb tree printing, useful only for debuging.
#function showtree(io::IO, t::IntervalBTree)
    #showtree(io, t.root, 0)
#end

#function showtree(t::IntervalBTree)
    #showtree(STDOUT, t)
#end

#function showtree(io::IO, t::InternalNode, indent::Int)
    #for (i, child) in enumerate(t.children)
        #showtree(io, child, indent+1)
        #if i <= length(t.keys)
            #print(io, repeat("  ", indent), t.keys[i], ":\n")
        #end
    #end
#end

#function showtree(io::IO, t::LeafNode, indent::Int)
    #for (k, v) in zip(t.keys, t.values)
        #print(io, repeat("  ", indent), k, ": ", v, "\n")
    #end
#end

include("map.jl")

if isdefined(Base, :iteratorsize)
    Base.iteratorsize(::IntervalFromIterator)         = Base.SizeUnknown()
    Base.iteratorsize(::IntervalIntersectionIterator) = Base.SizeUnknown()
    Base.iteratorsize(::IntersectionIterator)         = Base.SizeUnknown()
    Base.iteratorsize(::IntervalKeyIterator)          = Base.SizeUnknown()
    Base.iteratorsize(::IntervalValueIterator)        = Base.SizeUnknown()
end

end  # module IntervalTrees