text:
######################### Basic BED file tools #########################

## .bed file format ##################################################
#
# See http://pngu.mgh.harvard.edu/~purcell/plink/binary.shtml
# http://zzz.bwh.harvard.edu/plink/binary.shtml
# and https://www.cog-genomics.org/plink2/formats#bed
#
######################################################################

"""
    rawformat(snp::Integer, quartermap=quarterstohuman)

Returns RAW format from the BED format quarter-byte:

| BED    | RAW (default) | general         | Meaning                 |
|:------ |:------------- |:--------------- |:----------------------- |
| `0b00` | `0b10`        | `quartermap[1]` | homozygous minor allele |
| `0b11` | `0b00`        | `quartermap[2]` | homozygous major allele |
| `0b01` | `NA_byte`     | `quartermap[3]` | missing value           |
| `0b10` | `0b01`        | `quartermap[4]` | heterozygous            |

"""
@inline function rawformat(snp::Integer, quartermap=Consts.quarterstohuman)
    @inbounds return quartermap[snp + 1]
end

"""
    breakbyte(byte::UInt8, bytemap=Consts.defaultbytemap)

Return a `Tuple{Vararg{UInt8, 4}}` of the RAW-formatted SNP quarters, as
determined by `bytemap`, in `byte`.

"""
@inline function breakbyte(byte::UInt8, bytemap=Consts.defaultbytemap, flip::Bool=false)
    @inbounds return bytemap[byte + 1, flip + 1]
end

@inline function breakbyte(byte::UInt8, bytemap, snpind::Int, flip::Bool=false)
    @inbounds return breakbyte(byte, bytemap, flip)[snpind]
end

# Approach from Quantgen/BEDMatrix: shift by 2(n - 1) and then mask
# result.
#
# Performance-wise, the lookup table approach of breakbyte is a bit
# faster.
"""
    quarter(byte::UInt8, n::Integer, quartermap=quarterstohuman)

Returns the RAW-formatted `n`th quarter of `byte`, equivalent to
`breakbyte(byte)[n]`.

"""
function quarter(byte::UInt8, n::Integer, quartermap=Consts.quarterstohuman)
    @inbounds return quartermap[((byte >> 2(n - 1)) & 0b00000011) + 1]
end

@inline function flipBEDbyte(byte::UInt8)
    # 0b00 -> 0b11
    # 0b11 -> 0b00
    # 0b10 -> 0b10
    # 0b01 -> 0b01

    @inbounds return Consts.flipBEDbytes[byte + 1]
end

"""
    BEDsize(filebase::AbstractString)

Returns the number of rows and columns (as a tuple) of the BED matrix
in `filebase*".bed"` as determined by linecounts of `filebase*".fam"`
and `filebase*".bim"`.

"""
function BEDsize(filebase::AbstractString)
    famfile, bimfile = filebase*".fam", filebase*".bim"
    isfile(famfile) || error("missing FAM file \"$famfile\"")
    isfile(bimfile) || error("missing BIM file \"$bimfile\"")

    n = countlines(famfile)
    p = countlines(bimfile)
    return n, p
end

"""
    checkmagic(bytes)

Determine if first two "magic" bytes match plink format. If not throws
an error, else returns `true`.

If bytes are UInt16, then it should be in big endian format.

"""
# Puts first two bytes into big endian format
checkmagic(bytes::Vector{UInt8}) = checkmagic((convert(UInt16, bytes[1]) << 8) + bytes[2])

checkmagic(twobytes::UInt16) = (twobytes == Consts.plinkmagic ||
                                error("Bad magic: not a plink bed file"))

function checkmagic(bedstream::IO)
    seekstart(bedstream)

    # Note that I don't just read(bedstream, UInt16) to avoid system
    # endian dependence, although I believe julia currently only runs
    # on little endian machines.
    return checkmagic((convert(UInt16, read(bedstream, UInt8)) << 8) + read(bedstream, UInt8))
end

"""
    BEDmode(byte)

Returns the BED mode as determined by the third byte of the file:
either the current standard, `:SNPmajor`, or older plink formats,
`:SNPminor`.

"""
BEDmode(byte::UInt8) = Consts.modes[byte]

function BEDmode(bedstream::IO)
    seek(bedstream, 2)
    return BEDmode(read(bedstream, UInt8))
end

BEDmode(bytevector::Vector{UInt8}) = BEDmode(bytevector[3])

getbytemap{T}(navalue::T) = getbytemap((convert(T, 0b10), navalue,
                                        convert(T, 0b01), convert(T, 0b00)))
function getbytemap{T}(quartermap::Tuple{T, T, T, T})
    bytemap = hcat([(rawformat(byte & 0b00_00_00_11, quartermap),
                     rawformat((byte & 0b00_00_11_00) >> 2, quartermap),
                     rawformat((byte & 0b00_11_00_00) >> 4, quartermap),
                     rawformat((byte & 0b11_00_00_00) >> 6, quartermap)) for
                    byte in 0x00:0xff],
                   [(rawformat(byte & 0b00_00_00_11, quartermap),
                     rawformat((byte & 0b00_00_11_00) >> 2, quartermap),
                     rawformat((byte & 0b00_11_00_00) >> 4, quartermap),
                     rawformat((byte & 0b11_00_00_00) >> 6, quartermap)) for
                    byte in Consts.flipBEDbytes])

   bytemap
end

"""
    unsafe_breakbyte!(vector::AbstractArray, byte::UInt8,
                      bytemap=bytetoquarters, vecstart=1,
                      quarterstart=1, quarterstop=4,
                      flip::Bool=false)

Fills `vector` starting at `vecstart` with the `quarterstart`th snp
through the `quarterstop`th snp in `byte`, as determined by `bytemap`.

## (Unchecked) Constraints:
* 1 <= `quarterstart` <= `quarterstop` <= 4
* 1 <= `vecstart` <= `end - (quarterstop - quarterstart)`

"""
function unsafe_breakbyte!(vector::AbstractArray, byte::UInt8,
                           bytemap=Consts.bytetoquarters, flip::Bool=false,
                           vecstart::Integer=1, quarterstart::Integer=1, quarterstop::Integer=4)

    @inbounds copy!(vector, vecstart,
                    breakbyte(byte, bytemap, flip), quarterstart,
                    quarterstop - quarterstart + 1)
    vector
end

"""
    unsafe_copybytestosnps!(snparray::AbstractArray, bytearray::AbstractArray{UInt8},
                            bytestart::Integer, quarterstart::Integer,
                            bytestop::Integer, quarterstop::Integer,
                            deststart::Integer=1, bytemap=bytetoquarters,
                            flip::Bool=false)

Fills `snparray[deststart:(deststart + num_snps)]` with snps, where

    num_snps = 4*(bytestop - bytestart) + (quarterstop - quarterstart) + 1

The SNPs start with the `quarterstart`th snp in the `bytestart`th
byte, and end with `quarterstop` and `bytestop`.  The SNP
representation is determined by the `bytemap`. Currently only supports
unit stride for both `snparray` and `bytearray`, since at present
there is no compelling use-case involving non-unit strides.

Utility function used by more front-facing functions.

If `flip` is `true`, then minor and major allele are flipped when copying.

"""
function unsafe_copybytestosnps!(snparray::AbstractArray, bytearray::AbstractArray{UInt8},
                                 bytestart::Integer, quarterstart::Integer,
                                 bytestop::Integer, quarterstop::Integer,
                                 deststart::Integer=1, bytemap=Consts.bytetoquarters,
                                 flip::Bool=false)

    # First byte
    if quarterstart != 1
        stop = ifelse(bytestart == bytestop, quarterstop, 4)
        @inbounds unsafe_breakbyte!(snparray, bytearray[bytestart], bytemap, flip,
                                    deststart, quarterstart, stop)
        deststart += stop - quarterstart + 1
        bytestart += 1
    end

    if bytestart <= bytestop
        # Last byte
        if quarterstop != 4
            @inbounds unsafe_breakbyte!(snparray, bytearray[bytestop], bytemap, flip,
                                        deststart + 4*(bytestop - bytestart), 1, quarterstop)
            bytestop -= 1
        end

        # Main course
        @simd for bytej in bytestart:bytestop
            @inbounds unsafe_breakbyte!(snparray, bytearray[bytej], bytemap, flip,
                                        deststart + 4*(bytej - bytestart))
        end
    end

    return snparray
end

######################### Reading .bed into a Matrix #########################

"""
    BEDintomatrix{T<:Real}(bedfilename::AbstractString; datatype::Type{T}=UInt8,
                           navalue=NA_byte, n::Integer=0, p::Integer=0, use_mmap=true)

Returns a `Matrix{T}` representation of the BED file `bedfilename`. If
the number of rows or columns, `n` and `p`, are not provided, then
they will attempt to be determined from corresponding .fam and .bim
files. `use_mmap` determines whether to memory map the bedfile and
then read into the matrix for potential speedups. Missing values are
set to `navalue`.

"""
function BEDintomatrix{T<:Real}(bedfilename::AbstractString;
                                datatype::Type{T}=UInt8,
                                navalue=NA_byte,
                                n::Integer=0,
                                p::Integer=0,
                                use_mmap::Bool=true)
    if !isfile(bedfilename)
        isfile(bedfilename*".bed") || error("Cannot find file \"$bedfilename\"")
        bedfilename = bedfilename*".bed"
    end
    filebase = splitext(bedfilename)[1]

    if n == p == 0
        n, p = BEDsize(filebase)
    end

    A = Matrix{T}(n, p)

    fsz = filesize(bedfilename)
    if use_mmap && 0 < fsz < typemax(Int)
        bedfilevector = open(bedfilename, "r") do bedfile
            bvec = Mmap.mmap(bedfile, Vector{UInt8}, (Int(fsz),))
            checkmagic(bvec)
            BEDmode(bvec) == :SNPmajor || error("SNPminor mode not supported")
            bvec
        end
        BEDintomatrix!(A, bedfilevector, convert(T, navalue))
    else
        open(bedfilename, "r") do bedfin
            checkmagic(bedfin)
            BEDmode(bedfin) == :SNPmajor || error("SNPminor mode not supported")

            BEDintomatrix!(A, bedfin, convert(T, navalue))
        end
    end

    return A::Matrix{T}
end

function BEDintomatrix!{T<:Real}(A::AbstractMatrix{T}, bedvector::Vector{UInt8}, navalue=NA_byte)
    n, p = size(A)
    byteheight = ceil(Int, n/4)
    quarterstop = n - 4*(byteheight - 1)

    bytemap = getbytemap(convert(T, navalue))

    @inbounds for col in 1:p
        unsafe_copybytestosnps!(A, bedvector, byteheight*(col - 1) + 3 + 1, 1, byteheight*col + 3, quarterstop, n*(col - 1) + 1, bytemap)
    end

    return A
end

"""
    BEDintomatrix!{T<:Real}(A::AbstractMatrix{T}, bedstream::IO, navalue=NA_byte)

Fills `A` with type `T` representation of `bedstream` that corresponds
to .bed file format. `A` must have correct dimensions, as determined
via `BEDMatrices.BEDsize`, for example.

"""
function BEDintomatrix!{T<:Real}(A::AbstractMatrix{T}, bedstream::IO, navalue=NA_byte)
    n, p = size(A)
    bytestop = ceil(Int, n/4)
    quarterstop = n - 4*(bytestop - 1)

    bytemap = getbytemap(convert(T, navalue))

    bytecol = Vector{UInt8}(bytestop)
    @inbounds for col in 1:p
        read!(bedstream, bytecol)
        unsafe_copybytestosnps!(A, bytecol, 1, 1, bytestop, quarterstop, n*(col - 1) + 1, bytemap)
    end

    return A
end

######################### BEDMatrix type #########################
# See implementation of BitArray for inspiration julia/base/bitarray.jl
# in /usr/share/julia

## Things we get for "free", but may want to overload for better performance, etc.:
# * checkbounds(B, i...)
# * eltype(B)
# * length(B)
# * ndims(B)
# * strides(B)
# * view(B, i...) -- not sure if it is memory and computationally efficient...
# * matrix operations: v * B and B * v, norm(B), etc.

# Allowing for `X` to be a general subtype `S` of `AbstractMatrix`
# requires exposing the type as below. See discussion at
# http://docs.julialang.org/en/stable/manual/performance-tips/#type-declarations

immutable BEDMatrix{T, S<:AbstractMatrix} <: DenseArray{T, 2}
    n::Int                               # number of samples (rows)
    p::Int                               # number of SNPs (columns)
    X::S                                 # Matrix{UInt8} bed file
    navalue::T                           # representation of missing

    path::String                         # file location
    colnames::Vector{String}             # SNP labels, can be used for indexing
    rownames::Vector{String}             # sample labels, can be used for indexing
    colnamemap::Dict{String, Int}        # for quick name indexing by column
    rownamemap::Dict{String, Int}        # for quick name indexing by row

    _byteheight::Int                     # number of bytes in each column
    _lastrowSNPheight::Int               # number of SNPs in last byte of each column

    _bytemap::Matrix{Tuple{T, T, T, T}}  # quarters and flipped quarters for 0x00:0xff

    _flip::BitVector                     # whether to flip SNP major--minor
                                         # allele encoding for each SNP

    function BEDMatrix{T, S}(n::Integer, p::Integer, X::S, navalue::T,
                       path::AbstractString, colnames::AbstractVector, rownames::AbstractVector,
                       bytemap, flip::BitVector) where {T, S<:AbstractMatrix{UInt8}}
        byteheight = ceil(Int, n/4)
        lastrowheight = n - 4*(byteheight - 1)

        size(X) == (byteheight, p) || throw(DimensionMismatch("Matrix dimensions $(size(X)) do not agree with supplied BED dimensions (n = $n, p = $p)"))
        length(colnames) == p || throw(DimensionMismatch("colnames has incorrect length"))
        length(rownames) == n || throw(DimensionMismatch("rownames has incorrect length"))
        size(bytemap) == (256, 2) || error("bytemap must be of size (256, 2)")
        length(flip) == p || throw(DimensionMismatch("flip has length $(length(flip)); expected length of p = $p"))

        colnamemap = Dict{String, Int}()
        for (cidx, colname) in enumerate(colnames)
            colnamemap[colname] = cidx
        end

        rownamemap = Dict{String, Int}()
        for (ridx, rowname) in enumerate(rownames)
            rownamemap[rowname] = ridx
        end

        return new(n, p, X, navalue, path, colnames, rownames, colnamemap, rownamemap,
                   byteheight, lastrowheight, bytemap, flip)
    end
end

"""
    BEDMatrix(bedfilename::AbstractString;
              datatype::DataType=Int8, nsamples::Integer=0, nSNPs::Integer=0, navalue=NA_byte
              famfile::AbstractString="", bimfile::AbstractString="",
              quartermap::Tuple=Consts.quarterstohuman, flip::AbstractVector=[])

Create a `BEDMatrix` of type `datatype` using memory mapping of BED
file, `bedfilename`. Use `navalue` for missing values. `famfile` and
`bimfile` are only required if the .fam and .bim file are not in same
location as or do not share the same base name as the
`bedfilename`. `nsamples` and `nSNPs` are only used if .bim and .fam
file cannot be found and are not provided; in which case the
rownames/colnames will be generic.

`convert(datatype, x)` must work for `x` in `[0b00, 0b01, 0b10, navalue]`.

## Optional arguments
 - `datatype`: the return type for indexing.
 - `nsamples`: number of samples if .fam file cannot be found
 - `nSNPs`: the number of SNPs if .bim file cannot be found
 - `navalue`: the representation of missing values; this
   should not be used with quartermap.
 - `famfile`: the .fam file, only needed if the .fam file does
   not have the same base name as `bedfilename`.
 - `bimfile`: the .bim file, only needed if the .bim file does
   not have the same base name as `bedfilename`.
 - `quartermap`: a `Tuple` of the respective representation of
   homozygous minor, homozygous major, missing, and heterozygous
   values. Useful if you don't want to use standard RAW format
   encoding.
 - `flip`: a `Vector` of indices that should have the major--minor
   representation flipped. `flip` may either be a list of indices, a
   logical vector with `false` representing no flipping, or a list of
   columnnames. This may be recovered or changed later using
   `getflips` and `setflips!`, respectively.

## Examples
```julia
julia> bed = BEDMatrix("test/data/example.bed");

julia> bed[1:5, 1:5]
5×5 Array{Int8,2}:
 0  1  1  1  0
 1  1  1  1  3
 1  0  0  2  0
 2  0  0  0  1
 0  1  0  0  0

julia> bed = BEDMatrix("test/data/example", datatype=Nullable{Int}, navalue=Nullable{Int}());

julia> bed[1:5, 1:5]
5×5 Array{Nullable{Int64},2}:
 0  1  1  1  0
 1  1  1  1  #NULL
 1  0  0  2  0
 2  0  0  0  1
 0  1  0  0  0

julia> bed = BEDMatrix("test/data/example.bed", datatype=Float64, navalue=NaN);

julia> bed[1:5, 1:5]
5×5 Array{Float64,2}:
 0.0  1.0  1.0  1.0    0.0
 1.0  1.0  1.0  1.0  NaN
 1.0  0.0  0.0  2.0    0.0
 2.0  0.0  0.0  0.0    1.0
 0.0  1.0  0.0  0.0    0.0

julia> bed = BEDMatrix("test/data/example", famfile="badfilename", nsamples=50);

julia> rownames(bed)[1:5]
5-element Array{String,1}:
 "sample_1"
 "sample_2"
 "sample_3"
 "sample_4"
 "sample_5"
```

"""
function BEDMatrix(bedfilename::AbstractString;
                   datatype::DataType=Int8,
                   nsamples::Integer=0, nSNPs::Integer=0, navalue=NA_byte,
                   famfile::AbstractString="", bimfile::AbstractString="",
                   quartermap::Tuple=Consts.quarterstohuman, flip::AbstractVector=Int[])

    if !isfile(bedfilename)
        isfile(bedfilename*".bed") || error("Cannot find file \"$bedfilename\"")
        bedfilename = bedfilename*".bed"
    end
    filebase = splitext(bedfilename)[1]

    rownames = create_rownames(filebase, famfile, nsamples)
    colnames = create_colnames(filebase, bimfile, nSNPs)
    n, p = length(rownames), length(colnames)

    X = open(filebase*".bed", "r") do bedfile
        checkmagic(bedfile)
        BEDmode(bedfile) == :SNPmajor || error("Old-style SNP-minor mode bed"*
                                               " files not supported. Use plink"*
                                               " to convert to SNP-major format")

        Mmap.mmap(bedfile, Matrix{UInt8}, (ceil(Int, n/4), p))
    end

    # validate quartermap and navalue
    if quartermap == Consts.quarterstohuman
        bytemap = getbytemap(convert(datatype, navalue))
    else
        bytemap = getbytemap(quartermap)
        navalue = quartermap[2]
    end

    # Process and validate flip
    flip = process_flip(flip, colnames)

    return BEDMatrix{datatype, typeof(X)}(n, p, X, convert(datatype, navalue),
                                          abspath(bedfilename), colnames, rownames, bytemap, flip)
end

function parsefamline(line)
    fields = split(line)

    (fields[1], fields[2], parse(Int, fields[3]), parse(Int, fields[4]), parse(Int, fields[5]), parse(Int, fields[6]))
end

famrowname(line) = join(split(line)[1:2], '_')

readrownames(famfile::AbstractString) = map(famrowname, readlines(famfile))

function create_rownames(filebase::AbstractString, famfile::AbstractString="", nsamples::Integer=0)
    famfile = famfile == "" ? filebase*".fam" : famfile
    if isfile(famfile)
        rownames = readrownames(famfile)
    elseif nsamples > 0
        rownames = [string("sample_", j) for j in 1:nsamples]
    else
        error("Cannot find FAM file \"$famfile\" and nsamples not provided")
    end

    return rownames
end

function bimcolname(line)
    fields = split(line)
    string(fields[2], '_', fields[5])
end

readcolnames(bimfile::AbstractString) = map(bimcolname, readlines(bimfile))

function create_colnames(filebase::AbstractString, bimfile::AbstractString="", nSNPs::Integer=0)
    bimfile = bimfile == "" ? filebase*".bim" : bimfile

    if isfile(bimfile)
        colnames = readcolnames(bimfile)
    elseif nSNPs > 0
        colnames = [string("SNP_", j) for j in 1:nSNPs]
    else
        error("Cannot find BIM file \"$bimfile\" and nSNPs not provided")
    end
end

function process_flip(flip, colnames)
    p = length(colnames)

    if eltype(flip) == Bool
        length(flip) == p || throw(DimensionMismatch("flip has length $(length(flip)) expected length of p = $p"))
        flip = BitVector(flip)
    elseif eltype(flip) <: Integer
        checkbounds(Bool, 1:p, flip) || error("Invalid index in flip vector")
        flip = BitVector(map(x -> x in flip, 1:p))
    elseif eltype(flip) <: AbstractString
        for name in flip
            name in colnames || error("\"$name\" in flip is not a column_name")
        end
        flip = BitVector(map(x -> x in flip, colnames))
    else
        error("Invalid flip vector: expected logical, index, or column name vector")
    end

    return flip
end

######################## Getters, Properties, and Setters ########################

Base.size(B::BEDMatrix) = (B.n, B.p)
function Base.size(B::BEDMatrix, k::Integer)
    k > 0 || Base.throw_boundserror(size(B), k)
    # This is the same behavior as existing Base.size methods
    return k == 1 ? B.n : ifelse(k == 2, B.p, 1)
end

Base.IndexStyle{T<:BEDMatrix}(::Type{T}) = Base.IndexCartesian()

# not sure if this is the right definition; see
# https://github.com/JuliaLang/julia/issues/16614.
# This is the size of `B[:, :]`.
Base.sizeof{T, S}(B::BEDMatrix{T, S}) = B.n*B.p*sizeof(T)

"""
    path(B::BEDMatrix)

Returns the location of the .bed file used to create `X`.

"""
path(B::BEDMatrix) = B.path

"""
    rownames(B::BEDMatrix)

Returns `Vector{String}` of row names. These may be used for indexing.

"""
rownames(B::BEDMatrix) = B.rownames

"""
    colnames(B::BEDMatrix)

Returns `Vector{String}` of column names. These may be used for indexing.

"""
colnames(B::BEDMatrix) = B.colnames

"""
    NArep(B::BEDMatrix)

Returns value used for missing entries.

"""
NArep(B::BEDMatrix) = B.navalue

getrow(B::BEDMatrix, rowname::AbstractString) = B.rownamemap[rowname]
getcol(B::BEDMatrix, colname::AbstractString) = B.colnamemap[colname]

function _get_quartermap_from_bytemap(bytemap::Matrix)
    # 0b11_10_01_00 is the BED encoding for (homozygous minor,
    # homozygous major, missing, heterozygous). Recall the reverse
    # (little endian-like) order of the SNPs.
    bytemap[0b11_10_01_00 + 1, false + 1]
end

"""
    getquartermap(bed::BEDMatrix)

Returns the `quartermap` for `bed`. This is a `Tuple` of four values
giving the representation of respectively (homozygous minor,
homozygous major, missing value, heterozygous).

"""
function getquartermap(bed::BEDMatrix)
    _get_quartermap_from_bytemap(bed._bytemap)
end

"""
    getflips(B::BEDMatrix)

Returns a `Vector{Int}` of the columns which have a flipped major--minor allele
representation. Use `setflips!` to change this.

"""
getflips(B::BEDMatrix) = find(B._flip)

"""
    setflips!(B::BEDMatrix, indices)

(Re)sets which columns of the BEDMatrix use a flipped major--minor
representation. For example, `setflips!(bed, [])` would make all columns use the
standard encoding and `setflips!(bed, indices(bed, 2))` would make all columns
use a flipped representation.

"""
function setflips!(B::BEDMatrix, flip_logical::AbstractVector{Bool})
    length(flip_logical) == size(B, 2) || throw(DimensionsMismatch("flip had incorrect length"))
    copy!(B._flip, flip_logical)
    B._flip
end

function setflips!{K<:Integer}(B::BEDMatrix, flip_indices::AbstractVector{K})
    checkbounds(1:B.p, flip_indices)

    fill!(B._flip, false)

    @inbounds for col in flip_indices
        B._flip[col] = true
    end

    B._flip
end

function setflips!{S<:AbstractString}(B::BEDMatrix, flip_names::AbstractVector{S})
    # I don't do this in the same loop as below, since I want an
    # atomic operation.
    for name in flip
        name in colnames || error("\"$name\" in flip is not a column_name")
    end
    fill!(B._flip, false)

    for name in flip
        B._flip[getcol(B, name)] = true
    end

    B._flip
end

#################### Indexing ####################

Base.getindex{T<:AbstractString}(B::BEDMatrix, rownames::AbstractVector{T}, col) = B[[getrow(B, rowname) for rowname in rownames], col]
Base.getindex{T<:AbstractString}(B::BEDMatrix, row, colnames::AbstractVector{T}) = B[row, [getcol(B, colname) for colname in colnames]]
Base.getindex{T<:AbstractString, S<:AbstractString}(B::BEDMatrix, rownames::AbstractVector{S}, colnames::AbstractVector{T}) = [B[row, col] for row in rownames, col in colnames]

Base.getindex(B::BEDMatrix, rowname::AbstractString, colname::AbstractString) = B[getrow(B, rowname), getcol(B, colname)]
Base.getindex(B::BEDMatrix, rowname::AbstractString, col) = B[getrow(B, rowname), col]
Base.getindex(B::BEDMatrix, row, colname::AbstractString) = B[row, getcol(B, colname)]

Base.getindex{T<:AbstractString}(B::BEDMatrix, rownames::AbstractVector{T}, colname::AbstractString) = [B[row, colname] for row in rownames]
Base.getindex{T<:AbstractString}(B::BEDMatrix, rowname::AbstractString, colnames::AbstractVector{T}) = [B[rowname, col] for col in colnames]

# This is is the only getindex method we _need_, the other methods are
# provided for better performance, or convenience.
function Base.getindex(B::BEDMatrix, row::Integer, col::Integer)
    @boundscheck checkbounds(B, row, col)

    unsafe_getindex(B, row, col)
end

function Base.getindex(B::BEDMatrix, ::Colon, col::Integer)
    @boundscheck checkbounds(B, :, col)

    unsafe_getcol(B, col)
end

function Base.getindex{T, S, K<:Integer}(B::BEDMatrix{T, S}, ::Colon, cols::AbstractVector{K})
    @boundscheck checkbounds(B, :, cols)

    matrix = Matrix{T}(B.n, length(cols))

    for (cidx, col) in enumerate(cols)
        unsafe_getcol!(matrix, B, col, (cidx - 1)*B.n + 1)
    end
    return matrix
end

function Base.getindex(B::BEDMatrix, rrange::AbstractUnitRange, col::Integer)
    @boundscheck checkbounds(B, rrange, col)

    unsafe_getrowrange(B, rrange, col)
end

function Base.getindex{T, S, K <: Integer}(B::BEDMatrix{T, S}, rrange::AbstractUnitRange, cols::AbstractVector{K})
    @boundscheck checkbounds(B, rrange, cols)

    matrix = Matrix{T}(length(rrange), length(cols))
    for (cidx, col) in enumerate(cols)
        unsafe_getrowrange!(matrix, B, rrange, col, (cidx - 1)*length(rrange) + 1)
    end
    return matrix
end

function unsafe_getindex{T, S}(B::BEDMatrix{T, S}, row::Integer, col::Integer)
    byterow, snpind = rowtobytequarter(row)

    @inbounds return breakbyte(B.X[byterow, col], B._bytemap, snpind, B._flip[col])
end

function unsafe_getcol{T, S}(B::BEDMatrix{T, S}, col::Integer)
    column = Vector{T}(B.n)

    unsafe_getcol!(column, B, col)
    return column
end

function unsafe_getcol!{T, S}(column::AbstractArray{T}, B::BEDMatrix{T, S}, col::Integer, deststart=1)
    @inbounds unsafe_copybytestosnps!(column, B.X, B._byteheight*(col - 1) + 1, 1, B._byteheight*col,
                                      B._lastrowSNPheight, deststart, B._bytemap, B._flip[col])
    return column
end

function unsafe_getrowrange{T, S}(B::BEDMatrix{T, S}, rrange::AbstractUnitRange, col::Integer)
    vector = Vector{T}(length(rrange))
    unsafe_getrowrange!(vector, B, rrange, col)
    return vector
end

function unsafe_getrowrange!{T, S}(vector::AbstractArray{T}, B::BEDMatrix{T, S}, rrange::AbstractUnitRange, col::Integer, deststart=1)
    bytestart, quarterstart = rowtobytequarter(first(rrange))
    bytestart += B._byteheight*(col - 1)

    bytestop, quarterstop = rowtobytequarter(last(rrange))
    bytestop += B._byteheight*(col - 1)

    @inbounds unsafe_copybytestosnps!(vector, B.X, bytestart, quarterstart, bytestop, quarterstop, deststart, B._bytemap, B._flip[col])
    return vector
end

function rowtobytequarter(row::Integer)
    # Curse you julia for your foolish 1-indexing!
    return (((row - 1) >> 2) + 1, ((row - 1) & 3) + 1)
end

"""
    getquarterblock(B::BEDMatrix, row::Integer, col::Integer)

Returns `(fourquarters, snpind)`, where `fourquarters` is the block of
quarters that `B[row, col]` belongs to: `B[row, col] =
fourquarters[snpind]`.

"""
function getquarterblock(B::BEDMatrix, row::Integer, col::Integer)
    byterow, snpind = rowtobytequarter(row)

    @inbounds return (breakbyte(B.X[byterow, col], B._bytemap, B._flip[col]), byterow, snpind)
end


text:
#
# ccalls.jl -
#
# Calls to C-API of Nüvü Camēras.
#

if isfile(joinpath(@__DIR__,"..","deps","deps.jl"))
    include("../deps/deps.jl")
else
    error("NuvuCameras not properly installed.  Please run `Pkg.build(\"NuvuCameras\")` to create file \"",joinpath(@__DIR__,"..","deps","deps.jl"),"\"")
end

ext: jl


text:
module EllipsisNotation

  import Base.getindex, Base.setindex!
  const ..   = Val{:..}

  setindex!{T}(A::AbstractArray{T,1}, x, ::Type{Val{:..}}, n::Int) = A[n] = x
  setindex!{T}(A::AbstractArray{T,2}, x, ::Type{Val{:..}}, n::Int) = A[ :, n] = x
  setindex!{T}(A::AbstractArray{T,3}, x, ::Type{Val{:..}}, n::Int) = A[ :, :, n] =x
  setindex!{T}(A::AbstractArray{T,4}, x, ::Type{Val{:..}}, n::Int) = A[ :, :, :, n] =x
  setindex!{T}(A::AbstractArray{T,5}, x, ::Type{Val{:..}}, n::Int) = A[ :, :, :, :, n] =x
  setindex!{T}(A::AbstractArray{T,6}, x, ::Type{Val{:..}}, n::Int) = A[ :, :, :, :, :, n] =x
  setindex!{T}(A::AbstractArray{T,7}, x, ::Type{Val{:..}}, n::Int) = A[ :, :, :, :, :, :, n] =x
  setindex!{T}(A::AbstractArray{T,8}, x, ::Type{Val{:..}}, n::Int) = A[ :, :, :, :, :, :, :, n] =x
  setindex!{T}(A::AbstractArray{T,9}, x, ::Type{Val{:..}}, n::Int) = A[ :, :, :, :, :, :, :, :, n] =x
  setindex!{T}(A::AbstractArray{T,10}, x, ::Type{Val{:..}}, n::Int) = A[ :, :, :, :, :, :, :, :, :, n] =x

  getindex{T}(A::AbstractArray{T,1}, ::Type{Val{:..}}, n::Int) = A[n]
  getindex{T}(A::AbstractArray{T,2}, ::Type{Val{:..}}, n::Int) = A[ :, n]
  getindex{T}(A::AbstractArray{T,3}, ::Type{Val{:..}}, n::Int) = A[ :, :, n]
  getindex{T}(A::AbstractArray{T,4}, ::Type{Val{:..}}, n::Int) = A[ :, :, :, n]
  getindex{T}(A::AbstractArray{T,5}, ::Type{Val{:..}}, n::Int) = A[ :, :, :, :, n]
  getindex{T}(A::AbstractArray{T,6}, ::Type{Val{:..}}, n::Int) = A[ :, :, :, :, :, n]
  getindex{T}(A::AbstractArray{T,7}, ::Type{Val{:..}}, n::Int) = A[ :, :, :, :, :, :, n]
  getindex{T}(A::AbstractArray{T,8}, ::Type{Val{:..}}, n::Int) = A[ :, :, :, :, :, :, :, n]
  getindex{T}(A::AbstractArray{T,9}, ::Type{Val{:..}}, n::Int) = A[ :, :, :, :, :, :, :, :, n]
  getindex{T}(A::AbstractArray{T,10}, ::Type{Val{:..}}, n::Int) = A[ :, :, :, :, :, :, :, :, :, n]

  setindex!{T}(A::AbstractArray{T,1}, x, n::Int, ::Type{Val{:..}}) = A[n] = x
  setindex!{T}(A::AbstractArray{T,2}, x, n::Int, ::Type{Val{:..}}) = A[n, :] = x
  setindex!{T}(A::AbstractArray{T,3}, x, n::Int, ::Type{Val{:..}}) = A[n, :, :] = x
  setindex!{T}(A::AbstractArray{T,4}, x, n::Int, ::Type{Val{:..}}) = A[n, :, :, :] = x
  setindex!{T}(A::AbstractArray{T,5}, x, n::Int, ::Type{Val{:..}}) = A[n, :, :, :, :] = x
  setindex!{T}(A::AbstractArray{T,6}, x, n::Int, ::Type{Val{:..}}) = A[n, :, :, :, :, :] = x
  setindex!{T}(A::AbstractArray{T,7}, x, n::Int, ::Type{Val{:..}}) = A[n, :, :, :, :, :, :] = x
  setindex!{T}(A::AbstractArray{T,8}, x, n::Int, ::Type{Val{:..}}) = A[n, :, :, :, :, :, :, :] = x
  setindex!{T}(A::AbstractArray{T,9}, x, n::Int, ::Type{Val{:..}}) = A[n, :, :, :, :, :, :, :, :] = x
  setindex!{T}(A::AbstractArray{T,10}, x, n::Int, ::Type{Val{:..}}) = A[n, :, :, :, :, :, :, :, :, :] = x

#=
  setindex!{T}(A::AbstractArray{T,1}, x, n::Int, ::Type{Val{:...}}) = A[n] = x
  setindex!{T}(A::AbstractArray{T,2}, x, n::Int, ::Type{Val{:...}}) = A[n, :] = reshape(x,1,size(x,1)) #x
  setindex!{T}(A::AbstractArray{T,3}, x, n::Int, ::Type{Val{:...}}) = A[n, :, :] = reshape(x,1,size(x,1),size(x,2)) #x
  setindex!{T}(A::AbstractArray{T,4}, x, n::Int, ::Type{Val{:...}}) = A[n, :, :, :] = reshape(x,1,size(x,1),size(x,2),size(x,3)) #x
  setindex!{T}(A::AbstractArray{T,5}, x, n::Int, ::Type{Val{:...}}) = A[n, :, :, :, :] = reshape(x,1,size(x,1),size(x,2),size(x,3),size(x,4)) #x
  setindex!{T}(A::AbstractArray{T,6}, x, n::Int, ::Type{Val{:...}}) = A[n, :, :, :, :, :] = reshape(x,1,size(x,1),size(x,2),size(x,3),size(x,4),size(x,5)) #x
  setindex!{T}(A::AbstractArray{T,7}, x, n::Int, ::Type{Val{:...}}) = A[n, :, :, :, :, :, :] =
  reshape(x,1,size(x,1),size(x,2),size(x,3),size(x,4),size(x,5),size(x,6)) #x
  setindex!{T}(A::AbstractArray{T,8}, x, n::Int, ::Type{Val{:...}}) = A[n, :, :, :, :, :, :, :] =
  reshape(x,1,size(x,1),size(x,2),size(x,3),size(x,4),size(x,5),size(x,6),size(x,7)) #x
  setindex!{T}(A::AbstractArray{T,9}, x, n::Int, ::Type{Val{:...}}) = A[n, :, :, :, :, :, :, :, :] =
  reshape(x,1,size(x,1),size(x,2),size(x,3),size(x,4),size(x,5),size(x,6),size(x,7),size(x,8)) #x
  setindex!{T}(A::AbstractArray{T,10}, x, n::Int, ::Type{Val{:...}}) = A[n, :, :, :, :, :, :, :, :, :] =
  reshape(x,1,size(x,1),size(x,2),size(x,3),size(x,4),size(x,5),size(x,6),size(x,7),size(x,8),size(x,9)) #x
  =#

  getindex{T}(A::AbstractArray{T,1}, n::Int, ::Type{Val{:..}}) = A[n]
  getindex{T}(A::AbstractArray{T,2}, n::Int, ::Type{Val{:..}}) = A[n, :]
  getindex{T}(A::AbstractArray{T,3}, n::Int, ::Type{Val{:..}}) = A[n, :, :]
  getindex{T}(A::AbstractArray{T,4}, n::Int, ::Type{Val{:..}}) = A[n, :, :, :]
  getindex{T}(A::AbstractArray{T,5}, n::Int, ::Type{Val{:..}}) = A[n, :, :, :, :]
  getindex{T}(A::AbstractArray{T,6}, n::Int, ::Type{Val{:..}}) = A[n, :, :, :, :, :]
  getindex{T}(A::AbstractArray{T,7}, n::Int, ::Type{Val{:..}}) = A[n, :, :, :, :, :, :]
  getindex{T}(A::AbstractArray{T,8}, n::Int, ::Type{Val{:..}}) = A[n, :, :, :, :, :, :, :]
  getindex{T}(A::AbstractArray{T,9}, n::Int, ::Type{Val{:..}}) = A[n, :, :, :, :, :, :, :, :]
  getindex{T}(A::AbstractArray{T,10}, n::Int, ::Type{Val{:..}}) = A[n, :, :, :, :, :, :, : ,: ,:]

  export .., getindex, setindex!

end # module


text:
#############################################################################
#  Copyright 2016, Iain Dunning, Joey Huchette, Miles Lubin, and contributors
#  This Source Code Form is subject to the terms of the Mozilla Public
#  License, v. 2.0. If a copy of the MPL was not distributed with this
#  file, You can obtain one at http://mozilla.org/MPL/2.0/.
#
# Author: Louis Luangkesorn
# Date: Jan 30, 2015
#
# Allocation of passenger cars to trains to minimize cars required or car-miles run.
# Based on
#   Fourer, D.M. Gay and Brian W. Kernighan, A Modeling Language for Mathematical
#   Programming, http://www.ampl.com/REFS/amplmod.ps.gz
#   Appendix D
#############################################################################
using JuMP

ORIG = ["GARY", "CLEV", "PITT"];
DEST = ["FRA", "DET", "LAN", "WIN", "STL", "FRE", "LAF"];

supply = [1400 2600 2900];
demand = [ 900 1200  600 400 1700 1100 1000];

@assert sum(supply) == sum(demand);

cost = [
39   14   11   14   16   82    8;
27    9   12    9   26   95   17;
24   14   17   13   28   99   20
]

m = Model();

@variable(m, Trans[i=1:length(ORIG), j=1:length(DEST)] >= 0);

@objective(m, Min, sum{cost[i,j] * Trans[i,j], i=1:length(ORIG), j=1:length(DEST)});

@constraint(m, xyconstr[i=1:1:length(ORIG)], sum{Trans[i,j], j=1:length(DEST)} == supply[i]);

@constraint(m, xyconstr[j = 1:length(DEST)], sum{Trans[i,j], i=1:length(ORIG)} == demand[j]);

println("Solving original problem...")
status = solve(m);

if status == :Optimal
	@printf("Optimal!\n");
	@printf("Objective value: %d\n", getobjectivevalue(m));
	@printf("Transpotation:\n");
	for j = 1:length(DEST)
		@printf("\t%s", DEST[j]);
	end
	@printf("\n");
	for i = 1:length(ORIG)
		@printf("%s", ORIG[i]);
		for j = 1:length(DEST)
			@printf("\t%d", getvalue(Trans[i,j]));
		end
		@printf("\n");
	end
else
    @printf("No solution\n");
end


text:
function annotatesvg(p)
    ifstream = open(p[:outfile])
    iobuffer = IOBuffer()

    chr = p[:chromosome]
    features = p[:features]
    genes = @genes(chr, feature(gene) in features)

    ingenes = false
    stripedgene = false
    atgene = 0
    atline = 0
    for line in eachline(ifstream)
        atline += 1
        if atline == 2
            line *= "<title>" * p[:chromosome].name * "</title>"
        elseif atgene > length(genes)
            println(iobuffer, line)
            continue
        end
        if !stripedgene && isnewgene(line)
            ingenes = true
            stripedgene = false
            atgene += 1
            annotation = replace(replace(string(genes[atgene]), r"^ +" => ""), r"\n +" => "\n")
            line = line[1:end - 2] * ">\n<title>" * annotation * "</title></path>\n"
        elseif occursin(r"<path style=\"fill:none;stroke-width:[^;]+;stroke-linecap:[^;]+;stroke-linejoin:miter;stroke:rgb\([^)]+\);stroke-opacity:1;stroke-miterlimit:10;\" d=\"M \S+ \S+ L \S+ \S+ L \S+ \S+ L \S+ \S+ L \S+ \S+ L \S+ \S+ L \S+ \S+ \"", line)
            # Second half of spanning genes
            annotation = replace(replace(string(genes[atgene]), r"^ +" => ""), r"\n +" => "\n")
            line = line[1:end - 2] * ">\n<title>" * annotation * "</title></path>\n"
        elseif ingenes && occursin(r"<path", line)
            annotation = replace(replace(string(genes[atgene]), r"^ +" => ""), r"\n +" => "\n")
            line = line[1:end - 2] * ">\n<title>" * annotation * "</title></path>\n"
        end

        if ingenes && occursin(r"<g", line)
            stripedgene = true
        elseif ingenes && occursin(r"</g>", line)
            stripedgene = false
        end

        # leaving the genes section, no more annotations
        if ingenes && occursin(r"d=\"M \S+ \S+ C ", line)
            ingenes = false
        end

        println(iobuffer, line)
    end

    close(ifstream)

    ofstream = open(p[:outfile], "w")
    print(ofstream, String(take!(iobuffer)))
    close(ofstream)

    return nothing
end

function isnewgene(line)
    # Full, non-spanning
    if occursin(r"<path style=\"fill:none;stroke-width:[^;]+;stroke-linecap:[^;]+;stroke-linejoin:[^;]+;stroke:rgb\([^)]+\);stroke-opacity:1;stroke-miterlimit:10;\" d=\"M \S+ \S+ L \S+ \S+ L \S+ \S+ L \S+ \S+ L \S+ \S+ L \S+ \S+ L \S+ \S+ L \S+ \S+ \"", line)
        return true
    end
    # Short, non-spanning
    if occursin(r"<path style=\"fill:none;stroke-width:[^;]+;stroke-linecap:[^;]+;stroke-linejoin:[^;]+;stroke:rgb\([^)]+\);stroke-opacity:1;stroke-miterlimit:10;\" d=\"M \S+ \S+ L \S+ \S+ L \S+ \S+ L \S+ \S+ L \S+ \S+ \"", line)
        c = match(r"d=\"M (\S+) (\S+) L \S+ \S+ L \S+ \S+ L \S+ \S+ L (\S+) (\S+)", line).captures
        if c[1] == c[3] && c[2] == c[4]
           return true
        end
    end
    # First half, complement spanning, long
    if occursin(r"<path style=\"fill:none;stroke-width:[^;]+;stroke-linecap:[^;]+;stroke-linejoin:[^;]+;stroke:rgb\([^)]+\);stroke-opacity:1;stroke-miterlimit:10;\" d=\"M \S+ \S+ L \S+ \S+ L \S+ \S+ L \S+ \S+ L \S+ \S+ L \S+ \S+ L \S+ \S+ \"", line)
        c = match(r"d=\"M (\S+) \S+ L \S+ \S+ L \S+ \S+ L (\S+) \S+", line).captures
        if parse(Float64, c[1]) > parse(Float64, c[2])
            return true
        end
    end
    # First half complement spanning, short
    if occursin(r"<path style=\"fill:none;stroke-width:[^;]+;stroke-linecap:[^;]+;stroke-linejoin:[^;]+;stroke:rgb\([^)]+\);stroke-opacity:1;stroke-miterlimit:10;\" d=\"M \S+ \S+ L \S+ \S+ L \S+ \S+ L \S+ \S+ L \S+ \S+ \"", line)
        c = match(r"d=\"M (\S+) (\S+) L \S+ \S+ L (\S+) \S+ L \S+ \S+ L (\S+) (\S+)", line).captures
        if c[1] == c[4] && c[2] != c[5] && c[1] > c[3]
           return true
        end
    end
    # First half, non-compliment
    if occursin(r"<path style=\"fill:none;stroke-width:[^;]+;stroke-linecap:[^;]+;stroke-linejoin:[^;]+;stroke:rgb\([^)]+\);stroke-opacity:1;stroke-miterlimit:10;\" d=\"M \S+ \S+ L \S+ \S+ L \S+ \S+ L \S+ \S+ \"", line)
        c = match(r"d=\"M \S+ (\S+) L \S+ \S+ L \S+ \S+ L \S+ (\S+)", line).captures
        if parse(Float64, c[1]) > parse(Float64, c[2])
            return true
        end
    end
    return false
end


text:
using VectorSphericalWaves

"""
## Failure, doesn't support ^
using MultiFloat
MultiFloats.use_bigfloat_transcendentals()
"""
n_max, kr_r, kr_i, θ, ϕ = 2, 1e6, 1e2, 0.2, 0.3;
kr = complex(kr_r, kr_i);
M_N_wave_all_m_n(n_max, kr, θ, ϕ; kind="regular");

# BigFloat
n_max, kr_r, kr_i, θ, ϕ = 2, BigFloat(1e6), BigFloat(1e2), BigFloat(0.2), BigFloat(0.3);
kr = complex(kr_r, kr_i);
M_N_wave_all_m_n(n_max, kr, θ, ϕ; kind="regular");


text:
#=
Functions related to exchanging ghost data.
=#

"""
    ExchangeGroup

A description of information that needs to be communicated between processes during
a ghost exchange. It includes the ID of the involved processes, the included elements,
and the variable and component indices. The ghost exchange function will use this
to set up the data and path to communicate.
"""
mutable struct ExchangeGroup
    members::Array{Int,1} # ranks of the members of the group

    ghost_elements_only::Bool   # True if the exchange is for ghost element pairs
    exchange_all_data::Bool     # True if all variables/components are exchanged

    elements::Array{Int,2}  # 1xN list of global indices of elements to share. 2xN Pairs if ghost_elements_only.
    variables::Array{Int,1} # indices of variables to share (empty if exchange_all_data=true)
    components::Array{Array{Int,1},1} # components of each variable to share (empty if exchange_all_data=true)
end

"""

"""
function build_exchange_groups(grid; members=[], variables=[], components=[], ghost_elements_only=true, exchange_all_data=true)
    # This list of exchange groups will be returned
    ex_groups = [];

    if ghost_elements_only
        if exchange_all_data
            # If ghost_elements_only and exchange_all_data, this matches the ghost info from the grid
            ngroups = grid.num_neighbor_partitions;
            for gi=1:ngroups
                neighbor = grid.neighboring_partitions[gi];

                push!(ex_groups, ExchangeGroup([neighbor], ghost_elements_only, exchange_all_data, grid.ghost_index[gi], [], []));
            end
        else
            # If exchange_all_data = false, this matches the ghost info from the grid, but specifies some components
            ngroups = grid.num_neighbor_partitions;
            for gi=1:ngroups
                neighbor = grid.neighboring_partitions[gi];

                push!(ex_groups, ExchangeGroup([neighbor], ghost_elements_only, exchange_all_data, grid.ghost_index[gi], variables, components));
            end
        end
    else # all local elements are shared with ranks that have same partition
        if exchange_all_data

        else

        end
    end

    return ex_groups;
end

#=
Exchange ghosts.
This does both send and receive.
Values are sent/received directly from the variable.values arrays

NOTE: This is currently only set up for the case where each process has a distinct partition.
For other strategies, this needs to be modified.
=#
function exchange_ghosts(var, grid, tag=0)
    tag = config.num_procs + Int(round(tag));
    if config.solver_type == FV
        return exchange_ghosts_fv(var, grid, tag);
    else
        printerr("Ghost exchange not ready for "*config.solver_type, fatal=true)
    end
end

function exchange_ghosts_fv(var, grid, tag)
    n_neighbors = grid.num_neighbor_partitions;
    # How big does each buffer need to be?
    buffer_size = zeros(Int, n_neighbors);
    send_arrays = fill(zeros(0), n_neighbors);
    recv_arrays = fill(zeros(0), n_neighbors);
    requests = Vector{MPI.Request}(undef, n_neighbors * 2);
    if typeof(var) <: Array
        # multiple variables being solved for simultaneously
        dofs_per_node = 0;
        for vi=1:length(var)
            dofs_per_node += var[vi].total_components;
        end
    else
        # one variable
        dofs_per_node = var.total_components;
    end
    for ni=1:n_neighbors
        buffer_size[ni] = dofs_per_node * grid.ghost_counts[ni]; # number of values to send to neighbor ni
        send_arrays[ni] = zeros(buffer_size[ni]);
        recv_arrays[ni] = zeros(buffer_size[ni]);

        # Fill the send buffer
        next_ind = 1;
        if typeof(var) <: Array
            for vi=1:length(var)
                for ei=1:grid.ghost_counts[ni]
                    for comp_i=1:var[vi].total_components
                        send_arrays[ni][next_ind] = var[vi].values[comp_i,grid.ghost_index[ni][2,ei]];
                        next_ind += 1;
                    end
                end
            end
        else
            for ei=1:grid.ghost_counts[ni]
                for comp_i=1:var.total_components
                    send_arrays[ni][next_ind] = var.values[comp_i,grid.ghost_index[ni][2,ei]];
                    next_ind += 1;
                end
            end
        end
    end

    # send ghosts
    if config.use_mpi
        for ni=1:grid.num_neighbor_partitions
            requests[ni*2-1] = MPI.Isend(send_arrays[ni], grid.neighboring_partitions[ni], tag, MPI.COMM_WORLD);
        end
    end
    # receive ghosts
    if config.use_mpi
        for ni=1:grid.num_neighbor_partitions
            requests[ni*2] = MPI.Irecv!(recv_arrays[ni], grid.neighboring_partitions[ni], tag, MPI.COMM_WORLD);
        end

        # wait for it
        stats = MPI.Waitall!(requests);
    end

    # Put the received ghost values in place
    for ni=1:n_neighbors
        next_ind = 1;
        if typeof(var) <: Array
            for vi=1:length(var)
                for ei=1:grid.ghost_counts[ni]
                    for comp_i=1:var[vi].total_components
                        var[vi].values[comp_i,grid.ghost_index[ni][1,ei]] = recv_arrays[ni][next_ind];
                        next_ind += 1;
                    end
                end
            end
        else
            for ei=1:grid.ghost_counts[ni]
                for comp_i=1:var.total_components
                    var.values[comp_i,grid.ghost_index[ni][1,ei]] = recv_arrays[ni][next_ind];
                    next_ind += 1;
                end
            end
        end
    end
end

max_stars_repo_name: paralab/Finch

max_issues_repo_name: paralab/Finch

max_forks_repo_name: paralab/Finch


text:
# Copyright (c) Microsoft Corporation.
# Licensed under the MIT license.

module TestUtils
include("../../../../src/julia/shared/load.jl")
include("../../../../src/julia/runner/load.jl")
using ADPerfTest
import Benchmark

export can_objective_run_multiple_times!

"""
    can_objective_run_multiple_times!(context, objective)

Checks whether `objective` function of some `Test<Input, Output>` instance runs for different time when the supplied
times parameter is different. To do so this function uses `Benchmark.find_repeats_for_minimum_measurable_time`
function. It tries to find such a minimum_measurable_time that `find_repeats_for_minimum_measurable_time` will
return a number of repeats other than `1` or `measurable_time_not_achieved`.
If `objective` ignores its times parameter, we won't be able to find it.

# Arguments
- `context::Any`: `context` field of some `Test<Input, Output>` object.
- `objective::Function`: either `calculate_objective!` or `calculate_jacobian!` field of
the same `Test<Input, Output>` object.
"""
function can_objective_run_multiple_times!(context::Any, objective::Function)
    minimum_measurable_time = 0.05
    repeats, _, total_time = Benchmark.find_repeats_for_minimum_measurable_time!(context, minimum_measurable_time, objective)
    while (repeats == 1)
        # minimum_measurable_time * 2 ensures, that minimum_measurable_time * 2 will grow, while
        # result.total_time * 2 is a good guess for the time needed for at least 2 repeats
        minimum_measurable_time = max(minimum_measurable_time * 2, total_time * 2);
        repeats, _, total_time = Benchmark.find_repeats_for_minimum_measurable_time!(context, minimum_measurable_time, objective)
    end
    repeats ≠ Benchmark.measurable_time_not_achieved
end

end

max_stars_repo_name: dsyme/ADBench

max_issues_repo_name: dsyme/ADBench

max_forks_repo_name: dsyme/ADBench


text:
include("idl_types.jl")
include("common-funcs.jl")
include("common-macros.jl")

# RPC client
type RPCclient
    ptr::Ptr{Void}
    process::Union{Void, Base.Process}
end
RPCclient() = RPCclient(C_NULL, nothing)
RPCclient(ptr::Ptr{Void}) = RPCclient(ptr, nothing)

pclient = RPCclient()

function rpc_init()
    ccall((:IDL_RPCInit, idlrpc), Ptr{Void}, (Clong, Ptr{UInt8}), 0, C_NULL)
end

function rpc_cleanup()
    ecode = ccall((:IDL_RPCCleanup, idlrpc), Cint, (Ptr{Void}, Cint), pclient.ptr, 0)
    if ecode != 1
        error("IDL.exit: failed")
    end
    if pclient.process != nothing
        kill(pclient.process)
    end
    return
end

function init()
    olderr = STDERR
    (rd, wr) = redirect_stderr() #Redirect error messages
    ptr = rpc_init()
    if ptr != C_NULL #Check if idlrpc is already running
        global pclient = RPCclient(ptr)
    else #Start up idlrpc
        print("Initializing IDL ")
        proc = spawn(`idlrpc`)
        ptr = C_NULL
        cnt = 0
        while ptr == C_NULL && cnt < 60 #Allow for startup time
            ptr = rpc_init()
            cnt = cnt + 1
            sleep(1)
            print(".")
        end
        println("")
        ptr == C_NULL && error("IDL.init: IDLRPC init failed")
        global pclient = RPCclient(ptr, proc)
    end
    capture(true)
    redirect_stderr(olderr)
    #Register cleanup function to be called at exit
    atexit(rpc_cleanup)
end

function execute_converted(str::AbstractString)
    # does no conversion of interpolated vars, continuation chars, or newlines
    ecode = ccall((:IDL_RPCExecuteStr, idlrpc), Cint, (Ptr{Void},Ptr{UInt8}), pclient.ptr, str)
    if ecode != 1
        # since error get printed by IDL, we just reset error state
        ecode = ccall((:IDL_RPCExecuteStr, idlrpc), Cint, (Ptr{Void},Ptr{UInt8}), pclient.ptr,
                      "message, /RESET")
        flush()
        return false
    end
    flush()
    return true
end

function capture(flag::Bool)
    nlines = flag ? 5000 : 0
    ecode = ccall((:IDL_RPCOutputCapture, idlrpc), Cint, (Ptr{Void}, Cint), pclient.ptr, nlines)
    if ecode != 1
        error("IDL.capture: IDL_RPCOutputCapture failed")
    end
    return nothing
end

function get_output!(line_s::IDL_RPC_LINE_S)
    ecode = ccall((:IDL_RPCOutputGetStr, idlrpc), Cint, (Ptr{Void},Ref{IDL_RPC_LINE_S},Cint),
                  pclient.ptr, line_s, 0)
    ecode == 1 ? true : false
end

flag_set(x, flag) = (x & flag) == flag

function flush()
    line_s = IDL_RPC_LINE_S()
    # gc()
    while get_output!(line_s)
        println(unsafe_string(line_s.buf))
        flag_set(line_s.flags, IDL_TOUT_F_NLPOST) && print("\n")
    end
end

# no free_cb needed in idlrpc (I think)
free_cb = C_NULL

# NOTE: Put_var makes a copy of the data in the array when it is put into idlrpc process.
#       This is different than callable idl where the pointer to the data is copied.
#       I think the difference is because callable idl runs in the same process but
#       idlrpc does not. Thus, no free_cb is needed in idlrpc version.
function put_var{T,N}(arr::Array{T,N}, name::AbstractString)
    if !isbits(eltype(arr)) || (idl_type(arr) < 0)
        error("IDL.put_var: only works with some vars containing bits types")
    end
    dim = zeros(Int, IDL_MAX_ARRAY_DIM)
    dim[1:N] = [size(arr)...]
    vptr = ccall((:IDL_RPCImportArray, idlrpc), Ptr{IDL_Variable},
                 (Cint, IDL_ARRAY_DIM, Cint, Ptr{T}, IDL_ARRAY_FREE_CB),
                 N, dim, idl_type(arr), arr, free_cb)
    ecode = ccall((:IDL_RPCSetVariable, idlrpc), Cint,
                  (Ptr{Void}, Ptr{UInt8}, Ptr{IDL_Variable}),
                  pclient.ptr, name, vptr)
    if ecode != 1
        error("IDL.put_var: failed")
    end
    return
end

# there must be a slicker way to do this?
function uint_size(x)
    if sizeof(x) == 1
        UInt8
    elseif sizeof(x) == 2
        UInt16
    elseif sizeof(x) == 4
        UInt32
    elseif sizeof(x) == 8
        UInt64
    elseif sizeof(x) == 16
        UInt128
    end
end

function put_var(x, name::AbstractString)
    if !isbits(x) || (idl_type(x) < 0)
        error("IDL.put_var: only works with some vars containing bits types")
    end
    vptr = get_vptr(name)
    ccall((:IDL_RPCStoreScalar, idlrpc), Void,
          (Ptr{IDL_Variable}, Cint, Ref{UInt128}),
          vptr, idl_type(x), Ref{UInt128}(convert(UInt128,reinterpret(uint_size(x),x))))
    ecode = ccall((:IDL_RPCSetVariable, idlrpc), Cint,
                  (Ptr{Void}, Ptr{UInt8}, Ptr{IDL_Variable}),
                  pclient.ptr, name, vptr)
    return
end

function put_var{T}(x::Complex{T}, name::AbstractString)
    T == Float32 || T == Float64 || error("IDL.put_var: only floating point complex types allowed")
    if T == Float64
        y = (convert(UInt128, reinterpret(UInt64, imag(x))) << 64) + reinterpret(UInt64, real(x))
    else
        y = (convert(UInt128, reinterpret(UInt32, imag(x))) << 32) + reinterpret(UInt32, real(x))
    end
    vptr = get_vptr(name)
    ccall((:IDL_RPCStoreScalar, idlrpc), Void,
          (Ptr{IDL_Variable}, Cint, Ref{UInt128}),
          vptr, idl_type(x), Ref{UInt128}(y))
    ecode = ccall((:IDL_RPCSetVariable, idlrpc), Cint,
                  (Ptr{Void}, Ptr{UInt8}, Ptr{IDL_Variable}),
                  pclient.ptr, name, vptr)
    return
end

# function put_var(str::AbstractString, name::AbstractString)
#     # Sort of a HACK: do direcly since ImportNamedArray doesn't work
#     execute("$name = '$str'")
#     return
# end

function put_var(str::AbstractString, name::AbstractString)
    idl_string = IDL_String()
    ccall((:IDL_RPCStrStore, idlrpc), Void, (Ptr{IDL_String}, Ptr{Cchar}),
          idl_string, str)
    ecode = ccall((:IDL_RPCSetVariable, idlrpc), Cint,
                  (Ptr{Void}, Ptr{UInt8}, Ptr{IDL_Variable}),
                  pclient.ptr, name, vptr)
    if ecode != 1
        error("IDL.put_var: failed")
    end
end

function get_name(vptr::Ptr{IDL_Variable})
    # not implemented for RPC
    str = ccall((:IDL_RPCVarName, idlrpc), Ptr{Cchar}, (Ptr{IDL_Variable},), vptr)
    return unsafe_string(str)
end

function get_vptr(name::AbstractString)
    # returns C_NULL if name not in scope
    ccall((:IDL_RPCGetMainVariable, idlrpc), Ptr{IDL_Variable}, (Ptr{Void},Ptr{UInt8}),
          pclient.ptr, name)
end

function get_var(name::AbstractString)
    vptr = ccall((:IDL_RPCGetMainVariable, idlrpc), Ptr{IDL_Variable}, (Ptr{Void},Ptr{UInt8}),
                 pclient.ptr, name)
    # NOTE: IDL_RPCGetVariable never seems to return NULL in spite of docs
    if vptr == C_NULL
        error("IDL.get_var: variable $name does not exist")
    end
    var = get_var(vptr, name)
    # not sure if this is needed?
    vptr = ccall((:IDL_RPCDeltmp, idlrpc), Void, (Ptr{IDL_Variable},), vptr)
    var
end


text:
function comp_loc(x,y, dat_col)
    chr_start = split(x)[1]
    chr_start_2 = split(y)[1]
    if chr_start != chr_start_2
        return "Unique"
    end
    #println(x)
    pos_start = split(x)[2]
    #pos_end = replace(split(split(x)[8],";")[1],"END=" => "")
    pos_end = find_pos_val("END=",x)
    #svlen = replace(split(split(x)[8],";")[4],"SVLEN=" => "")
    svlen = find_pos_val("SVLEN=",x)
    svtype = find_pos_val("SVTYPE=",x)
    if occursin("INV",svtype)
    elseif occursin("DEL",svtype)#svtype != "DEL"
    elseif occursin("DUP:TANDEM",svtype)
    else
        pos_end = string(parse(Int, pos_end) + parse(Int,svlen))
    end
    #=if svtype != "DEL"
        pos_end = string(parse(Int, pos_end) + parse(Int,svlen))
    end=#

    pos_start_2 = split(split(split(split(y)[dat_col],":")[11],"-")[1],"_")[end]
    pos_end_2 = split(split(split(split(y)[dat_col],":")[11],"-")[2],"_")[end]
    if pos_start != pos_start_2 &&  pos_end != pos_end_2
        return "Unique"
    elseif pos_start == pos_start_2 && pos_end == pos_end_2
        return "Merged"
    else
        #println("error")
        #println(svtype*"\t"*pos_start*"\t"*pos_end*"\t"*pos_start_2*"\t"*pos_end_2)
        #println(x)
        return "singleton"
    end
end

function find_pos_val(val, instring)
    t= split(split(instring)[8],";")
    ind = [occursin(val,i) for i in t]
    y= replace(t[ind][1],val => "")
    return y
end

"""
    check_pos(x:String, y::String)::Bool

Compare two input SV lines to make sure that the position of the first line is lower than or equal to the other

"""
function check_pos(x,y)
    chr_start = split(x)[1]
    pos_1 = split(x)[2]
    chr_start_2 = split(y[1])[1]
    pos_2 = split(y[1])[2]
    if parse(Int,pos_1 ) <= parse(Int,pos_2) && chr_start == chr_start_2
        return true
    elseif parse(Int,pos_1 ) >  parse(Int,pos_2) && chr_start != chr_start_2
        return true
    else
        return false
    end
end

function compare_break_points(x,y,comp_col)
    pos_start_X = split(x)[2]
    #pos_end = replace(split(split(x)[8],";")[1],"END=" => "")
    pos_end_x = find_pos_val("END=",x)
    svlen_x = find_pos_val("SVLEN=",x)
    svtype_x = find_pos_val("SVTYPE=",x)
    if svtype_x != "INV"
        pos_end_x = string(parse(Int, pos_end_x) + parse(Int,svlen_x))
    end

    pos_start_y = split(split(split(split(y)[comp_col],":")[11],"-")[1],"_")[end]
    pos_end_y = split(split(split(split(y)[comp_col],":")[11],"-")[2],"_")[end]
    if pos_start_x == pos_start_y && pos_end_x == pos_end_y
        bk_pt_overlap ="Both"
    elseif pos_start_x == pos_start_y && pos_end_x != pos_end_y
        bk_pt_overlap = "Left"
    elseif pos_end_x == pos_end_y
        bk_pt_overlap = "Right"
    else
        bk_pt_overlap = "Ranged"
    end
    consensus = string(split(split(y)[11],":")[10] == split(split(y)[10],":")[10])
    bk_pt_overlap = bk_pt_overlap*"\t"*consensus*"\t"*replace(y,"\t" => "|")
    return bk_pt_overlap
end

function ann_sv(sv,l1, l2, comp_col)
    if sv =="Unique"
        return l1*"\t"*sv*"\tNA\tNA\tNA"
    else
        overlap = compare_break_points(l1,l2[1],comp_col)
        #println("Done3")
        return l1*"\t"*sv*"\t"*overlap
    end
end

function annotate_pindel(v1,v2,dat_col,comp_col)
    check_line = iterate(v2)
    v3=[]
    #i = 1
    #j = 1
    for line1 in v1
        if check_pos(line1,check_line)
            sv = comp_loc(line1,check_line[1],dat_col)
        elseif iterate(v2, check_line[2]) !== nothing
            check_line = iterate(v2,check_line[2])
            #j +=1
            sv= comp_loc(line1,check_line[1],dat_col)
        else
            sv= comp_loc(line1,check_line[1],dat_col)
        end
        val = ann_sv(sv,line1,check_line,comp_col)
        #println(val)
        push!(v3,val)
        #i +=1
        #println(string(i)*" " *string(j))
    end
    return v3
end

function print_res(filePath,vec)
    fheader = "CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\tFORMAT\tGT\tAD_REF\tAD_ALT"
    fheader = fheader*"\tStatus_of_merge\tbreak_point_kind\tSeq_overlap_status\tResults_of_merge"
    fheader = fheader*"\tSV_Type\tSV_length\tHompolymer_length\tfrac_ALT"
    f = open(filePath,"w")
    println(f,fheader)
    for line in vec
        tmpline = split(line)
        tmpline[10] =replace(tmpline[10],":" => "\t")
        tmpline[10] =replace(tmpline[10],"," => "\t")
        AD_ref = parse(Float64,split(tmpline[10])[2])
        AD_alt = parse(Float64,split(tmpline[10])[3])

        tmpline_to_print = tmpline
        tmpline_to_print[4] = first(tmpline_to_print[4],10)
        tmpline_to_print[5] = first(tmpline_to_print[5],10)
        tmpline_to_print[14] = first(tmpline_to_print[14],10)

        push!(tmpline_to_print, find_pos_val("SVTYPE=",line))
        push!(tmpline_to_print, find_pos_val("SVLEN=",line))
        push!(tmpline_to_print, find_pos_val("HOMLEN=",line))
        push!(tmpline_to_print, string(AD_alt/(AD_alt+AD_ref)))
        println(f,join(tmpline_to_print,"\t"))
    end
    close(f)
end

datdir = "/Users/aragaven/scratch/test_julia/test_pindel/pindel_vcf";
samp1 = "DSY303_OP";
samp2 = "DSY303_WH";
merged_file = "/Users/aragaven/Documents/Research/CBC/ene/DSY303_merged.vcf";
vec1 = readlines(joinpath(datdir,samp1*".vcf"));
vec2 = readlines(joinpath(datdir,samp2*".vcf"));
vec3 = readlines(merged_file);
vec1_1 = filter(x -> !startswith(x, "#"), vec1);
vec2_2 = filter(x -> !startswith(x, "#"), vec2);
vec3_3 = filter(x -> !startswith(x, "#"), vec3);
vec4 = annotate_pindel(vec1_1,vec3_3,10,11);
vec5 = annotate_pindel(vec2_2,vec3_3,11,10);

filep1 = "/Users/aragaven/Documents/Research/CBC/ene/DSY303_OP.results.tsv"
print_res(filep1,vec4)
filep2 = "/Users/aragaven/Documents/Research/CBC/ene/DSY303_WH.results.tsv"
print_res(filep2,vec5)

max_stars_repo_name: compbiocore/BioJulia-cookbook

max_issues_repo_name: compbiocore/BioJulia-cookbook

max_forks_repo_name: compbiocore/BioJulia-cookbook


text:
using SymPy

opmap = Dict(:(*) => "*",
             :(.*) => "*",
             :(/) => "/",
             :(./) => "/",
             :(+) => "+",
             :(-) => "-")

priomap = Dict(:(*) => 2, :(.*) => 2, :(/) => 2, :(./) => 2, :(+) => 1, :(-) => 1)

function js_code(ex::Expr, prio, pop, freevars)
  if ex.head == :call
    op = ex.args[1]
    args = ex.args[2:end]

    if op == :(-) && length(args) == 1
      return string("-", js_code(args[1], 3, 0, freevars))
    elseif haskey(opmap, op)
      nprio = priomap[op]

      ret = join([js_code(args[i], nprio, i == 1 ? 0 : op, freevars) for i in 1:length(args)], string(" ", opmap[op], " "))

      if nprio < prio || (pop != 0 && op != pop && nprio == prio)
        ret = string("(", ret, ")")
      end

      return ret
    elseif op == :(^) || op == :(.^)
      if isa(args[1], Symbol) && isa(args[2], Number)
        expanded = Expr(:call, :(*), [args[1] for i in 1:args[2]]...)
        return js_code(expanded, prio, :(^), freevars)
      else
        return string("Math.pow(", js_code(args[1], 0, 0, freevars), ", ", js_code(args[2], 0, 0, freevars), ")")
      end
    elseif op == :log
      return string("Math.log(", js_code(args[1], 0, 0, freevars) , ")")
    else
      return string(ex)
    end
  elseif ex.head == :if
    cond = js_code(ex.args[1], 0, 0, freevars)
    truepart = js_code(ex.args[2], 0, 0, freevars)
    falsepart = js_code(ex.args[3], 0, 0, freevars)

    return string(cond, " ? ", truepart, " : ", falsepart)
  elseif ex.head == :comparison
    lhs = js_code(ex.args[1], 3, 0, freevars)
    rhs = js_code(ex.args[3], 3, 0, freevars)

    return string(lhs, " ", ex.args[2], " ", rhs)
  else
    return string(ex)
  end
end

function js_code(sym::Symbol, prio, op, freevars)
  ret = string(sym)

  if ret == "ln2"
    return "Math.LN2"
  end

  if !(ret in freevars)
    push!(freevars, ret)
  end

  if ismatch(r"^x[0-9]+$", ret) || ismatch(r"_$", ret) || ret == "apexDistance"
    return ret
  elseif ret in ("ascensionFactor", "descensionFactor")
    return string("p.", ret)
  else
    return string("d.", ret)
  end
end

js_code(sym::Sym, prio, op, freevars) = js_code(parse(julia_code(sym)), prio, op, freevars)
js_code(ex, prio, op, freevars) = string(ex)

js_code(ex) = js_code(ex, 0, 0, [])

js_code(terms::Tuple) = js_code(fname, terms, [])

function js_code(terms::Tuple, paramvars)
  terms, exprs = terms
  expr = exprs[1]

  body = []
  termnames = []
  freevars = []

  for term in terms
    push!(termnames, string(term[1]))
    push!(body, string("  const ", term[1], " = ", js_code(term[2], 0, 0, freevars), ";"))
  end

  if length(terms) > 0
    push!(body, "")
  end

  push!(body, string("  return ", js_code(expr, 0, 0, freevars), ";"))

  freevars = filter(x -> !(x in termnames), freevars)

  for var in freevars
    if !(var in paramvars)
      push!(paramvars, var)
    end
  end

  return join(body, "\n")
end

max_stars_repo_name: rihorn/ky-libaries

max_issues_repo_name: rihorn/ky-libaries

max_forks_repo_name: rihorn/ky-libaries


text:
"""
    graph_small_world(n::Integer, m::Integer, p::Real)

Small-world network model of Watts and Strogatz (1998).

# Arguments
* `n` — number of nodes,
* `k` — number of neighbors for each node, should be even,
* `p` — probability of random rewiring of each edge.

# References
Watts, D. J., & Strogatz, S. H. (1998).
Collective dynamics of “small-world” networks. Nature, 393(6684), 440–442.
"""
function graph_small_world(n::Integer, k::Integer, p::Real; kwargs...)
    @assert(n > 5,         "number of nodes (n) should be > 5")
    @assert(0 < k < n / 2, "number of neighbors (k) is not in a valid range")
    @assert(k % 2 == 0,    "number of neighbors (k) should be even")
    @assert(0 < p < 1,     "rewiring probability (p) is not in (0, 1) interval")
    g = Graph(n; kwargs...)
    for i = 1:n, offset = 1:Int(k / 2)
        j = (i + offset - 1) % n + 1
        if rand() < p || hasedge(g, i, j)
            while true
                j = rand(1:n)
                if j != i && !hasedge(g, i, j)
                    break
                end
            end
        end
        addedge!(g, i, j)
    end
    return g
end

@deprecate(graph_smallworld, graph_small_world)


text:
using Documenter
using Literate
using TaskMaster

Literate.markdown(joinpath(@__DIR__, "../examples/introduction.jl"), joinpath(@__DIR__,"src/"); credit = false, name = "index")

Literate.markdown(joinpath(@__DIR__, "../examples/adaptive.jl"), joinpath(@__DIR__,"src/"); credit = false, name = "adaptive")

makedocs(sitename="TaskMaster.jl",pages = ["index.md","adaptive.md"])

deploydocs(
     repo = "github.com/akels/TaskMaster.jl.git",
 )


text:
using MPI

function runmpi(tests, file)
  MPI.Initialized() && !MPI.Finalized() &&
  error("runmpi does not work if MPI has been "*
        "Initialized but not Finalizd")

  # The code below was modified from the MPI.jl file runtests.jl
  #
  # Code coverage command line options; must correspond to src/julia.h
  # and src/ui/repl.c
  JL_LOG_NONE = 0
  JL_LOG_USER = 1
  JL_LOG_ALL = 2
  coverage_opts = Dict{Int, String}(JL_LOG_NONE => "none",
                                    JL_LOG_USER => "user",
                                    JL_LOG_ALL => "all")
  coverage_opt = coverage_opts[Base.JLOptions().code_coverage]
  testdir = dirname(file)

  if haskey(ENV, "SLURM_JOB_ID")
    oversubscribe = `--oversubscribe`
  else
    oversubscribe = ``
  end

  for (n, f) in tests
    cmd = `mpiexec $oversubscribe -n $n $(Base.julia_cmd()) --startup-file=no --project=$(Base.active_project()) --code-coverage=$coverage_opt $(joinpath(testdir, f))`

    @info "Running MPI test..." n f cmd
    # Running this way prevents:
    #   Balance Law Solver | No tests
    # since external tests are not returned as passed/fail
    @time @test (run(cmd); true)
  end
end

max_stars_repo_name: akshaysridhar/CLIMA

max_issues_repo_name: akshaysridhar/CLIMA

max_forks_repo_name: akshaysridhar/CLIMA


text:
using JSON
using BytePairEncoding
using BytePairEncoding: UnMap
using Transformers.Basic

abstract type GPT2 <: PretrainedTokenizer end

# wrapper for GPT2 Tokenizer with required functionalities
"""
struct GPT2Tokenizer <: GPT2
    encoder::Vocabulary{String}
    bpe_encode::GenericBPE
    bpe_decode::UnMap
    vocab::Dict{String, Any}
    unk_token::String
    unk_id::Int
    eos_token::String
    eos_token_id::Int
    pad_token::String
    pad_token_id::Int
end

Structure to hold all essential information / functions for GPT2 tokenizer

"""
struct GPT2Tokenizer <: GPT2
    encoder::Vocabulary{String}
    bpe_encode::GenericBPE
    bpe_decode::UnMap
    vocab::Dict{String, Any}
    unk_token::String
    unk_id::Int
    eos_token::String
    eos_token_id::Int
    pad_token::String
    pad_token_id::Int
end

"""
    load_pretrained_tokenizer(ty::Type{T}; unk_token="<|endoftext|>", eos_token="<|endoftext|>", pad_token="<|endoftext|>") where T<:PretrainedTokenizer

Load GPT2 tokenizer using Datadeps for pretrained bpe and vocab. Returns tokenizer as `GPT2Tokenizer` structure.

"""
function load_pretrained_tokenizer(ty::Type{T}; unk_token="<|endoftext|>", eos_token="<|endoftext|>", pad_token="<|endoftext|>") where T<:PretrainedTokenizer
    path_bpe = joinpath(datadep"BPE", readdir(datadep"BPE")[1])
    path_vocab = joinpath(datadep"Vocab", readdir(datadep"Vocab")[1])
    load_pretrained_tokenizer(path_bpe, path_vocab, unk_token, eos_token, pad_token)
end

"""
    load_pretrained_tokenizer(path_bpe, path_vocab, unk_token, eos_token, pad_token)

Load pretrained tokenizer for GPT2 from provided bpe and vocab file path. Initialises `unk_token`, `eos_token`, `pad_token` as provided with the function. Returns tokenizer as `GPT2Tokenizer` structure.

"""
function load_pretrained_tokenizer(path_bpe, path_vocab, unk_token, eos_token, pad_token)
    vocab = JSON.parsefile(path_vocab)
    labels = map(x->x[1], sort!(collect(pairs(vocab)), by=x->x[2]))
    encoder = Vocabulary(labels, unk_token)
    bpe_encode = ByteLevelBPE(path_bpe)
    bpe_decode = BytePairEncoding.UnMap(bpe_encode.codemap)
    GPT2Tokenizer(encoder, bpe_encode, bpe_decode, vocab, unk_token, encoder(unk_token),
            eos_token, encoder(eos_token), pad_token, encoder(eos_token))
end

"""
    tokenize(t::GPT2Tokenizer, text::AbstractString)

Function to tokenize given `text` with tokenizer bpe encoder (`t.bpe_encode`). Returns a string vector of tokens.
"""
function tokenize(t::GPT2Tokenizer, text::AbstractString)
    t.bpe_encode(text)
end

"""
    encode(t::GPT2Tokenizer, text::AbstractString; add_prefix_space=false)

Returns the encoded vector of tokens (mapping from vocab of Tokenizer) for `text`. If `add_prefix_space`=true, add space at the start of 'text' before tokenization.

# Example

For single text:

```julia
encode(tokenizer, text)
```

For vector of text:

```julia
map(x->encode(tokenizer, x), text_vector)
```
"""
function encode(t::GPT2Tokenizer, text::AbstractString; add_prefix_space=false)
    if add_prefix_space==true
        text = string(" ", text)
    end
    tokens = tokenize(t, text)
    t.encoder(tokens)
end

"""
    encode(t::GPT2Tokenizer, tokens::Vector{String})

Function to encode tokens vectors to their integer mapping from vocab of tokenizer.
"""
function encode(t::GPT2Tokenizer, tokens::Vector{String})
    t.encoder(tokens)
end

"""
    (t::GPT2Tokenizer)(text::AbstractString; add_prefix_space=false)

Encode the text with tokenizer and returns the encoded vector. If `add_prefix_space`=true, add space at the start of 'text' before tokenization.

# Examples:

For a single text:
```julia
tokenizer(text; add_prefix_space=true)
```

For vector of texts, use:

```julia
map(x->encode(tokenizer, x), text_vector)

# or

tokenizer.(text_vector)
```

Also checkout [`encode`](@ref PPLM.encode)
"""
function (t::GPT2Tokenizer)(text::AbstractString; add_prefix_space=false)
    encode(t, text; add_prefix_space=add_prefix_space)
end

decode(vocab::Vocabulary{T}, i::Int) where T = 0 <= i <= length(vocab) ? vocab.list[i] : vocab.unk

"""
    decode(vocab::Vocabulary{T}, is::Vector{Int}) where T

Return decoded vector of `string` tokens from the indices vector `is`, using the vocab.
"""
function decode(vocab::Vocabulary{T}, is::Vector{Int}) where T
    tokens = Vector{String}(undef, length(is))
    for (idx, i) ∈ enumerate(is)
        token = decode(vocab, i)
        tokens[idx] = token
    end
    tokens
end

"""
    decode(t::GPT2Tokenizer, tokens_ids::Vector{Int})

Return decoded vector of `string` tokens from the indices vector `tokens_ids`, using the tokenizer `t` encoder .
"""
function decode(t::GPT2Tokenizer, tokens_ids::Vector{Int})
    decode(t.encoder, tokens_ids)
end

"""
    detokenize(t::GPT2Tokenizer, tokens::Vector{String})

BPE Decode the vector of strings, using the tokenizer `t`.
"""
function detokenize(t::GPT2Tokenizer, tokens::Vector{String})
    t.bpe_decode(join(tokens))
end

# Example for detokenizing multiple examples: map(token_ids->detokenize(tokenizer, token_ids), tokens)
"""
    detokenize(t::GPT2Tokenizer, tokens_ids::Vector{Int})

Decode and Detokenize the vector of indices `token_ids`. Returns the final sentence after detokenization.

# Example

For single vector of token_ids:

```julia
detokenize(tokenizer, token_ids)
```

For vector of vector of `token_ids`, use:

``` julia
map(x->decode(tokenizer, x), tokens_id_vector_of_vector)
```
"""
function detokenize(t::GPT2Tokenizer, tokens_ids::Vector{Int})
    tokens_list = decode(t, tokens_ids)
    detokenize(t, tokens_list)
end


text:
using Luna

a = 13e-6
gas = :H2
pres = 5
flength = 200e-2

τfwhm = 20e-15
λ0 = 800e-9
energy = 1e-6

grid = Grid.RealGrid(flength, λ0, (400e-9, 2000e-9), 40e-12)

m = Capillary.MarcatiliMode(a, gas, pres, loss=false)
aeff = let m=m
    z -> Modes.Aeff(m, z=z)
end

densityfun = let dens0=PhysData.density(gas, pres)
    z -> dens0
end

responses = (Nonlinear.Kerr_field(PhysData.γ3_gas(gas)),
             Nonlinear.RamanPolarField(grid.to, Raman.raman_response(grid.to, gas)))

linop, βfun!, frame_vel, αfun = LinearOps.make_const_linop(grid, m, λ0)

inputs = Fields.GaussField(λ0=λ0, τfwhm=τfwhm, energy=energy)

Eω, transform, FT = Luna.setup(grid, densityfun, responses, inputs, βfun!, aeff)

statsfun = Stats.default(grid, Eω, m, linop, transform; gas=gas)
output = Output.MemoryOutput(0, grid.zmax, 201, statsfun)

Luna.run(Eω, grid, linop, transform, FT, output)

##
Plotting.pygui(true)
Plotting.stats(output)
Plotting.prop_2D(output)
Plotting.time_1D(output, range(0,flength,length=5), trange=(-0.5e-12, 2e-12))
Plotting.spec_1D(output, range(0,flength,length=5), λrange=(500e-9, 1400e-9))


text:
operator(pa, ::S)=pa.opS
operator(pa, ::G)=pa.opG

"""
A linear operator corresponding to convolution with either s or g
"""
function create_operator(pa, attrib)
	global to

	fw=(y,x)->@timeit to "F" F!(y, x, pa, attrib)
	bk=(y,x)->@timeit to "Fadj" Fadj!(y, x, pa, attrib)

	return LinearMap{collect(typeof(pa).parameters)...}(fw, bk,
		  length(pa.optm.cal.d),  # length of output
		  ninv(pa, attrib),
		  ismutating=true, isposdef=true)
end

function F!(y, x, pa, attrib)
	x_to_model!(x, pa, attrib) # modify pa.optm.cal.s
	Conv.mod!(pa.optm.cal, Conv.D()) # modify pa.optm.cal.d
	copyto!(y,pa.optm.cal.d)
	return nothing
end

function Fadj!(y, x, pa, attrib)
	copyto!(pa.optm.ddcal,x)
	Fadj!(pa, y, x, pa.optm.ddcal, attrib)
	return nothing
end

function filt!(y,x,pa::P_bandpass)
	copyto!(pa.x,x)
	mul!(pa.dfreq, pa.dfftp, pa.x)
	for i in eachindex(pa.dfreq)
		@inbounds pa.dfreq[i] *= pa.filter[i]
	end
	mul!(pa.y, pa.difftp, pa.dfreq)
	copyto!(y,pa.y)
end

function create_operator(pa::P_bandpass)
	fw=(y,x)->filt!(y, x, pa)
	F=LinearMap(fw, fw,# length of output,
                 pa.nt,
                 pa.nt,
                 ismutating=true, isposdef=true)

end
#=

function bandlimit_operator(n)

	function bandlimit!(y, x,)
		responsetype = Bandpass(0.1,0.4; fs=1)
		designmethod = Butterworth(16)
		b=digitalfilter(responsetype, designmethod)
		DSP.filt!(y, b, x)
	end

	function adjtest()
		x=randn(size(F,2))
		y=randn(size(F,1))
		a=LinearAlgebra.dot(y,F*x)
		b=LinearAlgebra.dot(x,adjoint(F)*y)
		c=LinearAlgebra.dot(x, transpose(F)*F*x)
		println("adjoint test: ", a, "\t", b)
		return isapprox(a,b,rtol=1e-6)
	end

	adjtest()

	=#


text:
@testset "Incremental problem tests" begin
    @test incrementalProblem([1, 2, 3], 0.5, [0, 1, 1], getKnapsackConstraints([1, 2, 3], 3), KnapsackProblemDescriptor(3)) == ([1, 1, 0], 3.0)
end


text:
include("../imports.jl")

@testset ExtendedTestSet "`throttle` steps" begin
    cb = CustomCallback(Events.StepBegin, TrainingPhase) do learner
            throw(CancelStepException("test"))
    end
    # only cancel every second batch
    throttledcb = throttle(cb, Events.StepBegin, freq = 2)

    learner = testlearner(Recorder(), throttledcb, coeff = 3)
    epoch!(learner, TrainingPhase())
    # should have made 8 steps instead of 16
    @test learner.cbstate.history[TrainingPhase()].steps == 8
end


text:
module Analysis

using DemeNet: Intent, Deme, ID, AbstractID, Notary, Contract, Certificate, Intent, Consensus
using PeaceVote: AbstractChain, AbstractProposal, AbstractVote, AbstractBraid

using ..Types: Proposal, Vote
using ..Ledgers: load, Ledger

import PeaceVote: attest, voters

function attest(chain::Vector{Union{Certificate,Contract}},notary::Notary)
    messages = []
    for record in chain
        if typeof(record) <: Certificate
            intent = Intent(record,notary)
            push!(messages,intent)
        elseif typeof(record) <: Contract
            braid = Consensus(record,notary)
            push!(messages,braid)
        else
            error("The type $(typeof(record)) is not valid for braidchain")
        end
    end

    return messages
end

function addvoters!(voters::Set{ID},input,output)
    for i in input
        @assert i in voters
    end

    for o in output
        push!(voters,o)
    end
end

function voters!(voters::Set{ID},input,output)
    addvoters!(voters,input,output)

    for i in input
        pop!(voters,i)
    end

end

function voters!(voters::Set{ID},messages::Vector)
    for msg in messages
        if typeof(msg)==Intent{ID}
            push!(voters,msg.document)
        elseif typeof(msg) <: Consensus{T} where T<:AbstractBraid
            input = unique(msg.references)
            output = unique(msg.document.ids)
            @assert length(input)==length(output)
            voters!(voters,input,output)
        end
    end
end

# """
#     Returns a parrent braid
# """
# function parrent(id,braids::Vector{Braid})
#     for b in braids
#         if id in b.outputs
#             return b
#         end
#     end
# end

# """
#     Returns a child braid
# """
# function child(id,braids::Vector{Braid})
#     for b in braids
#         if id in b.inputs
#             return b
#         end
#     end
# end

proposals(messages::Vector) = findall(msg -> typeof(msg) <: Intent{T} where T<:AbstractProposal,messages)
votes(messages::Vector) = filter(msg -> typeof(msg) <: Intent{T} where T<:AbstractVote,messages)
#tickets(messages::Vector) = filter(msg -> typeof(msg)==Ticket,messages)

#function voters(proposal::Proposal,messages)
# function voters(messages,proposal::Proposal)
#     index = findfirst(item -> item==proposal,messages)
#     vset = Set()
#     voters!(vset,messages[1:index])
#     return vset
# end

# function voters(messages,option::AbstractOption)
#     index = findfirst(x -> typeof(x)<:Proposal && x.uuid==option.pid,messages)
#     vset = Set()
#     voters!(vset,messages[1:index])
#     return vset
# end

function voters(messages,index::Integer)
    vset = Set{ID}()
    voters!(vset,messages[1:index])
    return vset
end

### Index is the one by which one can filter the votes of the proposal
#voters(messages,option::AbstractOption) = voters(messages,option.index)

function voters(braidchain)
    vset = Set{ID}()
    voters!(vset,braidchain)
    return vset
end

function members(braidchain)
    mset = Set{ID}()
    for msg in braidchain
        if typeof(msg) <: Intent{T} where T<:AbstractID
            push!(mset,msg.document.id)
        end
    end
    return mset
end

function members(braidchain,ca)
    mset = Set{ID}()
    for msg in braidchain
        if typeof(msg) <: Intent{T} where T<:AbstractID
            @assert msg.reference in ca "certificate with reference=$(msg.reference) is not valid"
            push!(mset,msg.document.id)
        end
    end
    return mset
end

function allvoters(braidchain)
    vset = Set{ID}()
    for msg in braidchain
        if typeof(msg) <: ID
            push!(vset,msg.document.id)
        elseif typeof(msg) <: Braid
            addvoters!(vset,msg)
        end
    end

    return vset
end

### Need to improve code so it would dispatch on proposal type at the particular index.
function normalcount(index::Int,deme::Deme,ledger::Ledger) ### Here I could have
    #loadedledger = load(chain.ledger) #BraidChain(deme).records
    loadedledger = load(ledger)
    #messages = attest(loadedledger,chain.deme.notary)
    messages = attest(loadedledger,deme.notary)
    intentproposal = messages[index]
    @assert typeof(intentproposal) <: Intent{T} where T<:AbstractProposal "Not a proposal"

    proposal = intentproposal.document

    voters = Set{ID}()
    voters!(voters,messages[1:index])

    ispvote(msg) = typeof(msg) <: Intent{T} where T<:Vote && msg.reference in voters && msg.document.pid==index

    tally = zeros(Int,length(proposal.options))

    for msg in messages[end:-1:index]
        if ispvote(msg)
            tally[msg.document.vote] += 1
            pop!(voters,msg.reference)
        end
    end

    return tally
end

# preferentialcount(proposal::AbstractProposal,deme::BraidChain) = error("Not yet implemented")

# quadraticcount(proposal::AbstractProposal,deme::BraidChain) = error("Not yet implemented")

end


text:
function sumlengths(points)
    T = eltype(points[1])
    result = zeros(T, length(points))
    for i=1:length(points)
        i0 = max(i-1,1)
        p1, p2 = points[i0], points[i]
        if !(any(map(isnan, p1)) || any(map(isnan, p2)))
            result[i] = result[i0] + norm(p1-p2)
        else
            result[i] = result[i0]
        end
    end
    result
end

intensity_convert(intensity, verts) = intensity
function intensity_convert(intensity::VecOrSignal{T}, verts) where T
    if length(value(intensity)) == length(value(verts))
        GLBuffer(intensity)
    else
        Texture(intensity)
    end
end

#TODO NaNMath.min/max?
dist(a, b) = abs(a-b)
mindist(x, a, b) = min(dist(a, x), dist(b, x))
function gappy(x, ps)
    n = length(ps)
    x <= first(ps) && return first(ps) - x
    for j=1:(n-1)
        p0 = ps[j]
        p1 = ps[min(j+1, n)]
        if p0 <= x && p1 >= x
            return mindist(x, p0, p1) * (isodd(j) ? 1 : -1)
        end
    end
    return last(ps) - x
end
function ticks(points, resolution)
    Float16[gappy(x, points) for x = linspace(first(points),last(points), resolution)]
end

function _default(position::Union{VecTypes{T}, MatTypes{T}}, s::style"lines", data::Dict) where T<:Point
    pv = value(position)
    p_vec = if isa(position, GPUArray)
        position
    else
        const_lift(position) do p
            pvv = vec(p)
            if length(pvv) < 4 # geometryshader doesn't work with less then 4
                return [pvv..., fill(T(NaN), 4-length(pvv))...]
            else
                return pvv
            end
        end
    end

    @gen_defaults! data begin
        dims::Vec{2, Int32} = const_lift(position) do p
            sz = ndims(p) == 1 ? (length(p), 1) : size(p)
            Vec{2, Int32}(sz)
        end
        vertex              = p_vec => GLBuffer
        intensity           = nothing
        color_map           = nothing => Texture
        color_norm          = nothing
        color               = (color_map == nothing ? default(RGBA, s) : nothing) => GLBuffer
        thickness::Float32  = 2f0
        pattern             = nothing
        fxaa                = false
        preferred_camera    = :orthographic_pixel
        boundingbox         = const_lift(x-> GLBoundingBox(to_cpu_mem(x)), p_vec)
        indices             = const_lift(length, p_vec) => to_index_buffer
        shader              = GLVisualizeShader("fragment_output.frag", "util.vert", "lines.vert", "lines.geom", "lines.frag")
        gl_primitive        = GL_LINE_STRIP_ADJACENCY
        startend            = const_lift(p_vec) do vec
            l = length(vec)
            map(1:l) do i
                (i == 1 || isnan(vec[max(i-1, 1)])) && return Float32(0) # start
                (i == l || isnan(vec[min(i+1, l)])) && return Float32(1) # end
                Float32(2) # segment
            end
        end => GLBuffer
    end
    if pattern != nothing
        if !isa(pattern, Texture)
            if !isa(pattern, Vector)
                error("Pattern needs to be a Vector of floats. Found: $(typeof(pattern))")
            end
            tex = GLAbstraction.Texture(ticks(pattern, 100), x_repeat = :repeat)
            data[:pattern] = tex
        end
        @gen_defaults! data begin
            pattern_length = Float32(last(pattern))
            lastlen   = const_lift(sumlengths, p_vec) => GLBuffer
            maxlength = const_lift(last, lastlen)
        end
    end
    data[:intensity] = intensity_convert(intensity, vertex)
    data
end

to_points(x::Vector{LineSegment{T}}) where {T} = reinterpret(T, x, (length(x)*2,))

_default(positions::VecTypes{LineSegment{T}}, s::Style, data::Dict) where {T <: Point} =
    _default(const_lift(to_points, positions), style"linesegment"(), data)

buffer_or_float(x::Signal{<: Number}) = map(Float32, x)
buffer_or_float(x::Number) = Float32(x)
buffer_or_float(x) = GLBuffer(x) # buffer should be able to handle signals

function _default(positions::VecTypes{T}, s::style"linesegment", data::Dict) where T <: Point
    @gen_defaults! data begin
        vertex              = positions => GLBuffer
        color               = default(RGBA, s, 1) => GLBuffer
        thickness           = 2f0 => buffer_or_float
        shape               = RECTANGLE
        pattern             = nothing
        fxaa                = false
        indices             = const_lift(length, positions) => to_index_buffer
        preferred_camera    = :orthographic_pixel
        # TODO update boundingbox
        boundingbox         = const_lift(x-> GLBoundingBox(to_cpu_mem(x)), positions)
        shader              = GLVisualizeShader("fragment_output.frag", "util.vert", "line_segment.vert", "line_segment.geom", "lines.frag")
        gl_primitive        = GL_LINES
    end
    if !isa(pattern, Texture) && pattern != nothing
        if !isa(pattern, Vector)
            error("Pattern needs to be a Vector of floats")
        end
        tex = GLAbstraction.Texture(ticks(pattern, 100), x_repeat = :repeat)
        data[:pattern] = tex
        data[:pattern_length] = Float32(last(pattern))
    end
    data
end

function _default(positions::Vector{T}, range::Range, s::style"lines", data::Dict) where T <: AbstractFloat
    length(positions) != length(range) && throw(
        DimensionMismatsch("length of $(typeof(positions)) $(length(positions)) and $(typeof(range)) $(length(range)) must match")
    )
    _default(points2f0(positions, range), s, data)
end

function _default(
        geometry::TOrSignal{G}, s::style"lines", data::Dict
    ) where G<:GeometryPrimitive{3}
    wireframe(geometry, data)
end
function _default(
        geometry::TOrSignal{GLNormalMesh}, s::style"lines", data::Dict
    )
    wireframe(geometry, data)
end
function wireframe(
        geometry, data::Dict
    )
    points = const_lift(geometry) do g
        decompose(Point3f0, g) # get the point representation of the geometry
    end
    # Get line index representation
    indices = decompose(Face{2, GLIndex}, value(geometry))
    data[:indices] = reinterpret(GLuint, indices)
    _default(points, style"linesegment"(), data)
end

struct GridPreRender end

function (::GridPreRender)()
    glEnable(GL_CULL_FACE)
    glDepthMask(GL_FALSE)
    glCullFace(GL_BACK)
end

function _default(c::TOrSignal{T}, ::Style{:grid}, data) where T<:AABB
    @gen_defaults! data begin
        primitive::GLPlainMesh = c
        bg_color = RGBA{Float32}(0.99,0.99,0.99,1)
        grid_color = RGBA{Float32}(0.8,0.8,0.8,1)
        grid_thickness = Vec3f0(0.999)
        gridsteps = Vec3f0(5)
        shader = GLVisualizeShader("fragment_output.frag", "grid.vert", "grid.frag")
        boundingbox = c
        prerender = GridPreRender()
        postrender = () -> glDisable(GL_CULL_FACE);
    end
end

function line_indices(array)
    len = length(array)
    result = Array(GLuint, len*2)
    idx = 1
    for i=0:(len-3), j=0:1
        result[idx] = i+j
        idx += 1
    end
    result
    #GLuint[i+j for i=0:(len-3) for j=0:1] # on 0.5
end
"""
Fast, non anti aliased lines
"""
function _default(position::VecTypes{T}, s::style"speedlines", data::Dict) where T <: Point
    @gen_defaults! data begin
        vertex       = position => GLBuffer
        color_map    = nothing  => Vec2f0
        indices      = const_lift(line_indices, position) => to_index_buffer
        color        = (color_map == nothing ? default(RGBA{Float32}, s) : nothing) => GLBuffer
        color_norm   = nothing  => Vec2f0
        intensity    = nothing  => GLBuffer
        shader       = GLVisualizeShader("fragment_output.frag", "dots.vert", "dots.frag")
        gl_primitive = GL_LINES
    end
end


text:
export reset_default_options

########################### OptionsSparse ###########################
"""
    OptionsSparse

Options for sparse linear algebra.
"""
mutable struct OptionsSparse
    auto_reorder::Bool
    solver::String
    OptionsSparse() = new(true, "SparseLU")
end

########################### OptionsNewtonRaphson ###########################
"""
    OptionsNewtonRaphson_LineSearch

Options for the line search in the Newton Raphson algorithm.
"""
mutable struct OptionsNewtonRaphson_LineSearch
    c1::Float64
    ρ_hi::Float64
    ρ_lo::Float64
    iterations::Int64
    maxstep::Int64
    αinitial::Float64
    OptionsNewtonRaphson_LineSearch() = new(1e-4, 0.5, 0.1, 1000, 9999999, 1.0)
end

"""
    OptionsNewtonRaphson

Options for Newton Raphson related algorithms.
"""
mutable struct OptionsNewtonRaphson
    max_iter::Int64
    verbose::Bool
    rtol::Float64
    tol::Float64
    LM::Float64
    linesearch::Bool
    linesearch_options::OptionsNewtonRaphson_LineSearch
    OptionsNewtonRaphson() = new(100, false, 1e-12, 1e-12, 0.0, false, OptionsNewtonRaphson_LineSearch())
end

mutable struct OptionsTraining
    training::PyObject
    verbose::Bool
    OptionsTraining() = new(placeholder(false), true)
end

######################################################
"""
    Options

The container of all options in ADCME.
"""
mutable struct Options
    sparse::OptionsSparse
    newton_raphson::OptionsNewtonRaphson
    training::OptionsTraining
    Options() = new(OptionsSparse(), OptionsNewtonRaphson(), OptionsTraining())
end

"""
    reset_default_options()

Results the ADCME options to default.
"""
function reset_default_options()
    global options
    options = Options()
end