text:
# BigBed
# ======
#
# The BigBed format is documented in
#   Kent, W. James, et al. "BigWig and BigBed: # enabling browsing of large
#   distributed datasets." Bioinformatics 26.17 (2010): 2204-2207.
#
# The low level details are documented in a series of tables in the supplement
# of that paper. The immutable types defined below are exactly the data layout
# described in those tables. They are labeled with the table they correspand to.
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/Bio.jl/blob/master/LICENSE.md

const BIGBED_MAGIC = 0x8789F2EB
const BIGWIG_MAGIC = 0x888FFC26
const BIGBED_BTREE_MAGIC = 0x78CA8C91
const BIGBED_RTREE_MAGIC = 0x2468ACE0

const BIGBED_MAX_ZOOM_LEVELS = 10
const BIGBED_CURRENT_VERSION = 4

# See Supplemental Table 5
immutable BigBedHeader
    magic::UInt32
    version::UInt16
    zoom_levels::UInt16
    chromosome_tree_offset::UInt64
    full_data_offset::UInt64
    full_index_offset::UInt64
    field_count::UInt16
    defined_field_count::UInt16
    auto_sql_offset::UInt64
    total_summary_offset::UInt64
    uncompress_buf_size::UInt32
    reserved::UInt64
end

function Base.read(io::IO, ::Type{BigBedHeader})
    return BigBedHeader(
        read(io, UInt32), read(io, UInt16), read(io, UInt16),
        read(io, UInt64), read(io, UInt64), read(io, UInt64),
        read(io, UInt16), read(io, UInt16), read(io, UInt64),
        read(io, UInt64), read(io, UInt32), read(io, UInt64))
end

function Base.write(io::IO, header::BigBedHeader)
    write(io, header.magic)
    write(io, header.version)
    write(io, header.zoom_levels)
    write(io, header.chromosome_tree_offset)
    write(io, header.full_data_offset)
    write(io, header.full_index_offset)
    write(io, header.field_count)
    write(io, header.defined_field_count)
    write(io, header.auto_sql_offset)
    write(io, header.total_summary_offset)
    write(io, header.uncompress_buf_size)
    write(io, header.reserved)
end

# See Supplemental Table 6
immutable BigBedZoomHeader
    reduction_level::UInt32
    reserved::UInt32
    data_offset::UInt64
    index_offset::UInt64
end

function Base.read(io::IO, ::Type{BigBedZoomHeader})
    return BigBedZoomHeader(
        read(io, UInt32), read(io, UInt32), read(io, UInt64), read(io, UInt64))
end

function Base.write(io::IO, header::BigBedZoomHeader)
    write(io, header.reduction_level)
    write(io, header.reserved)
    write(io, header.data_offset)
    write(io, header.index_offset)
end

# Supplemental Table 7
immutable BigBedTotalSummary
    bases_covered::UInt64
    min_val::Float64
    max_val::Float64
    sum_data::Float64
    sum_squares::Float64
end

function Base.read(io::IO, ::Type{BigBedTotalSummary})
    return BigBedTotalSummary(
        read(io, UInt64), read(io, UInt64), read(io, UInt64),
        read(io, UInt64), read(io, UInt64))
end

function Base.write(io::IO, summary::BigBedTotalSummary)
    write(io, summary.bases_covered)
    write(io, summary.min_val)
    write(io, summary.max_val)
    write(io, summary.sum_data)
    write(io, summary.sum_squares)
end

# Supplemental Table 19
immutable BigBedZoomData
    chrom_id::UInt32
    chrom_start::UInt32
    chrom_end::UInt32
    valid_count::UInt32
    min_val::Float32
    max_val::Float32
    sum_data::Float32
    sum_squares::Float32
end

function Base.read(io::IO, ::Type{BigBedZoomData})
    return BigBedZoomData(
        read(io, UInt32), read(io, UInt32), read(io, UInt32), read(io, UInt32),
        read(io, Float32), read(io, Float32), read(Float32), read(Float32))
end

function Base.write(io::IO, data::BigBedZoomData)
    write(io, data.chrom_id)
    write(io, data.chrom_start)
    write(io, data.chrom_end)
    write(io, data.valid_count)
    write(io, data.min_val)
    write(io, data.max_val)
    write(io, data.sum_data)
    write(io, data.sum_squares)
end

# Supplemental Table 8
immutable BigBedBTreeHeader
    magic::UInt32
    block_size::UInt32
    key_size::UInt32
    val_size::UInt32
    item_count::UInt64
    reserved::UInt64
end

function Base.read(io::IO, ::Type{BigBedBTreeHeader})
    return BigBedBTreeHeader(
        read(io, UInt32), read(io, UInt32), read(io, UInt32),
        read(io, UInt32), read(io, UInt64), read(io, UInt64))
end

function Base.write(io::IO, header::BigBedBTreeHeader)
    write(io, header.magic, header.block_size, header.key_size,
          header.val_size, header.item_count, header.reserved)
end

# Supplementary Table 9
immutable BigBedBTreeNode
    isleaf::UInt8
    reserved::UInt8
    count::UInt16
end

function Base.read(io::IO, ::Type{BigBedBTreeNode})
    return BigBedBTreeNode(read(io, UInt8), read(io, UInt8), read(io, UInt16))
end

function Base.write(io::IO, node::BigBedBTreeNode)
    write(io, node.isleaf)
    write(io, node.reserved)
    write(io, node.count)
end

# Supplemental Table 10
type BigBedBTreeLeafNode
    key::Vector{UInt8}
    chrom_id::UInt32
    chrom_size::UInt32
end

function BigBedBTreeLeafNode(keysize::Integer)
    BigBedBTreeLeafNode(Array(UInt8, keysize), 0, 0)
end

function Base.read!(io::IO, node::BigBedBTreeLeafNode)
    nb = readbytes!(io, node.key, length(node.key))
    if nb < length(node.key)
        error("Unexpected end of input.")
    end
    node.chrom_id = read(io, UInt32)
    node.chrom_size = read(io, UInt32)
    return nb + 8
end

# Supplemental Table 11
type BigBedBTreeInternalNode
    key::Vector{UInt8}
    child_offset::UInt64
end

function BigBedBTreeInternalNode(keysize::Integer)
    BigBedBTreeInternalNode(Array(UInt8, keysize), 0)
end

function Base.read!(io::IO, node::BigBedBTreeInternalNode)
    if readbytes!(io, node.key, length(node.key)) < length(node.key)
        error("Unexpected end of input.")
    end
    node.child_offset = read(io, UInt64)
end

# Supplemental Table 13
immutable BigWigSectionHeader
    chrom_id::UInt32
    chrom_start::UInt32
    chrom_end::UInt32
    item_step::UInt32
    item_span::UInt32
    data_type::UInt8
    reserved::UInt8
    item_count::UInt16
end

const BIGWIG_DATATYPE_BEDGRAPH  = UInt8(1)
const BIGWIG_DATATYPE_VARSTEP   = UInt8(2)
const BIGWIG_DATATYPE_FIXEDSTEP = UInt8(3)

function Base.write(io::IO, header::BigWigSectionHeader)
    write(io, header.chrom_id)
    write(io, header.chrom_start)
    write(io, header.chrom_end)
    write(io, header.item_step)
    write(io, header.item_span)
    write(io, header.data_type)
    write(io, header.reserved)
    write(io, header.item_count)
end

immutable BedGraphItem
    chrom_start::UInt32
    chrom_end::UInt32
    val::Float32
end

function Base.write(io::IO, item::BedGraphItem)
    write(io, item.chrom_start)
    write(io, item.chrom_end)
    write(io, item.val)
end

# Supplemental Table 14
immutable BigBedRTreeHeader
    magic::UInt32
    block_size::UInt32
    item_count::UInt64
    start_chrom_ix::UInt32
    start_base::UInt32
    end_chrom_ix::UInt32
    end_base::UInt32
    end_file_offset::UInt64
    items_per_slot::UInt32
    reserved::UInt32
end

function Base.read(io::IO, ::Type{BigBedRTreeHeader})
    return BigBedRTreeHeader(
        read(io, UInt32), read(io, UInt32), read(io, UInt64), read(io, UInt32),
        read(io, UInt32), read(io, UInt32), read(io, UInt32), read(io, UInt64),
        read(io, UInt32), read(io, UInt32))
end

function Base.write(io::IO, header::BigBedRTreeHeader)
    write(io, header.magic)
    write(io, header.block_size)
    write(io, header.item_count)
    write(io, header.start_chrom_ix)
    write(io, header.start_base)
    write(io, header.end_chrom_ix)
    write(io, header.end_base)
    write(io, header.end_file_offset)
    write(io, header.items_per_slot)
    write(io, header.reserved)
end

# Supplemental Table 15
immutable BigBedRTreeNode
    isleaf::UInt8
    reserved::UInt8
    count::UInt16
end

function Base.read(io::IO, ::Type{BigBedRTreeNode})
    return BigBedRTreeNode(read(io, UInt8), read(io, UInt8), read(io, UInt16))
end

function Base.write(io::IO, node::BigBedRTreeNode)
    write(io, node.isleaf)
    write(io, node.reserved)
    write(io, node.count)
end

# Supplemental Table 16
immutable BigBedRTreeLeafNode
    start_chrom_ix::UInt32
    start_base::UInt32
    end_chrom_ix::UInt32
    end_base::UInt32
    data_offset::UInt64
    data_size::UInt64
end

function Base.read(io::IO, ::Type{BigBedRTreeLeafNode})
    return BigBedRTreeLeafNode(
        read(io, UInt32), read(io, UInt32), read(io, UInt32), read(io, UInt32),
        read(io, UInt64), read(io, UInt64))
end

function Base.write(io::IO, node::BigBedRTreeLeafNode)
    write(io, node.start_chrom_ix)
    write(io, node.start_base)
    write(io, node.end_chrom_ix)
    write(io, node.end_base)
    write(io, node.data_offset)
    write(io, node.data_size)
end

# Supplemental Table 17
immutable BigBedRTreeInternalNode
    start_chrom_ix::UInt32
    start_base::UInt32
    end_chrom_ix::UInt32
    end_base::UInt32
    data_offset::UInt64
end

function Base.read(io::IO, ::Type{BigBedRTreeInternalNode})
    return BigBedRTreeInternalNode(
        read(io, UInt32), read(io, UInt32), read(io, UInt32), read(io, UInt32),
        read(io, UInt64))
end

function Base.write(io::IO, node::BigBedRTreeInternalNode)
    write(io, node.start_chrom_ix)
    write(io, node.start_base)
    write(io, node.end_chrom_ix)
    write(io, node.end_base)
    write(io, node.data_offset)
end

immutable BigBed <: FileFormat end
immutable BigWig <: FileFormat end

type BigBedData <: IntervalStream{BEDMetadata}
    stream::BufferedInputStream
    header::BigBedHeader
    zoom_headers::Vector{BigBedZoomHeader}
    autosql::AbstractString
    summary::BigBedTotalSummary
    btree_header::BigBedBTreeHeader
    rtree_header::BigBedRTreeHeader
    data_count::UInt32

    # preallocated space for reading and searching the B-tree
    btree_internal_nodes::Vector{BigBedBTreeInternalNode}
    btree_leaf_nodes::Vector{BigBedBTreeLeafNode}
    key::Vector{UInt8}
    node_keys::Vector{Vector{UInt8}}
    uncompressed_data::Vector{UInt8}
end

type BigBedDataParser <: AbstractParser
    state::Ragel.State

    # intermediate values used during parsing
    chrom_id::UInt32
    red::Float32
    green::Float32
    blue::Float32
    block_size_idx::Int
    block_first_idx::Int
    seq_names::Nullable{Vector{StringField}}
    assumed_seqname::Nullable{StringField}

    function BigBedDataParser(input::BufferedInputStream;
                              seq_names::Nullable{Vector{StringField}}=Nullable{Vector{StringField}}(),
                              assumed_seqname::Nullable{StringField}=Nullable{StringField}())
        cs = _bigbedparser_start
        return new(Ragel.State(cs, input), 0, 0.0, 0.0, 0.0, 1, 1, seq_names, assumed_seqname)
    end
end

include("bigbed-parser.jl")

# Open a BigBed file for reading.  Once opened, entries can be read from the
# file either by iterating over it, or by indexing into it with an interval.
function Base.open(stream::BufferedInputStream, ::Type{BigBed})
    # header
    header = read(stream, BigBedHeader)
    if header.magic != BIGBED_MAGIC
        error("Input is not a valid BigBed file.")
    end

    if header.version < 3
        error("Input is an older unsupported version of BigBed.")
    end

    # zoom headers
    zoom_headers = Array(BigBedZoomHeader, header.zoom_levels)
    for i in 1:header.zoom_levels
        zoom_headers[i] = read(stream, BigBedZoomHeader)
    end

    # autosql
    if header.auto_sql_offset == 0
        autosql = ""
    else
        seek(stream, header.auto_sql_offset)
        autosql_buf = IOBuffer()
        while (c = read(stream, UInt8)) != 0x00
            write(autosql_buf, c)
        end
        # TODO: eventually we should parse this and do something useful with it
        autosql = takebuf_string(autosql_buf)
    end

    # total summary
    seek(stream, header.total_summary_offset)
    summary = read(stream, BigBedTotalSummary)

    # b-tree header
    seek(stream, header.chromosome_tree_offset)
    btree_header = read(stream, BigBedBTreeHeader)
    if btree_header.magic != BIGBED_BTREE_MAGIC
        error("BigBed B-Tree magic number was incorrect. File may be corrupt or malformed.")
    end

    # data_count
    seek(stream, header.full_data_offset)
    data_count = read(stream, UInt32)

    # r-tree header
    seek(stream, header.full_index_offset)
    rtree_header = read(stream, BigBedRTreeHeader)
    if rtree_header.magic != BIGBED_RTREE_MAGIC
        error("BigBed R-Tree magic number was incorrect. File may be corrupt or malformed.")
    end

    return BigBedData(stream, header, zoom_headers, autosql, summary,
                      btree_header, rtree_header, data_count,
                      BigBedBTreeInternalNode[BigBedBTreeInternalNode(btree_header.key_size)
                                              for _ in 1:btree_header.block_size],
                      BigBedBTreeLeafNode[BigBedBTreeLeafNode(btree_header.key_size)
                                          for _ in 1:btree_header.block_size],
                      Array(UInt8, btree_header.key_size),
                      Array(Vector{UInt8}, btree_header.block_size),
                      Array(UInt8, header.uncompress_buf_size))
end

function memisless(a::Vector{UInt8}, b::Vector{UInt8})
    if length(a) != length(b)
        return length(a) < length(b)
    end
    i = 1
    while i ≤ length(a)
        if a[i] != b[i]
            return a[i] < b[i]
        end
        i += 1
    end
    return false
end

# Return all sequence (name, id, size) tuples in a BigBed B-tree.
function first_btree_leaf_position(bb::BigBedData)
    # find the first leaf-node in the b-tree
    offset = bb.header.chromosome_tree_offset + sizeof(BigBedBTreeHeader)
    seek(bb.stream, offset)
    leafpos = 0
    while true
        node = read(bb.stream, BigBedBTreeNode)
        if node.isleaf != 0
            leafpos = offset
            break
        else
            for i in 1:bb.btree_header.block_size
                read!(bb.stream, bb.btree_internal_nodes[i])
                bb.node_keys[i] = bb.btree_internal_nodes[i].key
            end
            offset = bb.btree_internal_nodes[1].child_offset
            seek(bb.stream, offset)
        end
    end

    if leafpos == 0
        error("Malformed BigBed file: no leaf nodes found in index.")
    end

    return leafpos
end

"""
An iterator over all entries in a BigBed file.
"""
type BigBedIteratorState
    seq_names::Vector{StringField}
    data_count::Int
    data_num::Int
    data_offset::UInt
    parser::BigBedDataParser
    parser_isdone::Bool
    next_interval::Interval{BEDMetadata}
end

function Base.start(bb::BigBedData)
    # read sequence names
    leafpos = first_btree_leaf_position(bb)
    seq_names = Array(StringField, bb.btree_header.item_count)
    leafnode = BigBedBTreeLeafNode(bb.btree_header.key_size)
    seek(bb.stream, leafpos)
    i = 1
    while i <= bb.btree_header.item_count
        node = read(bb.stream, BigBedBTreeNode)
        @assert node.isleaf != 0
        for j in 1:node.count
            read!(bb.stream, leafnode)
            p = findfirst(leafnode.key, 0) - 1
            if p == -1
                p = length(leafnode.key)
            end
            seq_names[i] = StringField(leafnode.key[1:p])
            i += 1
        end
    end

    # read the first data block
    seek(bb.stream, bb.header.full_data_offset)
    data_count = read(bb.stream, UInt64)
    zlib_stream = ZlibInflateInputStream(bb.stream, reset_on_end=false)
    unc_block_size = readbytes!(zlib_stream, bb.uncompressed_data,
                               length(bb.uncompressed_data))

    parser = BigBedDataParser(
        BufferedInputStream(bb.uncompressed_data, unc_block_size),
        seq_names=Nullable(seq_names))

    next_interval = BEDInterval()
    parser_isdone = isnull(tryread!(parser, next_interval))

    return BigBedIteratorState(seq_names, data_count, 1,
                               zlib_stream.source.zstream.total_in,
                               parser, parser_isdone, next_interval)
end

function Base.next(bb::BigBedData, state::BigBedIteratorState)
    value = copy(state.next_interval)

    state.data_num += 1
    if state.data_num > state.data_count
        return value, state
    end

    state.parser_isdone = isnull(tryread!(state.parser, state.next_interval))

    if state.parser_isdone
        seek(bb.stream, bb.header.full_data_offset + state.data_offset + sizeof(UInt64))
        zlib_stream = ZlibInflateInputStream(bb.stream, reset_on_end=false)

        unc_block_size = readbytes!(zlib_stream, bb.uncompressed_data,
                                    length(bb.uncompressed_data))
        state.parser = BigBedDataParser(
             BufferedInputStream(bb.uncompressed_data, unc_block_size),
             seq_names=Nullable(state.seq_names))
        state.data_offset += zlib_stream.source.zstream.total_in

        state.parser_isdone = isnull(tryread!(state.parser, state.next_interval))
        @assert !state.parser_isdone
    end

    return value, state
end

function Base.done(bb::BigBedData, state::BigBedIteratorState)
    return state.data_num > state.data_count
end

"""
An iterator over entries in a BigBed file that intersect a given interval.

Constructed by indexing into `BigBedData` with an interval.
"""
type BigBedIntersectIterator
    bb::BigBedData

    query_seqname::StringField
    query_first::Int64
    query_last::Int64
    query_chrom_id::UInt32
    query_chrom_size::UInt32

    # (offset, size) pairs giving the extents of data blocks that may contain
    # (overlapping intervals
    blocks::Vector{Tuple{UInt64, UInt64}}

    # Index of block currently being parsed
    block_num::Int

    parser::Nullable{BigBedDataParser}
    nextinterval::BEDInterval
    done::Bool
end

# Find the given seqname in the BigBed file's index and read the corresponding
# sequence id and length.
function lookup_seqname(bb::BigBedData, seqname::AbstractString)
    seek(bb.stream, bb.header.chromosome_tree_offset + sizeof(BigBedBTreeHeader))

    fill!(bb.key, 0)
    copy!(bb.key, seqname)

    while true
        node = read(bb.stream, BigBedBTreeNode)
        if node.isleaf == 0
            for i in 1:bb.btree_header.block_size
                read!(bb.stream, bb.btree_internal_nodes[i])
                bb.node_keys[i] = bb.btree_internal_nodes[i].key
            end
            i = searchsortedfirst(bb.node_keys, bb.key, Base.Order.Lt(memisless))
            if !(1 <= i <= bb.btree_header.block_size)
                break
            end

            seek(bb.stream, bb.btree_internal_nodes[i].child_offset)
        else
            for i in 1:bb.btree_header.block_size
                read!(bb.stream, bb.btree_leaf_nodes[i])
                bb.node_keys[i] = bb.btree_leaf_nodes[i].key
            end
            i = searchsortedfirst(bb.node_keys, bb.key, Base.Order.Lt(memisless))
            if !(1 <= i <= bb.btree_header.block_size) ||
                bb.btree_leaf_nodes[i].key != bb.key
                break
            end

            return (bb.btree_leaf_nodes[i].chrom_id,
                    bb.btree_leaf_nodes[i].chrom_size)
        end
    end
    error(string("Seqname \"", seqname, "\" is not present in the BigBed file."))
end

function Base.intersect(bb::BigBedData, query::Interval)
    seqname, first, last = query.seqname, query.first, query.last
    chrom_id, chrom_size = lookup_seqname(bb, seqname)

    # stack of nodes yet to visit
    root_position = bb.header.full_index_offset + sizeof(BigBedRTreeHeader)
    node_positions = UInt64[root_position]

    # stack of data blocks that may contain overlapping intervals
    blocks = Tuple{UInt64, UInt64}[]

    seek(bb.stream, bb.header.full_index_offset + sizeof(BigBedRTreeHeader))
    while !isempty(node_positions)
        pos = pop!(node_positions)
        seek(bb.stream, pos)
        node = read(bb.stream, BigBedRTreeNode)
        if node.isleaf == 0
            for i in 1:node.count
                internal_node = read(bb.stream, BigBedRTreeInternalNode)
                if internal_node.start_chrom_ix <= chrom_id <= internal_node.end_chrom_ix &&
                    (chrom_id < internal_node.end_chrom_ix || first <= internal_node.end_base) &&
                    (chrom_id > internal_node.start_chrom_ix || last - 1 >= internal_node.start_base)
                    push!(node_positions, internal_node.data_offset)
                end
            end
        else
            for i in 1:node.count
                leaf_node = read(bb.stream, BigBedRTreeLeafNode)
                if leaf_node.start_chrom_ix <= chrom_id <= leaf_node.end_chrom_ix &&
                    (chrom_id < leaf_node.end_chrom_ix || first <= leaf_node.end_base) &&
                    (chrom_id > leaf_node.start_chrom_ix || last - 1 >= leaf_node.start_base)
                    push!(blocks, (leaf_node.data_offset, leaf_node.data_size))
                end
            end
        end
    end

    return BigBedIntersectIterator(bb, seqname, first, last, chrom_id,
                                   chrom_size, blocks, 0,
                                   Nullable{BigBedDataParser}(),
                                   BEDInterval(), false)
end

function find_next_intersection!(it::BigBedIntersectIterator)
    it.done = true
    while !isnull(it.parser) || it.block_num < length(it.blocks)
        it.done = true
        while !isnull(it.parser)
            parser = get(it.parser)
            # Run the parser until we find an intersection
            it.done = isnull(tryread!(parser, it.nextinterval))
            if it.done
                it.parser = Nullable{BigBedDataParser}()
                break
            end

            if parser.chrom_id == it.query_chrom_id &&
               it.nextinterval.first <= it.query_last &&
               it.nextinterval.last >= it.query_first
                it.done = false
                return
            end
        end

        if it.block_num < length(it.blocks)
            # advance to the next block of interest ant initialize a new parser
            it.block_num += 1
            block_offset, block_size = it.blocks[it.block_num]

            seek(it.bb.stream, block_offset)
            @assert block_size <= length(it.bb.uncompressed_data)
            unc_block_size = readbytes!(ZlibInflateInputStream(it.bb.stream, reset_on_end=false),
                                        it.bb.uncompressed_data,
                                        length(it.bb.uncompressed_data))
            it.parser = BigBedDataParser(
                BufferedInputStream(it.bb.uncompressed_data, unc_block_size),
                assumed_seqname=Nullable(it.query_seqname))
        else
            it.done = true
            return
        end
    end
end

function Base.start(it::BigBedIntersectIterator)
    find_next_intersection!(it)
    return nothing
end

function Base.next(it::BigBedIntersectIterator, ::Void)
    value = copy(it.nextinterval)
    find_next_intersection!(it)
    return value, nothing
end

function Base.done(it::BigBedIntersectIterator, ::Void)
    return it.done
end

# BigBed Output
# -------------

function write_zeros(out::IO, n::Integer)
    for i in 1:n
        write(out, 0x00)
    end
end

immutable BigBedChromInfo
    name::StringField   # chromosome name
    id::UInt32          # unique id for chromosome
    size::UInt32        # size of chromosome
    item_count::UInt64  # number of items
end

immutable BigBedBounds
    offset::UInt64
    chrom_ix::UInt32
    start::UInt32
    stop::UInt32
end

function bigbed_chrom_info(intervals::IntervalCollection, chrom_sizes::Dict)
    info = Array(BigBedChromInfo, length(intervals.trees))
    seqnames = collect(StringField, keys(intervals.trees))
    sort!(seqnames)
    for (i, seqname) in enumerate(seqnames)
        info[i] = BigBedChromInfo(
            seqname,
            i - 1,
            chrom_sizes[seqname],
            length(intervals.trees[seqname]))
    end
    return info
end

function bigbed_btree_count_levels(max_block_size, item_count)
    levels = 1
    while item_count > max_block_size
        item_count = div(item_count + max_block_size - 1, max_block_size)
        levels += 1
    end
    return levels
end

function bigbed_btree_write_index_level(out::IO, block_size,
                                        items::Vector{BigBedChromInfo},
                                        key_size, index_offset, level)
    # See writeIndexLevel in bPlusTree.c

    # calculate number of nodes to write at this level
    slot_size_per = block_size ^ level
    node_size_per = slot_size_per * block_size
    node_count = div(length(items) + node_size_per - 1, node_size_per)

    # calculate sizes and offsets
    block_header_size = 4
    value_size = 8
    bytes_in_index_block = block_header_size + block_size * (key_size + sizeof(UInt64))
    bytes_in_leaf_block = block_header_size + block_size * (key_size + value_size)
    bytes_in_next_level_block = level == 1 ? bytes_in_leaf_block : bytes_in_index_block
    level_size = node_count * bytes_in_index_block
    end_level = index_offset + level_size
    next_child = end_level

    for i in 1:node_size_per:length(items)
        # calculate size of this block
        count_one = min(block_size, div(length(items) - (i-1) + slot_size_per - 1, slot_size_per))

        # write block header
        write(out, BigBedBTreeNode(0, 0, count_one))

        # write used slots
        slots_used = 0
        end_ix = min(length(items), i + node_size_per)
        for j in i:slot_size_per:end_ix
            item = items[j]
            write(out, item.name)
            # pad name
            for k in length(item.name)+1:key_size
                write(out, 0x00)
            end
            write(out, convert(UInt64, next_child))
            next_child += bytes_in_next_level_block
            slots_used += 1
        end
        @assert slots_used == count_one "Incorrect number of items written to index"

        # write out empty slots as all zero
        slot_size = key_size + 8
        for j in count_one:(block_size-1)
            for k in 1:slot_size
                write(out, 0x00)
            end
        end
    end

    return end_level
end

function bigbed_btree_write_leaf_level(out::IO, block_size,
                                       items::Vector{BigBedChromInfo},
                                       key_size, val_size)
    # See writeLeafLevel in bPlusTree.c

    count_one = 0
    count_left = length(items)
    i = 0
    while i < length(items)
        count_one = min(count_left, block_size)

        # write block header
        write(out, BigBedBTreeNode(1, 0, count_one))

        # write out position in genome and in file for each item
        for j in 1:count_one
            @assert i + j <= length(items) "Incorrect item count when writing B-tree leaves"
            item = items[i + j]
            write(out, item.name)
            # pad name
            for k in length(item.name)+1:key_size
                write(out, 0x00)
            end
            write(out, convert(UInt32, item.id))
            write(out, convert(UInt32, item.size))
        end

        # pad out any unused bits of last block with zeroes
        slot_size = key_size + val_size
        for j in count_one:(block_size-1)
            for k in 1:slot_size
                write(out, 0x00)
            end
        end

        count_left -= count_one
        i += count_one
    end
end

function bigbed_write_chrom_tree(out::IO, intervals::IntervalCollection,
                                 block_size,
                                 chrom_sizes::Dict=Dict{AbstractString, Int64}())
    # See bbiWriteChromInfo in bbiWrite.c
    length(intervals.trees)
    for seqname in keys(intervals.trees)
        if !haskey(chrom_sizes, seqname)
            chrom_sizes[seqname] =
                length(intervals.trees[seqname]) > 0 ?
                    intervals.trees[seqname].root.maxend : 0
        end
    end

    max_chrom_name_size = 0
    for seqname in keys(intervals.trees)
        max_chrom_name_size = max(max_chrom_name_size, length(seqname))
    end

    info = bigbed_chrom_info(intervals, chrom_sizes)
    item_count = length(info)
    chrom_block_size = min(block_size, item_count)
    value_size = 8; # sizeof(chrom_id) + sizeof(chrom_size)

    # write header. See bptFileBulkIndex in bPlusTree.c
    header = BigBedBTreeHeader(
        BIGBED_BTREE_MAGIC,
        chrom_block_size,
        max_chrom_name_size,
        value_size,
        item_count,
        0)
    write(out, header)
    index_offset = position(out)

    # write non-leaf nodes.
    levels = bigbed_btree_count_levels(chrom_block_size, item_count)
    for i in (levels-1):-1:1
        end_level_offset = bigbed_btree_write_index_level(out, chrom_block_size,
                                                          info, max_chrom_name_size,
                                                          index_offset, i)
        index_offset = position(out)
        @assert index_offset == end_level_offset "Incorrect file offsets reported."
    end

    # write leaf nodes
    bigbed_btree_write_leaf_level(out, chrom_block_size, info,
                                  max_chrom_name_size, value_size)
    return info
end

function bigbed_count_sections_needed(items::Vector{BigBedChromInfo}, items_per_slot)
    count = 0
    for item in items
        count_one = div(item.item_count + items_per_slot - 1, items_per_slot)
        count += count_one
    end
    return count
end

function bigbed_write_blocks(out::IO, intervals::IntervalCollection,
                             chrom_info::Vector{BigBedChromInfo},
                             items_per_slot, bounds::Vector{BigBedBounds},
                             section_count, compressed::Bool,
                             res_try_count, res_scales, res_sizes)
    # See writeBlocks in bedToBigBed.c

    section_ix = 0
    buf = IOBuffer()
    res_ends = zeros(Int, res_try_count)
    max_block_size = 0

    for info in chrom_info
        tree = intervals.trees[info.name]
        it_state = start(tree)
        fill!(res_ends, 0)

        while !done(tree, it_state)
            section_ix += 1
            item_ix = 1
            block_offset = position(out)
            first_interval = first(intervals.trees[info.name])
            start_pos = -1
            end_pos  = first_interval.last

            while item_ix <= items_per_slot && !done(tree, it_state)
                interval, it_state = next(tree, it_state)

                interval_start = interval.first - 1 # bed is 0-based
                interval_end   = interval.last
                if start_pos < 0
                    start_pos = interval_start
                end
                end_pos = max(end_pos, interval_end)

                write(buf, convert(UInt32, info.id))
                write(buf, convert(UInt32, interval_start))
                write(buf, convert(UInt32, interval_end))
                write_optional_fields(buf, interval, false)
                write(buf, '\0')

                item_ix += 1

                for res_try in 1:res_try_count
                    res_end = res_ends[res_try]
                    if interval_start >= res_end
                        res_sizes[res_try] += 1
                        res_ends[res_try] = res_end = interval_start + res_scales[res_try]
                    end
                    while interval_end > res_end
                        res_sizes[res_try] += 1
                        res_ends[res_try] = res_end = res_end + res_scales[res_try]
                    end
                end
            end

            block = takebuf_array(buf)
            max_block_size = max(max_block_size, length(block))

            if compressed
                writer = ZlibDeflateOutputStream(out)
                write(writer, block)
                flush(writer)
            else
                write(out, block)
            end

            bounds[section_ix] = BigBedBounds(block_offset, info.id,
                                              start_pos, end_pos)
        end
    end

    @assert section_ix == section_count "Incorrect number of sections when writing BigBed data"

    return max_block_size
end

function bigwig_write_blocks{T<:Number}(out::IO, intervals::IntervalCollection{T},
                                        chrom_info::Vector{BigBedChromInfo},
                                        items_per_slot, bounds::Vector{BigBedBounds},
                                        section_count, compressed::Bool,
                                        res_try_count, res_scales, res_sizes)
    section_ix = 0
    res_ends = zeros(Int, res_try_count)
    items = Array(BedGraphItem, items_per_slot)
    max_section_size = 0

    # See writeSections in BedGraphToBigWig.c
    for info in chrom_info
        tree = intervals.trees[info.name]
        it_state = start(tree)
        fill!(res_ends, 0)

        while !done(tree, it_state)
            section_ix += 1
            item_ix = 1

            while item_ix <= items_per_slot && !done(tree, it_state)
                interval, it_state = next(tree, it_state)
                interval_start = interval.first - 1
                interval_end = interval.last
                items[item_ix] = BedGraphItem(
                    interval_start, interval_end, interval.metadata)
                item_ix += 1

                for res_try in 1:res_try_count
                    res_end = res_ends[res_try]
                    if interval_start >= res_end
                        res_sizes[res_try] += 1
                        res_ends[res_try] = res_end = interval_start + res_scales[res_try]
                    end
                    while interval_end > res_end
                        res_sizes[res_try] += 1
                        res_ends[res_try] = res_end = res_end + res_scales[res_try]
                    end
                end
            end

            # write a section
            chrom_id = info.id
            section_start = items[1].chrom_start
            section_end = items[end].chrom_end

            section_offset = position(out)
            header = BigWigSectionHeader(chrom_id, section_start, section_end,
                                         0, 0, BIGWIG_DATATYPE_BEDGRAPH, 0, item_ix - 1)

            if compressed
                writer = ZlibDeflateOutputStream(out)
                write(writer, header)
                for i in 1:item_ix-1
                    write(writer, items[i])
                end
                flush(writer)
            else
                write(out, header)
                for i in 1:item_ix-1
                    write(out, items[i])
                end
            end
            section_end_offset = position(out)
            max_section_size = max(max_section_size, section_end_offset - section_offset)

            bounds[section_ix] = BigBedBounds(section_offset, info.id,
                                              section_start, section_end)
        end
    end

    @assert section_ix == section_count "Incorrect number of sections when writing BigWig data"

    return max_section_size
end

type BigBedRTree
    next::Nullable{BigBedRTree}
    children::Nullable{BigBedRTree}
    parent::Nullable{BigBedRTree}
    start_chrom_ix::UInt32
    start_base::UInt32
    end_chrom_ix::UInt32
    end_base::UInt32
    start_file_offset::UInt64
    end_file_offset::UInt64

    function BigBedRTree(next, children, parent, start_chrom_ix, start_base,
                         end_chrom_ix, end_base, start_file_offset, end_file_offset)
        return new(next, children, parent, start_chrom_ix, end_chrom_ix,
                   start_base, end_base, start_file_offset, end_file_offset)
    end

    function BigBedRTree()
        return new()
    end
end

function Base.copy(tree::BigBedRTree)
    return BigBedRTree(tree.next, tree.children, tree.parent,
                       tree.start_chrom_ix, tree.start_base,
                       tree.end_chrom_ix, tree.end_base,
                       tree.start_file_offset, tree.end_file_offset)
end

function Base.length(list::Nullable{BigBedRTree})
    count = 0
    while !isnull(list)
        count += 1
        list = get(list).next
    end
    return count
end

function bigbed_calc_rtree_level_sizes!(level_sizes, tree::Nullable{BigBedRTree},
                                        level, max_level)
    # See calcLevelSizes in cirTree.c

    el = tree
    while !isnull(el)
        level_sizes[level] += 1
        if level <= max_level
            bigbed_calc_rtree_level_sizes!(level_sizes, get(el).children,
                                           level + 1, max_level)
        end
        el = get(el).next
    end
end

function bigbed_rtree_write_index_level(out::IO, block_size, child_node_size,
                                        tree::BigBedRTree, offset_of_first_child,
                                        cur_level, dest_level)
    offset = offset_of_first_child
    if cur_level == dest_level
        # we've reached the right level, write out a node header
        count_one = length(tree.children)
        write(out, BigBedRTreeNode(0, 0, count_one))

        # write out elements of this node
        nel = tree.children
        while !isnull(nel)
            el = get(nel)
            write(out, BigBedRTreeInternalNode(el.start_chrom_ix, el.start_base,
                                               el.end_chrom_ix, el.end_base, offset))
            offset += child_node_size
            nel = el.next
        end

        # write out zeroes for empty slots in node
        for i in count_one:(block_size-1)
            for j in 1:sizeof(BigBedRTreeInternalNode)
                write(out, 0x00)
            end
        end
    else
        # otherwise recurse on children
        nel = tree.children
        while !isnull(nel)
            el = get(nel)
            offset = bigbed_rtree_write_index_level(out, block_size,
                                                    child_node_size, el,
                                                    cur_level + 1,
                                                    dest_level,
                                                    offset)
            nel = el.next
        end
    end
    return offset
end

function bigbed_rtree_write_leaf_level(out::IO, items_per_slot, leaf_node_size,
                                       tree::BigBedRTree, cur_level, leaf_level)
    if cur_level == leaf_level
        # we've reached the right level, write out node header
        count_one = length(tree.children)
        write(out, BigBedRTreeNode(1, 0, count_one))

        nel = tree.children
        while !isnull(nel)
            el = get(nel)
            write(out, BigBedRTreeLeafNode(el.start_chrom_ix, el.start_base,
                                           el.end_chrom_ix, el.end_base,
                                           el.start_file_offset,
                                           el.end_file_offset - el.start_file_offset))
            nel = el.next
        end

        # write out zeroes for empty slots in node
        for i in count_one:(items_per_slot-1)
            for j in 1:sizeof(BigBedRTreeInternalNode)
                write(out, 0x00)
            end
        end
    else
        # otherwise recurse on children
        nel = tree.children
        while !isnull(nel)
            el = get(nel)
            bigbed_rtree_write_leaf_level(out, items_per_slot, leaf_node_size,
                                          el, cur_level + 1, leaf_level)
            nel = el.next
        end
    end
end

function Base.write(out::IO, tree::BigBedRTree, block_size, level_count)
    # See writeTreeToOpenFile in cirTree.c

    level_sizes = zeros(Int, level_count)
    bigbed_calc_rtree_level_sizes!(level_sizes, Nullable{BigBedRTree}(tree), 1, level_count)

    level_offsets = Array(UInt64, level_count)

    index_node_size = sizeof(BigBedRTreeNode) + block_size * sizeof(BigBedRTreeInternalNode)
    leaf_node_size = sizeof(BigBedRTreeNode) + block_size * sizeof(BigBedRTreeLeafNode)

    offset = position(out)
    for i in 1:level_count
        level_offsets[i] = offset
        offset += level_sizes[i] * index_node_size
    end

    final_level = level_count - 2
    for i in 1:final_level
        child_node_size = i == final_level ? leaf_node_size : index_node_size
        bigbed_rtree_write_index_level(out, block_size, child_node_size, tree,
                                       level_offsets[i+1], 1, i)
        @assert position(out) == level_offsets[i+1]
    end

    leaf_level = level_count - 1
    bigbed_rtree_write_leaf_level(out, block_size, leaf_node_size, tree,
                                  1, leaf_level)
end

function Base.reverse!(tree::Nullable{BigBedRTree})
    if isnull(tree)
        return tree
    end

    nodes = BigBedRTree[]
    node = tree
    while !isnull(node)
        push!(nodes, get(node))
        node = get(node).next
    end

    for i in length(nodes):-1:2
        nodes[i].next = nodes[i-1]
    end
    nodes[1].next = Nullable{BigBedRTree}()
    return Nullable{BigBedRTree}(nodes[end])
end

function bigbed_build_rtree(bounds::Vector{BigBedBounds}, block_size,
                            items_per_slot, end_file_offset)
    # See rTreeFromChromRangeArray in cirTree.c

    list  = Nullable{BigBedRTree}()
    next_offset = isempty(bounds) ? 0 : bounds[1].offset
    for i in 1:items_per_slot:length(bounds)
        final_iteration = false
        one_size = length(bounds) - i + 1
        if one_size > items_per_slot
            one_size = items_per_slot
        else
            final_iteration = true
        end

        bound = bounds[i]
        el = BigBedRTree()
        el.next = Nullable{BigBedRTree}()
        el.children = Nullable{BigBedRTree}()
        el.parent = Nullable{BigBedRTree}()
        el.start_chrom_ix = el.end_chrom_ix = bound.chrom_ix
        el.start_base = bound.start
        el.end_base = bound.stop
        el.start_file_offset = next_offset

        # Figure out end of element from offset of next element (or file size
        # for final element.)
        if final_iteration
            next_offset = end_file_offset
        else
            next_offset = bounds[i + one_size].offset
        end
        el.end_file_offset = next_offset

        # expand area spanned to include all items in block
        for j in 2:one_size
            bound = bounds[i + j - 1]
            if bound.chrom_ix < el.start_chrom_ix
                el.start_chrom_ix = bound.chrom_ix
                el.start_base = bound.start
            elseif bound.chrom_ix == el.start_chrom_ix
                el.start_base = min(el.start_base, bound.start)
            elseif bound.chrom_ix > el.end_chrom_ix
                el.end_chrom_ix = bound.chrom_ix
                el.end_base = bound.stop
            elseif bound.chrom_ix == el.end_chrom_ix
                el.end_base = max(el.end_base, bound.stop)
            end
        end

        el.next = list
        list = Nullable{BigBedRTree}(el)
    end
    list = reverse!(list)

    # Now iterate through making more and more condensed versions until have
    # just one
    tree = list
    level_count = 1
    while (!isnull(tree) && !isnull(get(tree).next)) || level_count < 2
        list = Nullable{BigBedRTree}()
        slots_used = block_size
        parent = Nullable{BigBedRTree}()
        nel = tree
        while !isnull(nel)
            el = get(nel)
            next_el = el.next
            if slots_used >= block_size
                slots_used = 1
                parent = Nullable{BigBedRTree}(copy(el))
                get(parent).children = el
                el.parent = parent
                el.next = Nullable{BigBedRTree}()

                get(parent).next = list
                list = parent
            else
                slots_used += 1
                el.next = get(parent).children
                get(parent).children = el
                el.parent = parent
                if el.start_chrom_ix < get(parent).start_chrom_ix
                    get(parent).start_chrom_ix = el.start_chrom_ix
                    get(parent).start_base = el.start_base
                elseif el.start_chrom_ix == get(parent).start_chrom_ix
                    get(parent).start_base = min(el.start_base, get(parent).start_base)
                end

                if el.end_chrom_ix > get(parent).end_chrom_ix
                    get(parent).end_chrom_ix = el.end_chrom_ix
                    get(parent).end_base = el.end_base
                elseif el.end_chrom_ix == get(parent).end_chrom_ix
                    get(parent).end_base  = max(el.end_base, get(parent).end_base)
                end
            end
            nel = next_el
        end

        list = reverse!(list)
        el = list
        while !isnull(el)
            get(el).children = reverse!(get(el).children)
            el = get(el).next
        end
        tree = list
        level_count += 1
    end

    @assert !isnull(tree)

    (get(tree), level_count)
end

function bigbed_write_index(out::IO, bounds::Vector{BigBedBounds},
                            block_size, items_per_slot, end_file_offset)
    # See cirTreeFileBulkIndexToOpenFile in cirTree.c

    rtree, level_count = bigbed_build_rtree(bounds, block_size, items_per_slot,
                                            end_file_offset)
    write(out, BigBedRTreeHeader(BIGBED_RTREE_MAGIC, block_size, length(bounds),
                                 rtree.start_chrom_ix, rtree.start_base,
                                 rtree.end_chrom_ix, rtree.end_base,
                                 end_file_offset, items_per_slot, 0))
    write(out, rtree, block_size, level_count)
end

# Write out sum to file, keeping track of minimal info on it in, and also adding
# it to second level summary
function output_one_summary_further_reduce(out::IO, sum::BigBedZoomData,
                                           twice_reduced_list::Vector{BigBedZoomData},
                                           double_reduction_size,
                                           bounds_array::Vector{BigBedBounds},
                                           bounds_array_idx,
                                           chrom_size)
    # See bbiOutputOneSummaryFurtherReduce in bbiWrite.c

    @assert bounds_array_idx <= length(bounds_array)

    bounds_array[bounds_array_idx] = BigBedBounds(position(out), sum.chrom_id,
                                                  sum.chrom_start, sum.chrom_end)
    write(out, sum)

    # Fold summary info into twice_reduced_list
    if isempty(twice_reduced_list) || twice_reduced_list[end].chrom_id != sum.chrom_id ||
       twice_reduced_list[end].chrom_start + double_reduction_size < sum.chrom_end
        push!(twice_reduced_list, sum)
    else
        twice_reduced_list[end] = BigBedZoomData(
            twice_reduced_list[end].chrom_id,
            twice_reduced_list[end].chrom_start,
            sum.chrom_end,
            twice_reduced_list[end].valid_count + sum.valid_count,
            min(sum.min_val, twice_reduced_list[end].min_val),
            max(sum.max_val, twice_reduced_list[end].max_val),
            twice_reduced_list[end].sum_data + sum.sum_data,
            twice_reduced_list[end].sum_squares + sum.sum_squares)
    end

    return bounds_array_idx + 1
end

# Write out data reduced by factor of initialReduction.  Also calculate and keep
# in memory next reduction level.  This is more work than some ways, but it
# keeps us from having to keep the first reduction entirely in memory.
function bigbed_write_reduced_once_return_reduced_twice(
                out::IO, intervals::IntervalCollection,
                interval_tree_transform::Function,
                chrom_info::Vector{BigBedChromInfo},
                initial_reduction, initial_reduced_count,
                zoom_increment, block_size, items_per_slot, compressed)

    # See writeReducedOnceReturnReducedTwice in bedToBigBed.c

    double_reduction_size = initial_reduction * zoom_increment
    bounds_array = Array(BigBedBounds, initial_reduced_count)
    bounds_array_idx = 1
    twice_reduced_list = BigBedZoomData[]

    ret_data_start = position(out)
    write(out, convert(UInt32, initial_reduced_count))

    # This gets a little complicated I'm afraid.  The strategy is to:
    #  1) Build up a range tree that represents coverage depth on that chromosome
    #     This also has the nice side effect of getting rid of overlaps.
    #  2) Stream through the range tree, outputting the initial summary level and
    #     further reducing.
    #

    valid_count = UInt64(0)
    min_val = 0.0
    max_val = 0.0
    sum_data = 0.0
    sum_squares = 0.0
    first_time = true

    sum = Nullable{BigBedZoomData}()

    for info in chrom_info
        sum = Nullable{BigBedZoomData}()
        for range in interval_tree_transform(intervals.trees[info.name])
            val = convert(Float64, range.metadata)
            chrom_start = range.first - 1
            chrom_end = range.last
            size = chrom_end - chrom_start
            @assert size >= 0 "Bad size: $chrom_start, $chrom_end"
            if first_time
                valid_count = size
                min_val = max_val = val
                sum_data = val * size
                sum_squares = val * val * size
                first_time = false
            else
                valid_count += size
                min_val = min(min_val, val)
                max_val = max(max_val, val)
                sum_data += val * size
                sum_squares += val * val * size
            end

            # If start past existing block then output it
            if !isnull(sum) && get(sum).chrom_end <= chrom_start
                bounds_array_idx = output_one_summary_further_reduce(
                        out, get(sum), twice_reduced_list,
                        double_reduction_size, bounds_array,
                        bounds_array_idx, info.size)
                sum = Nullable{BigBedZoomData}()
            end

            # If don't have a summary we're working on now, make one
            if isnull(sum)
                sum = Nullable{BigBedZoomData}(
                    BigBedZoomData(info.id, chrom_start,
                                   min(info.size, chrom_start + initial_reduction),
                                   0, val, val, 0.0, 0.0))
            end

            # Deal with case where might have to split an item between multiple
            # summaries. This loop handles all but the final affected summary in
            # that case
            s = get(sum)
            while chrom_end > s.chrom_end
                # Fold in bits that ovelap with existing summary and output
                overlap = min(chrom_end, s.chrom_end) - max(chrom_start, s.chrom_start)
                @assert overlap > 0
                @assert overlap <= size

                s = BigBedZoomData(
                        s.chrom_id, s.chrom_start, s.chrom_end,
                        s.valid_count + overlap,
                        min(s.min_val, val), max(s.max_val, val),
                        s.sum_data + val * overlap,
                        s.sum_squares + val * val * overlap)

                bounds_array_idx = output_one_summary_further_reduce(
                        out, s, twice_reduced_list, double_reduction_size,
                        bounds_array, bounds_array_idx, info.size)
                size -= overlap

                # Move summary to next part
                chrom_start = s.chrom_end
                #chrom_end = min(info.size, chrom_start + initial_reduction)
                s = BigBedZoomData(
                        s.chrom_id, chrom_start,
                        min(info.size, chrom_start + initial_reduction),
                        0, val, val, 0.0, 0.0)
            end

            # Add to summary
            s = BigBedZoomData(
                    s.chrom_id, s.chrom_start, s.chrom_end,
                    s.valid_count + size,
                    min(s.min_val, val), max(s.max_val, val),
                    s.sum_data + val * size,
                    s.sum_data + val * val * size)
            sum = Nullable(s)
        end

        if !isnull(sum)
            bounds_array_idx = output_one_summary_further_reduce(
                    out, get(sum), twice_reduced_list, double_reduction_size,
                    bounds_array, bounds_array_idx, info.size)
        end
    end

    # Write out 1st zoom index
    index_offset = position(out)
    @assert bounds_array_idx == length(bounds_array) + 1
    bigbed_write_index(out, bounds_array, block_size, items_per_slot,
                       index_offset)
    return twice_reduced_list, BigBedTotalSummary(valid_count, min_val, max_val,
                                                  sum_data, sum_squares)
end

function bigbed_write_summary_and_index(
        out::IO, summary_list::Vector{BigBedZoomData}, block_size,
        items_per_slot, compressed)
    # See: bbiWriteSummaryAndIndex in bbiWrite.c
    if compressed
        bigbed_write_summary_and_index_comp(out, summary_list, block_size, items_per_slot)
    else
        bigbed_write_summary_and_index_unc(out, summary_list, block_size, items_per_slot)
    end
end

function bigbed_write_summary_and_index_unc(
        out::IO, summary_list::Vector{BigBedZoomData}, block_size,
        items_per_slot)
    # See: bbiWriteSummaryAndIndexUnc in bbiWrite.c
    count = length(summary_list)
    write(out, convert(UInt32, count))
    bounds_array = Array(BigBedBounds, count)
    for (i, summary) in enumerate(summary_list)
        bounds_array[i] = BigBedBounds(
            position(out), summary.chrom_id, summary.chrom_start,
            summary.chrom_end)
        write(out, summary)
    end
    index_offset = position(out)
    bigbed_write_index(out, bounds_array, block_size, items_per_slot, false)
    return index_offset
end

function bigbed_write_summary_and_index_comp(
        out::IO, summary_list::Vector{BigBedZoomData}, block_size,
        items_per_slot)
    # See: bbiWriteSummaryAndIndexComp in bbiWrite.c
    count = length(summary_list)
    write(out, convert(UInt32, count))
    bounds_array = Array(BigBedBounds, count)

    items_left = count
    sum_ix = 1
    i = 1
    while items_left > 0
        items_in_slot = min(items_per_slot, items_left)
        writer = ZlibDeflateOutputStream(out)
        file_pos = position(out)
        for _ in 1:items_in_slot
            summary = summary_list[sum_ix]
            bounds_array[sum_ix] = BigBedBounds(file_pos, summary.chrom_id,
                                                summary.chrom_start, summary.chrom_end)
            write(writer, summary_list[sum_ix])
            sum_ix += 1
            if sum_ix > count
                break
            end
        end
        items_left -= items_in_slot
        flush(writer)
    end
    index_offset = position(out)
    bigbed_write_index(out, bounds_array, block_size, items_per_slot, false)
    return index_offset
end

function bigbed_summary_simple_reduce(summary_list::Vector{BigBedZoomData}, reduction)
    new_summary_list = BigBedZoomData[]
    for summary in summary_list
        if isempty(new_summary_list) ||
           new_summary_list[end].chrom_id != summary.chrom_id ||
           summary.chrom_end > new_summary_list[end].chrom_start + reduction
            push!(new_summary_list, summary)
        else
            @assert new_summary_list[end].chrom_end < summary.chrom_end
            new_summary_list[end] = BigBedZoomData(
                new_summary_list[end].chrom_id,
                new_summary_list[end].chrom_start,
                summary.chrom_end,
                new_summary_list[end].valid_count + summary.valid_count,
                min(new_summary_list[end].min_val, summary.min_val),
                max(new_summary_list[end].max_val, summary.max_val),
                new_summary_list[end].sum_data + summary.sum_data,
                new_summary_list[end].sum_squares + summary.sum_squares)
       end
    end
    return new_summary_list
end

function bed_field_count(intervals::IntervalCollection{BEDMetadata})
    num_fields = 0
    for interval in intervals
        interval_num_fields = 3 + interval.metadata.used_fields
        if num_fields == 0 || num_fields > interval_num_fields
            num_fields = interval_num_fields
        end
    end
    return num_fields
end

function Base.write(out::IO, ::Type{BigBed}, intervals::IntervalCollection{BEDMetadata};
                    block_size::Int=256, items_per_slot::Int=512,
                    compressed::Bool=true)
    write_bigbed_bigwig(out, BigBed, intervals, block_size, items_per_slot, compressed)
end

function Base.write{T<:Number}(out::IO, ::Type{BigWig}, intervals::IntervalCollection{T};
                               block_size::Int=256, items_per_slot::Int=512,
                               compressed::Bool=true)
    write_bigbed_bigwig(out, BigWig, intervals, block_size, items_per_slot, compressed)
end

# Unified write function for bigwig and bigbed
function write_bigbed_bigwig(out::IO, fmt::Union{Type{BigBed}, Type{BigWig}},
                             intervals::IntervalCollection,
                             block_size::Int=256, items_per_slot::Int=512,
                             compressed::Bool=true)

    # See the function bbFileCreate in bedToBigBed.c in kent to see the only
    # other implementation of this I'm aware of. Also bedGraphToBigWig in
    # bedGraphToBigWig.c

    if !applicable(seek, out, 1)
        error("BigBed can only be written to seekable output streams.")
    end

    field_count = fmt == BigBed ? bed_field_count(intervals) : 0

    # write dummy headers, come back and fill them later
    write_zeros(out, sizeof(BigBedHeader))
    write_zeros(out, BIGBED_MAX_ZOOM_LEVELS * sizeof(BigBedZoomHeader))

    # TODO: optional autoSql specification
    as_offset = 0

    total_summary_offset = position(out)
    total_summary = BigBedTotalSummary(0, 0.0, 0.0, 0.0, 0.0)
    write(out, total_summary)

    chrom_tree_offset = position(out)
    chrom_info = bigbed_write_chrom_tree(out, intervals, block_size)

    # set up to keep track of possible initial reduction levels
    total_bases = 0
    for interval in intervals
        total_bases += interval.last - interval.first + 1
    end
    ave_span = total_bases / length(intervals)

    res_try_count = 10
    res_increment = 4
    res_scales = Array(Int, res_try_count)
    res_sizes = Array(Int, res_try_count)
    min_zoom = 10
    res = max(round(Int, ave_span), min_zoom)
    for res_try in 1:res_try_count
        res_sizes[res_try] = 0
        res_scales[res_try] = res
        if res > 1000000000
            res_try_count = res_try
            break
        end
        res *= res_increment
    end

    # write out primary full resolution data in sections, collect stats to
    # use for reductions
    data_offset = position(out)
    # TODO: The BigBed/BigWig specification differs from Kent's impementation
    # here. The spec says this count is a 32bit int, but the implmentation uses
    # a 64bit int. We follow the implementation for practical reasons.
    write(out, convert(UInt64, length(intervals)))
    block_count = bigbed_count_sections_needed(chrom_info, items_per_slot)
    bounds = Array(BigBedBounds, block_count)

    if fmt == BigBed
        max_block_size = bigbed_write_blocks(out, intervals, chrom_info, items_per_slot,
                                             bounds, block_count, compressed,
                                             res_try_count, res_scales, res_sizes)
    else
        max_block_size = bigwig_write_blocks(out, intervals, chrom_info, items_per_slot,
                                             bounds, block_count, compressed,
                                             res_try_count, res_scales, res_sizes)
    end

    # write out primary data index
    index_offset = position(out)
    bigbed_write_index(out, bounds, block_size, 1, index_offset)

    # declare arrays and vars that track the zoom levels we actually output
    zoom_levels = 0
    zoom_amounts = Array(UInt32, BIGBED_MAX_ZOOM_LEVELS)
    zoom_data_offsets = Array(UInt64, BIGBED_MAX_ZOOM_LEVELS)
    zoom_index_offsets = Array(UInt64, BIGBED_MAX_ZOOM_LEVELS)
    zoom_levels = 0

    if ave_span > 0
        data_size = index_offset - data_offset
        max_reduced_size = div(data_size, 2)
        initial_reduction = 0
        initial_reduced_count = 0

        for res_try in 1:res_try_count
            reduced_size = res_sizes[res_try] * sizeof(BigBedZoomData)
            if compressed
                reduced_size = div(reduced_size, 2) # estimate!
            end

            if reduced_size <= max_reduced_size
                initial_reduction = res_scales[res_try]
                initial_reduced_count = res_sizes[res_try]
                break
            end
        end

        if initial_reduction > 0
            zoom_increment = 4
            rezoomed_list, total_summary = bigbed_write_reduced_once_return_reduced_twice(
                    out, intervals,
                    fmt == BigBed ? coverage : identity,
                    chrom_info, initial_reduction,
                    initial_reduced_count, zoom_increment, block_size,
                    items_per_slot, compressed)

            zoom_amounts[1] = initial_reduction
            zoom_levels = 1

            zoom_count = initial_reduced_count
            reduction = initial_reduction * zoom_increment
            while zoom_levels < BIGBED_MAX_ZOOM_LEVELS
                rezoom_count = length(rezoomed_list)
                if rezoom_count >= zoom_count
                    break
                end

                zoom_count = rezoom_count
                zoom_data_offsets[zoom_levels] = position(out)
                zoom_index_offsets[zoom_levels] =
                        bigbed_write_summary_and_index(out, rezoomed_list, block_size,
                                                       items_per_slot, compressed)
                zoom_amounts[zoom_levels] = reduction
                zoom_levels += 1
                reduction *= zoom_increment

                rezoomed_list = bigbed_summary_simple_reduce(rezoomed_list, reduction)
            end
        end
    end

    uncompress_buf_size = 0
    if compressed
        max_zoom_uncomp_size = items_per_slot * sizeof(BigBedZoomData)
        uncompress_buf_size = max(max_block_size, max_zoom_uncomp_size)
    end

    # go back and rewrite header
    seek(out, 0)
    write(out,
        BigBedHeader(fmt == BigBed ? BIGBED_MAGIC : BIGWIG_MAGIC,
                     BIGBED_CURRENT_VERSION, zoom_levels,
                     chrom_tree_offset, data_offset, index_offset,
                     field_count, min(field_count, 12), as_offset,
                     total_summary_offset, uncompress_buf_size, 0))

    # write summary headers with data
    for i in 1:zoom_levels
        write(out, BigBedZoomHeader(zoom_amounts[i], 0, zoom_data_offsets[i],
                                    zoom_index_offsets[i]))
    end

    # write rest of summary header with no data
    for i in zoom_levels:(BIGBED_MAX_ZOOM_LEVELS-1)
        write(out, BigBedZoomHeader(0, 0, 0, 0))
    end

    # write total summary
    seek(out, total_summary_offset)
    write(out, total_summary)

    # write end signature
    seekend(out)
    write(out, convert(UInt32, fmt == BigBed ? BIGBED_MAGIC : BIGWIG_MAGIC))
end


text:
import Knet
using FileIO, ImageCore, ImageMagick, ImageTransformations

# main
function imagenet_preprocess(img::Matrix{<:RGB}; normalize=nothing, resolution=224, format="whcn", atype=Knet.atype)
    if normalize === nothing; normalize = identity; end
    minsize = resolution * 8 ÷ 7                          # 224 => 256
    img = imresize(img, ratio=minsize/minimum(size(img))) # min(h,w)=256
    hcenter,vcenter = size(img) .>> 1                     # h÷2, w÷2
    half = resolution ÷ 2
    img = img[hcenter-half+1:hcenter+half, vcenter-half+1:vcenter+half] # hw=224,224
    img = Float32.(channelview(img))                                    # chw=3,224,224
    img = normalize(img)
    img = reshape(img, (1, size(img)...)) # nchw=1,3,224,224
    if format != "nchw"
        perm = findfirst.(collect(format), "nchw")
        img = permutedims(img, (perm...,))
    end
    return atype(img)
end

# fallback for arbitrary array types
imagenet_preprocess(x; atype=Knet.atype, o...) = atype(x)

# grayscale images
function imagenet_preprocess(img::Matrix{<:Gray}; o...)
    imagenet_preprocess(RGB.(img); o...)
end

# files and urls
function imagenet_preprocess(file::String; o...)
    if occursin(r"^https?://", file)
        mktemp() do fn,io
            imagenet_preprocess(load(download(url,fn)); o...)
        end
    else
        imagenet_preprocess(load(file); o...)
    end
end

ext: jl


text:
function _DiscontinuousFESpace(
  vector_type::Type,
  trian::Triangulation,
  cell_fe::CellFE)

  cell_shapefuns, cell_dof_basis = compute_cell_space(cell_fe,trian)

  cell_dof_ids, nfree = compute_discontinuous_cell_dofs(cell_fe.cell_ctype,cell_fe.ctype_num_dofs)

  ndirichlet = 0
  dirichlet_dof_tag = Int8[]
  dirichlet_cells = Int32[]
  cell_is_dirichlet = Fill(false,num_cells(trian))
  ntags = 0

  UnconstrainedFESpace(
    vector_type,
    nfree,
    ndirichlet,
    cell_dof_ids,
    cell_shapefuns,
    cell_dof_basis,
    cell_is_dirichlet,
    dirichlet_dof_tag,
    dirichlet_cells,
    ntags)
end

function compute_discontinuous_cell_dofs(cell_to_ctype,ctype_to_nldofs)

  ncells = length(cell_to_ctype)
  ptrs = zeros(Int32,ncells+1)
  for (cell, ctype) in enumerate(cell_to_ctype)
    nldofs = ctype_to_nldofs[ctype]
    ptrs[cell+1] = nldofs
  end

  length_to_ptrs!(ptrs)

  ndata = ptrs[end]-1
  data = collect(Int32,1:ndata)

  (Table(data,ptrs), ndata)

end


text:
# Autogenerated wrapper script for libwebp_jll for x86_64-w64-mingw32
export cwebp, dwebp, gif2webp, img2webp, libwebp, libwebpdecoder, libwebpdemux, libwebpmux, webpinfo, webpmux

using Giflib_jll
using JpegTurbo_jll
using libpng_jll
using Libtiff_jll
using Libglvnd_jll
JLLWrappers.@generate_wrapper_header("libwebp")
JLLWrappers.@declare_executable_product(cwebp)
JLLWrappers.@declare_executable_product(dwebp)
JLLWrappers.@declare_executable_product(gif2webp)
JLLWrappers.@declare_executable_product(img2webp)
JLLWrappers.@declare_library_product(libwebp, "libwebp-7.dll")
JLLWrappers.@declare_library_product(libwebpdecoder, "libwebpdecoder-3.dll")
JLLWrappers.@declare_library_product(libwebpdemux, "libwebpdemux-2.dll")
JLLWrappers.@declare_library_product(libwebpmux, "libwebpmux-3.dll")
JLLWrappers.@declare_executable_product(webpinfo)
JLLWrappers.@declare_executable_product(webpmux)
function __init__()
    JLLWrappers.@generate_init_header(Giflib_jll, JpegTurbo_jll, libpng_jll, Libtiff_jll, Libglvnd_jll)
    JLLWrappers.@init_executable_product(
        cwebp,
        "bin\\cwebp.exe",
    )

    JLLWrappers.@init_executable_product(
        dwebp,
        "bin\\dwebp.exe",
    )

    JLLWrappers.@init_executable_product(
        gif2webp,
        "bin\\gif2webp.exe",
    )

    JLLWrappers.@init_executable_product(
        img2webp,
        "bin\\img2webp.exe",
    )

    JLLWrappers.@init_library_product(
        libwebp,
        "bin\\libwebp-7.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libwebpdecoder,
        "bin\\libwebpdecoder-3.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libwebpdemux,
        "bin\\libwebpdemux-2.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libwebpmux,
        "bin\\libwebpmux-3.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        webpinfo,
        "bin\\webpinfo.exe",
    )

    JLLWrappers.@init_executable_product(
        webpmux,
        "bin\\webpmux.exe",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
# AVIO

import Base: read, read!, show, close, eof, isopen, seekstart

export read, read!, pump, openvideo, opencamera, playvideo, viewcam, play

type StreamInfo
    stream_index0::Int             # zero-based
    stream::AVStream
    codec_ctx::AVCodecContext
end

abstract StreamContext

if isdefined(Main, :Color)
    typealias EightBitTypes Union(UInt8, Main.FixedPointNumbers.Ufixed8, Main.Color.RGB{Main.FixedPointNumbers.Ufixed8})
elseif isdefined(Main, :FixedPointNumbers)
    typealias EightBitTypes Union(UInt8, Main.FixedPointNumbers.Ufixed8)
else
    typealias EightBitTypes UInt8
end

# An audio-visual input stream/file
type AVInput{I}
    io::I
    apFormatContext::Vector{Ptr{AVFormatContext}}
    apAVIOContext::Vector{Ptr{AVIOContext}}
    avio_ctx_buffer_size::Uint
    aPacket::Vector{AVPacket}           # Reusable packet

    unknown_info::Vector{StreamInfo}
    video_info::Vector{StreamInfo}
    audio_info::Vector{StreamInfo}
    data_info::Vector{StreamInfo}
    subtitle_info::Vector{StreamInfo}
    attachment_info::Vector{StreamInfo}

    listening::Set{Int}
    stream_contexts::Array{StreamContext}

    isopen::Bool
end

function show(io::IO, avin::AVInput)
    println(io, "AVInput(", avin.io, ", ...), with")
    (len = length(avin.video_info))      > 0 && println(io, "  $len video stream(s)")
    (len = length(avin.audio_info))      > 0 && println(io, "  $len audio stream(s)")
    (len = length(avin.data_info))       > 0 && println(io, "  $len data stream(s)")
    (len = length(avin.subtitle_info))   > 0 && println(io, "  $len subtitle stream(s)")
    (len = length(avin.attachment_info)) > 0 && println(io, "  $len attachment stream(s)")
    (len = length(avin.unknown_info))    > 0 && println(io, "  $len unknown stream(s)")
end

type VideoTranscodeContext
    transcode_context::Ptr{SwsContext}

    transcode_interp::Cint
    target_pix_fmt::Cint
    target_bits_per_pixel::Int
    width::Int
    height::Int

    target_buf::Array{UInt8,3}  # TODO: this needs to be more general
    aTargetVideoFrame::Vector{AVFrame}
    apTargetDataBuffers::Vector{Ptr{UInt8}}
    apTargetLineSizes::Vector{Cint}
end

const TRANSCODE=true
const NO_TRANSCODE=false

type VideoReader{transcode} <: StreamContext
    avin::AVInput
    stream_info::StreamInfo

    stream_index0::Int
    pVideoCodecContext::Ptr{AVCodecContext}
    pVideoCodec::Ptr{AVCodec}
    aVideoFrame::Vector{AVFrame}   # Reusable frame
    aFrameFinished::Vector{Int32}

    format::Cint
    width::Cint
    height::Cint
    framerate::Rational
    aspect_ratio::Rational

    frame_queue::Vector{Vector{UInt8}}
    transcodeContext::VideoTranscodeContext
end

show(io::IO, vr::VideoReader) = print(io, "VideoReader(...)")

# type AudioContext <: StreamContext
#     stream_index0::Int             # zero-based
#     stream::AVStream
#     codec_ctx::AVCodecContext

#     sample_format::Int
#     sample_rate::Int
#     #sample_bits::Int
#     channels::Int
# end

# type SubtitleContext <: StreamContext
#     stream_index0::Int             # zero-based
#     stream::AVStream
#     codec_ctx::AVCodecContext
# end

# type DataContext <: StreamContext
#     stream_index0::Int             # zero-based
#     stream::AVStream
#     codec_ctx::AVCodecContext
# end

# type AttachmentContext <: StreamContext
#     stream_index0::Int             # zero-based
#     stream::AVStream
#     codec_ctx::AVCodecContext
# end

# type UnknownContext <: StreamContext
#     stream_index0::Int             # zero-based
#     stream::AVStream
#     codec_ctx::AVCodecContext
# end

# Pump input for data
function pump(c::AVInput)
    pFormatContext = c.apFormatContext[1]

    while true
        !c.isopen && break

        Base.sigatomic_begin()
        av_read_frame(pFormatContext, pointer(c.aPacket)) < 0 && break
        Base.sigatomic_end()

        packet = c.aPacket[1]
        stream_index = packet.stream_index

        # If we're not listening to this stream, skip it
        if stream_index in c.listening
            # Decode the packet, and check if the frame is complete
            frameFinished = decode_packet(c.stream_contexts[stream_index+1], c.aPacket)
            av_free_packet(pointer(c.aPacket))

            # If the frame is complete, we're done
            frameFinished && return stream_index
        else
            av_free_packet(pointer(c.aPacket))
        end
    end

    return -1
end

pump(r::StreamContext) = pump(r.avin)

function _read_packet(pavin::Ptr{AVInput}, pbuf::Ptr{UInt8}, buf_size::Cint)
    avin = unsafe_pointer_to_objref(pavin)
    out = pointer_to_array(pbuf, (buf_size,))
    convert(Cint, readbytes!(avin.io, out))
end

const read_packet = cfunction(_read_packet, Cint, (Ptr{AVInput}, Ptr{UInt8}, Cint))

function open_avinput(avin::AVInput, io::IO, input_format=C_NULL)

    !isreadable(io) && error("IO not readable")

    # These allow control over how much of the stream to consume when
    # determining the stream type
    # TODO: Change these defaults if necessary, or allow user to set
    #av_opt_set(avin.apFormatContext[1], "probesize", "100000000", 0)
    #av_opt_set(avin.apFormatContext[1], "analyzeduration", "1000000", 0)

    # Allocate the io buffer used by AVIOContext
    # Must be done with av_malloc, because it could be reallocated
    if (avio_ctx_buffer = av_malloc(avin.avio_ctx_buffer_size)) == C_NULL
        error("unable to allocate AVIOContext buffer (out of memory)")
    end

    # Allocate AVIOContext
    if (avin.apAVIOContext[1] = avio_alloc_context(avio_ctx_buffer, avin.avio_ctx_buffer_size,
                                                   0, pointer_from_objref(avin),
                                                   read_packet, C_NULL, C_NULL)) == C_NULL
        abuffer = [avio_ctx_buffer]
        av_freep(abuffer)
        error("Unable to allocate AVIOContext (out of memory)")
    end

    # pFormatContext->pb = pAVIOContext
    av_setfield(avin.apFormatContext[1], :pb, avin.apAVIOContext[1])

    # "Open" the input
    if avformat_open_input(avin.apFormatContext, C_NULL, input_format, C_NULL) != 0
        error("Unable to open input")
    end

    nothing
end

function open_avinput(avin::AVInput, source::String, input_format=C_NULL)
    if avformat_open_input(avin.apFormatContext,
                           source,
                           input_format,
                           C_NULL)    != 0
        error("Could not open file $source")
    end

    nothing
end

function AVInput{T<:Union(IO, String)}(source::T, input_format=C_NULL; avio_ctx_buffer_size=65536)

    # Register all codecs and formats
    av_register_all()
    av_log_set_level(AVUtil.AV_LOG_ERROR)

    aPacket = [AVPacket()]
    apFormatContext = Ptr{AVFormatContext}[avformat_alloc_context()]
    apAVIOContext = Ptr{AVIOContext}[C_NULL]

    # Allocate this object (needed to pass into AVIOContext in open_avinput)
    avin = AVInput{T}(source, apFormatContext, apAVIOContext, avio_ctx_buffer_size,
                      aPacket, [StreamInfo[] for i=1:6]..., Set(Int[]), StreamContext[], false)

    # Make sure we deallocate everything on exit
    finalizer(avin, close)

    # Set up the format context and open the input, based on the type of source
    open_avinput(avin, source, input_format)
    avin.isopen = true

    # Get the stream information
    if avformat_find_stream_info(avin.apFormatContext[1], C_NULL) < 0
        error("Unable to find stream information")
    end

    # Load streams, codec_contexts
    formatContext = unsafe_load(avin.apFormatContext[1]);

    for i = 1:formatContext.nb_streams
        pStream = unsafe_load(formatContext.streams,i)
        stream = unsafe_load(pStream)
        codec_ctx = unsafe_load(stream.codec)
        codec_type = codec_ctx.codec_type

        stream_info = StreamInfo(i-1, stream, codec_ctx)

        if codec_type == AVMEDIA_TYPE_VIDEO
            push!(avin.video_info, stream_info)
        elseif codec_type == AVMEDIA_TYPE_AUDIO
            push!(avin.audio_info, stream_info)
        elseif codec_type == AVMEDIA_TYPE_DATA
            push!(avin.data_info, stream_info)
        elseif codec_type == AVMEDIA_TYPE_SUBTITLE
            push!(avin.subtitle_info, stream_info)
        elseif codec_type == AVMEDIA_TYPE_ATTACHMENT
            push!(avin.attachment_info, stream_info)
        elseif codec_type == AVMEDIA_TYPE_UNKNOWN
            push!(avin.unknown_info, stream_info)
        end
    end

    resize!(avin.stream_contexts, formatContext.nb_streams)

    avin
end

function VideoReader(avin::AVInput, video_stream=1;
                     transcode::Bool=true,
                     transcode_interpolation=SWS_BILINEAR,
                     target_format=PIX_FMT_RGB24)

    1 <= video_stream <= length(avin.video_info) || error("video stream $video_stream not found")

    stream_info = avin.video_info[video_stream]

    # Get basic stream info
    pVideoCodecContext = stream_info.stream.codec
    codecContext = stream_info.codec_ctx

    width, height = codecContext.width, codecContext.height
    pix_fmt = codecContext.pix_fmt
    pix_fmt < 0 && error("Unknown pixel format")

    framerate = codecContext.time_base.den // codecContext.time_base.num

    aspect_num = codecContext.sample_aspect_ratio.num
    aspect_den = codecContext.sample_aspect_ratio.den
    aspect_ratio = (aspect_den == 0) ? 0 // 1 : aspect_num // aspect_den

    # Find the decoder for the video stream
    pVideoCodec = avcodec_find_decoder(codecContext.codec_id)
    pVideoCodec == C_NULL && error("Unsupported Video Codec")

    # Open the decoder
    avcodec_open2(pVideoCodecContext, pVideoCodec, C_NULL) < 0 && error("Could not open codec")

    aVideoFrame = [AVFrame()]
    aFrameFinished = Int32[0]

    aTargetVideoFrame = [AVFrame()]

    # Set up transcoding
    # TODO: this should be optional
    pFmtDesc = av_pix_fmt_desc_get(target_format)
    bits_per_pixel = av_get_bits_per_pixel(pFmtDesc)

    if transcode && bits_per_pixel % 8 != 0
        error("Unsupported format (bits_per_pixel = $bits_per_pixel)")
    end

    N = @compat Int64(bits_per_pixel >> 3)
    target_buf = Array(UInt8, bits_per_pixel>>3, width, height)

    sws_context = sws_getContext(width, height, pix_fmt,
                                 width, height, target_format,
                                 transcode_interpolation, C_NULL, C_NULL, C_NULL)

    avpicture_fill(pointer(aTargetVideoFrame), target_buf, target_format, width, height)
    apTargetDataBuffers = reinterpret(Ptr{UInt8}, [aTargetVideoFrame[1].data])
    apTargetLineSizes   = reinterpret(Cint,       [aTargetVideoFrame[1].linesize])

    transcodeContext = VideoTranscodeContext(sws_context,
                                             transcode_interpolation,
                                             target_format,
                                             bits_per_pixel,
                                             width,
                                             height,
                                             target_buf,
                                             aTargetVideoFrame,
                                             apTargetDataBuffers,
                                             apTargetLineSizes)

    vr = VideoReader{transcode}(avin,
                                stream_info,

                                stream_info.stream_index0,
                                pVideoCodecContext,
                                pVideoCodec,
                                aVideoFrame,
                                aFrameFinished,

                                pix_fmt,
                                width,
                                height,
                                framerate,
                                aspect_ratio,

                                Vector{UInt8}[],
                                transcodeContext)

    idx0 = stream_info.stream_index0
    push!(avin.listening, idx0)
    avin.stream_contexts[idx0+1] = vr

    vr
end

VideoReader{T<:Union(IO, String)}(s::T, args...; kwargs...) = VideoReader(AVInput(s), args...; kwargs... )

function decode_packet(r::VideoReader, aPacket)
    # Do we already have a complete frame that hasn't been consumed?
    if have_decoded_frame(r)
        apSourceDataBuffers = reinterpret(Ptr{UInt8}, [r.aVideoFrame[1].data])
        sz = avpicture_get_size(r.format, r.width, r.height)
        push!(r.frame_queue, copy(pointer_to_array(apSourceDataBuffers[1], sz)))
        reset_frame_flag!(r)
    end

    avcodec_decode_video2(r.pVideoCodecContext,
                          r.aVideoFrame,
                          r.aFrameFinished,
                          aPacket)

    return have_decoded_frame(r)
end

# Converts a grabbed frame to the correct format (RGB by default)
# TODO: type stability?
function retrieve(r::VideoReader{TRANSCODE}) # true=transcode

    while !have_frame(r)
        idx = pump(r.avin)
        idx == r.stream_index0 && break
        idx == -1 && throw(EOFError())
    end

    t = r.transcodeContext

    if t.target_bits_per_pixel % 8 != 0
        error("Unsupported video retrieval format.  Please allocate the buffer yourself and call retrieve!()")
    end

    if t.target_bits_per_pixel == 8
        buf = Array(UInt8, r.width, r.height)
    else
        buf = Array(UInt8, t.target_bits_per_pixel>>3, r.width, r.height)
    end

    retrieve!(r, buf)
end

function retrieve(r::VideoReader{NO_TRANSCODE}) # false=don't transcode
    while !have_frame(r)
        idx = pump(r.avin)
        idx == r.stream_index0 && break
        idx == -1 && throw(EOFError())
    end

    # TODO: set actual dimensions ?
    buf_sz = avpicture_get_size(r.format, r.width, r.height)
    buf = Array(UInt8, buf_sz)

    retrieve!(r, buf)
end

# Converts a grabbed frame to the correct format (RGB by default)
function retrieve!{T<:EightBitTypes}(r::VideoReader{TRANSCODE}, buf::Array{T})
    while !have_frame(r)
        idx = pump(r.avin)
        idx == r.stream_index0 && break
        idx == -1 && throw(EOFError())
    end

    if !bufsize_check(r, buf)
        error("Buffer is the wrong size")
    end

    t = r.transcodeContext

    # Set up target buffer
    if pointer(buf) != t.apTargetDataBuffers[1]
        avpicture_fill(pointer(t.aTargetVideoFrame),
                       buf,
                       t.target_pix_fmt,
                       r.width, r.height)
        t.apTargetDataBuffers = reinterpret(Ptr{UInt8}, [t.aTargetVideoFrame[1].data])
        t.apTargetLineSizes   = reinterpret(Cint,       [t.aTargetVideoFrame[1].linesize])

    end

    apSourceDataBuffers = isempty(r.frame_queue) ? reinterpret(Ptr{UInt8}, [r.aVideoFrame[1].data]) :
                                                   reinterpret(Ptr{UInt8}, [shift!(r.frame_queue)])
    apSourceLineSizes   = reinterpret(Cint, [r.aVideoFrame[1].linesize])

    Base.sigatomic_begin()
    sws_scale(t.transcode_context,
              apSourceDataBuffers,
              apSourceLineSizes,
              zero(Int32),
              r.height,
              t.apTargetDataBuffers,
              t.apTargetLineSizes)
    Base.sigatomic_end()

    reset_frame_flag!(r)

    return buf
end

function retrieve!{T<:EightBitTypes}(r::VideoReader{NO_TRANSCODE}, buf::Array{T})
    while !have_frame(r)
        idx = pump(r.avin)
        idx == r.stream_index0 && break
        idx == -1 && throw(EOFError())
    end

    if !bufsize_check(r, buf)
        error("Buffer is the wrong size")
    end
end

# Utility functions

# Not exported
open(filename::String) = AVInput(filename)
openvideo(args...; kwargs...) = VideoReader(args...; kwargs...)

read(r::VideoReader) = retrieve(r)
read!{T<:EightBitTypes}(r::VideoReader, buf::AbstractArray{T}) = retrieve!(r, buf)

isopen{I<:IO}(avin::AVInput{I}) = isopen(avin.io)
isopen(avin::AVInput) = avin.isopen
isopen(r::VideoReader) = isopen(r.avin)

bufsize_check{T<:EightBitTypes}(r::VideoReader{NO_TRANSCODE}, buf::Array{T}) = (length(buf)*sizeof(T) == avpicture_get_size(r.format, r.width, r.height))
bufsize_check{T<:EightBitTypes}(r::VideoReader{TRANSCODE}, buf::Array{T}) = bufsize_check(r.transcodeContext, buf)
bufsize_check{T<:EightBitTypes}(t::VideoTranscodeContext, buf::Array{T}) = (length(buf)*sizeof(T) == avpicture_get_size(t.target_pix_fmt, t.width, t.height))

have_decoded_frame(r) = r.aFrameFinished[1] > 0  # TODO: make sure the last frame was made available
have_frame(r::StreamContext) = !isempty(r.frame_queue) || have_decoded_frame(r)
have_frame(avin::AVInput) = any(Bool[have_frame(avin.stream_contexts[i+1]) for i in avin.listening])

reset_frame_flag!(r) = (r.aFrameFinished[1] = 0)

function seekstart(s::VideoReader, video_stream=1)
    !isopen(s) && throw(ErrorException("Video input stream is not open!"))

    pCodecContext = s.pVideoCodecContext # AVCodecContext

    seekstart(s.avin, video_stream)
    avcodec_flush_buffers(pCodecContext)

    return s
end

function seekstart{T<:String}(avin::AVInput{T}, video_stream = 1)
    # AVFormatContext
    fc = avin.apFormatContext[1]

    # Get stream information
    stream_info = avin.video_info[video_stream]
    seek_stream_index = stream_info.stream_index0
    stream = stream_info.stream
    first_dts = stream.first_dts

    # Seek
    ret = avformat_seek_file(fc, seek_stream_index, first_dts, first_dts, first_dts, AVSEEK_FLAG_BACKWARD)

    ret < 0 && throw(ErrorException("Could not seek to start of stream"))

    return avin
end

## This doesn't work...
#seekstart{T<:IO}(avin::AVInput{T}, video_stream = 1) = seekstart(avin.io)
seekstart{T<:IO}(avin::AVInput{T}, video_stream = 1) = throw(ErrorException("Sorry, Seeking is not supported from IO streams"))

function eof(avin::AVInput)
    !isopen(avin) && return true
    have_frame(avin) && return false
    got_frame = (pump(avin) != -1)
    return !got_frame
end

function eof{I<:IO}(avin::AVInput{I})
    !isopen(avin) && return true
    have_frame(avin) && return false
    return eof(avin.io)
end

eof(r::VideoReader) = eof(r.avin)

close(r::VideoReader) = close(r.avin)
_close(r::VideoReader) = avcodec_close(r.pVideoCodecContext)

# Free AVIOContext object when done
function close(avin::AVInput)
    Base.sigatomic_begin()
    isopen = avin.isopen
    avin.isopen = false
    Base.sigatomic_end()

    !isopen && return

    for i in avin.listening
        _close(avin.stream_contexts[i+1])
    end

    # Fix for segmentation fault issue #44
    empty!(avin.listening)

    Base.sigatomic_begin()
    if avin.apFormatContext[1] != C_NULL
        avformat_close_input(avin.apFormatContext)
    end
    Base.sigatomic_end()

    Base.sigatomic_begin()
    if avin.apAVIOContext[1] != C_NULL
        p = av_pointer_to_field(avin.apAVIOContext[1], :buffer)
        p != C_NULL && av_freep(p)
        av_freep(avin.apAVIOContext)
    end
    Base.sigatomic_end()
end

### Camera Functions

if have_avdevice()
    import AVDevice
    AVDevice.avdevice_register_all()

    function get_camera_devices(ffmpeg, idev, idev_name)
        CAMERA_DEVICES = UTF8String[]

        read_vid_devs = false
        out,err = readall_stdout_stderr(`$ffmpeg -list_devices true -f $idev -i $idev_name`)
        buf = length(out) > 0 ? out : err
        for line in eachline(IOBuffer(buf))
            if contains(line, "video devices")
                read_vid_devs = true
                continue
            elseif contains(line, "audio devices") || contains(line, "exit") || contains(line, "error")
                read_vid_devs = false
                continue
            end

            if read_vid_devs
                m = match(r"""\[.*"(.*)".?""", line)
                if m != nothing
                    push!(CAMERA_DEVICES, m.captures[1])
                end

                # Alternative format (TODO: could be combined with the regex above)
                m = match(r"""\[.*\] \[[0-9]\] (.*)""", line)
                if m != nothing
                    push!(CAMERA_DEVICES, m.captures[1])
                end
            end
        end

        return CAMERA_DEVICES
    end

    @windows_only begin
        ffmpeg = joinpath(Pkg.dir("VideoIO"), "deps", "ffmpeg-2.2.3-win$WORD_SIZE-shared", "bin", "ffmpeg.exe")

        DEFAULT_CAMERA_FORMAT = AVFormat.av_find_input_format("dshow")
        CAMERA_DEVICES = get_camera_devices(ffmpeg, "dshow", "dummy")
        DEFAULT_CAMERA_DEVICE = length(CAMERA_DEVICES) > 0 ? CAMERA_DEVICES[1] : "0"

    end

    @linux_only begin
        import Glob
        DEFAULT_CAMERA_FORMAT = AVFormat.av_find_input_format("video4linux2")
        CAMERA_DEVICES = Glob.glob("video*", "/dev")
        DEFAULT_CAMERA_DEVICE = length(CAMERA_DEVICES) > 0 ? CAMERA_DEVICES[1] : ""
    end

    @osx_only begin
        ffmpeg = joinpath(INSTALL_ROOT, "bin", "ffmpeg")

        DEFAULT_CAMERA_FORMAT = AVFormat.av_find_input_format("avfoundation")
        global CAMERA_DEVICES = UTF8String[]
        try
            CAMERA_DEVICES = get_camera_devices(ffmpeg, "avfoundation", "\"\"")
        catch
            try
                CAMERA_DEVICES = get_camera_devices(ffmpeg, "qtkit", "\"\"")
            end
        end

        DEFAULT_CAMERA_DEVICE = length(CAMERA_DEVICES) > 0 ? CAMERA_DEVICES[1] : "FaceTime"
        #DEFAULT_CAMERA_DEVICE = "Integrated"
    end

    function opencamera(device=DEFAULT_CAMERA_DEVICE, format=DEFAULT_CAMERA_FORMAT, args...; kwargs...)
        camera = AVInput(device, format)
        VideoReader(camera, args...; kwargs...)
    end
end

try
    if isa(Main.Images, Module)
        # Define read and retrieve for Images
        global retrieve, retrieve!, read!, read
        for r in [:read, :retrieve]
            r! = symbol("$(r)!")

            @eval begin
                # read!, retrieve!
                $r!(c::VideoReader, img::Main.Images.Image) = ($r!(c, Main.Images.data(img)); img)

                # read, retrieve
                function $r(c::VideoReader, ::Type{Main.Images.Image}, colorspace="RGB", colordim=1, spatialorder=["x","y"])
                    img = Main.Images.colorim($r(c::VideoReader))
                end
            end
        end
    end
end

try
    if isa(Main.ImageView, Module)
        # Define read and retrieve for Images
        global playvideo, viewcam, play

        function play(f, flip=false)
            img = read(f, Main.Images.Image)
            canvas, _ = Main.ImageView.view(img, flipx=flip, interactive=false)
            buf = Main.Images.data(img)

            while !eof(f)
                read!(f, buf)
                Main.ImageView.view(canvas, img, flipx=flip, interactive=false)
                sleep(1/f.framerate)
            end
        end

        function playvideo(video)
            f = VideoIO.openvideo(video)
            play(f)
        end

        if have_avdevice()
            function viewcam(device=DEFAULT_CAMERA_DEVICE, format=DEFAULT_CAMERA_FORMAT)
                camera = opencamera(device, format)
                play(camera, true)
            end
        else
            function viewcam()
                error("libavdevice not present")
            end
        end
    end
catch
    global playvideo, viewcam, play
    no_imageview() = error("Please load ImageView before VideoIO to enable play(...), playvideo(...) and viewcam()")
    play() = no_imageview()
    playvideo() = no_imageview()
    viewcam() = no_imageview()
end


text:
results =  []
which = 0
mypush!(a) = let
    global which
    which += 1
    push!(results, (which, a))
end

# Shin, Lee
# mypush!([100.46, 98.71, 98.74])
# mypush!([90.75, 92.90, 95.06])
mypush!([73.45, 81.47, 87.36])
mypush!([108.3, 103.04, 101.39])
mypush!([101.85, 100.79, 100.43])
mypush!([102.31, 100.84, 100.43])

# Sangtarash et al
mypush!([104.5, 101.0, 100.3])
# mypush!([0.94, 0.98, 1.01] .* 100)
mypush!([104.8, 100.5, 99.7])
mypush!([94.3, 97.2, 98.8])
mypush!([103.6, 100.2, 99.6])
mypush!([104.7, 100.5, 99.7])
mypush!([116.8, 102.8, 100.8])
# mypush!([77.1, 83.0, 94.0])
mypush!([94.3, 97.2, 98.9])
mypush!([106.6, 104.4, 103.4])
# mypush!([99.5, 98.6, 99.3])
mypush!([93.7, 97.4, 99.0])
mypush!([94.3, 97.3, 98.9])
# mypush!([])
mypush!([99.4, 99.7, 99.9])

# Rezaiee
mypush!([58.08, 90.94, 98.84])
mypush!([52.86, 88.31, 98.09])
mypush!([44.63, 84.28, 96.68])
mypush!([62.04, 92.25, 99.04])
# mypush!([# ])
mypush!([94.37, 97.25, 98.90])
mypush!([104.7, 101.1, 99.6])

# Briassoulis
mypush!([0.964, 0.984, 0.998] .*100)
mypush!([120.1, 104.6, 101.0])
# mypush!([94.0, 97.0, 100.0]
# mypush!([])
mypush!([108.3, 101.5, 100.0])
# mypush!([94.0, 98.0, 101.0])
mypush!([95.0, 97.4, 98.7])

# Argyris
# mypush!([69.7, 90.2, 100.1])

# Cen, LI
mypush!([103.74, 100.29, 99.61])
mypush!([100.14, 99.56, 99.51])
mypush!([103.56, 100.19, 99.66])
mypush!([116.8, 102.8, 100.8])
# mypush!([92.84, 96.09, 99.08])
# mypush!([77.12, 83.07, 94.31])
mypush!([119.12, 104.2, 100.63])
# mypush!([])
# mypush!([94.2, 100.8, 100.5])

# Ko et al
mypush!([95.50, 98.51, 99.32])
mypush!([94.32, 97.26, 98.86])
mypush!([104.8, 100.5, 99.73])
mypush!([108.3, 101.5, 100])
mypush!([104.8, 100.5, 99.6])
mypush!([104.7, 100.5, 99.74])

# Jun, et al
# mypush!([73.12, 87.43, 95.93])
# mypush!([66.77, 85.58, 95.40])
# mypush!([72.78, 90.15, 100.43])
mypush!([96.10, 99.31, 99.83])
mypush!([89.22, 97.62, 99.50])
mypush!([94.31, 97.29, 98.88])
# mypush!([94.01, 98.02, 100.99])
# mypush!([104.99, 100.78])
mypush!([121.89, 105.39, 101.69])

# Nguyen-Van et al
mypush!([119.12, 104.2, 100.63])
mypush!([119.25, 104.22, 100.66])
mypush!([120.1, 104.6, 101.0])
mypush!([97.6, 98.6, 99.3])
mypush!([104.7, 100.5, 99.7])
mypush!([121.9, 105.4, 101.7])

# Mostafa
mypush!([104.9, 101.3, 99.6])
mypush!([104.7, 101.1, 99.6])
# mypush!([99.7, ])
mypush!([98.6, 99.3, 99.6])
mypush!([96.0, 98.4, 99.9])
mypush!([102.9, 100.1, 99.2])

# Coupling Effect
mypush!([.326802, .316916, .313979] ./ (0.3024) .* 100)
mypush!([0.329038, 0.317741, 0.314225] ./ (0.3024) .* 100)
# # mypush!([0.095964, ] ./ (0.3024) .* 100)
# mypush!([.220731, .274365, .307401] ./ (0.3024) .* 100)
# mypush!([.284311, .293325, .305336] ./ (0.3024) .* 100)
mypush!([.316522, .303931, .301618] ./ (0.3024) .* 100)

mypush!(-0.23851E+00, -2.994369E-01, -3.021210E-01] ./ (-0.3024) .* 100)

# The extrapolated value with Richardson formula is 0.3022447

differences = []
for r in results
    @show r[1]
    @show q1, q2, q3 = r[2]
    @show qtrue = (q2^2 - q1 * q3) / (2*q2 - q1 - q3)
    push!(differences, (r[1], qtrue - 100))
end

differences = sort(differences, lt = (x, y) -> x[2] < y[2])
@show differences
# @show [(i, v) in  for (i, v) in enumerate(differences)]

using PGFPlotsX

objects = []

@pgf p = PGFPlotsX.Plot(
{
color = "black",
line_width  = 0.7,
style = "solid",
mark = "*"
},
Coordinates([(i, v[2]) for (i, v) in enumerate(differences)])
)
push!(objects, p)
# mypush!(objects, LegendEntry("$(r[1])"))

@pgf ax = Axis(
{
xlabel = "Study [ND]",
ylabel = "Departure from 0.3024 [\\%]",
            # xmin = range[1],
            # xmax = range[2],
xmode = "linear",
ymode = "linear",
yminorgrids = "true",
grid = "both",
legend_style = {
at = Coordinate(0.5, 1.05),
anchor = "south",
legend_columns = -1
},
},
objects...
)

display(ax)
pgfsave("Convergence-metaanalysis.pdf", ax)


text:
module Artificial

using CUDA
using LinearAlgebra
using MatrixOptInterface
using SparseArrays

import ..PhaseOne
import Simplex
import Simplex: nrows, ncols

const MatOI = MatrixOptInterface
const MOI = MatOI.MOI

"""
    This creates a reformulation to begin with all artificial basis.
    This is the simplest way to construct the initial basis.

    The canonical form

    minimize c x
    subject to
        A x == b
        xl <= x <= xu

    is reformulated to

    minimize a1 + a2
    subject to
        A1 x - a1 == b - A1 xN
        A2 x + a2 == b - A2 xN
        xl <= x <= xu
        0 <= a1 <= Inf
        0 <= a2 <= Inf,

    where
        if xl[j] > -Inf
            xN[j] = xl[j]
        elseif xu[j] < Inf
            xN[j] = xu[j]
        else
            xN[j] = 0.0
        end
"""
function reformulate(lp::MatOI.LPSolverForm{T, AT, VT}; basis::Array{Int,1} = Int[]) where {T, AT, VT}

    artif_idx = Int[]
    if length(basis) > 0
        # Add the artifical variables that are identified as basis
        for j in 1:length(basis)
            # if basis indicates an artificial variable,
            if basis[j] > ncols(lp)
                if basis[j] - ncols(lp) > nrows(lp)
                    @error "Invalid basis provided for phase one"
                end
                push!(artif_idx, basis[j])
                basis[j] = ncols(lp) + length(artif_idx)
            end
        end
    else
        # if no basis is provided, return (m+1,...,m+n)
        artif_idx = collect(1:nrows(lp)) .+ ncols(lp)
    end

    nartif = length(artif_idx)
    lp_nrows = nrows(lp)
    lp_ncols = ncols(lp) + nartif

    # objective coefficient
    c = append!(zeros(ncols(lp)), ones(nartif))
    c = VT(c)

    # column bounds
    # TODO: build thing directly on GPU
    xl = append!(deepcopy(Array(lp.v_lb)), zeros(nartif))
    xu = append!(deepcopy(Array(lp.v_ub)), fill(Inf, nartif))
    xl = VT(xl)
    xu = VT(xu)

    x = VT(undef, lp_ncols)

    for j = 1:lp_ncols
        if xl[j] > -Inf && abs(xl[j]) <= abs(xu[j])
            x[j] = xl[j]
        elseif xu[j] < Inf && abs(xl[j]) > abs(xu[j])
            x[j] = xu[j]
        else
            x[j] = 0.0
        end
    end

    # rhs
    b = VT(undef, nrows(lp))
    copyto!(b, lp.b .- lp.A * x[1:ncols(lp)])

    # create the submatrix for slack
    I = Array{Int64}(undef, nartif)
    J = collect(1:nartif)
    V = VT(undef, nartif)
    for i = 1:nartif
        I[i] = artif_idx[i] - ncols(lp)
        if b[I[i]] < 0
            V[i] = -1
            x[ncols(lp)+i] = -b[I[i]]
        else
            V[i] = 1
            x[ncols(lp)+i] = b[I[i]]
        end
    end
    S = sparse(I, J, V, lp_nrows, nartif)

    if VT <: Array
        # create the matrix in canonical form
        A = [sparse(Matrix(lp.A)) S]
    elseif VT <: CuArray
        # on GPU
        A = [lp.A Matrix(S)]
    else
        error("Not supported vector type $(VT)")
    end

    artif_col_idx = Int[]
    for j = 1:nartif
        push!(artif_col_idx, ncols(lp)+j)
    end

    senses = fill(MatOI.EQUAL_TO, nrows(lp))
    return MatOI.LPSolverForm{T, typeof(A), typeof(c)}(
        MOI.MIN_SENSE,
        c, A, lp.b, senses, xl, xu
    ), artif_col_idx
end

function run(prob::MatOI.LPSolverForm, Tv::Type; kwargs...)::Simplex.SpxData

    if !haskey(kwargs, :pivot_rule)
        @error "Argument :pivot_rule is not provided."
    end
    pivot_rule = kwargs[:pivot_rule]

    # convert to phase-one form
    p1lp, artif = reformulate(prob)

    # load the problem
    spx = Simplex.SpxData(p1lp, Tv)
    spx.pivot_rule = deepcopy(pivot_rule)

    # set basis
    basic = collect((ncols(p1lp)-nrows(p1lp)+1):ncols(p1lp))
    Simplex.set_basis(spx, basic)

    # Run simplex method
    Simplex.run_core(spx)

    for j in artif
        spx.lpdata.v_lb[j] = 0.0
        spx.lpdata.v_ub[j] = 0.0
    end

    return spx
end

end # module