text:
struct LevelDataset
    contents::Vector{String}
end

function get_level(dataset::LevelDataset, n::Int)
    start = 12 * n + 2
    finish = 12 * n + 11
    return dataset.contents[start:finish]
end

"""
The dataset (000.txt) used is [boxoban-levels](https://github.com/deepmind/boxoban-levels/blob/master/medium/train/000.txt) (LICENSE file included along with 000.txt)

The following is the citation of the dataset:

@misc{boxobanlevels,
author = {Arthur Guez, Mehdi Mirza, Karol Gregor, Rishabh Kabra, Sebastien Racaniere, Theophane Weber, David Raposo, Adam Santoro, Laurent Orseau, Tom Eccles, Greg Wayne, David Silver, Timothy Lillicrap, Victor Valdes},
title = {An investigation of Model-free planning: boxoban levels},
howpublished= {https://github.com/deepmind/boxoban-levels/},
year = "2018",
}

The dataset file can be updated by creating suitable [hook](https://github.com/JuliaReinforcementLearning/ReinforcementLearningCore.jl/blob/master/src/core/hooks.jl)
"""
mutable struct SokobanDirected{T, R} <: AbstractGridWorld
    world::GridWorldBase{Tuple{Agent, Wall, Box, Target}}
    agent_pos::CartesianIndex{2}
    agent_dir::AbstractDirection
    reward::T
    rng::R
    dataset::LevelDataset
    box_pos::Vector{CartesianIndex{2}}
    target_pos::Vector{CartesianIndex{2}}
    done::Bool
end

@generate_getters(SokobanDirected)
@generate_setters(SokobanDirected)

mutable struct SokobanUndirected{T, R} <: AbstractGridWorld
    world::GridWorldBase{Tuple{Agent, Wall, Box, Target}}
    agent_pos::CartesianIndex{2}
    reward::T
    rng::R
    dataset::LevelDataset
    box_pos::Vector{CartesianIndex{2}}
    target_pos::Vector{CartesianIndex{2}}
    done::Bool
end

@generate_getters(SokobanUndirected)
@generate_setters(SokobanUndirected)

#####
# Directed
#####

function SokobanDirected(; T = Float32, file = joinpath(dirname(pathof(@__MODULE__)), "envs/sokoban/boxoban-levels/medium/train/000.txt"), rng = Random.GLOBAL_RNG)
    dataset = LevelDataset(readlines(file))
    level = get_level(dataset, 0)
    height = length(level)
    width = length(level[1])

    objects = (AGENT, WALL, BOX, TARGET)
    world = GridWorldBase(objects, height, width)

    agent_pos = CartesianIndex(2, 2)
    world[AGENT, agent_pos] = true

    agent_dir = RIGHT
    reward = zero(T)

    box_pos = CartesianIndex{2}[]
    target_pos = CartesianIndex{2}[]
    done = false

    env = SokobanDirected(world, agent_pos, agent_dir, reward, rng, dataset, box_pos, target_pos, done)

    RLBase.reset!(env)

    return env
end

RLBase.StateStyle(env::SokobanDirected) = RLBase.InternalState{Any}()

RLBase.state_space(env::SokobanDirected, ::RLBase.Observation, ::RLBase.DefaultPlayer) = nothing
const SOKOBAN_DIRECTED_LAYERS = SA.SVector(2, 3, 4)
RLBase.state(env::SokobanDirected, ::RLBase.Observation, ::RLBase.DefaultPlayer) = get_grid(get_world(env), get_agent_pos(env), get_agent_dir(env), get_half_size(env), SOKOBAN_DIRECTED_LAYERS)

RLBase.state_space(env::SokobanDirected, ::RLBase.InternalState, ::RLBase.DefaultPlayer) = nothing
RLBase.state(env::SokobanDirected, ::RLBase.InternalState, ::RLBase.DefaultPlayer) = (copy(get_grid(env)), get_agent_dir(env))

RLBase.action_space(env::SokobanDirected, ::RLBase.DefaultPlayer) = DIRECTED_NAVIGATION_ACTIONS
RLBase.reward(env::SokobanDirected, ::RLBase.DefaultPlayer) = get_reward(env)
RLBase.is_terminated(env::SokobanDirected) = get_done(env)

function RLBase.reset!(env::SokobanDirected{T}) where {T}
    rng = get_rng(env)
    world = get_world(env)

    world[:, :, :] .= false
    env.box_pos = CartesianIndex{2}[]
    env.target_pos = CartesianIndex{2}[]

    level_num = rand(rng, 0:999)
    level = get_level(env.dataset, level_num)

    set_level!(env, level)

    set_agent_dir!(env, RIGHT)

    set_reward!(env, zero(T))
    set_done!(env, false)

    return nothing
end

function (env::SokobanDirected{T})(action::AbstractTurnAction) where {T}
    new_dir = turn(action, get_agent_dir(env))
    set_agent_dir!(env, new_dir)
    world = get_world(env)

    set_done!(env, all(pos -> get_world(env)[TARGET, pos], env.box_pos))
    set_reward!(env, zero(T))

    return nothing
end

function (env::SokobanDirected{T})(action::AbstractMoveAction) where {T}
    world = get_world(env)

    r1 = sum(pos -> world[TARGET, pos], env.box_pos) # can just use get_reward(env) here?
    agent_pos = get_agent_pos(env)
    dir = get_agent_dir(env)
    dest = move(action, dir, agent_pos)
    if !world[WALL, dest]
        beyond_dest = move(action, dir, dest)
        if !world[BOX, dest]
            world[AGENT, agent_pos] = false
            set_agent_pos!(env, dest)
            world[AGENT, dest] = true
        else
            if !world[BOX, beyond_dest] && !world[WALL, beyond_dest]
                world[BOX, dest] = false
                world[BOX, beyond_dest] = true

                idx = findfirst(pos -> pos == dest, env.box_pos)
                env.box_pos[idx] = beyond_dest
                world[AGENT, agent_pos] = false
                set_agent_pos!(env, dest)
                world[AGENT, dest] = true
            end
        end
    end

    set_done!(env, all(pos -> get_world(env)[TARGET, pos], env.box_pos))
    r2 = sum(pos -> world[TARGET, pos], env.box_pos)
    set_reward!(env, convert(T, r2 - r1))

    return nothing
end

function set_level!(env::SokobanDirected, level::Vector{String})
    world = get_world(env)
    for i in 1:length(level), j in 1:length(level[1])
        pos = CartesianIndex(i, j)
        char = level[i][j]

        if char === '@'
            set_agent_pos!(env, pos)
            world[AGENT, pos] = true
            set_agent_dir!(env, RIGHT)
        elseif char === '#'
            world[WALL, pos] = true
        elseif char === '$'
            push!(env.box_pos, pos)
            world[BOX, pos] = true
        elseif char === '.'
            push!(env.target_pos, pos)
            world[TARGET, pos] = true
        end
    end

    return nothing
end

#####
# Undirected
#####

function SokobanUndirected(; T = Float32, file = joinpath(dirname(pathof(@__MODULE__)), "envs/sokoban/boxoban-levels/medium/train/000.txt"), rng = Random.GLOBAL_RNG)
    dataset = LevelDataset(readlines(file))
    level = get_level(dataset, 0)
    height = length(level)
    width = length(level[1])

    objects = (AGENT, WALL, BOX, TARGET)
    world = GridWorldBase(objects, height, width)

    agent_pos = CartesianIndex(2, 2)
    world[AGENT, agent_pos] = true

    reward = zero(T)

    box_pos = CartesianIndex{2}[]
    target_pos = CartesianIndex{2}[]
    done = false

    env = SokobanUndirected(world, agent_pos, reward, rng, dataset, box_pos, target_pos, done)

    RLBase.reset!(env)

    return env
end

RLBase.StateStyle(env::SokobanUndirected) = RLBase.InternalState{Any}()

RLBase.state_space(env::SokobanUndirected, ::RLBase.Observation, ::RLBase.DefaultPlayer) = nothing
const SOKOBAN_UNDIRECTED_LAYERS = SA.SVector(2, 3, 4)
RLBase.state(env::SokobanUndirected, ::RLBase.Observation, ::RLBase.DefaultPlayer) = get_grid(get_world(env), get_agent_pos(env), get_half_size(env), SOKOBAN_UNDIRECTED_LAYERS)

RLBase.state_space(env::SokobanUndirected, ::RLBase.InternalState, ::RLBase.DefaultPlayer) = nothing
RLBase.state(env::SokobanUndirected, ::RLBase.InternalState, ::RLBase.DefaultPlayer) = copy(get_grid(env))

RLBase.action_space(env::SokobanUndirected, player::RLBase.DefaultPlayer) = UNDIRECTED_NAVIGATION_ACTIONS
RLBase.reward(env::SokobanUndirected, ::RLBase.DefaultPlayer) = get_reward(env)
RLBase.is_terminated(env::SokobanUndirected) = get_done(env)

function RLBase.reset!(env::SokobanUndirected{T}) where {T}
    rng = get_rng(env)
    world = get_world(env)

    world[:, :, :] .= false
    env.box_pos = CartesianIndex{2}[]
    env.target_pos = CartesianIndex{2}[]

    level_num = rand(rng, 0:999)
    level = get_level(env.dataset, level_num)

    set_level!(env, level)

    set_reward!(env, zero(T))
    set_done!(env, false)

    return nothing
end

function (env::SokobanUndirected{T})(action::AbstractMoveAction) where {T}
    world = get_world(env)

    r1 = sum(pos -> world[TARGET, pos], env.box_pos) # can just use get_reward(env) here?
    agent_pos = get_agent_pos(env)
    dest = move(action, agent_pos)
    if !world[WALL, dest]
        beyond_dest = move(action, dest)
        if !world[BOX, dest]
            world[AGENT, agent_pos] = false
            set_agent_pos!(env, dest)
            world[AGENT, dest] = true
        else
            if !world[BOX, beyond_dest] && !world[WALL, beyond_dest]
                world[BOX, dest] = false
                world[BOX, beyond_dest] = true

                idx = findfirst(pos -> pos == dest, env.box_pos)
                env.box_pos[idx] = beyond_dest
                world[AGENT, agent_pos] = false
                set_agent_pos!(env, dest)
                world[AGENT, dest] = true
            end
        end
    end

    set_done!(env, all(pos -> get_world(env)[TARGET, pos], env.box_pos))
    r2 = sum(pos -> world[TARGET, pos], env.box_pos)
    set_reward!(env, convert(T, r2 - r1))

    return nothing
end

function set_level!(env::SokobanUndirected, level::Vector{String})
    world = get_world(env)
    for i in 1:length(level), j in 1:length(level[1])
        pos = CartesianIndex(i, j)
        char = level[i][j]

        if char === '@'
            set_agent_pos!(env, pos)
            world[AGENT, pos] = true
        elseif char === '#'
            world[WALL, pos] = true
        elseif char === '$'
            push!(env.box_pos, pos)
            world[BOX, pos] = true
        elseif char === '.'
            push!(env.target_pos, pos)
            world[TARGET, pos] = true
        end
    end

    return nothing
end


text:
Sensor = Class(100)
sensor = Object(Sensor)
Action = Class(120)
action = Object(Action)
act = DFunction([Sensor, Sensor, Action], Action)
apply(act, [sensor, sensor, action])

Sensor = Class(100)
sensor = Object(Sensor)
Action = Class(120)
action = Object(Action)
act = DFunction([Sensor, Action], Action)
tree = Tree(act, Dict(), [Tree(sensor, Dict(), []), Tree(action, Dict(), [])])
init_tree!(tree)
eval_tree!(tree).value[:value] - apply(act, [sensor, action]).value

Sensor = Class(100)
sensor = Object(Sensor)
Action = Class(120)
action = Object(Action)
act = DFunction([Sensor, Action], Action)
tree = Tree(act, Dict(), [Tree(sensor, Dict(), []), Tree(action, Dict(), [])])
init_tree!(tree)
d = ones(length(eval_tree!(tree).value[:value])) - eval_tree!(tree).value[:value]
for i in 1:5000
    tree.value[:d] .= ones(length(eval_tree!(tree).value[:value])) - eval_tree!(tree).value[:value]
    bp_tree!(tree)
end
maximum(abs(ones(length(eval_tree!(tree).value[:value])) - eval_tree!(tree).value[:value]))

Sensor = Class(100)
sensor = Object(Sensor)
Action = Class(120)
action = Object(Action)
act = DFunction([Sensor], Action)
invact = DFunction([Action], Sensor)
axiom = Axiom([1, [2, 3]], 3)

@time train!(axiom, [invact, act, sensor], 10000)
@time foreach(i -> eval_tree!(axiom.tree1), 1:10000)
maximum(abs(eval_tree!(axiom.tree2).value[:value] - eval_tree!(axiom.tree1).value[:value]))

Profile.clear()
@profile train!(axiom, [invact, act, sensor], 1000)
## Profile.print()
Profile.print(format = :flat)

ext: jl

max_stars_repo_name: Non-Contradiction/NextGeneration

max_issues_repo_name: Non-Contradiction/NextGeneration

max_forks_repo_name: Non-Contradiction/NextGeneration


text:
### Calculate the line-of-sight (nautical) miles for each flight and the total travelled by each passenger
reload("nmr")

# this parses records in AComp_Passenger_data.cs into JSON, handling errors
nmr.runjob(nmr.NMR(4, ["AComp_Passenger_data.csv"], nmr.mapper_parserecordacomp, "acomplineofsight.json"))

# this gives us list of flights with [flightid, [srcairport, dstairport]]
nmr.runjob(nmr.NMR(4, ["acomplineofsight.json"], nmr.reducer_lineofsightpassenger, "lineofsightpassenger.json"))

# combine results
c = nmr.runcombiner("lineofsightpassenger.json", nmr.combiner_parsejson)

# sanity checks
@test typeof(c) == Vector{String}
@test length(c) == 17

# we need to join the above with airport coordinates to calculate line of sight for each flight
# we must then create a mr job to calculate passengers and their flights with src and dst airport
# .. to calculate total miles travelled by each passenger

# join the two output files
nmr.runjob(nmr.NMR(4, ["acomplineofsight.json", "lineofsightpassenger.json"], nmr.joiner_lineofsight, "joinedlineofsight.json"))

max_stars_repo_name: deyandyankov/nmr

max_issues_repo_name: deyandyankov/nmr

max_forks_repo_name: deyandyankov/nmr


text:
struct SpectralTransform{NF<:AbstractFloat}
    # SIZE OF SPECTRAL GRID
    trunc::Int      # Spectral truncation
    mx::Int         # Number of zonal wavenumbers
    nx::Int         # Number of total wavenumbers

    # FFT plans
    # rfft_plan::FFTW.rFFTWPlan
    # irfft_plan::AbstractFFTs.ScaledPlan

    # LEGENDRE ARRAYS
    leg_weight::Array{NF,1}          # Legendre weights
    nsh2::Array{Int,1}              # What's this?
    leg_poly::Array{Complex{NF},3}   # Legendre polynomials

    # HARMONIC AND BIHARMONIC DIFFUSION
    ∇²::Array{NF,2}          # Laplacian = l*(l+1)/(R_earth^2)
    ∇⁻²::Array{NF,2}         # inverse Laplacian
    ∇⁴::Array{NF,2}          # Laplacian squared, for biharmonic diffusion

    # Quantities required by functions grad, uvspec, and vds
    gradx::Array{NF,1}
    uvdx::Array{NF,2}
    uvdym::Array{NF,2}
    uvdyp::Array{NF,2}
    gradym::Array{NF,2}
    gradyp::Array{NF,2}
    vddym::Array{NF,2}
    vddyp::Array{NF,2}
end

struct GeoSpectral{NF<:AbstractFloat}
    geometry::Geometry{NF}
    spectral::SpectralTransform{NF}
end

function GeoSpectral(P::Parameters)
    G = Geometry(P)
    S = SpectralTransform(P,G)
    return GeoSpectral{P.NF}(G,S)
end

function SpectralTransform(P::Parameters,G::Geometry)

    @unpack nlon, nlat, nlon_half, nlat_half = G
    @unpack coslat_NH = G
    @unpack R_earth, trunc = P

    # SIZE OF SPECTRAL GRID
    mx = trunc+1
    nx = trunc+1

    # PLAN THE FFTs
    # rfft_plan = plan_rfft(rand(NF,nlon))
    # irfft_plan = plan_irfft(rand(Complex{NF},nlon_half+1),nlon)

    # LEGENDRE WEIGHTS from pole to equator (=first half or array)
    leg_weight = gausslegendre(nlat)[2][1:nlat_half]

    # Spectral packing of speedy is m',n', where m' = m, n' = m+n with m,n being the
    # conventional wavenumbers. Due to Julia's 1-based indexing subtract two, as
    # the Legendre polynomials start with
    nsh2 = zeros(Int, nx)
    for n in 1:nx
        for m in 1:mx
            if m + n - 2 <= trunc + 1
                nsh2[n] = nsh2[n] + 1
            end
        end
    end

    # Epsilon-factors for the recurrence relation of the associated Legendre polynomials.
    ε,ε⁻¹ = ε_recurrence(mx,nx)

    # Generate associated Legendre polynomials
    # get_legendre_poly computes the polynomials at a particular latitiude
    leg_poly = zeros(mx, nx, nlat_half)
    for j in 1:nlat_half
        leg_poly[:,:,j] = legendre_polynomials(j,ε,ε⁻¹,mx,nx,G)
    end

    # leg_poly = zeros(mx, nx, nlat_half)
    # for j in 1:nlat_half
    #     leg_poly[:,:,j] = AssociatedLegendrePolynomials.λlm(0:mx-1,0:nx-1,coslat_NH[j])
    # end

    # LAPLACIANS for harmonic & biharmonic diffusion
    ∇²,∇⁻²,∇⁴ = Laplacians(mx,nx,R_earth)

    gradx   = zeros(mx)          #TODO what's this?
    uvdx    = zeros(mx, nx)
    uvdym   = zeros(mx, nx)
    uvdyp   = zeros(mx, nx)
    gradym  = zeros(mx, nx)
    gradyp  = zeros(mx, nx)
    vddym   = zeros(mx, nx)
    vddyp   = zeros(mx, nx)

    for m in 1:mx
        for n in 1:nx
            m1  = m - 1
            m2  = m1 + 1
            el1 = m + n - 2
            if n == 1
                gradx[m]   = m1/R_earth
                uvdx[m,1]  = -R_earth/(m1 + 1)
                uvdym[m,1] = 0.0
                vddym[m,1] = 0.0
            else
                uvdx[m,n]   = -R_earth*m1/(el1*(el1 + 1.0))
                gradym[m,n] = (el1 - 1.0)*ε[m2,n]/R_earth
                uvdym[m,n]  = -R_earth*ε[m2,n]/el1
                vddym[m,n]  = (el1 + 1.0)*ε[m2,n]/R_earth
            end
            gradyp[m,n] = (el1 + 2.0)*ε[m2,n+1]/R_earth
            uvdyp[m,n]  = -R_earth*ε[m2,n+1]/(el1 + 1.0)
            vddyp[m,n]  = el1*ε[m2,n+1]/R_earth
        end
    end

    SpectralTransform{P.NF}(    trunc,mx,nx,
                                # rfft_plan,irfft_plan,
                                leg_weight,nsh2,leg_poly,∇²,∇⁻²,∇⁴,
                                gradx,uvdx,uvdym,uvdyp,gradym,gradyp,vddym,vddyp)
end

"""
Laplacian operator in spectral space via element-wise matrix-matrix multiplication.
"""
function ∇²(A::Array{Complex{NF},2},
            G::GeoSpectral{NF}) where {NF<:AbstractFloat}
    return -G.spectral.∇².*A
end

"""
In-place version of ∇².
"""
function ∇²!(   Out::Array{Complex{NF},2},
                In::Array{Complex{NF},2},
                G::GeoSpectral{NF}) where {NF<:AbstractFloat}

    @unpack ∇² = G.spectral

    mx,nx = size(In)
    @boundscheck (mx,nx) == size(Out) || throw(BoundsError())
    @boundscheck (mx,nx) == size(∇²) || throw(BoundsError())

    for n in 1:nx
        for m in 1:mx
            @inbounds Out[m,n] = -∇²[m,n]*In[m,n]
        end
    end
end

"""
Inverse Laplacian in spectral space via element-wise matrix-matrix multiplication.
"""
function ∇⁻²(   A::Array{Complex{NF},2},
                G::GeoSpectral{NF}) where {NF<:AbstractFloat}
    return -G.spectral.∇⁻².*A
end

"""
Transform a spectral array into grid-point space.
"""
function gridded(  input::Array{Complex{NF},2},
                   G::GeoSpectral{NF}) where {NF<:AbstractFloat}
    return fourier_inverse(legendre_inverse(input,G),G)
end

"""
Transform a gridded array into spectral space.
"""
function spectral(  input::Array{NF,2},
                    G::GeoSpectral{NF}) where {NF<:AbstractFloat}
    return legendre(fourier(input,G),G)
end

function grad!( ψ::Array{NF,2},
                psdx::Array{Complex{NF},2},
                psdy::Array{NF,2},
                G::GeoSpectral{NF}) where {NF<:AbstractFloat}

    #TODO boundscheck

    @unpack trunc, mx, nx = G.spectral
    @unpack gradx, gradyp, gradym = G.spectral

    for n in 1:nx
        psdx[:,n] = gradx.*ψ[:,n]*im
    end

    for m in 1:mx
        psdy[m,1]  =  gradyp[m,1]*ψ[m,2]
        psdy[m,nx] = -gradym[m,nx]*ψ[m,trunc+1]
    end

    for n in 2:trunc+1
        for m in 1:mx
            psdy[m,n] = -gradym[m,n]*ψ[m,n-1] + gradyp[m,n]*ψ[m,n+1]
        end
    end
end

function vds!(  ucosm::Array{NF,2},
                vcosm::Array{NF,2},
                vorm::Array{NF,2},
                divm::Array{NF,2},
                G::GeoSpectral{NF}) where {NF<:AbstractFloat}

    #TODO boundscheck

    @unpack trunc, mx, nx = G.spectral
    @unpack gradx, vddym, vddyp = G.spectral

    #TODO preallocate in a diagnosticvars struct
    zp = zeros(Complex{NF}, mx,nx)
    zc = zeros(Complex{NF}, mx,nx)

    for n in 1:nx
        zp[:,n] = gradx.*ucosm[:,n]*im
        zc[:,n] = gradx.*vcosm[:,n]*im
    end

    for m in 1:mx
        #TODO this has an implicit conversion to complex{NF}, issue?
        vorm[m,1]  = zc[m,1] - vddyp[m,1]*ucosm[m,2]
        vorm[m,nx] = vddym[m,nx]*ucosm[m,trunc+1]
        divm[m,1]  = zp[m,1] + vddyp[m,1]*vcosm[m,2]
        divm[m,nx] = -vddym[m,nx]*vcosm[m,trunc+1]
    end

    for n in 2:trunc+1
        for m in 1:mx
            #TODO same here
            vorm[m,n] =  vddym[m,n]*ucosm[m,n-1] - vddyp[m,n]*ucosm[m,n+1] + zc[m,n]
            divm[m,n] = -vddym[m,n]*vcosm[m,n-1] + vddyp[m,n]*vcosm[m,n+1] + zp[m,n]
        end
    end
end

function uvspec!(   vorm::Array{NF,2},
                    divm::Array{NF,2},
                    ucosm::Array{NF,2},
                    vcosm::Array{NF,2},
                    G::GeoSpectral{NF}) where {NF<:AbstractFloat}

    #TODO boundscheck

    @unpack trunc, mx, nx = G.spectral
    @unpack uvdx, uvdyp, uvdym = G.spectral

    #TODO preallocate elsewhere
    zp = uvdx.*vorm*im
    zc = uvdx.*divm*im

    for m in 1:mx
        ucosm[m,1]  =  zc[m,1] - uvdyp[m,1]*vorm[m,2]
        ucosm[m,nx] =  uvdym[m,nx]*vorm[m,trunc+1]
        vcosm[m,1]  =  zp[m,1] + uvdyp[m,1]*divm[m,2]
        vcosm[m,nx] = -uvdym[m,nx]*divm[m,trunc+1]
    end

    for n in 2:trunc+1
        for m in 1:mx
          vcosm[m,n] = -uvdym[m,n]*divm[m,n-1] + uvdyp[m,n]*divm[m,n+1] + zp[m,n]
          ucosm[m,n] =  uvdym[m,n]*vorm[m,n-1] - uvdyp[m,n]*vorm[m,n+1] + zc[m,n]
        end
    end
end

function vdspec!(   ug::Array{NF,2},
                    vg::Array{NF,2},
                    vorm::Array{NF,2},
                    divm::Array{NF,2},
                    kcos::Bool,
                    G::GeoSpectral{NF}) where {NF<:AbstractFloat}

    #TODO boundscheck

    @unpack nlat, nlon, cosgr, cosgr2 = G.geometry

    #TODO preallocate elsewhere
    ug1 = zeros(NF, nlon, nlat)
    vg1 = zeros(NF, nlon, nlat)

    # either cosgr or cosgr2
    cosgr = kcos ? cosgr : cosgr2

    for j in 1:nlat
        for i in 1:nlon
            ug1[i,j] = ug[i,j]*cosgr[j]
            vg1[i,j] = vg[i,j]*cosgr[j]
        end
    end

    #TODO add spectral_trans and geometry as arguments
    specu = spectral(ug1,G)
    specv = spectral(vg1,G)
    vds!(specu, specv, vorm, divm)
end

"""Truncate spectral field by seting the spectral coefficients of the lower right triangle to zero. """
function spectral_truncation!(  A::AbstractArray{Complex{NF},2},    # spectral field to be truncated
                                trunc::Int                          # truncate to total wave number `trunc`
                                ) where {NF<:AbstractFloat}         # number format NF
    m,n = size(A)
    zeero = zero(Complex{NF})   # zero element of elements in A

    @inbounds for j in 1:n
        for i in 1:m
            if i+j-2 > trunc    # if total wavenumber larger than trunc
                A[i,j] = zeero
            end
        end
    end
end

"""Spectral truncation with unpacking Geospectral struct."""
function spectral_truncation!(  A::AbstractArray{Complex{NF},2},
                                G::GeoSpectral{NF}
                                ) where NF
    spectral_truncation!(A,G.spectral.trunc)    # unpack GeoSpectral struct
end

"""Spectral truncation of a grid-point field with memory allocation."""
function spectral_truncation(   input::AbstractArray{NF,2},
                                G::GeoSpectral{NF}
                                ) where NF
    input_spectral = spectral(input,G)          # allocates memory
    spectral_truncation!(input_spectral,G)      # in-place truncation

    # allocates memory to return spectrally truncated gridded field
    return gridded(input_spectral, G)
end

"""In-place version of spectral trunction of a grid-point field."""
function spectral_truncation!(  input::AbstractArray{NF,2},
                                input_spectral::AbstractArray{Complex{NF},2},
                                G::GeoSpectral{NF}
                                ) where NF
    spectral!(input_spectral,input,G)       # in-place spectral transform from input to input_spectral
    spectral_truncation!(input_spectral,G)  # in-place truncation
    gridded!(input,input_spectral,G)        # in-place backtransform
end

function spectral_truncation!(  input::Array{NF,2},
                                G::GeoSpectral{NF}
                                ) where NF
    spectral_truncation!(input,spectral(input,G),G)
end

"""
Computes Laplacian matrix-operators (for element-wise multiplication).
Laplacian, Laplacian squared and inverse Laplacian in spectral space
correspond to a multiplication with the total wavenumber:

    ∇²_n^m = N*(N+1)/R_earth^2

with N = m+n-2 the total wavenumber -2 due to 1-based indexing. The biharmonic
operator is ∇⁴ = (∇²)², the inverse Laplacian is ∇⁻² = 1 ./ ∇².
"""
function Laplacians(mx::Integer,nx::Integer,R_earth::Real)
    ∇²   = zeros(mx, nx)
    ∇⁻²  = zeros(mx, nx)
    ∇⁴   = zeros(mx, nx)

    for n in 1:nx
        for m in 1:mx
            # total wavenumber is m+n, -2 due to Julia's 1-based indexing
            N = m+n-2
            ∇²[m,n] = N*(N+1)/R_earth^2
            ∇⁴[m,n] = ∇²[m,n]^2
        end
    end

    # inverse Laplacian, the first coefficient being zero corresponds
    # to some (?, TODO) boundary conditions
    ∇⁻²[1,1] = 0.0                  # don't divide by 0
    ∇⁻²[2:end] = 1 ./ ∇⁻²[2:end]    # all other elements in matrix

    return ∇²,∇⁻²,∇⁴
end


text:
module ImageStack

using Defer, FileIO, Images, LinearAlgebra, VideoIO

export Frames, frames, load
export Point, Box, Circle, FreeForm
export width, height, area, box

include("imagestack/dataset.jl")
include("imagestack/accumulator.jl")
include("imagestack/features.jl")

end


text: # TODO: Test HyperLogLog


text:
# Creates a ground constraint

using SparseArrays

include("../constraints/constraints.jl")

"""
    makeGroundConstraint(NSteps::Int64, nDim::Int64, index = 2)

Creates an [`AL_AffineInequality`](@ref) constraint such that the rocket stays
above the ground. (Constraint has form Ax - b ≤ 0). Ignores the last point
in the trajectory.

Usage: `NSteps` states the number of timesteps in the trajectory. `nDim` states
the dimensionality of the problem (e.g. 1D, 2D, or 3D). Index describes which
direction is up (e.g. index = 2 means "y" is up).
"""
function makeGroundConstraint(NSteps::Int64, nDim::Int64, index = 2)
    # Safegaurd to assue the input is reasonable
    if index > nDim
        print("Invalid index. Changed from $index to ")
        index = min(nDim, index)
        println(index)
    end

    # Build up the sparse matrices
    sizeTraj = 3 * NSteps * nDim + 2 * nDim
    AMat = spzeros(NSteps, sizeTraj)
    bVec = spzeros(NSteps)

    # AMat = zeros(NSteps, sizeTraj)
    # bVec = zeros(NSteps)

    # Iterate through and select each index that needs to be flipped to -1
    currInd = index
    currRow = 1
    while currRow <= NSteps
        AMat[currRow, currInd] = -1

        currRow += 1
        currInd += 3 * nDim
    end

    # return  the result as a constraint struct
    return AL_AffineInequality(AMat, bVec)

end

max_stars_repo_name: danineamati/TrajOptSOCPs


text:
function _update_till_convergence!( pcs::PotentialCalculationSetup{T, S, 3},
                                    convergence_limit,
                                    via_KernelAbstractions::Bool;
                                    n_iterations_between_checks = 500,
                                    depletion_handling::Val{depletion_handling_enabled} = Val{false}(),
                                    only2d::Val{only_2d} = Val{false}(),
                                    is_weighting_potential::Val{_is_weighting_potential} = Val{false}(),
                                    use_nthreads::Int = Base.Threads.nthreads(),
                                    max_n_iterations::Int = 10_000, # -1
                                    verbose::Bool = true
                                ) where {T, S, depletion_handling_enabled, only_2d, _is_weighting_potential, DAT} #  <: GPUArrays.AbstractGPUArray
    device = KernelAbstractions.get_device(pcs.potential)
    ndrange = size(pcs.potential)[1:3] .- 2
    kernel = get_sor_kernel(S, device, Val(via_KernelAbstractions))
    c_limit = _is_weighting_potential ? convergence_limit : abs(convergence_limit * pcs.bias_voltage)
    c = (one(c_limit) + c_limit) * 10 # Has to be larger than c_limit at the beginning
    n_performed_iterations = 0
    tmp_potential = similar(pcs.potential, ndrange)
    inner_ranges = broadcast(i -> 2:size(tmp_potential, i) + 1, (1, 2, 3))
    is_logging(io) = isa(io, Base.TTY) == false || (get(ENV, "CI", nothing) == "true")
    cs = fill(c, 4) # 4 is chosen by testing
    if verbose prog = ProgressThresh(c_limit; dt = 0.1, desc = "Convergence: ", output = stderr, enabled = !is_logging(stderr)) end
    while c > c_limit
        for _ in 1:n_iterations_between_checks-1
            update!(pcs, kernel, ndrange; use_nthreads, depletion_handling, is_weighting_potential, only2d)
            n_performed_iterations += 1
        end
        begin
            tmp_potential[:, :, :] .= view(pcs.potential, inner_ranges..., 1)
            update!(pcs, kernel, ndrange; use_nthreads, depletion_handling, is_weighting_potential, only2d)
            tmp_potential[:, :, :] .-= view(pcs.potential, inner_ranges..., 1)
            n_performed_iterations += 1
            c = maximum(abs.(tmp_potential))
            if verbose ProgressMeter.update!(prog, c) end
            cs = circshift(cs, -1)
            cs[end] = c
            cs_μ = mean(cs)
            cs_σ = std(cs, mean = cs_μ)
            if cs_σ < c_limit
                # Convergence limit not reached but the value of c does not change anymore
                # Especially needed in case of undepleted detectors (grid points switching between depleted and undepleted)
                break
            end
        end
        if n_performed_iterations >= max_n_iterations break end
    end
    if verbose ProgressMeter.finish!(prog) end
    return c
end


text:
using CSV
using DataFrames
using DelimitedFiles
using YAML
using Query

"""
    share_labor!(d::Dict, set::Dict)
"""
function share_labor!(d::Dict, set::Dict)
    println("  Calculating labor share")
    !(:region in collect(keys(d))) && share_region!(d, set)
    _share_lshr0!(d, set)

    cols = [:yr,:r,:s,:value]
    df = copy(d[:gsp])
    df[!,:value] .= df[:,:cmp] ./ df[:,:comp]
    df = dropnan(df[:,cols])

    df = SLiDE._join_to_operate(df, d[:region], d[:lshr0])
    df = edit_with(df, Rename.([:x1,:x2,:x3], [:value,:region,:lshr0]))

    # Use the national average labor share (calculated when partitioning)
    # in cases where the labor share is zero (e.g. banking, finance, etc.).
    ii = .&(df[:,:value] .== 0.0, df[:,:region] .> 0.0)
    df[ii,:value] .= df[ii,:lshr0]
    df = fill_zero(df)

    # Save labor to input in functions to calculate annual and regional average wage shares.
    d[:labor] = df[:,cols]
    d[:labor_temp] = copy(d[:labor])

    df[!,:wg] .= _share_wg!(d, set)[:,:value]
    df[!,:hw] .= _share_hw!(d, set)[:,:value]
    df[!,:avg_wg] = _share_avg_wg!(d)[:,:value]
    df[!,:sec_labor] = _share_sec_labor!(d)[:,:value]

    df[!,:labor_temp] .= df[:,:value]

    wg = .&(df[:,:wg], df[:,:avg_wg] .!== 0.0)
    df[wg,:value] .= df[wg,:avg_wg]
    df[df[:,:hw],:value] .= df[df[:,:hw],:sec_labor]
    d[:labor] = dropnan(df[:,cols])
    d[:labor_calc] = df
end

"""
    share_region!(d::Dict, set::Dict)
`region`: Regional share of value added
"""
function share_region!(d::Dict, set::Dict)
    println("  Calculating regional share of value added")
    :gdpcat in propertynames(d[:gsp]) && _share_gsp!(d, set)

    cols = [:yr,:r,:s,:value]
    df = edit_with(copy(d[:gsp]), Rename(:gdp,:value))[:,cols]
    df = df / transform_over(df, :r)

    # Let the used and scrap sectors be an average of other sectors.
    # These are the only sectors that have NaN values.
    df = df[.!isnan.(df[:,:value]),:]
    df_s  = combine_over(df, :s)
    df_s /= transform_over(df_s, :r)

    df_s = crossjoin(DataFrame(s = set[:oth,:use]), df_s)

    d[:region] = dropnan(dropmissing(sort([df; df_s])))
    verify_over(d[:region],:r) !== true && @error("Regional shares don't sum to 1.")
    return d[:region]
end

"`gsp`: Calculated gross state product."
function _share_gsp!(d::Dict, set::Dict)
    df = copy(d[:gsp])
    :gdpcat in propertynames(d[:gsp]) && (df = unstack(dropzero(d[:gsp]), :gdpcat, :value))

    # df = edit_with(df, Replace.(Symbol.(set[:gdpcat]),missing,0.0))
    df = edit_with(df, Replace.(find_oftype(df, AbstractFloat), missing, 0.0))

    df[!,:sudo] .= df[:,:gdp] - df[:,:taxsbd]
    df[!,:comp] .= df[:,:cmp] + df[:,:gos]

    df[!,:calc] .= df[:,:cmp] + df[:,:gos] + df[:,:taxsbd]
    df[!,:diff] .= df[:,:calc] - df[:,:gdp]

    d[:gsp] = df
end

"`lshr0`: Labor share of value added"
function _share_lshr0!(d::Dict, set::Dict)
    va0 = edit_with(unstack(copy(d[:va0]), :va, :value),
        [Replace.(Symbol.(set[:va]), missing, 0.0); Drop(:units,"all","==")])

    d[:lshr0]  = va0[:,[:yr,:s,:compen]]
    d[:lshr0] /= (va0[:,[:yr,:s,:compen]] + va0[:,[:yr,:s,:surplus]])

    # !!!!! _partition_lshr0 needs to come after calibration.
    # Order is: io, calibrate, share, disagg.
    d[:lshr0][va0[:,:surplus] .< 0,:value] .= 1.0
    dropmissing!(d[:lshr0])
end

# " `netval`: Factor totals"
# function _share_netval!(d::Dict)
#     cols = [:yr,:r,:s,:sudo,:comp]
#     df = copy(d[:gsp])
#     df[!,:sudo] .= df[:,:gdp] - df[:,:taxsbd]
#     df[!,:comp] .= df[:,:cmp] + df[:,:gos]
#     d[:netval] = df[:,cols]
# end

"`wg(yr,r,s)`: Index pairs with high wage shares (>1)"
function _share_wg!(d::Dict, set::Dict)
    d[:wg] = copy(d[:labor])
    d[:wg][!,:value] .= d[:wg][:,:value] .> 1
    set[:wg] = values.(eachrow(d[:wg][d[:wg][:,:value], find_oftype(d[:wg],Not(Bool))]))
    return d[:wg] # must return the DataFrame
end

"`hw(r,s)`: (region,sector) pairings with ALL wage shares > 1"
function _share_hw!(d::Dict, set::Dict)
    !(:wg in collect(keys(d))) && _share_wg!(d, set)

    df = copy(d[:wg])
    d[:hw] = combine_over(df, :yr; fun = prod)
    set[:hw] = values.(eachrow(d[:hw][d[:hw][:,:value], find_oftype(d[:hw],String)]))
    return transform_over(df, :yr; fun = prod) # must return the DataFrame
end

function _share_avg_wg!(d::Dict)
    # Here, WiNDC uses .!wg, which DOES include labor = 1.
    !(:wg in collect(keys(d))) && _share_wg!(d, set)
    not_wg = edit_with(copy(d[:wg]), Replace(:value,Bool,Not))

    df = copy(d[:labor]) * not_wg
    d[:avg_wg] = edit_with(combine_over(df, :yr) / combine_over(not_wg, :yr), Replace(:value,NaN,0.0))
    return edit_with(transform_over(df, :yr) / transform_over(not_wg, :yr), Replace(:value,NaN,0.0))
end

function _share_sec_labor!(d::Dict)
    # Cannot use .!wg because this will include labor = 1.
    not_wg = copy(d[:labor])
    not_wg[!,:value] .= not_wg[:,:value] .< 1

    df = copy(d[:labor]) * not_wg
    d[:sec_labor] = combine_over(df, :r) / combine_over(not_wg, :r)
    return transform_over(df, :r) / transform_over(not_wg, :r)
end

max_stars_repo_name: matbesancon/SLiDE

max_issues_repo_name: matbesancon/SLiDE

max_forks_repo_name: matbesancon/SLiDE


text:
using Test
using QuantumOpticsBase
using LinearAlgebra, SparseArrays, Random

mutable struct test_operators{BL<:Basis,BR<:Basis} <: AbstractOperator{BL,BR}
  basis_l::BL
  basis_r::BR
  data::Matrix{ComplexF64}
  test_operators(b1::Basis, b2::Basis, data) = length(b1) == size(data, 1) && length(b2) == size(data, 2) ? new{typeof(b1),typeof(b2)}(b1, b2, data) : throw(DimensionMismatch())
end

@testset "operators" begin

Random.seed!(0)

b1 = GenericBasis(3)
b2 = GenericBasis(2)
b = b1 ⊗ b2
op1 = randoperator(b1)
op2 = randoperator(b2)
op = randoperator(b, b)
op_test = test_operators(b, b, op.data)
op_test2 = test_operators(b1, b, randoperator(b1, b).data)
op_test3 = test_operators(b1 ⊗ b2, b2 ⊗ b1, randoperator(b, b).data)
ψ = randstate(b)
ρ = randoperator(b)

@test basis(op1) == b1
@test length(op1) == length(op1.data) == length(b1)^2

@test_throws ArgumentError op_test*op_test
@test_throws ArgumentError -op_test

@test_throws ArgumentError 1 + op_test
@test_throws ArgumentError op_test + 1
@test_throws ArgumentError 1 - op_test
@test_throws ArgumentError op_test - 1

@test_throws ArgumentError dagger(op_test)
@test_throws ArgumentError op_test'
@test_throws ArgumentError identityoperator(test_operators, b, b)
@test_throws ArgumentError tr(op_test)
@test_throws ArgumentError ptrace(op_test, [1])
@test_throws ArgumentError ishermitian(op_test)
@test_throws MethodError dense(op_test)
@test_throws ArgumentError sparse(op_test)
@test_throws ArgumentError transpose(op_test)

@test expect(1, op1, ρ) ≈ expect(embed(b, 1, op1), ρ)
@test expect(1, op1, ψ) ≈ expect(embed(b, 1, op1), ψ)
@test expect(op, [ρ, ρ]) == [expect(op, ρ) for i=1:2]
@test expect(1, op1, [ρ, ψ]) == [expect(1, op1, ρ), expect(1, op1, ψ)]

@test variance(1, op1, ρ) ≈ variance(embed(b, 1, op1), ρ)
@test variance(1, op1, ψ) ≈ variance(embed(b, 1, op1), ψ)
@test variance(op, [ρ, ρ]) == [variance(op, ρ) for i=1:2]
@test variance(1, op1, [ρ, ψ]) == [variance(1, op1, ρ), variance(1, op1, ψ)]

@test tensor(op_test) === op_test
@test_throws ArgumentError tensor(op_test, op_test)
@test_throws ArgumentError permutesystems(op_test, [1, 2])

@test embed(b, b, [1,2], op) == embed(b, [1,2], op)
@test embed(b, Dict{Vector{Int}, SparseOpType}()) == identityoperator(b)
@test_throws QuantumOpticsBase.IncompatibleBases embed(b1⊗b2, [2], [op1])

b_comp = b⊗b
@test isapprox(embed(b_comp, [1,[3,4]], [op1,op]), dense(op1 ⊗ one(b2) ⊗ op))
@test isapprox(embed(b_comp, [[1,2],4], [op,op2]), dense(op ⊗ one(b1) ⊗ op2))
@test_throws QuantumOpticsBase.IncompatibleBases embed(b_comp, [[1,2],3], [op,op2])
@test_throws QuantumOpticsBase.IncompatibleBases embed(b_comp, [[1,3],4], [op,op2])

function basis_vec(n, N)
    x = zeros(Complex{Float64}, N)
    x[n+1] = 1
    return x
end
function basis_maker(dims...)
    function bm(ns...)
        bases = [basis_vec(n, dim) for (n, dim) in zip(ns, dims)][end:-1:1]
        return reduce(kron, bases)
    end
end

embed_op = embed(b_comp, [1,4], op)
bv = basis_maker(3,2,3,2)
all_idxs = [(idx, jdx) for (idx, jdx) in [Iterators.product(0:1, 0:2)...]]

m11 = reshape([Bra(b_comp, bv(0,idx,jdx,0)) * embed_op * Ket(b_comp, bv(0,kdx,ldx,0))
              for ((idx, jdx), (kdx, ldx)) in Iterators.product(all_idxs, all_idxs)], (6,6))
@test isapprox(m11 / op.data[1, 1], diagm(0=>ones(Complex{Float64}, 6)))

m21 = reshape([Bra(b_comp, bv(1,idx,jdx,0)) * embed_op * Ket(b_comp, bv(0,kdx,ldx,0))
              for ((idx, jdx), (kdx, ldx)) in Iterators.product(all_idxs, all_idxs)], (6,6))
@test isapprox(m21 / op.data[2,1], diagm(0=>ones(Complex{Float64}, 6)))

m12 = reshape([Bra(b_comp, bv(0,idx,jdx,0)) * embed_op * Ket(b_comp, bv(1,kdx,ldx,0))
              for ((idx, jdx), (kdx, ldx)) in Iterators.product(all_idxs, all_idxs)], (6,6))
@test isapprox(m12 / op.data[1,2], diagm(0=>ones(Complex{Float64}, 6)))

b_comp = b_comp⊗b_comp
OP_test1 = dense(tensor([op1,one(b2),op,one(b1),one(b2),op1,one(b2)]...))
OP_test2 = embed(b_comp, [1,[3,4],7], [op1,op,op1])
@test isapprox(OP_test1, OP_test2)

b8 = b2⊗b2⊗b2
cnot = [1 0 0 0; 0 1 0 0; 0 0 0 1; 0 0 1 0]
op_cnot = DenseOperator(b2⊗b2, cnot)
OP_cnot = embed(b8, [1,3], op_cnot)
@test ptrace(OP_cnot, [2])/2. == op_cnot
@test_throws AssertionError embed(b2⊗b2, [1,1], op_cnot)

@test_throws ArgumentError exp(sparse(op1))

@test one(b1).data == Diagonal(ones(b1.shape[1]))
@test one(op1).data == Diagonal(ones(b1.shape[1]))

@test_throws ArgumentError conj!(op_test)

bfock = FockBasis(2)
bspin = SpinBasis(1//2)
b = bfock⊗bspin⊗bspin

op12 = destroy(bfock)⊗sigmap(bspin)

@test embed(b, [1,2], op12) == destroy(bfock)⊗sigmap(bspin)⊗one(bspin)
@test embed(b, [1,3], op12) == destroy(bfock)⊗one(bspin)⊗sigmap(bspin)

end # testset


text:
import Base.Ordering
import Base.Forward
import Base.Reverse
import DataStructures.eq
import Base.lt
import Base.ForwardOrdering
import DataStructures.IntSemiToken

## Function fulldump dumps the entire tree; helpful for debugging.

function fulldump(t::DataStructures.BalancedTree23)
    thislevstack = Int[]
    rl = t.rootloc
    dpth = t.depth
    push!(thislevstack, rl)
    println("  rootloc = $rl depth = $dpth")
    nextlevstack = Int[]
    levcount = 0
    for mydepth = 1 : dpth
        isleaf = (dpth == mydepth)
        sz = size(thislevstack,1)
        if sz == 0
            break
        end
        levcount += 1
        println("-------------\n------LEVEL $levcount")
        for i = 1 : sz
            ii = thislevstack[i]
            p = t.tree[ii].parent
            sk1 = t.tree[ii].splitkey1
            sk2 = t.tree[ii].splitkey2
            c1 = t.tree[ii].child1
            c2 = t.tree[ii].child2
            c3 = t.tree[ii].child3
            dt = (mydepth == dpth)? "child(data)" : "child(tree)"
            if c3 == 0
                println("ii = $ii splitkey1 = /$sk1/ $dt.1 = $c1 $dt.2 = $c2 parent = $p")
                if !isleaf
                    push!(nextlevstack,c1)
                    push!(nextlevstack,c2)
                end
            else
                println("ii = $ii splitkey1 = /$sk1/ splitkey2 = /$sk2/ $dt.1 = $c1 $dt.2 = $c2 $dt.3 =$c3 parent = $p")
                if !isleaf
                    push!(nextlevstack,c1)
                    push!(nextlevstack,c2)
                    push!(nextlevstack,c3)
                end
            end
        end
        thislevstack = nextlevstack
        nextlevstack = Int[]
    end
    println("----- FREE TREE CELLS -----")
    for i = 1 : size(t.freetreeinds,1)
        println(t.freetreeinds[i])
    end

    println("----- DATA ARRAY ---")
    for j = 1 : size(t.data,1)
        k = t.data[j].k
        d = t.data[j].d
        p = t.data[j].parent
        if j > 2
            println("j = $j k = /$k/ d = /$d/ parent = /$p/")
        else
            println("j = $j (  k = /$k/ d = /$d/ ) parent = /$p/")
        end
    end
    println("----- FREE DATA CELLS -----")
    for i = 1 : size(t.freedatainds,1)
        println(t.freedatainds[i])
    end
end

push_test!(m, k, v) = push!(m, k=>v)

## Function checkcorrectness checks a balanced tree for correctness.

function checkcorrectness{K,D,Ord <: Ordering}(t::DataStructures.BalancedTree23{K,D,Ord},
                                               allowdups=false)
    dsz = size(t.data, 1)
    tsz = size(t.tree, 1)
    r = t.rootloc
    bfstreenodes = Vector{Int}(0)
    tdpth = t.depth
    intree = IntSet()
    levstart = Vector{Int}(tdpth)
    push!(bfstreenodes, r)
    levstart[1] = 1
    push!(intree, r)
    for curdepth = 2 : tdpth
        levstart[curdepth] = size(bfstreenodes, 1) + 1
        for l = levstart[curdepth - 1] : levstart[curdepth] - 1
            anc = bfstreenodes[l]
            c1 = t.tree[anc].child1
            if in(c1, intree)
                println("anc = $anc c1 = $c1")
                throw(ErrorException("Tree contains loops 1"))
            end
            push!(bfstreenodes, c1)
            push!(intree, c1)
            c2 = t.tree[anc].child2
            if in(c2, intree)
                throw(ErrorException("Tree contains loops 2"))
            end
            push!(bfstreenodes, c2)
            push!(intree, c2)
            c3 = t.tree[anc].child3
            if c3 > 0
                if in(c3, intree)
                    throw(ErrorException("Tree contains loops 3"))
                end
                push!(bfstreenodes, c3)
                push!(intree, c3)
            end
        end
    end
    bfstreesize = size(bfstreenodes, 1)
    dataused = IntSet()
    minkeys = Vector{K}(bfstreesize)
    maxkeys = Vector{K}(bfstreesize)
    for s = levstart[tdpth] : bfstreesize
        anc = bfstreenodes[s]
        c1 = t.tree[anc].child1
        if s == levstart[tdpth]
            if c1 != 1
                throw(ErrorException("Leftmost data descendant should be node 1"))
            end
        else
            minkeys[s] = t.data[c1].k
        end
        c2 = t.tree[anc].child2
        c3 = t.tree[anc].child3
        lastchild = c3 > 0? c3 : c2
        if s == bfstreesize
            if lastchild != 2
                throw(ErrorException("Rightmost data descendant should be node 2"))
            end
        else
            maxkeys[s] = t.data[lastchild].k
        end
        if s > levstart[tdpth] &&
            (lt(t.ord, minkeys[s], maxkeys[s - 1]) ||
             (!lt(t.ord, maxkeys[s-1],minkeys[s]) && !allowdups))
            println("tdpth = ", tdpth, " s = ", s,
                    " maxkeys[s-1] = ", maxkeys[s-1],
                    " minkeys[s] = ", minkeys[s])
            throw(ErrorException("Data nodes out of order"))
        end
        if s < bfstreesize || c3 > 0
            if !eq(t.ord, t.tree[anc].splitkey1, t.data[c2].k)
                throw(ErrorException("Splitkey1 of leaf should match key of 2nd child"))
            end
        end
        if s < bfstreesize && c3 > 0
            if !eq(t.ord, t.tree[anc].splitkey2, t.data[c3].k)
                throw(ErrorException("Splitkey2 of leaf should match key of 1st child"))
            end
        end
        if t.data[c1].parent != anc || t.data[c2].parent != anc ||
            (c3 > 0 && t.data[c3].parent != anc)
            println("c1 = $c1 c2 = $c2 c3 = $c3 anc = $anc")
            println("t.data[c1].parent = $(t.data[c1].parent) t.data[c2].parent = $(t.data[c2].parent)")
            if c3 > 0
                println("t.data[c3].parent = $(t.data[c3].parent)")
            end
            throw(ErrorException("Incorrect parent node for data child"))
        end
        push!(dataused, c1)
        push!(dataused, c2)
        if c3 > 0
            push!(dataused, c3)
        end
    end
    for curdepth = tdpth - 1 : -1 : 1
        cp = levstart[curdepth + 1]
        for s = levstart[curdepth] : levstart[curdepth + 1] - 1
            anc = bfstreenodes[s]
            c1 = t.tree[anc].child1
            @assert(c1 == bfstreenodes[cp])
            if s > levstart[curdepth]
                mk1 = minkeys[cp]
            end
            cp += 1
            if t.tree[c1].parent != anc
                throw(ErrorException("Parent/child1 links do not match"))
            end
            c2 = t.tree[anc].child2
            @assert(c2 == bfstreenodes[cp])
            mk2 = minkeys[cp]
            cp += 1
            if t.tree[c2].parent != anc
                throw(ErrorException("Parent/child2 links do not match"))
            end
            c3 = t.tree[anc].child3
            @assert(s == levstart[curdepth] ||
                    lt(t.ord,mk1,mk2) || (!lt(t.ord,mk2,mk1) && allowdups))
            if c3 > 0
                if t.tree[c3].parent != anc
                    throw(ErrorException("Parent/child3 links do not match"))
                end
                mk3 = minkeys[cp]
                cp += 1
                @assert(lt(t.ord,mk2, mk3) ||
                        !lt(t.ord,mk3,mk2) && allowdups)
            end
            if s > levstart[curdepth]
                minkeys[s] = mk1
            end
            if !eq(t.ord, t.tree[anc].splitkey1, mk2)
                throw(ErrorException("Minkey2 not equal to minimum key among descendants of child2"))
            end
            if c3 > 0 && !eq(t.ord, t.tree[anc].splitkey2, mk3)
                throw(ErrorException("Minkey3 not equal to minimum key among descendants of child3"))
            end
        end
    end
    freedata = IntSet()
    for i = 1 : size(t.freedatainds,1)
        fdi = t.freedatainds[i]
        if in(fdi, freedata)
            throw(ErrorException("t.freedatainds has repeated element $i"))
        end
        if fdi < 1 || fdi > dsz
            throw(ErrorException("t.freedatainds entry out of range"))
        end
        push!(freedata, fdi)
    end
    if last(t.useddatacells) > dsz
        throw(ErrorException("t.useddatacells has indices larger than t.data size"))
    end
    for i = 1 : dsz
        if (in(i, dataused) && !in(i, t.useddatacells)) ||
            (!in(i,dataused) && in(i, t.useddatacells))
            throw(ErrorException("Mismatch between actual data cells used and useddatacells array"))
        end
        if (in(i, freedata) && in(i, dataused)) ||
            (!in(i,freedata) && !in(i, dataused))
            throw(ErrorException("Mismatch between t.freedatainds and t.useddatacells"))
        end
    end
    freetree = IntSet()
    for i = 1 : size(t.freetreeinds,1)
        tfi = t.freetreeinds[i]
        if in(tfi, freetree)
            throw(ErrorException("Free tree index repeated twice"))
        end
        if tfi < 1 || tfi > tsz
            throw(ErrorException("Free tree index out of range"))
        end
        push!(freetree, tfi)
    end
    for i = 1 : tsz
        if (!in(i, intree) && !in(i, freetree)) ||
            (in(i, intree) && in(i, freetree))
            throw(ErrorException("Mismatch between t.freetreeinds and actual cells used"))
        end
    end
end

function test1()
    # a few basic tests of SortedDict to start
    m1 = SortedDict((Dict{Compat.ASCIIString,Compat.ASCIIString}()), Forward)
    kdarray = ["hello", "jello", "alpha", "beta", "fortune", "random",
               "july", "wednesday"]
    checkcorrectness(m1.bt, false)
    for i = 1 : div(size(kdarray,1), 2)
        k = kdarray[i*2-1]
        d = kdarray[i*2]
        #println("- inserting: k = $k d = $d")
        m1[k] = d
        checkcorrectness(m1.bt, false)
        # fulldump(m1.bt)
    end
    i1 = startof(m1)
    count = 0
    while i1 != pastendsemitoken(m1)
        count += 1
        k,d = deref((m1,i1))
        #println("+ reading: k = $k, d = $d")
        i1 = advance((m1,i1))
        checkcorrectness(m1.bt, false)
    end
    @test count == 4
    nothing
end

function test2()
    # test all methods of SortedDict here except loops
    m0 = SortedDict(Dict{Int, Float64}())
    m1 = SortedDict(Dict(8=>32.0, 12=>33.1, 6=>18.2))
    expected = ([6,8,12], [18.2, 32.0, 33.1])
    checkcorrectness(m1.bt, false)
    ii = startof(m1)
    m2 = packdeepcopy(m1)
    m3 = packcopy(m1)
    p = first(m1)
    @test p[1] == 6 && p[2] == 18.2
    @test in(Pair(8,32.0),m3)
    @test !in(Pair(8,32.1),m3)
    push_test!(m1, 12, 33.2)
    checkcorrectness(m1.bt, false)
    @test length(m1) == 3
    @test m1[12] == 33.2
    push_test!(m1, 12, 33.1)
    @test length(m1) == 3
    @test m1[12] == 33.1
    for j = 1 : 3
        @test ii != pastendsemitoken(m1)
        pr = deref((m1,ii))
        @test pr[1] == expected[1][j] && pr[2] == expected[2][j]
        checkcorrectness(m1.bt, false)
        oldii = ii
        ii = advance((m1,ii))
        delete!((m1,oldii))
    end
    checkcorrectness(m1.bt, false)
    checkcorrectness(m2.bt, false)
    @test length(m2) == 3
    ii = startof(m2)
    for j = 1 : 3
        pr = deref((m2,ii))
        @test pr[1] == expected[1][j] && pr[2] == expected[2][j]
        ii = advance((m2,ii))
    end

    checkcorrectness(m3.bt, false)
    @test length(m3) == 3
    ii = startof(m3)
    for j = 1 : 3
        pr = deref((m3,ii))
        @test pr[1] == expected[1][j] && pr[2] == expected[2][j]
        ii = advance((m3,ii))
    end

    @test isempty(m1)
    @test length(m1) == 0
    N = 5000
    for i = N : -1 : 2
        m1[i] = convert(Float64,i) ^ 2
        if i % 50 == 0
            checkcorrectness(m1.bt, false)
        end
    end
    @test !isempty(m1)
    assert(length(m1) == N - 1)
    for i = 2 : N
        d = pop!(m1, i)
        @test d == convert(Float64,i)^2
        if i % 50 == 0
            checkcorrectness(m1.bt, false)
        end
    end
    @test isempty(m1)
    @test length(m1) == 0
    for i = N : -1 : 2
        m1[i] = convert(Float64,i) ^ 2
        if i % 50 == 0
            checkcorrectness(m1.bt, false)
        end
    end
    ii = endof(m1)
    for i = 1 : N - 1
        pr = deref((m1,ii))
        @test pr[1] == N + 1 - i && pr[2] == convert(Float64,pr[1]) ^ 2
        ii = regress((m1,ii))
    end
    lastprime = 1
    while true
        ii = searchsortedafter(m1, lastprime)
        if ii == pastendsemitoken(m1)
            break
        end
        j = deref_key((m1,ii))
        for k = j * 2 : j : N
            p = find(m1, k)
            if p != pastendsemitoken(m1)
                delete!((m1,p))
            end
        end
        lastprime = j
    end
    checkcorrectness(m1.bt, false)
    @test ii == pastendsemitoken(m1)
    @test status((m1,ii)) == 3
    @test status((m1,SDSemiToken(-1))) == 0
    t = 0
    u = 0.0
    for pr in m1
        t += pr[1]
        u += pr[2]
    end
    numprimes = length(m1)
    pn = primes(N)
    @test t == sum(pn)
    @test u == sum(pn.^2)
    ij = endof(m1)
    @test deref_key((m1,ij)) == last(pn) &&
       convert(Float64, last(pn)^2) ==  deref_value((m1,ij))
    m1[6] = 49.0
    @test length(m1) == numprimes + 1
    @test m1[6] == 49.0
    b, i6 = insert!(m1, 6, 50.0)
    @test length(m1) == numprimes + 1
    @test !b
    p = deref((m1,i6))
    @test p[1] == 6 && p[2] == 50.0
    m1[i6] = 9.0
    p = deref((m1,i6))
    @test p[1] == 6 && p[2] == 9.0
    @test m1[i6] == 9.0
    b2, i7 = insert!(m1, 8, 51.0)
    @test b2
    @test length(m1) == numprimes + 2
    p = deref((m1,i7))
    @test p[1] == 8 && p[2] == 51.0
    delete!((m1,i7))
    z = pop!(m1, 6)
    checkcorrectness(m1.bt, false)
    @test z == 9.0
    i8 = startof(m1)
    p = deref((m1,i8))
    @test p[1] == 2 && p[2] == 4.0
    @test i8 != beforestartsemitoken(m1)
    @test status((m1,i8)) == 1
    i9 = regress((m1,i8))
    @test i9 == beforestartsemitoken(m1)
    @test status((m1,i9)) == 2
    i10 = find(m1, 17)
    i11 = regress((m1,i10))
    @test deref_key((m1,i11)) == 13
    i12 = searchsortedfirst(m1, 47)
    i13 = searchsortedfirst(m1, 48)
    @test deref_key((m1,i12)) == 47
    @test deref_key((m1,i13)) == 53
    i14 = searchsortedafter(m1, 47)
    i15 = searchsortedafter(m1, 48)
    @test deref_key((m1,i14)) == 53
    @test deref_key((m1,i15)) == 53
    i16 = searchsortedlast(m1, 47)
    i17 = searchsortedlast(m1, 48)
    @test deref_key((m1,i16)) == 47
    @test deref_key((m1,i17)) == 47
    ww = primes(N)
    cc = last(m1)
    @test cc[1] == last(ww)
    wwx = first(m1)
    @test wwx[1] == 2
    tpr = eltype(m1)
    @test tpr == Pair{Int,Float64}
    tpr2 = eltype(typeof(m1))
    @test tpr2 == Pair{Int,Float64}
    kt = keytype(m1)
    @test kt == Int
    kt2 = keytype(typeof(m1))
    @test kt2 == Int
    vt = valtype(m1)
    @test vt == Float64
    vt2 = valtype(typeof(m1))
    @test vt2 == Float64
    @test ordtype(m1) == ForwardOrdering
    @test ordtype(typeof(m1)) == ForwardOrdering

    co = orderobject(m1)
    @test co == Forward
    @test haskey(m1, 71)
    @test !haskey(m1, 77)
    @test get(m1, 70, -45.2) == -45.2
    @test get(m1, 83, -45.2) == convert(Float64,83)^2
    h = get!(m1, 5, 27.0)
    @test h == 25.0
    h = get!(m1, 6, 27.0)
    @test h == 27.0
    @test m1[6] == 27.0
    @test length(m1) == length(ww) + 1
    pop!(m1, 6)
    @test getkey(m1,7, 9) == 7
    @test getkey(m1,7.0, 9) == 7
    @test getkey(m1,12, 9) == 9
    delete!(m1, 17)
    @test length(m1) == length(ww) - 1
    @test deref_key((m1,advance((m1,find(m1,13))))) == 19
    empty!(m1)
    checkcorrectness(m1.bt, false)
    @test isempty(m1)
    c1 = SortedDict(Dict("Eggplants"=>3,
                        "Figs"=>9,
                        "Apples"=>7))
    c2 = SortedDict(Dict("Eggplants"=>6,
                        "Honeydews"=>19,
                        "Melons"=>11))
    @test !isequal(c1,c2)
    c3 = merge(c1, c2)
    checkcorrectness(c3.bt, false)
    c4 = SortedDict(Dict("Apples"=>7,
                        "Figs"=>9,
                        "Eggplants"=>6,
                        "Melons"=>11,
                        "Honeydews"=>19))
    @test isequal(c3,c4)
    c5 = SortedDict(Dict("Apples"=>7))
    @test !isequal(c4,c5)
    merge!(c1,c2)
    checkcorrectness(c1.bt, false)
    @test isequal(c3,c1)
    merge!(c3,c3)
    @test isequal(c3,c1)
    checkcorrectness(c3.bt, false)
    nothing
end

function bitreverse(i)
    zeroi = zero(i)
    onei = one(i)
    twoi = onei + onei
    r = zeroi
    for j = 1 : 32
        r *= twoi
        r += (i & onei)
        i = div(i,twoi)
    end
    r
end

function test3{T}(z::T)
    ## Test the loops
    zero1 = zero(z)
    one1 = one(z)
    two1 = one1 + one1
    m1 = SortedDict(Dict{T,T}())
    N = 1000
    for l = 1 : N
        lUi = convert(T, l)
        m1[bitreverse(lUi)] = lUi
    end
    count = 0
    for (stok,k,v) in semitokens(inclusive(m1, startof(m1), endof(m1)))
        for (stok2,k2,v2) in semitokens(exclusive(m1, startof(m1), pastendsemitoken(m1)))
            c = compare(m1,stok,stok2)
            if c < 0
                @test deref_key((m1,stok)) < deref_key((m1,stok2))
            elseif c == 0
                @test deref_key((m1,stok)) == deref_key((m1,stok2))
            else
                @test deref_key((m1,stok)) > deref_key((m1,stok2))
            end
            count += 1
        end
    end
    @test eltype(semitokens(exclusive(m1, startof(m1), pastendsemitoken(m1)))) ==
    Tuple{IntSemiToken, T, T}
    @test count == N^2
    N = 10000
    sk = zero1
    sv = zero1
    skhalf = zero1
    svhalf = zero1
    for l = 1 : N
        lUi = convert(T, l)
        brl = bitreverse(lUi)
        sk += brl
        m1[brl] = lUi
        sv += lUi
        if brl < div(N, 2)
            skhalf += brl
            svhalf += lUi
        end
    end
    count = 0
    sk2 = zero1
    sv2 = zero1
    for (k,v) in m1
        sk2 += k
        sv2 += v
        count += 1
    end
    @test count == N
    @test sk2 == sk
    @test sv == sv2
    sk2 = zero1
    for k in keys(m1)
        sk2 += k
    end
    @test eltype(keys(m1)) == T
    @test sk2 == sk

    sk2a = zero1

    for k in eachindex(m1)
        sk2a += k
    end
    @test eltype(eachindex(m1)) == T
    @test sk2a == sk

    sk2b = zero1
    for st in onlysemitokens(m1)
        sk2b += deref_key((m1,st))
    end
    @test sk2b == sk

    sv2 = zero1
    for v in values(m1)
        sv2 += v
    end
    @test eltype(values(m1)) == T

    @test sv == sv2
    count = 0
    for (st,k) in semitokens(keys(m1))
        @test deref_key((m1,st)) == k
        count += 1
    end
    @test eltype(semitokens(keys(m1))) == Tuple{IntSemiToken, T}

    @test count == N
    count = 0
    for (st,v) in semitokens(values(m1))
        @test deref_value((m1,st)) == v
        count += 1
    end
    @test count == N
    @test eltype(semitokens(values(m1))) == Tuple{IntSemiToken, T}

    pos1 = searchsortedfirst(m1, div(N,2))
    sk2 = zero1
    for k in keys(exclusive(m1, startof(m1), pos1))
        sk2 += k
    end
    @test sk2 == skhalf
    @test eltype(keys(exclusive(m1, startof(m1), pos1))) == T

    sk2a = zero1
    for k in eachindex(exclusive(m1, startof(m1), pos1))
        sk2a += k
    end
    @test sk2a == skhalf
    @test eltype(eachindex(exclusive(m1, startof(m1), pos1))) == T

    sv2 = zero1
    for v in values(exclusive(m1, startof(m1), pos1))
        sv2 += v
    end
    @test sv2 == svhalf
    count = 0
    for (k,v) in exclusive(m1, pastendsemitoken(m1), pastendsemitoken(m1))
        count += 1
    end
    @test eltype(keys(exclusive(m1, startof(m1), pos1))) == T
    @test count == 0

    count = 0
    for (k,v) in inclusive(m1, startof(m1), beforestartsemitoken(m1))
        count += 1
    end
    @test count == 0
    @test eltype(keys(inclusive(m1, startof(m1), beforestartsemitoken(m1)))) ==
       T

    count = 0
    sk5 = zero1
    for k in eachindex(inclusive(m1, startof(m1), startof(m1)))
        sk5 += k
        count += 1
    end
    @test count == 1 && sk5 == deref_key((m1,startof(m1)))
    @test eltype(eachindex(inclusive(m1, startof(m1), startof(m1)))) == T

    factors = SortedMultiDict(Int[], Int[])
    N = 1000
    len = 0
    sum1 = 0
    sum2 = 0
    for factor = 1 : N
        for multiple = factor : factor : N
            insert!(factors, multiple, factor)
            sum1 += multiple
            sum2 += factor
            len += 1
        end
    end

    sum1a = 0
    sum2a = 0

    for (k,v) in factors
        sum1a += k
        sum2a += v
    end
    @test sum1a == sum1 && sum2a == sum2

    sum1a = 0
    sum2a = 0
    for st in eachindex(factors)
        (k,v) = deref((factors,st))
        sum1a += k
        sum2a += v
    end
    @test sum1a == sum1 && sum2a == sum2

    sum1a = 0
    sum2a = 0
    for st in onlysemitokens(factors)
        (k,v) = deref((factors,st))
        sum1a += k
        sum2a += v
    end
    @test sum1a == sum1 && sum2a == sum2
    @test eltype(onlysemitokens(factors)) == IntSemiToken

    sum2 = 0
    for (k,v) in inclusive(factors,
                           searchsortedfirst(factors,70),
                           searchsortedlast(factors,70))
        sum2 += v
    end

    @test sum2 == 1 + 2 + 5 + 7 + 10 + 14 + 35 + 70
    @test eltype(inclusive(factors,
                           searchsortedfirst(factors,70),
                           searchsortedlast(factors,70))) == Pair{Int,Int}

    sum2 = 0
    for st in eachindex(inclusive(factors,
                                  searchsortedfirst(factors,70),
                                  searchsortedlast(factors,70)))
        v = deref_value((factors,st))
        sum2 += v
    end

    @test sum2 == 1 + 2 + 5 + 7 + 10 + 14 + 35 + 70
    @test eltype(eachindex(inclusive(factors,
                                     searchsortedfirst(factors,70),
                                     searchsortedlast(factors,70)))) == IntSemiToken

    sum3 = 0
    for (k,v) in exclusive(factors,
                           searchsortedfirst(factors,60),
                           searchsortedfirst(factors,61))
        sum3 += v
    end
    @test sum3 == 1 + 2 + 3 + 4 + 5 + 6 + 10 + 12 + 15 + 20 + 30 + 60
    @test eltype(exclusive(factors,
                           searchsortedfirst(factors,70),
                           searchsortedlast(factors,70))) == Pair{Int,Int}

    sum3 = 0
    for st in eachindex(exclusive(factors,
                                  searchsortedfirst(factors,60),
                                  searchsortedfirst(factors,61)))
        v = deref_value((factors,st))
        sum3 += v
    end
    @test sum3 == 1 + 2 + 3 + 4 + 5 + 6 + 10 + 12 + 15 + 20 + 30 + 60
    @test eltype(eachindex(exclusive(factors,
                                     searchsortedfirst(factors,70),
                                     searchsortedlast(factors,70)))) == IntSemiToken

    sum4 = 0
    for k in keys(factors)
        sum4 += k
    end
    @test sum4 == sum1
    @test eltype(keys(factors)) == Int

    sum5 = 0
    for v in values(factors)
        sum5 += v
    end

    @test sum5 == sum2a
    @test eltype(values(factors)) == Int

    sum2 = 0
    for k in keys(inclusive(factors,
                            searchsortedfirst(factors,70),
                            searchsortedlast(factors,70)))
        sum2 += k
    end
    @test sum2 == 70 * 8
    @test eltype(keys(inclusive(factors,
                                searchsortedfirst(factors,70),
                                searchsortedlast(factors,70)))) ==  Int

    sum3 = 0
    for k in keys(exclusive(factors,
                            searchsortedfirst(factors,60),
                            searchsortedfirst(factors,61)))
        sum3 += k
    end
    @test sum3 == 60 * 12
    @test eltype(keys(exclusive(factors,
                                searchsortedfirst(factors,60),
                                searchsortedfirst(factors,61)))) == Int

    sum2 = 0
    for v in values(inclusive(factors,
                              searchsortedfirst(factors,70),
                              searchsortedlast(factors,70)))
        sum2 += v
    end
    @test sum2 == 1 + 2 + 5 + 7 + 10 + 14 + 35 + 70
    @test eltype(values(inclusive(factors,
                                  searchsortedfirst(factors,60),
                                  searchsortedfirst(factors,61)))) == Int

    sum3 = 0
    for v in values(exclusive(factors,
                              searchsortedfirst(factors,60),
                              searchsortedfirst(factors,61)))
        sum3 += v
    end
    @test sum3 == 1 + 2 + 3 + 4 + 5 + 6 + 10 + 12 + 15 + 20 + 30 + 60
    @test eltype(values(exclusive(factors,
                                  searchsortedfirst(factors,60),
                                  searchsortedfirst(factors,61)))) == Int

    sum1b = 0
    sum2b = 0
    for (st,k,v) in semitokens(factors)
        @test deref_value((factors,st)) == v
        sum1b += k
        sum2b += v
    end
    @test sum1b == sum1a && sum2b == sum2a
    @test eltype(semitokens(factors)) == Tuple{IntSemiToken, Int, Int}

    sum2 = 0
    for (st,k,v) in semitokens(inclusive(factors,
                                         searchsortedfirst(factors,70),
                                         searchsortedlast(factors,70)))
        @test deref_value((factors,st)) == v
        sum2 += v
    end
    @test sum2 == 1 + 2 + 5 + 7 + 10 + 14 + 35 + 70
    @test eltype(semitokens(inclusive(factors,
                                         searchsortedfirst(factors,70),
                                         searchsortedlast(factors,70)))) ==
        Tuple{IntSemiToken, Int, Int}

    sum3 = 0
    for (st,k,v) in semitokens(exclusive(factors,
                                         searchsortedfirst(factors,60),
                                         searchsortedfirst(factors,61)))
        @test deref_value((factors,st)) == v
        sum3 += v
    end
    @test sum3 == 1 + 2 + 3 + 4 + 5 + 6 + 10 + 12 + 15 + 20 + 30 + 60
    @test eltype(semitokens(exclusive(factors,
                                      searchsortedfirst(factors,60),
                                      searchsortedfirst(factors,61)))) ==
         Tuple{IntSemiToken, Int, Int}

    sum4 = 0
    for (st,k) in semitokens(keys(factors))
        @test deref_key((factors,st)) == k && mod(k,deref_value((factors,st))) == 0
        sum4 += k
    end
    @test sum4 == sum1
    @test eltype(semitokens(keys(factors))) == Tuple{IntSemiToken,Int}

    sum5 = 0
    for (st,v) in semitokens(values(factors))
        @test deref_value((factors,st)) == v
        sum5 += v
    end
    @test sum5 == sum2a
    @test eltype(semitokens(values(factors))) == Tuple{IntSemiToken,Int}

    sum2 = 0
    for (st,k) in semitokens(keys(inclusive(factors,
                                           searchsortedfirst(factors,70),
                                           searchsortedlast(factors,70))))
        @test deref_key((factors,st)) == k && mod(k,deref_value((factors,st))) == 0
        sum2 += k
    end
    @test sum2 == 70 * 8
    @test eltype(semitokens(keys(inclusive(factors,
                                           searchsortedfirst(factors,70),
                                           searchsortedlast(factors,70))))) ==
        Tuple{IntSemiToken, Int}

    sum3 = 0
    for (st,k) in semitokens(keys(inclusive(factors,
                                            searchsortedfirst(factors,60),
                                            searchsortedlast(factors,60))))
        @test deref_key((factors,st)) == k && mod(k,deref_value((factors,st))) == 0
        sum3 += k
    end
    @test sum3 == 60 * 12
    @test eltype(semitokens(keys(inclusive(factors,
                                            searchsortedfirst(factors,60),
                                            searchsortedlast(factors,60))))) ==
        Tuple{IntSemiToken, Int}

    sum2 = 0
    for (st,v) in semitokens(values(inclusive(factors,
                                              searchsortedfirst(factors,70),
                                              searchsortedlast(factors,70))))
        @test deref_value((factors,st)) == v
        sum2 += v
    end
    @test sum2 == 1 + 2 + 5 + 7 + 10 + 14 + 35 + 70
    @test eltype(semitokens(values(inclusive(factors,
                                              searchsortedfirst(factors,70),
                                              searchsortedlast(factors,70))))) ==
      Tuple{IntSemiToken, Int}

    sum3 = 0
    for (st,v) in semitokens(values(exclusive(factors,
                                              searchsortedfirst(factors,60),
                                              searchsortedfirst(factors,61))))
        @test deref_value((factors,st)) == v
        sum3 += v
    end
    @test sum3 == 1 + 2 + 3 + 4 + 5 + 6 + 10 + 12 + 15 + 20 + 30 + 60
    @test eltype(semitokens(values(exclusive(factors,
                                              searchsortedfirst(factors,60),
                                              searchsortedfirst(factors,61))))) ==
       Tuple{IntSemiToken, Int}

    s = SortedSet([39, 24, 2, 14, 45, 107, 66])
    sum1 = 0
    for k in s
        sum1 += k
    end
    @test sum1 == sum([39, 24, 2, 14, 45, 107, 66])

    sum1 = 0
    for (st,k) in semitokens(s)
        @test deref((s,st)) == k
        sum1 += k
    end
    @test sum1 == sum([39, 24, 2, 14, 45, 107, 66])
    @test eltype(semitokens(s)) == Tuple{IntSemiToken, Int}

    sum1 = 0
    for st in onlysemitokens(s)
        k = deref((s,st))
        sum1 += k
    end
    @test sum1 == sum([39, 24, 2, 14, 45, 107, 66])
    @test eltype(onlysemitokens(s)) == IntSemiToken

    sum1 = 0
    for st in eachindex(s)
        k = deref((s,st))
        sum1 += k
    end
    @test sum1 == sum([39, 24, 2, 14, 45, 107, 66])
    @test eltype(eachindex(s)) == IntSemiToken

    sum2 = 0
    for k in inclusive(s,
                       searchsortedfirst(s, 24),
                       searchsortedfirst(s, 66))
        sum2 += k
    end
    @test sum2 == 24 + 39 + 45 + 66
    @test eltype(inclusive(s,
                       searchsortedfirst(s, 24),
                       searchsortedfirst(s, 66))) == Int

    sum2 = 0
    for st in eachindex(inclusive(s,
                                  searchsortedfirst(s, 24),
                                  searchsortedfirst(s, 66)))
        sum2 += deref((s,st))
    end
    @test sum2 == 24 + 39 + 45 + 66
    @test eltype(eachindex(inclusive(s,
                                     searchsortedfirst(s, 24),
                                     searchsortedfirst(s, 66)))) == IntSemiToken

    sum2 = 0
    for (st,k) in semitokens(inclusive(s,
                                       searchsortedfirst(s, 24),
                                       searchsortedfirst(s, 66)))
        @test deref((s,st)) == k
        sum2 += k
    end
    @test sum2 == 24 + 39 + 45 + 66
    @test eltype(semitokens(inclusive(s,
                                       searchsortedfirst(s, 24),
                                       searchsortedfirst(s, 66)))) ==
      Tuple{IntSemiToken, Int}

    sum2 = 0
    for st in onlysemitokens(inclusive(s,
                                       searchsortedfirst(s, 24),
                                       searchsortedfirst(s, 66)))
        sum2 += deref((s,st))
    end
    @test sum2 == 24 + 39 + 45 + 66
    @test eltype(onlysemitokens(inclusive(s,
                                          searchsortedfirst(s, 24),
                                          searchsortedfirst(s, 66)))) == IntSemiToken

    sum3 = 0
    for k in exclusive(s,
                       searchsortedfirst(s, 24),
                       searchsortedfirst(s, 66))
        sum3 += k
    end
    @test sum3 == 24 + 39 + 45
    @test eltype(exclusive(s,
                       searchsortedfirst(s, 24),
                       searchsortedfirst(s, 66))) == Int

    sum3 = 0
    for (st,k) in semitokens(exclusive(s,
                                       searchsortedfirst(s, 24),
                                       searchsortedfirst(s, 66)))
        @test deref((s,st)) == k
        sum3 += k
    end
    @test sum3 == 24 + 39 + 45
    @test eltype(semitokens(exclusive(s,
                                       searchsortedfirst(s, 24),
                                       searchsortedfirst(s, 66)))) ==
     Tuple{IntSemiToken, Int}

    sum3 = 0
    for st in eachindex(exclusive(s,
                                  searchsortedfirst(s, 24),
                                  searchsortedfirst(s, 66)))
        sum3 += deref((s,st))
    end
    @test sum3 == 24 + 39 + 45
    @test eltype(eachindex(exclusive(s,
                                     searchsortedfirst(s, 24),
                                     searchsortedfirst(s, 66)))) == IntSemiToken
    nothing
end

function test4()
    # test all the errors of sorted containers
    m = SortedDict(Dict("a" => 6, "bb" => 9))
    @test_throws KeyError println(m["b"])
    m2 = SortedDict(Dict{Compat.ASCIIString, Int}())
    @test_throws BoundsError println(first(m2))
    @test_throws BoundsError println(last(m2))
    state1 = start(m2)
    @test_throws BoundsError next(m2, state1)
    i1 = find(m,"a")
    delete!((m,i1))
    i2 = find(m,"bb")
    @test_throws BoundsError start(inclusive(m,i1,i2))
    @test_throws BoundsError start(exclusive(m,i1,i2))
    @test_throws KeyError delete!(m,"a")
    @test_throws KeyError pop!(m,"a")
    m3 = SortedDict((Dict{Compat.ASCIIString, Int}()), Reverse)
    @test_throws ArgumentError isequal(m2, m3)
    @test_throws BoundsError m[i1]
    @test_throws BoundsError regress((m,beforestartsemitoken(m)))
    @test_throws BoundsError advance((m,pastendsemitoken(m)))
    m1 = SortedMultiDict(Int[], Int[])
    @test_throws ArgumentError m3 = SortedMultiDict(["a", "b"], [1,2,3])
    @test_throws ArgumentError isequal(SortedMultiDict(["a"],[1]), SortedMultiDict(["b"], [1.0]))
    @test_throws ArgumentError isequal(SortedMultiDict(["a"],[1],Reverse), SortedMultiDict(["b"], [1]))
    @test_throws BoundsError first(m1)
    @test_throws BoundsError last(m1)
    s = SortedSet([3,5])
    @test_throws KeyError delete!(s,7)
    @test_throws KeyError pop!(s, 7)
    pop!(s)
    pop!(s)
    @test_throws BoundsError pop!(s)
    @test_throws BoundsError first(s)
    @test_throws BoundsError last(s)
    @test_throws ArgumentError isequal(SortedSet(["a"]), SortedSet([1]))
    @test_throws ArgumentError isequal(SortedSet(["a"]), SortedSet(["b"],Reverse))
    @test_throws ArgumentError (("a",6) in m)
    @test_throws ArgumentError ((2,5) in m1)
    nothing
end

function seekfile(fname)
    #fullname = joinpath(Pkg.dir("DataStructures"), "test", fname)
    fname
end

immutable CaseInsensitive <: Ordering
end

lt(::CaseInsensitive, a, b) = isless(lowercase(a), lowercase(b))
eq(::CaseInsensitive, a, b) = isequal(lowercase(a), lowercase(b))

function test5()
    ## Test use of alternative orderings in test5
    keylist = ["Apple", "aPPle", "berry", "CHerry", "Dairy", "diary"]
    vallist = [6,9,-4,2,1,8]
    m = SortedDict(Dict{Compat.ASCIIString,Int}())
    for j = 1:6
        m[keylist[j]] = vallist[j]
    end
    checkcorrectness(m.bt, false)
    expectedord1 = [1,4,5,2,3,6]
    count = 0
    for p in m
        count += 1
        @test p[1] == keylist[expectedord1[count]] &&
                p[2] == vallist[expectedord1[count]]
    end
    @test count == 6
    m2 = SortedDict((Dict{Compat.ASCIIString, Int}()), Reverse)
    for j = 1 : 6
        m2[keylist[j]] = vallist[j]
    end
    checkcorrectness(m2.bt, false)
    expectedord2 = [6,3,2,5,4,1]
    count = 0
    for p in m2
        count += 1
        @test p[1] == keylist[expectedord2[count]] &&
                p[2] == vallist[expectedord2[count]]
    end
    @test count == 6
    m3 = SortedDict((Dict{Compat.ASCIIString, Int}()), CaseInsensitive())
    for j = 1 : 6
        m3[keylist[j]] = vallist[j]
    end
    @test "BERRY" in keys(m3)
    @test !("BERRY" in collect(keys(m3)))
    checkcorrectness(m3.bt, false)
    expectedord3 = [2,3,4,5,6]
    count = 0
    for p in m3
        count += 1
        @test p[1] == keylist[expectedord3[count]] &&
                p[2] == vallist[expectedord3[count]]
    end
    @test count == 5
    m3empty = similar(m3)
    @test eltype(m3empty) == Pair{Compat.ASCIIString, Int} &&
       orderobject(m3empty) == CaseInsensitive() &&
       length(m3empty) == 0 && ordtype(m3empty) == CaseInsensitive &&
       ordtype(typeof(m3empty)) == CaseInsensitive
    m4 = SortedDict((Dict{Compat.ASCIIString,Int}()), Lt((x,y) -> isless(lowercase(x),lowercase(y))))
    for j = 1 : 6
        m4[keylist[j]] = vallist[j]
    end
    checkcorrectness(m4.bt, false)
    count = 0
    for p in m4
        count += 1
        @test p[1] == keylist[expectedord3[count]] &&
                p[2] == vallist[expectedord3[count]]
    end
    @test count == 5
    nothing
end

## test6, test6a and test6b are not run by package
## testing but are included for timing tests.
## They are identical except for their usage of ordering objects.

function test6(numtrial::Int, expectedk::String, expectedd::String)
    NSTRINGPAIR = 50000
    m1 = SortedDict(Dict{Compat.ASCIIString,Compat.ASCIIString}())
    strlist = Compat.ASCIIString[]
    open(seekfile("wordsScram.txt"), "r") do inio
        for j = 1 : NSTRINGPAIR * 2
            push!(strlist, chomp(readline(inio)))
        end
    end
    for trial = 1 : numtrial
        empty!(m1)
        count = 1
        for j = 1 : NSTRINGPAIR
            k = strlist[count]
            d = strlist[count + 1]
            m1[k] = d
            count += 2
        end
        delct = div(NSTRINGPAIR, 6)
        for j = 1 : delct
            l = find(m1, strlist[j * 2 + 1])
            delete!((m1,l))
        end
        spec = div(NSTRINGPAIR * 3,4)
        l = startof(m1)
        for j = 1 : spec
            l = advance((m1,l))
        end
        k,d = deref((m1,l))
        ekn = expectedk
        edn = expectedd
        @test k == ekn && d == edn
    end
end

function test6a(numtrial::Int, expectedk::String, expectedd::String)
    NSTRINGPAIR = 50000
    m1 = SortedDict((Dict{Compat.ASCIIString, Compat.ASCIIString}()), Lt(isless))
    strlist = Compat.ASCIIString[]
    open(seekfile("wordsScram.txt"), "r") do inio
        for j = 1 : NSTRINGPAIR * 2
            push!(strlist, chomp(readline(inio)))
        end
    end
    for trial = 1 : numtrial
        empty!(m1)
        count = 1
        for j = 1 : NSTRINGPAIR
            k = strlist[count]
            d = strlist[count + 1]
            m1[k] = d
            count += 2
        end
        delct = div(NSTRINGPAIR, 6)
        for j = 1 : delct
            l = find(m1, strlist[j * 2 + 1])
            delete!((m1,l))
        end
        spec = div(NSTRINGPAIR * 3,4)
        l = startof(m1)
        for j = 1 : spec
            l = advance((m1,l))
        end
        k,d = deref((m1,l))
        ekn = expectedk
        edn = expectedd
        @test k == ekn && d == edn
    end
end

function SDConstruct(a::Associative; lt::Function=isless, by::Function=identity)
    if by == identity
        return SortedDict(a, Lt(lt))
    elseif lt == isless
        return SortedDict(a, By(by))
    else
        throw(ErrorException("having both by and lt not both implemented"))
    end
end

function test6b(numtrial::Int, expectedk::String, expectedd::String)
    NSTRINGPAIR = 50000
    m1 = SDConstruct((Dict{Compat.ASCIIString,Compat.ASCIIString}()), lt=isless)
    strlist = Compat.ASCIIString[]
    open(seekfile("wordsScram.txt"), "r") do inio
        for j = 1 : NSTRINGPAIR * 2
            push!(strlist, chomp(readline(inio)))
        end
    end
    for trial = 1 : numtrial
        empty!(m1)
        count = 1
        for j = 1 : NSTRINGPAIR
            k = strlist[count]
            d = strlist[count + 1]
            m1[k] = d
            count += 2
        end
        delct = div(NSTRINGPAIR, 6)
        for j = 1 : delct
            l = find(m1, strlist[j * 2 + 1])
            delete!((m1,l))
        end
        spec = div(NSTRINGPAIR * 3,4)
        l = startof(m1)
        for j = 1 : spec
            l = advance((m1,l))
        end
        k,d = deref((m1,l))
        ekn = expectedk
        edn = expectedd
        @test k == ekn && d == edn
    end
end

function test7()
    # Test all methods of SortedMultiDict except loops
    factors = SortedMultiDict(Int[], Int[])
    N = 1000
    checkcorrectness(factors.bt, true)
    len = 0
    for factor = 1 : N
        for multiple = factor : factor : N
            insert!(factors, multiple, factor)
            len += 1
        end
    end
    @test length(factors) == len
    @test Pair(70,2) in factors
    @test Pair(70,14) in factors
    @test !(Pair(70,15) in factors)
    @test !(Pair(N+1,15) in factors)
    @test eltype(factors) == Pair{Int,Int}
    @test eltype(typeof(factors)) == Pair{Int,Int}

    @test keytype(factors) == Int
    @test keytype(typeof(factors)) == Int
    @test valtype(factors) == Int
    @test valtype(typeof(factors)) == Int
    @test ordtype(factors) == ForwardOrdering
    @test ordtype(typeof(factors)) == ForwardOrdering

    push_test!(factors, 70,3)
    @test length(factors) == len+1
    @test Pair(70,3) in factors
    i = searchsortedfirst(factors, 70)
    dcount = 0
    for (s,k,v) in semitokens(inclusive(factors, i, endof(factors)))
        @test k == 70
        if v == 3
            delete!((factors,s))
            dcount += 1
            break
        end
    end
    @test dcount == 1

    @test orderobject(factors) == Forward
    @test haskey(factors, 60)
    @test !haskey(factors, -1)
    @test 60 in keys(factors)
    @test !(-1 in keys(factors))
    checkcorrectness(factors.bt, true)
    i = startof(factors)
    i = advance((factors,i))
    @test deref((factors,i)) == Pair(2,1)
    @test deref_key((factors,i)) == 2
    @test deref_value((factors,i)) == 1
    @test factors[i] == 1
    factors[i] = 7
    @test deref((factors,i)) == Pair(2,7)
    factors[i] = 1
    i = regress((factors,i))
    i = regress((factors,i))
    @test i == beforestartsemitoken(factors)
    pr = first(factors)
    @test pr == Pair(1,1)
    pr2 = last(factors)
    @test pr2 == Pair(N,N)
    i = searchsortedfirst(factors,77)
    @test deref((factors,i)) == Pair(77,1)
    i = searchsortedlast(factors,77)
    @test deref((factors,i)) == Pair(77,77)
    i = searchsortedafter(factors,77)
    @test deref((factors,i)) == Pair(78,1)
    expected = [1,2,4,5,8,10,16,20,40,80]
    i1,i2 = searchequalrange(factors, 80)
    i = i1
    for e in expected
        @test deref_value((factors,i)) == e
        i = advance((factors,i))
    end
    @test compare(factors,i,i2) != 0
    @test compare(factors,regress((factors,i)),i2) == 0
    @test compare(factors,i,i1) != 0
    insert!(factors, 80, 6)
    @test length(factors) == len + 1
    checkcorrectness(factors.bt, true)
    expected1 = deepcopy(expected)
    push!(expected1, 6)
    i1,i2 = searchequalrange(factors, 80)
    i = i1
    for e in expected1
        @test deref_value((factors,i)) == e
        i = advance((factors,i))
    end
    @test compare(factors,i2,regress((factors,i))) == 0
    @test compare(factors,i,i1) != 0
    delete!((factors,i2))
    @test length(factors) == len
    checkcorrectness(factors.bt, true)
    i1,i2 = searchequalrange(factors, 80)
    i = i1
    for e in expected
        @test deref_value((factors,i)) == e
        i = advance((factors,i))
    end
    @test compare(factors,regress((factors,i)),i2) == 0
    @test !isempty(factors)
    empty!(factors)
    checkcorrectness(factors.bt, true)
    @test length(factors) == 0
    @test isempty(factors)
    i = startof(factors)
    @test i == pastendsemitoken(factors)
    i = endof(factors)
    @test i == beforestartsemitoken(factors)
    i1,i2 = searchequalrange(factors, N + 2)
    @test i1 == pastendsemitoken(factors)
    @test i2 == beforestartsemitoken(factors)
    m1 = SortedMultiDict(["apples", "apples", "bananas"], [2.0, 1.0,1.5])
    checkcorrectness(m1.bt, true)
    m2 = SortedMultiDict(["bananas","apples", "apples"], [1.5, 2.0, 1.0])
    checkcorrectness(m2.bt, true)
    m3 = SortedMultiDict(["apples", "apples", "bananas"], [1.0, 2.0, 1.5])
    checkcorrectness(m3.bt, true)
    @test isequal(m1,m2)
    @test !isequal(m1,m3)
    @test !isequal(m1, SortedMultiDict(["apples"], [2.0]))
    stok = insert!(m2, "cherries", 6.1)
    checkcorrectness(m2.bt, true)
    @test !isequal(m1,m2)
    delete!((m2,stok))
    checkcorrectness(m2.bt, true)
    @test isequal(m1,m2)
    m4 = deepcopy(m2)
    checkcorrectness(m4.bt, true)
    @test isequal(m1,m4)
    m5 = packcopy(m2)
    checkcorrectness(m5.bt, true)
    @test isequal(m1,m5)
    m6 = packdeepcopy(m2)
    checkcorrectness(m6.bt, true)
    @test isequal(m1,m6)

    m1 = SortedMultiDict(["bananas", "apples", "cherries", "cherries", "oranges"],
                         [1.0, 2.0, 3.0, 4.0, 5.0])
    m2 = SortedMultiDict(["apples", "cherries", "cherries", "bananas", "plums"],
                         [6.0, 7.0, 8.0, 9.0, 10.0])
    m3 = SortedMultiDict(["apples", "apples", "bananas", "bananas",
                          "cherries", "cherries", "cherries", "cherries",
                          "oranges", "plums"],
                         [2.0, 6.0, 1.0, 9.0, 3.0, 4.0, 7.0, 8.0, 5.0, 10.0])
    m3empty = similar(m3)
    @test (eltype(m3empty) == Pair{Compat.ASCIIString, Float64}) &&
        orderobject(m3empty) == Forward &&
        length(m3empty) == 0
    m4 = merge(m1, m2)
    @test isequal(m3, m4)
    m5 = merge(m2, m1)
    @test !isequal(m3, m5)
    merge!(m1, m2)
    @test isequal(m1, m3)
    m7 = SortedMultiDict(Int[], Int[])
    n1 = 10000
    for k = 1 : n1
        insert!(m7, k, k+1)
    end
    for k = 1 : n1
        insert!(m7, k, k+2)
    end
    for k = 1 : n1
        i1, i2 = searchequalrange(m7, k)
        count = 0
        for (key,v) in inclusive(m7, i1, i2)
            count += 1
            @test key == k
            @test v == k + count
        end
        @test count == 2
        count = 0
        for (key,v) in inclusive(m7, searchequalrange(m7,k))
            count += 1
            @test key == k
            @test v == k + count
        end
        @test count == 2
    end
    nothing
end

function test8()
    # Test SortedSet
    N = 1000
    sm = 0.0
    m = SortedSet(Float64[])
    smallest = 10.0
    largest = -10.0
    for j = 1 : N
        u = j * exp(1)
        ui = u - floor(u)
        push!(m, ui)
        sm += ui
        smallest = min(smallest,ui)
        largest = max(largest,ui)
    end
    isnew,st = insert!(m, 72.5)
    @test isnew
    @test deref((m,st)) == 72.5
    delete!((m,st))
    isnew,st = insert!(m, 73.5)
    @test isnew
    @test deref((m,st)) == 73.5
    delete!(m, 73.5)
    checkcorrectness(m.bt, false)
    count = 0
    sm2 = 0.0
    prev = -1.0
    for k in m
        sm2 += k
        count += 1
        @test k >= prev
    end
    @test abs(sm2 - sm) <= 1e-10
    @test count == N
    @test length(m) == N
    ii2 = searchsortedfirst(m, 0.5)
    i3 = startof(m)
    v = first(m)
    @test v == smallest
    @test deref((m,i3)) == v
    i4 = endof(m)
    w = last(m)
    @test w == largest
    @test deref((m,i4)) == w
    i5 = beforestartsemitoken(m)
    @test advance((m,i5)) == i3
    i6 = pastendsemitoken(m)
    @test regress((m,i6)) == i4
    @test advance((m,i5)) != i4
    @test regress((m,i6)) != i3
    j1 = searchsortedfirst(m,0.5)
    j2 = searchsortedlast(m,0.5)
    j3 = searchsortedafter(m,0.5)
    @test deref((m,j1)) > 0.5
    @test deref((m,j2)) < 0.5
    @test advance((m,j2)) == j1
    @test j1 == j3
    k1 = searchsortedfirst(m,smallest)
    k2 = searchsortedlast(m,smallest)
    k3 = searchsortedafter(m,smallest)
    @test deref((m,k1)) == smallest
    @test deref((m,k2)) == smallest
    @test deref((m,regress((m,k3)))) == smallest
    secondsmallest = deref((m,k3))
    sk = searchsortedfirst(m,0.4)
    ek = searchsortedfirst(m,0.6)
    dcount = 0
    for (st,k) in semitokens(exclusive(m,sk,ek))
        delete!((m,st))
        dcount += 1
        if dcount % 20 == 0
            checkcorrectness(m.bt, false)
        end
    end
    newcount = 0
    for k in m
        newcount += 1
        @test k < 0.4 || k > 0.6
    end
    @test newcount == N - dcount
    @test smallest in m
    @test haskey(m,smallest)
    @test !(0.5 in m)
    @test !haskey(m,0.5)
    @test eltype(m) == Float64
    @test eltype(typeof(m)) == Float64
    @test keytype(m) == Float64
    @test keytype(typeof(m)) == Float64
    @test orderobject(m) == Forward
    @test ordtype(m) == ForwardOrdering
    @test ordtype(typeof(m)) == ForwardOrdering
    pop!(m, smallest)
    checkcorrectness(m.bt, false)
    @test length(m) == N - dcount - 1
    key1 = pop!(m)
    @test key1 == secondsmallest
    @test length(m) == N - dcount - 2
    checkcorrectness(m.bt, false)
    @test !isempty(m)
    empty!(m)
    @test isempty(m)
    m1 = SortedSet(["blue", "orange", "red"])
    m2 = SortedSet(["orange", "blue", "red"])
    m3 = SortedSet(["orange", "yellow", "red"])
    m3empty = similar(m3)
    @test eltype(m3empty) == Compat.ASCIIString &&
       length(m3empty) == 0
    @test isequal(m1,m2)
    @test !isequal(m1,m3)
    @test !isequal(m1, SortedSet(["blue"]))
    m4 = packcopy(m3)
    @test isequal(m3,m4)
    m5 = packdeepcopy(m4)
    @test isequal(m3,m4)
    m6 = deepcopy(m5)
    @test isequal(m3,m5)
    checkcorrectness(m1.bt, false)
    checkcorrectness(m2.bt, false)
    checkcorrectness(m3.bt, false)
    checkcorrectness(m4.bt, false)
    checkcorrectness(m5.bt, false)
    checkcorrectness(m5.bt, false)
    m7 = union(m1, ["yellow"])
    m8 = union(m3, SortedSet(["blue"]))
    @test isequal(m7,m8)
    @test !isequal(m1,m8)
    union!(m1, ["yellow"])
    @test isequal(m1,m8)
    m8a = intersect(m8)
    @test isequal(m8a,m8)
    m9 = intersect(m8, SortedSet(["yellow", "red", "white"]))
    @test isequal(m9, SortedSet(["red", "yellow"]))
    m9a = intersect(m8, SortedSet(["yellow", "red", "white"]), m8)
    @test isequal(m9a, SortedSet(["red", "yellow"]))
    m10 = symdiff(m8,  SortedSet(["yellow", "red", "white"]))
    @test isequal(m10, SortedSet(["white", "blue", "orange"]))
    m11 = symdiff(m8, SortedSet(["yellow", "red", "blue", "orange",
                                 "zinc"]))
    @test isequal(m11, SortedSet(["zinc"]))
    m12 = symdiff(SortedSet(["yellow", "red", "blue", "orange",
                                 "zinc"]), m8)
    @test isequal(m12, SortedSet(["zinc"]))
    m13 = setdiff(m8, SortedSet(["yellow", "red", "white"]))
    @test isequal(m13, SortedSet(["blue", "orange"]))
    m14 = setdiff(m8, SortedSet(["blue"]))
    @test isequal(m14, SortedSet(["orange", "yellow", "red"]))
    @test issubset(["yellow", "blue"], m8)
    @test !issubset(["blue", "green"], m8)
    setdiff!(m8, SortedSet(["yellow", "red", "white"]))
    @test isequal(m8, SortedSet(["blue", "orange"]))
    nothing
end

# test the constructors of SortedDict and SortedMultiDict

function test9()

    sd1 = SortedDict("w" => 64, "p" => 12)
    @test length(sd1) == 2 && first(sd1) == ("p"=>12) &&
        last(sd1) == ("w"=>64)
    sd2 = SortedDict(Reverse, "w" => 64, "p" => 12)
    @test length(sd2) == 2 && last(sd2) == ("p"=>12) &&
        first(sd2) == ("w"=>64)
    sd3 = SortedDict(("w"=>64, "p"=>12))
    @test length(sd3) == 2 && first(sd3) == ("p"=>12) &&
        last(sd3) == ("w"=>64)
    sd4 = SortedDict(("w"=>64, "p"=>12), Reverse)
    @test length(sd4) == 2 && last(sd4) == ("p"=>12) &&
        first(sd4) == ("w"=>64)
    sm1 = SortedMultiDict("w" => 64, "p" => 12, "p" => 9)
    @test length(sm1) == 3 && first(sm1) == ("p"=>12) &&
        last(sm1) == ("w"=>64)
    sm2 = SortedMultiDict(Reverse, "w" => 64, "p" => 12, "p" => 9)
    @test length(sm2) == 3 && last(sm2) == ("p"=>9) &&
        first(sm2) == ("w"=>64)
    sm3 = SortedMultiDict(("w"=>64, "p"=>12, "p"=> 9))
    @test length(sm3) == 3 && first(sm3) == ("p"=>12) &&
        last(sm3) == ("w"=>64)
    sm4 = SortedMultiDict(("w"=> 64, "p"=>12, "p"=>9), Reverse)
    @test length(sm4) == 3 && last(sm4) == ("p"=>9) &&
        first(sm4) == ("w"=>64)
    nothing
end

test1()
test2()
test3(0x00000000)
test4()
test5()
#test6(2, "soothingly", "compere")
test7()
test8()
test9()


text:
"""
	PenaltyWrapper

It is an inprocessing algorithm that wraps a Probabilistic classifier.
Optimal thresholds for each group in protected attribute is found to minimise accuracy+ alpha*fairness_measure^2.
Gradients are used to find optimal threshold values. alpha controls the fairness accuracy tradeoff.
"""
struct PenaltyWrapper{M<:MLJBase.Model} <: DeterministicComposite
	grp::Symbol
	classifier::M
	measure::MLJBase.Measure
	alpha::Real
	n_iters::Int
	lr::Float64
end

"""
	PenaltyWrapper(classifier=nothing, grp=:class, measure, alpha=1, n_iters=1000, lr=0.01)

Instantiates PenaltyWrapper which wraps the classifier and contains measure to be optimised.
It also provides various hyperparameters to the wrapper.
"""
function PenaltyWrapper(;
	grp::Symbol = :class,
    classifier::MLJBase.Model = nothing,
    measure::MLJBase.Measure,
    alpha::Real = 1, # Fairess-Accuracy tradeoff parameter in loss function
    n_iters = 10^3, # Number of gradient steps to perform per group
	lr::Float64 = 10^-2 # Learning Rate
)
    model = PenaltyWrapper(grp, classifier, measure, alpha, n_iters, lr)
    message = MLJBase.clean!(model)
    isempty(message) || @warn message
    return model
end

function MLJBase.clean!(model::PenaltyWrapper)
	warning = ""
    model.classifier!=nothing || (warning *= "No classifier specified in model\n")
    target_scitype(model.classifier) <: AbstractVector{<:Finite} || (warning *= "Only Binary Classifiers are supported\n")
    return warning
end

function MMI.fit(model::PenaltyWrapper, verbosity::Int,
	X, y)
	grps = X[:, model.grp]
	labels = levels(y)
	favLabel = labels[2]
	unfavLabel = labels[1]
	n_grps = length(levels(grps)) # Number of different values for sensitive attribute
	n = length(y)

	y = categorical(y.==favLabel)

	grp_idx = Dict()
	levels_ = levels(grps)
	for i in 1:n_grps
	    grp_idx[levels_[i]] = i
	end

	λ = fill(0.5, n_grps) #Thresholds for various groups
	mat = Array{Any, 3}(undef, n_grps, 2, 2)
	mat .= 0.25
	ft = FairTensor(mat, levels(grps))

	score(x) = x.prob_given_ref[2]
	all_indices = convert(Array, 1:n)
	for itr in 1:model.n_iters
		train_indices = shuffle!(all_indices)[1:Int(round(0.7n))]
		tune_indices = all_indices[Int(round(0.7n))+1:n]
		mch = machine(model.classifier, X[train_indices, :], y[train_indices])
		fit!(mch, verbosity=0)
		ŷ = MMI.predict(mch, X[tune_indices, :])
		for i in 1:n_grps
			indices = grps[tune_indices].==levels_[i]
			y_tune_indices = convert(Array, y[tune_indices][indices])
			if all(y_tune_indices.==unfavLabel) || all(y_tune_indices.==favLabel) continue end
			scores = score.(ŷ[indices])
			dis0 = Distributions.fit(Distributions.Normal, scores[y_tune_indices .== unfavLabel])
			dis1 = Distributions.fit(Distributions.Normal, scores[y_tune_indices .== favLabel])
			function lambdafair(λ)
				pr = Array{Any, 2}(undef, 2, 2)
				pr[1, 1] = 1-Distributions.cdf(dis1, λ) #TruePositive
				pr[1, 2] = 1-Distributions.cdf(dis0, λ) #FalsePositive
				pr[2, 1] = Distributions.cdf(dis1, λ) #FalseNegative
				pr[2, 2] = Distributions.cdf(dis0, λ) #TrueNegative
				pr /= sum(pr)
				ft[i, :, :] = pr
				ft
			end
			compositeloss(ft::FairTensor) =  accuracy(ft) +
				model.alpha*(model.measure(ft, grp=levels_[i])-model.measure(ft, grp=levels_[1]))^2
			loss(λ) = compositeloss(lambdafair(λ[1]))

			λᵢ_grad = ForwardDiff.gradient(loss, [λ[i]])
			λ[i] -= λᵢ_grad[1]*model.lr
		end
	end

	# Make y and ŷ as an attribute. And then predict whether that is true or false
	mach = machine(model.classifier, X, y)
	fit!(mach, verbosity=verbosity)
	fitresult = [mach.fitresult, λ]
	return fitresult, nothing, nothing
end

function MMI.predict(model::PenaltyWrapper, fitresult, Xnew)
	classifier_fitresult, λ = fitresult
	n = size(Xnew)[1]
	grps = Xnew[:, model.grp]
	levels_ = levels(grps)
	n_grps = length(levels(grps))
	grp_idx = Dict() # Maps group_name => index of group in levels_
	for i in 1:n_grps
		grp_idx[levels_[i]] = i
	end

	score(x) = x.prob_given_ref[2]
	preds = MMI.predict(model.classifier, classifier_fitresult, Xnew) # preds are UnivariateFinite
	scores = score.(preds)
	ŷ = zeros(Int, n)
	for i in 1:n
		ŷ[i] = scores[i]>λ[grp_idx[Xnew[i, model.grp]]]
	end
	return ŷ
end

MMI.input_scitype(::Type{<:PenaltyWrapper{M}}) where M = input_scitype(M)
MMI.target_scitype(::Type{<:PenaltyWrapper{M}}) where M = AbstractVector{<:Finite{2}}


text:
using MLDatasets: CIFAR10

struct CIFAR10Dataset{T, D} <: ImageDataset{T, D}
    X
    y
    Y
    Xt
    yt
    Yt
end

n_display(::CIFAR10Dataset) = 64

function CIFAR10Dataset(
    n_data::Int=50_000;
    seed::Int=1,
    test_ratio=1/5,
    n_test::Int=ratio2num(n_data, test_ratio),
    is_flatten::Bool=false,
    alpha::T=5f-2,
    is_link::Bool=true,
    rng=MersenneTwister(seed),
) where {T<:AbstractFloat}
    if n_data > 50_000
        @warn "Attempts to access more than 50,000 train points in CIFAR10; clamped to 50,000."
        n_data = 50_000
    end
    if n_test > 10_000
        @warn "Attempts to access more than 10,000 test points in CIFAR10; clamped to 10,000."
        n_test = 10_000
    end
    X, y, Y, Xt, yt, Yt = get_image_data(
        CIFAR10, n_data, n_test, is_flatten, alpha, is_link, 10, (2, 1, 3, 4); rng=rng
    )
    @info "Oh you just get the CIFAR10 dataset" n_data=n_data n_test=n_test is_flatten=is_flatten alpha=alpha is_link=is_link
    return CIFAR10Dataset{Val{is_link}, 3072}(X, y, Y, Xt, yt, Yt)
end