text:
#
#
# Encapulates arbitrary solvers by a debug solver, which is induced by Decoration Options
#
#
#
#
# A general framework for solvers of problems on Manifolds
#
# This file introduces fallbacks for not yet implemented parts and the general
# function to run the solver
export initializeSolver!, doSolverStep!, getSolverResult, stopSolver!
"""
    initializeSolver!(p,o)
Initialize the solver to the optimization [`Problem`](@ref) by initializing all
values in the [`DebugOptions`](@ref)` o`.
"""
function initializeSolver!(p::P,o::O) where {P <: Problem, O <: DebugOptions}
    initializeSolver!(p,o.options)
    get(o.debugDictionary,:Start,DebugDivider(""))(p,getOptions(o),0)
    get(o.debugDictionary,:All,DebugDivider(""))(p,getOptions(o),0)
end
"""
    doSolverStep!(p,o,iter)
Do one iteration step (the `iter`th) for [`Problem`](@ref)` p` by modifying
the values in the [`Options`](@ref)` o.options` and print `Debug`.
"""
function doSolverStep!(p::P,o::O, i) where {P <: Problem, O <: DebugOptions}
    doSolverStep!(p,o.options,i)
    get(o.debugDictionary,:Step,DebugDivider(""))(p,getOptions(o),i)
    get(o.debugDictionary,:All,DebugDivider(""))(p,getOptions(o),i)
end
"""
    getSolverResult(p,o)
Return the final result after all iterations that is stored within the
(modified during the iterations) [`Options`](@ref) `o`.
"""
function getSolverResult(p::P,o::O) where {P <: Problem, O <: DebugOptions}
    return getSolverResult(p, o.options)
end
"""
    stopSolver!(p,o,i)

determine whether the solver for [`Problem`](@ref) `p` and the [`DebugOptions`](@ref) `o`
should stop at iteration `i`. If so, print all debug from `:All` and `:Final`.
"""
function stopSolver!(p::P,o::O,i::Int) where {P <: Problem, O <: DebugOptions}
    s = stopSolver!(p,o.options,i)
    if s
        get(o.debugDictionary,:Stop,DebugDivider(""))(p,getOptions(o),typemin(Int))
        get(o.debugDictionary,:All,DebugDivider(""))(p,getOptions(o),typemin(Int))
    end
    return s
end

ext: jl


text:
using NCDatasets

"""
Run a coarse thermal bubble simulation and save the output to NetCDF at the
10th time step. Then read back the output and test that it matches the model's
state.
"""
function run_thermal_bubble_netcdf_tests(arch)
    Nx, Ny, Nz = 16, 16, 16
    Lx, Ly, Lz = 100, 100, 100

    grid = RegularCartesianGrid(size=(Nx, Ny, Nz), length=(Lx, Ly, Lz))
    closure = ConstantIsotropicDiffusivity(ν=4e-2, κ=4e-2)
    model = IncompressibleModel(architecture=arch, grid=grid, closure=closure)
    simulation = Simulation(model, Δt=6, stop_iteration=10)

    # Add a cube-shaped warm temperature anomaly that takes up the middle 50%
    # of the domain volume.
    i1, i2 = round(Int, Nx/4), round(Int, 3Nx/4)
    j1, j2 = round(Int, Ny/4), round(Int, 3Ny/4)
    k1, k2 = round(Int, Nz/4), round(Int, 3Nz/4)
    model.tracers.T.data[i1:i2, j1:j2, k1:k2] .+= 0.01

    outputs = Dict(
        "v" => model.velocities.v,
        "u" => model.velocities.u,
        "w" => model.velocities.w,
        "T" => model.tracers.T,
        "S" => model.tracers.S
    )

    nc_writer = NetCDFOutputWriter(model, outputs, filename="dump_test.nc", frequency=10)
    push!(simulation.output_writers, nc_writer)

    xC_slice = 1:10
    xF_slice = 2:11
    yC_slice = 10:15
    yF_slice = 1
    zC_slice = 10
    zF_slice = 9:11

    nc_sliced_writer =
        NetCDFOutputWriter(model, outputs, filename="dump_test_sliced.nc", frequency=10,
                           xC=xC_slice, xF=xF_slice, yC=yC_slice,
                           yF=yF_slice, zC=zC_slice, zF=zF_slice)

    push!(simulation.output_writers, nc_sliced_writer)

    run!(simulation)
    close(nc_writer)
    close(nc_sliced_writer)

    ds3 = Dataset("dump_test.nc")
    u = ds3["u"][:, :, :, end]
    v = ds3["v"][:, :, :, end]
    w = ds3["w"][:, :, :, end]
    T = ds3["T"][:, :, :, end]
    S = ds3["S"][:, :, :, end]
    close(ds3)

    @test all(u .≈ Array(interiorparent(model.velocities.u)))
    @test all(v .≈ Array(interiorparent(model.velocities.v)))
    @test all(w .≈ Array(interiorparent(model.velocities.w)))
    @test all(T .≈ Array(interiorparent(model.tracers.T)))
    @test all(S .≈ Array(interiorparent(model.tracers.S)))

    ds2 = Dataset("dump_test_sliced.nc")
    u_sliced = ds2["u"][:, :, :, end]
    v_sliced = ds2["v"][:, :, :, end]
    w_sliced = ds2["w"][:, :, :, end]
    T_sliced = ds2["T"][:, :, :, end]
    S_sliced = ds2["S"][:, :, :, end]
    close(ds2)

    @test all(u_sliced .≈ Array(interiorparent(model.velocities.u))[xF_slice, yC_slice, zC_slice])
    @test all(v_sliced .≈ Array(interiorparent(model.velocities.v))[xC_slice, yF_slice, zC_slice])
    @test all(w_sliced .≈ Array(interiorparent(model.velocities.w))[xC_slice, yC_slice, zF_slice])
    @test all(T_sliced .≈ Array(interiorparent(model.tracers.T))[xC_slice, yC_slice, zC_slice])
    @test all(S_sliced .≈ Array(interiorparent(model.tracers.S))[xC_slice, yC_slice, zC_slice])
end

function run_netcdf_function_output_tests(arch)
    N = 16
    L = 1
    model = IncompressibleModel(grid=RegularCartesianGrid(size=(N, N, N), length=(L, 2L, 3L)))
    simulation = Simulation(model, Δt=1.25, stop_iteration=3)

    # Define scalar, vector, 2D slice, and 3D field outputs
    f(model) = model.clock.time^2
    g(model) = @. model.clock.time * exp(model.grid.zC)
    h(model) = @. model.clock.time * sin(model.grid.xC) * cos(model.grid.yC')

    outputs = Dict("scalar" => f, "profile" => g, "slice" => h)
    dims = Dict("scalar" => (), "profile" => ("zC",), "slice" => ("xC", "yC"))

    output_attributes = Dict(
        "scalar"  => Dict("longname" => "Some scalar", "units" => "bananas"),
        "profile" => Dict("longname" => "Some vertical profile", "units" => "watermelons"),
        "slice"   => Dict("longname" => "Some slice", "units" => "mushrooms")
    )

    global_attributes = Dict("location" => "Bay of Fundy", "onions" => 7)

    simulation.output_writers[:fruits] =
        NetCDFOutputWriter(
            model, outputs; frequency=1, filename="test_function_outputs.nc", dimensions=dims,
            global_attributes=global_attributes, output_attributes=output_attributes)

    run!(simulation)
    close(simulation.output_writers[:fruits])

    ds = Dataset("test_function_outputs.nc", "r")

    @test ds.attrib["location"] == "Bay of Fundy"
    @test ds.attrib["onions"] == 7

    @test length(ds["time"]) == 4
    @test ds["time"][:] == [1.25i for i in 0:3]

    @test ds["scalar"].attrib["longname"] == "Some scalar"
    @test ds["scalar"].attrib["units"] == "bananas"
    @test ds["scalar"][:] == [(1.25i)^2 for i in 0:3]
    @test dimnames(ds["scalar"]) == ("time",)

    @test ds["profile"].attrib["longname"] == "Some vertical profile"
    @test ds["profile"].attrib["units"] == "watermelons"
    @test ds["profile"][:, end] == @. 3.75 * exp(model.grid.zC)
    @test size(ds["profile"]) == (N, 4)
    @test dimnames(ds["profile"]) == ("zC", "time")

    @test ds["slice"].attrib["longname"] == "Some slice"
    @test ds["slice"].attrib["units"] == "mushrooms"
    @test ds["slice"][:, :, end] == @. 3.75 * sin(model.grid.xC) * cos(model.grid.yC')
    @test size(ds["slice"]) == (N, N, 4)
    @test dimnames(ds["slice"]) == ("xC", "yC", "time")

    close(ds)
    return nothing
end

function run_jld2_file_splitting_tests(arch)
    model = IncompressibleModel(grid=RegularCartesianGrid(size=(16, 16, 16), length=(1, 1, 1)))
    simulation = Simulation(model, Δt=1, stop_iteration=10)

    u(model) = Array(model.velocities.u.data.parent)
    fields = Dict(:u => u)

    function fake_bc_init(file, model)
        file["boundary_conditions/fake"] = π
    end

    ow = JLD2OutputWriter(model, fields; dir=".", prefix="test", frequency=1,
                          init=fake_bc_init, including=[:grid],
                          max_filesize=200KiB, force=true)

    push!(simulation.output_writers, ow)

    # 531 KiB of output will be written which should get split into 3 files.
    run!(simulation)

    # Test that files has been split according to size as expected.
    @test filesize("test_part1.jld2") > 200KiB
    @test filesize("test_part2.jld2") > 200KiB
    @test filesize("test_part3.jld2") < 200KiB

    for n in string.(1:3)
        filename = "test_part" * n * ".jld2"
        jldopen(filename, "r") do file
            # Test to make sure all files contain structs from `including`.
            @test file["grid/Nx"] == 16

            # Test to make sure all files contain info from `init` function.
            @test file["boundary_conditions/fake"] == π
        end

        # Leave test directory clean.
        rm(filename)
    end
end

"""
Run two coarse rising thermal bubble simulations and make sure that when
restarting from a checkpoint, the restarted simulation matches the non-restarted
simulation numerically.
"""
function run_thermal_bubble_checkpointer_tests(arch)
    Nx, Ny, Nz = 16, 16, 16
    Lx, Ly, Lz = 100, 100, 100
    Δt = 6

    grid = RegularCartesianGrid(size=(Nx, Ny, Nz), length=(Lx, Ly, Lz))
    closure = ConstantIsotropicDiffusivity(ν=4e-2, κ=4e-2)
    true_model = IncompressibleModel(architecture=arch, grid=grid, closure=closure)

    # Add a cube-shaped warm temperature anomaly that takes up the middle 50%
    # of the domain volume.
    i1, i2 = round(Int, Nx/4), round(Int, 3Nx/4)
    j1, j2 = round(Int, Ny/4), round(Int, 3Ny/4)
    k1, k2 = round(Int, Nz/4), round(Int, 3Nz/4)
    true_model.tracers.T.data[i1:i2, j1:j2, k1:k2] .+= 0.01

    checkpointed_model = deepcopy(true_model)

    true_simulation = Simulation(true_model, Δt=Δt, stop_iteration=9)
    run!(true_simulation)

    checkpointed_simulation = Simulation(checkpointed_model, Δt=Δt, stop_iteration=5)
    checkpointer = Checkpointer(checkpointed_model, frequency=5, force=true)
    push!(checkpointed_simulation.output_writers, checkpointer)

    # Checkpoint should be saved as "checkpoint5.jld" after the 5th iteration.
    run!(checkpointed_simulation)

    # Remove all knowledge of the checkpointed model.
    checkpointed_model = nothing

    # model_kwargs = Dict{Symbol, Any}(:boundary_conditions => SolutionBoundaryConditions(grid))
    restored_model = restore_from_checkpoint("checkpoint_iteration5.jld2")

    for n in 1:4
        time_step!(restored_model, Δt, euler=false)
    end

    rm("checkpoint_iteration0.jld2", force=true)
    rm("checkpoint_iteration5.jld2", force=true)

    # Now the true_model and restored_model should be identical.
    @test all(restored_model.velocities.u.data     .≈ true_model.velocities.u.data)
    @test all(restored_model.velocities.v.data     .≈ true_model.velocities.v.data)
    @test all(restored_model.velocities.w.data     .≈ true_model.velocities.w.data)
    @test all(restored_model.tracers.T.data        .≈ true_model.tracers.T.data)
    @test all(restored_model.tracers.S.data        .≈ true_model.tracers.S.data)
    @test all(restored_model.timestepper.Gⁿ.u.data .≈ true_model.timestepper.Gⁿ.u.data)
    @test all(restored_model.timestepper.Gⁿ.v.data .≈ true_model.timestepper.Gⁿ.v.data)
    @test all(restored_model.timestepper.Gⁿ.w.data .≈ true_model.timestepper.Gⁿ.w.data)
    @test all(restored_model.timestepper.Gⁿ.T.data .≈ true_model.timestepper.Gⁿ.T.data)
    @test all(restored_model.timestepper.Gⁿ.S.data .≈ true_model.timestepper.Gⁿ.S.data)
end

@testset "Output writers" begin
    @info "Testing output writers..."

    for arch in archs
         @testset "NetCDF [$(typeof(arch))]" begin
             @info "  Testing NetCDF output writer [$(typeof(arch))]..."
             run_thermal_bubble_netcdf_tests(arch)
             run_netcdf_function_output_tests(arch)
         end

        @testset "JLD2 [$(typeof(arch))]" begin
            @info "  Testing JLD2 output writer [$(typeof(arch))]..."
            run_jld2_file_splitting_tests(arch)
        end

        @testset "Checkpointer [$(typeof(arch))]" begin
            @info "  Testing Checkpointer [$(typeof(arch))]..."
            run_thermal_bubble_checkpointer_tests(arch)
        end
    end
end


text:
##########################################################################
# The RandomLargeScaleNetwork is a network whose geography changes in
# each user drop. Both BSs and MSs are uniformly distributed over
# a rectangle.

##########################################################################
# Network definition
type RandomLargeScaleNetwork{MS_t <: PhysicalMS, BS_t <: PhysicalBS, System_t <: System, PropagationEnvironment_t <: PropagationEnvironment} <: PhysicalNetwork
    MSs::Vector{MS_t}
    BSs::Vector{BS_t}

    system::System_t
    num_MSs_per_cell::Int
    propagation_environment::PropagationEnvironment_t
    geography_size::@compat Tuple{Float64, Float64}
    MS_serving_BS_distance::Nullable{Float64} # null signifies random placement of MSs. Otherwise they are placed on a circle from the BS.
    aux_network_params::AuxNetworkParams

    assignment::Assignment
end

# Convenience constructor without network params and assignments
RandomLargeScaleNetwork(MSs, BSs, system, num_MSs_per_cell, propagation_environment, geography_size, MS_serving_BS_distance) =
    RandomLargeScaleNetwork(MSs, BSs, system, num_MSs_per_cell, propagation_environment, geography_size, MS_serving_BS_distance, AuxNetworkParams(), Assignment())

Base.show(io::IO, x::RandomLargeScaleNetwork) =
    print(io, "RandomLargeScaleNetwork(I = $(length(x.BSs)), Kc = $(x.num_MSs_per_cell), geography_size = $(x.geography_size), MS_serving_BS_distance = $(x.MS_serving_BS_distance))")
Base.showcompact(io::IO, x::RandomLargeScaleNetwork) =
    print(io, "RandomLargeScaleNetwork($(length(x.BSs)), $(x.num_MSs_per_cell), $(x.geography_size), $(x.MS_serving_BS_distance))")

function setup_random_large_scale_network(
    num_BSs, num_MSs_per_cell, num_MS_antennas, num_BS_antennas;
    system = SinglecarrierSystem(2e9, 15e3),
    propagation_environment = SimpleLargescaleFadingEnvironment(37.6, 15.3, 10, 8),
    geography_size = (1500., 1500.),
    MS_serving_BS_distance = Nullable(150.), # Nullable{Float64}() for random MS placement
    transmit_power = 10^(18.2/10), transmit_powers = transmit_power*ones(Float64, num_BSs),
    BS_antenna_gain_params = [ OmnidirectionalAntennaParams(0) for idx = 1:num_BSs ],
    user_priority = 1., user_priorities = user_priority*ones(Float64, num_BSs*num_MSs_per_cell),
    num_streams = 1, num_streamss = num_streams*ones(Int, num_BSs*num_MSs_per_cell),
    MS_antenna_gain_dB = 0., MS_antenna_gains_dB = MS_antenna_gain_dB*ones(Float64, num_BSs*num_MSs_per_cell),
    receiver_noise_figure = 9., receiver_noise_figures = receiver_noise_figure*ones(Float64, num_BSs*num_MSs_per_cell))

    isa(num_MS_antennas, Vector) || (num_MS_antennas = num_MS_antennas*ones(Int, num_BSs*num_MSs_per_cell))
    isa(num_BS_antennas, Vector) || (num_BS_antennas = num_BS_antennas*ones(Int, num_BSs))

    BSs = [ PhysicalBS(num_BS_antennas[i], Position(0, 0), transmit_powers[i], BS_antenna_gain_params[i]) for i = 1:num_BSs ]
    MSs = [ PhysicalMS(num_MS_antennas[k], Position(0, 0), Velocity(0, 0), user_priorities[k], num_streamss[k], MS_antenna_gains_dB[k], receiver_noise_figures[k], SimpleLargescaleFadingEnvironmentState(zeros(Float64, num_BSs), falses(num_BSs))) for k = 1:num_BSs*num_MSs_per_cell ]

    RandomLargeScaleNetwork(MSs, BSs, system, num_MSs_per_cell, propagation_environment, geography_size, MS_serving_BS_distance)
end

##########################################################################
# SNR functions (only work if MS_serving_BS_distance is not null)

function get_average_SNRs_dB{MS_t <: PhysicalMS, BS_t <: PhysicalBS, System_t <: System}(network::RandomLargeScaleNetwork{MS_t, BS_t, System_t, SimpleLargescaleFadingEnvironment})
    isnull(network.MS_serving_BS_distance) && error("Can only call get_average_SNRs_dB if BS-MS distance for served MSs is specified.")

    I = get_num_BSs(network); Kc = network.num_MSs_per_cell
    pathloss_alpha = network.propagation_environment.pathloss_alpha
    pathloss_beta = network.propagation_environment.pathloss_beta

    SNRs_dB = zeros(Float64, I*Kc)
    for i = 1:I
        Pt_dBm = get_transmit_power_dBm(network.BSs[i])
        for k = (i-1)*Kc+1:i*Kc
            sigma2_dBm = get_receiver_noise_power_dBm(network.MSs[k], network)
            SNRs_dB[k] = Pt_dBm - sigma2_dBm - pathloss_beta - pathloss_alpha*log10(get(network.MS_serving_BS_distance))
        end
    end

    return SNRs_dB
end

get_average_SNRs{MS_t <: PhysicalMS, BS_t <: PhysicalBS, System_t <: System}(network::RandomLargeScaleNetwork{MS_t, BS_t, System_t, SimpleLargescaleFadingEnvironment}) =
    10.^(get_average_SNRs_dB(network)/10)

function set_average_SNRs_dB!{MS_t <: PhysicalMS, BS_t <: PhysicalBS, System_t <: System}(network::RandomLargeScaleNetwork{MS_t, BS_t, System_t, SimpleLargescaleFadingEnvironment}, SNR_dB)
    isnull(network.MS_serving_BS_distance) && error("Can only call set_average_SNRs_dB! if BS-MS distance for served MSs is specified.")

    I = get_num_BSs(network); Kc = network.num_MSs_per_cell
    pathloss_alpha = network.propagation_environment.pathloss_alpha
    pathloss_beta = network.propagation_environment.pathloss_beta

    # Set transmit power of each BS based on the noise characteristics of its first receiving MS
    for i = 1:I
        k1 = (i-1)*Kc+1
        sigma2_dBm = get_receiver_noise_power_dBm(network.MSs[k1], network)
        Pt_dBm = SNR_dB + sigma2_dBm + pathloss_beta + pathloss_alpha*log10(get(network.MS_serving_BS_distance))
        set_transmit_power_dBm!(network.BSs[i], Pt_dBm)
    end
end

set_average_SNRs!{MS_t <: PhysicalMS, BS_t <: PhysicalBS, System_t <: System}(network::RandomLargeScaleNetwork{MS_t, BS_t, System_t, SimpleLargescaleFadingEnvironment}, SNR) =
    set_average_SNR_dB!(network, 10*log10(SNR))

##########################################################################
# Cell assignment functions

# Standard cell assignment functions, only applicable if MS_serving_BS_distance is not null
function IDCellAssignment!(channel, network::RandomLargeScaleNetwork)
    isnull(network.MS_serving_BS_distance) && error("Can only call IDCellAssignment! if BS-MS distance for served MSs is specified.")

    I = get_num_BSs(network); Kc = network.num_MSs_per_cell
    cell_assignment = Array(Int, I*Kc)

    for i = 1:I
        cell_assignment[(i-1)*Kc+1:i*Kc] = i
    end

    network.assignment = Assignment(cell_assignment, I)

    return AssignmentResults()
end

function LargeScaleFadingCellAssignment!(channel, network::RandomLargeScaleNetwork)
    if isnull(network.MS_serving_BS_distance)
        # Greedy scheduler based on the large scale fading realizations
        I = get_num_BSs(network); K = get_num_MSs(network)

        aux_params = get_aux_assignment_params(network)
        @defaultize_param! aux_params "max_num_MSs_per_BS" 1

        # Scheduling matrix
        cell_assignment_matrix = zeros(Int, K, I)

        # User selection metric
        F = (channel.large_scale_fading_factor.^2)*Diagonal(get_transmit_powers(network))
        Fsize = size(F)

        # Do greedy scheduling
        while !all(F .== 0.)
            _, idx = findmax(F)
            k, l = ind2sub(Fsize, idx)

            if sum(cell_assignment_matrix[:,l]) < aux_params["max_num_MSs_per_BS"]
                cell_assignment_matrix[k,l] = 1
                F[k,:] = 0.
            else
                F[:,l] = 0.
            end
        end

        network.assignment = Assignment(cell_assignment_matrix)

        return AssignmentResults()
    else
        # MSs served by their pre-specified BS
        return IDCellAssignment!(channel, network)
    end
end

##########################################################################
# Simulation functions
function draw_user_drop!{MS_t <: PhysicalMS, BS_t <: PhysicalBS, System_t <: System}(network::RandomLargeScaleNetwork{MS_t, BS_t, System_t, SimpleLargescaleFadingEnvironment})
    I = get_num_BSs(network); K = get_num_MSs(network)
    Kc = network.num_MSs_per_cell

    if isnull(network.MS_serving_BS_distance)
        # Drop BSs uniformly at random
        for i = 1:I
            BS_x = network.geography_size[1]*rand(); BS_y = network.geography_size[2]*rand()
            network.BSs[i].position = Position(BS_x, BS_y)
        end
        # Drop MSs uniformly at random
        for k = 1:K
            MS_x = network.geography_size[1]*rand(); MS_y = network.geography_size[2]*rand()
            network.MSs[k].position = Position(MS_x, MS_y)
        end
    else
        # Drop BSs uniformly at random. Also drop the correspondingly served MSs.
        for i = 1:I
            BS_x = network.geography_size[1]*rand(); BS_y = network.geography_size[2]*rand()
            network.BSs[i].position = Position(BS_x, BS_y)

            for k_idx = 1:Kc
                k = (i-1)*Kc + k_idx

                # Position
                theta = 2*pi*rand()
                dx = get(network.MS_serving_BS_distance)*cos(theta)
                dy = get(network.MS_serving_BS_distance)*sin(theta)
                network.MSs[k].position = Position(BS_x + dx, BS_y + dy)

                # Shadow fading (uncorrelated between BSs)
                network.MSs[k].propagation_environment_state =
                    SimpleLargescaleFadingEnvironmentState(network.propagation_environment.shadow_sigma_dB*randn(I), falses(I)) # LoS not used in this model
            end
        end
    end
end

function draw_channel{MS_t <: PhysicalMS, BS_t <: PhysicalBS}(network::RandomLargeScaleNetwork{MS_t, BS_t, SinglecarrierSystem, SimpleLargescaleFadingEnvironment})
    K = get_num_MSs(network); I = get_num_BSs(network)
    Ns = get_num_MS_antennas(network); Ms = get_num_BS_antennas(network)

    coefs = Array(Matrix{Complex128}, K, I)
    large_scale_fading_factor = Array(Float64, K, I)

    distances = get_distances(network)
    angles = get_angles(network)

    for k = 1:K
        for i = 1:I
            # Small scale fading
            coefs[k,i] = (1/sqrt(2))*(randn(Ns[k], Ms[i])
                                 + im*randn(Ns[k], Ms[i]))

            # Pathloss
            pathloss_factor = sqrt(10^(-(network.propagation_environment.pathloss_beta + network.propagation_environment.pathloss_alpha*log10(distances[k,i]))/10))

            # Shadow fading
            shadow_factor = sqrt(10^(network.MSs[k].propagation_environment_state.shadow_realizations_dB[i]/10))

            # Penetration loss
            penetration_loss_factor = sqrt(10^(-network.propagation_environment.penetration_loss_dB/10))

            # BS antenna gain
            bs_antenna_gain = sqrt(get_antenna_gain(network.BSs[i].antenna_params, angles[k,i]))

            # MS antenna gain
            ms_antenna_gain = sqrt(10^(network.MSs[k].antenna_gain_dB/10))

            # Apply scale factors
            large_scale_fading_factor[k,i] = pathloss_factor*shadow_factor*penetration_loss_factor*bs_antenna_gain*ms_antenna_gain
            coefs[k,i] *= large_scale_fading_factor[k,i]
        end
    end

    return SinglecarrierChannel(coefs, Ns, Ms, K, I, large_scale_fading_factor)
end

##########################################################################
# Visualization functions
function plot_network_layout(network::RandomLargeScaleNetwork)
    I = get_num_BSs(network); K = get_num_MSs(network)

    fig = PyPlot.figure()
    ax = fig[:add_subplot](1, 1, 1)

    # Rectangle edges
    Mx = network.geography_size[1]; My = network.geography_size[2]
    ax[:plot]([0, Mx], [0, 0], color="b", linestyle="-")
    ax[:plot]([0, Mx], [My, My], color="b", linestyle="-")
    ax[:plot]([0, 0], [0, My], color="b", linestyle="-")
    ax[:plot]([Mx, Mx], [0, My], color="b", linestyle="-")

    # MSs
    for i = 1:I
        pos = network.BSs[i].position
        ax[:plot](pos.x, pos.y; marker="x", color="b", markersize=6)
        ax[:text](pos.x, pos.y, "BS $i")
    end

    # MSs
    for k = 1:K
        pos = network.MSs[k].position
        ax[:plot](pos.x, pos.y; marker="o", color="r", markersize=6)
        ax[:text](pos.x, pos.y, "MS $k")
    end

    display(fig)
end


text:
module TestCompositionsBase

using CompositionsBase
using CompositionsBase: @var_str
using Test

@testset "binary" begin
    @test tuple ∘ inv ===
          compose(tuple, inv) ===
          var"⨟"(inv, tuple) ===
          opcompose(inv, tuple)
end

if VERSION >= v"1.5"  # TODO: remove this once Compact adds it
@testset "unary" begin
    @test ∘(tuple) === compose(tuple) === var"⨟"(tuple) === opcompose(tuple) === tuple
end
end  # if

if VERSION >= v"1.5"
    include("test_julia15.jl")
end

end  # module


text:
using CFITSIO
using Test

# `create_test_file` : Create a simple FITS file for testing, with the
# given header string added after the required keywords. The length of
# `header` must be a multiple of 80.  The purpose of creating such
# files is to test the parsing of non-standard FITS keyword records
# (non-standard files can't be created with cfitsio).

function create_test_file(fname::AbstractString, header::String)
    if length(header) % 80 != 0
        error("length of header must be multiple of 80")
    end

    f = open(fname, "w")

    stdhdr = "SIMPLE  =                    T / file does conform to FITS standard             BITPIX  =                  -64 / number of bits per data pixel                  NAXIS   =                    2 / number of data axes                            NAXIS1  =                   10 / length of data axis 1                          NAXIS2  =                   10 / length of data axis 2                          EXTEND  =                    T / FITS dataset may contain extensions            "
    endline = "END                                                                             "
    data = fill(0., (10, 10))  # 10x10 array of big-endian Float64 zeros

    # write header
    write(f, stdhdr)
    write(f, header)
    write(f, endline)

    # add padding
    block_position = (length(stdhdr) + length(header) + length(endline)) % 2880
    padding = (block_position == 0) ? 0 : 2880 - block_position
    write(f, " "^padding)

    # write data
    write(f, data)

    # add padding
    block_position = sizeof(data) % 2880
    padding = (block_position == 0) ? 0 : 2880 - block_position
    write(f, fill(0x00, (padding,)))

    close(f)
end

function writehealpix(filename, pixels, nside, ordering, coordsys)
    if eltype(pixels) == Float32
        tform = "1E"
    elseif eltype(pixels) == Float64
        tform = "1D"
    end

    file = fits_create_file("!"*filename)
    try
        fits_create_img(file, Int16, Int[])
        fits_write_date(file)
        fits_movabs_hdu(file, 1)
        fits_create_binary_tbl(file, length(pixels), [("SIGNAL", tform, "")], "BINTABLE")
        fits_write_key(file, "PIXTYPE",  "HEALPIX", "HEALPIX pixelization")
        fits_write_key(file, "ORDERING", ordering,  "Pixel ordering scheme (either RING or NESTED)")
        fits_write_key(file, "NSIDE",    nside,     "Resolution parameter for HEALPIX")
        fits_write_key(file, "COORDSYS", coordsys,  "Pixelization coordinate system")
        fits_write_comment(file, "G = galactic, E = ecliptic, C = celestial = equatorial")
        fits_write_col(file, 1, 1, 1, pixels)
    finally
        fits_close_file(file)
    end
end

function readhealpix(filename)
    file = fits_open_file(filename)
    try
        hdutype = fits_movabs_hdu(file, 2)
        tform, tform_comment = fits_read_key_str(file, "TFORM1")
        if tform == "1E"
            T = Float32
        elseif tform == "1D"
            T = Float64
        end

        naxes, naxis_comment = fits_read_key_lng(file, "NAXIS")
        naxis, nfound = fits_read_keys_lng(file, "NAXIS", 1, naxes)
        nside, nside_comment = fits_read_key_lng(file, "NSIDE")
        npix = 12*nside*nside

        ordering, ordering_comment = fits_read_key_str(file, "ORDERING")
        coordsys, coordsys_comment = fits_read_key_str(file, "COORDSYS")

        pixels = zeros(T, npix)
        fits_read_col(file, 1, 1, 1, pixels)

        return pixels, nside, ordering, coordsys

    finally
        fits_close_file(file)
    end
end

@testset "CFITSIO.jl" begin

    @testset "types" begin
        for (T, code) in (
            (UInt8, 11),
            (Int8, 12),
            (Bool, 14),
            (String, 16),
            (Cushort, 20),
            (Cshort, 21),
            (Cuint, 30),
            (Cint, 31),
            (Int64, 81),
            (Float32, 42),
            (Float64, 82),
            (ComplexF32, 83),
            (ComplexF64, 163),
        )
            @test cfitsio_typecode(T) == Cint(code)
        end

        for (T, code) in ((UInt8,     8), # BYTE_IMG
                        (Int16,    16), # SHORT_IMG
                        (Int32,    32), # LONG_IMG
                        (Int64,    64), # LONGLONG_IMG
                        (Float32, -32), # FLOAT_IMG
                        (Float64, -64), # DOUBLE_IMG
                        (Int8,     10), # SBYTE_IMG
                        (UInt16,   20), # USHORT_IMG
                        (UInt32,   40)) # ULONG_IMG
            @test bitpix_from_type(T) == code
            @test type_from_bitpix(Cint(code)) == type_from_bitpix(Val(Cint(code))) == T
            @test type_from_bitpix(Int16(code)) == T
        end
        @test_throws MethodError type_from_bitpix(7)
        @test_throws MethodError type_from_bitpix("BITPIX")
    end

    # test reading/writing Healpix maps as FITS binary tables using the Libcfitsio interface
    for T in (Float32, Float64)
        nside = 4
        npix = 12*nside*nside
        pixels = rand(T, npix)
        ordering = "NESTED"
        coordsys = "G"

        filename = tempname() * ".fits"
        try
            writehealpix(filename, pixels, nside, ordering, coordsys)
            @test readhealpix(filename) == (pixels, nside, ordering, coordsys)
        finally
            ispath(filename) && rm(filename)
        end
    end

    @testset "Miscellaneous" begin
        # test that this function works and returns the right type.
        @test typeof(libcfitsio_version()) === VersionNumber
        # test it parses a number as intended.
        @test libcfitsio_version(3.341)  === VersionNumber(3, 34, 1)
        @test libcfitsio_version(3.41f0) === VersionNumber(3, 41, 0)
    end

end


text:
module KrokiTest

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

using Kroki:
  @mermaid_str,
  @plantuml_str,
  Diagram,
  InvalidDiagramSpecificationError,
  InvalidOutputFormatError,
  StatusError, # Imported from HTTP through Kroki
  render
using Kroki.Service: setEndpoint!

testRenderError(
  title::AbstractString,
  diagram_type::Symbol,
  content::AbstractString,
  output_format::AbstractString,
  expected_error_type::DataType,
) = @testset "$(title)" begin
  diagram = Diagram(diagram_type, content)

  @test_throws(expected_error_type, render(diagram, output_format))

  @testset "rendering" begin
    service_response = "$(title) response"

    rendered_error = sprint(showerror, expected_error_type(service_response, diagram))

    @test occursin("Kroki service responded with:", rendered_error)
    @test occursin(content, rendered_error)
    @test occursin(service_response, rendered_error)
  end
end

function testShowMethodRenders(
  diagram::Diagram,
  mime_type::AbstractString,
  render_output_format::AbstractString,
)
  @test sprint(show, mime_type, diagram) == String(render(diagram, render_output_format))
end

@testset "Kroki" begin
  @testset "`render`" begin
    # This is not an exhaustive list of supported diagram types or output
    # formats, but serves to verify generic rendering logic is available for at
    # least the specified types and output formats
    DIAGRAM_EXAMPLES = [
      :PlantUML => "Kroki -> Julia: Hello",
      :Graphviz => "digraph D { Hello->World }",
      :ditaa => """
      +--------+     +---------+
      |        |     |         |
      |  User  | <-> | Another |
      |        |     |    User |
      |        |     |         |
      +--------+     +---------+
      """,
      # The specification of the diagram type is case-insensitive
      :plantuml => "Kroki <- Julia: Hello",
    ]

    # The PNG specification defines the structure of PNG files (see
    # http://www.libpng.org/pub/png/spec/1.2/PNG-Structure.html). The most
    # straight-forward parts to test for are the header and the end-of-file
    # marker (i.e. the 'IEND image trailer'), both consisting of 8 bytes
    PNG_HEADER = [0x89, 0x50, 0x4e, 0x47, 0x0d, 0x0a, 0x1a, 0x0a]
    PNG_EOF = [0x0, 0x49, 0x45, 0x4e, 0x44, 0xae, 0x42, 0x60, 0x82]
    @testset "$(diagram_format) to PNG" for (diagram_format, specification) in
                                            DIAGRAM_EXAMPLES
      rendered = render(Diagram(diagram_format, specification), "png")

      @test rendered[1:length(PNG_HEADER)] == PNG_HEADER
      @test rendered[(end - 8):end] == PNG_EOF
    end

    @testset "$(diagram_format) to SVG" for (diagram_format, specification) in
                                            DIAGRAM_EXAMPLES
      rendered = String(render(Diagram(diagram_format, specification), "svg"))

      @test startswith(rendered, "<?xml")
      # Some renderers (e.g. Graphviz) include additional whitespace/newlines
      # after the render, these should be ignored when matching
      @test endswith(rendered, r"</svg>\s?")
    end

    @testset "`RenderError`" begin
      # `RenderError`s are thrown from the `render` function whenever an error
      # occurs. Some of these benefit from rewrites into more descriptive
      # errors
      @testset "transforms rewritable `HTTP.ExceptionRequest.StatusError`s" begin
        # Errors in the Kroki service are thrown as
        # `HTTP.ExceptionRequest.StatusError` and should be rewritten in more
        # descriptive errors
        testRenderError(
          "invalid diagram specification",
          :PlantUML,
          # A missing `>` and message cause the specification to be invalid
          "Julia - Kroki:",
          "svg",
          InvalidDiagramSpecificationError,
        )

        testRenderError(
          "invalid output format",
          :Mermaid,
          # The Mermaid renderer does not support PNG output
          "graph TD; A --> B;",
          "png",
          InvalidOutputFormatError,
        )
      end

      @testset "passes unknown `HTTP.ExceptionRequest.StatusError`s as-is" begin
        # Any HTTP related errors that are not due to rendering errors in the
        # Kroki service (e.g. unknown endpoints), should be thrown from
        # `RenderError` as-is
        @test_throws(StatusError, render(Diagram(:non_existent_diagram_type, ""), "svg"))
      end

      @testset "passes other errors as-is" begin
        # Non-`StatusError`s (e.g. `IOError`s due to incorrect hostnames should
        # be thrown/returned as-is
        expected_service_host = "http://localhost:1"
        expected_diagram_type = :plantuml
        setEndpoint!(expected_service_host)

        try
          render(Diagram(expected_diagram_type, "A -> B: C"), "svg")
        catch exception
          rendered_buffer = sprint(showerror, exception)

          @test occursin("ECONNREFUSED", rendered_buffer)
          @test occursin(
            "$(expected_service_host)/$(expected_diagram_type)",
            rendered_buffer,
          )
        end

        setEndpoint!()
      end
    end
  end

  @testset "string literal shorthand for diagram types" begin
    # String literal macros should be defined as a convenient way for
    # specifying diagrams.
    #
    # This is not an exhaustive test as these are dynamically generated. The
    # basic functionality is only verified for key diagram types.
    shorthand_plantuml = plantuml"A -> B: C"
    longhand_plantuml = Diagram(:PlantUML, "A -> B: C")
    # The leading newlines make sure the alignment of the plain text
    # representations is identical across both calling methods
    @test "\n$shorthand_plantuml" == "\n$longhand_plantuml"

    shorthand_mermaid = mermaid"graph TD; A --> B"
    longhand_mermaid = Diagram(:mermaid, "graph TD; A --> B")
    @test shorthand_mermaid == longhand_mermaid
  end

  @testset "`Base.show`" begin
    # Svgbob diagrams only support SVG output. Any other formats should throw
    # `InvalidOutputFormatError`s when called directly.
    #
    # To prevent compatible `AbstractDisplay`s from trying to render
    # incompatible diagram types in certain formats resulting in errors,
    # `Base.showable` should be overridden to indicate the diagram cannot be
    # rendered in the specified MIME type
    svgbob_diagram = Diagram(:svgbob, "-->[_...__... ]")
    @test_throws(InvalidOutputFormatError, sprint(show, "image/png", svgbob_diagram))
    @test !showable("image/png", svgbob_diagram)
    testShowMethodRenders(svgbob_diagram, "image/svg+xml", "svg")

    plantuml_diagram = Diagram(:PlantUML, "A -> B: C")
    testShowMethodRenders(plantuml_diagram, "image/png", "png")
    testShowMethodRenders(plantuml_diagram, "image/svg+xml", "svg")

    @testset "`text/plain`" begin
      # PlantUML diagrams can be rendered nicely in text/plain based
      # environments
      testShowMethodRenders(plantuml_diagram, "text/plain", "utxt")

      # Other diagram types should simply display their `specification`
      @test sprint(show, "text/plain", svgbob_diagram) == svgbob_diagram.specification
    end
  end

  include("./kroki/service_test.jl")
end

end


text:
# -------------------------------- #
# Custom JSON output for our types #
# -------------------------------- #
JSON.lower(a::HasFields) = a.fields
JSON.lower(c::Cycler) = c.vals
_json_lower(x) = JSON.lower(x)
_json_lower(x::Union{Bool,String,Number,Nothing,Missing}) = x
_json_lower(x::Union{Tuple,AbstractArray}) = _json_lower.(x)
_json_lower(d::Dict) = Dict{Any,Any}(k => _json_lower(v) for (k, v) in pairs(d))
_json_lower(a::HasFields) = Dict{Any,Any}(k => _json_lower(v) for (k, v) in pairs(a.fields))
_json_lower(c::Cycler) = c.vals
function _json_lower(c::ColorScheme)::Vector{Tuple{Float64,String}}
    N = length(c.colors)
    map(ic -> ((ic[1] - 1) / (N - 1), _json_lower(ic[2])), enumerate(c.colors))
end

_maybe_set_attr!(hf::HasFields, k::Symbol, v::Any) =
    get(hf, k, nothing) == nothing && setindex!(hf, v, k)

# special case for associative to get nested application
function _maybe_set_attr!(hf::HasFields, k1::Symbol, v::AbstractDict)
    for (k2, v2) in v
        _maybe_set_attr!(hf, Symbol(k1, "_", k2), v2)
    end
end

function _maybe_set_attr!(p::Plot, k1::Symbol, v::AbstractDict)
    for (k2, v2) in v
        _maybe_set_attr!(p, Symbol(k1, "_", k2), v2)
    end
end

function _maybe_set_attr!(p::Plot, k::Symbol, v)
    foreach(t -> _maybe_set_attr!(t, k, v), p.data)
end

function _maybe_set_attr!(p::Plot, k::Symbol, v::Cycler)
    ix = 0
    for t in p.data
        if t[k] == Dict()  # was empty
            t[k] = v[ix += 1]
        end
    end
end

function JSON.lower(p::Plot)
    out = Dict(
        :data => _json_lower(p.data),
        :layout => _json_lower(p.layout),
        :frames => _json_lower(p.frames),
        :config => _json_lower(p.config)
    )

    if templates.default !== "none" && _isempty(get(out[:layout], :template, Dict()))
        out[:layout][:template] = _json_lower(templates[templates.default])
    end
    out
end

# Let string interpolation stringify to JSON format
Base.print(io::IO, a::Union{Shape,GenericTrace,PlotlyAttribute,Layout,Plot}) = print(io, JSON.json(a))
Base.print(io::IO, a::Vector{T}) where {T <: GenericTrace} = print(io, JSON.json(a))

GenericTrace(d::AbstractDict{Symbol}) = GenericTrace(pop!(d, :type, "scatter"), d)
GenericTrace(d::AbstractDict{T}) where {T <: AbstractString} = GenericTrace(_symbol_dict(d))
Layout(d::AbstractDict{T}) where {T <: AbstractString} = Layout(_symbol_dict(d))

function JSON.parse(::Type{Plot}, str::AbstractString)
    d = JSON.parse(str)
    data = GenericTrace[GenericTrace(tr) for tr in d["data"]]
    layout = Layout(d["layout"])
    Plot(data, layout)
end

JSON.parsefile(::Type{Plot}, fn) =
    open(fn, "r") do f; JSON.parse(Plot, String(read(f))) end


text:
const catpool_seed = UInt === UInt32 ? 0xe3cf1386 : 0x356f2c715023f1a5

hashlevels(levs::AbstractVector) = foldl((h, x) -> hash(x, h), levs, init=catpool_seed)

CategoricalPool{T, R, V}(ordered::Bool=false) where {T, R, V} =
    CategoricalPool{T, R, V}(T[], ordered)
CategoricalPool{T, R}(ordered::Bool=false) where {T, R} =
    CategoricalPool{T, R}(T[], ordered)
CategoricalPool{T}(ordered::Bool=false) where {T} =
    CategoricalPool{T, DefaultRefType}(T[], ordered)

CategoricalPool{T, R}(levels::Vector, ordered::Bool=false) where {T, R} =
    CategoricalPool{T, R, CategoricalValue{T, R}}(convert(Vector{T}, levels), ordered)
CategoricalPool(levels::Vector{T}, ordered::Bool=false) where {T} =
    CategoricalPool{T, DefaultRefType}(convert(Vector{T}, levels), ordered)

CategoricalPool(invindex::Dict{T, R}, ordered::Bool=false) where {T, R <: Integer} =
    CategoricalPool{T, R, CategoricalValue{T, R}}(invindex, ordered)

Base.convert(::Type{T}, pool::T) where {T <: CategoricalPool} = pool

Base.convert(::Type{CategoricalPool{S}}, pool::CategoricalPool{T, R}) where {S, T, R <: Integer} =
    convert(CategoricalPool{S, R}, pool)

function Base.convert(::Type{CategoricalPool{T, R}}, pool::CategoricalPool) where {T, R <: Integer}
    if length(levels(pool)) > typemax(R)
        throw(LevelsException{T, R}(levels(pool)[typemax(R)+1:end]))
    end

    levelsT = convert(Vector{T}, pool.levels)
    invindexT = convert(Dict{T, R}, pool.invindex)
    return CategoricalPool{T, R, CategoricalValue{T, R}}(levelsT, invindexT, pool.ordered)
end

Base.copy(pool::CategoricalPool{T, R, V}) where {T, R, V} =
    CategoricalPool{T, R, V}(copy(pool.levels), copy(pool.invindex),
                             pool.ordered, pool.hash)

function Base.show(io::IO, pool::CategoricalPool{T, R}) where {T, R}
    @static if VERSION >= v"1.6.0"
        @printf(io, "%s{%s, %s}([%s])", CategoricalPool, T, R,
                join(map(repr, levels(pool)), ", "))
    else
        @printf(io, "%s{%s,%s}([%s])", CategoricalPool, T, R,
                join(map(repr, levels(pool)), ", "))
    end

    pool.ordered && print(io, " with ordered levels")
end

Base.length(pool::CategoricalPool) = length(pool.levels)

Base.getindex(pool::CategoricalPool, i::Integer) = CategoricalValue(pool, i)
Base.get(pool::CategoricalPool, level::Any) = pool.invindex[level]
Base.get(pool::CategoricalPool, level::Any, default::Any) = get(pool.invindex, level, default)

"""
after calling this method, add the returned value to pool.invindex:
it doesn't do this itself to avoid doing a dict lookup twice
"""
@inline function push_level!(pool::CategoricalPool{T, R}, level) where {T, R}
    x = convert(T, level)
    n = length(pool)
    if n >= typemax(R)
        throw(LevelsException{T, R}([level]))
    end

    i = R(n + 1)
    push!(pool.levels, x)
    if pool.hash !== nothing
        pool.hash = hash(x, pool.hash)
    end
    pool.equalto = C_NULL
    pool.subsetof = C_NULL
    i
end

"""
    mergelevels(ordered::Bool, levels::AbstractVector...) -> (vec::Vector, ordered′::Bool)

Merge vectors of values `levels` and return:
- `vec`: a superset of all values in `levels`, respecting orders of values
  in each vector of levels if possible
- `ordered′`: if `ordered=true`, whether order comparisons between all pairs
  of levels in `vec` have a defined result based on orders of values in input `levels`
"""
function mergelevels(ordered::Bool, levels::AbstractVector...)
    T = cat_promote_eltype(levels...)
    res = Vector{T}(undef, 0)

    nonempty_lv = findfirst(!isempty, levels)
    if nonempty_lv === nothing
        # no levels
        return res, ordered
    elseif all(l -> isempty(l) || l == levels[nonempty_lv], levels)
        # Fast path if all non-empty levels are equal
        append!(res, levels[nonempty_lv])
        return res, ordered
    end

    union!(res, levels...)

    # Ensure that relative orders of levels are preserved if possible,
    # giving priority to the first sets of levels in case of conflicts
    for levs in reverse(levels)
        levelsmap = indexin(res, levs)

        # Do not touch levels from other sets at the end
        n = length(res)
        @inbounds for i in length(levelsmap):-1:1
            levelsmap[i] === nothing || break
            levelsmap[i] = n
        end
        j = 1
        @inbounds for i in 1:length(levelsmap)
            if levelsmap[i] === nothing
                levelsmap[i] = j
            else
                j = levelsmap[i] + 1
            end
        end
        permute!(res, sortperm(levelsmap, alg=Base.Sort.MergeSort))
    end

    # Check that result is ordered
    if ordered
        levelsmaps = [indexin(res, l) for l in levels]

        # Check that each original order is preserved
        for m in levelsmaps
            issorted(Iterators.filter(x -> x != nothing, m)) || return res, false
        end

        # Check that all order relations between pairs of subsequent elements
        # are defined in at least one set of original levels
        pairs = fill(false, length(res)-1)
        for m in levelsmaps
            @inbounds for i in eachindex(pairs)
                pairs[i] |= (m[i] != nothing) & (m[i+1] != nothing)
            end
            all(pairs) && return res, true
        end
    end

    res, false
end

@inline function Base.get!(pool::CategoricalPool, level::Any)
    get!(pool.invindex, level) do
        ordered!(pool, false)
        push_level!(pool, level)
    end
end

@inline function Base.get!(pool::CategoricalPool, level::CategoricalValue)
    if pool === level.pool || pool == level.pool
        return refcode(level)
    end
    if level.pool ⊈ pool
        newlevs, ordered = merge_pools(pool, level.pool)
        levels!(pool, newlevs)
        ordered!(pool, ordered)
    end
    get!(pool, unwrap(level))
end

@inline function Base.push!(pool::CategoricalPool, level)
    get!(pool.invindex, level) do
        push_level!(pool, level)
    end
    return pool
end

# TODO: optimize for multiple additions
function Base.append!(pool::CategoricalPool, levels)
    for level in levels
        push!(pool, level)
    end
    return pool
end

# Do not override Base.merge as for internal use we need to use the type and orderedness
# of the first pool rather than promoting both pools
function merge_pools(a::CategoricalPool{T}, b::CategoricalPool) where {T}
    if length(a) == 0 && length(b) == 0
        newlevs = T[]
        ordered = isordered(a)
    elseif length(a) == 0
        newlevs = Vector{T}(levels(b))
        ordered = isordered(b)
    elseif length(b) == 0
        newlevs = copy(levels(a))
        ordered = isordered(a)
    else
        ordered = isordered(a) && (isordered(b) || b ⊆ a)
        nl, ordered = mergelevels(ordered, a.levels, b.levels)
        newlevs = convert(Vector{T}, nl)
    end
    newlevs, ordered
end

@inline function Base.hash(pool::CategoricalPool, h::UInt)
    if pool.hash === nothing
        pool.hash = hashlevels(levels(pool))
    end
    hash(pool.hash, h)
end

@inline function Base.:(==)(a::CategoricalPool, b::CategoricalPool)
    pa = pointer_from_objref(a)
    pb = pointer_from_objref(b)
    # Checking both ways is needed to detect changes to a or b
    if a === b || (a.equalto == pb && b.equalto == pa)
        return true
    else
        if hash(a) == hash(b) && isequal(a.levels, b.levels)
            a.equalto = pb
            b.equalto = pa
            return true
        else
            return false
        end
    end
end

# Efficient equivalent of issubset(levels(a), levels(b)), i.e. ignoring order
function Base.issubset(a::CategoricalPool, b::CategoricalPool)
    pa = pointer_from_objref(a)
    pb = pointer_from_objref(b)
    # Checking both ways is needed to detect changes to a or b
    if a === b || (a.equalto == pb && b.equalto == pa) || a.subsetof == pb
        return true
    else
        # Equality check is faster than subset check so do it first
        if hash(a) == hash(b) && isequal(a.levels, b.levels)
            a.equalto = pb
            b.equalto = pa
            return true
        elseif issubset(a.levels, keys(b.invindex))
            a.subsetof = pb
            return true
        else
            return false
        end
    end
end

# Contrary to the CategoricalArray one, this method only allows adding new levels at the end
# so that existing CategoricalValue objects still point to the same value
function levels!(pool::CategoricalPool{S, R}, newlevels::Vector;
                 checkunique::Bool=true) where {S, R}
    levs = convert(Vector{S}, newlevels)
    if checkunique && !allunique(levs)
        throw(ArgumentError(string("duplicated levels found in levs: ",
                                   join(unique(filter(x->sum(levs.==x)>1, levs)), ", "))))
    elseif length(levs) < length(pool) || view(levs, 1:length(pool)) != pool.levels
        throw(ArgumentError("removing or reordering levels of existing CategoricalPool is not allowed"))
    end

    n = length(levs)

    if n > typemax(R)
        throw(LevelsException{S, R}(setdiff(levs, levels(pool))[typemax(R)-length(levels(pool))+1:end]))
    end

    empty!(pool.invindex)
    resize!(pool.levels, n)
    pool.hash = nothing
    pool.equalto = C_NULL
    pool.subsetof = C_NULL
    for i in 1:n
        v = levs[i]
        pool.levels[i] = v
        pool.invindex[v] = i
    end

    return pool
end

DataAPI.levels(pool::CategoricalPool) = pool.levels

isordered(pool::CategoricalPool) = pool.ordered
ordered!(pool::CategoricalPool, ordered) = (pool.ordered = ordered; pool)

# LevelsException
function Base.showerror(io::IO, err::LevelsException{T, R}) where {T, R}
    levs = join(repr.(err.levels), ", ", " and ")
    print(io, "cannot store level(s) $levs since reference type $R can only hold $(typemax(R)) levels. Use the decompress function to make room for more levels.")
end


text:
@testset "Transverse Mercator projection and inverse" begin

    # File has 5000 entries of lon between ±50°, lat between ±90° (with UTM scaling 0.9996)
    # Adapted from Charles Karney's test data (http://zenodo.org/record/32470#.VzF-Krp97CI)
    f = open("TMcoords_lite.bin") # Binary file, Float64: lat, lon, x, y, γ, k
    tmdat = read(f, Float64, (5000, 6))
    close(f)

    tm = Geodesy.TransverseMercator(wgs84)

    for i = 1:size(tmdat,1)
        lat = tmdat[i,1]
        lon = tmdat[i,2]
        x = tmdat[i,3]
        y = tmdat[i,4]
        γ = tmdat[i,5]
        k = tmdat[i,6]

        (x2, y2, γ2, k2) = Geodesy.transverse_mercator_forward(0.0, lat, lon, 0.9996, tm)

        @test (x2, y2, γ2, k2) ≈ (x, y, γ, k)

        (lat2, lon2, γ2, k2) = Geodesy.transverse_mercator_reverse(0.0, x, y, 0.9996, tm)

        @test (lat2, lon2, γ2, k2) ≈ (lat, lon, γ, k)
    end
end


text:
"""
`PairwiseListMatrix{T, diagonal}` is a (squared) symmetric matrix that stores a `list` of values of type `T` for the pairwise comparison/evaluation of `nelements`.
If `diagonal` is `true` the first element of the list is `1, 1` otherwise is `1, 2`.
If `diagonal` is `false` the diagonal values are stored in a vector on the `diag` field.
Labels can be stored on the field `labels` as an `IndexedArray`.
"""
type PairwiseListMatrix{T, diagonal} <: AbstractArray{T, 2}
  list::AbstractVector{T} # i.e. It can be a DataArray
  diag::AbstractVector{T}
  labels::IndexedArray
  nelements::Int
end

"Returns `true` if the list has diagonal values."
@inline hasdiagonal{T, diagonal}(::PairwiseListMatrix{T, diagonal}) = diagonal

"Retuns the `list` vector."
@inline getlist(plm::PairwiseListMatrix) = plm.list

"Retuns the `diag` vector (which contains the diagonal values if `diagonal` is `true`)."
@inline getdiag(plm::PairwiseListMatrix) = plm.diag

# Creation
# ========

function _test_nelements(vector, nelements::Int, name::ASCIIString)
  if length(vector) != 0 && length(vector) != nelements
    throw(DimensionMismatch(string(name, " must have ", nelements, " names!")))
  end
end

# listlength
# ==========

"Returns the length of the `list` field"
lengthlist(plm::PairwiseListMatrix) = length(getlist(plm))

lengthlist{diagonal}(nelements::Int, ::Type{Val{diagonal}}) = diagonal ? div(nelements*(nelements+1),2) : div(nelements*(nelements-1),2)

"""Returns the list length needed for a pairwise measures or comparisons of `nelements`.
If `diagonal` is `true`, diagonal values are included in the list.

```
julia> using PairwiseListMatrices

julia> PLM = PairwiseListMatrix([1, 2, 3, 4, 5, 6], false)
4x4 PairwiseListMatrices.PairwiseListMatrix{Int64,false}:
 0  1  2  3
 1  0  4  5
 2  4  0  6
 3  5  6  0

julia> lengthlist(4, false)
6

julia> lengthlist(PLM)
6

```
"""
lengthlist(nelements::Int, diagonal::Bool) = lengthlist(nelements, Val{diagonal})

# Indexes: ij2k
# =============

ij2k{diagonal}(i::Int, j::Int, nelements::Int, ::Type{Val{diagonal}}) = div(diagonal ?  (nelements*(nelements+1))-((nelements-i)*(nelements-i+1)) :
                                                                              (nelements*(nelements-1))-((nelements-i)*(nelements-i-1)), 2) - nelements + j

"""Returns the `k` index of the `list` from the indixes `i` and `j` with `i<j` from a matrix of `nelements` by `nelements`.
`diagonal` should be `true` or `Val{true}` if the diagonal values are on the `list`. You must not use it with `i>j`.

```
julia> PLM = PairwiseListMatrix([10,20,30,40,50,60], true)
3x3 PairwiseListMatrices.PairwiseListMatrix{Int64,true}:
 10  20  30
 20  40  50
 30  50  60

julia> ij2k(1, 2, 3, true)
2

julia> getlist(PLM)[2]
20

```
"""
ij2k(i::Int, j::Int, nelements::Int, diagonal::Bool) = ij2k(i, j, nelements, Val{diagonal})

# Empties
# -------

"""
An empty `PairwiseListMatrix` can be created for a given `Type` and a number of elements `nelements`.
Optionally, the vector or IndexedArray can be used to declare names/labels to each element.
The `diagonal` (default to `false`) can be declared as `true` to indicate that the list needs space for the diagonal elements.
If `diagonal` is `false` the diagonal values are stored in a vector on the `diag` field instead of being on the list.
This vector can be filled with the optional `diagonalvalue` arguments (default to `0`).
"""
function PairwiseListMatrix{T}(::Type{T}, nelements::Int, labels::IndexedArray,
                               diagonal::Bool=false, diagonalvalue::T=zero(T))
  _test_nelements(labels, nelements, "labels")
  if diagonal
    return PairwiseListMatrix{T, diagonal}(Array(T, lengthlist(nelements, Val{true})), T[], labels, nelements)
  else
    return PairwiseListMatrix{T, diagonal}(Array(T, lengthlist(nelements, Val{false})), fill!(Array(T, nelements), diagonalvalue), labels, nelements)
  end
end

PairwiseListMatrix{T}(::Type{T}, nelements::Int, labels::AbstractVector, diagonal::Bool=false,
                      diagonalvalue::T=zero(T)) = PairwiseListMatrix(T, nelements, IndexedArray(labels), diagonal, diagonalvalue)

PairwiseListMatrix{T}(::Type{T}, nelements::Int, diagonal::Bool=false,
                      diagonalvalue::T=zero(T)) = PairwiseListMatrix(T, nelements, IndexedArray{Any}(), diagonal, diagonalvalue)

# Convert
# -------

convert{T, F, D}(::Type{PairwiseListMatrix{T, D}},
                 mat::PairwiseListMatrix{F, D}) = PairwiseListMatrix{T, D}(convert(Vector{T}, mat.list),
                                                                           convert(Vector{T}, mat.diag), copy(mat.labels), copy(mat.nelements))

function convert{T, F}(::Type{PairwiseListMatrix{T, true}}, mat::PairwiseListMatrix{F, false})
  N = copy(mat.nelements)
  ncol = size(mat, 2)
  list = Array(T, lengthlist(N, Val{true}))
  k = 0
  for i in 1:ncol
    for j in i:ncol
      @inbounds list[k += 1] = mat[i,j]
    end
  end
  PairwiseListMatrix{T, true}(list, T[], copy(mat.labels), N)
end

function convert{T, F}(::Type{PairwiseListMatrix{T, false}}, mat::PairwiseListMatrix{F, true})
  N = copy(mat.nelements)
  ncol = size(mat, 2)
  diag = Array(T, N)
  matlist = mat.list
  list = Array(T, length(matlist) - N)
  l = 0
  d = 0
  m = 0
  for i in 1:ncol
    for j in i:ncol
      if i != j
        @inbounds list[l += 1] = matlist[m += 1]
      else
        @inbounds diag[d += 1] = matlist[m += 1]
      end
    end
  end
  PairwiseListMatrix{T, false}(list, diag, copy(mat.labels), N)
end

function convert{T, F, D}(::Type{PairwiseListMatrix{T, D}}, mat::Matrix{F})
  N = size(mat, 1)
  if N != size(mat, 2)
    throw(ErrorException("It should be a squared matrix."))
  end
  plm = PairwiseListMatrix(T, N, D)
  for i in 1:N
    @inbounds for j in i:N
      value = mat[i,j]
      if i != j && value != mat[j,i]
        throw(ErrorException("It should be a symmetric matrix."))
      end
      plm[i,j] = value
    end
  end
  plm
end

# convert{T <: PairwiseListMatrix, F}(::Type{Array{F,2}}, mat::T) = convert(Matrix{F}, full(mat))

# This is faster than list comprehension (2.4 x)
function convert{T, F}(::Type{Array{F,2}}, lm::PairwiseListMatrix{T, true})
  N = lm.nelements
  complete = Array(F, N, N)
  list = lm.list
  k = 0
  for col in 1:N
    @inbounds for row in col:N
      complete[row, col] = complete[col, row] = F(list[k += 1])
    end
  end
  complete
end

function convert{T, F}(::Type{Array{F,2}}, lm::PairwiseListMatrix{T, false})
  N = lm.nelements
  complete = Array(F, N, N)
  list = lm.list
  diag = lm.diag
  k = 0
  l = 0
  for col in 1:(N-1)
    @inbounds for row in (col+1):N
      complete[row, col] = F(list[k += 1])
    end
    @inbounds for row in (col+1):N
      complete[col, row] = F(list[l += 1])
    end
  end
  @inbounds for i in 1:N
    complete[i, i] = F(diag[i])
  end
  complete
end

# full
# ====

"""Returns a full dense matrix.
This converts a `PairwiseListMatrix{T, D}` into a `Matrix{T}`"""
full{T, D}(lm::PairwiseListMatrix{T, D}) = convert(Array{T, 2}, lm)

full{T, S <: PairwiseListMatrix}(m::Symmetric{T, S}) = full(m.data)

# From a list
# -----------

@inline _nelements(len::Int) = @fastmath div(1+Int(sqrt(1+8*len)),2)
@inline _nelements_with_diagonal(len::Int) = @fastmath div(Int(sqrt(1+8*len)-1),2)

"""
A `PairwiseListMatrix` can be created from a `list`. Optionally, the vector or IndexedArray can be used for declare `labels` to each element.
The `diagonal` (default to `false`) could be declared as `true` to indicate that the list has the diagonal elements.
If `diagonal` is `false`, the diagonal values are stored in a vector on the `diag` field instead of being on the list.
This vector can be filled with the optional `diagonalvalue` arguments (default to `0`).
"""
function PairwiseListMatrix{T}(list::AbstractVector{T}, labels::IndexedArray,
                               diagonal::Bool=false, diagonalvalue::T=zero(T))
  if diagonal
    nelements = _nelements_with_diagonal(length(list))
    _test_nelements(labels, nelements, "labels")
    return PairwiseListMatrix{T, diagonal}(list, T[], labels, nelements)
  else
    nelements = _nelements(length(list))
    _test_nelements(labels, nelements, "labels")
    return PairwiseListMatrix{T, diagonal}(list, fill!(Array(T, nelements), diagonalvalue), labels, nelements)
  end
end

PairwiseListMatrix{T}(list::AbstractVector{T}, labels::AbstractVector, diagonal::Bool=false,
                      diagonalvalue::T=zero(T)) = PairwiseListMatrix(list, IndexedArray(labels), diagonal, diagonalvalue)

PairwiseListMatrix{T}(list::AbstractVector{T}, diagonal::Bool=false,
                      diagonalvalue::T=zero(T)) = PairwiseListMatrix(list, IndexedArray{Any}(), diagonal, diagonalvalue)

# AbstractArray methods
# =====================

size(m::PairwiseListMatrix) = (m.nelements, m.nelements)
length(m::PairwiseListMatrix) = m.nelements * m.nelements

eltype{T, diagonal}(m::PairwiseListMatrix{T, diagonal}) = T

similar{T, diagonal}(m::PairwiseListMatrix{T, diagonal}) = PairwiseListMatrix{T, diagonal}(similar(m.list), copy(m.diag),
                                                                                copy(m.labels), copy(m.nelements))
similar{T, diagonal, S}(m::PairwiseListMatrix{T, diagonal}, ::Type{S}) = PairwiseListMatrix{S, diagonal}(similar(m.list, S), similar(m.diag, S),
                                                                                copy(m.labels), copy(m.nelements))

copy{T, diagonal}(m::PairwiseListMatrix{T, diagonal}) = PairwiseListMatrix{T, diagonal}(copy(m.list), copy(m.diag), copy(m.labels), copy(m.nelements))

zeros{T, diagonal}(m::PairwiseListMatrix{T, diagonal}) = PairwiseListMatrix{T, diagonal}(zeros(m.list), zeros(m.diag), copy(m.labels), copy(m.nelements))
ones{T, diagonal}(m::PairwiseListMatrix{T, diagonal})  = PairwiseListMatrix{T, diagonal}(ones(m.list),  ones(m.diag),  copy(m.labels), copy(m.nelements))

# TODO:  p = zeros(PairwiseListMatrix{Float64,false}, 3)

# Indexing (getindex)
# ===================

function getindex{T}(lm::PairwiseListMatrix{T, true}, i::Int, j::Int)
  if i <= j
    return(lm.list[ij2k(i, j, lm.nelements, Val{true})])
  else
    return(lm.list[ij2k(j, i, lm.nelements, Val{true})])
  end
end

function getindex{T}(lm::PairwiseListMatrix{T, false}, i::Int, j::Int)
  if i < j
    return(lm.list[ij2k(i, j, lm.nelements, Val{false})])
  elseif i > j
    return(lm.list[ij2k(j, i, lm.nelements, Val{false})])
  else
    return(lm.diag[i])
  end
end

Base.linearindexing(m::PairwiseListMatrix) = Base.LinearFast()

function getindex{T}(lm::PairwiseListMatrix{T, true}, i::Int)
  n = lm.nelements
  row = Int(ceil(i/n))
  col = i - n*(row-1)
  if row <= col
    @inbounds return lm.list[ij2k(row, col, n, Val{true})]
  else
    @inbounds return lm.list[ij2k(col, row, n, Val{true})]
  end
end

function getindex{T}(lm::PairwiseListMatrix{T, false}, i::Int)
  n = lm.nelements
  row = Int(ceil(i/n))
  col = i - n*(row-1)
  if row < col
    @inbounds return lm.list[ij2k(row, col, n, Val{false})]
  elseif row > col
    @inbounds return lm.list[ij2k(col, row, n, Val{false})]
  else
    @inbounds return lm.diag[row]
  end
end

# Labels
# ======

"Get the labels/names of the row/columns of the matrix"
labels(lm::PairwiseListMatrix) = lm.labels

"You can use labels! for add labels/names to the row/columns of the matrix"
function labels!(lm::PairwiseListMatrix, labels::IndexedArray)
  _test_nelements(labels, lm.nelements, "labels")
  lm.labels = labels
end

labels!(lm::PairwiseListMatrix, labels::AbstractVector) = labels!(lm, IndexedArray(labels))

# Indexing using labels (getlabel)
# ================================

"Like `getindex`, but using the labels/names instead of `Int` numbers."
function getlabel(lm::PairwiseListMatrix, i, j)
  if isempty(lm.labels)
    throw(ErrorException("There are not labels in the matrix. You can use labels!(...) for add them."))
  else
    return getindex(lm, findfirst(lm.labels, i), findfirst(lm.labels, j))
  end
end

# Set values (setindex!)
# ======================

function setindex!{T}(lm::PairwiseListMatrix{T, true}, v, i::Int, j::Int)
  if i <= j
    return setindex!(lm.list, v, ij2k(i, j, lm.nelements, Val{true}))
  else
    return setindex!(lm.list, v, ij2k(j, i, lm.nelements, Val{true}))
  end
end

function setindex!{T}(lm::PairwiseListMatrix{T, false}, v, i::Int, j::Int)
  if i < j
    return setindex!(lm.list, v, ij2k(i, j, lm.nelements, Val{false}))
  elseif i > j
    return setindex!(lm.list, v, ij2k(j, i, lm.nelements, Val{false}))
  else
    return setindex!(lm.diag, v, i)
  end
end

function setindex!{T}(lm::PairwiseListMatrix{T, true}, v, i::Int)
  n = lm.nelements
  row = Int(ceil(i/n))
  col = i - n*(row-1)
  if row <= col
    return setindex!(lm.list, v, ij2k(row, col, n, Val{true}))
  else
    return setindex!(lm.list, v, ij2k(col, row, n, Val{true}))
  end
end

function setindex!{T}(lm::PairwiseListMatrix{T, false}, v, i::Int)
  n = lm.nelements
  row = Int(ceil(i/n))
  col = i - n*(row-1)
  if row < col
    return setindex!(lm.list, v, ij2k(row, col, n, Val{false}))
  elseif row > col
    return setindex!(lm.list, v, ij2k(col, row, n, Val{false}))
  else
    return setindex!(lm.diag, v, row)
  end
end

# Set values using labels (setlabel!)
# ===================================

"Like `setindex!`, but using the labels/names instead of `Int` numbers."
function setlabel!(lm::PairwiseListMatrix, v, i, j)
  if isempty(lm.labels)
    throw(ErrorException("There are not labels in the matrix. You can use labels!(...) for add them."))
  else
    return setindex!(lm, v, findfirst(lm.labels, i), findfirst(lm.labels, j))
  end
end

# Transpose
# =========

transpose(lm::PairwiseListMatrix) = lm
transpose!(lm::PairwiseListMatrix) = lm

ctranspose(lm::PairwiseListMatrix) = lm
ctranspose!(lm::PairwiseListMatrix) = lm

# diag
# ====

diag{T}(lm::PairwiseListMatrix{T, false}) = lm.diag

diag{T}(lm::PairwiseListMatrix{T, true}) = T[ lm[i,i] for i in 1:lm.nelements ]

# Unary operations
# ================

sqrt{T <: Real}(lm::PairwiseListMatrix{T, true}) = PairwiseListMatrix{Float64, true }(sqrt(lm.list), lm.diag, lm.labels, lm.nelements)
sqrt{T <: Real}(lm::PairwiseListMatrix{T, false})= PairwiseListMatrix{Float64, false}(sqrt(lm.list), sqrt(lm.diag), lm.labels, lm.nelements)

for una in (:abs, :-)
  @eval begin
    $(una){T}(lm::PairwiseListMatrix{T, true}) = PairwiseListMatrix{T, true }($(una)(lm.list), lm.diag, lm.labels, lm.nelements)
    $(una){T}(lm::PairwiseListMatrix{T, false})= PairwiseListMatrix{T, false}($(una)(lm.list), $(una)(lm.diag), lm.labels, lm.nelements)
  end
end

svd(m::PairwiseListMatrix) = svd(full(m))

# Binary operations
# =================

for bin in ( :-, :+, :.*, :./, :.+, :.- )

  @eval begin

    function $(bin){T}(A::PairwiseListMatrix{T, true}, B::PairwiseListMatrix{T, true})
      if A.labels != B.labels || A.nelements != B.nelements
        return($(bin)(full(A), full(B)))
      end
      PairwiseListMatrix{T, true}($(bin)(A.list, B.list), copy(A.diag), copy(A.labels), copy(A.nelements))
    end

    function $(bin){T}(A::PairwiseListMatrix{T, false}, B::PairwiseListMatrix{T, false})
      if A.labels != B.labels || A.nelements != B.nelements
        return($(bin)(full(A), full(B)))
      end
      PairwiseListMatrix{T, false}($(bin)(A.list, B.list), $(bin)(A.diag, B.diag), copy(A.labels), copy(A.nelements))
    end

    $(bin)(A::PairwiseListMatrix, B::PairwiseListMatrix) = $(bin)(full(A), full(B))

  end

end

for bin in (:/, :./)

  @eval begin

    $(bin){T <: AbstractFloat}(A::PairwiseListMatrix{T, true},  B::T) = PairwiseListMatrix{T, true }($(bin)(A.list, B), copy(A.diag), copy(A.labels), copy(A.nelements))
    $(bin){T <: AbstractFloat}(A::PairwiseListMatrix{T, false}, B::T) = PairwiseListMatrix{T, false}($(bin)(A.list, B), $(bin)(A.diag, B), copy(A.labels), copy(A.nelements))

    $(bin){T}(A::PairwiseListMatrix{T, true},  B::Integer) = PairwiseListMatrix{Float64, true }($(bin)(A.list, B), copy(A.diag), copy(A.labels), copy(A.nelements))
    $(bin){T}(A::PairwiseListMatrix{T, false}, B::Integer) = PairwiseListMatrix{Float64, false}($(bin)(A.list, B), $(bin)(A.diag, B), copy(A.labels), copy(A.nelements))

  end

end

for bin in (:+, :-)
  @eval begin
    $(bin)(A::PairwiseListMatrix{Bool, true},  B::Bool) = PairwiseListMatrix{Int, true }($(bin)(A.list, B), copy(A.diag), copy(A.labels), copy(A.nelements))
    $(bin)(A::PairwiseListMatrix{Bool, false}, B::Bool) = PairwiseListMatrix{Int, false}($(bin)(A.list, B), $(bin)(A.diag, B), copy(A.labels), copy(A.nelements))
  end
end

for bin in (:.+, :.-, :.*, :-, :+)

  @eval begin

    $(bin){T <: Number}(A::PairwiseListMatrix{T, true},  B::T) = PairwiseListMatrix{T, true }($(bin)(A.list, B), copy(A.diag), copy(A.labels), copy(A.nelements))
    $(bin){T <: Number}(A::PairwiseListMatrix{T, false}, B::T) = PairwiseListMatrix{T, false}($(bin)(A.list, B), $(bin)(A.diag, B), copy(A.labels), copy(A.nelements))

  end

end

for bin in (:*, :/)

  @eval $(bin)(A::PairwiseListMatrix, B::PairwiseListMatrix) = $(bin)(full(A), full(B))

end

# Faster mean and sum
# ===================

sum{T}(m::PairwiseListMatrix{T, false}) = 2*sum(m.list) + sum(m.diag)
sum{T}(m::PairwiseListMatrix{T, true}) =  2*sum(m.list) - sum(diag(m))

function _sum_kernel!(sum_i, list, N)
  k = 0
  l = 0
  for i in 1:N
    l += 1
    @inbounds @simd for j in i:N
      sum_i[i] += list[k += 1]
    end
    @inbounds for j in (i+1):N
      sum_i[j] += list[l += 1]
    end
  end
  sum_i
end

function _sum_kernel!(sum_i, diag, list, N)
  k = 0
  l = 0
  for i in 1:(N-1)
    @inbounds @simd for j in (i+1):N
      sum_i[i] += list[k += 1]
    end
    @inbounds for j in (i+1):N
      sum_i[j] += list[l += 1]
    end
  end
  @inbounds for i in 1:N
    sum_i[i] += diag[i]
  end
  sum_i
end

function _test_and_set_sum{T}(::Type{T}, region, N)
  if region == 1
    sum_i = zeros(T, 1, N)
  elseif region == 2
    sum_i = zeros(T, N, 1)
  else
    throw(ErrorException("region should be 1 or 2"))
  end
  sum_i
end

function sum{T}(lm::PairwiseListMatrix{T, true}, region::Int)
  N = lm.nelements
  sum_i = _test_and_set_sum(T, region, N)
  _sum_kernel!(sum_i, lm.list, N)
end

function sum{T}(lm::PairwiseListMatrix{T, false}, region::Int)
  N = lm.nelements
  sum_i = _test_and_set_sum(T, region, N)
  _sum_kernel!(sum_i, lm.diag, lm.list, N)
end

mean(m::PairwiseListMatrix) = sum(m)/length(m)
mean(m::PairwiseListMatrix, region::Int) = sum(m, region) ./ m.nelements

# Sum/Mean without diagonal: sum/mean_nodiag
# ------------------------------------------

"Sum the values outside the diagonal"
sum_nodiag{T}(m::PairwiseListMatrix{T, false}) = T(2) * sum(m.list)
sum_nodiag{T}(m::PairwiseListMatrix{T, true}) = T(2) * sum(m.list) - sum(diag(m))

function _sum_nodiag_kernel!{T}(sum_i, lm::PairwiseListMatrix{T, true}, N)
  list = lm.list
  k = 0
  for i in 1:N
    for j in i:N
      k += 1
      if i != j
        @inbounds value = list[k]
        @inbounds sum_i[i] += value
        @inbounds sum_i[j] += value
      end
    end
  end
  sum_i
end

function _sum_nodiag_kernel!{T}(sum_i, lm::PairwiseListMatrix{T, false}, N)
  list = lm.list
  k = 0
  l = 0
  for i in 1:(N-1)
    @inbounds @simd for j in (i+1):N
      sum_i[i] += list[k += 1]
    end
    @inbounds for j in (i+1):N
      sum_i[j] += list[l += 1]
    end
  end
  sum_i
end

function sum_nodiag{T, diagonal}(lm::PairwiseListMatrix{T, diagonal}, region::Int)
  N = lm.nelements
  sum_i = _test_and_set_sum(T, region, N)
  _sum_nodiag_kernel!(sum_i, lm, N)
end

"Mean of the values outside the diagonal"
mean_nodiag(m::PairwiseListMatrix) = sum_nodiag(m) / (length(m) - m.nelements)
mean_nodiag(m::PairwiseListMatrix, region::Int) = sum_nodiag(m, region) ./ (m.nelements-1)

# Operations on Vector{PairwiseListMatrix}
# ========================================

# sum
# ---

@inline _has_diagonal{T, diagonal}(x::PairwiseListMatrix{T, diagonal}) = diagonal

function sum{T <: PairwiseListMatrix}(list::Vector{T})
  samples = length(list)
  if samples == 0
    throw(ErrorException("Empty list"))
  end
  if samples == 1
    return list[1]
  else
    start = copy(list[1])
    i = 2
    start_list = start.list
    N = length(start_list)
    i = 2
    while i <= samples
      @inbounds ylist = list[i].list
      if length(ylist) != N
        throw(ErrorException("Different number of elements"))
      end
      for k in 1:N
        @inbounds start_list[k] += ylist[k]
      end
      i += 1
    end
    if !_has_diagonal(start)
      i = 2
      start_diag = start.diag
      while i <= samples
        @inbounds ydiag = list[i].diag
        for k in 1:length(start_diag)
          @inbounds start_diag[k] += ydiag[k]
        end
        i += 1
      end
    end
    return start
  end
end

# std
# ---

function varm{T <: PairwiseListMatrix}(list::Vector{T}, mean::PairwiseListMatrix)
  samples = length(list)
  if samples < 2
    throw(ErrorException("You need at least 2 samples."))
  end
  out = zeros(list[1])
  out_list = out.list
  mean_list = mean.list
  N = length(mean_list)
  @inbounds for sample in list
    sample_list = sample.list
    if length(sample_list) != N
      throw(ErrorException("Different number of elements"))
    end
    for k in 1:N
      @inbounds out_list[k] += abs2( sample_list[k] - mean_list[k])
    end
  end
  if !_has_diagonal(out)
    out_diag = out.diag
    mean_diag = mean.diag
    @inbounds for sample in list
      sample_diag = sample.diag
      for k in 1:length(out_diag)
        @inbounds out_diag[k] += abs2( sample_diag[k] - mean_diag[k])
      end
    end
  end
  out ./ samples
end

function var{T <: PairwiseListMatrix}(list::Vector{T}; mean=nothing)
  mean === nothing ? varm(list, Base.mean(list)) : varm(list, mean)
end

std{T <: PairwiseListMatrix}(list::Vector{T}; mean=nothing) = sqrt(var(list, mean=mean))

# zscore
# ======

function zscore!{T <: PairwiseListMatrix, E <: AbstractFloat, D}(list::Vector{T}, mat::PairwiseListMatrix{E, D})
  list_mean = mean(list)
  if size(mat) != size(list_mean)
    throw(ErrorException("PairwiseListMatrices must have the same size"))
  end
  if D != _has_diagonal(list[1])
    throw(ErrorException("PairwiseListMatrices must have the same diagonal parameter"))
  end
  list_std  = std(list, mean=list_mean)

  mat_list = mat.list
  std_list = list_std.list
  mean_list = list_mean.list

  ElementType = eltype(mat_list)

  @inbounds for i in 1:length(mat_list)
    s = std_list[i]
    m = mean_list[i]
    value = mat_list[i]
    if (m - value + one(ElementType)) ≉ one(ElementType)
      if (s + one(ElementType)) ≉ one(ElementType)
        mat_list[i] = (value - m)/s
      else
        mat_list[i] = NaN
      end
    else
      mat_list[i] = zero(ElementType)
    end
  end

  if _has_diagonal(mat)
    mat_diag = mat.diag
    std_diag = list_std.diag
    mean_diag = list_mean.diag
    @inbounds for i in 1:length(mat_diag)
      s = std_diag[i]
      m = mean_diag[i]
      value = mat_diag[i]
      if (m - value + one(ElementType)) ≉ one(ElementType)
        if (s + one(ElementType)) ≉ one(ElementType)
          mat_diag[i] = (value - m)/s
        else
          mat_diag[i] = NaN
        end
      else
        mat_diag[i] = zero(ElementType)
      end
    end
  end

  mat
end

zscore!{T <: PairwiseListMatrix, E <: Integer, D}(list::Vector{T}, mat::PairwiseListMatrix{E, D}) = zscore!(list, convert(PairwiseListMatrix{Float64, D}, mat))
zscore(list, mat) = zscore!(list, copy(mat))

# Tables
# ======

"""
Creates a `Matrix{Any}` is useful for `writedlm` and/or `writecsv`.
Labels are stored in the columns 1 and 2, and the values in the column 3.
Diagonal values are included by default.

```
julia> list
3x3 PairwiseListMatrices.PairwiseListMatrix{Int64,false}:
 0  1  2
 1  0  3
 2  3  0

julia> to_table(list)
6x3 Array{Any,2}:
 "A"  "A"  0
 "A"  "B"  1
 "A"  "C"  2
 "B"  "B"  0
 "B"  "C"  3
 "C"  "C"  0

julia> to_table(list, false)
3x3 Array{Any,2}:
 "A"  "B"  1
 "A"  "C"  2
 "B"  "C"  3

```
"""
function to_table(plm::PairwiseListMatrix, diagonal::Bool=true)
  N = plm.nelements
  labels = length(plm.labels) != 0 ? plm.labels : IndexedArray(collect(1:N))
  table = Array(Any, diagonal ? div(N*(N+1),2) : div(N*(N-1),2), 3)
  t = 0
  @iterateupper plm diagonal begin
    t += 1
    table[t, 1] = labels[i]
    table[t, 2] = labels[j]
    table[t, 3] = list[k]
  end
  table
end

"""
Creation of a `PairwiseListMatrix` from a `Matrix`.
By default the columns with the labels for i (slow) and j (fast) are 1 and 2.
Values are taken from the column 3 by default.

```
julia> data = readcsv("example.csv")
3x3 Array{Any,2}:
 "A"  "B"  10
 "A"  "C"  20
 "B"  "C"  30

julia> from_table(data, Int, false)
3x3 PairwiseListMatrices.PairwiseListMatrix{Int64,false}:
  0  10  20
 10   0  30
 20  30   0

```
"""
function from_table{T}(table::Matrix, value::Type{T}, diagonal::Bool, labelcols::Vector{Int}=[1,2], valuecol::Int=3)
  PairwiseListMatrix(convert(Vector{T}, table[:,valuecol]), unique(table[:,labelcols]), diagonal)
end

# write...
# ========

"""
This function takes the filename as first argument and a `PairwiseListMatrix` as second argument.
If the third positional argument is `true` (default) the diagonal is included in the output.
The keyword argument `delim` (by default is `'\t’`) allows to modified the character used as delimiter.

julia> PLM  = PairwiseListMatrix(collect(1:6), true)
3x3 PairwiseListMatrices.PairwiseListMatrix{Int64,true}:
 1  2  3
 2  4  5
 3  5  6

julia> writedlm("example.txt", PLM, false, delim=' ')

shell> cat example.txt
1 2 2
1 3 3
2 3 5

"""
function writedlm(filename::ByteString, plm::PairwiseListMatrix, diagonal::Bool=true; delim::Char='\t')
  open(filename, "w") do fh
    if length(plm.labels) != 0
      labels = plm.labels
      @iterateupper plm diagonal println(fh, labels[i], delim, labels[j], delim, list[k])
    else
      @iterateupper plm diagonal println(fh, i, delim, j, delim, list[k])
    end
  end
end

"""
This function takes the filename as first argument and a `PairwiseListMatrix` as second argument.
If the third positional argument is `true` (default) the diagonal is included in the output.

julia> PLM  = PairwiseListMatrix(collect(1:6), true)
3x3 PairwiseListMatrices.PairwiseListMatrix{Int64,true}:
 1  2  3
 2  4  5
 3  5  6

julia> writecsv("example.csv", PLM)

shell> cat example.csv
1,1,1
1,2,2
1,3,3
2,2,4
2,3,5
3,3,6

julia> writecsv("example.csv", PLM, false)

shell> cat example.csv
1,2,2
1,3,3
2,3,5

"""
writecsv(filename::ByteString, plm::PairwiseListMatrix, diagonal::Bool=true) = writedlm(filename, plm, diagonal, delim=',')

# triu! & triu
# ============

triu!(::PairwiseListMatrix, args...) = throw(ErrorException("PairwiseListMatrix must be Symmetric, use triu instead of triu!"))
triu(mat::PairwiseListMatrix, args...) = triu(full(mat), args...)

# join
# ====

"""
This function join two PairwiseListMatrices by their labels, returning two PairwiseListMatrices with same size and labels.
There are 4 `kind` of joins:
  - `:inner` : Intersect. The output matrices only include the labels that are in both PairwiseListMatrices
  - `:outer` : Union. Include the labels of the two PairwiseListMatrices.
  - `:left` : Only use labels form the first argument.
  - `:right` : Only use labels form the second argument.
`NaN`s are filled in where needed to complete joins. The default value for missing values can be changed passing a tuple to `missing`.

```
julia> l = PairwiseListMatrix([1.,2.,3.], ['a', 'b', 'c'], false) # a b c
3x3 PairwiseListMatrices.PairwiseListMatrix{Float64,false}:
 0.0  1.0  2.0
 1.0  0.0  3.0
 2.0  3.0  0.0

julia> r = PairwiseListMatrix([1.,2.,3.], ['b', 'c', 'd'], false) #   b c d
3x3 PairwiseListMatrices.PairwiseListMatrix{Float64,false}:
 0.0  1.0  2.0
 1.0  0.0  3.0
 2.0  3.0  0.0

julia> join(l, r, kind=:inner)                                    #   b c
(
2x2 PairwiseListMatrices.PairwiseListMatrix{Float64,false}:
 0.0  3.0
 3.0  0.0,

2x2 PairwiseListMatrices.PairwiseListMatrix{Float64,false}:
 0.0  1.0
 1.0  0.0)

julia> join(l, r, kind=:outer)                                    # a b c d
(
4x4 PairwiseListMatrices.PairwiseListMatrix{Float64,false}:
   0.0    1.0    2.0  NaN
   1.0    0.0    3.0  NaN
   2.0    3.0    0.0  NaN
 NaN    NaN    NaN    NaN,

4x4 PairwiseListMatrices.PairwiseListMatrix{Float64,false}:
 NaN  NaN    NaN    NaN
 NaN    0.0    1.0    2.0
 NaN    1.0    0.0    3.0
 NaN    2.0    3.0    0.0)

```
"""
function join{L <: AbstractFloat, R <: AbstractFloat, DL, DR}(left::PairwiseListMatrix{L, DL}, right::PairwiseListMatrix{R, DR}; kind::Symbol = :inner ,
                                                              missing::Tuple{L,R} = (L(NaN), R(NaN)) )

  labels_left  = labels(left)
  labels_right = labels(right)

  if labels_left == labels_right
    return(left, right)
  end

  if     kind == :inner
    out_labels = intersect(labels_left, labels_right)
    N = length(out_labels)
    out_L = PairwiseListMatrix(L, N, out_labels, DL)
    out_R = PairwiseListMatrix(R, N, out_labels, DR)
    for i in 1:N
      li = out_labels[i]
      for j in i:N
        lj = out_labels[j]
        out_L[i,j] = getlabel(left,  li, lj)
        out_R[i,j] = getlabel(right, li, lj)
      end
    end
    return(out_L, out_R)
  elseif kind == :left
    out_labels = labels_left
    N = length(out_labels)
    out_R = PairwiseListMatrix(R, N, out_labels, DR)
    for i in 1:N
      li = out_labels[i]
      flag_i_r = li in labels_right
      for j in i:N
        lj = out_labels[j]
        out_R[i,j] = (flag_i_r && (lj in labels_right)) ? getlabel(right, li, lj) : missing[2]
      end
    end
    return(left, out_R)
  elseif kind == :right
    out_labels = labels_right
    N = length(out_labels)
    out_L = PairwiseListMatrix(L, N, out_labels, DL)
    for i in 1:N
      li = out_labels[i]
      flag_i_l = li in labels_left
      for j in i:N
        lj = out_labels[j]
        out_L[i,j] = (flag_i_l  && (lj in labels_left)) ? getlabel(left,  li, lj) : missing[1]
      end
    end
    return(out_L, right)
  elseif kind == :outer
    out_labels = union(labels_left, labels_right)
    N = length(out_labels)
    out_L = PairwiseListMatrix(L, N, out_labels, DL)
    out_R = PairwiseListMatrix(R, N, out_labels, DR)
    for i in 1:N
      li = out_labels[i]
      flag_i_l = li in labels_left
      flag_i_r = li in labels_right
      for j in i:N
        lj = out_labels[j]
        out_L[i,j] = (flag_i_l && (lj in labels_left))  ? getlabel(left,  li, lj) : missing[1]
        out_R[i,j] = (flag_i_r && (lj in labels_right)) ? getlabel(right, li, lj) : missing[2]
      end
    end
    return(out_L, out_R)
  else
    throw(ArgumentError("Unknown kind of join requested: use :inner, :left, :right or :outer"))
  end

end


text:
# Autogenerated wrapper script for HELICS_jll for i686-linux-musl-cxx11
export libhelicsSharedLib

using ZeroMQ_jll
JLLWrappers.@generate_wrapper_header("HELICS")
JLLWrappers.@declare_library_product(libhelicsSharedLib, "libhelicsSharedLib.so.2")
function __init__()
    JLLWrappers.@generate_init_header(ZeroMQ_jll)
    JLLWrappers.@init_library_product(
        libhelicsSharedLib,
        "lib/libhelicsSharedLib.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using Pkg
Pkg.add("JSON")
Pkg.add("JLD")
Pkg.add("Plots")
Pkg.add("ArgParse")
Pkg.add("JuliaFormatter")
Pkg.add("PlotThemes")
Pkg.add("Formatting")
Pkg.add("Colors")
Pkg.add("Base64")
Pkg.add("Dates")

using ArgParse
using Plots
using JSON
using JLD
using Dates
using Base64
using PlotThemes
using Formatting
using Colors

max_stars_repo_name: Agapanthus/server-diary

max_issues_repo_name: Agapanthus/server-diary

max_forks_repo_name: Agapanthus/server-diary


text:
export OptimizationProblem, primalDualSolve
######
# Simple Primal-Dual IP for the problem
#
# min_x f(x)
#   Cx -h <= 0.0

immutable OptimizationProblem
  d :: Int64
  m :: Int64
  f_g_h :: Function
  #R^d |----> R, returns a tuple with function value, gradient, hessian
  C :: Matrix{Float64} #C \in R^(d,m)
  h :: Vector{Float64}
end

# BV page 610
# The newton system is
# [\nabla^2 f_0, C' ]
# [-diag(\lambda) *C, -diag(C*x -h ) ]

function primalDualSolve(p :: OptimizationProblem, x :: Vector{Float64}, lambda :: Vector{Float64}, mu :: Float64, e_feas :: Float64, e :: Float64, max_iters :: Int64, alpha :: Float64, beta :: Float64)
  assert(mu > 1)
  assert(0.0 < beta < 1.0)
  assert(p.d == length(x))
  assert(p.m == length(lambda))
  assert(all(p.C*x .<= p.h))
  assert(all(lambda .> 0.0))
  assert(e > 0.0 && e_feas > 0.0)
  eta_hat(x,lambda) = -dot(lambda,p.C*x-p.h)

  function compute_r_t(f,g,x,lambda,t)
    [g + p.C'*lambda;
    -lambda.*(p.C*x - p.h) - (1.0/t)*ones(p.m)]
  end

  iter = 1
  for iter = 1:max_iters
    # println("ITER :", iter, ". Printing values of x: ", x)
    e_hat = eta_hat(x,lambda)
    t = (p.m*mu)/e_hat
    (f,g,h) = p.f_g_h(x)
    assert(length(g) == p.d)
    assert(size(h) == (p.d,p.d))

    #check termination conditions... Here?
    r_old = norm(compute_r_t(f,g,x,lambda,t))
    if (iter == max_iters)
      println(norm(g + p.C'*lambda))
      println(e_hat)
      println()
    end
    if (norm(g + p.C'*lambda) < e_feas && e_hat < e)
      break;
    end
    #form newton system...
    newton_system = [h p.C';
      -diagm(lambda)*p.C -diagm(p.C*x - p.h)]
    #target
    r_dual = g + p.C'*lambda
    r_cent = -lambda.*(p.C*x - p.h) - (1.0/t)*ones(p.m)
    #solve
    search_direction = -(newton_system\[r_dual; r_cent])
    sd_x = search_direction[1:p.d]
    # println(" Printing values of sd_x:  ", sd_x)
    sd_lambda = search_direction[p.d+1:end]
    #linesearch
    neg_inds = find(sd_lambda .< 0.0) #bullshit
    s_max = 1.0
    if length(neg_inds) > 0
      s_max = min(1.0, minimum(-lambda[neg_inds]./sd_lambda[neg_inds]))
    end
    #hardcoded in BV
    s = 0.995*s_max
    #backtracking for feasibility
    loopingIndex2 = 1
    while(!all(p.C*(x + s*sd_x) .<= p.h ))
      s = beta*s
      if loopingIndex2 ==1000
	  println("########################################################")
	  println("########################################################")
	  println("### WARINING ENTERING INFINITE LOOP, EXIT WITH ERROR ###")
    println("########################################################")
	  println("########################################################")
	  println("########################################################")
	  error()
      end
      loopingIndex2 = loopingIndex2 + 1
    end

    #backtracking for some other thing
    loopingIndex = 1
    while(true)
      if mod(loopingIndex,100) == 0
        println("### loopin' ",loopingIndex," !########")
      end
      loopingIndex = loopingIndex + 1
      x_plus = x + s*sd_x
      lambda_plus = lambda + s*sd_lambda
      f_t,g_t,h_t = p.f_g_h(x_plus) #only need f_g
      #compute norm of r_t
      r_t = compute_r_t(f_t,g_t,x_plus,lambda_plus,t)
      if (vecnorm(r_t) <= (1-alpha*s)*r_old)
        break;
      end
      s = beta*s
    end
    x = x + s*sd_x
    lambda = lambda + s*sd_lambda
  end
  if iter == max_iters
    warn("Hit max iters in interior point method!")
  end
  return x,lambda
end


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md

using JuliaFEM
using JuliaFEM.Testing

@testset "test continuum nonlinear elasticity with surface load" begin

    nodes = Dict{Int64, Node}(
    1 => [0.0, 0.0, 0.0],
    2 => [1.0, 0.0, 0.0],
    3 => [1.0, 1.0, 0.0],
    4 => [0.0, 1.0, 0.0],
    5 => [0.0, 0.0, 1.0],
    6 => [1.0, 0.0, 1.0],
    7 => [1.0, 1.0, 1.0],
    8 => [0.0, 1.0, 1.0])

    element1 = Element(Hex8, [1, 2, 3, 4, 5, 6, 7, 8])
    element2 = Element(Quad4, [5, 6, 7, 8])
    update!([element1, element2], "geometry", nodes)
    update!([element1], "youngs modulus", 900.0)
    update!([element1], "poissons ratio", 0.25)
    update!([element2], "displacement traction force", Vector{Float64}[[0.0, 0.0, -100.0] for i=1:4])

    elasticity_problem = Problem(Elasticity, "solve continuum block", 3)
    elasticity_problem.properties.finite_strain = true
    push!(elasticity_problem, element1)
    push!(elasticity_problem, element2)

    symxy = Element(Quad4, [1, 2, 3, 4])
    symxz = Element(Quad4, [1, 2, 6, 5])
    symyz = Element(Quad4, [1, 4, 8, 5])
    update!([symxy, symxz, symyz], "geometry", nodes)
    symxy["displacement 3"] = 0.0
    symxz["displacement 2"] = 0.0
    symyz["displacement 1"] = 0.0
    boundary_problem = Problem(Dirichlet, "symmetry boundary conditions", 3, "displacement")
    push!(boundary_problem, symxy, symxz, symyz)

    solver = NonlinearSolver("solve 3d block")
    push!(solver, elasticity_problem)
    push!(solver, boundary_problem)
    solver()

    disp = element1("displacement", [1.0, 1.0, 1.0], 0.0)
    info("displacement at tip: $disp")
    # verified using Code Aster.
    # 2015-12-12-continuum-elasticity/vim c3d_grot_gdep_traction_force.comm
    @test isapprox(disp, [3.17431158889468E-02, 3.17431158889468E-02, -1.38591518927826E-01]; rtol=1.0e-4)
end


text:
function test_num_constraints(bridged_mock, F, S, n)
    @test MOI.get(bridged_mock, MOI.NumberOfConstraints{F, S}()) == n
    @test length(MOI.get(bridged_mock, MOI.ListOfConstraintIndices{F, S}())) == n
    @test ((F, S) in MOI.get(bridged_mock, MOI.ListOfConstraints())) == !iszero(n)
end

function warn_incomplete_list_num_constraints(BT, list_num_constraints)
    for (S,) in MOIB.added_constrained_variable_types(BT)
        F = MOIU.variable_function_type(S)
        if !any(c -> c[1] == F && c[2] == S, list_num_constraints)
            error("Bridges of type $BT add constrained variable in $S but their number is not tested.")
        end
    end
    for (F, S) in MOIB.added_constraint_types(BT)
        if !any(c -> c[1] == F && c[2] == S, list_num_constraints)
            error("Bridges of type $BT add $F-in-$S constraints but their number is not tested.")
        end
    end
end

"""
    test_delete_bridge(
        m::MOIB.AbstractBridgeOptimizer, ci::MOI.ConstraintIndex{F, S}, nvars::Int,
        list_num_constraints::Tuple; used_bridges = 1, num_bridged = 1) where {F, S}

Test deletion of the constraint `ci` in model `m`.

* The number of variables added in `m` is `nvars` (ignoring the variables added by the bridges).
* For each `(Fi, Si, n)` in `list_num_constraints`, the number of `Fi`-in-`Si` constraints added
  to `m` is `n` (ignoring the constraints added by the bridges).
  All constraint types added by the bridge bridging `ci` should appear in `list_num_constraints`.
  Otherwise a error is thrown to ensure good test coverage.
* The number of bridges that will be deleted when deleted `ci` is `used_bridges`.
* The number of `F`-in-`S` constraints is `num_bridged`.
"""
function test_delete_bridge(
        m::MOIB.AbstractBridgeOptimizer, ci::MOI.ConstraintIndex{F, S}, nvars::Int,
        list_num_constraints::Tuple; used_bridges = 1, num_bridged = 1) where {F, S}
    warn_incomplete_list_num_constraints(typeof(MOIB.bridge(m, ci)), list_num_constraints)
    function num_bridges()
        return count(bridge -> true, values(MOIB.Constraint.bridges(m)))
    end
    start_num_bridges = num_bridges()
    @test MOI.get(m, MOI.NumberOfVariables()) == nvars
    @test length(MOI.get(m, MOI.ListOfVariableIndices())) == nvars
    test_num_constraints(m, F, S, num_bridged)
    for num_constraints in list_num_constraints
        test_num_constraints(m, num_constraints...)
    end
    @test MOI.is_valid(m, ci)
    MOI.delete(m, ci)
    @test_throws MOI.InvalidIndex(ci) MOI.delete(m, ci)
    @test !MOI.is_valid(m, ci)
    @test num_bridges() == start_num_bridges - used_bridges
    test_num_constraints(m, F, S, num_bridged - 1)
    # As the bridge has been removed, if the constraints it has created where not removed, it wouldn't be there to decrease this counter anymore
    @test MOI.get(m, MOI.NumberOfVariables()) == nvars
    @test length(MOI.get(m, MOI.ListOfVariableIndices())) == nvars
    for num_constraints in list_num_constraints
        test_num_constraints(m, num_constraints...)
    end
end
# Test deletion of variable bridge used for variable `vi`
function test_delete_bridged_variable(
    m::MOIB.AbstractBridgeOptimizer, vi::MOI.VariableIndex, S::Type,
    nvars::Int, list_num_constraints::Tuple; used_bridges = 1, num_bridged = 1, used_constraints = 1)
    warn_incomplete_list_num_constraints(typeof(MOIB.bridge(m, vi)), list_num_constraints)
    function num_bridges()
        return count(bridge -> true, values(MOIB.Variable.bridges(m)))
    end
    start_num_bridges = num_bridges()
    @test MOI.get(m, MOI.NumberOfVariables()) == nvars
    @test length(MOI.get(m, MOI.ListOfVariableIndices())) == nvars
    if S != MOI.Reals
        F = S <: MOI.AbstractScalarSet ? MOI.SingleVariable : MOI.VectorOfVariables
        test_num_constraints(m, F, S, num_bridged)
    end
    for num_constraints in list_num_constraints
        test_num_constraints(m, num_constraints...)
    end
    @test MOI.is_valid(m, vi)
    MOI.delete(m, vi)
    @test_throws MOI.InvalidIndex(vi) MOI.delete(m, vi)
    @test !MOI.is_valid(m, vi)
    @test num_bridges() == start_num_bridges - used_bridges
    if S != MOI.Reals
        test_num_constraints(m, F, S, num_bridged - used_constraints)
    end
    @test MOI.get(m, MOI.NumberOfVariables()) == nvars - 1
    @test length(MOI.get(m, MOI.ListOfVariableIndices())) == nvars - 1
    for num_constraints in list_num_constraints
        test_num_constraints(m, num_constraints...)
    end
end
# Test deletion of variable bridge used for vector of variables `vis`
function test_delete_bridged_variables(
    m::MOIB.AbstractBridgeOptimizer, vis::Vector{MOI.VariableIndex}, S::Type,
    nvars::Int, list_num_constraints::Tuple; used_bridges = 1, num_bridged = 1)
    warn_incomplete_list_num_constraints(typeof(MOIB.bridge(m, vis[1])), list_num_constraints)
    function num_bridges()
        return count(bridge -> true, values(MOIB.Variable.bridges(m)))
    end
    start_num_bridges = num_bridges()
    @test MOI.get(m, MOI.NumberOfVariables()) == nvars
    @test length(MOI.get(m, MOI.ListOfVariableIndices())) == nvars
    if S != MOI.Reals
        F = S <: MOI.AbstractScalarSet ? MOI.SingleVariable : MOI.VectorOfVariables
        test_num_constraints(m, F, S, num_bridged)
    end
    for num_constraints in list_num_constraints
        test_num_constraints(m, num_constraints...)
    end
    @test all(vi -> MOI.is_valid(m, vi), vis)
    MOI.delete(m, vis)
    @test_throws MOI.InvalidIndex(vis[1]) MOI.delete(m, vis)
    @test all(vi -> !MOI.is_valid(m, vi), vis)
    @test num_bridges() == start_num_bridges - used_bridges
    if S != MOI.Reals
        test_num_constraints(m, F, S, num_bridged - 1)
    end
    @test MOI.get(m, MOI.NumberOfVariables()) == nvars - length(vis)
    @test length(MOI.get(m, MOI.ListOfVariableIndices())) == nvars - length(vis)
    for num_constraints in list_num_constraints
        test_num_constraints(m, num_constraints...)
    end
end
function test_delete_objective(
    m::MOIB.AbstractBridgeOptimizer, nvars::Int,
    list_num_constraints::Tuple; used_bridges = 1)
    warn_incomplete_list_num_constraints(typeof(MOIB.Objective.root_bridge(MOIB.Objective.bridges(m))), list_num_constraints)
    function num_bridges()
        return length(MOIB.Objective.bridges(m))
    end
    start_num_bridges = num_bridges()
    @test MOI.get(m, MOI.NumberOfVariables()) == nvars
    @test length(MOI.get(m, MOI.ListOfVariableIndices())) == nvars
    for num_constraints in list_num_constraints
        test_num_constraints(m, num_constraints...)
    end
    MOI.set(m, MOI.ObjectiveSense(), MOI.FEASIBILITY_SENSE)
    @test MOI.get(m, MOI.ObjectiveSense()) == MOI.FEASIBILITY_SENSE
    @test num_bridges() == start_num_bridges - used_bridges
    @test MOI.get(m, MOI.NumberOfVariables()) == nvars
    @test length(MOI.get(m, MOI.ListOfVariableIndices())) == nvars
    for num_constraints in list_num_constraints
        test_num_constraints(m, num_constraints...)
    end
end


text:
using SnoopCompile
using Test

uncompiled(x) = x + 1
if VERSION >= v"1.2.0-DEV.573"
    include_string(Main, """
    @testset "snoopi" begin
        timing_data = @snoopi uncompiled(2)
        @test any(td->td[2].def.name == :uncompiled, timing_data)
        # Ensure older methods can be tested
        a = rand(Float16, 5)
        timing_data = @snoopi sum(a)
        @test any(td->td[2].def.name == :sum, timing_data)
    end
    """)
end

# issue #26
@snoopc "/tmp/anon.log" begin
    map(x->x^2, [1,2,3])
end
data = SnoopCompile.read("/tmp/anon.log")
pc = SnoopCompile.parcel(reverse!(data[2]))
@test length(pc[:Base]) <= 1

#=
# Simple call
let str = "sum"
    keep, pcstring, topmod = SnoopCompile.parse_call("Foo.any($str)")
    @test keep
    @test pcstring == "Tuple{$str}"
    @test topmod == :Main
end

# Operator
let str = "Base.:*, Int, Int"
    keep, pcstring, topmod = SnoopCompile.parse_call("Foo.any($str)")
    @test keep
    @test pcstring == "Tuple{$str}"
    @test topmod == :Base
end

# Function as argument
let str = "typeof(Base.identity), Array{Bool, 1}"
    keep, pcstring, topmod = SnoopCompile.parse_call("Foo.any($str, Vararg{Any, N} where N)")
    @test keep
    @test pcstring == "Tuple{$str, Int}"
    @test topmod == :Base
end

# Anonymous function closure in a new module as argument
let func = (@eval Main module SnoopTestTemp
            func = () -> (y = 2; (x -> x > y))
        end).func
    str = "getfield(SnoopTestTemp, Symbol(\"$(typeof(func()))\")), Array{Float32, 1}"
    keep, pcstring, topmod = SnoopCompile.parse_call("Foo.any($str)")
    @test keep
    @test pcstring == "Tuple{$str}"
    @test topmod == :SnoopTestTemp
end

# Function as a type
let str = "typeof(Base.Sort.sort!), Array{Any, 1}, Base.Sort.MergeSortAlg, Base.Order.By{typeof(Base.string)}"
    keep, pcstring, topmod = SnoopCompile.parse_call("Foo.Bar.sort!($str)")
    @test keep
    @test pcstring == "Tuple{$str}"
    @test topmod == :Base
end
=#

include("colortypes.jl")