text:
"""
    This provides a few plotting routines, for example, for CreepLaws
"""

using LaTeXStrings
using Unitful
using Parameters
using ..Units
using ..MaterialParameters
using ..MeltingParam
using .Plots

using GeoParams: AbstractMaterialParam, AbstractMaterialParamsStruct
using .MaterialParameters.CreepLaw: CreepLawVariables, computeCreepLaw_TauII, AbstractCreepLaw
using .MaterialParameters.HeatCapacity: AbstractHeatCapacity, compute_heatcapacity
using .MaterialParameters.Conductivity: AbstractConductivity, compute_conductivity
using .MeltingParam: AbstractMeltingParam, compute_meltfraction

export
    PlotStressStrainrate_CreepLaw,
    PlotHeatCapacity,
    PlotConductivity,
    PlotMeltFraction,
    PlotPhaseDiagram

"""
    PlotStressStrainrate_CreepLaw(x::AbstractCreepLaw; p::CreepLawParams=nothing, Strainrate=(1e-18,1e-12), CreatePlot::Bool=false)

Plots deviatoric stress versus deviatoric strain rate for a single creeplaw.
    Note: if you want to create plots or use the `CreatePlot=true` option you need to install the `Plots.jl` package in julia
    which is not added as a dependency here (as it is a rather large dependency).

# Example 1
```julia-repl
julia> x=LinearViscous()
Linear viscosity: η=1.0e20 Pa s
julia> Tau_II, Eps_II,  = PlotStressStrainrate_CreepLaw(x);
```
Next you can plot this with
```julia-repl
julia> using Plots;
julia> plot(ustrip(Eps_II),ustrip(Tau_II), xaxis=:log, yaxis=:log,xlabel="strain rate [1/s]",ylabel="Dev. Stress [MPa]")
```
Note that `ustrip` removes the units of the arrays, as many of the plotting packages don't know how to deal with that.

You could also have done:
```julia-repl
julia> using Plots;
julia> Tau_II, Eps_II, pl = PlotStressStrainrate_CreepLaw(x,CreatePlot=true);
```
which will generate the following plot
![subet1](./assets/img/Stress_Strainrate_LinearViscous.png)

The plot can be customized as
```julia-repl
julia> plot(pl, title="Linear viscosity", linecolor=:red)
```
See the [Plots.jl](https://github.com/JuliaPlots/Plots.jl) package for more options.

"""
function PlotStressStrainrate_CreepLaw(x::AbstractCreepLaw; p=nothing, Strainrate=(1e-18/s,1e-12/s), CreatePlot::Bool=false)

    if isnothing(p); p = CreepLawParams();  end

    if isDimensional(x)==false
       error("The struct with Creep Law parameters: $(typeof(x)) should be in dimensional units for plotting. You can use Dimensionalize! to do that.")
    end

    # Define strainrate
    Eps_II = range(ustrip(Strainrate[1])/s, stop=ustrip(Strainrate[2])/s, length=101)/s
    Tau_II = computeCreepLaw_TauII(Eps_II, x, p)                  # deviatoric stress

    # Transfer to GeoUnits
    Eps_II = GeoUnit(Eps_II);
    Tau_II = GeoUnit(Tau_II/1e6);

    if CreatePlot
        try
            pl = plot(ustrip(Eps_II), ustrip(Tau_II),
                xaxis=:log, xlabel=L"\textrm{deviatoric strain rate  } \dot{\varepsilon}_{II} \textrm{    [1/s]}",
                yaxis=:log, ylabel=L"\textrm{deviatoric stress  }\tau_{II} \textrm{    [MPa]}",
                legend=false,show = true)
        catch
            error("It seems that you did not install, or did not load Plots.jl. For plotting, please add that with `add Plots` in the package manager and type `using Plots` before running this.")
        end

        return Tau_II, Eps_II, pl
    else
        return Tau_II, Eps_II
    end

end

"""
    T,Cp,plt = PlotHeatCapacity(cp::AbstractHeatCapacity; T=nothing, plt=nothing, lbl=nothing)

Creates a plot of temperature `T` vs. heat capacity, as specified in cp (which can be temperature-dependent).

# Optional parameters
- T: temperature range
- plt: a previously generated plotting object
- lbl: label of the curve

# Example
```
julia> cp = T_HeatCapacity_Whittacker()
julia> T,Cp,plt = PlotHeatCapacity(cp)
```
you can now save the figure to disk with:
```
julia> using Plots
julia> savefig(plt,"Tdependent_heatcapacity.png")
```

"""
function PlotHeatCapacity(cp::AbstractHeatCapacity; T=nothing, plt=nothing, lbl=nothing)

    if isnothing(T)
        T = (273:10:1250)*K
    end

    Cp       =   compute_heatcapacity(T,cp)
    if length(Cp) == 1
        Cp = ones(size(T))*Cp
    end

    if isnothing(plt)
        plt = plot(ustrip(T), ustrip(Cp), label=lbl)
    else
        plt = plot!(ustrip(T), ustrip(Cp), label=lbl)
    end
    plot!(plt,   xlabel="Temperature [$(unit(T[1]))]",
                 ylabel="Cp [$(unit(Cp[1]))]")
    gui(plt)

    return T,Cp, plt
end

"""
    T,Kk,plt = PlotConductivity(cp::AbstractConductivity; T=nothing, plt=nothing, lbl=nothing)

Creates a plot of temperature `T` vs. thermal conductivity, as specified in `k` (which can be temperature-dependent).

# Optional parameters
- `T`: temperature range
- `plt`: a previously generated plotting object
- `lbl`: label of the curve

# Example
```
julia> k = T_Conductivity_Whittacker()
julia> T,KK,plt = PlotConductivity(k)
```
you can now save the figure to disk with:
```
julia> using Plots
julia> savefig(plt,"Tdependent_conductivity.png")
```

"""
function PlotConductivity(k::AbstractConductivity; T=nothing, P=nothing, plt=nothing, lbl=nothing)

    if isnothing(T)
        T = (273:10:1250)*K
    end
    if isnothing(P)
        P = 1e6Pa*ones(size(T))
    end

    Cond       =   compute_conductivity(P,T,k)
    if length(Cond) == 1
        Cond = ones(size(T))*Cond
    end

    if isnothing(plt)
        plt = plot(ustrip(T), ustrip(Cond), label=lbl)
    else
        plt = plot!(ustrip(T), ustrip(Cond), label=lbl)
    end
    plot!(plt,   xlabel="Temperature [$(unit(T[1]))]",
                 ylabel="Thermal conductivity [$(unit(Cond[1]))]")
    gui(plt)

    return T,Cond, plt
end

"""
    T,phi,plt = PlotMeltFraction(p::AbstractMeltingParam; T=nothing, plt=nothing, lbl=nothing)

Creates a plot of temperature `T` vs. melt fraction, as specified in `p`.
The 1D curve can be evaluated at a specific pressure `P` which can be given as a scalar or as an array of the same size as `T`

# Optional parameters
- `T`: temperature range
- `P`: pressure
- `plt`: a previously generated plotting object
- `lbl`: label of the curve

# Example
```
julia> p        =  MeltingParam_Caricchi()
julia> T,phi,plt = PlotMeltFraction(p)
```
you can now save the figure to disk with:
```
julia> using Plots
julia> savefig(plt,"MeltFraction.png")
```

"""
function PlotMeltFraction(p::AbstractMeltingParam; T=nothing, P=nothing, plt=nothing, lbl=nothing)

    if isnothing(T)
        T = (500:10:1500)*K
    end
    T_C = ustrip(T) .- 273.15

    if isnothing(P)
        P = 1e6Pa*ones(size(T))
    end
    if length(P) == 1
        P = P*ones(size(T))
    end

    phi       =   compute_meltfraction(P,T,p)
    if length(phi) == 1
        phi = ones(size(T))*phi
    end

    if isnothing(plt)
       plt = plot(T_C, ustrip(phi), label=lbl)
    else
       plt = plot!(T_C, ustrip(phi), label=lbl)
    end
    plot!(plt,   xlabel="Temperature [C]",
                 ylabel="Melt Fraction \\Phi")
    gui(plt)

    return T,phi, plt
end

"""
    plt, data, Tvec, Pvec = PlotPhaseDiagram(p::PhaseDiagram_LookupTable; fieldname::Symbol, Tvec=nothing, Pvec=nothing)

Plots a phase diagram as a function of `T` (x-axis) and `P` (y-axis).
We either use the default ranges of the diagram, or you can specify the temperature and pressure ranges (while specifying units).
The return arguments are the plotting object `plt` (so you can modify properties) as well as the data that is being plotted

Example
=======
```julia
julia> PD_data =  Read_LaMEM_Perple_X_Diagram("Peridotite.in")
Perple_X/LaMEM Phase Diagram Lookup Table:
                      File    :   Peridotite.in
                      T       :   293.0 - 1573.000039
                      P       :   1.0e7 - 2.9999999944e9
                      fields  :   :meltRho, :meltRho, :meltFrac, :rockRho, :Rho, :rockVp
                                  :rockVs, :rockVpVs, :meltVp, :meltVs, :meltVpVs
                                  :Vp, :Vs, :VpVs, :cpxFrac
julia> PlotPhaseDiagram(PD_data,:meltFrac, Tvec=(100:1:1400).*C, Pvec=(.1:.1:30).*kbar )
```
This will generate the following plot
![subet2](./assets/img/PhaseDiagram.png)

You can also use the default pressure/temperature ranges in the diagrams:
```julia
julia> PlotPhaseDiagram(PD_data,:Rho)
```

"""
function PlotPhaseDiagram(p::PhaseDiagram_LookupTable, fieldn::Symbol; Tvec=nothing, Pvec=nothing)

    data = getfield( p, fieldn)
    if isnothing(Tvec)
        Tvec_K  =   data.itp.knots[1]
        Tvec    =   Tvec_K;
    else
        Tvec_K  =   Float64.(uconvert.(K,Tvec))
    end
    if isnothing(Pvec)
        Pvec_Pa =   data.itp.knots[2]
        Pvec    =   Pvec_Pa
    else
        Pvec_Pa =   Float64.(uconvert.(Pa,Pvec))
    end

    data_scalar = data(ustrip.(Tvec_K), ustrip.(Pvec_Pa) )

    plt = heatmap(ustrip.(Tvec), ustrip.(Pvec), data_scalar',
                title  = string(fieldn),
                xlabel = "T [$(unit(Tvec[1]))]",
                ylabel = "P [$(unit(Pvec[1]))]", c=:batlow)

    display(plt)

    return  plt, data_scalar, Tvec, Pvec
end


text:
using zfp_jll

# zfp constants
const HEADER_MAGIC = 1
const HEADER_META = 2
const HEADER_MODE = 4
const HEADER_FULL = 7

# zfp types declaration and size
zfp_type(::Type{Int32}) = 1
zfp_type(::Type{Int64}) = 2
zfp_type(::Type{Float32}) = 3
zfp_type(::Type{Float64}) = 4
zfp_type(::Type) = 0

@enum ZfpType begin
    zfp_type_none = 0
    zfp_type_int32 = 1
    zfp_type_int64 = 2
    zfp_type_float = 3
    zfp_type_double = 4
end

#TODO mapping from C to Julia seems to be inconsistent here
@enum ZfpExecPolicy begin
    zfp_exec_serial = 0     # serial execution (default)
    zfp_exec_omp    = 1     # OpenMP multi-threaded execution
    zfp_exec_cuda   = 2     # CUDA parallel execution
end

function zfp_type(i::Int)
    i == 1 && return Int32
    i == 2 && return Int64
    i == 3 && return Float32
    i == 4 && return Float64
    throw(TypeError())
end

"""Size of zfp types (Int32,Int64,Float32,Float64) in bytes."""
zfp_type_size(i::Int64) = ccall((:zfp_type_size,libzfp),Int64,(Int64,),i)
zfp_type_size(::Type{T}) where T = zfp_type_size(zfp_type(T))

# READ IN ARRAYS
struct ZfpField
    type::ZfpType
    nx::Cuint
    ny::Cuint
    nz::Cuint
    nw::Cuint
    sx::Cint
    sy::Cint
    sz::Cint
    sw::Cint
    data::Ptr{Cvoid}
end

ZfpField(field::Ptr) = unsafe_load(Ptr{ZfpField}(field))

"""Pass a 1-D array into a zfp_field in C, in Julia only as Ptr{Cvoid}."""
function zfp_field(A::Array{T,1}) where T
    n = length(A)
    field = ccall((:zfp_field_1d,libzfp),Ptr{Cvoid},
        (Ptr{Cvoid},Cint,Cuint),A,zfp_type(T),n)
    sx = strides(A)[1]
    ccall((:zfp_field_set_stride_1d, libzfp), Cvoid, (Ptr{Cvoid}, Cint), field, sx)
    return field
end

"""Pass a 2-D array into a zfp_field in C, in Julia only as Ptr{Cvoid}."""
function zfp_field(A::Array{T,2}) where T
    nx,ny = size(A)
    field = ccall((:zfp_field_2d,libzfp),Ptr{Cvoid},
        (Ptr{Cvoid},Cint,Cuint,Cuint),A,zfp_type(T),nx,ny)
    sx,sy = strides(A)
    ccall((:zfp_field_set_stride_2d, libzfp), Cvoid, (Ptr{Cvoid}, Cint, Cint), field, sx, sy)
    return field
end

"""Pass a 3-D array into a zfp_field in C, in Julia only as Ptr{Cvoid}."""
function zfp_field(A::Array{T,3}) where T
    nx,ny,nz = size(A)
    field = ccall((:zfp_field_3d,libzfp),Ptr{Cvoid},
        (Ptr{Cvoid},Cint,Cuint,Cuint,Cuint),A,zfp_type(T),nx,ny,nz)
    sx,sy,sz = strides(A)
    ccall((:zfp_field_set_stride_3d, libzfp), Cvoid, (Ptr{Cvoid}, Cint, Cint, Cint),
        field, sx, sy, sz)
    return field
end

"""Pass a 4-D array into a zfp_field in C, in Julia only as Ptr{Cvoid}."""
function zfp_field(A::Array{T,4}) where T
    nx,ny,nz,nw = size(A)
    field = ccall((:zfp_field_4d,libzfp),Ptr{Cvoid},
        (Ptr{Cvoid},Cint,Cuint,Cuint,Cuint,Cuint),A,zfp_type(T),nx,ny,nz,nw)
    sx,sy,sz,sw = strides(A)
    ccall((:zfp_field_set_stride_4d, libzfp), Cvoid, (Ptr{Cvoid}, Cint, Cint, Cint, Cint),
        field, sx, sy, sz, sw)
    return field
end

"""Allocate an empty zfp field."""
zfp_field_alloc() = ccall((:zfp_field_alloc,libzfp),Ptr{Cvoid},(Ptr,),C_NULL)

"""Free the zfp field."""
function zfp_field_free(field::Ptr{Cvoid})
    ccall((:zfp_field_free,libzfp),Cvoid,(Ptr{Cvoid},),field)
end

"""Return type of a zfp field."""
zfp_field_type(field::Ptr{Cvoid}) = zfp_type(
    ccall((:zfp_field_type,libzfp),Int,(Ptr{Cvoid},),field))

"""Return the dimensionality (1,2,3 or 4) of the zfp field."""
zfp_field_dimensionality(field::Ptr{Cvoid}) = ccall((:zfp_field_dimensionality,libzfp),
    Int,(Ptr{Cvoid},),field)

"""Return the dimensionality (1,2,3 or 4) of the zfp field."""
zfp_field_pointer(field::Ptr{Cvoid}) = ccall((:zfp_field_pointer,libzfp),
    Ptr{Cvoid},(Ptr{Cvoid},),field)

"""Associate a zfp field with a data pointer."""
zfp_field_set_pointer(field::Ptr{Cvoid},ptr::Ptr) = ccall(
    (:zfp_field_set_pointer,libzfp),Cvoid,(Ptr{Cvoid},Ptr{Cvoid}),field,ptr)

# COMPRESSION OPTIONS
"""Open a stream (=object that holds the (de)compression settings) for zfp."""
zfp_stream_open() = ccall((:zfp_stream_open,libzfp),Ptr{Cvoid},(Ptr,),C_NULL)

"""Open a stream (=object that holds the (de)compression settings) for zfp
from an exisiting bitstream used for storing the compressed array."""
function zfp_stream_open(bitstream::Ptr{Cvoid})
    ccall((:zfp_stream_open,libzfp),Ptr{Cvoid},(Ptr{Cvoid},),bitstream)
end

"""Set bitrate (=bits per value) to set the compression rate directly.
Should not be larger than the bits per value of the uncompressed array."""
function zfp_stream_set_rate(stream::Ptr{Cvoid},rate::Real,type::Type,dims::Integer)
    ccall((:zfp_stream_set_rate,libzfp),Cdouble,
        (Ptr{Cvoid},Cdouble,Cuint,Cuint),stream,Float64(rate),zfp_type(type),dims)
end

"""Set the precision (≈ mantissa bits per value) for compression."""
function zfp_stream_set_precision(stream::Ptr{Cvoid},precision::Integer)
    ccall((:zfp_stream_set_precision,libzfp),Cuint,
        (Ptr{Cvoid},Cuint),stream,UInt(precision))
end

"""Set the accuracy (>=max abs error) for compression."""
function zfp_stream_set_accuracy(stream::Ptr{Cvoid},tol::AbstractFloat)
    ccall((:zfp_stream_set_accuracy,libzfp),Cdouble,
        (Ptr{Cvoid},Cdouble),stream,Float64(tol))
end

"""Set the zfp compression to lossless = reversible."""
function zfp_stream_set_reversible(stream::Ptr{Cvoid})
    ccall((:zfp_stream_set_reversible,libzfp),Cvoid,(Ptr{Cvoid},),stream)
end

"""Initialize a zfp stream C struct holding the compression settings.
    Only a Ptr{Cvoid} is returned to Julia."""
function zfp_stream(::Type{T},
                    ndims::Int;
                    tol::Real=0,
                    precision::Real=0,
                    rate::Int=0) where T

    stream = zfp_stream_open()

    # set the compression options
    if tol > 0
        zfp_stream_set_accuracy(stream,tol)
    elseif precision > 0
        zfp_stream_set_precision(stream,precision)
    elseif rate > 0
        maxrate = 8*zfp_type_size(T)
        rate <= maxrate || @warn "Rate was set to $rate-bit > $maxrate-bit for type $T"
        zfp_stream_set_rate(stream,rate,T,ndims)
    else  # lossless
        zfp_stream_set_reversible(stream)
    end

    return stream
end

# BUFFER
"""Retrieve the max size of the compressed stream in bytes pre-compression."""
function zfp_stream_maximum_size(stream::Ptr{Cvoid},field::Ptr{Cvoid})
    ccall((:zfp_stream_maximum_size,libzfp),Int,
        (Ptr{Cvoid},Ptr{Cvoid}),stream,field)
end

"""Retrieve the actual size of the compressed stream in bytes post-compression."""
function zfp_stream_compressed_size(stream::Ptr{Cvoid})
    ccall((:zfp_stream_compressed_size,libzfp),Int,(Ptr{Cvoid},),stream)
end

"""Open a buffer (= pointer for a preallocated data array) into a bitstream
    for zfp to flush the compressed array into."""
function stream_open(buffer::Ptr,bufsize::Int)
    ccall((:stream_open,libzfp),Ptr{Cvoid},(Ptr{Cvoid},Int),buffer,bufsize)
end

"""Close the compressed bitstream."""
function stream_close(bitstream::Ptr{Cvoid})
    ccall((:stream_close,libzfp),Cvoid,(Ptr{Cvoid},),bitstream)
end

"""Connect the zfp stream (=object that holds the compression settings) to the
bitstream that will contain the compressed array."""
function zfp_stream_set_bit_stream(stream::Ptr{Cvoid},bitstream::Ptr{Cvoid})
    ccall((:zfp_stream_set_bit_stream,libzfp),Cvoid,(Ptr{Cvoid},Ptr{Cvoid}),
        stream,bitstream)
end

"""Rewind the data stream."""
function zfp_stream_rewind(stream::Ptr{Cvoid})
    ccall((:zfp_stream_rewind,libzfp),Cvoid,(Ptr{Cvoid},),stream)
end

"""Rewind the data stream."""
function zfp_stream_flush(stream::Ptr{Cvoid})
    ccall((:zfp_stream_flush,libzfp),Cvoid,(Ptr{Cvoid},),stream)
end

"""Close the zfp stream."""
function zfp_stream_close(stream::Ptr{Cvoid})
    ccall((:zfp_stream_close,libzfp),Cvoid,(Ptr{Cvoid},),stream)
end

"""Write the header into the stream, which includes the compression parameters."""
function zfp_write_header(stream::Ptr{Cvoid},field::Ptr{Cvoid},HEADER::Int)
    ccall((:zfp_write_header,libzfp),Int,(Ptr{Cvoid},Ptr{Cvoid},Cuint),
        stream,field,HEADER)
end

"""Write the header into the stream, which includes the compression parameters."""
function zfp_read_header(stream::Ptr{Cvoid},field::Ptr{Cvoid},HEADER::Int)
    ccall((:zfp_read_header,libzfp),Int,(Ptr{Cvoid},Ptr{Cvoid},Cuint),
        stream,field,HEADER)
end

# SET OPENMP NUMBER OF THREADS
"""Set the number of OpenMP threads for compression, also switches ZfpExecPolicy to OpenMP."""
function zfp_stream_set_omp_threads(stream::Ptr{Cvoid},nthreads::Integer)
    success = ccall((:zfp_stream_set_omp_threads,libzfp),Cuint,
        (Ptr{Cvoid},Cuint),stream,UInt(nthreads))
    success == 0 && throw("Enabling OpenMP failed.")
end

"""Return the current execution policy (serial/OpenMP/CUDA)."""
function zfp_stream_execution(stream::Ptr{Cvoid})
    ccall((:zfp_stream_execution,libzfp),ZfpExecPolicy,(Ptr{Cvoid},),stream)
end

"""Set the current execution policy to serial, OpenMP or CUDA."""
function zfp_stream_set_execution(stream::Ptr{Cvoid},execution::Symbol)
    if execution == :serial
        exec_policy = ZfpExecPolicy(0)
    elseif execution == :openmp
        exec_policy = ZfpExecPolicy(1)
    elseif execution == :cuda
        # exec_policy = ZfpExecPolicy(2)
        throw("CUDA currently unsupported for ZfpCompression.jl.")
    else
        throw("Execution $execution unsupported.")
    end

    success = ccall((:zfp_stream_set_execution,libzfp),Int,
        (Ptr{Cvoid},ZfpExecPolicy),stream,exec_policy)

    success == 0 && throw("Enabling $execution failed.")
end

# COMPRESSION AND DECOMPRESSION
"""Low-level C call to run the compression."""
function zfp_compress(stream::Ptr{Cvoid},field::Ptr{Cvoid})
    ccall((:zfp_compress,libzfp),Int,(Ptr{Cvoid},Ptr{Cvoid}),stream,field)
end

"""Low-level C call to run the decompression."""
function zfp_decompress(stream::Ptr{Cvoid},field::Ptr{Cvoid})
    ccall((:zfp_decompress,libzfp),Int,(Ptr{Cvoid},Ptr{Cvoid}),stream,field)
end

function zfp_compress(  src::AbstractArray{T};
                        write_header::Bool=true,
                        nthreads::Int=1,
                        kws...) where {T<:Union{Int32,Int64,Float32,Float64}}

    ndims = length(size(src))
    ndims in [1,2,3,4] || throw(DimensionMismatch("Zfp compression only for 1-4D array."))

    zfpstream = zfp_stream(T,ndims;kws...)  # initialize the compression
    field = zfp_field(src)                  # turn src array into zfp field

    # preallocate the compressed array
    bufsize = zfp_stream_maximum_size(zfpstream,field)
    dest = Vector{UInt8}(undef,bufsize)             # allocate as UInt8
    bitstream = stream_open(pointer(dest),bufsize)  # turn array into zfp pointer
    zfp_stream_set_bit_stream(zfpstream,bitstream)  # connect bitstream pointer to zfp struct
    zfp_stream_rewind(zfpstream)

    # write header
    if write_header && zfp_write_header(zfpstream,field,HEADER_FULL) == 0
        throw(error("Writing header failed."))
    end

    # Enable OpenMP multi-threading
    if nthreads > 1
        zfp_stream_set_omp_threads(zfpstream,nthreads)
    end

    # perform compression
    success = zfp_compress(zfpstream,field)
    success == 0 && throw(error("Zfp compression failed."))
    zfp_stream_flush(zfpstream)
    compressed_size = zfp_stream_compressed_size(zfpstream)

    # free and close
    zfp_field_free(field)
    zfp_stream_close(zfpstream)
    stream_close(bitstream)

    return dest[1:compressed_size]
end

function zfp_decompress!(   dest::AbstractArray{T},
                            src::Vector{UInt8};
                            kws...) where {T<:Union{Int32,Int64,Float32,Float64}}

    ndims = length(size(dest))
    ndims in [1,2,3,4] || throw(DimensionMismatch("Zfp compression only for 1-4D array."))

    zfpstream = zfp_stream(T,ndims;kws...)    # initialize decompression
    field = zfp_field(dest)             # turn destination array into zfp pointer

    # declare src as the bitstream to decompress and connect to zfp struct
    bufsize = zfp_stream_maximum_size(zfpstream,field)
    bitstream = stream_open(pointer(src),bufsize)
    zfp_stream_set_bit_stream(zfpstream,bitstream)
    zfp_stream_rewind(zfpstream)

    # perform decompression
    compressed_size = zfp_decompress(zfpstream,field)

    # free and close
    zfp_field_free(field)
    zfp_stream_close(zfpstream)
    stream_close(bitstream)

    # check for failure
    compressed_size == 0 && throw(error("Zfp decompression failed."))
    return nothing
end

function zfp_decompress(src::Vector{UInt8})

    field = zfp_field_alloc()
    bitstream = stream_open(pointer(src),length(src))
    zfpstream = zfp_stream_open(bitstream)
    if zfp_read_header(zfpstream,field,HEADER_FULL) == 0
        throw(error("Reading header failed."))
    end

    # read header via ccalls
    T = zfp_field_type(field)
    ndims = zfp_field_dimensionality(field)

    # convert field pointer into Julia struct to access nx,ny,nz,nw
    ZF = ZfpField(field)
    n = filter(!=(0),(ZF.nx,ZF.ny,ZF.nz,ZF.nw))

    output = Array{T,ndims}(undef,n...)
    zfp_field_set_pointer(field,pointer(output))

    compressed_size = zfp_decompress(zfpstream,field)

    # free and close
    zfp_field_free(field)
    zfp_stream_close(zfpstream)
    stream_close(bitstream)

    # check for failure
    compressed_size == 0 && throw(error("Zfp decompression failed."))
    return output
end

ext: jl


text:
using CSV

struct RPLHeader
    width::Int
    height::Int
    depth::Int
    offset::Int
    dataType::Type{<:Real}
    byteOrder::Symbol # :bigendian, :littleendian
    recordedBy::Symbol # :vector, :image
    function RPLHeader(
        w::Int,
        h::Int,
        d::Int,
        o::Int,
        dt::Type{<:Real},
        bo::Symbol,
        rb::Symbol,
    )
        @assert w > 0
        @assert h > 0
        @assert d > 0
        @assert o >= 0
        @assert dt in [
            Int8,
            Int16,
            Int32,
            Int64,
            UInt8,
            UInt16,
            UInt32,
            UInt64,
            Float16,
            Float32,
            Float64,
        ]
        @assert bo in [:littleendian, :bigendian]
        @assert rb in [:vector, :image]
        new(w, h, d, o, dt, bo, rb)
    end
end

readrplraw(rplio::IO, rawio::IO) = readraw(rawio, readrpl(rplio))

"""
    readrplraw(rplfilename::AbstractString, rawfilename::AbstractString)
    readrplraw(rplio::IO, rawio::IO)

Read a RPL/RAW file pair from IO or filename into an Array obect.  The reader supports :bigendian, :littleendian
ordering and :vector or :image alignment.  Since the Array can be very large it is beneficial to release and collected
the associated memory when you are done with the data by assigning `nothing` to the variable (or allowing it to go
out of scope) and then calling `GC.gc()` to force a garbage collection.

    * Ordering: The individual data items may be in `:littleendian` or `:bigendian`.
    ** `:littleendian` => Intel/AMD and others
    ** `:bigendian` => ARM/PowerPC/Motorola
    * Alignment:  The data in the file can be organized as `:vector` or `:image`.  However, the data will be
reorganized into 'vector' format when returned as a Array.
    **`:vector` => Spectrum/data vector as contiguous blocks by pixel
    ** `:image` => Each channel of data organized in image planes
    * Data types: signed/unsigned 8/16/32-bit integers or 16-bit/32-bit/64-bit floats

    readrplraw(filenamebase::AbstractString)::Array{<:Real}

Read the files filenamebase*".rpl" and filenamebase*".raw" into an Array.  Maintains the data type
of the values in the RAW file.

#### Standard LISPIX Parameters in .rpl File

.rpl files consist of 9 lines.  Each line consists of a 'key'<tab>'value' where there is one and only one tab and
possibly other space between the parameter name and parameter value. Parameter names are case-insensitive.
The first line in the files is "key<tab>value".  Subsequent lines contain the keys and values described in this
table.

| **key**      |  **value**    | **description**                          |
| :------------|---------------|:-----------------------------------------|
| width        | 849           |  pixels per row       integer            |
| height       | 846           |  rows                 integer            |
| depth        | 4096          |  images or spec pts   integer            |
| offset       | 0             |  bytes to skip        integer            |
| data-length  | 1             |  bytes per pixel      1, 2, 4, or 8      |
| data-type    | unsigned      |  signed, unsigned, or float              |
| byte-order   | dont-care     |  big-endian, little-endian, or dont-care |
| record-by    | vector        |  image, vector, or dont-care             |

This .rpl file indicates the image is 849 pixels wide and 846 pixels high, with 4096 levels in the depth dimension.
"""
function readrplraw(rplfilename::AbstractString, rawfilename::AbstractString)::Array{<:Real}
    open(rplfilename, read = true) do rplio
        open(rawfilename, read = true) do rawio
            return readrplraw(rplio, rawio)
        end
    end
end

readrplraw(filenamebase::AbstractString)::Array{<:Real} =
    readrplraw(String(filenamebase) * ".rpl", String(filenamebase) * ".raw")

function readrpl(io::IO)::RPLHeader
    w, h, d, o = -1, -1, -1, -1
    dl, dtv, bo, rb = -1, "", :unknown, :unknown
    for row in CSV.Rows(read(io))
        key, val = uppercase(strip(row[1])), strip(row[2])
        if key == "WIDTH"
            w = parse(Int, val)
        elseif key == "HEIGHT"
            h = parse(Int, val)
        elseif key == "DEPTH"
            d = parse(Int, val)
        elseif key == "OFFSET"
            o = parse(Int, val)
        elseif key == "DATA-LENGTH"
            dl = parse(Int, val)
        elseif key == "DATA-TYPE"
            dtv = uppercase(val)
        elseif key == "BYTE-ORDER"
            bo = uppercase(val) == "BIG-ENDIAN" ? :bigendian : :littleendian
        elseif key == "RECORD-BY"
            rb = uppercase(val) == "IMAGE" ? :image : :vector # :dontcare => :vector
        end
    end
    dt = UInt8
    if dtv == "SIGNED"
        if dl == 1
            dt = Int8
        elseif dl == 2
            dt = Int16
        elseif dl == 4
            dt = Int32
        elseif dl == 8
            dt = Int64
        else
            error("Unexpected data-length $dl")
        end
    elseif dtv == "UNSIGNED"
        if dl == 1
            dt = UInt8
        elseif dl == 2
            dt = UInt16
        elseif dl == 4
            dt = UInt32
        elseif dl == 8
            dt = UInt64
        else
            error("Unexpected data-length $dl")
        end
    elseif dtv == "FLOAT"
        if dl == 2
            dt = Float16
        elseif dl == 4
            dt = Float32
        elseif dl == 8
            dt = Float64
        else
            error("Unexpected data-length $dl")
        end
    else
        error("Unexepected data-type $dt.")
    end
    return RPLHeader(w, h, d, o, dt, bo, rb)
end

"""
    readraw(ios::IOStream, rpl::RPLHeader)::Array{<:Real}

Construct an Array from the data in binary from the specified stream.  Always returns the data in [d,w,h] order
regardless of how the data is organized on disk.
"""
function readraw(ios::IOStream, rpl::RPLHeader)::Array{<:Real}
    if rpl.recordedBy == :vector
        size = ( rpl.depth, rpl.height, rpl.width )
    else
        size = ( rpl.height, rpl.width, rpl.depth )
    end
    seek(ios, rpl.offset)
    data = Array{rpl.dataType}(undef, size...)
    # Data in a RPL file is organized ch -> row -> col
    for r in Base.OneTo(rpl.height), c in Base.OneTo(rpl.width)
        read!(ios, @view data[:, r, c])
    end
    nativeendian = (ENDIAN_BOM == 0x04030201 ? :littleendian : :bigendian)
    if rpl.byteOrder != nativeendian
        if nativeendian == :littleendian
            foreach(i -> data[i] = ntoh(data[i]), eachindex(data))
        else
            foreach(i -> data[i] = hton(data[i]), eachindex(data))
        end
    end
    res = data
    if rpl.recordedBy == :image
        # Reorder as :vector
        res = Array{rpl.dataType}(undef, rpl.depth, rpl.width, rpl.height)
        for d in 1:rpl.depth, h in 1:rpl.height, w in 1:rpl.width
            res[d, h, w] = data[h, w, d]
        end
    end
    return res
end

function writerplraw(rplbasefile::String, arr::AbstractArray{<:Real})
    open(rplbasefile * ".rpl", write = true) do rplio
        println(rplio, "key\tvalue")
        println(rplio, "width\t$(size(arr,2))")
        println(rplio, "height\t$(size(arr,3))")
        println(rplio, "depth\t$(size(arr,1))")
        println(rplio, "offset\t0")
        println(rplio, "data-length\t$(sizeof(eltype(arr)))")
        if eltype(arr) in (Int8, Int16, Int32, Int64)
            println(rplio, "data-type\tsigned")
        elseif eltype(arr) in (UInt8, UInt16, UInt32, UInt64)
            println(rplio, "data-type\tunsigned")
        elseif eltype(arr) in (Float16, Float32, Float64)
            println(rplio, "data-type\tfloat")
        else
            error("Unsupported type $(eltype(arr)) in write RPL/RAW.")
        end
        println(rplio, "byte-order\tlittle-endian")
        println(rplio, "record-by\tvector")
    end
    open(rplbasefile * ".raw", write = true) do rawio
        write(rawio, arr)
    end
end

max_stars_repo_name: NicholasWMRitchie/NeXLSpectrum

max_issues_repo_name: NicholasWMRitchie/NeXLSpectrum

max_forks_repo_name: NicholasWMRitchie/NeXLSpectrum


text:
module Glucose

include("data.jl")
include("model_bayes.jl")
include("model.jl")
include("test.jl")

end

max_stars_repo_name: zenna/expect

max_issues_repo_name: zenna/expect

max_forks_repo_name: zenna/expect


text:
# Benchmark the counting of alanine residues in a PDB file

using Bio.Structure

pdb_filepath = "pdbs/1AKE.pdb"
struc = read(pdb_filepath, PDB)

function count()
    alanineselector(res::AbstractResidue) = resnameselector(res, ["ALA"])
    return countresidues(struc, alanineselector)
end

# Run to JIT compile
count()

elapsed = @elapsed count()

println(elapsed)

max_stars_repo_name: chemlove/pdb-benchmarks

max_issues_repo_name: chemlove/pdb-benchmarks

max_forks_repo_name: chemlove/pdb-benchmarks


text:
using Preprocessing, Gencode

function vis(visfields,app,mesh)

nt = size(visfields,4);

if app.viselem == []
    ne = size(mesh.t,2);
    app.viselem = 1:ne;
end

if app.porder>1
    visorder = min(app.porder,8);
else
    visorder = app.porder;
end
xpe,telem,~,~,~ = Preprocessing.masternodes(visorder,app.nd,app.elemtype);
shape = Preprocessing.mkshape(app.porder,mesh.xpe,xpe,app.elemtype);
shape = shape[:,:,1]';

dgnodes = Preprocessing.createdgnodes(mesh.p,mesh.t[:,app.viselem],mesh.f[:,app.viselem],mesh.curvedboundary,mesh.curvedboundaryexpr,visorder);
cgnodes, cgelcon, cgcells, celltype = createcggrid(dgnodes,telem);
dgnodes = Preprocessing.createdgnodes(mesh.p,mesh.t[:,app.viselem],mesh.f[:,app.viselem],mesh.curvedboundary,mesh.curvedboundaryexpr,app.porder);

# find paraview executable
app.paraview = Preprocessing.findexec(app.paraview, app.version);

if nt==1
    tm = shape*reshape(visfields[:,:,app.viselem],(size(shape,2), size(visfields,2)*length(app.viselem)));
    tm = reshape(tm,(size(shape,1), size(visfields,2), length(app.viselem)));
    vtuwrite(app.visfilename, cgnodes, cgelcon, cgcells, celltype, app.visscalars, app.visvectors, tm);
    if length(app.paraview)>0
        str = app.paraview * " --data="* app.visfilename * ".vtu";
        run(Gencode.string2cmd(str), wait=false);
    end
else
    if length(app.visdt)==0
        app.visdt = 1;
    end
    pvdwrite(app.visfilename, cgnodes, cgelcon, cgcells, celltype, app.visscalars, app.visvectors, visfields[:,:,app.viselem,:],app.visdt,shape);
    if length(app.paraview)>0
        str = app.paraview * " --data=" * app.visfilename * ".pvd";
        run(Gencode.string2cmd(str), wait=false);
    end
end

return dgnodes

end

max_stars_repo_name: glwagner/Exasim

max_issues_repo_name: glwagner/Exasim

max_forks_repo_name: glwagner/Exasim


text:
abstract type PopulationDynamicsModel end

"""
    PopState{T<:Real}
        P::Matrix{T}

Hold the population state at a given time.
"""
struct PopState{T<:Real} <: Any
    P::Matrix{T}
end

vecstate(P::PopState) = vec(P.P)
sum(P::PopState) = sum(P.P)

getindex(P::PopState, i) = P.P[i]
getindex(P::PopState, i, j) = P.P[i, j]
setindex!(P::PopState, x, i) = setindex!(P.P, x, i)
setindex!(P::PopState, x, i, j) = setindex!(P.P, x, i, j)
copy(P::PopState) = PopState(copy(P.P))

length(P::PopState) = length(P.P)
size(P::PopState) = size(P.P)
size(P::PopState, n::Integer) = size(P.P, n)

"""
    Schaefer{T<:Real}
        r::T
        K::T

Parameters of a Schaefer population dynamics model.
"""
struct Schaefer{T<:Real} <: PopulationDynamicsModel
    r::T
    K::T
end

"""
    (S::Schaefer)(p::PopState)

Step a population state forward by one time increment. Population change is
based on a region-wide carrying capacity, but each cell steps forward
individually.
 """
function (S::Schaefer)(P::PopState)
    Ptot = sum(P)
    Pnew = copy(P)
    real_r = S.r * (1 - Ptot / S.K)
    # Here and below, can't use @. macro because then Pnew.P is lowered to
    # getfield.(Pnew, :P), which doesn't work (or make sense).
    Pnew.P .= P.P .+ real_r .* P.P
    Pnew
end

"""
    PellaTomlinson{T<:Real}
        r::T
        K::T
        m::T

A Pella-Tomlinson population dynamics model with growth rate `r`, carrying
capacity `K`, and shape parameter `m`. Population will change over time
following

```math
P_{t+1} = P_t + \\frac{r}{m-1} P_t \\left(1 - \\left(\\frac{\\sum P_t}{K}\\right)^{m-1}\\right)
```

Note that the Schaefer model is a special case of this model where ``m = 2``,
and the Fox model is a limiting case where ``m \\to 1``.
"""
struct PellaTomlinson{T<:Real} <: PopulationDynamicsModel
    r::T
    K::T
    m::T
end

"""
    (PT::PellaTomlinson)(P::PopState)

Step a population state forward by one time increment. Population change is
based on a region-wide carrying capacity but each cell steps forward
independently.
"""
function (PT::PellaTomlinson)(P::PopState)
    Ptot = sum(P)
    Pnew = copy(P)
    real_r = PT.r / (PT.m - 1) * (1 - (Ptot / PT.K) ^ (PT.m - 1))
    @. Pnew.P = P.P + real_r * P.P
    Pnew
end

"""
    StochasticProduction{M<:PopulationDynamicsModel, D<:Distribution}
        dynmod::M
        dist::D

Use `dynmod` population dynamics, and then apply multiplicative noise in the
resulting population.
"""
struct StochasticProduction{M<:PopulationDynamicsModel,
                            D<:Union{Distribution,DomainDistribution}} <:PopulationDynamicsModel
    dynmod::M
    dist::D
end

function (sp::StochasticProduction{M,D})(P::PopState) where {M,D<:DomainDistribution}
    Pnew = sp.dynmod(P)
    Pnew.P .= Pnew.P .* rand(sp.dist)
    Pnew
end

function (sp::StochasticProduction{M,D})(P::PopState) where {M,D<:UnivariateDistribution}
    Pnew = sp.dynmod(P)
    Pnew.P .= Pnew.P .* rand(sp.dist, size(Pnew)...)
    Pnew
end

function (sp::StochasticProduction{M,D})(P::PopState) where {M,D<:MultivariateDistribution}
    Pnew = sp.dynmod(P)
    Pnew.P .= Pnew.P .* reshape(rand(sp.dist), size(Pnew)...)
    Pnew
end


text:
"""
    F = JopAtan(spc, c)

where `F` is the shifted arc tangent operator F = arctan( x(t)/c ) + π/2

x(t) = 0 corresponds to F = π/2
if |x(t)| ≤ c, then π/4 ≤ F ≤ 3π/4
"""
JopAtan(spc::JetSpace{T}, c = 1) where {T} = JopNl(dom = spc, rng = spc, f! = JopAtan_f!, df! = JopAtan_df!, df′! = JopAtan_df′!, s = (c=T(c),))
export JopAtan

function JopAtan_f!(d::AbstractArray, m::AbstractArray; c)
    d .= atan.(m ./ c) .+ (π/2)
end

function JopAtan_df!(δd::AbstractArray, δm::AbstractArray; mₒ, c)
    δd .= (c ./ (c^2 .+ mₒ.^2)) .* δm
end

function JopAtan_df′!(δm::AbstractArray, δd::AbstractArray; mₒ, c)
    δm .= conj(c ./ (c^2 .+ mₒ.^2)) .* δd
end


text:
function run(args::Dict)

  cur_task = args["<task>"]

  if cur_task == nothing
    cur_task = "main"
  end

  task_params = args["<field>"]

  run(cur_task, task_params)

end

function run(cur_task::String, task_params::AbstractArray=[])

  target_name = bump()

  if !isdefined(Symbol(target_name))
    eval(parse("using $target_name"))
  end

  eval(parse("$target_name.$cur_task($(join(task_params,",")))"))

end


text:
function simulation_init()
    A = Ne
    C = Ms+1
    D = Ms-Ne+1
    return Simulation{A, C, D}(
    zeros(Float64, N+2, 3), #x
    zeros(Float64, N+2, 3), #v
    MVector{3, Float64}(0.0, 0.0, 0.0),    #Δ_no
    zeros(Int64, N, 12), #nn_arr
    zeros(Float64, N+2, 3), #dhp_neutral
    zeros(Float64, N+2, 3), #dhp_ion
    zeros(Float64, N+2, 3), #dhp_coup
    MVector{3, Float64}(0.0, 0.0, 0.0), #hp
    zeros(Float64, N+2, 3), #F
    MVector{A, Int16}(zeros(Int16, Ne)), #surfp
    MVector{C, Int16}(zeros(Int16, Ms+1)), #occnum
    MVector{D, Int16}(zeros(Int16, Ms-Ne+1)), #surfh
    MVector{A, Int16}(zeros(Int16, Ne)), #surfpinit
    MVector{A, Int16}(zeros(Int16, Ne)), #surfpnew
    MVector{1, Float64}(10.0*Å), # trajzmin
    MVector{1, Float64}(0.0), #trajtheta
    zeros(Float64, Ms+1, Ms+1), #H
    zeros(ComplexF64, Ms+1, Ne), #psi
    zeros(ComplexF64, Ms+1, Ne), #phi
    zeros(Float64, Ms+1, Ms+1), #eigvec_H
    MVector{C, Float64}(zeros(Float64, Ms+1)), #eigval_H
    zeros(Float64, Ms+1, Ms+1), #dhdea
    zeros(Float64, Ms+1, Ms+1), #dhdv
    MVector{1, ComplexF64}(zero(ComplexF64)), #akl
    MVector{1, Float64}(0.0), #akk
    zeros(Float64, Ms+1,Ms+1), #DM
    zeros(Float64, Ne, Ms+1-Ne), #blk
    zeros(Float64, Ms+1, Ms+1), #blk2
    zeros(Float64, Ne*(Ms+1-Ne)), #Pb
    zeros(Float64, tsteps), #storage_aop
    zeros(Float64, Ms+1, tsteps), #storage_op,
    zeros(Float64, 5, tsteps), #storage_e
    zeros(Float64, 6, tsteps), #storage_xno
    zeros(Float64, 6, tsteps), #storage_vno
    zeros(Float64, tsteps), #storage_temp
    zeros(Float64, tsteps), #storage_phop
    MVector{1, ComplexF64}(zero(ComplexF64)), #"phipsi"
    MVector{1, Float64}(0.0), #"Pbmaxest"
    zeros(Float64, N+2, 3),#vdot
    zeros(Float64, N+2, 3),#vtemp
    MVector{1, Float64}(0.0), #vscale
    zeros(Float64, tsteps), #KEt
    zeros(Float64, tsteps), #PEt
    zeros(Float64, Int(round(tsteps/2))),
    zeros(Float64, Int(round(tsteps/2))),
    zeros(Float64, Int(round(tsteps/2))),
    zeros(Int64, Int(round(tsteps/2))),
    MVector{1, Int64}(0),
    MVector{1, Int64}(0),
    MVector{1, Int64}(0)
    )
end

function simulation_constructor_x_v!(s::Simulation)

    vvib_no = sqrt(e_vib_i*2.0/(mass_arr[1]*mass_arr[2]/(mass_arr[1] + mass_arr[2])))		# Initial vibrational velocity of NO
    vz_no = sqrt(e_trans_i*2.0/(mass_arr[1] + mass_arr[2]))				# Initial z velocity of NO

    #initial position of NO
    x0no = zeros(Float64, 2, 3)

    xno = rand() #temporary variable
    #xno = 0.1#rand() #temporary variable
    xno = @. xno * [cell[1], cell[2], 1.0*Å]
    #xno[3] = xno[3] + 4.5*Å # Initial z position of NO = 10
    xno[3] = xno[3] + 10.5*Å
    theta_no = rand() # Initial orientation of NO, 0 = O-down, PI = N-down
    #theta_no = 0.1#rand() # Initial orientation of NO, 0 = O-down, PI = N-down
    theta_no = asin(2.0 * theta_no - 1.0) + π/2.0
    phi_no = rand()
    #phi_no = 0.1#rand()
    phi_no = 2.0*π*phi_no

    x0no[1, :] = xno .+ r_no/2.0 * [sin(phi_no)*sin(theta_no), cos(phi_no)*sin(theta_no), cos(theta_no)]
    x0no[2, :] = xno .- r_no/2.0 * [sin(phi_no)*sin(theta_no), cos(phi_no)*sin(theta_no), cos(theta_no)]
    #x0no[1, 3] = x0no[1, 3] +1Å
    xno = x0no[1,:] - x0no[2, :]
    xno = xno / norm(xno)

    #initial velocity of no
    v0no = zeros(Float64, 2, 3)
    phi_no = 0.0
    v0no[1, :] = vz_no * [sin(phi_no)*sin(phi_inc), -cos(phi_no)*sin(phi_inc), -cos(phi_inc)]
    v0no[2, :] = v0no[1, :]
    v0no[1, :] = v0no[1, :] + vvib_no * xno * mass_arr[2] / (mass_arr[1] + mass_arr[2])
    v0no[2, :] = v0no[2, :] - vvib_no * xno * mass_arr[1] / (mass_arr[1] + mass_arr[2])

    #add rotational velocity
    theta_no = rand()
    #theta_no = 0.1#rand()
    theta_no = asin(2.0 * theta_no - 1.0) + π/2.0
    vrot_no = @. [-xno[3], 0.0, xno[1]] * cos(theta_no)/sqrt(xno[1]^2 + xno[3]^2)
    vrot_no = @. vrot_no + [xno[1], -(xno[1]^2 + xno[3]^2)/xno[2], xno[3]] *
                sin(theta_no)/sqrt((xno[1]^2 + xno[3]^2) * (1 + (xno[1]^2 + xno[3]^2)/xno[2]^2))

    vvib_no = sqrt(2.0 * e_rot_i/(mass_arr[1] + mass_arr[2]))
    v0no[1, :] = @. v0no[1, :] + vvib_no * vrot_no * sqrt(mass_arr[2]/mass_arr[1])
    v0no[2, :] = @. v0no[2, :] - vvib_no * vrot_no * sqrt(mass_arr[1]/mass_arr[2])
    #stack information to x_au0 and v_au0
    v0 = [v0no; v_au0]
    x0 = [x0no; x_au0]
    s.x .= x0
    s.v .= v0
    s.Δ_no .= x0[1, :] - x0[2, :]
end

function simulation_constructor_nn!(s::Simulation)
    nn_lattice!(s)
end

function simulation_constructor_energy(s::Simulation)
    s.hp .= get_E_all(s)
end

function simulation_constructor_force(s::Simulation)
    get_F_all(s)
end


text:
using LightGraphs
using LinearAlgebra
using SparseArrays
using IterativeSolvers
using Plots
using GraphPlot
using ColorSchemes

function network_automorphism_classes(g::SimpleGraph{<:Integer}; tol=1e-4)
    r = lobpcg(float(adjacency_matrix((g))), true, 1, tol = 1e-10)
    x = copy(r.X[:,1]) # get the dominant eigenvector
    p = sortperm(x)
    x .= x[p]
    k = 1
    classes = [[p[1]]]
    for i in 2:length(x)
        if abs(x[i]-x[i-1]) < tol
            push!(classes[k], p[i])
        else
            k += 1
            push!(classes, [p[i]])
        end
    end
    return classes, r.X[:,1]
end

A = [0 1 1 0 0 0 0 0;
     1 0 1 0 0 0 0 0;
     1 1 0 1 1 1 0 0;
     0 0 1 0 0 0 1 1;
     0 0 1 0 0 0 1 1;
     0 0 1 0 0 0 1 1;
     0 0 0 1 1 1 0 0;
     0 0 0 1 1 1 0 0
    ]

g = SimpleGraph(A)
gplot(g)

cs = network_automorphism_classes(g)

function pairpermutation_matrix(n,i,j)
    P = spzeros(Int64, n,n)
    for k in 1:n
        if k != i || k != j
            P[k,k] = 1
        end
    end
    P[i,j] = 1
    P[j,i] = 1
    return P
end

function checkclasses(g::SimpleGraph{<:Integer}, cs)
    A = adjacency_matrix(g)
    N = nv(g)
    for c in cs
        n = length(c)
        if n > 1
            for i in 2:n # NB incorrect test
                P = pairpermutation_matrix(N,c[1],c[i])
                if A != P*A*P
                    return false
                end
            end
        end
    end
    return true
end

checkclasses(g, cs)

g = graphfamous("karate")
cs, v = network_automorphism_classes(g)
classes = filter(c->length(c)>1, cs)
allclasses = vcat(classes...)
ps = [(k,v[k]) for k in allclasses]
colors = ColorSchemes.Set1_4.colors;
nodecols = fill(colorant"lightgray", nv(g));
for i in 1:length(classes), v in classes[i]
    nodecols[v] = colors[i];
end
gplot(g, nodelabel=1:34, nodefillc = nodecols)

A = float(adjacency_matrix(g))
λ, V = eigen(Array(A))
[V[allclasses,i] for i in 1:34 if any(abs.(V[:,i]) .< 1e-6) ]


text:
using FileIO
using FCSFiles
using GeometricalPredicates

#= Gets flow data for a dictionary of tubes.

Given a dictionary of index => fcs file name, returns a dictionary of
index => filtered flow values.
=#
function get_facs_data(tubedict)
    # want to use B1-A for information -- FITC
    facs_data = Dict(i => get_filtered_flow_data(tubedict[i]) for i in collect(keys(tubedict)))

    if debug
        println("got facs data")
    end

    return facs_data
end

#= Filters flow data for live cells

Takes arguments for FSC-A, SSC-A, and B1A -- uses a polygon gate made in
Cytoflow to filter B1A values if the corresponding FSC-A and SSC-A values
indicate that the cell is a live singlet.
=#
function filter_flow_data(fscs, sscs, mchs, b1a)
    min_coords = [minimum(fscs), minimum(sscs)]
    max_coords = [maximum(fscs), minimum(sscs)]

    # GeometricalPredicates requires all coordinates to be in [1, 2)
    fscs_scaled = [(f-min_coords[1])/(max_coords[1]-min_coords[1]) + 1.0 for f in fscs]
    sscs_scaled = [(s-min_coords[2])/(max_coords[2]-min_coords[2]) + 1.0 for s in sscs]

    point1 = (50000,  10000)
    point2 = (200000, 10000)
    point3 = (200000, 100000)
    point4 = (50000,  100000)
    point5 = (50000,  10000)

    #  point1 = (113097.94767349683, 188528.74478919798)
    #  point2 = (158699.69036009547, 77818.35897898073)
    #  point3 = (86319.3037719132,   32617.497996363905)
    #  point4 = (63789.55549038625,  44560.87239599226)
    #  point5 = (63789.55549038625,  104693.81475858176)
    #  point6 = (90599.38825827425,  179128.20042318373)
    #  point7 = (113097.94767349683, 188528.74478919798) # close the polygon
    #  points = [point1, point2, point3, point4, point5, point6]

    # Scales a point to be in the domain [1, 2) for use with GeometricalPredicates
    function scale_point(point, min_coords, max_coord)
        return GeometricalPredicates.Point(
                    (point[1]-min_coords[1])/(max_coords[1]-min_coords[1]) + 1.0,
                    (point[2]-min_coords[2])/(max_coords[2]-min_coords[2]) + 1.0)
    end

    point1 = scale_point(point1, min_coords, max_coords)
    point2 = scale_point(point2, min_coords, max_coords)
    point3 = scale_point(point3, min_coords, max_coords)
    point4 = scale_point(point4, min_coords, max_coords)
    point5 = scale_point(point5, min_coords, max_coords)
    #  point6 = scale_point(point6, min_coords, max_coords)
    #  point7 = scale_point(point7, min_coords, max_coords)

    livegate = Polygon(point1, point2, point3, point4, point5)
    #  livegate = Polygon(point1, point2, point3, point4, point5, point6, point7)

    b1a_filtered = Float64[]

    for i in 1:length(b1a)
        point = GeometricalPredicates.Point(fscs_scaled[i], sscs_scaled[i])
        if inpolygon(livegate, point)>0
            if (mchs[i] > 1000) # gate on mCherry-positive cells
                append!(b1a_filtered, b1a[i])
            end
        end
    end

    pct = (1.0 * length(b1a_filtered))/(1.0 * length(fscs_scaled))

    return b1a_filtered
end

#= Returns filtered flow data given path to FCS file

Accepts string path to FCS file, and returns the filtered flow data -- in this
case, we're only really concerned with B1-A values gated on FSC-A/SSC-A
=#
function get_filtered_flow_data(fcsfile)
    filedata = load(fcsfile)
    b1a = filedata["B1-A"]
    fsc = filedata["FSC-A"]
    ssc = filedata["SSC-A"]
    mchs = filedata["Y2-A"]

    return filter_flow_data(fsc, ssc, mchs, b1a)
end

#= Given a dictionary of input => FCS filepath, returns an array of [inputs. outputs]

Note here that each output is a specific value/read from the FCS file
=#
function get_data_array(facsdict)
    data_arr = []
    for input in collect(keys(facsdict))
        yvals = facsdict[input]
        append!(data_arr, [[input y] for y in yvals])
    end

    return [[[Float32(x[1]) for x in data_arr]]
            [[Float32(x[2]) for x in data_arr]]]
end

#= Given a dictionary of input => FCS files, returns a dict of input => samples =#
function get_samples(tubedict)
    samples_dict = Dict(
        i => get_filtered_flow_data(tubedict[i])
                    for i in collect(keys(tubedict)))
end


text:
using CUDArt
import CUDAdrv
using Base.Test

dev = CUDAdrv.CuDevice(0)
ctx = CUDAdrv.CuContext(dev)

CUDArt.clean_cache()   # for deterministic testing purposes

## basic compilation & execution

let
    CUDArt.@compile dev kernel """
    __global__ void kernel()
    {
    }
    """

    CUDAdrv.cudacall(kernel, 1, 1, ())
end

@test_throws CompileError let
    CUDArt.@compile dev kernel """
    __global__ void kernel()
    {
        invalid code
    }
    """
end

@test_throws CUDAdrv.CuError let
    CUDArt.@compile dev wrongname """
    __global__ void kernel()
    {
    }
    """
end

## argument passing

dims = (16, 16)
len = prod(dims)

CUDArt.@compile dev kernel_copy """
__global__ void kernel_copy(const float *input, float *output)
{
    int i = blockIdx.x * blockDim.x + threadIdx.x;

    output[i] = input[i];
}
"""

let
    input = round.(rand(Cfloat, dims) * 100)

    input_dev = CUDAdrv.CuArray(input)
    output_dev = CUDAdrv.CuArray{Cfloat}(dims)

    CUDAdrv.cudacall(kernel_copy, 1, len,
                     Tuple{Ptr{Cfloat}, Ptr{Cfloat}},
                     input_dev, output_dev)
    output = Array(output_dev)
    @test input ≈ output
end

CUDAdrv.destroy!(ctx)


text:
using Test, FinancialMonteCarlo;
@show "Black Scholes Model"
S0 = 100.0;
K = 100.0;
r = 0.02;
T = 1.0;
d = 0.01;
D = 90.0;

Nsim = 10000;
Nstep = 30;
sigma = 0.2;
mc = MonteCarloConfiguration(Nsim, Nstep);
mc1 = MonteCarloConfiguration(Nsim, Nstep, FinancialMonteCarlo.AntitheticMC());
toll = 0.8

rfCurve = ZeroRate(r);

FwdData = Forward(T)
EUData = EuropeanOption(T, K)
EUBin = BinaryEuropeanOption(T, K)
AMData = AmericanOption(T, K)
AmBin = BinaryEuropeanOption(T, K)
BarrierData = BarrierOptionDownOut(T, K, D)
AsianFloatingStrikeData = AsianFloatingStrikeOption(T)
AsianFixedStrikeData = AsianFixedStrikeOption(T, K)
Model = BlackScholesProcess(sigma, Underlying(S0, d));

display(Model)
ST = FinancialMonteCarlo.distribution(Model, rfCurve, mc, T);
@show FwdPrice = exp(-r * T) * sum(ST) / length(ST);

@test abs(FwdPrice - 99.1078451563562) < toll


text:
using Amb
using Test

@testset "Pythagorean triples" begin
    intbetween(lo, hi) = (require(lo <= hi); @amb lo intbetween(lo+1, hi))

    function triple(lo, hi)
        i = intbetween(lo, hi)
        j = intbetween(i, hi)
        k = intbetween(j, hi)
        require(i*i + j*j == k*k)
        (i, j, k)
    end
    @test collect(ambiter(()->triple(1,20))) == [(3, 4, 5), (5, 12, 13),
                                                 (6, 8, 10), (8, 15, 17),
                                                 (9, 12, 15), (12, 16, 20)]
end


text:
module IntervalTrees

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

using Base: notnothing

include("slice.jl")

"""
Types deriving from `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 type AbstractInterval{T} end

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

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

function Base.isless(u::AbstractInterval, v::AbstractInterval)
    return basic_isless(u, v)
end

"""
A basic interval.
"""
struct Interval{T} <: AbstractInterval{T}
    first::T
    last::T
end
Base.convert(::Type{Interval{T}}, range::AbstractRange{T}) where T =
    Interval(first(range), last(range))
Interval(range::AbstractRange{T}) where T = convert(Interval{T}, range)

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

IntervalValue(range::AbstractRange{K}, value::V) where {K, V} =
    IntervalValue(first(range), last(range), value)

valtype(::Type{IntervalValue{K,V}}) where {K, V} = V

value(i::IntervalValue{K, V}) where {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 type Node{K, V, B} end

mutable struct 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 end ofd this node
    maxend::K

    # Maximum child subtree end-points
    maxends::Slice{K, B}

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

    parent::Union{Nothing, InternalNode{K, V, B}}

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

    function InternalNode{K,V,B}() where {K,V,B}
        return new{K,V,B}(
            Slice{Interval{K}, B}(), zero(K), Slice{K, B}(), Slice{Node{K, V, B}, B}(),
            nothing,
            nothing,
            nothing)
    end
end

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

mutable struct LeafNode{K, V, B} <: Node{K, V, B}
    entries::Vector{V}
    keys::Vector{Interval{K}}
    count::UInt32 # number of stored entries

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

    parent::Union{Nothing, InternalNode{K, V, B}}

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

    function LeafNode{K,V,B}() where {K,V,B}
        return new{K,V,B}(
            Array{V}(undef, B), Array{Interval{K}}(undef, B), 0,
            zero(K),
            nothing,
            nothing,
            nothing)
    end
end

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

function Base.resize!(leaf::LeafNode{K, V, B}, n) where {K, V, B}
    @assert 0 <= n <= B
    leaf.count = n
end

function Base.splice!(leaf::LeafNode{K, V, B}, i::Integer) where {K, V, B}
    slice_splice!(leaf.keys, leaf.count, i)
    x, leaf.count = slice_splice!(leaf.entries, leaf.count, i)
    return x
end

function Base.insert!(leaf::LeafNode{K, V, B}, i::Integer, value::V) where {K, V, B}
    slice_insert!(leaf.keys, leaf.count, i, Interval{K}(first(value), last(value)))
    slice_insert!(leaf.entries, leaf.count, i, value)
    leaf.count += 1
end

function Base.push!(leaf::LeafNode{K, V, B}, value::V) where {K, V, B}
    leaf.count += 1
    leaf.entries[leaf.count] = value
    leaf.keys[leaf.count] = Interval{K}(first(value), last(value))
end

function Base.pop!(leaf::LeafNode{K, V, B}) where {K, V, B}
    x = leaf.entries[leaf.count]
    leaf.count -= 1
    return x
end

function Base.iterate(leaf::LeafNode, i::Int=1)
    return i ≤ length(leaf) ? (leaf.entries[i], i + 1) : nothing
end

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

    function IntervalBTree{K,V,B}() where {K,V,B}
        return new{K,V,B}(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{K,V,B}(entries::AbstractVector{V}) where {K,V,B}
        if !issorted(entries, lt=basic_isless)
            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{K,V,B}(LeafNode{K, V, B}(), 0)
        end

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

        maxends = Vector{K}(undef, numleaves)
        minkeys = Vector{Interval{K}}(undef, 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[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])
                    push!(parents[i].maxends, maxends[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{K,V,B}(children[1], n)
    end
end

# Default B-tree order
const 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 t.count
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 t.count == 0
end

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

mutable struct IntervalBTreeIteratorState{K, V, B}
    leaf::Union{Nothing, LeafNode{K, V, B}}
    i::Int
end

function Base.iterate(
        t::IntervalBTree{K, V, B},
        state::IntervalBTreeIteratorState{K, V, B}=iterinitstate(t)) where {K, V, B}
    if state.leaf === nothing || isempty(state.leaf)
        return nothing
    end
    leaf = notnothing(state.leaf)
    entry = leaf.entries[state.i]
    if state.i < length(leaf)
        state.i += 1
    else
        state.leaf = leaf.right
        state.i = 1
    end
    return entry, state
end

function iterinitstate(t::IntervalBTree{K,V,B}) where {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(node, 1)
end

# Iterate from a given starting from the first interval that intersects a given point.
struct IntervalFromIterator{K, V, B}
    t::IntervalBTree{K, V, B}
    p::K
end

Base.eltype(::Type{IntervalFromIterator{K, V, B}}) where {K, V, B} = V
Base.IteratorSize(::Type{IntervalFromIterator{K, V, B}}) where {K, V, B} = Base.SizeUnknown()

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

function Base.iterate(
        it::IntervalFromIterator{K, V, B},
        state::IntervalBTreeIteratorState{K, V, B}=iterinitstate(it)) where {K, V, B}
    return iterate(it.t, state)
end

function iterinitstate(it::IntervalFromIterator{K, V, B}) where {K, V, B}
    node, i = firstfrom(it.t, it.p)
    if i == 0
        return IntervalBTreeIteratorState{K, V, B}(nothing, i)
    else
        return IntervalBTreeIteratorState{K, V, B}(node, i)
    end
end

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

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

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

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

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

    return left, right
end

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

    m = length(left)

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

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

    resize!(left.children, div(m, 2))
    resize!(left.maxends, 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(t::InternalNode{K, V, B}) where {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(t::LeafNode{K, V, B}) where {K, V, B}
    maxend = zero(K)
    for i in 1:length(t)
        maxend = max(maxend, last(t.keys[i]))
    end
    return maxend
end

sameparent(u::Node{K, V, B}, v::Node{K, V, B}) where {K, V, B} = (u.parent == v.parent)

# 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

# Find the root node
function root(t::Node)
    while t.parent !== nothing
        t = notnothing(t.parent)
    end
    return t
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(t::LeafNode{K, V, B}) where {K, V, B}
    return t.keys[1]
end

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

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

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

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

    insert!(node, 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 = rightleaf.keys[1]

        # travel back up the tree setting maxend values and splitting as needed
        parent = node.parent
        leftnode = leftleaf
        child = leftleaf
        rightnode = rightleaf
        hassplit = true
        while parent !== nothing
            p = notnothing(parent)
            if hassplit
                i = findidx(p, entry) # key index
                j = i <= length(p) - 1 && entry < p.keys[i] ? i : i + 1 # child index
                p.maxends[j] = p.children[j].maxend
                insert!(p.children, j + 1, rightnode)
                insert!(p.maxends, j + 1, p.children[j+1].maxend)
                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))
                ifind = findfirst(isequal(child), p.children)
                if ifind === nothing
                    ifind = 0
                end
                p.maxends[ifind] = child.maxend
            end
            child = p
            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
            push!(t.root.maxends, leftnode.maxend)
            push!(t.root.maxends, rightnode.maxend)
            leftnode.parent = t.root
            rightnode.parent = t.root
        end
    else
        # travel back up the tree setting maxend values
        parent = node.parent
        child = node
        while parent !== nothing
            p = notnothing(parent)
            p.maxend = max(p.maxend, last(entry))
            ifind = findfirst(isequal(child), p.children)
            if ifind === nothing
                ifind = 0
            end
            p.maxends[ifind] = child.maxend
            parent = p.parent
            child = p
        end
    end

    return entry
end

# Deleting
# --------

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

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

function deletefirst!(t::IntervalBTree{K, V, B}, key::Interval{K}) where {K, V, B}
    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 = nothing
    end

    return t
end

# Indicate what steps are need to account for an updated child.
struct 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
struct 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 if present. Return a DeleteResult object.
function _deletefirst!(t::InternalNode{K, V, B}, key::Interval{K}) where {K, V, B}
    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)
        t.maxends[j] = t.children[j].maxend
        if j > 1
            t.maxends[j-1] = t.children[j-1].maxend
            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)
        t.maxends[j] = t.children[j].maxend
        if j < length(t)
            t.maxends[j+1] = t.children[j+1].maxend
            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)
        splice!(t.maxends, deleteidx)
        t.maxend = nodemaxend(t)
        if ans.fate == KEYFATE_DELETE
            t.maxends[j-1] = t.children[j-1].maxend
        elseif ans.fate == KEYFATE_DELETE_RIGHT
            t.maxends[j] = t.children[j].maxend
        end

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

        # borrow right
        if t.right !== nothing
            right = notnothing(t.right)
            if sameparent(right, t) && length(right) > minsize
                splice!(right.keys, 1)
                push!(t.children, splice!(right.children, 1))
                push!(t.maxends, splice!(right.maxends, 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 t.left !== nothing
            left = notnothing(t.left)
            if sameparent(left, t) && length(left) > minsize
                insert!(t.children, 1, pop!(left.children))
                insert!(t.maxends, 1, pop!(left.maxends))
                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 t.left !== nothing
            left = notnothing(t.left)
            if sameparent(left, t)
                merge!(left, t)
                left.right = t.right
                if t.right !== nothing
                    t.right = notnothing(t.right)
                    t.right.left = left
                end

                return DeleteResult(true, KEYFATE_DELETE)
            end
        end

        # merge with right
        if t.right !== nothing
            right = notnothing(t.right)
            if sameparent(right, t)
                merge!(t, right)
                if right.right !== nothing
                    right.right = notnothing(right.right)
                    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 # KEYFATE_NONE
        t.maxend = nodemaxend(t)
        t.maxends[j] = t.children[j].maxend
        return DeleteResult(true, KEYFATE_NONE)
    end
end

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

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

    splice!(t, 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 t.right !== nothing
        right = notnothing(t.right)
        if sameparent(right, t) && length(right) > minsize
            push!(t, splice!(right, 1))
            t.maxend = max(t.maxend, last(t.keys[t.count]))
            right.maxend = nodemaxend(right)

            return DeleteResult(true, KEYFATE_UPDATE_RIGHT)
        end
    end

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

            return DeleteResult(true, KEYFATE_UPDATE_LEFT)
        end
    end

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

    # merge with right
    if t.right !== nothing
        right = notnothing(t.right)
        if sameparent(right, t)
            merge!(t, right)
            if right.right !== nothing
                right.right = notnothing(right.right)
                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!(left::LeafNode{K, V, B}, right::LeafNode{K, V, B}) where {K, V, B}
    leftlen, rightlen = length(left), length(right)
    @assert leftlen + rightlen <= B
    resize!(left, leftlen + rightlen)
    copyto!(left.entries, leftlen+1, right.entries, 1, length(right))
    copyto!(left.keys, leftlen+1, right.keys, 1, length(right))
    left.maxend = nodemaxend(left)
    return
end

# Join two internal nodes into one
function merge!(left::InternalNode{K, V, B}, right::InternalNode{K, V, B}) where {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)
    resize!(left.maxends, leftlen + rightlen)
    copyto!(left.children, leftlen+1, right.children, 1, length(right.children))
    copyto!(left.maxends, leftlen+1, right.maxends, 1, length(right.maxends))
    for i in leftlen+1:length(left.children)
        left.children[i].parent = left
    end
    left.keys[leftlen] = minkey(left.children[leftlen+1])
    copyto!(left.keys, leftlen+1, right.keys, 1, length(right.keys))
    left.maxend = nodemaxend(left)
    return
end

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

function findidx(t::LeafNode{K, V, B}, key::AbstractInterval{K}) where {K, V, B}
    if isempty(t)
        return 0
    end
    # TODO: search keys
    return searchsortedfirst(t.entries, key, 1, Int(t.count),
                             Base.ord(basic_isless, identity, false))
end

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

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

function Base.haskey(t::InternalNode{K, V, B}, key::AbstractInterval{K}) where {K, V, B}
    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(t::IntervalBTree{K, V, B}, key0::Tuple{Any, Any}) where {K, V, B}
    key = Interval{K}(key0[1], key0[2])
    return haskey(t.root, key)
end

function Base.findfirst(t::LeafNode{K, V, B}, key::AbstractInterval{K}, f) where {K, V, B}
    i = findidx(t, key)
    while 1 <= i <= length(t) &&
          first(t.keys[i]) == first(key) &&
          last(t.keys[i]) == last(key)
        if f(t.entries[i], key)
            return t.entries[i]
        end
        i += 1
        if i > length(t) && t.right !== nothing
            t = notnothing(t.right)
            i = 1
        end
    end

    return nothing
end

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

true_cmp(a, b) = true

function Base.findfirst(t::IntervalBTree{K, V, B}, key::AbstractInterval{K},
                        f=true_cmp) where {K, V, B}
    return findfirst(t.root, key, f)
end

function Base.findfirst(t::IntervalBTree{K, V, B}, key0::Tuple{Any, Any},
                        f=true_cmp) where {K, V, B}
    return findfirst(t, Interval{K}(key0[1], key0[2]), f)
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(key1::AbstractInterval{K}, key2::AbstractInterval{K}) where 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(t::IntervalBTree{K, V, B}, query) where {K, V, B}
    return hasintersection(t.root, convert(K, query))
end

function hasintersection(t::InternalNode{K, V, B}, query::K) where {K, V, B}
    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(t::LeafNode{K, V, B}, query::K) where {K, V, B}
    if isempty(t) || t.maxend < query
        return false
    end

    for i in 1:length(t)
        if first(t.keys[i]) <= query <= last(t.keys[i])
            return true
        elseif query < first(t.keys[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.
mutable struct Intersection{K, V, B}
    index::Int
    node::LeafNode{K, V, B}

    Intersection{K,V,B}(index, node) where {K,V,B} = new{K,V,B}(index, node)
    Intersection{K,V,B}() where {K,V,B} = new{K,V,B}(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!(t::IntervalBTree{K, V, B},
                            query::AbstractInterval{K},
                            lower::Union{Nothing, V},
                            out::Intersection{K, V, B},
                            filter::F) where {F, K, V, B}
    return firstintersection!(t.root, query, lower, out, filter)
end

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

    (query_first, query_last) = (first(query), last(query))

    i = lower === nothing ? 1 : searchsortedfirst(t.keys, notnothing(lower))

    @inbounds while i <= length(t.children)
        if i > 1 && t.keys[i-1].first > query_last
            break
        end
        if t.maxends[i] >= query_first
            firstintersection!(t.children[i], query, lower, out, filter)
            if out.index > 0
                return
            end
        end
        i += 1
    end

    out.index = 0
    return
end

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

    # TODO: search keys
    i = lower === nothing ? 1 :
        searchsortedfirst(t.entries, notnothing(lower), 1, Int(t.count),
                          Base.ord(basic_isless, identity, false))

    while i <= length(t)
        if intersects(t.entries[i], query) && filter(t.entries[i], query)
            out.index = i
            out.node = t
            return
        elseif last(query) < first(t.keys[i])
            break
        end
        i += 1
    end

    out.index = 0
    return
end

function firstintersection(t::InternalNode{K, V, B},
                           query::Interval{K},
                           lower::Union{Nothing, Interval{K}}) where {K, V, B}
    if isempty(t) || t.maxend < first(query)
        return (nothing, 1)
    end

    i = lower === nothing ? 1 : searchsortedfirst(t.keys, notnothing(lower))

    @inbounds while i <= length(t.children)
        if i > 1 && t.keys[i-1].first > last(query)
            break
        end
        if t.maxends[i] >= first(query)
            w, k = firstintersection(t.children[i], query, lower)
            w !== nothing && return notnothing(w), k
        end
        i += 1
    end

    return (nothing, 1)
end

function firstintersection(t::LeafNode{K, V, B}, query::Interval{K},
                           lower::Union{Nothing, Interval{K}}) where {K, V, B}

    if isempty(t) || t.maxend < first(query)
        return (nothing, 1)
    end

    i = lower === nothing ? 1 : searchsortedfirst(t.keys, notnothing(lower))

    while i <= length(t)
        if intersects(t.keys[i], query) &&
           (lower === nothing || !basic_isless(t.keys[i], notnothing(lower)))
            return (t, i)
        elseif last(query) < first(t.keys[i])
            break
        end
        i += 1
    end

    return (nothing, 1)
end

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

function firstfrom(t::InternalNode{K, V, B}, query::K) where {K, V, B}
    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(t::LeafNode{K, V, B}, query::K) where {K, V, B}
    if isempty(t)
        return (t, 0)
    end

    for i in 1:length(t)
        if last(t.keys[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!(t::LeafNode{K, V, B}, i::Integer,
                           query::AbstractInterval{K},
                           out::Intersection{K, V, B},
                           filter::F) where {F, K, V, B}
    j = i + 1
    while true
        while j <= length(t)
            if intersects(t.keys[j], query) && filter(t.entries[j], query)
                out.index = j
                out.node = t
                return
            end
            j += 1
        end
        j = 1

        t.right === nothing && break

        t = notnothing(t.right)

        if minstart(t) > last(query)
            break
        end

        # When we hit a leaf node with no possible intersections, we start over
        # from the top of the leaf. This is a heuristic to avoid linear search
        # behavior in the presence of extremely long intervals.
        if t.maxend < first(query)
            return firstintersection!(root(t), query,
                                      t.entries[t.count], out, filter)
        end
    end

    out.index = 0
    return
end

mutable struct IntervalIntersectionIterator{F, K, V, B}
    filter::F
    intersection::Intersection{K, V, B}
    t::IntervalBTree{K, V, B}
    query::AbstractInterval{K}

    function IntervalIntersectionIterator{F,K,V,B}(filter, intersection, t, query) where {F,K,V,B}
        return new{F,K,V,B}(filter, intersection, t, query)
    end

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

Base.eltype(::Type{IntervalIntersectionIterator{F,K,V,B}}) where {F,K,V,B} = V
Base.IteratorSize(::Type{IntervalIntersectionIterator{F,K,V,B}}) where {F,K,V,B} = Base.SizeUnknown()

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

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

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

function Base.iterate(it::IntervalIntersectionIterator{F, K, V, B}) where {F, K, V, B}
    firstintersection!(it.t, it.query, nothing, it.intersection, it.filter)
    return iterate(it, nothing)
end

function Base.iterate(it::IntervalIntersectionIterator{F, K, V, B}, _) where {F, K, V, B}
    intersection = it.intersection
    if intersection.index == 0
        return nothing
    end
    entry = intersection.node.entries[intersection.index]
    nextintersection!(intersection.node, intersection.index,
                      it.query, intersection, it.filter)
    return entry, nothing
end

mutable struct IntersectionIterator{F, K, V1, B1, V2, B2}
    filter::F

    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

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

    function IntersectionIterator{F,K,V1,B1,V2,B2}(filter::F, t1, t2, successive::Bool) where {F,K,V1,B1,V2,B2}
        return new{F,K,V1,B1,V2,B2}(filter, t1, t2, successive)
    end
end

Base.eltype(::Type{IntersectionIterator{F,K,V1,B1,V2,B2}}) where {F,K,V1,B1,V2,B2} = Tuple{V1,V2}
Base.IteratorSize(::Type{IntersectionIterator{F,K,V1,B1,V2,B2}}) where {F,K,V1,B1,V2,B2} = Base.SizeUnknown()

function Base.iterate(it::IntersectionIterator, _=iterinitstate(it))
    if it.isdone
        return nothing
    end
    if it.successive
        u = it.u
        w = it.w
        value = (u.entries[it.i], w.entries[it.k])
        @nextleafkey(w, it.w, it.k)
        successive_nextintersection!(it)
    else
        u = it.u
        w = it.w
        value = (u.entries[it.i], w.entries[it.k])
        @nextleafkey(w, it.w, it.k)
        iterative_nextintersection!(it)
    end
    return value, nothing
end

function iterinitstate(it::IntersectionIterator)
    it.isdone = true
    if it.successive
        it.u = isempty(it.t1) ? nothing : firstleaf(it.t1)
        it.w = isempty(it.t2) ? nothing : firstleaf(it.t2)
        it.i, it.k = 1, 1
        successive_nextintersection!(it)
    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) ? nothing : firstleaf(it.t1)
        it.v = isempty(it.t2) ? nothing : firstleaf(it.t2)
        it.w = it.v
        it.i, it.j, it.k = 1, 1, 1
        iterative_nextintersection!(it)
    end
    return nothing  # iterator is stateful
end

function iterative_nextintersection!(
    it::IntersectionIterator{F, K, V1, B1, V2, B2}) where {F, 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 true
        u === nothing && return
        unode = notnothing(u)
        ukey = unode.keys[i]

        v === nothing && return
        vnode = notnothing(v)
        vkey = vnode.keys[j]

        # find next intersection w
        while w !== nothing
            wnode = notnothing(w)
            wkey = wnode.keys[k]

            if first(wkey) <= last(ukey)
                if first(ukey) <= last(wkey) && it.filter(unode.entries[i], wnode.entries[k])
                    it.isdone = false
                    break
                end
                @nextleafkey(wnode, w, k)
            else
                break
            end
        end

        # intersection was found
        if !it.isdone
            break
        end

        # if no intersection found, advance u and start again
        @nextleafkey(unode, u, i)

        u === nothing && break

        # find new v corresponding to new u
        unode = notnothing(u)
        ukey = unode.keys[i]
        while v !== nothing
            vnode = notnothing(v)
            vkey = vnode.keys[j]

            if last(vkey) < first(ukey)
                @nextleafkey(vnode, v, j)
            else
                break
            end
        end

        w, k = v, j
    end

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

function successive_nextintersection!(
    it::IntersectionIterator{F, K, V1, B1, V2, B2}) where {F, K, V1, B1, V2, B2}

    u, w, i, k = it.u, it.w, it.i, it.k

    it.isdone = true

    while true
        u === nothing && return
        unode = notnothing(u)
        ukey = unode.keys[i]

        # find next intersection w
        while w !== nothing
            wnode = notnothing(w)
            wkey = wnode.keys[k]

            if first(wkey) <= last(ukey)
                if first(ukey) <= last(wkey) && it.filter(unode.entries[i], wnode.entries[k])
                    it.isdone = false
                    break
                end

                # next w key
                if k < length(wnode)
                    k += 1
                else
                    k = 1
                    w = wnode.right

                    # When we hit a leaf node with no possible intersections, we
                    # start over from the top of the leaf. This is a heuristic
                    # to avoid linear search behavior in the presence of
                    # extremely long intervals.
                    if w !== nothing && w.maxend < first(ukey)
                        wnode = notnothing(w)
                        w, k = firstintersection(it.t2.root, ukey,
                                                 wnode.keys[wnode.count])
                    end
                end
            else
                break
            end
        end

        # intersection was found
        if !it.isdone
            break
        end

        # if no intersection found, advance u and start again
        @nextleafkey(unode, u, i)

        u === nothing && break

        # find u's first intersection
        unode = notnothing(u)
        ukey = unode.keys[i]

        w, k = firstintersection(it.t2.root, ukey, nothing)
    end

    it.u = u
    it.w = w
    it.i = i
    it.k = k
    return
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(t1::IntervalBTree{K, V1, B1},
                        t2::IntervalBTree{K, V2, B2},
                        filter::F=true_cmp; method=:auto) where {F, K, V1, B1, V2, B2}
    # We decide heuristically which intersection algorithm to use.
    m = length(t1)
    n = length(t2)

    if method == :auto
        iterative_cost  = n + m
        successive_cost = 0.25 * m * log(1 + n) + 1e5
        if iterative_cost < successive_cost
            return IntersectionIterator{F, K, V1, B1, V2, B2}(filter, t1, t2, false)
        else
            return IntersectionIterator{F, K, V1, B1, V2, B2}(filter, t1, t2, true)
        end
    elseif method == :successive
        return IntersectionIterator{F, K, V1, B1, V2, B2}(filter, t1, t2, true)
    elseif method == :iterative
        return IntersectionIterator{F, K, V1, B1, V2, B2}(filter, 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")

end  # module IntervalTrees


text:
# using MLDatasets: MNIST
using Base.Iterators
import Pkg; Pkg.add("Augmentor"); using Augmentor

function test_image_loader()
    augm = CropSize(28,28)
    trn, vld = mk_image_minibatch("../data/flowers",
                4; split=true, fr=0.2,
                balanced=true, shuffle=true,
                train=true,
                aug_pipl=augm, pre_proc=preproc_imagenet)

    lab = get_class_labels(trn)

    return size(first(trn)[1]) == (28,28,3,4) &&
           size(lab) == (3,)
end

function test_image_preload()
    tst = mk_image_minibatch("../data/flowers",
                4; split=false, fr=0.2,
                balanced=false, shuffle=false,
                train=false, pre_load=true,
                aug_pipl=nothing, pre_proc=nothing)

    return size(first(tst)) == (256, 256, 3, 4)
end

function test_image2arr()
        img = load("../data/elecat/cat.jpg")
        x = image2array(img)
        return size(x) == (1024, 1001, 3) && eltype(x) == Float32
end

function test_array2image()
        img1 = array2image(rand(256,256,3))
        img2 = array2image(rand(256,256))
        img3 = array2RGB(rand(256,256))

        return eltype(img1) <: RGB &&
               eltype(img2) <: Gray &&
               eltype(img3) <: RGB
end

# imagenet:
#
function test_preproc_imagenet()

    img = rand(32,32,3)
    img = preproc_imagenet(img)
    return size(img) == (32,32,3)
end

function test_in_classes()
    c = get_imagenet_classes()
    return length(c) == 1000
end